Chaste  Release::2018.1
MutableVertexMesh.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "MutableVertexMesh.hpp"
37 
38 #include "LogFile.hpp"
39 #include "UblasCustomFunctions.hpp"
40 #include "Warnings.hpp"
41 
42 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
44  std::vector<VertexElement<ELEMENT_DIM,SPACE_DIM>*> vertexElements,
45  double cellRearrangementThreshold,
46  double t2Threshold,
47  double cellRearrangementRatio,
48  double protorosetteFormationProbability,
49  double protorosetteResolutionProbabilityPerTimestep,
50  double rosetteResolutionProbabilityPerTimestep)
51  : mCellRearrangementThreshold(cellRearrangementThreshold),
52  mCellRearrangementRatio(cellRearrangementRatio),
53  mT2Threshold(t2Threshold),
54  mProtorosetteFormationProbability(protorosetteFormationProbability),
55  mProtorosetteResolutionProbabilityPerTimestep(protorosetteResolutionProbabilityPerTimestep),
56  mRosetteResolutionProbabilityPerTimestep(rosetteResolutionProbabilityPerTimestep),
57  mCheckForInternalIntersections(false),
58  mDistanceForT3SwapChecking(5.0)
59 {
60  // Threshold parameters must be strictly positive
61  assert(cellRearrangementThreshold > 0.0);
62  assert(t2Threshold > 0.0);
63  assert(protorosetteFormationProbability >= 0.0);
64  assert(protorosetteFormationProbability <= 1.0);
65  assert(protorosetteResolutionProbabilityPerTimestep >= 0.0);
66  assert(protorosetteResolutionProbabilityPerTimestep <= 1.0);
67  assert(rosetteResolutionProbabilityPerTimestep >= 0.0);
68  assert(rosetteResolutionProbabilityPerTimestep <= 1.0);
69 
70  // Reset member variables and clear mNodes and mElements
71  Clear();
72 
73  // Populate mNodes and mElements
74  for (unsigned node_index=0; node_index<nodes.size(); node_index++)
75  {
76  Node<SPACE_DIM>* p_temp_node = nodes[node_index];
77  this->mNodes.push_back(p_temp_node);
78  }
79  for (unsigned elem_index=0; elem_index<vertexElements.size(); elem_index++)
80  {
81  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_temp_vertex_element = vertexElements[elem_index];
82  this->mElements.push_back(p_temp_vertex_element);
83  }
84 
85  // If in 3D, then also populate mFaces
86  if (SPACE_DIM == 3)
87  {
88  // Use a std::set to keep track of which faces have been added to mFaces
89  std::set<unsigned> faces_counted;
90 
91  // Loop over mElements
92  for (unsigned elem_index=0; elem_index<this->mElements.size(); elem_index++)
93  {
94  // Loop over faces of this element
95  for (unsigned face_index=0; face_index<this->mElements[elem_index]->GetNumFaces(); face_index++)
96  {
97  VertexElement<ELEMENT_DIM-1, SPACE_DIM>* p_face = this->mElements[elem_index]->GetFace(face_index);
98 
99  // If this face is not already contained in mFaces, then add it and update faces_counted
100  if (faces_counted.find(p_face->GetIndex()) == faces_counted.end())
101  {
102  this->mFaces.push_back(p_face);
103  faces_counted.insert(p_face->GetIndex());
104  }
105  }
106  }
107  }
108 
109  // Register elements with nodes
110  for (unsigned index=0; index<this->mElements.size(); index++)
111  {
112  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_temp_vertex_element = this->mElements[index];
113  for (unsigned node_index=0; node_index<p_temp_vertex_element->GetNumNodes(); node_index++)
114  {
115  Node<SPACE_DIM>* p_temp_node = p_temp_vertex_element->GetNode(node_index);
116  p_temp_node->AddElement(p_temp_vertex_element->GetIndex());
117  }
118  }
119 
120  this->mMeshChangesDuringSimulation = true;
121 }
122 
123 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
125  : mCellRearrangementThreshold(0.01),
126  mCellRearrangementRatio(1.5),
127  mT2Threshold(0.001),
128  mProtorosetteFormationProbability(0.0),
129  mProtorosetteResolutionProbabilityPerTimestep(0.0),
130  mRosetteResolutionProbabilityPerTimestep(0.0),
131  mCheckForInternalIntersections(false),
132  mDistanceForT3SwapChecking(5.0)
133 {
134  // Note that the member variables initialised above will be overwritten as soon as archiving is complete
135  this->mMeshChangesDuringSimulation = true;
136  Clear();
137 }
138 
139 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
141 {
142  Clear();
143 }
144 
145 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
147 {
148  return mCellRearrangementThreshold;
149 }
150 
151 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
153 {
154  return mT2Threshold;
155 }
156 
157 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
159 {
160  return mCellRearrangementRatio;
161 }
162 
163 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
165 {
166  return this->mProtorosetteFormationProbability;
167 }
168 
169 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
171 {
172  return this->mProtorosetteResolutionProbabilityPerTimestep;
173 }
174 
175 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
177 {
178  return this->mRosetteResolutionProbabilityPerTimestep;
179 }
180 
181 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
183 {
184  mDistanceForT3SwapChecking = distanceForT3SwapChecking;
185 }
186 
187 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
189 {
190  return mDistanceForT3SwapChecking;
191 }
192 
193 
194 
195 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
197 {
198  return mCheckForInternalIntersections;
199 }
200 
201 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
203 {
204  mCellRearrangementThreshold = cellRearrangementThreshold;
205 }
206 
207 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
209 {
210  mT2Threshold = t2Threshold;
211 }
212 
213 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
215 {
216  mCellRearrangementRatio = cellRearrangementRatio;
217 }
218 
219 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
221 {
222  // Check that the new value is in [0,1]
223  if (protorosetteFormationProbability < 0.0)
224  {
225  EXCEPTION("Attempting to assign a negative probability.");
226  }
227  if (protorosetteFormationProbability > 1.0)
228  {
229  EXCEPTION("Attempting to assign a probability greater than one.");
230  }
231 
232  // Assign the new value
233  mProtorosetteFormationProbability = protorosetteFormationProbability;
234 }
235 
236 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
238 {
239  // Check that the new value is in [0,1]
240  if (protorosetteResolutionProbabilityPerTimestep < 0.0)
241  {
242  EXCEPTION("Attempting to assign a negative probability.");
243  }
244  if (protorosetteResolutionProbabilityPerTimestep > 1.0)
245  {
246  EXCEPTION("Attempting to assign a probability greater than one.");
247  }
248 
249  // Assign the new value
250  mProtorosetteResolutionProbabilityPerTimestep = protorosetteResolutionProbabilityPerTimestep;
251 }
252 
253 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
255 {
256  // Check that the new value is in [0,1]
257  if (rosetteResolutionProbabilityPerTimestep < 0.0)
258  {
259  EXCEPTION("Attempting to assign a negative probability.");
260  }
261  if (rosetteResolutionProbabilityPerTimestep > 1.0)
262  {
263  EXCEPTION("Attempting to assign a probability greater than one.");
264  }
265 
266  // Assign the new value
267  mRosetteResolutionProbabilityPerTimestep = rosetteResolutionProbabilityPerTimestep;
268 }
269 
270 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
272 {
273  mCheckForInternalIntersections = checkForInternalIntersections;
274 }
275 
276 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
278 {
279  mDeletedNodeIndices.clear();
280  mDeletedElementIndices.clear();
281 
283 }
284 
285 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
287 {
288  return this->mNodes.size() - mDeletedNodeIndices.size();
289 }
290 
291 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
293 {
294  return this->mElements.size() - mDeletedElementIndices.size();
295 }
296 
297 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
298 std::vector< c_vector<double, SPACE_DIM> > MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLocationsOfT1Swaps()
299 {
300  return mLocationsOfT1Swaps;
301 }
302 
303 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
305 {
306  return mLastT2SwapLocation;
307 }
308 
309 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
310 std::vector< c_vector<double, SPACE_DIM> > MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLocationsOfT3Swaps()
311 {
312  return mLocationsOfT3Swaps;
313 }
314 
315 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
317 {
318  mLocationsOfT1Swaps.clear();
319 }
320 
321 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
323 {
324  mLocationsOfT3Swaps.clear();
325 }
326 
327 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
329 {
330  if (mDeletedNodeIndices.empty())
331  {
332  pNewNode->SetIndex(this->mNodes.size());
333  this->mNodes.push_back(pNewNode);
334  }
335  else
336  {
337  unsigned index = mDeletedNodeIndices.back();
338  pNewNode->SetIndex(index);
339  mDeletedNodeIndices.pop_back();
340  delete this->mNodes[index];
341  this->mNodes[index] = pNewNode;
342  }
343  return pNewNode->GetIndex();
344 }
345 
346 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
348 {
349  unsigned new_element_index = pNewElement->GetIndex();
350 
351  if (new_element_index == this->mElements.size())
352  {
353  this->mElements.push_back(pNewElement);
354  }
355  else
356  {
357  this->mElements[new_element_index] = pNewElement;
358  }
359  pNewElement->RegisterWithNodes();
360  return pNewElement->GetIndex();
361 }
362 
363 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
365 {
366  this->mNodes[nodeIndex]->SetPoint(point);
367 }
368 
369 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
371  c_vector<double, SPACE_DIM> axisOfDivision,
372  bool placeOriginalElementBelow)
373 {
374  assert(SPACE_DIM == 2); // LCOV_EXCL_LINE
375  assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
376 
377  // Get the centroid of the element
378  c_vector<double, SPACE_DIM> centroid = this->GetCentroidOfElement(pElement->GetIndex());
379 
380  // Create a vector perpendicular to the axis of division
381  c_vector<double, SPACE_DIM> perp_axis;
382  perp_axis(0) = -axisOfDivision(1);
383  perp_axis(1) = axisOfDivision(0);
384 
385  /*
386  * Find which edges the axis of division crosses by finding any node
387  * that lies on the opposite side of the axis of division to its next
388  * neighbour.
389  */
390  unsigned num_nodes = pElement->GetNumNodes();
391  std::vector<unsigned> intersecting_nodes;
392  bool is_current_node_on_left = (inner_prod(this->GetVectorFromAtoB(pElement->GetNodeLocation(0), centroid), perp_axis) >= 0);
393  for (unsigned i=0; i<num_nodes; i++)
394  {
395  bool is_next_node_on_left = (inner_prod(this->GetVectorFromAtoB(pElement->GetNodeLocation((i+1)%num_nodes), centroid), perp_axis) >= 0);
396  if (is_current_node_on_left != is_next_node_on_left)
397  {
398  intersecting_nodes.push_back(i);
399  }
400  is_current_node_on_left = is_next_node_on_left;
401  }
402 
403  // If the axis of division does not cross two edges then we cannot proceed
404  if (intersecting_nodes.size() != 2)
405  {
406  EXCEPTION("Cannot proceed with element division: the given axis of division does not cross two edges of the element");
407  }
408 
409  std::vector<unsigned> division_node_global_indices;
410  unsigned nodes_added = 0;
411 
412  // Find the intersections between the axis of division and the element edges
413  for (unsigned i=0; i<intersecting_nodes.size(); i++)
414  {
415  /*
416  * Get pointers to the nodes forming the edge into which one new node will be inserted.
417  *
418  * Note that when we use the first entry of intersecting_nodes to add a node,
419  * we change the local index of the second entry of intersecting_nodes in
420  * pElement, so must account for this by moving one entry further on.
421  */
422  Node<SPACE_DIM>* p_node_A = pElement->GetNode((intersecting_nodes[i]+nodes_added)%pElement->GetNumNodes());
423  Node<SPACE_DIM>* p_node_B = pElement->GetNode((intersecting_nodes[i]+nodes_added+1)%pElement->GetNumNodes());
424 
425  // Find the indices of the elements owned by each node on the edge into which one new node will be inserted
426  std::set<unsigned> elems_containing_node_A = p_node_A->rGetContainingElementIndices();
427  std::set<unsigned> elems_containing_node_B = p_node_B->rGetContainingElementIndices();
428 
429  c_vector<double, SPACE_DIM> position_a = p_node_A->rGetLocation();
430  c_vector<double, SPACE_DIM> position_b = p_node_B->rGetLocation();
431  c_vector<double, SPACE_DIM> a_to_b = this->GetVectorFromAtoB(position_a, position_b);
432 
433  c_vector<double, SPACE_DIM> intersection;
434 
435  if (norm_2(a_to_b) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
436  {
437  WARNING("Edge is too small for normal division; putting node in the middle of a and b. There may be T1 swaps straight away.");
439  intersection = position_a + 0.5*a_to_b;
440  }
441  else
442  {
443  // Find the location of the intersection
444  double determinant = a_to_b[0]*axisOfDivision[1] - a_to_b[1]*axisOfDivision[0];
445 
446  // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
447  c_vector<double, SPACE_DIM> moved_centroid;
448  moved_centroid = position_a + this->GetVectorFromAtoB(position_a, centroid);
449 
450  double alpha = (moved_centroid[0]*a_to_b[1] - position_a[0]*a_to_b[1]
451  -moved_centroid[1]*a_to_b[0] + position_a[1]*a_to_b[0])/determinant;
452 
453  intersection = moved_centroid + alpha*axisOfDivision;
454 
455  /*
456  * If then new node is too close to one of the edge nodes, then reposition it
457  * a distance mCellRearrangementRatio*mCellRearrangementThreshold further along the edge.
458  */
459  c_vector<double, SPACE_DIM> a_to_intersection = this->GetVectorFromAtoB(position_a, intersection);
460  if (norm_2(a_to_intersection) < mCellRearrangementThreshold)
461  {
462  intersection = position_a + mCellRearrangementRatio*mCellRearrangementThreshold*a_to_b/norm_2(a_to_b);
463  }
464 
465  c_vector<double, SPACE_DIM> b_to_intersection = this->GetVectorFromAtoB(position_b, intersection);
466  if (norm_2(b_to_intersection) < mCellRearrangementThreshold)
467  {
468  assert(norm_2(a_to_intersection) > mCellRearrangementThreshold); // to prevent moving intersection back to original position
469 
470  intersection = position_b - mCellRearrangementRatio*mCellRearrangementThreshold*a_to_b/norm_2(a_to_b);
471  }
472  }
473 
474  /*
475  * The new node is boundary node if the 2 nodes are boundary nodes and the elements don't look like
476  * ___A___
477  * | | |
478  * |___|___|
479  * B
480  */
481  bool is_boundary = false;
482  if (p_node_A->IsBoundaryNode() && p_node_B->IsBoundaryNode())
483  {
484  if (elems_containing_node_A.size() != 2 ||
485  elems_containing_node_B.size() != 2 ||
486  elems_containing_node_A != elems_containing_node_B)
487  {
488  is_boundary = true;
489  }
490  }
491 
492  // Add a new node to the mesh at the location of the intersection
493  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, is_boundary, intersection[0], intersection[1]));
494  nodes_added++;
495 
496  // Now make sure the new node is added to all neighbouring elements
497 
498  // Find common elements
499  std::set<unsigned> shared_elements;
500  std::set_intersection(elems_containing_node_A.begin(),
501  elems_containing_node_A.end(),
502  elems_containing_node_B.begin(),
503  elems_containing_node_B.end(),
504  std::inserter(shared_elements, shared_elements.begin()));
505 
506  // Iterate over common elements
507  unsigned node_A_index = p_node_A->GetIndex();
508  unsigned node_B_index = p_node_B->GetIndex();
509  for (std::set<unsigned>::iterator iter = shared_elements.begin();
510  iter != shared_elements.end();
511  ++iter)
512  {
513  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(*iter);
514 
515  // Find which node has the lower local index in this element
516  unsigned local_indexA = p_element->GetNodeLocalIndex(node_A_index);
517  unsigned local_indexB = p_element->GetNodeLocalIndex(node_B_index);
518 
519  unsigned index = local_indexB;
520 
521  // If node B has a higher index then use node A's index...
522  if (local_indexB > local_indexA)
523  {
524  index = local_indexA;
525 
526  // ...unless nodes A and B share the element's last edge
527  if ((local_indexA == 0) && (local_indexB == p_element->GetNumNodes()-1))
528  {
529  index = local_indexB;
530  }
531  }
532  else if ((local_indexB == 0) && (local_indexA == p_element->GetNumNodes()-1))
533  {
534  // ...otherwise use node B's index, unless nodes A and B share the element's last edge
535  index = local_indexA;
536  }
537 
538  // Add new node to this element
539  this->GetElement(*iter)->AddNode(this->GetNode(new_node_global_index), index);
540  }
541 
542  // Store index of new node
543  division_node_global_indices.push_back(new_node_global_index);
544  }
545 
546  // Now call DivideElement() to divide the element using the new nodes
547  unsigned new_element_index = DivideElement(pElement,
548  pElement->GetNodeLocalIndex(division_node_global_indices[0]),
549  pElement->GetNodeLocalIndex(division_node_global_indices[1]),
550  placeOriginalElementBelow);
551 
552  return new_element_index;
553 }
554 
555 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
557  bool placeOriginalElementBelow)
558 {
559  assert(SPACE_DIM == 2); // LCOV_EXCL_LINE
560  assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
561 
562  c_vector<double, SPACE_DIM> short_axis = this->GetShortAxisOfElement(pElement->GetIndex());
563 
564  unsigned new_element_index = DivideElementAlongGivenAxis(pElement, short_axis, placeOriginalElementBelow);
565  return new_element_index;
566 }
567 
568 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
570  unsigned nodeAIndex,
571  unsigned nodeBIndex,
572  bool placeOriginalElementBelow)
573 {
574  assert(SPACE_DIM == 2); // LCOV_EXCL_LINE
575  assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
576 
577  // Sort nodeA and nodeB such that nodeBIndex > nodeAindex
578  assert(nodeBIndex != nodeAIndex);
579  unsigned node1_index = (nodeAIndex < nodeBIndex) ? nodeAIndex : nodeBIndex; // low index
580  unsigned node2_index = (nodeAIndex < nodeBIndex) ? nodeBIndex : nodeAIndex; // high index
581 
582  // Store the number of nodes in the element (this changes when nodes are deleted from the element)
583  unsigned num_nodes = pElement->GetNumNodes();
584 
585  // Copy the nodes in this element
586  std::vector<Node<SPACE_DIM>*> nodes_elem;
587  for (unsigned i=0; i<num_nodes; i++)
588  {
589  nodes_elem.push_back(pElement->GetNode(i));
590  }
591 
592  // Get the index of the new element
593  unsigned new_element_index;
594  if (mDeletedElementIndices.empty())
595  {
596  new_element_index = this->mElements.size();
597  }
598  else
599  {
600  new_element_index = mDeletedElementIndices.back();
601  mDeletedElementIndices.pop_back();
602  delete this->mElements[new_element_index];
603  }
604 
605  // Add the new element to the mesh
606  AddElement(new VertexElement<ELEMENT_DIM,SPACE_DIM>(new_element_index, nodes_elem));
607 
614  // Find lowest element
616  double height_midpoint_1 = 0.0;
617  double height_midpoint_2 = 0.0;
618  unsigned counter_1 = 0;
619  unsigned counter_2 = 0;
620 
621  for (unsigned i=0; i<num_nodes; i++)
622  {
623  if (i>=node1_index && i<=node2_index)
624  {
625  height_midpoint_1 += pElement->GetNode(i)->rGetLocation()[1];
626  counter_1++;
627  }
628  if (i<=node1_index || i>=node2_index)
629  {
630  height_midpoint_2 += pElement->GetNode(i)->rGetLocation()[1];
631  counter_2++;
632  }
633  }
634  height_midpoint_1 /= (double)counter_1;
635  height_midpoint_2 /= (double)counter_2;
636 
637  for (unsigned i=num_nodes; i>0; i--)
638  {
639  if (i-1 < node1_index || i-1 > node2_index)
640  {
641  if (height_midpoint_1 < height_midpoint_2)
642  {
643  if (placeOriginalElementBelow)
644  {
645  pElement->DeleteNode(i-1);
646  }
647  else
648  {
649  this->mElements[new_element_index]->DeleteNode(i-1);
650  }
651  }
652  else
653  {
654  if (placeOriginalElementBelow)
655  {
656  this->mElements[new_element_index]->DeleteNode(i-1);
657  }
658  else
659  {
660  pElement->DeleteNode(i-1);
661  }
662  }
663  }
664  else if (i-1 > node1_index && i-1 < node2_index)
665  {
666  if (height_midpoint_1 < height_midpoint_2)
667  {
668  if (placeOriginalElementBelow)
669  {
670  this->mElements[new_element_index]->DeleteNode(i-1);
671  }
672  else
673  {
674  pElement->DeleteNode(i-1);
675  }
676  }
677  else
678  {
679  if (placeOriginalElementBelow)
680  {
681  pElement->DeleteNode(i-1);
682  }
683  else
684  {
685  this->mElements[new_element_index]->DeleteNode(i-1);
686  }
687  }
688  }
689  }
690 
691  return new_element_index;
692 }
693 
694 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
696 {
697  assert(SPACE_DIM == 2); // LCOV_EXCL_LINE
698 
699  // Mark any nodes that are contained only in this element as deleted
700  for (unsigned i=0; i<this->mElements[index]->GetNumNodes(); i++)
701  {
702  Node<SPACE_DIM>* p_node = this->mElements[index]->GetNode(i);
703 
704  if (p_node->rGetContainingElementIndices().size() == 1)
705  {
706  DeleteNodePriorToReMesh(p_node->GetIndex());
707  }
708 
709  // Mark all the nodes contained in the removed element as boundary nodes
710  p_node->SetAsBoundaryNode(true);
711  }
712 
713  // Mark this element as deleted
714  this->mElements[index]->MarkAsDeleted();
715  mDeletedElementIndices.push_back(index);
716 }
717 
718 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
720 {
721  this->mNodes[index]->MarkAsDeleted();
722  mDeletedNodeIndices.push_back(index);
723 }
724 
725 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
727 {
728  // Find the indices of the elements owned by each node
729  std::set<unsigned> elements_containing_nodeA = pNodeA->rGetContainingElementIndices();
730  std::set<unsigned> elements_containing_nodeB = pNodeB->rGetContainingElementIndices();
731 
732  // Find common elements
733  std::set<unsigned> shared_elements;
734  std::set_intersection(elements_containing_nodeA.begin(),
735  elements_containing_nodeA.end(),
736  elements_containing_nodeB.begin(),
737  elements_containing_nodeB.end(),
738  std::inserter(shared_elements, shared_elements.begin()));
739 
740  // Check that the nodes have a common edge and not more than 2
741  assert(!shared_elements.empty());
742  assert(shared_elements.size()<=2u);
743 
744  // Specify if it's a boundary node
745  bool is_boundary_node = false;
746  if (shared_elements.size()==1u)
747  {
748  // If only one shared element then must be on the boundary.
749  assert((pNodeA->IsBoundaryNode()) && (pNodeB->IsBoundaryNode()));
750  is_boundary_node = true;
751  }
752 
753  // Create a new node (position is not important as it will be changed)
754  Node<SPACE_DIM>* p_new_node = new Node<SPACE_DIM>(GetNumNodes(), is_boundary_node, 0.0, 0.0);
755 
756  // Update the node location
757  c_vector<double, SPACE_DIM> new_node_position = pNodeA->rGetLocation() + 0.5*this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation());
758  ChastePoint<SPACE_DIM> point(new_node_position);
759  p_new_node->SetPoint(new_node_position);
760 
761  // Add node to mesh
762  this->mNodes.push_back(p_new_node);
763 
764  // Iterate over common elements
765  unsigned node_A_index = pNodeA->GetIndex();
766  unsigned node_B_index = pNodeB->GetIndex();
767  for (std::set<unsigned>::iterator iter = shared_elements.begin();
768  iter != shared_elements.end();
769  ++iter)
770  {
771  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(*iter);
772 
773  // Find which node has the lower local index in this element
774  unsigned local_indexA = p_element->GetNodeLocalIndex(node_A_index);
775  unsigned local_indexB = p_element->GetNodeLocalIndex(node_B_index);
776 
777  unsigned index = local_indexB;
778 
779  // If node B has a higher index then use node A's index...
780  if (local_indexB > local_indexA)
781  {
782  index = local_indexA;
783 
784  // ...unless nodes A and B share the element's last edge
785  if ((local_indexA == 0) && (local_indexB == p_element->GetNumNodes()-1))
786  {
787  index = local_indexB;
788  }
789  }
790  else if ((local_indexB == 0) && (local_indexA == p_element->GetNumNodes()-1))
791  {
792  // ...otherwise use node B's index, unless nodes A and B share the element's last edge
793  index = local_indexA;
794  }
795 
796  // Add new node to this element
797  this->GetElement(*iter)->AddNode(p_new_node, index);
798  }
799 }
800 
801 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
803 {
804  // Make sure the map is big enough. Each entry will be set in the loop below.
805  rElementMap.Resize(this->GetNumAllElements());
806 
807  // Remove any elements that have been marked for deletion and store all other elements in a temporary structure
808  std::vector<VertexElement<ELEMENT_DIM, SPACE_DIM>*> live_elements;
809  for (unsigned i=0; i<this->mElements.size(); i++)
810  {
811  if (this->mElements[i]->IsDeleted())
812  {
813  delete this->mElements[i];
814  rElementMap.SetDeleted(i);
815  }
816  else
817  {
818  live_elements.push_back(this->mElements[i]);
819  rElementMap.SetNewIndex(i, (unsigned)(live_elements.size()-1));
820  }
821  }
822 
823  // Sanity check
824  assert(mDeletedElementIndices.size() == this->mElements.size() - live_elements.size());
825 
826  // Repopulate the elements vector and reset the list of deleted element indices
827  mDeletedElementIndices.clear();
828  this->mElements = live_elements;
829 
830  // Finally, reset the element indices to run from zero
831  for (unsigned i=0; i<this->mElements.size(); i++)
832  {
833  this->mElements[i]->ResetIndex(i);
834  }
835 
836  // Remove deleted nodes
837  RemoveDeletedNodes();
838 }
839 
840 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
842 {
843  // Remove any nodes that have been marked for deletion and store all other nodes in a temporary structure
844  std::vector<Node<SPACE_DIM>*> live_nodes;
845  for (unsigned i=0; i<this->mNodes.size(); i++)
846  {
847  if (this->mNodes[i]->IsDeleted())
848  {
849  delete this->mNodes[i];
850  }
851  else
852  {
853  live_nodes.push_back(this->mNodes[i]);
854  }
855  }
856 
857  // Sanity check
858  assert(mDeletedNodeIndices.size() == this->mNodes.size() - live_nodes.size());
859 
860  // Repopulate the nodes vector and reset the list of deleted node indices
861  this->mNodes = live_nodes;
862  mDeletedNodeIndices.clear();
863 
864  // Finally, reset the node indices to run from zero
865  for (unsigned i=0; i<this->mNodes.size(); i++)
866  {
867  this->mNodes[i]->SetIndex(i);
868  }
869 }
870 
871 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
873 {
874  // Make sure that we are in the correct dimension - this code will be eliminated at compile time
875  assert(SPACE_DIM==2 || SPACE_DIM==3); // LCOV_EXCL_LINE
876  assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
877 
878 
879  if (SPACE_DIM == 2)
880  {
881  // Make sure the map is big enough
882  rElementMap.Resize(this->GetNumAllElements());
883 
884  /*
885  * To begin the remeshing process, we do not need to call Clear() and remove all current data,
886  * since cell birth, rearrangement and death result only in local remeshing of a vertex-based
887  * mesh. Instead, we just remove any deleted elements and nodes.
888  */
889  RemoveDeletedNodesAndElements(rElementMap);
890  bool recheck_mesh = true;
891  while (recheck_mesh == true)
892  {
893  // We check for any short edges and perform swaps if necessary and possible.
894  recheck_mesh = CheckForSwapsFromShortEdges();
895  }
896 
897  // Check for element intersections
898  recheck_mesh = true;
899  while (recheck_mesh == true)
900  {
901  // Check mesh for intersections, and perform T3 swaps where required
902  recheck_mesh = CheckForIntersections();
903  }
904 
905  RemoveDeletedNodes();
906 
907  /*
908  * This is handled in a separate method to allow child classes to implement additional ReMeshing functionality
909  * (see #2664).
910  */
911  this->CheckForRosettes();
912  }
913  else // 3D
914  {
915 // LCOV_EXCL_START
916  EXCEPTION("Remeshing has not been implemented in 3D (see #827 and #860)\n");
917 // LCOV_EXCL_STOP
919  }
920 }
921 
922 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
924 {
925  VertexElementMap map(GetNumElements());
926  ReMesh(map);
927 }
928 
929 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
931 {
932  // Loop over elements to check for T1 swaps
933  for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
934  elem_iter != this->GetElementIteratorEnd();
935  ++elem_iter)
936  {
938 
939  unsigned num_nodes = elem_iter->GetNumNodes();
940  assert(num_nodes > 0);
941 
942  // Loop over the nodes contained in this element
943  for (unsigned local_index=0; local_index<num_nodes; local_index++)
944  {
945  // Find locations of the current node and anticlockwise node
946  Node<SPACE_DIM>* p_current_node = elem_iter->GetNode(local_index);
947  unsigned local_index_plus_one = (local_index+1)%num_nodes;
948  Node<SPACE_DIM>* p_anticlockwise_node = elem_iter->GetNode(local_index_plus_one);
949 
950  // Find distance between nodes
951  double distance_between_nodes = this->GetDistanceBetweenNodes(p_current_node->GetIndex(), p_anticlockwise_node->GetIndex());
952 
953  // If the nodes are too close together...
954  if (distance_between_nodes < mCellRearrangementThreshold)
955  {
956  // ...then check if any triangular elements are shared by these nodes...
957  std::set<unsigned> elements_of_node_a = p_current_node->rGetContainingElementIndices();
958  std::set<unsigned> elements_of_node_b = p_anticlockwise_node->rGetContainingElementIndices();
959 
960  std::set<unsigned> shared_elements;
961  std::set_intersection(elements_of_node_a.begin(), elements_of_node_a.end(),
962  elements_of_node_b.begin(), elements_of_node_b.end(),
963  std::inserter(shared_elements, shared_elements.begin()));
964 
965  bool both_nodes_share_triangular_element = false;
966  for (std::set<unsigned>::const_iterator it = shared_elements.begin();
967  it != shared_elements.end();
968  ++it)
969  {
970  if (this->GetElement(*it)->GetNumNodes() <= 3)
971  {
972  both_nodes_share_triangular_element = true;
973  break;
974  }
975  }
976 
977  // ...and if none are, then perform the required type of swap and halt the search, returning true
978  if (!both_nodes_share_triangular_element)
979  {
980  IdentifySwapType(p_current_node, p_anticlockwise_node);
981  return true;
982  }
983  }
984  }
985  }
986 
987  return false;
988 }
989 
990 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
992 {
993  // Loop over elements to check for T2 swaps
994  for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
995  elem_iter != this->GetElementIteratorEnd();
996  ++elem_iter)
997  {
998  // If this element is triangular...
999  if (elem_iter->GetNumNodes() == 3)
1000  {
1001  // ...and smaller than the threshold area...
1002  if (this->GetVolumeOfElement(elem_iter->GetIndex()) < GetT2Threshold())
1003  {
1004  // ...then perform a T2 swap and break out of the loop
1005  PerformT2Swap(*elem_iter);
1007  rElementMap.SetDeleted(elem_iter->GetIndex());
1008  return true;
1009  }
1010  }
1011  }
1012  return false;
1013 }
1014 
1015 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1017 {
1018  // If checking for internal intersections as well as on the boundary, then check that no nodes have overlapped any elements...
1019  if (mCheckForInternalIntersections)
1020  {
1022  for (typename AbstractMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator node_iter = this->GetNodeIteratorBegin();
1023  node_iter != this->GetNodeIteratorEnd();
1024  ++node_iter)
1025  {
1026  assert(!(node_iter->IsDeleted()));
1027 
1028  for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
1029  elem_iter != this->GetElementIteratorEnd();
1030  ++elem_iter)
1031  {
1032  unsigned elem_index = elem_iter->GetIndex();
1033 
1034  // Check that the node is not part of this element
1035  if (node_iter->rGetContainingElementIndices().count(elem_index) == 0)
1036  {
1037  if (this->ElementIncludesPoint(node_iter->rGetLocation(), elem_index))
1038  {
1039  PerformIntersectionSwap(&(*node_iter), elem_index);
1040  return true;
1041  }
1042  }
1043  }
1044  }
1045  }
1046  else
1047  {
1048  // ...otherwise, just check that no boundary nodes have overlapped any boundary elements
1049  // First: find all boundary element and calculate their centroid only once
1050  std::vector<unsigned> boundary_element_indices;
1051  std::vector< c_vector<double, SPACE_DIM> > boundary_element_centroids;
1052  for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
1053  elem_iter != this->GetElementIteratorEnd();
1054  ++elem_iter)
1055  {
1056  if (elem_iter->IsElementOnBoundary())
1057  {
1058  unsigned element_index = elem_iter->GetIndex();
1059  boundary_element_indices.push_back(element_index);
1060  // should be a map but I am too lazy to look up the syntax
1061  boundary_element_centroids.push_back(this->GetCentroidOfElement(element_index));
1062  }
1063  }
1064 
1065  // Second: Check intersections only for those nodes and elements within
1066  // mDistanceForT3SwapChecking within each other (node<-->element centroid)
1067  for (typename AbstractMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator node_iter = this->GetNodeIteratorBegin();
1068  node_iter != this->GetNodeIteratorEnd();
1069  ++node_iter)
1070  {
1071  if (node_iter->IsBoundaryNode())
1072  {
1073  assert(!(node_iter->IsDeleted()));
1074 
1075  // index in boundary_element_centroids and boundary_element_indices
1076  unsigned boundary_element_index = 0;
1077  for (std::vector<unsigned>::iterator elem_iter = boundary_element_indices.begin();
1078  elem_iter != boundary_element_indices.end();
1079  ++elem_iter)
1080  {
1081  // Check that the node is not part of this element
1082  if (node_iter->rGetContainingElementIndices().count(*elem_iter) == 0)
1083  {
1084  c_vector<double, SPACE_DIM> node_location = node_iter->rGetLocation();
1085  c_vector<double, SPACE_DIM> element_centroid = boundary_element_centroids[boundary_element_index];
1086  double node_element_distance = norm_2(this->GetVectorFromAtoB(node_location, element_centroid));
1087 
1088  if ( node_element_distance < mDistanceForT3SwapChecking )
1089  {
1090  if (this->ElementIncludesPoint(node_iter->rGetLocation(), *elem_iter))
1091  {
1092  this->PerformT3Swap(&(*node_iter), *elem_iter);
1093  return true;
1094  }
1095  }
1096  }
1097  // increment the boundary element index
1098  boundary_element_index +=1u;
1099  }
1100  }
1101  }
1102 
1103  }
1104  return false;
1105 }
1106 
1107 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1109 {
1110  // Find the sets of elements containing nodes A and B
1111  std::set<unsigned> nodeA_elem_indices = pNodeA->rGetContainingElementIndices();
1112  std::set<unsigned> nodeB_elem_indices = pNodeB->rGetContainingElementIndices();
1113 
1114  // Form the set union
1115  std::set<unsigned> all_indices, temp_union_set;
1116  std::set_union(nodeA_elem_indices.begin(), nodeA_elem_indices.end(),
1117  nodeB_elem_indices.begin(), nodeB_elem_indices.end(),
1118  std::inserter(temp_union_set, temp_union_set.begin()));
1119  all_indices.swap(temp_union_set); // temp_set will be deleted, all_indices now contains all the indices of elements
1120  // that touch the potentially swapping nodes
1121 
1122  if ((nodeA_elem_indices.size()>3) || (nodeB_elem_indices.size()>3))
1123  {
1124  /*
1125  * Looks like
1126  *
1127  * \
1128  * \ A B
1129  * ---o---o---
1130  * /
1131  * /
1132  *
1133  */
1134 
1135  /*
1136  * This case is handled in a separate method to allow child classes to implement different
1137  * functionality for high-order-junction remodelling events (see #2664).
1138  */
1139  this->HandleHighOrderJunctions(pNodeA, pNodeB);
1140  }
1141  else // each node is contained in at most three elements
1142  {
1143  switch (all_indices.size())
1144  {
1145  case 1:
1146  {
1147  /*
1148  * Each node is contained in a single element, so the nodes must lie on the boundary
1149  * of the mesh, as shown below. In this case, we merge the nodes and tidy up node
1150  * indices through calls to PerformNodeMerge() and RemoveDeletedNodes().
1151  *
1152  * A B
1153  * ---o---o---
1154  */
1155  assert(pNodeA->IsBoundaryNode());
1156  assert(pNodeB->IsBoundaryNode());
1157 
1158  PerformNodeMerge(pNodeA, pNodeB);
1159  RemoveDeletedNodes();
1160  break;
1161  }
1162  case 2:
1163  {
1164  if (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==2)
1165  {
1166  if (pNodeA->IsBoundaryNode() && pNodeB->IsBoundaryNode())
1167  {
1168  /*
1169  * The node configuration is as shown below, with voids on either side. In this case
1170  * we perform a T1 swap, which separates the elements.
1171  *
1172  * \ /
1173  * \ / Node A
1174  * (1) | (2) (element number in brackets)
1175  * / \ Node B
1176  * / \
1177  */
1178  PerformT1Swap(pNodeA, pNodeB,all_indices);
1179  }
1180  else if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
1181  {
1182  /*
1183  * The node configuration is as shown below, with a void on one side. We should not
1184  * be able to reach this case at present, since we allow only for three-way junctions
1185  * or boundaries, so we throw an exception.
1186  *
1187  * \ /
1188  * \ / Node A
1189  * (1) | (2) (element number in brackets)
1190  * x Node B
1191  * |
1192  */
1193  EXCEPTION("There is a non-boundary node contained only in two elements; something has gone wrong.");
1194  }
1195  else
1196  {
1197  /*
1198  * Each node is contained in two elements, so the nodes lie on an internal edge, as shown below.
1199  * We should not be able to reach this case at present, since we allow only for three-way junctions
1200  * or boundaries, so we throw an exception.
1201  *
1202  * A B
1203  * ---o---o---
1204  */
1205  EXCEPTION("There are non-boundary nodes contained only in two elements; something has gone wrong.");
1206  }
1207  }// from [if (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==2)]
1208  else
1209  {
1210  /*
1211  * The node configuration looks like that shown below. In this case, we merge the nodes
1212  * and tidy up node indices through calls to PerformNodeMerge() and RemoveDeletedNodes().
1213  *
1214  * Outside
1215  * /
1216  * --o--o (2)
1217  * (1) \
1218  *
1219  * ///\todo this should be a T1 swap (see #1263 and #2401)
1220  * Referring to the todo: this should probably stay a node-merge. If this is a T1 swap then
1221  * the single boundary node will travel from element 1 to element 2, but still remain a single node.
1222  * I.e. we would not reduce the total number of nodes in this situation.
1223  */
1224  PerformNodeMerge(pNodeA, pNodeB);
1225  RemoveDeletedNodes();
1226  }
1227  break;
1228  }
1229  case 3:
1230  {
1231  if (nodeA_elem_indices.size()==1 || nodeB_elem_indices.size()==1)
1232  {
1233  /*
1234  * One node is contained in one element and the other node is contained in three elements.
1235  * We should not be able to reach this case at present, since we allow each boundary node
1236  * to be contained in at most two elements, so we throw an exception.
1237  *
1238  * A B
1239  *
1240  * empty /
1241  * / (3)
1242  * ---o---o----- (element number in brackets)
1243  * (1) \ (2)
1244  * \
1245  */
1246  assert(pNodeA->IsBoundaryNode());
1247  assert(pNodeB->IsBoundaryNode());
1248 
1249  EXCEPTION("There is a boundary node contained in three elements something has gone wrong.");
1250  }
1251  else if (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==2)
1252  {
1253  // The short edge must be at the boundary. We need to check whether this edge is
1254  // adjacent to a triangular void before we swap. If it is a triangular void, we perform a T2-type swap.
1255  // If not, then we perform a normal T1 swap. I.e. in detail we need to check whether the
1256  // element in nodeA_elem_indices which is not in nodeB_elem_indices contains a shared node
1257  // with the element in nodeB_elem_indices which is not in nodeA_elem_indices.
1258 
1259  std::set<unsigned> element_A_not_B, temp_set;
1260  std::set_difference(all_indices.begin(), all_indices.end(), nodeB_elem_indices.begin(),
1261  nodeB_elem_indices.end(), std::inserter(temp_set, temp_set.begin()));
1262  element_A_not_B.swap(temp_set);
1263 
1264  // There must be only one such element
1265  assert(element_A_not_B.size() == 1);
1266 
1267  std::set<unsigned> element_B_not_A;
1268  std::set_difference(all_indices.begin(), all_indices.end(), nodeA_elem_indices.begin(),
1269  nodeA_elem_indices.end(), std::inserter(temp_set, temp_set.begin()));
1270  element_B_not_A.swap(temp_set);
1271 
1272  // There must be only one such element
1273  assert(element_B_not_A.size() == 1);
1274 
1275  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_A_not_B = this->mElements[*element_A_not_B.begin()];
1276  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_B_not_A = this->mElements[*element_B_not_A.begin()];
1277 
1278  unsigned local_index_1 = p_element_A_not_B->GetNodeLocalIndex(pNodeA->GetIndex());
1279  unsigned next_node_1 = p_element_A_not_B->GetNodeGlobalIndex((local_index_1 + 1)%(p_element_A_not_B->GetNumNodes()));
1280  unsigned previous_node_1 = p_element_A_not_B->GetNodeGlobalIndex(
1281  (local_index_1 + p_element_A_not_B->GetNumNodes() - 1)%(p_element_A_not_B->GetNumNodes()));
1282  unsigned local_index_2 = p_element_B_not_A->GetNodeLocalIndex(pNodeB->GetIndex());
1283  unsigned next_node_2 = p_element_B_not_A->GetNodeGlobalIndex(
1284  (local_index_2 + 1)%(p_element_B_not_A->GetNumNodes()));
1285  unsigned previous_node_2 = p_element_B_not_A->GetNodeGlobalIndex(
1286  (local_index_2 + p_element_B_not_A->GetNumNodes() - 1)%(p_element_B_not_A->GetNumNodes()));
1287 
1288  if (next_node_1 == previous_node_2 || next_node_2 == previous_node_1)
1289  {
1290  /*
1291  * The node configuration looks like that shown below, and both nodes must be on the boundary.
1292  * In this case we remove the void through a call to PerformVoidRemoval().
1293  *
1294  * A C B A B
1295  * /\ \ /
1296  * /v \ \ (1) /
1297  * (3)o----o (1) or (2) o----o (3) (element number in brackets, v is a void)
1298  * / (2) \ \v /
1299  * / \ \/
1300  * C
1301  */
1302  assert(pNodeA->IsBoundaryNode());
1303  assert(pNodeB->IsBoundaryNode());
1304 
1305  // Get the third node in the triangular void
1306 
1307  unsigned nodeC_index;
1308  if (next_node_1 == previous_node_2 && next_node_2 != previous_node_1)
1309  {
1310  nodeC_index = next_node_1;
1311  }
1312  else if (next_node_2 == previous_node_1 && next_node_1 != previous_node_2)
1313  {
1314  nodeC_index = next_node_2;
1315  }
1316  else
1317  {
1318  assert(next_node_1 == previous_node_2 && next_node_2 == previous_node_1);
1326  EXCEPTION("Triangular element next to triangular void, not implemented yet.");
1327  }
1328 
1329  if (p_element_A_not_B->GetNumNodes() == 3u || p_element_B_not_A->GetNumNodes() == 3u)
1330  {
1338  EXCEPTION("Triangular element next to triangular void, not implemented yet.");
1339  }
1340 
1341  PerformVoidRemoval(pNodeA, pNodeB, this->mNodes[nodeC_index]);
1342  }
1343  else
1344  {
1345  /*
1346  * The node configuration looks like that below, and both nodes must lie on the boundary.
1347  * In this case we perform a T1 swap.
1348  *
1349  * A B A B
1350  * \ empty/ \ /
1351  * \ / \(1) /
1352  * (3) o--o (1) or (2) o--o (3) (element number in brackets)
1353  * / (2)\ / \
1354  * / \ /empty \
1355  */
1356  assert(pNodeA->IsBoundaryNode());
1357  assert(pNodeB->IsBoundaryNode());
1358 
1359  PerformT1Swap(pNodeA, pNodeB, all_indices);
1360  }
1361  } // from else if (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==2)
1362  else
1363  {
1364  // In this case, one node must be contained in two elements and the other in three elements.
1365  assert ( (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==3)
1366  || (nodeA_elem_indices.size()==3 && nodeB_elem_indices.size()==2) );
1367 
1368  // They can't both be boundary nodes
1369  assert(!(pNodeA->IsBoundaryNode() && pNodeB->IsBoundaryNode()));
1370 
1371  if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
1372  {
1373  /*
1374  * The node configuration looks like that shown below. We perform a T1 swap in this case.
1375  *
1376  * A B A B
1377  * \ / \ /
1378  * \ (1)/ \(1) /
1379  * (3) o--o (empty) or (empty) o--o (3) (element number in brackets)
1380  * / (2)\ /(2) \
1381  * / \ / \
1382  */
1383  PerformT1Swap(pNodeA, pNodeB, all_indices);
1384  }
1385  else
1386  {
1387  /*
1388  * The node configuration looks like that shown below. We should not be able to reach this case
1389  * at present, since we allow only for three-way junctions or boundaries, so we throw an exception.
1390  *
1391  * A B A B
1392  * \ /
1393  * \ (1) (1) /
1394  * (3) o--o--- or ---o--o (3) (element number in brackets)
1395  * / (2) (2) \
1396  * / \
1397  */
1398  EXCEPTION("There are non-boundary nodes contained only in two elements; something has gone wrong.");
1399  }
1400  }
1401  break;
1402  }
1403  case 4:
1404  {
1405  /*
1406  * The node configuration looks like that shown below. We perform a T1 swap in this case.
1407  *
1408  * \(1)/
1409  * \ / Node A
1410  * (2) | (4) (element number in brackets)
1411  * / \ Node B
1412  * /(3)\
1413  */
1414 
1415  /*
1416  * This case is handled in a separate method to allow child classes to implement different
1417  * functionality for junction remodelling events (see #2664).
1418  */
1419  if (mProtorosetteFormationProbability > RandomNumberGenerator::Instance()->ranf())
1420  {
1421  this->PerformNodeMerge(pNodeA, pNodeB);
1422  this->RemoveDeletedNodes();
1423  }
1424  else
1425  {
1426  this->PerformT1Swap(pNodeA, pNodeB, all_indices);
1427  }
1428  break;
1429  }
1430  default:
1431  // This can't happen
1432  NEVER_REACHED;
1433  }
1434  }
1435 }
1436 
1437 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1439 {
1440  // Find the sets of elements containing each of the nodes, sorted by index
1441  std::set<unsigned> nodeA_elem_indices = pNodeA->rGetContainingElementIndices();
1442  std::set<unsigned> nodeB_elem_indices = pNodeB->rGetContainingElementIndices();
1443 
1444  // Move node A to the mid-point
1445  pNodeA->rGetModifiableLocation() += 0.5 * this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation());
1446 
1447  // Update the elements previously containing node B to contain node A
1448  unsigned node_B_index = pNodeB->GetIndex();
1449  for (std::set<unsigned>::const_iterator it = nodeB_elem_indices.begin(); it != nodeB_elem_indices.end(); ++it)
1450  {
1451  // Find the local index of node B in this element
1452  unsigned node_B_local_index = this->mElements[*it]->GetNodeLocalIndex(node_B_index);
1453  assert(node_B_local_index < UINT_MAX); // this element contains node B
1454 
1455  /*
1456  * If this element already contains node A, then just remove node B.
1457  * Otherwise replace it with node A in the element and remove it from mNodes.
1458  */
1459  if (nodeA_elem_indices.count(*it) != 0)
1460  {
1461  this->mElements[*it]->DeleteNode(node_B_local_index);
1462  }
1463  else
1464  {
1465  // Replace node B with node A in this element
1466  this->mElements[*it]->UpdateNode(node_B_local_index, pNodeA);
1467  }
1468  }
1469 
1470  assert(!(this->mNodes[node_B_index]->IsDeleted()));
1471  this->mNodes[node_B_index]->MarkAsDeleted();
1472  mDeletedNodeIndices.push_back(node_B_index);
1473 }
1474 
1475 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1477  Node<SPACE_DIM>* pNodeB,
1478  std::set<unsigned>& rElementsContainingNodes)
1479 {
1480  // First compute and store the location of the T1 swap, which is at the midpoint of nodes A and B
1481  double distance_between_nodes_CD = mCellRearrangementRatio*mCellRearrangementThreshold;
1482 
1483  c_vector<double, SPACE_DIM> nodeA_location = pNodeA->rGetLocation();
1484  c_vector<double, SPACE_DIM> nodeB_location = pNodeB->rGetLocation();
1485  c_vector<double, SPACE_DIM> vector_AB = this->GetVectorFromAtoB(nodeA_location, nodeB_location);
1486  mLocationsOfT1Swaps.push_back(nodeA_location + 0.5*vector_AB);
1487 
1488  double distance_AB = norm_2(vector_AB);
1489  if (distance_AB < 1e-10)
1490  {
1491  EXCEPTION("Nodes are too close together, this shouldn't happen");
1492  }
1493 
1494  /*
1495  * Compute the locations of two new nodes C, D, placed on either side of the
1496  * edge E_old formed by nodes A and B, such that the edge E_new formed by the
1497  * new nodes is the perpendicular bisector of E_old, with |E_new| 'just larger'
1498  * (mCellRearrangementRatio) than mThresholdDistance.
1499  *
1500  * We implement the following changes to the mesh:
1501  *
1502  * The element whose index was in nodeA_elem_indices but not nodeB_elem_indices,
1503  * and the element whose index was in nodeB_elem_indices but not nodeA_elem_indices,
1504  * should now both contain nodes A and B.
1505  *
1506  * The element whose index was in nodeA_elem_indices and nodeB_elem_indices, and which
1507  * node C lies inside, should now only contain node A.
1508  *
1509  * The element whose index was in nodeA_elem_indices and nodeB_elem_indices, and which
1510  * node D lies inside, should now only contain node B.
1511  *
1512  * Iterate over all elements involved and identify which element they are
1513  * in the diagram then update the nodes as necessary.
1514  *
1515  * \(1)/
1516  * \ / Node A
1517  * (2) | (4) elements in brackets
1518  * / \ Node B
1519  * /(3)\
1520  */
1521 
1522  // Move nodes A and B to C and D respectively
1523  c_vector<double, SPACE_DIM> vector_CD;
1524  vector_CD(0) = -vector_AB(1) * distance_between_nodes_CD / distance_AB;
1525  vector_CD(1) = vector_AB(0) * distance_between_nodes_CD / distance_AB;
1526 
1527  c_vector<double, SPACE_DIM> nodeC_location = nodeA_location + 0.5*vector_AB - 0.5*vector_CD;
1528  c_vector<double, SPACE_DIM> nodeD_location = nodeC_location + vector_CD;
1529 
1530  pNodeA->rGetModifiableLocation() = nodeC_location;
1531  pNodeB->rGetModifiableLocation() = nodeD_location;
1532 
1533  // Find the sets of elements containing nodes A and B
1534  std::set<unsigned> nodeA_elem_indices = pNodeA->rGetContainingElementIndices();
1535  std::set<unsigned> nodeB_elem_indices = pNodeB->rGetContainingElementIndices();
1536 
1537  for (std::set<unsigned>::const_iterator it = rElementsContainingNodes.begin();
1538  it != rElementsContainingNodes.end();
1539  ++it)
1540  {
1541  // If, as in element 3 above, this element does not contain node A (now C)...
1542  if (nodeA_elem_indices.find(*it) == nodeA_elem_indices.end())
1543  {
1544  // ...then add it to the element just after node B (now D), going anticlockwise
1545  unsigned nodeB_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeB->GetIndex());
1546  assert(nodeB_local_index < UINT_MAX);
1547 
1548  this->mElements[*it]->AddNode(pNodeA, nodeB_local_index);
1549  }
1550  else if (nodeB_elem_indices.find(*it) == nodeB_elem_indices.end())
1551  {
1552  // Do similarly if the element does not contain node B (now D), as in element 4 above
1553  unsigned nodeA_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeA->GetIndex());
1554  assert(nodeA_local_index < UINT_MAX);
1555 
1556  this->mElements[*it]->AddNode(pNodeB, nodeA_local_index);
1557  }
1558  else
1559  {
1560  // If the element contains both nodes A and B (now C and D respectively)...
1561  unsigned nodeA_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeA->GetIndex());
1562  unsigned nodeB_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeB->GetIndex());
1563 
1564  assert(nodeA_local_index < UINT_MAX);
1565  assert(nodeB_local_index < UINT_MAX);
1566 
1567  /*
1568  * Locate local index of nodeA and nodeB and use the ordering to
1569  * identify the element, if nodeB_index > nodeA_index then element 4
1570  * and if nodeA_index > nodeB_index then element 2
1571  */
1572  unsigned nodeB_local_index_plus_one = (nodeB_local_index + 1)%(this->mElements[*it]->GetNumNodes());
1573 
1574  if (nodeA_local_index == nodeB_local_index_plus_one)
1575  {
1576  /*
1577  * In this case the local index of nodeA is the local index of
1578  * nodeB plus one so we are in element 2 so we remove nodeB
1579  */
1580  this->mElements[*it]->DeleteNode(nodeB_local_index);
1581  }
1582  else
1583  {
1584  assert(nodeB_local_index == (nodeA_local_index + 1)%(this->mElements[*it]->GetNumNodes())); // as A and B are next to each other
1585  /*
1586  * In this case the local index of nodeA is the local index of
1587  * nodeB minus one so we are in element 4 so we remove nodeA
1588  */
1589  this->mElements[*it]->DeleteNode(nodeA_local_index);
1590  }
1591  }
1592  }
1593 
1594  // Sort out boundary nodes
1595  if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
1596  {
1597  if (pNodeA->GetNumContainingElements() == 3)
1598  {
1599  pNodeA->SetAsBoundaryNode(false);
1600  }
1601  else
1602  {
1603  pNodeA->SetAsBoundaryNode(true);
1604  }
1605  if (pNodeB->GetNumContainingElements() == 3)
1606  {
1607  pNodeB->SetAsBoundaryNode(false);
1608  }
1609  else
1610  {
1611  pNodeB->SetAsBoundaryNode(true);
1612  }
1613  }
1614 }
1615 
1616 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1618 {
1619  assert(SPACE_DIM == 2); // LCOV_EXCL_LINE
1620  assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
1621 
1622 
1623  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elementIndex);
1624  unsigned num_nodes = p_element->GetNumNodes();
1625 
1626  std::set<unsigned> elements_containing_intersecting_node;
1627 
1628  for (unsigned node_local_index=0; node_local_index<num_nodes; node_local_index++)
1629  {
1630  unsigned node_global_index = p_element->GetNodeGlobalIndex(node_local_index);
1631 
1632  std::set<unsigned> node_elem_indices = this->GetNode(node_global_index)->rGetContainingElementIndices();
1633 
1634  for (std::set<unsigned>::const_iterator elem_iter = node_elem_indices.begin();
1635  elem_iter != node_elem_indices.end();
1636  ++elem_iter)
1637  {
1638  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_neighbouring_element = this->GetElement(*elem_iter);
1639  unsigned num_nodes_in_neighbouring_element = p_neighbouring_element->GetNumNodes();
1640 
1641  // Check if element contains the intersecting node
1642  for (unsigned node_index_2 = 0; node_index_2 < num_nodes_in_neighbouring_element; node_index_2++)
1643  {
1644  if (p_neighbouring_element->GetNodeGlobalIndex(node_index_2) == pNode->GetIndex())
1645  {
1646  elements_containing_intersecting_node.insert(p_neighbouring_element->GetIndex());
1647  }
1648  }
1649  }
1650  }
1651  /*
1652  * If there are not two elements containing the intersecting node then the node is coming from the other side of the element
1653  * and there is no way to fix it unless you want to make two new elements.
1654  */
1655  assert(elements_containing_intersecting_node.size() == 2);
1656 
1657  std::set<unsigned> all_elements_containing_intersecting_node = pNode->rGetContainingElementIndices();
1658 
1659  std::set<unsigned> intersecting_element;
1660 
1661  std::set_difference(all_elements_containing_intersecting_node.begin(), all_elements_containing_intersecting_node.end(),
1662  elements_containing_intersecting_node.begin(), elements_containing_intersecting_node.end(),
1663  std::inserter(intersecting_element, intersecting_element.begin()));
1664 
1665  /*
1666  * Identify nodes and elements to perform switch on
1667  * Intersecting node is node A
1668  * Other node is node B
1669  *
1670  * Element 1 only contains node A
1671  * Element 2 has nodes B and A (in that order)
1672  * Element 3 only contains node B
1673  * Element 4 has nodes A and B (in that order)
1674  */
1675  unsigned node_A_index = pNode->GetIndex();
1676  unsigned node_B_index;
1677  unsigned element_1_index = *(intersecting_element.begin());
1678  unsigned element_2_index;
1679  unsigned element_3_index = elementIndex;
1680  unsigned element_4_index;
1681 
1682  std::set<unsigned>::iterator iter = elements_containing_intersecting_node.begin();
1683  unsigned element_a_index = *(iter);
1684  iter++;
1685  unsigned element_b_index = *(iter);
1686 
1687  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_a = this->GetElement(element_a_index);
1688  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_b = this->GetElement(element_b_index);
1689 
1690  std::set<unsigned> element_a_nodes;
1691  for (unsigned node_index = 0; node_index < p_element_a->GetNumNodes(); node_index++)
1692  {
1693  element_a_nodes.insert(p_element_a->GetNodeGlobalIndex(node_index));
1694  }
1695 
1696  std::set<unsigned> element_b_nodes;
1697  for (unsigned node_index = 0; node_index < p_element_b->GetNumNodes(); node_index++)
1698  {
1699  element_b_nodes.insert(p_element_b->GetNodeGlobalIndex(node_index));
1700  }
1701 
1702  std::set<unsigned> switching_nodes;
1703  std::set_intersection(element_a_nodes.begin(), element_a_nodes.end(),
1704  element_b_nodes.begin(), element_b_nodes.end(),
1705  std::inserter(switching_nodes, switching_nodes.begin()));
1706 
1707  assert(switching_nodes.size() == 2);
1708 
1709  // Check intersecting node is this set
1710  assert(switching_nodes.find(node_A_index) != switching_nodes.end());
1711  switching_nodes.erase(node_A_index);
1712 
1713  assert(switching_nodes.size() == 1);
1714 
1715  node_B_index = *(switching_nodes.begin());
1716 
1717  // Now identify elements 2 and 4
1718  unsigned node_A_local_index_in_a = p_element_a->GetNodeLocalIndex(node_A_index);
1719  unsigned node_B_local_index_in_a = p_element_a->GetNodeLocalIndex(node_B_index);
1720 
1721  if ((node_B_local_index_in_a+1)%p_element_a->GetNumNodes() == node_A_local_index_in_a)
1722  {
1723  assert((p_element_b->GetNodeLocalIndex(node_A_index)+1)%p_element_b->GetNumNodes()
1724  == p_element_b->GetNodeLocalIndex(node_B_index));
1725 
1726  // Element 2 is element a, element 4 is element b
1727  element_2_index = element_a_index;
1728  element_4_index = element_b_index;
1729  }
1730  else
1731  {
1732  assert((p_element_b->GetNodeLocalIndex(node_B_index)+1)%p_element_b->GetNumNodes()
1733  == p_element_b->GetNodeLocalIndex(node_A_index));
1734 
1735  // Element 2 is element b, element 4 is element a
1736  element_2_index = element_b_index;
1737  element_4_index = element_a_index;
1738  }
1739 
1740  unsigned intersected_edge = this->GetLocalIndexForElementEdgeClosestToPoint(pNode->rGetLocation(), elementIndex);
1741 
1742  unsigned node_A_local_index_in_1 = this->GetElement(element_1_index)->GetNodeLocalIndex(node_A_index);
1743 
1744  unsigned node_A_local_index_in_2 = this->GetElement(element_2_index)->GetNodeLocalIndex(node_A_index);
1745  unsigned node_B_local_index_in_2 = this->GetElement(element_2_index)->GetNodeLocalIndex(node_B_index);
1746 
1747  unsigned node_B_local_index_in_3 = this->GetElement(elementIndex)->GetNodeLocalIndex(node_B_index);
1748 
1749  unsigned node_A_local_index_in_4 = this->GetElement(element_4_index)->GetNodeLocalIndex(node_A_index);
1750  unsigned node_B_local_index_in_4 = this->GetElement(element_4_index)->GetNodeLocalIndex(node_B_index);
1751 
1752  if (intersected_edge==node_B_local_index_in_3)
1753  {
1754  /*
1755  * Add node B to element 1 after node A
1756  * Add node A to element 3 after node B
1757  *
1758  * Remove node B from element 2
1759  * Remove node A from element 4
1760  */
1761  this->mElements[element_1_index]->AddNode(this->mNodes[node_B_index], node_A_local_index_in_1);
1762  this->mElements[element_3_index]->AddNode(this->mNodes[node_A_index], node_B_local_index_in_3);
1763 
1764  this->mElements[element_2_index]->DeleteNode(node_B_local_index_in_2);
1765  this->mElements[element_4_index]->DeleteNode(node_A_local_index_in_4);
1766  }
1767  else
1768  {
1769  assert((intersected_edge+1)%num_nodes==node_B_local_index_in_3);
1770 
1771  // Add node B to element 1 before node A and add node A to element 3 before node B
1772  unsigned node_before_A_in_1 = (node_A_local_index_in_1 - 1)%this->GetElement(element_1_index)->GetNumNodes();
1773  unsigned node_before_B_in_3 = (node_B_local_index_in_3 - 1)%this->GetElement(element_3_index)->GetNumNodes();
1774  this->mElements[element_1_index]->AddNode(this->mNodes[node_B_index], node_before_A_in_1);
1775  this->mElements[element_3_index]->AddNode(this->mNodes[node_A_index], node_before_B_in_3);
1776 
1777  // Remove node A from element 2 and remove node B from element 4
1778  this->mElements[element_2_index]->DeleteNode(node_A_local_index_in_2);
1779  this->mElements[element_4_index]->DeleteNode(node_B_local_index_in_4);
1780  }
1781 }
1782 
1783 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1785 {
1786  // The given element must be triangular for us to be able to perform a T2 swap on it
1787  assert(rElement.GetNumNodes() == 3);
1788 
1789  // Note that we define this vector before setting it, as otherwise the profiling build will break (see #2367)
1790  c_vector<double, SPACE_DIM> new_node_location;
1791  new_node_location = this->GetCentroidOfElement(rElement.GetIndex());
1792  mLastT2SwapLocation = new_node_location;
1793 
1794  // Create a new node at the element's centroid; this will be a boundary node if any existing nodes were on the boundary
1795  bool is_node_on_boundary = false;
1796  for (unsigned i=0; i<3; i++)
1797  {
1798  if (rElement.GetNode(i)->IsBoundaryNode())
1799  {
1800  is_node_on_boundary = true;
1801  break;
1802  }
1803  }
1804  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(GetNumNodes(), new_node_location, is_node_on_boundary));
1805  Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
1806 
1807  // Loop over each of the three nodes contained in rElement
1808  for (unsigned i=0; i<3; i++)
1809  {
1810  // For each node, find the set of other elements containing it
1811  Node<SPACE_DIM>* p_node = rElement.GetNode(i);
1812  std::set<unsigned> containing_elements = p_node->rGetContainingElementIndices();
1813  containing_elements.erase(rElement.GetIndex());
1814 
1815  // For each of these elements...
1816  for (std::set<unsigned>::iterator elem_iter = containing_elements.begin(); elem_iter != containing_elements.end(); ++elem_iter)
1817  {
1818  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_this_elem = this->GetElement(*elem_iter);
1819 
1820  // ...throw an exception if the element is triangular...
1821  if (p_this_elem->GetNumNodes() < 4)
1822  {
1823  EXCEPTION("One of the neighbours of a small triangular element is also a triangle - dealing with this has not been implemented yet");
1824  }
1825 
1826  // ...otherwise, replace p_node with p_new_node unless this has already happened (in which case, delete p_node from the element)
1827  if (p_this_elem->GetNodeLocalIndex(new_node_global_index) == UINT_MAX)
1828  {
1829  p_this_elem->ReplaceNode(p_node, p_new_node);
1830  }
1831  else
1832  {
1833  p_this_elem->DeleteNode(p_this_elem->GetNodeLocalIndex(p_node->GetIndex()));
1834  }
1835  }
1836  }
1837 
1838  // We also have to mark pElement, pElement->GetNode(0), pElement->GetNode(1), and pElement->GetNode(2) as deleted
1839  mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(0));
1840  mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(1));
1841  mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(2));
1842 
1843  rElement.GetNode(0)->MarkAsDeleted();
1844  rElement.GetNode(1)->MarkAsDeleted();
1845  rElement.GetNode(2)->MarkAsDeleted();
1846 
1847  mDeletedElementIndices.push_back(rElement.GetIndex());
1848  rElement.MarkAsDeleted();
1849 }
1850 
1851 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1853 {
1854  assert(SPACE_DIM == 2); // LCOV_EXCL_LINE - code will be removed at compile time
1855  assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE - code will be removed at compile time
1856 
1857  assert(pNode->IsBoundaryNode());
1858 
1859  // Store the index of the elements containing the intersecting node
1860  std::set<unsigned> elements_containing_intersecting_node = pNode->rGetContainingElementIndices();
1861 
1862  // Get the local index of the node in the intersected element after which the new node is to be added
1863  unsigned node_A_local_index = this->GetLocalIndexForElementEdgeClosestToPoint(pNode->rGetLocation(), elementIndex);
1864 
1865  // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
1866  c_vector<double, SPACE_DIM> node_location;
1867  node_location = pNode->rGetModifiableLocation();
1868 
1869  // Get element
1870  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elementIndex);
1871  unsigned num_nodes = p_element->GetNumNodes();
1872 
1873  // Get the nodes at either end of the edge to be divided
1874  unsigned vertexA_index = p_element->GetNodeGlobalIndex(node_A_local_index);
1875  unsigned vertexB_index = p_element->GetNodeGlobalIndex((node_A_local_index+1)%num_nodes);
1876 
1877  // Check these nodes are also boundary nodes if this fails then the elements have become concave and you need a smaller timestep
1878  if (!this->mNodes[vertexA_index]->IsBoundaryNode() || !this->mNodes[vertexB_index]->IsBoundaryNode())
1879  {
1880  EXCEPTION("A boundary node has intersected a non-boundary edge; this is because the boundary element has become concave. You need to rerun the simulation with a smaller time step to prevent this.");
1881  }
1882 
1883  // Get the nodes at either end of the edge to be divided and calculate intersection
1884  c_vector<double, SPACE_DIM> vertexA = p_element->GetNodeLocation(node_A_local_index);
1885  c_vector<double, SPACE_DIM> vertexB = p_element->GetNodeLocation((node_A_local_index+1)%num_nodes);
1886  c_vector<double, SPACE_DIM> vector_a_to_point = this->GetVectorFromAtoB(vertexA, node_location);
1887 
1888  c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
1889 
1890  c_vector<double, SPACE_DIM> edge_ab_unit_vector = vector_a_to_b/norm_2(vector_a_to_b);
1891  c_vector<double, SPACE_DIM> intersection = vertexA + edge_ab_unit_vector*inner_prod(vector_a_to_point, edge_ab_unit_vector);
1892 
1893  // Store the location of the T3 swap, the location of the intersection with the edge
1895  // is called (see #2401) - we should correct this in these cases!
1896 
1897  mLocationsOfT3Swaps.push_back(intersection);
1898 
1899  if (pNode->GetNumContainingElements() == 1)
1900  {
1901  // Get the index of the element containing the intersecting node
1902  unsigned intersecting_element_index = *elements_containing_intersecting_node.begin();
1903 
1904  // Get element
1905  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_intersecting_element = this->GetElement(intersecting_element_index);
1906 
1907  unsigned local_index = p_intersecting_element->GetNodeLocalIndex(pNode->GetIndex());
1908  unsigned next_node = p_intersecting_element->GetNodeGlobalIndex((local_index + 1)%(p_intersecting_element->GetNumNodes()));
1909  unsigned previous_node = p_intersecting_element->GetNodeGlobalIndex((local_index + p_intersecting_element->GetNumNodes() - 1)%(p_intersecting_element->GetNumNodes()));
1910 
1911  // Check to see if the nodes adjacent to the intersecting node are contained in the intersected element between vertices A and B
1912  if (next_node == vertexA_index || previous_node == vertexA_index || next_node == vertexB_index || previous_node == vertexB_index)
1913  {
1914  unsigned common_vertex_index;
1915 
1916  if (next_node == vertexA_index || previous_node == vertexA_index)
1917  {
1918  common_vertex_index = vertexA_index;
1919  }
1920  else
1921  {
1922  common_vertex_index = vertexB_index;
1923  }
1924 
1925  assert(this->mNodes[common_vertex_index]->GetNumContainingElements()>1);
1926 
1927  std::set<unsigned> elements_containing_common_vertex = this->mNodes[common_vertex_index]->rGetContainingElementIndices();
1928  std::set<unsigned>::const_iterator it = elements_containing_common_vertex.begin();
1929  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*it);
1930  it++;
1931  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*it);
1932 
1933  // Find the number and indices of common vertices between element_1 and element_2
1934  unsigned num_common_vertices = 0;
1935  std::vector<unsigned> common_vertex_indices;
1936  for (unsigned i=0; i<p_element_common_1->GetNumNodes(); i++)
1937  {
1938  for (unsigned j=0; j<p_element_common_2->GetNumNodes(); j++)
1939  {
1940  if (p_element_common_1->GetNodeGlobalIndex(i)==p_element_common_2->GetNodeGlobalIndex(j))
1941  {
1942  num_common_vertices++;
1943  common_vertex_indices.push_back(p_element_common_1->GetNodeGlobalIndex(i));
1944  }
1945  }
1946  }
1947 
1948  if (num_common_vertices == 1 || this->mNodes[common_vertex_index]->GetNumContainingElements() > 2)
1949  {
1950  /*
1951  * This is the situation here.
1952  *
1953  * From To
1954  * _ _
1955  * | <--- |
1956  * | /\ |\
1957  * | / \ | \
1958  * _|/____\ _|__\
1959  *
1960  * The edge goes from vertexA--vertexB to vertexA--pNode--vertexB
1961  */
1962 
1963  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
1964  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
1965 
1966  // Move original node
1967  pNode->rGetModifiableLocation() = intersection;
1968 
1969  // Add the moved nodes to the element (this also updates the node)
1970  this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
1971 
1972  // Check the nodes are updated correctly
1973  assert(pNode->GetNumContainingElements() == 2);
1974  }
1975  else if (num_common_vertices == 2)
1976  {
1977  // The two elements must have an edge in common. Find whether the common edge is the same as the
1978  // edge that is merged onto.
1979 
1980  if ((common_vertex_indices[0]==vertexA_index && common_vertex_indices[1]==vertexB_index) ||
1981  (common_vertex_indices[1]==vertexA_index && common_vertex_indices[0]==vertexB_index))
1982  {
1983  /*
1984  * Due to a previous T3 swap the situation looks like this.
1985  *
1986  * pNode
1987  * \ |\ /
1988  * \ | \ /
1989  * \_______|__\/
1990  * /A | B
1991  * / \
1992  *
1993  * A T3 Swap would merge pNode onto an edge of its own element.
1994  * We prevent this by just removing pNode. By doing this we also avoid the
1995  * intersecting element to be concave.
1996  */
1997 
1998  // Delete pNode in the intersecting element
1999  unsigned p_node_local_index = this->
2000  GetElement(intersecting_element_index)->GetNodeLocalIndex(pNode->GetIndex());
2001  this->GetElement(intersecting_element_index)->DeleteNode(p_node_local_index);
2002 
2003  // Mark all three nodes as deleted
2004  pNode->MarkAsDeleted();
2005  mDeletedNodeIndices.push_back(pNode->GetIndex());
2006  }
2007  else
2008  {
2009  /*
2010  * This is the situation here.
2011  *
2012  * C is common_vertex D is the other one.
2013  *
2014  * From To
2015  * _ D _
2016  * | <--- |
2017  * | /\ |\
2018  * C|/ \ | \
2019  * _|____\ _|__\
2020  *
2021  * The edge goes from vertexC--vertexB to vertexC--pNode--vertexD
2022  * then vertex B is removed as it is no longer needed.
2023  */
2024 
2025  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2026  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2027 
2028  // Move original node
2029  pNode->rGetModifiableLocation() = intersection;
2030 
2031  // Replace common_vertex with the the moved node (this also updates the nodes)
2032  this->GetElement(elementIndex)->ReplaceNode(this->mNodes[common_vertex_index], pNode);
2033 
2034  // Remove common_vertex
2035  unsigned common_vertex_local_index = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(common_vertex_index);
2036  this->GetElement(intersecting_element_index)->DeleteNode(common_vertex_local_index);
2037  assert(this->mNodes[common_vertex_index]->GetNumContainingElements() == 0);
2038 
2039  this->mNodes[common_vertex_index]->MarkAsDeleted();
2040  mDeletedNodeIndices.push_back(common_vertex_index);
2041 
2042  // Check the nodes are updated correctly
2043  assert(pNode->GetNumContainingElements() == 2);
2044  }
2045  }
2046  else if (num_common_vertices == 4)
2047  {
2048  /*
2049  * The two elements share edges CA and BD due to previous swaps but not the edge AB
2050  *
2051  * From To
2052  * D___ D___
2053  * | |
2054  * B|\ |
2055  * | \ |
2056  * | / |
2057  * A|/ |
2058  * C_|__ C_|__
2059  *
2060  * We just remove the intersecting node as well as vertices A and B.
2061  */
2062 
2063  // Delete node A and B in the intersected element
2064  this->GetElement(elementIndex)->DeleteNode(node_A_local_index);
2065  unsigned node_B_local_index = this->
2066  GetElement(elementIndex)->GetNodeLocalIndex(vertexB_index);
2067  this->GetElement(elementIndex)->DeleteNode(node_B_local_index);
2068 
2069  // Delete nodes A and B in the intersecting element
2070  unsigned node_A_local_index_intersecting_element = this->
2071  GetElement(intersecting_element_index)->GetNodeLocalIndex(vertexA_index);
2072  this->GetElement(intersecting_element_index)->DeleteNode(node_A_local_index_intersecting_element);
2073  unsigned node_B_local_index_intersecting_element = this->
2074  GetElement(intersecting_element_index)->GetNodeLocalIndex(vertexB_index);
2075  this->GetElement(intersecting_element_index)->DeleteNode(node_B_local_index_intersecting_element);
2076 
2077  // Delete pNode in the intersecting element
2078  unsigned p_node_local_index = this->
2079  GetElement(intersecting_element_index)->GetNodeLocalIndex(pNode->GetIndex());
2080  this->GetElement(intersecting_element_index)->DeleteNode(p_node_local_index);
2081 
2082  // Mark all three nodes as deleted
2083  pNode->MarkAsDeleted();
2084  mDeletedNodeIndices.push_back(pNode->GetIndex());
2085  this->mNodes[vertexA_index]->MarkAsDeleted();
2086  mDeletedNodeIndices.push_back(vertexA_index);
2087  this->mNodes[vertexB_index]->MarkAsDeleted();
2088  mDeletedNodeIndices.push_back(vertexB_index);
2089  }
2090  else
2091  {
2092  // This can't happen as nodes can't be on the internal edge of 2 elements.
2093  NEVER_REACHED;
2094  }
2095  }
2096  else
2097  {
2098  /*
2099  * From To
2100  * ____ _______
2101  * / \
2102  * /\ ^ / \
2103  * / \ |
2104  *
2105  * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
2106  */
2107 
2108  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2109  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2110  edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index+1)%num_nodes);
2111 
2112  // Move original node
2113  pNode->rGetModifiableLocation() = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2114 
2115  // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2116  c_vector<double, SPACE_DIM> new_node_location;
2117  new_node_location = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2118 
2119  // Add new node which will always be a boundary node
2120  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
2121 
2122  // Add the moved and new nodes to the element (this also updates the node)
2123  this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2124  this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
2125 
2126  // Add the new node to the original element containing pNode (this also updates the node)
2127  this->GetElement(intersecting_element_index)->AddNode(this->mNodes[new_node_global_index], this->GetElement(intersecting_element_index)->GetNodeLocalIndex(pNode->GetIndex()));
2128 
2129  // The nodes must have been updated correctly
2130  assert(pNode->GetNumContainingElements() == 2);
2131  assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
2132  }
2133  }
2134  else if (pNode->GetNumContainingElements() == 2)
2135  {
2136  // Find the nodes contained in elements containing the intersecting node
2137  std::set<unsigned>::const_iterator it = elements_containing_intersecting_node.begin();
2138 
2139  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_1 = this->GetElement(*it);
2140  unsigned num_nodes_elem_1 = p_element_1->GetNumNodes();
2141  it++;
2142 
2143  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_2 = this->GetElement(*it);
2144  unsigned num_nodes_elem_2 = p_element_2->GetNumNodes();
2145 
2146  unsigned node_global_index = pNode->GetIndex();
2147 
2148  unsigned local_index_1 = p_element_1->GetNodeLocalIndex(node_global_index);
2149  unsigned next_node_1 = p_element_1->GetNodeGlobalIndex((local_index_1 + 1)%num_nodes_elem_1);
2150  unsigned previous_node_1 = p_element_1->GetNodeGlobalIndex((local_index_1 + num_nodes_elem_1 - 1)%num_nodes_elem_1);
2151 
2152  unsigned local_index_2 = p_element_2->GetNodeLocalIndex(node_global_index);
2153  unsigned next_node_2 = p_element_2->GetNodeGlobalIndex((local_index_2 + 1)%num_nodes_elem_2);
2154  unsigned previous_node_2 = p_element_2->GetNodeGlobalIndex((local_index_2 + num_nodes_elem_2 - 1)%num_nodes_elem_2);
2155 
2156  // Check to see if the nodes adjacent to the intersecting node are contained in the intersected element between vertices A and B
2157  if ((next_node_1 == vertexA_index || previous_node_1 == vertexA_index || next_node_2 == vertexA_index || previous_node_2 == vertexA_index) &&
2158  (next_node_1 == vertexB_index || previous_node_1 == vertexB_index || next_node_2 == vertexB_index || previous_node_2 == vertexB_index))
2159  {
2160  /*
2161  * Here we have
2162  * __
2163  * /| /
2164  * __ / | --> ___/
2165  * \ | \
2166  * \|__ \
2167  *
2168  * Where the node on the left has overlapped the edge A B
2169  *
2170  * Move p_node to the intersection on A B and merge AB and p_node
2171  */
2172 
2173  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2174  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2175  edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index+1)%num_nodes);
2176 
2177  // Check they are all boundary nodes
2178  assert(pNode->IsBoundaryNode());
2179  assert(this->mNodes[vertexA_index]->IsBoundaryNode());
2180  assert(this->mNodes[vertexB_index]->IsBoundaryNode());
2181 
2182  // Move p_node to the intersection with the edge AB
2183  pNode->rGetModifiableLocation() = intersection;
2184  pNode->SetAsBoundaryNode(false);
2185 
2186  // Add pNode to the intersected element
2187  this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2188 
2189  // Remove vertex A from elements
2190  std::set<unsigned> elements_containing_vertex_A = this->mNodes[vertexA_index]->rGetContainingElementIndices();
2191  for (std::set<unsigned>::const_iterator iter = elements_containing_vertex_A.begin();
2192  iter != elements_containing_vertex_A.end();
2193  iter++)
2194  {
2195  this->GetElement(*iter)->DeleteNode(this->GetElement(*iter)->GetNodeLocalIndex(vertexA_index));
2196  }
2197 
2198  // Remove vertex A from the mesh
2199  assert(this->mNodes[vertexA_index]->GetNumContainingElements() == 0);
2200  this->mNodes[vertexA_index]->MarkAsDeleted();
2201  mDeletedNodeIndices.push_back(vertexA_index);
2202 
2203  // Remove vertex B from elements
2204  std::set<unsigned> elements_containing_vertex_B = this->mNodes[vertexB_index]->rGetContainingElementIndices();
2205  for (std::set<unsigned>::const_iterator iter = elements_containing_vertex_B.begin();
2206  iter != elements_containing_vertex_B.end();
2207  iter++)
2208  {
2209  this->GetElement(*iter)->DeleteNode(this->GetElement(*iter)->GetNodeLocalIndex(vertexB_index));
2210  }
2211 
2212  // Remove vertex B from the mesh
2213  assert(this->mNodes[vertexB_index]->GetNumContainingElements()==0);
2214  this->mNodes[vertexB_index]->MarkAsDeleted();
2215  mDeletedNodeIndices.push_back(vertexB_index);
2216  }
2217  else
2218  {
2219  if (next_node_1 == vertexA_index || previous_node_1 == vertexA_index || next_node_2 == vertexA_index || previous_node_2 == vertexA_index)
2220  {
2221  // Get elements containing vertexA_index (the common vertex)
2222 
2223  assert(this->mNodes[vertexA_index]->GetNumContainingElements() > 1);
2224 
2225  std::set<unsigned> elements_containing_vertex_A = this->mNodes[vertexA_index]->rGetContainingElementIndices();
2226  std::set<unsigned>::const_iterator iter = elements_containing_vertex_A.begin();
2227  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*iter);
2228  iter++;
2229  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*iter);
2230 
2231  // Calculate the number of common vertices between element_1 and element_2
2232  unsigned num_common_vertices = 0;
2233  for (unsigned i=0; i<p_element_common_1->GetNumNodes(); i++)
2234  {
2235  for (unsigned j=0; j<p_element_common_2->GetNumNodes(); j++)
2236  {
2237  if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
2238  {
2239  num_common_vertices++;
2240  }
2241  }
2242  }
2243 
2244  if (num_common_vertices == 1 || this->mNodes[vertexA_index]->GetNumContainingElements() > 2)
2245  {
2246  /*
2247  * From To
2248  * _ B _ B
2249  * | <--- |
2250  * | /|\ |\
2251  * | / | \ | \
2252  * | / | \ |\ \
2253  * _|/___|___\ _|_\_\
2254  * A A
2255  *
2256  * The edge goes from vertexA--vertexB to vertexA--pNode--new_node--vertexB
2257  */
2258 
2259  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2260  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2261  edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index+1)%num_nodes);
2262 
2263  // Move original node and change to non-boundary node
2264  pNode->rGetModifiableLocation() = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2265  pNode->SetAsBoundaryNode(false);
2266 
2267  // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2268  c_vector<double, SPACE_DIM> new_node_location;
2269  new_node_location = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2270 
2271  // Add new node, which will always be a boundary node
2272  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
2273 
2274  // Add the moved nodes to the element (this also updates the node)
2275  this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
2276  this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2277 
2278  // Add the new nodes to the original elements containing pNode (this also updates the node)
2279  if (next_node_1 == previous_node_2)
2280  {
2281  p_element_1->AddNode(this->mNodes[new_node_global_index], (local_index_1 + p_element_1->GetNumNodes() - 1)%(p_element_1->GetNumNodes()));
2282  }
2283  else
2284  {
2285  assert(next_node_2 == previous_node_1);
2286 
2287  p_element_2->AddNode(this->mNodes[new_node_global_index], (local_index_2 + p_element_2->GetNumNodes() - 1)%(p_element_2->GetNumNodes()));
2288  }
2289 
2290  // Check the nodes are updated correctly
2291  assert(pNode->GetNumContainingElements() == 3);
2292  assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
2293  }
2294  else if (num_common_vertices == 2)
2295  {
2296  /*
2297  * From To
2298  * _ B _ B
2299  * |<--- |
2300  * | /|\ |\
2301  * |/ | \ | \
2302  * | | \ |\ \
2303  * _|__|___\ _|_\_\
2304  * A A
2305  *
2306  * The edge goes from vertexA--vertexB to vertexA--pNode--new_node--vertexB
2307  * then vertexA is removed
2308  */
2309 
2310  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2311  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2312  edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index+1)%num_nodes);
2313 
2314  // Move original node and change to non-boundary node
2315  pNode->rGetModifiableLocation() = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2316  pNode->SetAsBoundaryNode(false);
2317 
2318  // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2319  c_vector<double, SPACE_DIM> new_node_location;
2320  new_node_location = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2321 
2322  // Add new node, which will always be a boundary node
2323  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
2324 
2325  // Add the moved nodes to the element (this also updates the node)
2326  this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
2327  this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2328 
2329  // Add the new nodes to the original elements containing pNode (this also updates the node)
2330  if (next_node_1 == previous_node_2)
2331  {
2332  p_element_1->AddNode(this->mNodes[new_node_global_index], (local_index_1 + p_element_1->GetNumNodes() - 1)%(p_element_1->GetNumNodes()));
2333  }
2334  else
2335  {
2336  assert(next_node_2 == previous_node_1);
2337  p_element_2->AddNode(this->mNodes[new_node_global_index], (local_index_2 + p_element_2->GetNumNodes() - 1)%(p_element_2->GetNumNodes()));
2338  }
2339 
2340  // Remove vertex A from the mesh
2341  p_element_common_1->DeleteNode(p_element_common_1->GetNodeLocalIndex(vertexA_index));
2342  p_element_common_2->DeleteNode(p_element_common_2->GetNodeLocalIndex(vertexA_index));
2343 
2344  assert(this->mNodes[vertexA_index]->GetNumContainingElements()==0);
2345 
2346  this->mNodes[vertexA_index]->MarkAsDeleted();
2347  mDeletedNodeIndices.push_back(vertexA_index);
2348 
2349  // Check the nodes are updated correctly
2350  assert(pNode->GetNumContainingElements() == 3);
2351  assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
2352  }
2353  else
2354  {
2355  // This can't happen as nodes can't be on the internal edge of two elements
2356  NEVER_REACHED;
2357  }
2358  }
2359  else if (next_node_1 == vertexB_index || previous_node_1 == vertexB_index || next_node_2 == vertexB_index || previous_node_2 == vertexB_index)
2360  {
2361  // Get elements containing vertexB_index (the common vertex)
2362 
2363  assert(this->mNodes[vertexB_index]->GetNumContainingElements()>1);
2364 
2365  std::set<unsigned> elements_containing_vertex_B = this->mNodes[vertexB_index]->rGetContainingElementIndices();
2366  std::set<unsigned>::const_iterator iter = elements_containing_vertex_B.begin();
2367  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*iter);
2368  iter++;
2369  VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*iter);
2370 
2371  // Calculate the number of common vertices between element_1 and element_2
2372  unsigned num_common_vertices = 0;
2373  for (unsigned i=0; i<p_element_common_1->GetNumNodes(); i++)
2374  {
2375  for (unsigned j=0; j<p_element_common_2->GetNumNodes(); j++)
2376  {
2377  if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
2378  {
2379  num_common_vertices++;
2380  }
2381  }
2382  }
2383 
2384  if (num_common_vertices == 1 || this->mNodes[vertexB_index]->GetNumContainingElements() > 2)
2385  {
2386  /*
2387  * From To
2388  * _B_________ _B____
2389  * |\ | / | / /
2390  * | \ | / |/ /
2391  * | \ | / | /
2392  * | \|/ |/
2393  * _| <--- _|
2394  * A
2395  *
2396  * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
2397  */
2398 
2399  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2400  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2401  edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index+1)%num_nodes);
2402 
2403  // Move original node and change to non-boundary node
2404  pNode->rGetModifiableLocation() = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2405  pNode->SetAsBoundaryNode(false);
2406 
2407  // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2408  c_vector<double, SPACE_DIM> new_node_location;
2409  new_node_location = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2410 
2411  // Add new node which will always be a boundary node
2412  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
2413 
2414  // Add the moved nodes to the element (this also updates the node)
2415  this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2416  this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
2417 
2418  // Add the new nodes to the original elements containing pNode (this also updates the node)
2419  if (next_node_1 == previous_node_2)
2420  {
2421  p_element_2->AddNode(this->mNodes[new_node_global_index], local_index_2);
2422  }
2423  else
2424  {
2425  assert(next_node_2 == previous_node_1);
2426  p_element_1->AddNode(this->mNodes[new_node_global_index], local_index_1);
2427  }
2428 
2429  // Check the nodes are updated correctly
2430  assert(pNode->GetNumContainingElements() == 3);
2431  assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
2432  }
2433  else if (num_common_vertices == 2)
2434  {
2435  /*
2436  * From To
2437  * _B_______ _B____
2438  * | | / | / /
2439  * | | / |/ /
2440  * |\ | / | /
2441  * | \|/ |/
2442  * _| <--- _|
2443  * A
2444  *
2445  * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
2446  * then vertexB is removed
2447  */
2448 
2449  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2450  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2451  edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index+1)%num_nodes);
2452 
2453  // Move original node and change to non-boundary node
2454  pNode->rGetModifiableLocation() = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2455  pNode->SetAsBoundaryNode(false);
2456 
2457  // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2458  c_vector<double, SPACE_DIM> new_node_location;
2459  new_node_location = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2460 
2461  // Add new node which will always be a boundary node
2462  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
2463 
2464  // Add the moved nodes to the element (this also updates the node)
2465  this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2466  this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
2467 
2468  // Add the new nodes to the original elements containing pNode (this also updates the node)
2469  if (next_node_1 == previous_node_2)
2470  {
2471  p_element_2->AddNode(this->mNodes[new_node_global_index], local_index_2);
2472  }
2473  else
2474  {
2475  assert(next_node_2 == previous_node_1);
2476  p_element_1->AddNode(this->mNodes[new_node_global_index], local_index_1);
2477  }
2478 
2479  // Remove vertex B from the mesh
2480  p_element_common_1->DeleteNode(p_element_common_1->GetNodeLocalIndex(vertexB_index));
2481  p_element_common_2->DeleteNode(p_element_common_2->GetNodeLocalIndex(vertexB_index));
2482 
2483  assert(this->mNodes[vertexB_index]->GetNumContainingElements()==0);
2484 
2485  this->mNodes[vertexB_index]->MarkAsDeleted();
2486  mDeletedNodeIndices.push_back(vertexB_index);
2487 
2488  // Check the nodes are updated correctly
2489  assert(pNode->GetNumContainingElements() == 3);
2490  assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
2491  }
2492  else
2493  {
2494  // This can't happen as nodes can't be on the internal edge of two elements
2495  NEVER_REACHED;
2496  }
2497  }
2498  else
2499  {
2500  /*
2501  * From To
2502  * _____ _______
2503  * / | \
2504  * /|\ ^ / | \
2505  * / | \ |
2506  *
2507  * The edge goes from vertexA--vertexB to vertexA--new_node_1--pNode--new_node_2--vertexB
2508  */
2509 
2510  // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2511  intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2512  edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index+1)%num_nodes);
2513 
2514  // Move original node and change to non-boundary node
2515  pNode->rGetModifiableLocation() = intersection;
2516  pNode->SetAsBoundaryNode(false);
2517 
2518  c_vector<double, SPACE_DIM> new_node_1_location;
2519  new_node_1_location = intersection - mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2520  c_vector<double, SPACE_DIM> new_node_2_location;
2521  new_node_2_location = intersection + mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
2522 
2523  // Add new nodes which will always be boundary nodes
2524  unsigned new_node_1_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_1_location[0], new_node_1_location[1]));
2525  unsigned new_node_2_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_2_location[0], new_node_2_location[1]));
2526 
2527  // Add the moved and new nodes to the element (this also updates the node)
2528  this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_2_global_index], node_A_local_index);
2529  this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2530  this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_1_global_index], node_A_local_index);
2531 
2532  // Add the new nodes to the original elements containing pNode (this also updates the node)
2533  if (next_node_1 == previous_node_2)
2534  {
2535  p_element_1->AddNode(this->mNodes[new_node_2_global_index], (local_index_1 + p_element_1->GetNumNodes() - 1)%(p_element_1->GetNumNodes()));
2536  p_element_2->AddNode(this->mNodes[new_node_1_global_index], local_index_2);
2537  }
2538  else
2539  {
2540  assert(next_node_2 == previous_node_1);
2541 
2542  p_element_1->AddNode(this->mNodes[new_node_1_global_index], local_index_1);
2543  p_element_2->AddNode(this->mNodes[new_node_2_global_index], (local_index_2 + p_element_2->GetNumNodes() - 1)%(p_element_2->GetNumNodes()));
2544  }
2545 
2546  // Check the nodes are updated correctly
2547  assert(pNode->GetNumContainingElements() == 3);
2548  assert(this->mNodes[new_node_1_global_index]->GetNumContainingElements() == 2);
2549  assert(this->mNodes[new_node_2_global_index]->GetNumContainingElements() == 2);
2550  }
2551  }
2552  }
2553  else
2554  {
2555  EXCEPTION("Trying to merge a node, contained in more than 2 elements, into another element, this is not possible with the vertex mesh.");
2556  }
2557 }
2558 
2559 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2561 {
2562  // Calculate void centroid
2563  c_vector<double, SPACE_DIM> nodes_midpoint = pNodeA->rGetLocation()
2564  + this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation()) / 3.0
2565  + this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeC->rGetLocation()) / 3.0;
2566 
2567  /*
2568  * In two steps, merge nodes A, B and C into a single node. This is implemented in such a way that
2569  * the ordering of their indices does not matter.
2570  */
2571 
2572  PerformNodeMerge(pNodeA, pNodeB);
2573 
2574  Node<SPACE_DIM>* p_merged_node = pNodeB;
2575 
2576  if (pNodeB->IsDeleted())
2577  {
2578  p_merged_node = pNodeA;
2579  }
2580 
2581  PerformNodeMerge(pNodeC, p_merged_node);
2582 
2583  if (p_merged_node->IsDeleted())
2584  {
2585  p_merged_node = pNodeC;
2586  }
2587 
2588  p_merged_node->rGetModifiableLocation() = nodes_midpoint;
2589 
2590  // Tag remaining node as non-boundary
2591  p_merged_node->SetAsBoundaryNode(false);
2592 
2593  // Remove the deleted nodes and re-index
2594  RemoveDeletedNodes();
2595 }
2596 
2597 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2599 {
2600  unsigned node_a_rank = pNodeA->rGetContainingElementIndices().size();
2601  unsigned node_b_rank = pNodeB->rGetContainingElementIndices().size();
2602 
2603  if ((node_a_rank > 3) && (node_b_rank > 3))
2604  {
2605  // The code can't handle this case
2606  EXCEPTION("Both nodes involved in a swap event are contained in more than three elements");
2607  }
2608  else // the rosette degree should increase in this case
2609  {
2610  assert(node_a_rank > 3 || node_b_rank > 3);
2611  this->PerformRosetteRankIncrease(pNodeA, pNodeB);
2612  this->RemoveDeletedNodes();
2613  }
2614 }
2615 
2616 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2618 {
2619  /*
2620  * One of the nodes will have 3 containing element indices, the other
2621  * will have at least four. We first identify which node is which.
2622  */
2623 
2624  unsigned node_a_index = pNodeA->GetIndex();
2625  unsigned node_b_index = pNodeB->GetIndex();
2626 
2627  unsigned node_a_rank = pNodeA->rGetContainingElementIndices().size();
2628  unsigned node_b_rank = pNodeB->rGetContainingElementIndices().size();
2629 
2630  unsigned lo_rank_index = (node_a_rank < node_b_rank) ? node_a_index : node_b_index;
2631  unsigned hi_rank_index = (node_a_rank < node_b_rank) ? node_b_index : node_a_index;
2632 
2633  // Get pointers to the nodes, sorted by index
2634  Node<SPACE_DIM>* p_lo_rank_node = this->GetNode(lo_rank_index);
2635  Node<SPACE_DIM>* p_hi_rank_node = this->GetNode(hi_rank_index);
2636 
2637  // Find the sets of elements containing each of the nodes, sorted by index
2638  std::set<unsigned> lo_rank_elem_indices = p_lo_rank_node->rGetContainingElementIndices();
2639  std::set<unsigned> hi_rank_elem_indices = p_hi_rank_node->rGetContainingElementIndices();
2640 
2667  for (std::set<unsigned>::const_iterator it = lo_rank_elem_indices.begin();
2668  it != lo_rank_elem_indices.end();
2669  ++it)
2670  {
2671  // Find the local index of lo_rank_node in this element
2672  unsigned lo_rank_local_index = this->mElements[*it]->GetNodeLocalIndex(lo_rank_index);
2673  assert(lo_rank_local_index < UINT_MAX); // double check this element contains lo_rank_node
2674 
2675  /*
2676  * If this element already contains the hi_rank_node, we are in the situation of elements
2677  * C and D above, so we just remove lo_rank_node.
2678  *
2679  * Otherwise, we are in element E, so we must replace lo_rank_node with high-rank node,
2680  * and remove it from mNodes.
2681  *
2682  * We can check whether hi_rank_node is in this element using the set::count() method.
2683  */
2684 
2685  if (hi_rank_elem_indices.count(*it) > 0)
2686  {
2687  // Delete lo_rank_node from current element
2688  this->mElements[*it]->DeleteNode(lo_rank_local_index);
2689  }
2690  else
2691  {
2692  // Update lo_rank_node with all information (including index and location) of hi_rank_node
2693  this->mElements[*it]->UpdateNode(lo_rank_local_index, p_hi_rank_node);
2694  }
2695  }
2696 
2697  // Tidy up the mesh by ensuring the global instance of lo_rank_node is deleted
2698  assert(!(this->mNodes[lo_rank_index]->IsDeleted()));
2699  this->mNodes[lo_rank_index]->MarkAsDeleted();
2700  this->mDeletedNodeIndices.push_back(lo_rank_index);
2701 }
2702 
2703 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2705 {
2706  // Double check we are dealing with a protorosette
2707  assert(pProtorosetteNode->rGetContainingElementIndices().size() == 4);
2708 
2709  // Get random number (0, 1, 2 or 3), as the resolution axis is assumed to be random
2710  unsigned random_elem_increment = RandomNumberGenerator::Instance()->randMod(4);
2711 
2712  // Find global indices of elements around the protorosette node
2713  std::set<unsigned> protorosette_node_containing_elem_indices = pProtorosetteNode->rGetContainingElementIndices();
2714 
2715  // Select random element by advancing iterator a random number times
2716  std::set<unsigned>::const_iterator elem_index_iter(protorosette_node_containing_elem_indices.begin());
2717  advance(elem_index_iter, random_elem_increment);
2718 
2746  /*
2747  * We need to find the global indices of elements B, C and D. We do this with set intersections.
2748  */
2749 
2750  unsigned elem_a_idx = *elem_index_iter;
2751  unsigned elem_b_idx = UINT_MAX;
2752  unsigned elem_c_idx = UINT_MAX;
2753  unsigned elem_d_idx = UINT_MAX;
2754 
2755  // Get pointer to element we've chosen at random (element A)
2756  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_elem_a = this->GetElement(elem_a_idx);
2757 
2758  // Get all necessary info about element A and the protorosette node
2759  unsigned num_nodes_elem_a = p_elem_a->GetNumNodes();
2760  unsigned protorosette_node_global_idx = pProtorosetteNode->GetIndex();
2761  unsigned protorosette_node_local_idx = p_elem_a->GetNodeLocalIndex(protorosette_node_global_idx);
2762 
2763  // Find global indices of previous (cw) and next (ccw) nodes, locally, from the protorosette node, in element A
2764  unsigned prev_node_global_idx = p_elem_a->GetNodeGlobalIndex((protorosette_node_local_idx + num_nodes_elem_a - 1) % num_nodes_elem_a);
2765  unsigned next_node_global_idx = p_elem_a->GetNodeGlobalIndex((protorosette_node_local_idx + 1) % num_nodes_elem_a);
2766 
2767  // Get the set of elements the previous and next nodes are contained in
2768  Node<SPACE_DIM>* p_prev_node = this->GetNode(prev_node_global_idx);
2769  Node<SPACE_DIM>* p_next_node = this->GetNode(next_node_global_idx);
2770  std::set<unsigned> prev_node_elem_indices = p_prev_node->rGetContainingElementIndices();
2771  std::set<unsigned> next_node_elem_indices = p_next_node->rGetContainingElementIndices();
2772 
2773  // Perform set intersections with the set of element indices which the protorosette node is contained in
2774  std::set<unsigned> intersection_with_prev;
2775  std::set<unsigned> intersection_with_next;
2776 
2777  // This intersection should contain just global indices for elements A and B
2778  std::set_intersection(protorosette_node_containing_elem_indices.begin(),
2779  protorosette_node_containing_elem_indices.end(),
2780  prev_node_elem_indices.begin(),
2781  prev_node_elem_indices.end(),
2782  std::inserter(intersection_with_prev, intersection_with_prev.begin()));
2783 
2784  // This intersection should contain just global indices for elements A and D
2785  std::set_intersection(protorosette_node_containing_elem_indices.begin(),
2786  protorosette_node_containing_elem_indices.end(),
2787  next_node_elem_indices.begin(),
2788  next_node_elem_indices.end(),
2789  std::inserter(intersection_with_next, intersection_with_next.begin()));
2790 
2791  assert(intersection_with_prev.size() == 2);
2792  assert(intersection_with_next.size() == 2);
2793 
2794  // Get global index of element B
2795  if (*intersection_with_prev.begin() != elem_a_idx)
2796  {
2797  elem_b_idx = *(intersection_with_prev.begin());
2798  }
2799  else
2800  {
2801  elem_b_idx = *(++(intersection_with_prev.begin()));
2802  }
2803  assert(elem_b_idx < UINT_MAX);
2804 
2805  // Get global index of element D
2806  if (*intersection_with_next.begin() != elem_a_idx)
2807  {
2808  elem_d_idx = *(intersection_with_next.begin());
2809  }
2810  else
2811  {
2812  elem_d_idx = *(++(intersection_with_next.begin()));
2813  }
2814  assert(elem_d_idx < UINT_MAX);
2815 
2816  // By elimination, the remaining unassigned index in the original set must be global index of element C
2817  for (elem_index_iter = protorosette_node_containing_elem_indices.begin();
2818  elem_index_iter != protorosette_node_containing_elem_indices.end();
2819  ++elem_index_iter)
2820  {
2821  if ((*elem_index_iter != elem_a_idx) && (*elem_index_iter != elem_b_idx) && (*elem_index_iter != elem_d_idx))
2822  {
2823  elem_c_idx = *elem_index_iter;
2824  }
2825  }
2826  assert(elem_c_idx < UINT_MAX);
2827 
2842  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_elem_b = this->GetElement(elem_b_idx);
2843  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_elem_c = this->GetElement(elem_c_idx);
2844  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_elem_d = this->GetElement(elem_d_idx);
2845 
2846  double swap_distance = (this->mCellRearrangementRatio) * (this->mCellRearrangementThreshold);
2847 
2848  // Get normalized vectors to centre of elements A and B from protorosette node
2849  c_vector<double, SPACE_DIM> node_to_elem_a_centre = this->GetCentroidOfElement(elem_a_idx) - pProtorosetteNode->rGetLocation();
2850  node_to_elem_a_centre /= norm_2(node_to_elem_a_centre);
2851 
2852  c_vector<double, SPACE_DIM> node_to_elem_c_centre = this->GetCentroidOfElement(elem_c_idx) - pProtorosetteNode->rGetLocation();
2853  node_to_elem_c_centre /= norm_2(node_to_elem_c_centre);
2854 
2855  // Calculate new node locations
2856  c_vector<double, SPACE_DIM> new_location_of_protorosette_node = pProtorosetteNode->rGetLocation() + (0.5 * swap_distance) * node_to_elem_a_centre;
2857  c_vector<double, SPACE_DIM> location_of_new_node = pProtorosetteNode->rGetLocation() + (0.5 * swap_distance) * node_to_elem_c_centre;
2858 
2859  // Move protorosette node to new location
2860  pProtorosetteNode->rGetModifiableLocation() = new_location_of_protorosette_node;
2861 
2862  // Create new node in correct location
2863  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(this->GetNumNodes(), location_of_new_node, false));
2864  Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
2865 
2876  unsigned local_idx_elem_b = p_elem_b->GetNodeLocalIndex(protorosette_node_global_idx);
2877  local_idx_elem_b = (local_idx_elem_b + p_elem_b->GetNumNodes() - 1) % p_elem_b->GetNumNodes();
2878  unsigned local_idx_elem_c = p_elem_c->GetNodeLocalIndex(protorosette_node_global_idx);
2879  unsigned local_idx_elem_d = p_elem_d->GetNodeLocalIndex(protorosette_node_global_idx);
2880 
2881  p_elem_b->AddNode(p_new_node, local_idx_elem_b);
2882  p_elem_c->AddNode(p_new_node, local_idx_elem_c);
2883  p_elem_d->AddNode(p_new_node, local_idx_elem_d);
2884 
2885  // All that is left is to remove the original protorosette node from element C
2886  p_elem_c->DeleteNode(p_elem_c->GetNodeLocalIndex(protorosette_node_global_idx));
2887 }
2888 
2889 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2891 {
2892  unsigned rosette_rank = pRosetteNode->rGetContainingElementIndices().size();
2893 
2894  // Double check we're dealing with a rosette
2895  assert(rosette_rank > 4);
2896 
2897  // Get random number in [0, 1, ..., n) where n is rank of rosette, as the resolution axis is assumed to be random
2898  unsigned random_elem_increment = RandomNumberGenerator::Instance()->randMod(rosette_rank);
2899 
2900  // Find global indices of elements around the protorosette node
2901  std::set<unsigned> rosette_node_containing_elem_indices = pRosetteNode->rGetContainingElementIndices();
2902 
2903  // Select random element by advancing iterator a random number times
2904  std::set<unsigned>::const_iterator elem_index_iter(rosette_node_containing_elem_indices.begin());
2905  advance(elem_index_iter, random_elem_increment);
2906 
2945  /*
2946  * We need to find the global indices of elements N and P. We do this with set intersections.
2947  */
2948 
2949  // Get the vertex element S (which we randomly selected)
2950  unsigned elem_s_idx = *elem_index_iter;
2951  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_elem_s = this->GetElement(elem_s_idx);
2952 
2953  unsigned elem_n_idx = UINT_MAX;
2954  unsigned elem_p_idx = UINT_MAX;
2955 
2956  // Get all necessary info about element S and the rosette node
2957  unsigned num_nodes_elem_s = p_elem_s->GetNumNodes();
2958  unsigned rosette_node_global_idx = pRosetteNode->GetIndex();
2959  unsigned rosette_node_local_idx = p_elem_s->GetNodeLocalIndex(rosette_node_global_idx);
2960 
2961  // Find global indices of previous (cw) and next (ccw) nodes, locally, from the rosette node, in element S
2962  unsigned prev_node_global_idx = p_elem_s->GetNodeGlobalIndex((rosette_node_local_idx + num_nodes_elem_s - 1) % num_nodes_elem_s);
2963  unsigned next_node_global_idx = p_elem_s->GetNodeGlobalIndex((rosette_node_local_idx + 1) % num_nodes_elem_s);
2964 
2965  // Get the set of elements that the previous and next nodes are contained in
2966  Node<SPACE_DIM>* p_prev_node = this->GetNode(prev_node_global_idx);
2967  Node<SPACE_DIM>* p_next_node = this->GetNode(next_node_global_idx);
2968  std::set<unsigned> prev_node_elem_indices = p_prev_node->rGetContainingElementIndices();
2969  std::set<unsigned> next_node_elem_indices = p_next_node->rGetContainingElementIndices();
2970 
2971  // Perform set intersections with the set of element indices that the rosette node is contained in
2972  std::set<unsigned> intersection_with_prev;
2973  std::set<unsigned> intersection_with_next;
2974 
2975  // This intersection should contain just global indices for elements S and N
2976  std::set_intersection(rosette_node_containing_elem_indices.begin(),
2977  rosette_node_containing_elem_indices.end(),
2978  prev_node_elem_indices.begin(),
2979  prev_node_elem_indices.end(),
2980  std::inserter(intersection_with_prev, intersection_with_prev.begin()));
2981 
2982  // This intersection should contain just global indices for elements S and P
2983  std::set_intersection(rosette_node_containing_elem_indices.begin(),
2984  rosette_node_containing_elem_indices.end(),
2985  next_node_elem_indices.begin(),
2986  next_node_elem_indices.end(),
2987  std::inserter(intersection_with_next, intersection_with_next.begin()));
2988 
2989  assert(intersection_with_prev.size() == 2);
2990  assert(intersection_with_next.size() == 2);
2991 
2992  // Get global index of element N
2993  if (*intersection_with_prev.begin() != elem_s_idx)
2994  {
2995  elem_n_idx = *intersection_with_prev.begin();
2996  }
2997  else
2998  {
2999  elem_n_idx = *(++(intersection_with_prev.begin()));
3000  }
3001  assert(elem_n_idx < UINT_MAX);
3002 
3003  // Get global index of element P
3004  if (*intersection_with_next.begin() != elem_s_idx)
3005  {
3006  elem_p_idx = *intersection_with_next.begin();
3007  }
3008  else
3009  {
3010  elem_p_idx = *(++(intersection_with_next.begin()));
3011  }
3012  assert(elem_p_idx < UINT_MAX);
3013 
3024  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_elem_n = this->GetElement(elem_p_idx);
3025  VertexElement<ELEMENT_DIM,SPACE_DIM>* p_elem_p = this->GetElement(elem_n_idx);
3026 
3027  double swap_distance = (this->mCellRearrangementRatio) * (this->mCellRearrangementThreshold);
3028 
3029  // Calculate location of new node
3030  c_vector<double, 2> node_to_selected_elem = this->GetCentroidOfElement(elem_s_idx) - pRosetteNode->rGetLocation();
3031  node_to_selected_elem /= norm_2(node_to_selected_elem);
3032  c_vector<double, 2> new_node_location = pRosetteNode->rGetLocation() + (swap_distance * node_to_selected_elem);
3033 
3034  // Create new node in correct location
3035  unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(this->GetNumNodes(), new_node_location, false));
3036  Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
3037 
3049  // Add new node, and remove rosette node, from element S
3050  unsigned node_local_idx_in_elem_s = p_elem_s->GetNodeLocalIndex(rosette_node_global_idx);
3051  p_elem_s->AddNode(p_new_node, node_local_idx_in_elem_s);
3052  p_elem_s->DeleteNode(node_local_idx_in_elem_s);
3053 
3054  // Add new node to element N
3055  unsigned node_local_idx_in_elem_n = p_elem_n->GetNodeLocalIndex(rosette_node_global_idx);
3056  node_local_idx_in_elem_n = (node_local_idx_in_elem_n + p_elem_n->GetNumNodes() - 1) % p_elem_n->GetNumNodes();
3057  p_elem_n->AddNode(p_new_node, node_local_idx_in_elem_n);
3058 
3059  // Add new node to element P
3060  unsigned node_local_idx_in_elem_p = p_elem_p->GetNodeLocalIndex(rosette_node_global_idx);
3061  p_elem_p->AddNode(p_new_node, node_local_idx_in_elem_p);
3062 }
3063 
3064 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3066 {
3075  // Vectors to store the nodes that need resolution events
3076  std::vector<Node<SPACE_DIM>* > protorosette_nodes;
3077  std::vector<Node<SPACE_DIM>* > rosette_nodes;
3078 
3079  // First loop in which we populate these vectors
3080  unsigned num_nodes = this->GetNumAllNodes();
3081  for (unsigned node_idx = 0 ; node_idx < num_nodes ; node_idx++)
3082  {
3083  Node<SPACE_DIM>* current_node = this->GetNode(node_idx);
3084  unsigned node_rank = current_node->rGetContainingElementIndices().size();
3085 
3086  if (node_rank < 4)
3087  {
3088  // Nothing to do if the node is not high-rank
3089  continue;
3090  }
3091  else if (node_rank == 4)
3092  {
3093  // For protorosette nodes, we check against a random number to decide if resolution is necessary
3094  if (mProtorosetteResolutionProbabilityPerTimestep >= RandomNumberGenerator::Instance()->ranf())
3095  {
3096  protorosette_nodes.push_back(current_node);
3097  }
3098  }
3099  else // if (node_rank > 4)
3100  {
3101  // For rosette nodes, we check against a random number to decide if resolution is necessary
3102  if (mRosetteResolutionProbabilityPerTimestep >= RandomNumberGenerator::Instance()->ranf())
3103  {
3104  rosette_nodes.push_back(current_node);
3105  }
3106  }
3107  }
3108 
3116  // First, resolve any protorosettes
3117  for (unsigned node_idx = 0 ; node_idx < protorosette_nodes.size() ; node_idx++)
3118  {
3119  Node<SPACE_DIM>* current_node = protorosette_nodes[node_idx];
3120 
3121  // Verify that node has not been marked for deletion, and that it is still contained in four elements
3122  assert( !(current_node->IsDeleted()) );
3123  assert( current_node->rGetContainingElementIndices().size() == 4 );
3124 
3125  // Perform protorosette resolution
3126  this->PerformProtorosetteResolution(current_node);
3127  }
3128 
3129  // Finally, resolve any rosettes
3130  for (unsigned node_idx = 0 ; node_idx < rosette_nodes.size() ; node_idx++)
3131  {
3132  Node<SPACE_DIM>* current_node = rosette_nodes[node_idx];
3133 
3134  // Verify that node has not been marked for deletion, and that it is still contained in at least four elements
3135  assert( !(current_node->IsDeleted()) );
3136  assert( current_node->rGetContainingElementIndices().size() > 4 );
3137 
3138  // Perform protorosette resolution
3139  this->PerformRosetteRankDecrease(current_node);
3140  }
3141 }
3142 
3143 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3145  unsigned indexA, unsigned indexB, c_vector<double,2> intersection)
3146 {
3156  c_vector<double, SPACE_DIM> vertexA = this->GetNode(indexA)->rGetLocation();
3157  c_vector<double, SPACE_DIM> vertexB = this->GetNode(indexB)->rGetLocation();
3158  c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
3159 
3160  if (norm_2(vector_a_to_b) < 4.0*mCellRearrangementRatio*mCellRearrangementThreshold)
3161  {
3162  WARNING("Trying to merge a node onto an edge which is too small.");
3163 
3164  c_vector<double, SPACE_DIM> centre_a_and_b = vertexA + 0.5*vector_a_to_b;
3165 
3166  vertexA = centre_a_and_b - 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*vector_a_to_b/norm_2(vector_a_to_b);
3167  ChastePoint<SPACE_DIM> vertex_A_point(vertexA);
3168  SetNode(indexA, vertex_A_point);
3169 
3170  vertexB = centre_a_and_b + 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*vector_a_to_b/norm_2(vector_a_to_b);
3171  ChastePoint<SPACE_DIM> vertex_B_point(vertexB);
3172  SetNode(indexB, vertex_B_point);
3173 
3174  intersection = centre_a_and_b;
3175  }
3176 
3177  // Reset distances
3178  vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
3179  c_vector<double,2> edge_ab_unit_vector = vector_a_to_b/norm_2(vector_a_to_b);
3180 
3181  // Reset the intersection away from vertices A and B to allow enough room for new nodes
3190  if (norm_2(intersection - vertexA) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
3191  {
3192  intersection = vertexA + 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
3193  }
3194  if (norm_2(intersection - vertexB) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
3195  {
3196  intersection = vertexB - 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
3197  }
3198  return intersection;
3199 }
3200 
3201 // Explicit instantiation
3202 template class MutableVertexMesh<1,1>;
3203 template class MutableVertexMesh<1,2>;
3204 template class MutableVertexMesh<1,3>;
3205 template class MutableVertexMesh<2,2>;
3206 template class MutableVertexMesh<2,3>;
3207 template class MutableVertexMesh<3,3>;
3208 
3209 // Serialization for Boost >= 1.36
unsigned DivideElement(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned nodeAIndex, unsigned nodeBIndex, bool placeOriginalElementBelow=false)
void SetRosetteResolutionProbabilityPerTimestep(double rosetteResolutionProbabilityPerTimestep)
double GetRosetteResolutionProbabilityPerTimestep() const
void PerformIntersectionSwap(Node< SPACE_DIM > *pNode, unsigned elementIndex)
double GetCellRearrangementRatio() const
void SetDistanceForT3SwapChecking(double distanceForT3SwapChecking)
virtual void Clear()
Definition: VertexMesh.cpp:579
c_vector< double, 2 > WidenEdgeOrCorrectIntersectionLocationIfNecessary(unsigned indexA, unsigned indexB, c_vector< double, 2 > intersection)
std::vector< c_vector< double, SPACE_DIM > > GetLocationsOfT3Swaps()
void RemoveDeletedNodesAndElements(VertexElementMap &rElementMap)
void SetPoint(ChastePoint< SPACE_DIM > point)
Definition: Node.cpp:115
void PerformRosetteRankDecrease(Node< SPACE_DIM > *pRosetteNode)
unsigned DivideElementAlongGivenAxis(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, c_vector< double, SPACE_DIM > axisOfDivision, bool placeOriginalElementBelow=false)
unsigned randMod(unsigned base)
Definition: Node.hpp:58
void SetAsBoundaryNode(bool value=true)
Definition: Node.cpp:127
bool CheckForT2Swaps(VertexElementMap &rElementMap)
unsigned GetNodeGlobalIndex(unsigned localIndex) const
std::vector< c_vector< double, SPACE_DIM > > GetLocationsOfT1Swaps()
unsigned DivideElementAlongShortAxis(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, bool placeOriginalElementBelow=false)
void PerformNodeMerge(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
virtual bool CheckForSwapsFromShortEdges()
void SetProtorosetteResolutionProbabilityPerTimestep(double protorosetteResolutionProbabilityPerTimestep)
bool IsBoundaryNode() const
Definition: Node.cpp:164
std::vector< VertexElement< ELEMENT_DIM-1, SPACE_DIM > * > mFaces
Definition: VertexMesh.hpp:83
#define EXCEPTION(message)
Definition: Exception.hpp:143
double GetT2Threshold() const
unsigned AddNode(Node< SPACE_DIM > *pNewNode)
void ReplaceNode(Node< SPACE_DIM > *pOldNode, Node< SPACE_DIM > *pNewNode)
double GetCellRearrangementThreshold() const
void PerformVoidRemoval(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB, Node< SPACE_DIM > *pNodeC)
std::set< unsigned > & rGetContainingElementIndices()
Definition: Node.cpp:300
void SetIndex(unsigned index)
Definition: Node.cpp:121
#define NEVER_REACHED
Definition: Exception.hpp:206
bool GetCheckForInternalIntersections() const
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
void SetCellRearrangementThreshold(double cellRearrangementThreshold)
bool mMeshChangesDuringSimulation
double GetProtorosetteResolutionProbabilityPerTimestep() const
void PerformT3Swap(Node< SPACE_DIM > *pNode, unsigned elementIndex)
void AddNode(Node< SPACE_DIM > *pNode, const unsigned &rIndex)
std::vector< Node< SPACE_DIM > * > mNodes
virtual void SetNode(unsigned nodeIndex, ChastePoint< SPACE_DIM > point)
unsigned AddElement(VertexElement< ELEMENT_DIM, SPACE_DIM > *pNewElement)
void PerformProtorosetteResolution(Node< SPACE_DIM > *pProtorosetteNode)
void SetDeleted(unsigned index)
double GetDistanceForT3SwapChecking() const
unsigned GetNumNodes() const
void PerformRosetteRankIncrease(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
unsigned GetNumNodes() const
void DeleteNodePriorToReMesh(unsigned index)
double GetProtorosetteFormationProbability() const
void SetNewIndex(unsigned oldIndex, unsigned newIndex)
static RandomNumberGenerator * Instance()
void RegisterWithNodes()
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
void SetProtorosetteFormationProbability(double protorosetteFormationProbability)
std::vector< VertexElement< ELEMENT_DIM, SPACE_DIM > * > mElements
Definition: VertexMesh.hpp:80
virtual void HandleHighOrderJunctions(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
void DivideEdge(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
void SetCheckForInternalIntersections(bool checkForInternalIntersections)
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition: Node.cpp:139
void AddElement(unsigned index)
Definition: Node.cpp:268
void DeleteElementPriorToReMesh(unsigned index)
unsigned GetIndex() const
unsigned GetNumElements() const
void SetCellRearrangementRatio(double cellRearrangementRatio)
void DeleteNode(const unsigned &rIndex)
unsigned GetNumContainingElements() const
Definition: Node.cpp:312
double GetNodeLocation(unsigned localIndex, unsigned dimension) const
void PerformT1Swap(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB, std::set< unsigned > &rElementsContainingNodes)
unsigned GetIndex() const
Definition: Node.cpp:158
unsigned GetNodeLocalIndex(unsigned globalIndex) const
void SetT2Threshold(double t2Threshold)
virtual void IdentifySwapType(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
c_vector< double, SPACE_DIM > & rGetModifiableLocation()
Definition: Node.cpp:151
void MarkAsDeleted()
Definition: Node.cpp:406
bool IsDeleted() const
Definition: Node.cpp:412
c_vector< double, SPACE_DIM > GetLastT2SwapLocation()
void PerformT2Swap(VertexElement< ELEMENT_DIM, SPACE_DIM > &rElement)
void Resize(unsigned size)