Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
MutableVertexMesh.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF 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"
41template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
42MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::MutableVertexMesh(std::vector<Node<SPACE_DIM>*> nodes,
43 std::vector<VertexElement<ELEMENT_DIM, SPACE_DIM>*> vertexElements,
44 double cellRearrangementThreshold,
45 double t2Threshold,
46 double cellRearrangementRatio,
47 double protorosetteFormationProbability,
48 double protorosetteResolutionProbabilityPerTimestep,
49 double rosetteResolutionProbabilityPerTimestep)
50 : mCellRearrangementThreshold(cellRearrangementThreshold),
51 mCellRearrangementRatio(cellRearrangementRatio),
52 mT2Threshold(t2Threshold),
53 mProtorosetteFormationProbability(protorosetteFormationProbability),
54 mProtorosetteResolutionProbabilityPerTimestep(protorosetteResolutionProbabilityPerTimestep),
55 mRosetteResolutionProbabilityPerTimestep(rosetteResolutionProbabilityPerTimestep),
56 mCheckForInternalIntersections(false),
57 mCheckForT3Swaps(true),
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->GenerateEdgesFromElements(vertexElements);
121
122 this->mMeshChangesDuringSimulation = true;
123}
124
125template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
126MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::MutableVertexMesh()
127 : mCellRearrangementThreshold(0.01),
128 mCellRearrangementRatio(1.5),
129 mT2Threshold(0.001),
130 mProtorosetteFormationProbability(0.0),
131 mProtorosetteResolutionProbabilityPerTimestep(0.0),
132 mRosetteResolutionProbabilityPerTimestep(0.0),
133 mCheckForInternalIntersections(false),
134 mCheckForT3Swaps(true),
135 mDistanceForT3SwapChecking(5.0)
136{
137 // Note that the member variables initialised above will be overwritten as soon as archiving is complete
138 this->mMeshChangesDuringSimulation = true;
139 Clear();
140}
141
142template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
147
148template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
149double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetCellRearrangementThreshold() const
150{
151 return mCellRearrangementThreshold;
152}
153
154template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
155double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetT2Threshold() const
156{
157 return mT2Threshold;
158}
159
160template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
161double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetCellRearrangementRatio() const
162{
163 return mCellRearrangementRatio;
164}
165
166template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
167double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetProtorosetteFormationProbability() const
168{
169 return this->mProtorosetteFormationProbability;
170}
171
172template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
173double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetProtorosetteResolutionProbabilityPerTimestep() const
174{
175 return this->mProtorosetteResolutionProbabilityPerTimestep;
176}
177
178template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
179double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetRosetteResolutionProbabilityPerTimestep() const
180{
181 return this->mRosetteResolutionProbabilityPerTimestep;
182}
183
184template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
185void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetDistanceForT3SwapChecking(double distanceForT3SwapChecking)
186{
187 mDistanceForT3SwapChecking = distanceForT3SwapChecking;
188}
189
190template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
191double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetDistanceForT3SwapChecking() const
192{
193 return mDistanceForT3SwapChecking;
194}
195
196template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
197bool MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetCheckForInternalIntersections() const
198{
199 return mCheckForInternalIntersections;
200}
201
202template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
203bool MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetCheckForT3Swaps() const
204{
205 return mCheckForT3Swaps;
206}
207
208template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
209void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetCellRearrangementThreshold(double cellRearrangementThreshold)
210{
211 mCellRearrangementThreshold = cellRearrangementThreshold;
212}
213
214template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
215void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetT2Threshold(double t2Threshold)
216{
217 mT2Threshold = t2Threshold;
218}
219
220template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
221void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetCellRearrangementRatio(double cellRearrangementRatio)
222{
223 mCellRearrangementRatio = cellRearrangementRatio;
224}
225
226template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
227void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetProtorosetteFormationProbability(double protorosetteFormationProbability)
228{
229 // Check that the new value is in [0, 1]
230 if (protorosetteFormationProbability < 0.0)
231 {
232 EXCEPTION("Attempting to assign a negative probability.");
233 }
234 if (protorosetteFormationProbability > 1.0)
235 {
236 EXCEPTION("Attempting to assign a probability greater than one.");
237 }
238
239 // Assign the new value
240 mProtorosetteFormationProbability = protorosetteFormationProbability;
241}
242
243template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
244void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetProtorosetteResolutionProbabilityPerTimestep(double protorosetteResolutionProbabilityPerTimestep)
245{
246 // Check that the new value is in [0, 1]
247 if (protorosetteResolutionProbabilityPerTimestep < 0.0)
248 {
249 EXCEPTION("Attempting to assign a negative probability.");
250 }
251 if (protorosetteResolutionProbabilityPerTimestep > 1.0)
252 {
253 EXCEPTION("Attempting to assign a probability greater than one.");
254 }
255
256 // Assign the new value
257 mProtorosetteResolutionProbabilityPerTimestep = protorosetteResolutionProbabilityPerTimestep;
258}
259
260template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
261void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetRosetteResolutionProbabilityPerTimestep(double rosetteResolutionProbabilityPerTimestep)
262{
263 // Check that the new value is in [0, 1]
264 if (rosetteResolutionProbabilityPerTimestep < 0.0)
265 {
266 EXCEPTION("Attempting to assign a negative probability.");
267 }
268 if (rosetteResolutionProbabilityPerTimestep > 1.0)
269 {
270 EXCEPTION("Attempting to assign a probability greater than one.");
271 }
272
273 // Assign the new value
274 mRosetteResolutionProbabilityPerTimestep = rosetteResolutionProbabilityPerTimestep;
275}
276
277template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
278void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetCheckForInternalIntersections(bool checkForInternalIntersections)
279{
280 mCheckForInternalIntersections = checkForInternalIntersections;
281}
282
283template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
284void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetCheckForT3Swaps(bool checkForT3Swaps)
285{
286 mCheckForT3Swaps = checkForT3Swaps;
287}
288
289template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
290void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::Clear()
291{
292 mDeletedNodeIndices.clear();
293 mDeletedElementIndices.clear();
294
295 VertexMesh<ELEMENT_DIM, SPACE_DIM>::Clear();
296}
297
298template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
299unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNumNodes() const
300{
301 return this->mNodes.size() - mDeletedNodeIndices.size();
302}
303
304template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
305unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNumElements() const
306{
307 return this->mElements.size() - mDeletedElementIndices.size();
308}
309
310template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
311std::vector< c_vector<double, SPACE_DIM> > MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLocationsOfT1Swaps()
312{
313 std::vector<T1SwapInfo<SPACE_DIM> > swap_info = mOperationRecorder.GetT1SwapsInfo();
314 std::vector< c_vector<double, SPACE_DIM> > swap_locations;
315 for (unsigned i=0; i<swap_info.size(); ++i)
316 {
317 swap_locations.push_back(swap_info[i].mLocation);
318 }
319 return swap_locations;
320}
321
322template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
323c_vector<double, SPACE_DIM> MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLastT2SwapLocation()
324{
325 return mLastT2SwapLocation;
326}
327
328template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
329std::vector< c_vector<double, SPACE_DIM> > MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLocationsOfT3Swaps()
330{
331 std::vector<T3SwapInfo<SPACE_DIM> > swap_info = mOperationRecorder.GetT3SwapsInfo();
332 std::vector< c_vector<double, SPACE_DIM> > swap_locations;
333 for (unsigned i=0; i<swap_info.size(); ++i)
334 {
335 swap_locations.push_back(swap_info[i].mLocation);
336 }
337 return swap_locations;
338}
339
340template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
341std::vector< c_vector<double, SPACE_DIM> > MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLocationsOfIntersectionSwaps()
342{
343 return mLocationsOfIntersectionSwaps;
344}
345
346template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
347void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ClearLocationsOfT1Swaps()
348{
349 mOperationRecorder.ClearT1SwapsInfo();
350}
351
352template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
353void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ClearLocationsOfT3Swaps()
354{
355 mOperationRecorder.ClearT3SwapsInfo();
356}
357
358template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
359void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ClearLocationsOfIntersectionSwaps()
360{
361 mLocationsOfIntersectionSwaps.clear();
362}
363
364template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
365unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::AddNode(Node<SPACE_DIM>* pNewNode)
366{
367 if (mDeletedNodeIndices.empty())
368 {
369 pNewNode->SetIndex(this->mNodes.size());
370 this->mNodes.push_back(pNewNode);
371 }
372 else
373 {
374 unsigned index = mDeletedNodeIndices.back();
375 pNewNode->SetIndex(index);
376 mDeletedNodeIndices.pop_back();
377 delete this->mNodes[index];
378 this->mNodes[index] = pNewNode;
379 }
380 return pNewNode->GetIndex();
381}
382
383template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
384unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::AddElement(VertexElement<ELEMENT_DIM, SPACE_DIM>* pNewElement)
385{
386 unsigned new_element_index = pNewElement->GetIndex();
387
388 if (new_element_index == this->mElements.size())
389 {
390 this->mElements.push_back(pNewElement);
391 }
392 else
393 {
394 this->mElements[new_element_index] = pNewElement;
395 }
396
397 pNewElement->RegisterWithNodes();
398 pNewElement->SetEdgeHelper(&(this->mEdgeHelper));
399 pNewElement->BuildEdges();
400 return pNewElement->GetIndex();
401}
402
403template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
404void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetNode(unsigned nodeIndex, ChastePoint<SPACE_DIM> point)
405{
406 this->mNodes[nodeIndex]->SetPoint(point);
407}
408
409template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
410unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DivideElementAlongGivenAxis(
411 [[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>* pElement,
412 [[maybe_unused]] c_vector<double, SPACE_DIM> axisOfDivision,
413 [[maybe_unused]] bool placeOriginalElementBelow)
414{
415 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
416 {
417 assert(SPACE_DIM == 2); // LCOV_EXCL_LINE
418 assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
419
420 // Storing element edge map prior to division. This is needed for division recording.
421 std::vector<unsigned> edgeIds;
422 for (unsigned i = 0; i < pElement->GetNumEdges(); i++)
423 {
424 edgeIds.push_back(pElement->GetEdge(i)->GetIndex());
425 }
426
427 // Get the centroid of the element
428 c_vector<double, SPACE_DIM> centroid = this->GetCentroidOfElement(pElement->GetIndex());
429 // Create a vector perpendicular to the axis of division
430 c_vector<double, SPACE_DIM> perp_axis;
431 perp_axis(0) = -axisOfDivision(1);
432 perp_axis(1) = axisOfDivision(0);
433
434 /*
435 * Find which edges the axis of division crosses by finding any node
436 * that lies on the opposite side of the axis of division to its next
437 * neighbour.
438 */
439 unsigned num_nodes = pElement->GetNumNodes();
440 std::vector<unsigned> intersecting_nodes;
441 bool is_current_node_on_left = (inner_prod(this->GetVectorFromAtoB(pElement->GetNodeLocation(0), centroid), perp_axis) >= 0);
442 for (unsigned i = 0; i < num_nodes; i++)
443 {
444 bool is_next_node_on_left = (inner_prod(this->GetVectorFromAtoB(pElement->GetNodeLocation((i + 1) % num_nodes), centroid), perp_axis) >= 0);
445 if (is_current_node_on_left != is_next_node_on_left)
446 {
447 intersecting_nodes.push_back(i);
448 }
449 is_current_node_on_left = is_next_node_on_left;
450 }
451
452 // If the axis of division does not cross two edges then we cannot proceed
453 if (intersecting_nodes.size() != 2)
454 {
455 EXCEPTION("Cannot proceed with element division: the given axis of division does not cross two edges of the element");
456 }
457
458 std::vector<unsigned> division_node_global_indices;
459 unsigned nodes_added = 0;
460
461 // Keeps track of elements (first) whose edge (second) had been split
462 std::vector<std::pair<VertexElement<ELEMENT_DIM, SPACE_DIM>*, unsigned> > edge_split_pairs;
463 std::vector<double> relative_new_node;
464 // Find the intersections between the axis of division and the element edges
465 for (unsigned i = 0; i < intersecting_nodes.size(); i++)
466 {
467 /*
468 * Get pointers to the nodes forming the edge into which one new node will be inserted.
469 *
470 * Note that when we use the first entry of intersecting_nodes to add a node,
471 * we change the local index of the second entry of intersecting_nodes in
472 * pElement, so must account for this by moving one entry further on.
473 */
474 Node<SPACE_DIM>* p_node_A = pElement->GetNode((intersecting_nodes[i] + nodes_added) % pElement->GetNumNodes());
475 Node<SPACE_DIM>* p_node_B = pElement->GetNode((intersecting_nodes[i] + nodes_added + 1) % pElement->GetNumNodes());
476
477 // Find the indices of the elements owned by each node on the edge into which one new node will be inserted
478 std::set<unsigned> elems_containing_node_A = p_node_A->rGetContainingElementIndices();
479 std::set<unsigned> elems_containing_node_B = p_node_B->rGetContainingElementIndices();
480
481 c_vector<double, SPACE_DIM> position_a = p_node_A->rGetLocation();
482 c_vector<double, SPACE_DIM> position_b = p_node_B->rGetLocation();
483 c_vector<double, SPACE_DIM> a_to_b = this->GetVectorFromAtoB(position_a, position_b);
484
485 c_vector<double, SPACE_DIM> intersection;
486
487 if (norm_2(a_to_b) < 2.0 * mCellRearrangementRatio * mCellRearrangementThreshold)
488 {
489 WARNING("Edge is too small for normal division; putting node in the middle of a and b. There may be T1 swaps straight away.");
491 intersection = position_a + 0.5 * a_to_b;
492 }
493 else
494 {
495 // Find the location of the intersection
496 double determinant = a_to_b[0] * axisOfDivision[1] - a_to_b[1] * axisOfDivision[0];
497
498 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
499 c_vector<double, SPACE_DIM> moved_centroid;
500 moved_centroid = position_a + this->GetVectorFromAtoB(position_a, centroid);
501
502 double alpha = (moved_centroid[0] * a_to_b[1] - position_a[0] * a_to_b[1]
503 - moved_centroid[1] * a_to_b[0] + position_a[1] * a_to_b[0])
504 / determinant;
505
506 intersection = moved_centroid + alpha * axisOfDivision;
507
508 /*
509 * If then new node is too close to one of the edge nodes, then reposition it
510 * a distance mCellRearrangementRatio*mCellRearrangementThreshold further along the edge.
511 */
512 c_vector<double, SPACE_DIM> a_to_intersection = this->GetVectorFromAtoB(position_a, intersection);
513 if (norm_2(a_to_intersection) < mCellRearrangementThreshold)
514 {
515 intersection = position_a + mCellRearrangementRatio * mCellRearrangementThreshold * a_to_b / norm_2(a_to_b);
516 }
517
518 c_vector<double, SPACE_DIM> b_to_intersection = this->GetVectorFromAtoB(position_b, intersection);
519 if (norm_2(b_to_intersection) < mCellRearrangementThreshold)
520 {
521 assert(norm_2(a_to_intersection) > mCellRearrangementThreshold); // to prevent moving intersection back to original position
522
523 intersection = position_b - mCellRearrangementRatio * mCellRearrangementThreshold * a_to_b / norm_2(a_to_b);
524 }
525 }
526
527 /*
528 * The new node is boundary node if the 2 nodes are boundary nodes and the elements don't look like
529 * ___A___
530 * | | |
531 * |___|___|
532 * B
533 */
534 bool is_boundary = false;
535 if (p_node_A->IsBoundaryNode() && p_node_B->IsBoundaryNode())
536 {
537 if (elems_containing_node_A.size() != 2 || elems_containing_node_B.size() != 2 || elems_containing_node_A != elems_containing_node_B)
538 {
539 is_boundary = true;
540 }
541 }
542
543 // Add a new node to the mesh at the location of the intersection
544 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, is_boundary, intersection[0], intersection[1]));
545 nodes_added++;
546
547 // Now make sure the new node is added to all neighbouring elements
548
549 // Find common elements
550 std::set<unsigned> shared_elements;
551 std::set_intersection(elems_containing_node_A.begin(),
552 elems_containing_node_A.end(),
553 elems_containing_node_B.begin(),
554 elems_containing_node_B.end(),
555 std::inserter(shared_elements, shared_elements.begin()));
556
557 // Iterate over common elements, including the element to be divided
558 unsigned node_A_index = p_node_A->GetIndex();
559 unsigned node_B_index = p_node_B->GetIndex();
560 bool original_element = false;
561 for (std::set<unsigned>::iterator iter = shared_elements.begin();
562 iter != shared_elements.end();
563 ++iter)
564 {
565 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(*iter);
566 original_element = p_element->GetIndex() == pElement->GetIndex();
567
568 // Find which node has the lower local index in this element
569 unsigned local_indexA = p_element->GetNodeLocalIndex(node_A_index);
570 unsigned local_indexB = p_element->GetNodeLocalIndex(node_B_index);
571
572 unsigned index = local_indexB;
573
574 // If node B has a higher index then use node A's index...
575 if (local_indexB > local_indexA)
576 {
577 index = local_indexA;
578
579 // ...unless nodes A and B share the element's last edge
580 if ((local_indexA == 0) && (local_indexB == p_element->GetNumNodes() - 1))
581 {
582 index = local_indexB;
583 }
584 }
585 else if ((local_indexB == 0) && (local_indexA == p_element->GetNumNodes() - 1))
586 {
587 // ...otherwise use node B's index, unless nodes A and B share the element's last edge
588 index = local_indexA;
589 }
590
591 // Record the edge split in neighbouring elements
592 if (!original_element && mTrackMeshOperations)
593 {
594 const unsigned num_edges = p_element->GetNumEdges();
595 const unsigned next_index = (index + 1) % num_edges;
596
597 auto prev_node = p_element->GetNode(index)->rGetLocation();
598 auto next_node = p_element->GetNode(next_index)->rGetLocation();
599 auto curr_node = this->GetNode(new_node_global_index)->rGetLocation();
600
601 c_vector<double, SPACE_DIM> last_to_next = this->GetVectorFromAtoB(prev_node, next_node);
602 c_vector<double, SPACE_DIM> last_to_curr = this->GetVectorFromAtoB(prev_node, curr_node);
603
604 double old_distance = norm_2(last_to_next);
605 double prev_curr_distance = norm_2(last_to_curr);
606
607 // Assert that new node in split edge is between a pair old nodes
608 assert(old_distance >= prev_curr_distance); // LCOV_EXCL_LINE
609
610 double theta = prev_curr_distance / old_distance;
611 relative_new_node.push_back(theta);
612 edge_split_pairs.push_back(std::pair<VertexElement<ELEMENT_DIM, SPACE_DIM>*, unsigned>(p_element, index));
613 }
614
615 // Add new node to this element
616 p_element->AddNode(this->GetNode(new_node_global_index), index);
617 }
618
619 // Store index of new node
620 division_node_global_indices.push_back(new_node_global_index);
621 }
622
623 // Now call DivideElement() to divide the element using the new nodes
624 unsigned new_element_index = DivideElement(pElement,
625 pElement->GetNodeLocalIndex(division_node_global_indices[0]),
626 pElement->GetNodeLocalIndex(division_node_global_indices[1]),
627 placeOriginalElementBelow);
628
629 // Record cell division info
630 CellDivisionInfo<SPACE_DIM> division_info;
631 division_info.mLocation = centroid;
632 division_info.mDaughterLocation1 = this->GetCentroidOfElement(pElement->GetIndex());
633 c_vector<double, SPACE_DIM> long_axis = this->GetShortAxisOfElement(pElement->GetIndex());
634 division_info.mDaughterLongAxis1(0) = -long_axis(1);
635 division_info.mDaughterLongAxis1(1) = long_axis(0);
636
637 division_info.mDaughterLocation2 = this->GetCentroidOfElement(new_element_index);
638 long_axis = this->GetShortAxisOfElement(new_element_index);
639 division_info.mDaughterLongAxis2(0) = -long_axis(1);
640 division_info.mDaughterLongAxis2(1) = long_axis(0);
641
642 division_info.mDivisionAxis = axisOfDivision;
643 mOperationRecorder.RecordCellDivisionInfo(division_info);
644
645 if (mTrackMeshOperations)
646 {
647 // Record edge rearrangements in the daughter cells ...
648 mOperationRecorder.RecordCellDivideOperation(edgeIds, pElement, this->mElements[new_element_index]);
649
650 // ... and in each neighbouring cell
651 for (unsigned i = 0; i < edge_split_pairs.size(); ++i)
652 {
653 mOperationRecorder.RecordEdgeSplitOperation(edge_split_pairs[i].first,
654 edge_split_pairs[i].second,
655 relative_new_node[i], true);
656 }
657 }
658
659 return new_element_index;
660 }
661 else
662 {
663 NEVER_REACHED;
664 }
665}
666
667template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
668unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DivideElementAlongShortAxis(
669 [[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>* pElement,
670 [[maybe_unused]] bool placeOriginalElementBelow)
671{
672 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
673 {
674 c_vector<double, SPACE_DIM> short_axis = this->GetShortAxisOfElement(pElement->GetIndex());
675
676 unsigned new_element_index = DivideElementAlongGivenAxis(pElement, short_axis, placeOriginalElementBelow);
677 return new_element_index;
678 }
679 else
680 {
681 NEVER_REACHED;
682 }
683}
684
685template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
686unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DivideElement(
687 [[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>* pElement,
688 [[maybe_unused]] unsigned nodeAIndex,
689 [[maybe_unused]] unsigned nodeBIndex,
690 [[maybe_unused]] bool placeOriginalElementBelow)
691{
692 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
693 {
694 // Sort nodeA and nodeB such that nodeBIndex > nodeAindex
695 assert(nodeBIndex != nodeAIndex);
696 unsigned node1_index = (nodeAIndex < nodeBIndex) ? nodeAIndex : nodeBIndex; // low index
697 unsigned node2_index = (nodeAIndex < nodeBIndex) ? nodeBIndex : nodeAIndex; // high index
698
699 // Store the number of nodes in the element (this changes when nodes are deleted from the element)
700 unsigned num_nodes = pElement->GetNumNodes();
701
702 // Copy the nodes in this element
703 std::vector<Node<SPACE_DIM>*> nodes_elem;
704 for (unsigned i = 0; i < num_nodes; i++)
705 {
706 nodes_elem.push_back(pElement->GetNode(i));
707 }
708
709 // Get the index of the new element
710 unsigned new_element_index;
711 if (mDeletedElementIndices.empty())
712 {
713 new_element_index = this->mElements.size();
714 }
715 else
716 {
717 new_element_index = mDeletedElementIndices.back();
718 mDeletedElementIndices.pop_back();
719 delete this->mElements[new_element_index];
720 }
721
722 // Add the new element to the mesh
723 AddElement(new VertexElement<ELEMENT_DIM, SPACE_DIM>(new_element_index, nodes_elem));
724
731 // Find lowest element
733 double height_midpoint_1 = 0.0;
734 double height_midpoint_2 = 0.0;
735 unsigned counter_1 = 0;
736 unsigned counter_2 = 0;
737
738 for (unsigned i = 0; i < num_nodes; i++)
739 {
740 if (i >= node1_index && i <= node2_index)
741 {
742 height_midpoint_1 += pElement->GetNode(i)->rGetLocation()[1];
743 counter_1++;
744 }
745 if (i <= node1_index || i >= node2_index)
746 {
747 height_midpoint_2 += pElement->GetNode(i)->rGetLocation()[1];
748 counter_2++;
749 }
750 }
751 height_midpoint_1 /= (double)counter_1;
752 height_midpoint_2 /= (double)counter_2;
753
754 for (unsigned i = num_nodes; i > 0; i--)
755 {
756 if (i - 1 < node1_index || i - 1 > node2_index)
757 {
758 if (height_midpoint_1 < height_midpoint_2)
759 {
760 if (placeOriginalElementBelow)
761 {
762 pElement->DeleteNode(i - 1);
763 }
764 else
765 {
766 this->mElements[new_element_index]->DeleteNode(i - 1);
767 }
768 }
769 else
770 {
771 if (placeOriginalElementBelow)
772 {
773 this->mElements[new_element_index]->DeleteNode(i - 1);
774 }
775 else
776 {
777 pElement->DeleteNode(i - 1);
778 }
779 }
780 }
781 else if (i - 1 > node1_index && i - 1 < node2_index)
782 {
783 if (height_midpoint_1 < height_midpoint_2)
784 {
785 if (placeOriginalElementBelow)
786 {
787 this->mElements[new_element_index]->DeleteNode(i - 1);
788 }
789 else
790 {
791 pElement->DeleteNode(i - 1);
792 }
793 }
794 else
795 {
796 if (placeOriginalElementBelow)
797 {
798 pElement->DeleteNode(i - 1);
799 }
800 else
801 {
802 this->mElements[new_element_index]->DeleteNode(i - 1);
803 }
804 }
805 }
806 }
807 // Re-build edges when division is performed
808 this->mElements[new_element_index]->RebuildEdges();
809 pElement->RebuildEdges();
810 return new_element_index;
811 }
812 else
813 {
814 NEVER_REACHED;
815 }
816}
817
818template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
819void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DeleteElementPriorToReMesh([[maybe_unused]] unsigned index)
820{
821 if constexpr (SPACE_DIM == 2)
822 {
823 // Mark any nodes that are contained only in this element as deleted
824 for (unsigned i = 0; i < this->mElements[index]->GetNumNodes(); i++)
825 {
826 Node<SPACE_DIM>* p_node = this->mElements[index]->GetNode(i);
827
828 if (p_node->rGetContainingElementIndices().size() == 1)
829 {
830 DeleteNodePriorToReMesh(p_node->GetIndex());
831 }
832
833 // Mark all the nodes contained in the removed element as boundary nodes
834 p_node->SetAsBoundaryNode(true);
835 }
836
837 // Mark this element as deleted
838 this->mElements[index]->MarkAsDeleted();
839 mDeletedElementIndices.push_back(index);
840 }
841 else
842 {
843 NEVER_REACHED;
844 }
845}
846
847template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
848void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DeleteNodePriorToReMesh(unsigned index)
849{
850 this->mNodes[index]->MarkAsDeleted();
851 mDeletedNodeIndices.push_back(index);
852}
853
854template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
855void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DivideEdge(Node<SPACE_DIM>* pNodeA, Node<SPACE_DIM>* pNodeB)
856{
857 // Find the indices of the elements owned by each node
858 std::set<unsigned> elements_containing_nodeA = pNodeA->rGetContainingElementIndices();
859 std::set<unsigned> elements_containing_nodeB = pNodeB->rGetContainingElementIndices();
860
861 // Find common elements
862 std::set<unsigned> shared_elements;
863 std::set_intersection(elements_containing_nodeA.begin(),
864 elements_containing_nodeA.end(),
865 elements_containing_nodeB.begin(),
866 elements_containing_nodeB.end(),
867 std::inserter(shared_elements, shared_elements.begin()));
868
869 // Check that the nodes have a common edge and not more than 2
870 assert(!shared_elements.empty());
871 assert(shared_elements.size()<=2u);
872
873 // Specify if it's a boundary node
874 bool is_boundary_node = false;
875 if (shared_elements.size()==1u)
876 {
877 // If only one shared element then must be on the boundary.
878 assert((pNodeA->IsBoundaryNode()) && (pNodeB->IsBoundaryNode()));
879 is_boundary_node = true;
880 }
881
882 // Create a new node (position is not important as it will be changed)
883 Node<SPACE_DIM>* p_new_node = new Node<SPACE_DIM>(GetNumNodes(), is_boundary_node, 0.0, 0.0);
884
885 // Update the node location
886 c_vector<double, SPACE_DIM> new_node_position = pNodeA->rGetLocation() + 0.5*this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation());
887 ChastePoint<SPACE_DIM> point(new_node_position);
888 p_new_node->SetPoint(new_node_position);
889
890 // Add node to mesh
891 this->mNodes.push_back(p_new_node);
892
893 // Iterate over common elements
894 unsigned node_A_index = pNodeA->GetIndex();
895 unsigned node_B_index = pNodeB->GetIndex();
896 for (std::set<unsigned>::iterator iter = shared_elements.begin();
897 iter != shared_elements.end();
898 ++iter)
899 {
900 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(*iter);
901
902 // Find which node has the lower local index in this element
903 unsigned local_indexA = p_element->GetNodeLocalIndex(node_A_index);
904 unsigned local_indexB = p_element->GetNodeLocalIndex(node_B_index);
905
906 unsigned index = local_indexB;
907
908 // If node B has a higher index then use node A's index...
909 if (local_indexB > local_indexA)
910 {
911 index = local_indexA;
912
913 // ...unless nodes A and B share the element's last edge
914 if ((local_indexA == 0) && (local_indexB == p_element->GetNumNodes()-1))
915 {
916 index = local_indexB;
917 }
918 }
919 else if ((local_indexB == 0) && (local_indexA == p_element->GetNumNodes()-1))
920 {
921 // ...otherwise use node B's index, unless nodes A and B share the element's last edge
922 index = local_indexA;
923 }
924
925 // Add new node to this element
926 this->GetElement(*iter)->AddNode(p_new_node, index);
927 }
928}
929
930template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
931void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::RemoveDeletedNodesAndElements(VertexElementMap& rElementMap)
932{
933 // Make sure the map is big enough. Each entry will be set in the loop below.
934 rElementMap.Resize(this->GetNumAllElements());
935 // Remove any elements that have been marked for deletion and store all other elements in a temporary structure
936 std::vector<VertexElement<ELEMENT_DIM, SPACE_DIM>*> live_elements;
937 for (unsigned i=0; i<this->mElements.size(); i++)
938 {
939 if (this->mElements[i]->IsDeleted())
940 {
941 delete this->mElements[i];
942 rElementMap.SetDeleted(i);
943 }
944 else
945 {
946 live_elements.push_back(this->mElements[i]);
947 rElementMap.SetNewIndex(i, (unsigned)(live_elements.size()-1));
948 }
949 }
950
951 // Sanity check
952 assert(mDeletedElementIndices.size() == this->mElements.size() - live_elements.size());
953
954 // Repopulate the elements vector and reset the list of deleted element indices
955 mDeletedElementIndices.clear();
956 this->mElements = live_elements;
957
958 // Finally, reset the element indices to run from zero
959 for (unsigned i=0; i<this->mElements.size(); i++)
960 {
961 this->mElements[i]->ResetIndex(i);
962 }
963
964 // Remove deleted nodes
965 RemoveDeletedNodes();
966}
967
968template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
969void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::RemoveDeletedNodes()
970{
971 // Remove any nodes that have been marked for deletion and store all other nodes in a temporary structure
972 // Also mark edges associated with the deleted nodes
973 std::vector<Node<SPACE_DIM>*> live_nodes;
974 for (unsigned i=0; i<this->mNodes.size(); i++)
975 {
976 if (this->mNodes[i]->IsDeleted())
977 {
978 delete this->mNodes[i];
979 }
980 else
981 {
982 live_nodes.push_back(this->mNodes[i]);
983 }
984 }
985
986 // Sanity check
987 assert(mDeletedNodeIndices.size() == this->mNodes.size() - live_nodes.size());
988 // Repopulate the nodes vector and reset the list of deleted node indices
989 this->mNodes = live_nodes;
990 mDeletedNodeIndices.clear();
991
992 // Finally, reset the node indices to run from zero
993 for (unsigned i=0; i<this->mNodes.size(); i++)
994 {
995 this->mNodes[i]->SetIndex(i);
996 }
997
998 // Remove deleted edges
999 this->mEdgeHelper.RemoveDeletedEdges();
1000}
1001
1002template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1003void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ReMesh([[maybe_unused]] VertexElementMap& rElementMap)
1004{
1005 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1006 {
1007 // Make sure the map is big enough
1008 rElementMap.Resize(this->GetNumAllElements());
1009
1010 /*
1011 * To begin the remeshing process, we do not need to call Clear() and remove all current data,
1012 * since cell birth, rearrangement and death result only in local remeshing of a vertex-based
1013 * mesh. Instead, we just remove any deleted elements and nodes.
1014 */
1015 RemoveDeletedNodesAndElements(rElementMap);
1016 bool recheck_mesh = true;
1017 while (recheck_mesh == true)
1018 {
1019 // We check for any short edges and perform swaps if necessary and possible.
1020 recheck_mesh = CheckForSwapsFromShortEdges();
1021 }
1022
1023 // Check for element intersections
1024 recheck_mesh = true;
1025 while (recheck_mesh == true)
1026 {
1027 // Check mesh for intersections, and perform T3 swaps where required
1028 recheck_mesh = CheckForIntersections();
1029 }
1030
1031 RemoveDeletedNodes();
1032
1033 /*
1034 * This is handled in a separate method to allow child classes to implement additional ReMeshing functionality
1035 * (see #2664).
1036 */
1037 this->CheckForRosettes();
1038 }
1039 else if constexpr (ELEMENT_DIM == 3 && SPACE_DIM == 3)
1040 {
1041 EXCEPTION("Remeshing has not been implemented in 3D (see Trac tickets #827, #860, #1422)\n"); // LCOV_EXCL_LINE
1043 }
1044 else
1045 {
1046 NEVER_REACHED;
1047 }
1048}
1049
1050template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1051void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ReMesh()
1052{
1053 VertexElementMap map(GetNumElements());
1054 ReMesh(map);
1055}
1056
1057template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1058bool MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::CheckForSwapsFromShortEdges()
1059{
1060 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1061 {
1062 // Loop over elements to check for T1 swaps
1063 for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
1064 elem_iter != this->GetElementIteratorEnd();
1065 ++elem_iter)
1066 {
1068
1069 unsigned num_nodes = elem_iter->GetNumNodes();
1070 assert(num_nodes > 0);
1071
1072 // Loop over the nodes contained in this element
1073 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1074 {
1075 // Find locations of the current node and anticlockwise node
1076 Node<SPACE_DIM>* p_current_node = elem_iter->GetNode(local_index);
1077 unsigned local_index_plus_one = (local_index + 1) % num_nodes;
1078 Node<SPACE_DIM>* p_anticlockwise_node = elem_iter->GetNode(local_index_plus_one);
1079
1080 // Find distance between nodes
1081 double distance_between_nodes = this->GetDistanceBetweenNodes(p_current_node->GetIndex(), p_anticlockwise_node->GetIndex());
1082
1083 // If the nodes are too close together...
1084 if (distance_between_nodes < mCellRearrangementThreshold)
1085 {
1086 // ...then check if any triangular elements are shared by these nodes...
1087 std::set<unsigned> elements_of_node_a = p_current_node->rGetContainingElementIndices();
1088 std::set<unsigned> elements_of_node_b = p_anticlockwise_node->rGetContainingElementIndices();
1089
1090 std::set<unsigned> shared_elements;
1091 std::set_intersection(elements_of_node_a.begin(), elements_of_node_a.end(),
1092 elements_of_node_b.begin(), elements_of_node_b.end(),
1093 std::inserter(shared_elements, shared_elements.begin()));
1094
1095 bool both_nodes_share_triangular_element = false;
1096 for (std::set<unsigned>::const_iterator it = shared_elements.begin();
1097 it != shared_elements.end();
1098 ++it)
1099 {
1100 if (this->GetElement(*it)->GetNumNodes() <= 3)
1101 {
1102 both_nodes_share_triangular_element = true;
1103 break;
1104 }
1105 }
1106
1107 // ...and if none are, then perform the required type of swap and halt the search, returning true
1108 if (!both_nodes_share_triangular_element)
1109 {
1110 IdentifySwapType(p_current_node, p_anticlockwise_node);
1111 return true;
1112 }
1113 }
1114 }
1115 }
1116
1117 return false;
1118 }
1119 else
1120 {
1121 NEVER_REACHED;
1122 }
1123}
1124
1125template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1126bool MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::CheckForT2Swaps(VertexElementMap& rElementMap)
1127{
1128 // Loop over elements to check for T2 swaps
1129 for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
1130 elem_iter != this->GetElementIteratorEnd();
1131 ++elem_iter)
1132 {
1133 // If this element is triangular...
1134 if (elem_iter->GetNumNodes() == 3)
1135 {
1136 // ...and smaller than the threshold area...
1137 if (this->GetVolumeOfElement(elem_iter->GetIndex()) < GetT2Threshold())
1138 {
1139 // ...then perform a T2 swap and break out of the loop
1140 PerformT2Swap(*elem_iter);
1142 rElementMap.SetDeleted(elem_iter->GetIndex());
1143 return true;
1144 }
1145 }
1146 }
1147 return false;
1148}
1149
1150template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1151bool MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::CheckForIntersections()
1152{
1153 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1154 {
1155 // If checking for internal intersections, then check that no nodes have overlapped any elements...
1156 if (mCheckForInternalIntersections)
1157 {
1158 // First neighbours are elements that contain the node. Second
1159 // neighbours are elements that share a cell-cell boundary with first
1160 // neighbours, but do not contain the node. Nodes can only intersect
1161 // second neighbours.
1162 for (auto node_iter = this->GetNodeIteratorBegin();
1163 node_iter != this->GetNodeIteratorEnd();
1164 ++node_iter)
1165 {
1166 assert(!(node_iter->IsDeleted()));
1167
1168 // Get all nodes of first neighbours
1169 std::set<unsigned> first_neighbour_node_indices;
1170
1171 auto first_neighbour_indices = node_iter->rGetContainingElementIndices();
1172
1173 for (auto elem_iter = first_neighbour_indices.begin();
1174 elem_iter != first_neighbour_indices.end();
1175 ++elem_iter)
1176 {
1177 auto p_element = this->GetElement(*elem_iter);
1178
1179 for (auto local_node_index = 0u;
1180 local_node_index < p_element->GetNumNodes();
1181 ++local_node_index)
1182 {
1183 first_neighbour_node_indices.insert(p_element->GetNodeGlobalIndex(local_node_index));
1184 }
1185 }
1186
1187 // Get all first and second neighbours
1188 std::set<unsigned> all_neighbours;
1189
1190 for (auto second_node_iter = first_neighbour_node_indices.begin();
1191 second_node_iter != first_neighbour_node_indices.end();
1192 ++second_node_iter)
1193 {
1194 auto containing_element_indices = this->GetNode(*second_node_iter)->rGetContainingElementIndices();
1195 all_neighbours.insert(containing_element_indices.begin(),
1196 containing_element_indices.end());
1197 }
1198
1199 // Second neighbours are the difference between all neighbours and
1200 // first neighbours
1201 std::set<unsigned> second_neighbour_indices;
1202 std::set_difference(
1203 all_neighbours.begin(), all_neighbours.end(),
1204 first_neighbour_indices.begin(), first_neighbour_indices.end(),
1205 std::inserter(second_neighbour_indices, second_neighbour_indices.begin()));
1206
1207 // Loop over second neighbours only
1208 for (auto elem_iter = second_neighbour_indices.begin();
1209 elem_iter != second_neighbour_indices.end();
1210 ++elem_iter)
1211 {
1212 unsigned elem_index = *elem_iter;
1213
1214 // Node should not be part of this element
1215 assert(node_iter->rGetContainingElementIndices().count(elem_index) == 0);
1216
1217 if (this->ElementIncludesPoint(node_iter->rGetLocation(), elem_index))
1218 {
1219 PerformIntersectionSwap(&(*node_iter), elem_index);
1220 return true;
1221 }
1222 }
1223 }
1224 }
1225
1226 if (mCheckForT3Swaps)
1227 {
1228 // If checking for T3 swaps, check that no boundary nodes have overlapped any boundary elements
1229 // First: find all boundary element and calculate their centroid only once
1230 std::vector<unsigned> boundary_element_indices;
1231 std::vector<c_vector<double, SPACE_DIM> > boundary_element_centroids;
1232 for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
1233 elem_iter != this->GetElementIteratorEnd();
1234 ++elem_iter)
1235 {
1236 if (elem_iter->IsElementOnBoundary())
1237 {
1238 unsigned element_index = elem_iter->GetIndex();
1239 boundary_element_indices.push_back(element_index);
1240 // should be a map but I am too lazy to look up the syntax
1241 boundary_element_centroids.push_back(this->GetCentroidOfElement(element_index));
1242 }
1243 }
1244
1245 // Second: Check intersections only for those nodes and elements within
1246 // mDistanceForT3SwapChecking within each other (node<-->element centroid)
1247 for (typename AbstractMesh<ELEMENT_DIM, SPACE_DIM>::NodeIterator node_iter = this->GetNodeIteratorBegin();
1248 node_iter != this->GetNodeIteratorEnd();
1249 ++node_iter)
1250 {
1251 if (node_iter->IsBoundaryNode())
1252 {
1253 assert(!(node_iter->IsDeleted()));
1254
1255 // index in boundary_element_centroids and boundary_element_indices
1256 unsigned boundary_element_index = 0;
1257 for (std::vector<unsigned>::iterator elem_iter = boundary_element_indices.begin();
1258 elem_iter != boundary_element_indices.end();
1259 ++elem_iter)
1260 {
1261 // Check that the node is not part of this element
1262 if (node_iter->rGetContainingElementIndices().count(*elem_iter) == 0)
1263 {
1264 c_vector<double, SPACE_DIM> node_location = node_iter->rGetLocation();
1265 c_vector<double, SPACE_DIM> element_centroid = boundary_element_centroids[boundary_element_index];
1266 double node_element_distance = norm_2(this->GetVectorFromAtoB(node_location, element_centroid));
1267
1268 if (node_element_distance < mDistanceForT3SwapChecking)
1269 {
1270 if (this->ElementIncludesPoint(node_iter->rGetLocation(), *elem_iter))
1271 {
1272 this->PerformT3Swap(&(*node_iter), *elem_iter);
1273 return true;
1274 }
1275 }
1276 }
1277 // increment the boundary element index
1278 boundary_element_index += 1u;
1279 }
1280 }
1281 }
1282 }
1283 return false;
1284 }
1285 else
1286 {
1287 NEVER_REACHED;
1288 }
1289}
1290
1291template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1292void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::IdentifySwapType(
1293 [[maybe_unused]] Node<SPACE_DIM>* pNodeA, [[maybe_unused]] Node<SPACE_DIM>* pNodeB)
1294{
1295 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1296 {
1297 // Find the sets of elements containing nodes A and B
1298 std::set<unsigned> nodeA_elem_indices = pNodeA->rGetContainingElementIndices();
1299 std::set<unsigned> nodeB_elem_indices = pNodeB->rGetContainingElementIndices();
1300
1301 // Form the set union
1302 std::set<unsigned> all_indices, temp_union_set;
1303 std::set_union(nodeA_elem_indices.begin(), nodeA_elem_indices.end(),
1304 nodeB_elem_indices.begin(), nodeB_elem_indices.end(),
1305 std::inserter(temp_union_set, temp_union_set.begin()));
1306 all_indices.swap(temp_union_set); // temp_set will be deleted, all_indices now contains all the indices of elements
1307 // that touch the potentially swapping nodes
1308
1309 if ((nodeA_elem_indices.size() > 3) || (nodeB_elem_indices.size() > 3))
1310 {
1311 /*
1312 * Looks like
1313 *
1314 * \
1315 * \ A B
1316 * ---o---o---
1317 * /
1318 * /
1319 *
1320 */
1321
1322 /*
1323 * This case is handled in a separate method to allow child classes to implement different
1324 * functionality for high-order-junction remodelling events (see #2664).
1325 */
1326 this->HandleHighOrderJunctions(pNodeA, pNodeB);
1327 }
1328 else // each node is contained in at most three elements
1329 {
1330 switch (all_indices.size())
1331 {
1332 case 1:
1333 {
1334 /*
1335 * Each node is contained in a single element, so the nodes must lie on the boundary
1336 * of the mesh, as shown below. In this case, we merge the nodes and tidy up node
1337 * indices through calls to PerformNodeMerge() and RemoveDeletedNodes().
1338 *
1339 * A B
1340 * ---o---o---
1341 */
1342 assert(pNodeA->IsBoundaryNode());
1343 assert(pNodeB->IsBoundaryNode());
1344 PerformNodeMerge(pNodeA, pNodeB);
1345 RemoveDeletedNodes();
1346 break;
1347 }
1348 case 2:
1349 {
1350 if (nodeA_elem_indices.size() == 2 && nodeB_elem_indices.size() == 2)
1351 {
1352 if (pNodeA->IsBoundaryNode() && pNodeB->IsBoundaryNode())
1353 {
1354 /*
1355 * The node configuration is as shown below, with voids on either side. In this case
1356 * we perform a T1 swap, which separates the elements.
1357 *
1358 * \ /
1359 * \ / Node A
1360 * (1) | (2) (element number in brackets)
1361 * / \ Node B
1362 * / \
1363 */
1364 PerformT1Swap(pNodeA, pNodeB, all_indices);
1365 }
1366 else if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
1367 {
1368 /*
1369 * The node configuration is as shown below, with a void on one side. We should not
1370 * be able to reach this case at present, since we allow only for three-way junctions
1371 * or boundaries, so we throw an exception.
1372 *
1373 * \ /
1374 * \ / Node A
1375 * (1) | (2) (element number in brackets)
1376 * x Node B
1377 * |
1378 */
1379 EXCEPTION("There is a non-boundary node contained only in two elements; something has gone wrong.");
1380 }
1381 else
1382 {
1383 /*
1384 * Each node is contained in two elements, so the nodes lie on an internal edge, as shown below.
1385 * We should not be able to reach this case at present, since we allow only for three-way junctions
1386 * or boundaries, so we throw an exception.
1387 *
1388 * A B
1389 * ---o---o---
1390 */
1391 EXCEPTION("There are non-boundary nodes contained only in two elements; something has gone wrong.");
1392 }
1393 } // from [if (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==2)]
1394 else
1395 {
1396 /*
1397 * The node configuration either looks like that shown below. In this case, we merge the nodes
1398 * and tidy up node indices through calls to PerformNodeMerge() and RemoveDeletedNodes().
1399 *
1400 * Outside
1401 * /
1402 * --o--o (2)
1403 * (1) \
1404 *
1405 * Or its an internal triangular void near the boundary. Like this
1406 *
1407 * x
1408 * |
1409 * o-o
1410 * |/ Where the area inside the triangle is a void and the horizontal edge
1411 * o is the short edge all the nodes are therefore boundary nodes.
1412 * | Here we remove the void and merge all the nodes with one of the nodes at the ends
1413 * x
1414 *
1415 * Or it's a more complicated situation that's currently not identified.
1416 *
1417 * So we search to differentiate these cases
1418 */
1419 assert((nodeA_elem_indices.size() == 1 && nodeB_elem_indices.size() == 2)
1420 || (nodeA_elem_indices.size() == 2 && nodeB_elem_indices.size() == 1));
1421
1422 // Node one is the potential convex node. Node 2 is the one in both elements.
1423 Node<SPACE_DIM>* p_node_1 = pNodeA;
1424 Node<SPACE_DIM>* p_node_2 = pNodeB;
1425 std::set<unsigned> node1_elem_indices = nodeA_elem_indices;
1426 std::set<unsigned> node2_elem_indices = nodeB_elem_indices;
1427
1428 if (nodeB_elem_indices.size() == 1)
1429 {
1430 p_node_1 = pNodeB;
1431 node1_elem_indices = nodeB_elem_indices;
1432 p_node_2 = pNodeA;
1433 node2_elem_indices = nodeA_elem_indices;
1434 }
1435 assert(node1_elem_indices.size() == 1);
1436 unsigned unique_element_index = *node1_elem_indices.begin();
1437
1438 node2_elem_indices.erase(unique_element_index);
1439 assert(node2_elem_indices.size() == 1);
1440 unsigned common_element_index = *node2_elem_indices.begin();
1441
1442 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_uniqiue_element = this->mElements[unique_element_index];
1443 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_common_element = this->mElements[common_element_index];
1444
1445 unsigned local_index_1 = p_uniqiue_element->GetNodeLocalIndex(p_node_1->GetIndex());
1446 unsigned next_node_1 = p_uniqiue_element->GetNodeGlobalIndex((local_index_1 + 1) % (p_uniqiue_element->GetNumNodes()));
1447 unsigned previous_node_1 = p_uniqiue_element->GetNodeGlobalIndex(
1448 (local_index_1 + p_uniqiue_element->GetNumNodes() - 1) % (p_uniqiue_element->GetNumNodes()));
1449 unsigned previous_previous_node_1 = p_uniqiue_element->GetNodeGlobalIndex(
1450 (local_index_1 + p_uniqiue_element->GetNumNodes() - 2) % (p_uniqiue_element->GetNumNodes()));
1451 unsigned local_index_2 = p_common_element->GetNodeLocalIndex(p_node_2->GetIndex());
1452 unsigned next_node_2 = p_common_element->GetNodeGlobalIndex(
1453 (local_index_2 + 1) % (p_common_element->GetNumNodes()));
1454 unsigned previous_node_2 = p_common_element->GetNodeGlobalIndex(
1455 (local_index_2 + p_common_element->GetNumNodes() - 1) % (p_common_element->GetNumNodes()));
1456
1457 if (next_node_1 == previous_node_2 || next_node_2 == previous_node_1)
1458 {
1459 /*
1460 * Here we have an internal triangular void on an internal edge. Can happen when a void is shrinking.
1461 *
1462 * x
1463 * |
1464 * o-o
1465 * |/ Where the area inside the triangle is a void and the horizontal edge
1466 * o is the short edge all the nodes are therefore boundary nodes.
1467 * | Here we remove the void and merge all the nodes with one of the nodes at
1468 * x the adjoining edges (x's)
1469 *
1470 *
1471 */
1472
1473 // First remove void
1474 // Find all three nodes in void
1475 Node<SPACE_DIM>* p_node_C = this->mNodes[next_node_2]; // The other node in the triangular void
1476 if (next_node_1 == previous_node_2)
1477 {
1478 p_node_C = this->mNodes[next_node_1];
1479 }
1480
1481 /*
1482 * In two steps, merge nodes A, B and C into a single node. This is implemented in such a way that
1483 * the ordering of their indices does not matter.
1484 */
1485
1486 PerformNodeMerge(pNodeA, pNodeB);
1487
1488 Node<SPACE_DIM>* p_merged_node = pNodeB;
1489
1490 if (pNodeB->IsDeleted())
1491 {
1492 p_merged_node = pNodeA;
1493 }
1494
1495 PerformNodeMerge(p_node_C, p_merged_node);
1496
1497 if (p_merged_node->IsDeleted())
1498 {
1499 p_merged_node = p_node_C;
1500 }
1501
1502 // Tag remaining node as non-boundary
1503 p_merged_node->SetAsBoundaryNode(false);
1504
1505 // Now merge this node with one of the nearest vertices keeping that vertices location.
1506 Node<SPACE_DIM>* p_end_node; // The neighbouring vertex to merge to
1507
1508 std::set<unsigned> previous_previous_elem_indices = this->mNodes[previous_previous_node_1]->rGetContainingElementIndices();
1509
1510 std::set<unsigned> shared_elements;
1511 std::set_intersection(all_indices.begin(),
1512 all_indices.end(),
1513 previous_previous_elem_indices.begin(),
1514 previous_previous_elem_indices.end(),
1515 std::inserter(shared_elements, shared_elements.begin()));
1516
1517 assert(shared_elements.size() < 3);
1518
1519 if (shared_elements.size() == 2)
1520 {
1521 // This neighbouring node is in the same 2 elements so treat this as end node
1522 p_end_node = this->mNodes[previous_previous_node_1];
1523 }
1524 else
1525 {
1526 /*
1527 * If this trips then the adjacent node isn't in both elements and we currently dont deal with this.
1528 * See #3080
1529 */
1530 NEVER_REACHED;
1531 }
1532
1533 p_merged_node->rGetModifiableLocation() = p_end_node->rGetLocation();
1534 // We perform the merge in this order so the first node is kept as this has correct boundary information.
1535 PerformNodeMerge(p_end_node, p_merged_node);
1536
1537 // Remove the deleted nodes and re-index
1538 RemoveDeletedNodes();
1539 }
1540 else
1541 {
1542 /*
1543 * The node configuration looks like that shown below. In this case, we merge the nodes
1544 * and tidy up node indices through calls to PerformNodeMerge() and RemoveDeletedNodes().
1545 *
1546 * Outside
1547 * /
1548 * --o--o (2)
1549 * (1) \
1550 */
1551 PerformNodeMerge(pNodeA, pNodeB);
1552 RemoveDeletedNodes();
1553 }
1554 }
1555 break;
1556 }
1557 case 3:
1558 {
1559 if (nodeA_elem_indices.size() == 1 || nodeB_elem_indices.size() == 1)
1560 {
1561 /*
1562 * One node is contained in one element and the other node is contained in three elements.
1563 * We should not be able to reach this case at present, since we allow each boundary node
1564 * to be contained in at most two elements, so we throw an exception.
1565 *
1566 * A B
1567 *
1568 * empty /
1569 * / (3)
1570 * ---o---o----- (element number in brackets)
1571 * (1) \ (2)
1572 * \
1573 */
1574 assert(pNodeA->IsBoundaryNode());
1575 assert(pNodeB->IsBoundaryNode());
1576
1577 EXCEPTION("There is a boundary node contained in three elements something has gone wrong.");
1578 }
1579 else if (nodeA_elem_indices.size() == 2 && nodeB_elem_indices.size() == 2)
1580 {
1581 // The short edge must be at the boundary. We need to check whether this edge is
1582 // adjacent to a triangular void before we swap. If it is a triangular void, we perform a T2-type swap.
1583 // If not, then we perform a normal T1 swap. I.e. in detail we need to check whether the
1584 // element in nodeA_elem_indices which is not in nodeB_elem_indices contains a shared node
1585 // with the element in nodeB_elem_indices which is not in nodeA_elem_indices.
1586
1587 std::set<unsigned> element_A_not_B, temp_set;
1588 std::set_difference(all_indices.begin(), all_indices.end(), nodeB_elem_indices.begin(),
1589 nodeB_elem_indices.end(), std::inserter(temp_set, temp_set.begin()));
1590 element_A_not_B.swap(temp_set);
1591
1592 // There must be only one such element
1593 assert(element_A_not_B.size() == 1);
1594
1595 std::set<unsigned> element_B_not_A;
1596 std::set_difference(all_indices.begin(), all_indices.end(), nodeA_elem_indices.begin(),
1597 nodeA_elem_indices.end(), std::inserter(temp_set, temp_set.begin()));
1598 element_B_not_A.swap(temp_set);
1599
1600 // There must be only one such element
1601 assert(element_B_not_A.size() == 1);
1602
1603 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_A_not_B = this->mElements[*element_A_not_B.begin()];
1604 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_B_not_A = this->mElements[*element_B_not_A.begin()];
1605
1606 unsigned local_index_1 = p_element_A_not_B->GetNodeLocalIndex(pNodeA->GetIndex());
1607 unsigned next_node_1 = p_element_A_not_B->GetNodeGlobalIndex((local_index_1 + 1) % (p_element_A_not_B->GetNumNodes()));
1608 unsigned previous_node_1 = p_element_A_not_B->GetNodeGlobalIndex(
1609 (local_index_1 + p_element_A_not_B->GetNumNodes() - 1) % (p_element_A_not_B->GetNumNodes()));
1610 unsigned local_index_2 = p_element_B_not_A->GetNodeLocalIndex(pNodeB->GetIndex());
1611 unsigned next_node_2 = p_element_B_not_A->GetNodeGlobalIndex(
1612 (local_index_2 + 1) % (p_element_B_not_A->GetNumNodes()));
1613 unsigned previous_node_2 = p_element_B_not_A->GetNodeGlobalIndex(
1614 (local_index_2 + p_element_B_not_A->GetNumNodes() - 1) % (p_element_B_not_A->GetNumNodes()));
1615
1616 if (next_node_1 == previous_node_2 || next_node_2 == previous_node_1)
1617 {
1618 /*
1619 * The node configuration looks like that shown below, and both nodes must be on the boundary.
1620 * In this case we remove the void through a call to PerformVoidRemoval().
1621 *
1622 * A C B A B
1623 * /\ \ /
1624 * /v \ \ (1) /
1625 * (3)o----o (1) or (2) o----o (3) (element number in brackets, v is a void)
1626 * / (2) \ \v /
1627 * / \ \/
1628 * C
1629 */
1630 assert(pNodeA->IsBoundaryNode());
1631 assert(pNodeB->IsBoundaryNode());
1632
1633 // Get the third node in the triangular void
1634
1635 unsigned nodeC_index;
1636 if (next_node_1 == previous_node_2 && next_node_2 != previous_node_1)
1637 {
1638 nodeC_index = next_node_1;
1639 }
1640 else if (next_node_2 == previous_node_1 && next_node_1 != previous_node_2)
1641 {
1642 nodeC_index = next_node_2;
1643 }
1644 else
1645 {
1646 assert(next_node_1 == previous_node_2 && next_node_2 == previous_node_1);
1654 EXCEPTION("Triangular element next to triangular void, not implemented yet.");
1655 }
1656
1657 if (p_element_A_not_B->GetNumNodes() == 3u || p_element_B_not_A->GetNumNodes() == 3u)
1658 {
1666 EXCEPTION("Triangular element next to triangular void, not implemented yet.");
1667 }
1668 PerformVoidRemoval(pNodeA, pNodeB, this->mNodes[nodeC_index]);
1669 }
1670 else
1671 {
1672 /*
1673 * The node configuration looks like that below, and both nodes must lie on the boundary.
1674 * In this case we perform a T1 swap.
1675 *
1676 * A B A B
1677 * \ empty/ \ /
1678 * \ / \‍(1) /
1679 * (3) o--o (1) or (2) o--o (3) (element number in brackets)
1680 * / (2)\ / \
1681 * / \ /empty \
1682 */
1683 assert(pNodeA->IsBoundaryNode());
1684 assert(pNodeB->IsBoundaryNode());
1685 PerformT1Swap(pNodeA, pNodeB, all_indices);
1686 }
1687 } // from else if (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==2)
1688 else
1689 {
1690 // In this case, one node must be contained in two elements and the other in three elements.
1691 assert((nodeA_elem_indices.size() == 2 && nodeB_elem_indices.size() == 3)
1692 || (nodeA_elem_indices.size() == 3 && nodeB_elem_indices.size() == 2));
1693
1694 // They can't both be boundary nodes
1695 assert(!(pNodeA->IsBoundaryNode() && pNodeB->IsBoundaryNode()));
1696
1697 if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
1698 {
1699 /*
1700 * The node configuration looks like that shown below. We perform a T1 swap in this case.
1701 *
1702 * A B A B
1703 * \ / \ /
1704 * \ (1)/ \‍(1) /
1705 * (3) o--o (empty) or (empty) o--o (3) (element number in brackets)
1706 * / (2)\ /(2) \
1707 * / \ / \
1708 */
1709 PerformT1Swap(pNodeA, pNodeB, all_indices);
1710 }
1711 else
1712 {
1713 /*
1714 * The node configuration looks like that shown below. We should not be able to reach this case
1715 * at present, since we allow only for three-way junctions or boundaries, so we throw an exception.
1716 *
1717 * A B A B
1718 * \ /
1719 * \ (1) (1) /
1720 * (3) o--o--- or ---o--o (3) (element number in brackets)
1721 * / (2) (2) \
1722 * / \
1723 */
1724 EXCEPTION("There are non-boundary nodes contained only in two elements; something has gone wrong.");
1725 }
1726 }
1727 break;
1728 }
1729 case 4:
1730 {
1731 /*
1732 * The node configuration looks like that shown below. We perform a T1 swap in this case.
1733 *
1734 * \‍(1)/
1735 * \ / Node A
1736 * (2) | (4) (element number in brackets)
1737 * / \ Node B
1738 * /(3)\
1739 */
1740
1741 /*
1742 * This case is handled in a separate method to allow child classes to implement different
1743 * functionality for junction remodelling events (see #2664).
1744 */
1745 if (mProtorosetteFormationProbability > RandomNumberGenerator::Instance()->ranf())
1746 {
1747 this->PerformNodeMerge(pNodeA, pNodeB);
1748 this->RemoveDeletedNodes();
1749 }
1750 else
1751 {
1752 this->PerformT1Swap(pNodeA, pNodeB, all_indices);
1753 }
1754 break;
1755 }
1756 default:
1757 // This can't happen
1758 NEVER_REACHED;
1759 }
1760 }
1761 }
1762 else
1763 {
1764 NEVER_REACHED;
1765 }
1766}
1767
1768template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1769void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformNodeMerge(Node<SPACE_DIM>* pNodeA, Node<SPACE_DIM>* pNodeB)
1770{
1771 // Find the sets of elements containing each of the nodes, sorted by index
1772 std::set<unsigned> nodeA_elem_indices = pNodeA->rGetContainingElementIndices();
1773 std::set<unsigned> nodeB_elem_indices = pNodeB->rGetContainingElementIndices();
1774
1775 // Move node A to the mid-point
1776 pNodeA->rGetModifiableLocation() += 0.5 * this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation());
1777
1778 // Update the elements previously containing node B to contain node A
1779 unsigned node_B_index = pNodeB->GetIndex();
1780 // For rebuilding edges affected elements
1781 std::set<unsigned> rebuilt_elements;
1782 for (std::set<unsigned>::const_iterator it = nodeB_elem_indices.begin(); it != nodeB_elem_indices.end(); ++it)
1783 {
1784 // Find the local index of node B in this element
1785 unsigned node_B_local_index = this->mElements[*it]->GetNodeLocalIndex(node_B_index);
1786 assert(node_B_local_index < UINT_MAX); // this element contains node B
1787
1788 /*
1789 * If this element already contains node A, then just remove node B.
1790 * Otherwise replace it with node A in the element and remove it from mNodes.
1791 */
1792 if (nodeA_elem_indices.count(*it) != 0)
1793 {
1794 std::vector<unsigned> edgeIds;
1795 for (unsigned i = 0; i < this->mElements[*it]->GetNumEdges(); i++)
1796 {
1797 edgeIds.push_back(this->mElements[*it]->GetEdge(i)->GetIndex());
1798 }
1799 const unsigned node_A_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeA->GetIndex());
1800 this->mElements[*it]->DeleteNode(node_B_local_index);
1801
1802 if (mTrackMeshOperations)
1803 {
1804 mOperationRecorder.RecordNodeMergeOperation(edgeIds, this->mElements[*it],
1805 std::pair<unsigned, unsigned>(node_A_local_index, node_B_local_index));
1806 }
1807 }
1808 else
1809 {
1810 // Replace node B with node A in this element
1811 this->mElements[*it]->UpdateNode(node_B_local_index, pNodeA);
1812 }
1813 rebuilt_elements.insert(this->mElements[*it]->GetIndex());
1814 this->mElements[*it]->RebuildEdges();
1815 }
1816
1817 assert(!(this->mNodes[node_B_index]->IsDeleted()));
1818 this->mNodes[node_B_index]->MarkAsDeleted();
1819 mDeletedNodeIndices.push_back(node_B_index);
1820 for (unsigned i:rebuilt_elements)
1821 this->GetElement(i)->RebuildEdges();
1822}
1823
1824template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1825void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformT1Swap(
1826 [[maybe_unused]] Node<SPACE_DIM>* pNodeA,
1827 [[maybe_unused]] Node<SPACE_DIM>* pNodeB,
1828 [[maybe_unused]] std::set<unsigned>& rElementsContainingNodes)
1829{
1830 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
1831 {
1832 // First compute and store the location of the T1 swap, which is at the midpoint of nodes A and B
1833 double distance_between_nodes_CD = mCellRearrangementRatio * mCellRearrangementThreshold;
1834
1835 c_vector<double, SPACE_DIM> nodeA_location = pNodeA->rGetLocation();
1836 c_vector<double, SPACE_DIM> nodeB_location = pNodeB->rGetLocation();
1837 c_vector<double, SPACE_DIM> vector_AB = this->GetVectorFromAtoB(nodeA_location, nodeB_location);
1838
1839 double distance_AB = norm_2(vector_AB);
1840 if (distance_AB < 1e-10)
1841 {
1842 EXCEPTION("Nodes are too close together, this shouldn't happen");
1843 }
1844
1845 /*
1846 * Compute the locations of two new nodes C, D, placed on either side of the
1847 * edge E_old formed by nodes A and B, such that the edge E_new formed by the
1848 * new nodes is the perpendicular bisector of E_old, with |E_new| 'just larger'
1849 * (mCellRearrangementRatio) than mThresholdDistance.
1850 *
1851 * We implement the following changes to the mesh:
1852 *
1853 * The element whose index was in nodeA_elem_indices but not nodeB_elem_indices,
1854 * and the element whose index was in nodeB_elem_indices but not nodeA_elem_indices,
1855 * should now both contain nodes A and B.
1856 *
1857 * The element whose index was in nodeA_elem_indices and nodeB_elem_indices, and which
1858 * node C lies inside, should now only contain node A.
1859 *
1860 * The element whose index was in nodeA_elem_indices and nodeB_elem_indices, and which
1861 * node D lies inside, should now only contain node B.
1862 *
1863 * Iterate over all elements involved and identify which element they are
1864 * in the diagram then update the nodes as necessary.
1865 *
1866 * \‍(1)/
1867 * \ / Node A
1868 * (2) | (4) elements in brackets
1869 * / \ Node B
1870 * /(3)\
1871 */
1872
1873 // Move nodes A and B to C and D respectively
1874 c_vector<double, SPACE_DIM> vector_CD;
1875 vector_CD(0) = -vector_AB(1) * distance_between_nodes_CD / distance_AB;
1876 vector_CD(1) = vector_AB(0) * distance_between_nodes_CD / distance_AB;
1877
1878 // Record T1Swap
1879 T1SwapInfo<SPACE_DIM> swap_info;
1880 swap_info.mLocation = nodeA_location + 0.5 * vector_AB;
1881 swap_info.mPreSwapEdge = vector_AB;
1882 swap_info.mPostSwapEdge = vector_CD;
1883 mOperationRecorder.RecordT1Swap(swap_info);
1884
1885 c_vector<double, SPACE_DIM> nodeC_location = nodeA_location + 0.5 * vector_AB - 0.5 * vector_CD;
1886 c_vector<double, SPACE_DIM> nodeD_location = nodeC_location + vector_CD;
1887
1888 pNodeA->rGetModifiableLocation() = nodeC_location;
1889 pNodeB->rGetModifiableLocation() = nodeD_location;
1890
1891 // Find the sets of elements containing nodes A and B
1892 std::set<unsigned> nodeA_elem_indices = pNodeA->rGetContainingElementIndices();
1893 std::set<unsigned> nodeB_elem_indices = pNodeB->rGetContainingElementIndices();
1894
1895 // For rebuilding affected elements
1896 std::set<unsigned> rebuilt_elements;
1897 for (std::set<unsigned>::const_iterator it = rElementsContainingNodes.begin();
1898 it != rElementsContainingNodes.end();
1899 ++it)
1900 {
1901 std::vector<unsigned> old_ids;
1902 for (unsigned i=0; i<this->mElements[*it]->GetNumEdges(); ++i)
1903 {
1904 old_ids.push_back(this->mElements[*it]->GetEdge(i)->GetIndex());
1905 }
1906
1907 // If, as in element 3 above, this element does not contain node A (now C)...
1908 if (nodeA_elem_indices.find(*it) == nodeA_elem_indices.end())
1909 {
1910 // ...then add it to the element just after node B (now D), going anticlockwise
1911 unsigned nodeB_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeB->GetIndex());
1912 assert(nodeB_local_index < UINT_MAX);
1913
1914 std::vector<unsigned> old_ids;
1915 for (unsigned i=0; i<this->mElements[*it]->GetNumEdges(); ++i)
1916 {
1917 old_ids.push_back(this->mElements[*it]->GetEdge(i)->GetIndex());
1918 }
1919
1920 this->mElements[*it]->AddNode(pNodeA, nodeB_local_index);
1921 if (mTrackMeshOperations)
1922 {
1923 mOperationRecorder.RecordNewEdgeOperation(this->mElements[*it], nodeB_local_index);
1924 }
1925 }
1926 else if (nodeB_elem_indices.find(*it) == nodeB_elem_indices.end())
1927 {
1928 // Do similarly if the element does not contain node B (now D), as in element 1 above
1929 unsigned nodeA_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeA->GetIndex());
1930 assert(nodeA_local_index < UINT_MAX);
1931
1932 std::vector<unsigned> old_ids;
1933 for (unsigned i=0; i<this->mElements[*it]->GetNumEdges(); ++i)
1934 {
1935 old_ids.push_back(this->mElements[*it]->GetEdge(i)->GetIndex());
1936 }
1937
1938 this->mElements[*it]->AddNode(pNodeB, nodeA_local_index);
1939 if (mTrackMeshOperations)
1940 {
1941 mOperationRecorder.RecordNewEdgeOperation(this->mElements[*it], nodeA_local_index);
1942 }
1943 }
1944 else
1945 {
1946 // If the element contains both nodes A and B (now C and D respectively)...
1947 unsigned nodeA_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeA->GetIndex());
1948 unsigned nodeB_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeB->GetIndex());
1949
1950 assert(nodeA_local_index < UINT_MAX);
1951 assert(nodeB_local_index < UINT_MAX);
1952
1953 /*
1954 * Locate local index of nodeA and nodeB and use the ordering to
1955 * identify the element, if nodeB_index > nodeA_index then element 4
1956 * and if nodeA_index > nodeB_index then element 2
1957 */
1958 unsigned nodeB_local_index_plus_one = (nodeB_local_index + 1) % (this->mElements[*it]->GetNumNodes());
1959 /*
1960 * T1 swap and subsequent A-B edge shrinkage is recorded as node merging
1961 */
1962 std::pair<unsigned, unsigned> deleted_node_indices;
1963 if (nodeA_local_index == nodeB_local_index_plus_one)
1964 {
1965 /*
1966 * In this case the local index of nodeA is the local index of
1967 * nodeB plus one so we are in element 2 so we remove nodeB
1968 */
1969 this->mElements[*it]->DeleteNode(nodeB_local_index);
1970 deleted_node_indices.first = nodeA_local_index;
1971 deleted_node_indices.second = nodeB_local_index;
1972 }
1973 else
1974 {
1975 assert(nodeB_local_index == (nodeA_local_index + 1) % (this->mElements[*it]->GetNumNodes())); // as A and B are next to each other
1976 /*
1977 * In this case the local index of nodeA is the local index of
1978 * nodeB minus one so we are in element 4 so we remove nodeA
1979 */
1980 this->mElements[*it]->DeleteNode(nodeA_local_index);
1981 deleted_node_indices.first = nodeB_local_index;
1982 deleted_node_indices.second = nodeA_local_index;
1983 }
1984 if (mTrackMeshOperations)
1985 {
1986 mOperationRecorder.RecordNodeMergeOperation(old_ids, this->mElements[*it],
1987 deleted_node_indices);
1988 }
1989 }
1990 rebuilt_elements.insert(this->mElements[*it]->GetIndex());
1991 }
1992 // Sort out boundary nodes
1993 if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
1994 {
1995 if (pNodeA->GetNumContainingElements() == 3)
1996 {
1997 pNodeA->SetAsBoundaryNode(false);
1998 }
1999 else
2000 {
2001 pNodeA->SetAsBoundaryNode(true);
2002 }
2003 if (pNodeB->GetNumContainingElements() == 3)
2004 {
2005 pNodeB->SetAsBoundaryNode(false);
2006 }
2007 else
2008 {
2009 pNodeB->SetAsBoundaryNode(true);
2010 }
2011 }
2012
2013 for (unsigned i : rebuilt_elements)
2014 {
2015 this->GetElement(i)->RebuildEdges();
2016 }
2017 }
2018 else
2019 {
2020 NEVER_REACHED;
2021 }
2022}
2023
2024template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2025void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformIntersectionSwap(
2026 [[maybe_unused]] Node<SPACE_DIM>* pNode, [[maybe_unused]] unsigned elementIndex)
2027{
2028 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
2029 {
2030 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elementIndex);
2031 unsigned num_nodes = p_element->GetNumNodes();
2032
2033 std::set<unsigned> elements_containing_intersecting_node;
2034
2035 for (unsigned node_local_index = 0; node_local_index < num_nodes; node_local_index++)
2036 {
2037 unsigned node_global_index = p_element->GetNodeGlobalIndex(node_local_index);
2038
2039 std::set<unsigned> node_elem_indices = this->GetNode(node_global_index)->rGetContainingElementIndices();
2040
2041 for (std::set<unsigned>::const_iterator elem_iter = node_elem_indices.begin();
2042 elem_iter != node_elem_indices.end();
2043 ++elem_iter)
2044 {
2045 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_neighbouring_element = this->GetElement(*elem_iter);
2046 unsigned num_nodes_in_neighbouring_element = p_neighbouring_element->GetNumNodes();
2047
2048 // Check if element contains the intersecting node
2049 for (unsigned node_index_2 = 0; node_index_2 < num_nodes_in_neighbouring_element; node_index_2++)
2050 {
2051 if (p_neighbouring_element->GetNodeGlobalIndex(node_index_2) == pNode->GetIndex())
2052 {
2053 elements_containing_intersecting_node.insert(p_neighbouring_element->GetIndex());
2054 }
2055 }
2056 }
2057 }
2058
2059 std::set<unsigned> all_elements_containing_intersecting_node = pNode->rGetContainingElementIndices();
2060
2061 assert(elements_containing_intersecting_node.size() >= 1);
2062 assert(all_elements_containing_intersecting_node.size() >= 1);
2063
2064 /*
2065 * Identify nodes and elements to perform switch on
2066 * Intersecting node is node A
2067 * Other node is node B
2068 *
2069 * Element 1 only contains node A
2070 * Element 2 has nodes B and A (in that order)
2071 * Element 3 only contains node B
2072 * Element 4 has nodes A and B (in that order)
2073 *
2074 * If node A is a boundary node, then elements 1, 2, or 4 can be missing.
2075 */
2076 unsigned node_A_index = pNode->GetIndex();
2077 unsigned node_B_index = UINT_MAX;
2078
2079 unsigned element_1_index = UINT_MAX;
2080 unsigned element_2_index = UINT_MAX;
2081 unsigned element_3_index = elementIndex;
2082 unsigned element_4_index = UINT_MAX;
2083
2084 // Get element 1
2085 std::set<unsigned> intersecting_element;
2086
2087 std::set_difference(all_elements_containing_intersecting_node.begin(), all_elements_containing_intersecting_node.end(),
2088 elements_containing_intersecting_node.begin(), elements_containing_intersecting_node.end(),
2089 std::inserter(intersecting_element, intersecting_element.begin()));
2090
2091 if (intersecting_element.size() == 1)
2092 {
2093 element_1_index = *(intersecting_element.begin());
2094 }
2095
2096 // Get element A
2097 std::set<unsigned>::iterator iter = elements_containing_intersecting_node.begin();
2098 unsigned element_a_index = *(iter);
2099 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_a = this->GetElement(element_a_index);
2100
2101 // Get element B
2102 unsigned element_b_index = UINT_MAX;
2103 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_b = nullptr;
2104 if (elements_containing_intersecting_node.size() == 2)
2105 {
2106 iter++;
2107 element_b_index = *(iter);
2108 p_element_b = this->GetElement(element_b_index);
2109 }
2110
2111 // T1 swaps can sometimes result in a concave triangular element (#3067)
2112 if ((p_element_a->GetNumNodes() == 3) || ((p_element_b != nullptr) && p_element_b->GetNumNodes() == 3))
2113 {
2114 EXCEPTION("A triangular element has become concave. "
2115 "You need to rerun the simulation with a smaller time step to prevent this.");
2116 }
2117
2118 // Get node B index
2119 unsigned node_A_local_index_in_a = p_element_a->GetNodeLocalIndex(node_A_index);
2120
2121 unsigned node_before_A_in_a = (node_A_local_index_in_a + p_element_a->GetNumNodes() - 1) % p_element_a->GetNumNodes();
2122 unsigned node_after_A_in_a = (node_A_local_index_in_a + 1) % p_element_a->GetNumNodes();
2123
2124 unsigned global_node_before_A_in_a = p_element_a->GetNodeGlobalIndex(node_before_A_in_a);
2125 unsigned global_node_after_A_in_a = p_element_a->GetNodeGlobalIndex(node_after_A_in_a);
2126
2127 for (unsigned node_index = 0; node_index < num_nodes; ++node_index)
2128 {
2129 if (p_element->GetNodeGlobalIndex(node_index) == global_node_before_A_in_a)
2130 {
2131 node_B_index = global_node_before_A_in_a;
2132 break;
2133 }
2134 else if (p_element->GetNodeGlobalIndex(node_index) == global_node_after_A_in_a)
2135 {
2136 node_B_index = global_node_after_A_in_a;
2137 break;
2138 }
2139 }
2140 /*
2141 * If node B cannot be found then there are no nodes before or after node A
2142 * in element a that are contained in the intersected element, and there is
2143 * no way to fix it unless you want to make two new elements.
2144 */
2145 if (node_B_index == UINT_MAX)
2146 {
2147 EXCEPTION("Intersection cannot be resolved without splitting the element into two new elements.");
2148 }
2149
2150 // Store location of intersection swap, which is at midpoint of nodes A and B
2151 c_vector<double, SPACE_DIM> nodeA_location = pNode->rGetLocation();
2152 c_vector<double, SPACE_DIM> nodeB_location = this->GetNode(node_B_index)->rGetLocation();
2153 c_vector<double, SPACE_DIM> vector_AB = this->GetVectorFromAtoB(nodeA_location, nodeB_location);
2154 mLocationsOfIntersectionSwaps.push_back(nodeA_location + 0.5 * vector_AB);
2155
2156 // Now identify elements 2 and 4
2157 unsigned node_B_local_index_in_a = p_element_a->GetNodeLocalIndex(node_B_index);
2158
2159 if ((node_B_local_index_in_a + 1) % p_element_a->GetNumNodes() == node_A_local_index_in_a)
2160 {
2161#ifndef NDEBUG
2162 if (element_b_index != UINT_MAX)
2163 {
2164 assert(p_element_b != nullptr);
2165 assert((p_element_b->GetNodeLocalIndex(node_A_index) + 1) % p_element_b->GetNumNodes()
2166 == p_element_b->GetNodeLocalIndex(node_B_index));
2167 }
2168#endif
2169
2170 // Element 2 is element a, element 4 is element b
2171 element_2_index = element_a_index;
2172 element_4_index = element_b_index;
2173 }
2174 else
2175 {
2176#ifndef NDEBUG
2177 if (element_b_index != UINT_MAX)
2178 {
2179 assert(p_element_b != nullptr);
2180 assert((p_element_b->GetNodeLocalIndex(node_B_index) + 1) % p_element_b->GetNumNodes()
2181 == p_element_b->GetNodeLocalIndex(node_A_index));
2182 }
2183#endif
2184
2185 // Element 2 is element b, element 4 is element a
2186 element_2_index = element_b_index;
2187 element_4_index = element_a_index;
2188 }
2189
2190 // Get local indices
2191 unsigned intersected_edge = this->GetLocalIndexForElementEdgeClosestToPoint(pNode->rGetLocation(), elementIndex);
2192
2193 unsigned node_A_local_index_in_1 = UINT_MAX;
2194 if (element_1_index != UINT_MAX)
2195 {
2196 node_A_local_index_in_1 = this->GetElement(element_1_index)->GetNodeLocalIndex(node_A_index);
2197 }
2198
2199 unsigned node_A_local_index_in_2 = UINT_MAX;
2200 unsigned node_B_local_index_in_2 = UINT_MAX;
2201
2202 if (element_2_index != UINT_MAX)
2203 {
2204 node_A_local_index_in_2 = this->GetElement(element_2_index)->GetNodeLocalIndex(node_A_index);
2205 node_B_local_index_in_2 = this->GetElement(element_2_index)->GetNodeLocalIndex(node_B_index);
2206 }
2207
2208 unsigned node_B_local_index_in_3 = this->GetElement(elementIndex)->GetNodeLocalIndex(node_B_index);
2209
2210 unsigned node_A_local_index_in_4 = UINT_MAX;
2211 unsigned node_B_local_index_in_4 = UINT_MAX;
2212
2213 if (element_4_index != UINT_MAX)
2214 {
2215 node_A_local_index_in_4 = this->GetElement(element_4_index)->GetNodeLocalIndex(node_A_index);
2216 node_B_local_index_in_4 = this->GetElement(element_4_index)->GetNodeLocalIndex(node_B_index);
2217 }
2218
2219 // Switch nodes
2220 if (intersected_edge == node_B_local_index_in_3)
2221 {
2222 /*
2223 * Add node B to element 1 after node A
2224 * Add node A to element 3 after node B
2225 *
2226 * Remove node B from element 2
2227 * Remove node A from element 4
2228 */
2229 if (element_1_index != UINT_MAX)
2230 {
2231 assert(node_A_local_index_in_1 != UINT_MAX);
2232 this->mElements[element_1_index]->AddNode(this->mNodes[node_B_index], node_A_local_index_in_1);
2233 if (mTrackMeshOperations)
2234 mOperationRecorder.RecordNewEdgeOperation(this->mElements[element_1_index], node_A_local_index_in_1);
2235 }
2236 this->mElements[element_3_index]->AddNode(this->mNodes[node_A_index], node_B_local_index_in_3);
2237 c_vector<double, SPACE_DIM> vector_B_to_node
2238 = this->GetVectorFromAtoB(this->GetElement(elementIndex)->GetNode(intersected_edge)->rGetLocation(), this->mNodes[node_B_index]->rGetLocation());
2239 c_vector<double, SPACE_DIM> vector_A_to_B
2240 = this->GetVectorFromAtoB(this->mNodes[node_B_index]->rGetLocation(), this->mNodes[node_A_index]->rGetLocation());
2241 // Insertion of node A splits the intersected edge
2242 if (mTrackMeshOperations)
2243 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), intersected_edge,
2244 norm_2(vector_A_to_B) / norm_2(vector_B_to_node));
2245
2246 if (element_2_index != UINT_MAX)
2247 {
2248 assert(node_B_local_index_in_2 != UINT_MAX);
2249 this->mElements[element_2_index]->DeleteNode(node_B_local_index_in_2);
2250
2251 if (mTrackMeshOperations)
2252 mOperationRecorder.RecordEdgeMergeOperation(this->mElements[element_2_index], node_B_local_index_in_2);
2253 }
2254 if (element_4_index != UINT_MAX)
2255 {
2256 assert(node_A_local_index_in_4 != UINT_MAX);
2257 this->mElements[element_4_index]->DeleteNode(node_A_local_index_in_4);
2258 if (mTrackMeshOperations)
2259 mOperationRecorder.RecordEdgeMergeOperation(this->mElements[element_4_index], node_A_local_index_in_4);
2260 }
2261 }
2262 else
2263 {
2264 assert((intersected_edge + 1) % num_nodes == node_B_local_index_in_3);
2265
2266 // Add node B to element 1 before node A and add node A to element 3 before node B
2267 if (element_1_index != UINT_MAX)
2268 {
2269 assert(node_A_local_index_in_1 != UINT_MAX);
2270 unsigned node_before_A_in_1 = (node_A_local_index_in_1 + this->GetElement(element_1_index)->GetNumNodes() - 1) % this->GetElement(element_1_index)->GetNumNodes();
2271 this->mElements[element_1_index]->AddNode(this->mNodes[node_B_index], node_before_A_in_1);
2272 // Insertion of node B creates a new edge
2273 if (mTrackMeshOperations)
2274 mOperationRecorder.RecordNewEdgeOperation(this->mElements[element_1_index], node_before_A_in_1 + 1);
2275 }
2276
2277 unsigned node_before_B_in_3 = (node_B_local_index_in_3 + this->GetElement(element_3_index)->GetNumNodes() - 1) % this->GetElement(element_3_index)->GetNumNodes();
2278 this->mElements[element_3_index]->AddNode(this->mNodes[node_A_index], node_before_B_in_3);
2279 c_vector<double, SPACE_DIM> vector_B_to_node
2280 = this->GetVectorFromAtoB(this->GetElement(elementIndex)->GetNode(intersected_edge)->rGetLocation(), this->mNodes[node_B_index]->rGetLocation());
2281 c_vector<double, SPACE_DIM> vector_A_to_B
2282 = this->GetVectorFromAtoB(this->mNodes[node_B_index]->rGetLocation(), this->mNodes[node_A_index]->rGetLocation());
2283 // Insertion of node A splits the intersected edge
2284 if (mTrackMeshOperations)
2285 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), intersected_edge,
2286 norm_2(vector_A_to_B) / norm_2(vector_B_to_node));
2287
2288 // Remove node A from element 2 and remove node B from element 4
2289 if (element_2_index != UINT_MAX)
2290 {
2291 assert(node_A_local_index_in_2 != UINT_MAX);
2292 this->mElements[element_2_index]->DeleteNode(node_A_local_index_in_2);
2293 if (mTrackMeshOperations)
2294 mOperationRecorder.RecordEdgeMergeOperation(this->mElements[element_2_index], node_A_local_index_in_2);
2295 }
2296 if (element_4_index != UINT_MAX)
2297 {
2298 assert(node_B_local_index_in_4 != UINT_MAX);
2299 this->mElements[element_4_index]->DeleteNode(node_B_local_index_in_4);
2300 if (mTrackMeshOperations)
2301 mOperationRecorder.RecordEdgeMergeOperation(this->mElements[element_4_index], node_B_local_index_in_4);
2302 }
2303 }
2304 if (element_1_index != UINT_MAX)
2305 {
2306 this->mElements[element_1_index]->RebuildEdges();
2307 }
2308 if (element_2_index != UINT_MAX)
2309 {
2310 this->mElements[element_2_index]->RebuildEdges();
2311 }
2312 if (element_3_index != UINT_MAX)
2313 {
2314 this->mElements[element_3_index]->RebuildEdges();
2315 }
2316 if (element_4_index != UINT_MAX)
2317 {
2318 this->mElements[element_4_index]->RebuildEdges();
2319 }
2320 // Set as boundary nodes if intersecting node is a boundary node
2321 if (all_elements_containing_intersecting_node.size() == 2)
2322 {
2323 // Case 1: element 1 is missing
2324 if (elements_containing_intersecting_node.size() == 2)
2325 {
2326 assert(this->mNodes[node_A_index]->IsBoundaryNode() == true);
2327 assert(this->mNodes[node_B_index]->IsBoundaryNode() == false);
2328 this->mNodes[node_B_index]->SetAsBoundaryNode(true);
2329 }
2330 else if (elements_containing_intersecting_node.size() == 1)
2331 {
2332 assert(this->mNodes[node_A_index]->IsBoundaryNode() == true);
2333 assert(this->mNodes[node_B_index]->IsBoundaryNode() == true);
2334 // Case 2: element 2 is missing
2335 if (element_2_index == UINT_MAX)
2336 {
2337 if (intersected_edge == node_B_local_index_in_3)
2338 {
2339 this->mNodes[node_B_index]->SetAsBoundaryNode(false);
2340 }
2341 else
2342 {
2343 this->mNodes[node_A_index]->SetAsBoundaryNode(false);
2344 }
2345 }
2346 // Case 3: element 4 is missing
2347 if (element_4_index == UINT_MAX)
2348 {
2349 if (intersected_edge == node_B_local_index_in_3)
2350 {
2351 this->mNodes[node_A_index]->SetAsBoundaryNode(false);
2352 }
2353 else
2354 {
2355 this->mNodes[node_B_index]->SetAsBoundaryNode(false);
2356 }
2357 }
2358 }
2359 else
2360 {
2361 NEVER_REACHED;
2362 }
2363 }
2364#ifndef NDEBUG
2365 // Case 4: elements 1 and 2 are missing
2366 // Case 5: elements 1 and 4 are missing
2367 else if (all_elements_containing_intersecting_node.size() == 1)
2368 {
2369 if (elements_containing_intersecting_node.size() == 1)
2370 {
2371 assert(this->mNodes[node_A_index]->IsBoundaryNode() == true);
2372 assert(this->mNodes[node_B_index]->IsBoundaryNode() == true);
2373 }
2374 else
2375 {
2376 NEVER_REACHED;
2377 }
2378 }
2379 // Node A is not a boundary node
2380 else
2381 {
2382 assert(this->mNodes[node_A_index]->IsBoundaryNode() == false);
2383 assert(this->mNodes[node_B_index]->IsBoundaryNode() == false);
2384 }
2385#endif
2386 }
2387 else
2388 {
2389 NEVER_REACHED;
2390 }
2391}
2392
2393template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2394void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformT2Swap([[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>& rElement)
2395{
2396 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
2397 {
2398 // The given element must be triangular for us to be able to perform a T2 swap on it
2399 assert(rElement.GetNumNodes() == 3);
2400 // Note that we define this vector before setting it, as otherwise the profiling build will break (see #2367)
2401 c_vector<double, SPACE_DIM> new_node_location;
2402 new_node_location = this->GetCentroidOfElement(rElement.GetIndex());
2403 mLastT2SwapLocation = new_node_location;
2404
2405 T2SwapInfo<SPACE_DIM> swap_info;
2406 swap_info.mCellId = rElement.GetIndex();
2407 swap_info.mLocation = new_node_location;
2408 mOperationRecorder.RecordT2Swap(swap_info);
2409
2410 // If this element has no neighbours, delete the element, its nodes and exit function
2411 if (this->GetNeighbouringElementIndices(rElement.GetIndex()).size() == 0)
2412 {
2413 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(0));
2414 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(1));
2415 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(2));
2416
2417 rElement.GetNode(0)->MarkAsDeleted();
2418 rElement.GetNode(1)->MarkAsDeleted();
2419 rElement.GetNode(2)->MarkAsDeleted();
2420
2421 mDeletedElementIndices.push_back(rElement.GetIndex());
2422 rElement.MarkAsDeleted();
2423
2424 return;
2425 }
2426
2427 // Create a new node at the element's centroid; this will be a boundary node if any existing nodes were on the boundary
2428 bool is_node_on_boundary = false;
2429 for (unsigned i = 0; i < 3; i++)
2430 {
2431 if (rElement.GetNode(i)->IsBoundaryNode())
2432 {
2433 is_node_on_boundary = true;
2434 break;
2435 }
2436 }
2437 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(GetNumNodes(), new_node_location, is_node_on_boundary));
2438 Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
2439
2440 std::set<unsigned> neigh_indices;
2441
2442 // Loop over each of the three nodes contained in rElement
2443 for (unsigned i = 0; i < 3; i++)
2444 {
2445 // For each node, find the set of other elements containing it
2446 Node<SPACE_DIM>* p_node = rElement.GetNode(i);
2447
2448 std::set<unsigned> containing_elements = p_node->rGetContainingElementIndices();
2449 containing_elements.erase(rElement.GetIndex());
2450 // For each of these elements...
2451 for (std::set<unsigned>::iterator elem_iter = containing_elements.begin(); elem_iter != containing_elements.end(); ++elem_iter)
2452 {
2453 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_this_elem = this->GetElement(*elem_iter);
2454
2455 neigh_indices.insert(p_this_elem->GetIndex());
2456
2457 // ...throw an exception if the element is triangular...
2458 if (p_this_elem->GetNumNodes() < 4)
2459 {
2460 EXCEPTION("One of the neighbours of a small triangular element is also a triangle - dealing with this has not been implemented yet");
2461 }
2462
2463 // ...otherwise, replace p_node with p_new_node unless this has already happened (in which case, delete p_node from the element)
2464 if (p_this_elem->GetNodeLocalIndex(new_node_global_index) == UINT_MAX)
2465 {
2466 p_this_elem->ReplaceNode(p_node, p_new_node);
2467 }
2468 else
2469 {
2470 std::vector<unsigned> old_ids;
2471 std::pair<unsigned, unsigned> node_pair;
2472 if (mTrackMeshOperations)
2473 {
2474 for (unsigned k = 0; k < p_this_elem->GetNumEdges(); ++k)
2475 {
2476 old_ids.push_back(p_this_elem->GetEdge(k)->GetIndex());
2477 }
2478 node_pair.first = p_this_elem->GetNodeLocalIndex(new_node_global_index);
2479 node_pair.second = p_this_elem->GetNodeLocalIndex(p_node->GetIndex());
2480 }
2481 p_this_elem->DeleteNode(p_this_elem->GetNodeLocalIndex(p_node->GetIndex()));
2482 if (mTrackMeshOperations)
2483 {
2484 mOperationRecorder.RecordNodeMergeOperation(old_ids, p_this_elem, node_pair, true);
2485 }
2486 }
2487 }
2488 }
2489
2490 // We also have to mark pElement, pElement->GetNode(0), pElement->GetNode(1), and pElement->GetNode(2) as deleted
2491 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(0));
2492 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(1));
2493 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(2));
2494
2495 rElement.GetNode(0)->MarkAsDeleted();
2496 rElement.GetNode(1)->MarkAsDeleted();
2497 rElement.GetNode(2)->MarkAsDeleted();
2498
2499 mDeletedElementIndices.push_back(rElement.GetIndex());
2500 rElement.MarkAsDeleted();
2501
2502 for (unsigned i : neigh_indices)
2503 {
2504 this->GetElement(i)->RebuildEdges();
2505 }
2506 }
2507 else
2508 {
2509 NEVER_REACHED;
2510 }
2511}
2512
2513template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2514void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformT3Swap(
2515 [[maybe_unused]] Node<SPACE_DIM>* pNode, [[maybe_unused]] unsigned elementIndex)
2516{
2517 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
2518 {
2519
2520 assert(pNode->IsBoundaryNode());
2521
2522 // Get element
2523 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elementIndex);
2524 unsigned num_nodes = p_element->GetNumNodes();
2525
2526 // Store the index of the elements containing the intersecting node
2527 std::set<unsigned> elements_containing_intersecting_node = pNode->rGetContainingElementIndices();
2528
2529 // Get the local index of the node in the intersected element after which the new node is to be added
2530 unsigned node_A_local_index = this->GetLocalIndexForElementEdgeClosestToPoint(pNode->rGetLocation(), elementIndex);
2531
2532 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2533 c_vector<double, SPACE_DIM> node_location;
2534 node_location = pNode->rGetModifiableLocation();
2535
2536 // Get the nodes at either end of the edge to be divided
2537 unsigned vertexA_index = p_element->GetNodeGlobalIndex(node_A_local_index);
2538 unsigned vertexB_index = p_element->GetNodeGlobalIndex((node_A_local_index + 1) % num_nodes);
2539
2540 // Check these nodes are also boundary nodes if this fails then the elements have become concave and you need a smaller timestep
2541 if (!this->mNodes[vertexA_index]->IsBoundaryNode() || !this->mNodes[vertexB_index]->IsBoundaryNode())
2542 {
2543 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.");
2544 }
2545
2546 // Get the nodes at either end of the edge to be divided and calculate intersection
2547 c_vector<double, SPACE_DIM> vertexA = p_element->GetNodeLocation(node_A_local_index);
2548 c_vector<double, SPACE_DIM> vertexB = p_element->GetNodeLocation((node_A_local_index + 1) % num_nodes);
2549 c_vector<double, SPACE_DIM> vector_a_to_point = this->GetVectorFromAtoB(vertexA, node_location);
2550
2551 c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
2552
2553 c_vector<double, SPACE_DIM> edge_ab_unit_vector = vector_a_to_b / norm_2(vector_a_to_b);
2554 c_vector<double, SPACE_DIM> intersection = vertexA + edge_ab_unit_vector * inner_prod(vector_a_to_point, edge_ab_unit_vector);
2555
2556 // Store the location of the T3 swap, the location of the intersection with the edge
2558 // is called (see #2401) - we should correct this in these cases!
2559 T3SwapInfo<SPACE_DIM> swap_info;
2560 swap_info.mLocation = intersection;
2561 mOperationRecorder.RecordT3Swap(swap_info);
2562
2563 // For rebuilding edges
2564 std::set<unsigned> rebuilt_elements;
2565 if (pNode->GetNumContainingElements() == 1)
2566 {
2567 // Get the index of the element containing the intersecting node
2568 unsigned intersecting_element_index = *elements_containing_intersecting_node.begin();
2569
2570 // Get element
2571 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_intersecting_element = this->GetElement(intersecting_element_index);
2572
2573 //Edge map before swap
2574 const unsigned num_edges = p_intersecting_element->GetNumEdges();
2575 std::vector<unsigned> old_ids(num_edges);
2576 for (unsigned i=0; i<num_edges; ++i)
2577 {
2578 old_ids[i] = p_intersecting_element->GetEdge(i)->GetIndex();
2579 }
2580
2581 unsigned local_index = p_intersecting_element->GetNodeLocalIndex(pNode->GetIndex());
2582 unsigned next_node = p_intersecting_element->GetNodeGlobalIndex((local_index + 1) % (p_intersecting_element->GetNumNodes()));
2583 unsigned previous_node = p_intersecting_element->GetNodeGlobalIndex((local_index + p_intersecting_element->GetNumNodes() - 1) % (p_intersecting_element->GetNumNodes()));
2584
2585 // Check to see if the nodes adjacent to the intersecting node are contained in the intersected element between vertices A and B
2586 if (next_node == vertexA_index || previous_node == vertexA_index || next_node == vertexB_index || previous_node == vertexB_index)
2587 {
2588 unsigned common_vertex_index;
2589
2590 if (next_node == vertexA_index || previous_node == vertexA_index)
2591 {
2592 common_vertex_index = vertexA_index;
2593 }
2594 else
2595 {
2596 common_vertex_index = vertexB_index;
2597 }
2598
2599 assert(this->mNodes[common_vertex_index]->GetNumContainingElements() > 1);
2600
2601 std::set<unsigned> elements_containing_common_vertex = this->mNodes[common_vertex_index]->rGetContainingElementIndices();
2602 std::set<unsigned>::const_iterator it = elements_containing_common_vertex.begin();
2603 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*it);
2604 it++;
2605 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*it);
2606
2607 // Find the number and indices of common vertices between element_1 and element_2
2608 unsigned num_common_vertices = 0;
2609 std::vector<unsigned> common_vertex_indices;
2610 for (unsigned i = 0; i < p_element_common_1->GetNumNodes(); i++)
2611 {
2612 for (unsigned j = 0; j < p_element_common_2->GetNumNodes(); j++)
2613 {
2614 if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
2615 {
2616 num_common_vertices++;
2617 common_vertex_indices.push_back(p_element_common_1->GetNodeGlobalIndex(i));
2618 }
2619 }
2620 }
2621
2622 if (num_common_vertices == 1 || this->mNodes[common_vertex_index]->GetNumContainingElements() > 2)
2623 {
2624 /*
2625 * This is the situation here.
2626 *
2627 * From To
2628 * _ _
2629 * | <--- |
2630 * | /\ |\
2631 * | / \ | \
2632 * _|/____\ _|__\
2633 *
2634 * The edge goes from vertexA--vertexB to vertexA--pNode--vertexB
2635 */
2636
2637 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2638 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2639
2640 // Move original node
2641 pNode->rGetModifiableLocation() = intersection;
2642
2643 // Record the edge split data
2644 const double a_to_b_length = norm_2(vector_a_to_b);
2645 c_vector<double, SPACE_DIM> vector_a_to_node = this->GetVectorFromAtoB(vertexA, intersection);
2646 const double a_to_node_length = norm_2(vector_a_to_node);
2647
2648 // Add the moved nodes to the element (this also updates the node)
2649 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2650
2651 if (mTrackMeshOperations)
2652 mOperationRecorder.RecordEdgeSplitOperation(p_element, node_A_local_index, a_to_node_length / a_to_b_length);
2653 rebuilt_elements.insert(p_element->GetIndex());
2654 // Check the nodes are updated correctly
2655 assert(pNode->GetNumContainingElements() == 2);
2656 }
2657 else if (num_common_vertices == 2)
2658 {
2659 // The two elements must have an edge in common. Find whether the common edge is the same as the
2660 // edge that is merged onto.
2661
2662 if ((common_vertex_indices[0] == vertexA_index && common_vertex_indices[1] == vertexB_index) || (common_vertex_indices[1] == vertexA_index && common_vertex_indices[0] == vertexB_index))
2663 {
2664 /*
2665 * Due to a previous T3 swap the situation looks like this.
2666 *
2667 * pNode
2668 * \ |\ /
2669 * \ | \ /
2670 * \_______|__\/
2671 * /A | B
2672 * / \
2673 *
2674 * A T3 Swap would merge pNode onto an edge of its own element.
2675 * We prevent this by just removing pNode. By doing this we also avoid the
2676 * intersecting element to be concave.
2677 */
2678
2679 // Obtain necessary data for event recording
2680 const unsigned downstream_index = (local_index + p_intersecting_element->GetNumNodes() - 1) % (p_intersecting_element->GetNumNodes());
2681
2682 // Delete pNode in the intersecting element
2683 p_intersecting_element->DeleteNode(local_index);
2684
2685 // We record node deletion as node merging.
2686 if (mTrackMeshOperations)
2687 {
2688 mOperationRecorder.RecordNodeMergeOperation(old_ids, p_intersecting_element, std::pair<unsigned, unsigned>(downstream_index, local_index));
2689 }
2690 rebuilt_elements.insert(p_intersecting_element->GetIndex());
2691
2692 // Mark all three nodes as deleted
2693 pNode->MarkAsDeleted();
2694 mDeletedNodeIndices.push_back(pNode->GetIndex());
2695 }
2696 else
2697 {
2698 /*
2699 * This is the situation here.
2700 *
2701 * C is common_vertex D is the other one.
2702 *
2703 * From To
2704 * _ D _
2705 * | <--- |
2706 * | /\ |\
2707 * C|/ \ | \
2708 * _|____\ _|__\
2709 *
2710 * The edge goes from vertexC--vertexB to vertexC--pNode--vertexD
2711 * then vertex B is removed as it is no longer needed.
2712 */
2713
2714 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2715 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2716
2717 // Move original node
2718 pNode->rGetModifiableLocation() = intersection;
2719
2720 // Replace common_vertex with the moved node (this also updates the nodes)
2721 this->GetElement(elementIndex)->ReplaceNode(this->mNodes[common_vertex_index], pNode);
2722
2723 // Remove common_vertex
2724 unsigned common_vertex_local_index = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(common_vertex_index);
2725 this->GetElement(intersecting_element_index)->DeleteNode(common_vertex_local_index);
2726 assert(this->mNodes[common_vertex_index]->GetNumContainingElements() == 0);
2727
2728 // Record edge merging in the intersecting element
2729 if (mTrackMeshOperations)
2730 mOperationRecorder.RecordEdgeMergeOperation(p_intersecting_element, common_vertex_local_index);
2731 rebuilt_elements.insert(p_intersecting_element->GetIndex());
2732
2733 this->mNodes[common_vertex_index]->MarkAsDeleted();
2734 mDeletedNodeIndices.push_back(common_vertex_index);
2735
2736 // Check the nodes are updated correctly
2737 assert(pNode->GetNumContainingElements() == 2);
2738 }
2739 }
2740 else if (num_common_vertices == 4)
2741 {
2742 /*
2743 * The two elements share edges CA and BD due to previous swaps but not the edge AB
2744 *
2745 * From To
2746 * D___ D___
2747 * | |
2748 * B|\ |
2749 * | \ |
2750 * | / |
2751 * A|/ |
2752 * C_|__ C_|__
2753 *
2754 * We just remove the intersecting node as well as vertices A and B.
2755 */
2756
2757 // Delete node A and B in the intersected element
2758 this->GetElement(elementIndex)->DeleteNode(node_A_local_index);
2759 if (mTrackMeshOperations)
2760 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(elementIndex), node_A_local_index);
2761 unsigned node_B_local_index = this->GetElement(elementIndex)->GetNodeLocalIndex(vertexB_index);
2762 this->GetElement(elementIndex)->DeleteNode(node_B_local_index);
2763 if (mTrackMeshOperations)
2764 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(elementIndex), node_B_local_index);
2765
2766 // Delete nodes A and B in the intersecting element
2767 unsigned node_A_local_index_intersecting_element = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(vertexA_index);
2768 this->GetElement(intersecting_element_index)->DeleteNode(node_A_local_index_intersecting_element);
2769 if (mTrackMeshOperations)
2770 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(intersecting_element_index), node_A_local_index_intersecting_element);
2771 unsigned node_B_local_index_intersecting_element = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(vertexB_index);
2772 this->GetElement(intersecting_element_index)->DeleteNode(node_B_local_index_intersecting_element);
2773 if (mTrackMeshOperations)
2774 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(intersecting_element_index), node_B_local_index_intersecting_element);
2775
2776 // Delete pNode in the intersecting element
2777 unsigned p_node_local_index = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(pNode->GetIndex());
2778 this->GetElement(intersecting_element_index)->DeleteNode(p_node_local_index);
2779 if (mTrackMeshOperations)
2780 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(intersecting_element_index), p_node_local_index);
2781
2782 rebuilt_elements.insert(elementIndex);
2783 rebuilt_elements.insert(intersecting_element_index);
2784 // Mark all three nodes as deleted
2785 pNode->MarkAsDeleted();
2786 mDeletedNodeIndices.push_back(pNode->GetIndex());
2787 this->mNodes[vertexA_index]->MarkAsDeleted();
2788 mDeletedNodeIndices.push_back(vertexA_index);
2789 this->mNodes[vertexB_index]->MarkAsDeleted();
2790 mDeletedNodeIndices.push_back(vertexB_index);
2791 }
2792 else
2793 {
2794 // This can't happen as nodes can't be on the internal edge of 2 elements.
2795 NEVER_REACHED;
2796 }
2797 }
2798 else
2799 {
2800 /*
2801 * From To
2802 * ____ _______
2803 * / \
2804 * /\ ^ / \
2805 * / \ |
2806 *
2807 * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
2808 */
2809
2810 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2811 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2812 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
2813
2814 // Move original node
2815 pNode->rGetModifiableLocation() = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
2816
2817 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2818 c_vector<double, SPACE_DIM> new_node_location;
2819 new_node_location = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
2820
2821 // Add new node which will always be a boundary node
2822 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
2823
2824 // Add the moved node (this also updates the node)
2825 // and record the split of the edge vertexA--vertexB due to insertion of pNode
2826 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2827 // Recording edge split
2828 c_vector<double, SPACE_DIM> vector_a_to_node = this->GetVectorFromAtoB(pNode->rGetLocation(), vertexA);
2829 if (mTrackMeshOperations)
2830 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index,
2831 norm_2(vector_a_to_node) / norm_2(vector_a_to_b));
2832
2833 // Add the new node to the element (this also updates the node)
2834 // and record the split of the edge vertexA--pNode by insertion of the new node
2835 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
2836
2837 // Recording edge split
2838 c_vector<double, SPACE_DIM> vector_a_to_new_node = this->GetVectorFromAtoB(new_node_location, vertexA);
2839 if (mTrackMeshOperations)
2840 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index,
2841 norm_2(vector_a_to_new_node) / norm_2(vector_a_to_node));
2842 // New node must be between vertexA and pNode
2843 assert(norm_2(vector_a_to_new_node) / norm_2(vector_a_to_node) < 1);
2844
2845 // Add the new node to the original element containing pNode (this also updates the node)
2846 const unsigned new_node_local_index = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(pNode->GetIndex());
2847 this->GetElement(intersecting_element_index)->AddNode(this->mNodes[new_node_global_index], new_node_local_index);
2848 if (mTrackMeshOperations)
2849 {
2850 mOperationRecorder.RecordNewEdgeOperation(this->GetElement(intersecting_element_index), new_node_local_index);
2851 }
2852
2853 // The nodes must have been updated correctly
2854 assert(pNode->GetNumContainingElements() == 2);
2855 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
2856
2857 rebuilt_elements.insert(elementIndex);
2858 rebuilt_elements.insert(intersecting_element_index);
2859 }
2860 }
2861 else if (pNode->GetNumContainingElements() == 2)
2862 {
2863 // Find the nodes contained in elements containing the intersecting node
2864 std::set<unsigned>::const_iterator it = elements_containing_intersecting_node.begin();
2865
2866 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_1 = this->GetElement(*it);
2867 unsigned num_nodes_elem_1 = p_element_1->GetNumNodes();
2868 it++;
2869
2870 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_2 = this->GetElement(*it);
2871 unsigned num_nodes_elem_2 = p_element_2->GetNumNodes();
2872
2873 unsigned node_global_index = pNode->GetIndex();
2874
2875 unsigned local_index_1 = p_element_1->GetNodeLocalIndex(node_global_index);
2876 unsigned next_node_1 = p_element_1->GetNodeGlobalIndex((local_index_1 + 1) % num_nodes_elem_1);
2877 unsigned previous_node_1 = p_element_1->GetNodeGlobalIndex((local_index_1 + num_nodes_elem_1 - 1) % num_nodes_elem_1);
2878
2879 unsigned local_index_2 = p_element_2->GetNodeLocalIndex(node_global_index);
2880 unsigned next_node_2 = p_element_2->GetNodeGlobalIndex((local_index_2 + 1) % num_nodes_elem_2);
2881 unsigned previous_node_2 = p_element_2->GetNodeGlobalIndex((local_index_2 + num_nodes_elem_2 - 1) % num_nodes_elem_2);
2882
2883 // Check to see if the nodes adjacent to the intersecting node are contained in the intersected element between vertices A and B
2884 if ((next_node_1 == vertexA_index || previous_node_1 == vertexA_index || next_node_2 == vertexA_index || previous_node_2 == vertexA_index) && (next_node_1 == vertexB_index || previous_node_1 == vertexB_index || next_node_2 == vertexB_index || previous_node_2 == vertexB_index))
2885 {
2886 /*
2887 * Here we have
2888 * __
2889 * /| /
2890 * __ / | --> ___/
2891 * \ | \
2892 * \|__ \
2893 *
2894 * Where the node on the left has overlapped the edge A B
2895 *
2896 * Move p_node to the intersection on A B and merge AB and p_node
2897 */
2898
2899 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2900 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2901 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
2902
2903 // Check they are all boundary nodes
2904 assert(pNode->IsBoundaryNode());
2905 assert(this->mNodes[vertexA_index]->IsBoundaryNode());
2906 assert(this->mNodes[vertexB_index]->IsBoundaryNode());
2907
2908 // Move p_node to the intersection with the edge AB
2909 pNode->rGetModifiableLocation() = intersection;
2910 pNode->SetAsBoundaryNode(false);
2911
2912 // Add pNode to the intersected element
2913 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2914
2915 // Record this as edge split
2916 const double a_to_b_length = norm_2(vector_a_to_b);
2917 c_vector<double, SPACE_DIM> vector_a_to_p = this->GetVectorFromAtoB(intersection, vertexA);
2918 const double split_ratio = norm_2(vector_a_to_p) / a_to_b_length;
2919 if (mTrackMeshOperations)
2920 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, split_ratio);
2921 rebuilt_elements.insert(elementIndex);
2922
2923 // Remove vertex A from elements and record this as edge merge operation
2924 std::set<unsigned> elements_containing_vertex_A = this->mNodes[vertexA_index]->rGetContainingElementIndices();
2925 for (std::set<unsigned>::const_iterator iter = elements_containing_vertex_A.begin();
2926 iter != elements_containing_vertex_A.end();
2927 iter++)
2928 {
2929 const unsigned this_vertexA_local_index = this->GetElement(*iter)->GetNodeLocalIndex(vertexA_index);
2930 this->GetElement(*iter)->DeleteNode(this_vertexA_local_index);
2931 if (mTrackMeshOperations)
2932 {
2933 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(*iter), this_vertexA_local_index);
2934 }
2935 rebuilt_elements.insert(this->GetElement(*iter)->GetIndex());
2936 }
2937
2938 // Remove vertex A from the mesh
2939 assert(this->mNodes[vertexA_index]->GetNumContainingElements() == 0);
2940 this->mNodes[vertexA_index]->MarkAsDeleted();
2941 mDeletedNodeIndices.push_back(vertexA_index);
2942
2943 // Remove vertex B from elements and record this as edge merge operation
2944 std::set<unsigned> elements_containing_vertex_B = this->mNodes[vertexB_index]->rGetContainingElementIndices();
2945 for (std::set<unsigned>::const_iterator iter = elements_containing_vertex_B.begin();
2946 iter != elements_containing_vertex_B.end();
2947 iter++)
2948 {
2949 const unsigned this_vertexB_local_index = this->GetElement(*iter)->GetNodeLocalIndex(vertexB_index);
2950 this->GetElement(*iter)->DeleteNode(this_vertexB_local_index);
2951 if (mTrackMeshOperations)
2952 {
2953 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(*iter), this_vertexB_local_index);
2954 }
2955 rebuilt_elements.insert(this->GetElement(*iter)->GetIndex());
2956 }
2957
2958 // Remove vertex B from the mesh
2959 assert(this->mNodes[vertexB_index]->GetNumContainingElements() == 0);
2960 this->mNodes[vertexB_index]->MarkAsDeleted();
2961 mDeletedNodeIndices.push_back(vertexB_index);
2962 }
2963 else
2964 {
2965 if (next_node_1 == vertexA_index || previous_node_1 == vertexA_index || next_node_2 == vertexA_index || previous_node_2 == vertexA_index)
2966 {
2967 // Get elements containing vertexA_index (the common vertex)
2968
2969 assert(this->mNodes[vertexA_index]->GetNumContainingElements() > 1);
2970
2971 std::set<unsigned> elements_containing_vertex_A = this->mNodes[vertexA_index]->rGetContainingElementIndices();
2972 std::set<unsigned>::const_iterator iter = elements_containing_vertex_A.begin();
2973 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*iter);
2974 iter++;
2975 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*iter);
2976
2977 // Calculate the number of common vertices between element_1 and element_2
2978 unsigned num_common_vertices = 0;
2979 for (unsigned i = 0; i < p_element_common_1->GetNumNodes(); i++)
2980 {
2981 for (unsigned j = 0; j < p_element_common_2->GetNumNodes(); j++)
2982 {
2983 if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
2984 {
2985 num_common_vertices++;
2986 }
2987 }
2988 }
2989
2990 if (num_common_vertices == 1 || this->mNodes[vertexA_index]->GetNumContainingElements() > 2)
2991 {
2992 /*
2993 * From To
2994 * _ B _ B
2995 * | <--- |
2996 * | /|\ |\
2997 * | / | \ | \
2998 * | / | \ |\ \
2999 * _|/___|___\ _|_\_\
3000 * A A
3001 *
3002 * The edge goes from vertexA--vertexB to vertexA--pNode--new_node--vertexB
3003 */
3004
3005 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3006 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3007 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3008
3009 // Move original node and change to non-boundary node
3010 pNode->rGetModifiableLocation() = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3011 pNode->SetAsBoundaryNode(false);
3012
3013 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
3014 c_vector<double, SPACE_DIM> new_node_location;
3015 new_node_location = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3016
3017 // Add new node, which will always be a boundary node
3018 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
3019
3020 // Add the moved nodes to the element (this also updates the node)
3021 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
3022
3023 const double a_to_b_length = norm_2(vector_a_to_b);
3024 c_vector<double, SPACE_DIM> a_to_new_node_vector = this->GetVectorFromAtoB(vertexA, new_node_location);
3025 const double a_to_new_length = norm_2(a_to_new_node_vector);
3026 const double ratio_new_node = a_to_new_length / a_to_b_length;
3027 if (mTrackMeshOperations)
3028 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_new_node);
3029 rebuilt_elements.insert(elementIndex);
3030
3031 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3032
3033 c_vector<double, SPACE_DIM> a_to_p_vector = this->GetVectorFromAtoB(vertexA, pNode->rGetLocation());
3034 const double ratio_p_node = norm_2(a_to_p_vector) / a_to_new_length;
3035 assert(ratio_p_node <= 1);
3036 if (mTrackMeshOperations)
3037 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_p_node);
3038
3039 // Add the new nodes to the original elements containing pNode (this also updates the node)
3040 if (next_node_1 == previous_node_2)
3041 {
3042 const unsigned insertion_local_index = (local_index_1 + p_element_1->GetNumNodes() - 1) % (p_element_1->GetNumNodes());
3043 p_element_1->AddNode(this->mNodes[new_node_global_index], insertion_local_index);
3044 if (mTrackMeshOperations)
3045 {
3046 mOperationRecorder.RecordNewEdgeOperation(p_element_1, insertion_local_index);
3047 }
3048 rebuilt_elements.insert(p_element_1->GetIndex());
3049 }
3050 else
3051 {
3052 assert(next_node_2 == previous_node_1);
3053 const unsigned insertion_local_index = (local_index_2 + p_element_2->GetNumNodes() - 1) % (p_element_2->GetNumNodes());
3054 p_element_2->AddNode(this->mNodes[new_node_global_index], insertion_local_index);
3055 if (mTrackMeshOperations)
3056 {
3057 mOperationRecorder.RecordNewEdgeOperation(p_element_2, insertion_local_index);
3058 }
3059 rebuilt_elements.insert(p_element_2->GetIndex());
3060 }
3061
3062 // Check the nodes are updated correctly
3063 assert(pNode->GetNumContainingElements() == 3);
3064 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
3065 }
3066 else if (num_common_vertices == 2)
3067 {
3068 /*
3069 * From To
3070 * _ B _ B
3071 * |<--- |
3072 * | /|\ |\
3073 * |/ | \ | \
3074 * | | \ |\ \
3075 * _|__|___\ _|_\_\
3076 * A A
3077 *
3078 * The edge goes from vertexA--vertexB to vertexA--pNode--new_node--vertexB
3079 * then vertexA is removed
3080 */
3081
3082 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3083 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3084 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3085
3086 // Move original node and change to non-boundary node
3087 pNode->rGetModifiableLocation() = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3088 pNode->SetAsBoundaryNode(false);
3089
3090 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
3091 c_vector<double, SPACE_DIM> new_node_location;
3092 new_node_location = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3093
3094 // Add new node, which will always be a boundary node
3095 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
3096
3097 // Add the moved nodes to the element (this also updates the node)
3098 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
3099
3100 const double a_to_b_length = norm_2(vector_a_to_b);
3101 c_vector<double, SPACE_DIM> a_to_new_node_vector = this->GetVectorFromAtoB(vertexA, new_node_location);
3102 const double a_to_new_length = norm_2(a_to_new_node_vector);
3103 const double ratio_new_node = a_to_new_length / a_to_b_length;
3104 if (mTrackMeshOperations)
3105 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_new_node);
3106
3107 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3108
3109 c_vector<double, SPACE_DIM> a_to_p_vector = this->GetVectorFromAtoB(vertexA, pNode->rGetLocation());
3110 const double ratio_p_node = norm_2(a_to_p_vector) / a_to_new_length;
3111 assert(ratio_p_node <= 1);
3112 if (mTrackMeshOperations)
3113 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_p_node);
3114
3115 rebuilt_elements.insert(elementIndex);
3116 // Add the new nodes to the original elements containing pNode (this also updates the node)
3117 if (next_node_1 == previous_node_2)
3118 {
3119 const unsigned insertion_local_index = (local_index_1 + p_element_1->GetNumNodes() - 1) % (p_element_1->GetNumNodes());
3120 p_element_1->AddNode(this->mNodes[new_node_global_index], insertion_local_index);
3121 if (mTrackMeshOperations)
3122 {
3123 mOperationRecorder.RecordNewEdgeOperation(p_element_1, insertion_local_index);
3124 }
3125 rebuilt_elements.insert(p_element_1->GetIndex());
3126 }
3127 else
3128 {
3129 assert(next_node_2 == previous_node_1);
3130 const unsigned insertion_local_index = (local_index_2 + p_element_2->GetNumNodes() - 1) % (p_element_2->GetNumNodes());
3131 p_element_2->AddNode(this->mNodes[new_node_global_index], insertion_local_index);
3132 if (mTrackMeshOperations)
3133 {
3134 mOperationRecorder.RecordNewEdgeOperation(p_element_2, insertion_local_index);
3135 }
3136 rebuilt_elements.insert(p_element_2->GetIndex());
3137 }
3138
3139 // Remove vertex A from the mesh
3140 const unsigned local_index_1 = p_element_common_1->GetNodeLocalIndex(vertexA_index);
3141 const unsigned local_index_2 = p_element_common_2->GetNodeLocalIndex(vertexA_index);
3142 p_element_common_1->DeleteNode(local_index_1);
3143 if (mTrackMeshOperations)
3144 {
3145 mOperationRecorder.RecordEdgeMergeOperation(p_element_common_1, local_index_1);
3146 }
3147 p_element_common_2->DeleteNode(local_index_2);
3148 if (mTrackMeshOperations)
3149 {
3150 mOperationRecorder.RecordEdgeMergeOperation(p_element_common_2, local_index_2);
3151 }
3152 assert(this->mNodes[vertexA_index]->GetNumContainingElements() == 0);
3153 rebuilt_elements.insert(p_element_common_1->GetIndex());
3154 rebuilt_elements.insert(p_element_common_2->GetIndex());
3155
3156 this->mNodes[vertexA_index]->MarkAsDeleted();
3157 mDeletedNodeIndices.push_back(vertexA_index);
3158
3159 // Check the nodes are updated correctly
3160 assert(pNode->GetNumContainingElements() == 3);
3161 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
3162 }
3163 else
3164 {
3165 // This can't happen as nodes can't be on the internal edge of two elements
3166 NEVER_REACHED;
3167 }
3168 }
3169 else if (next_node_1 == vertexB_index || previous_node_1 == vertexB_index || next_node_2 == vertexB_index || previous_node_2 == vertexB_index)
3170 {
3171 // Get elements containing vertexB_index (the common vertex)
3172
3173 assert(this->mNodes[vertexB_index]->GetNumContainingElements() > 1);
3174
3175 std::set<unsigned> elements_containing_vertex_B = this->mNodes[vertexB_index]->rGetContainingElementIndices();
3176 std::set<unsigned>::const_iterator iter = elements_containing_vertex_B.begin();
3177 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*iter);
3178 iter++;
3179 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*iter);
3180
3181 // Calculate the number of common vertices between element_1 and element_2
3182 unsigned num_common_vertices = 0;
3183 for (unsigned i = 0; i < p_element_common_1->GetNumNodes(); i++)
3184 {
3185 for (unsigned j = 0; j < p_element_common_2->GetNumNodes(); j++)
3186 {
3187 if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
3188 {
3189 num_common_vertices++;
3190 }
3191 }
3192 }
3193
3194 if (num_common_vertices == 1 || this->mNodes[vertexB_index]->GetNumContainingElements() > 2)
3195 {
3196 /*
3197 * From To
3198 * _B_________ _B____
3199 * |\ | / | / /
3200 * | \ | / |/ /
3201 * | \ | / | /
3202 * | \|/ |/
3203 * _| <--- _|
3204 * A
3205 *
3206 * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
3207 */
3208
3209 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3210 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3211 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3212
3213 // Move original node and change to non-boundary node
3214 pNode->rGetModifiableLocation() = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3215 pNode->SetAsBoundaryNode(false);
3216
3217 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
3218 c_vector<double, SPACE_DIM> new_node_location;
3219 new_node_location = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3220
3221 // Add new node which will always be a boundary node
3222 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
3223
3224 // Add the moved nodes to the element (this also updates the node)
3225 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3226 const double a_to_b_length = norm_2(vector_a_to_b);
3227 c_vector<double, SPACE_DIM> a_to_p_vector = this->GetVectorFromAtoB(vertexA, pNode->rGetLocation());
3228 const double ratio_p_node = norm_2(a_to_p_vector) / a_to_b_length;
3229 assert(ratio_p_node <= 1);
3230 if (mTrackMeshOperations)
3231 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_p_node);
3232
3233 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
3234 c_vector<double, SPACE_DIM> a_to_new_node_vector = this->GetVectorFromAtoB(vertexA, new_node_location);
3235 const double a_to_new_length = norm_2(a_to_new_node_vector);
3236 const double ratio_new_node = a_to_new_length / norm_2(a_to_p_vector);
3237 assert(ratio_new_node <= 1);
3238 if (mTrackMeshOperations)
3239 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_new_node);
3240 rebuilt_elements.insert(elementIndex);
3241 // Add the new nodes to the original elements containing pNode (this also updates the node)
3242 if (next_node_1 == previous_node_2)
3243 {
3244 p_element_2->AddNode(this->mNodes[new_node_global_index], local_index_2);
3245 if (mTrackMeshOperations)
3246 mOperationRecorder.RecordNewEdgeOperation(p_element_2, local_index_2);
3247 rebuilt_elements.insert(p_element_2->GetIndex());
3248 }
3249 else
3250 {
3251 assert(next_node_2 == previous_node_1);
3252 p_element_1->AddNode(this->mNodes[new_node_global_index], local_index_1);
3253 if (mTrackMeshOperations)
3254 mOperationRecorder.RecordNewEdgeOperation(p_element_1, local_index_1);
3255 rebuilt_elements.insert(p_element_1->GetIndex());
3256 }
3257
3258 // Check the nodes are updated correctly
3259 assert(pNode->GetNumContainingElements() == 3);
3260 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
3261 }
3262 else if (num_common_vertices == 2)
3263 {
3264 /*
3265 * From To
3266 * _B_______ _B____
3267 * | | / | / /
3268 * | | / |/ /
3269 * |\ | / | /
3270 * | \|/ |/
3271 * _| <--- _|
3272 * A
3273 *
3274 * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
3275 * then vertexB is removed
3276 */
3277
3278 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3279 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3280 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3281
3282 // Move original node and change to non-boundary node
3283 pNode->rGetModifiableLocation() = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3284 pNode->SetAsBoundaryNode(false);
3285
3286 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
3287 c_vector<double, SPACE_DIM> new_node_location;
3288 new_node_location = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3289
3290 // Add new node which will always be a boundary node
3291 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
3292
3293 // Add the moved nodes to the element (this also updates the node)
3294 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3295 const double a_to_b_length = norm_2(vector_a_to_b);
3296 c_vector<double, SPACE_DIM> a_to_p_vector = this->GetVectorFromAtoB(vertexA, pNode->rGetLocation());
3297 const double ratio_p_node = norm_2(a_to_p_vector) / a_to_b_length;
3298 assert(ratio_p_node <= 1);
3299 if (mTrackMeshOperations)
3300 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_p_node);
3301
3302 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
3303 c_vector<double, SPACE_DIM> a_to_new_node_vector = this->GetVectorFromAtoB(vertexA, new_node_location);
3304 const double a_to_new_length = norm_2(a_to_new_node_vector);
3305 const double ratio_new_node = a_to_new_length / norm_2(a_to_p_vector);
3306 assert(ratio_new_node <= 1);
3307 if (mTrackMeshOperations)
3308 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_new_node);
3309 rebuilt_elements.insert(elementIndex);
3310 // Add the new nodes to the original elements containing pNode (this also updates the node)
3311 if (next_node_1 == previous_node_2)
3312 {
3313 p_element_2->AddNode(this->mNodes[new_node_global_index], local_index_2);
3314 if (mTrackMeshOperations)
3315 mOperationRecorder.RecordNewEdgeOperation(p_element_2, local_index_2);
3316 rebuilt_elements.insert(p_element_2->GetIndex());
3317 }
3318 else
3319 {
3320 assert(next_node_2 == previous_node_1);
3321 p_element_1->AddNode(this->mNodes[new_node_global_index], local_index_1);
3322 if (mTrackMeshOperations)
3323 mOperationRecorder.RecordNewEdgeOperation(p_element_1, local_index_1);
3324 rebuilt_elements.insert(p_element_1->GetIndex());
3325 }
3326
3327 // Remove vertex B from the mesh
3328 const unsigned local_index_1 = p_element_common_1->GetNodeLocalIndex(vertexB_index);
3329 const unsigned local_index_2 = p_element_common_2->GetNodeLocalIndex(vertexB_index);
3330 p_element_common_1->DeleteNode(local_index_1);
3331 if (mTrackMeshOperations)
3332 {
3333 mOperationRecorder.RecordEdgeMergeOperation(p_element_common_1, local_index_1);
3334 }
3335 p_element_common_2->DeleteNode(local_index_2);
3336 if (mTrackMeshOperations)
3337 {
3338 mOperationRecorder.RecordEdgeMergeOperation(p_element_common_2, local_index_2);
3339 }
3340 rebuilt_elements.insert(p_element_common_1->GetIndex());
3341 rebuilt_elements.insert(p_element_common_2->GetIndex());
3342
3343 assert(this->mNodes[vertexB_index]->GetNumContainingElements() == 0);
3344
3345 this->mNodes[vertexB_index]->MarkAsDeleted();
3346 mDeletedNodeIndices.push_back(vertexB_index);
3347
3348 // Check the nodes are updated correctly
3349 assert(pNode->GetNumContainingElements() == 3);
3350 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
3351 }
3352 else
3353 {
3354 // This can't happen as nodes can't be on the internal edge of two elements
3355 NEVER_REACHED;
3356 }
3357 }
3358 else
3359 {
3360 /*
3361 * From To
3362 * _____ _______
3363 * / | \
3364 * /|\ ^ / | \
3365 * / | \ |
3366 *
3367 * The edge goes from vertexA--vertexB to vertexA--new_node_1--pNode--new_node_2--vertexB
3368 */
3369
3370 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3371 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3372 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3373
3374 // Move original node and change to non-boundary node
3375 pNode->rGetModifiableLocation() = intersection;
3376 pNode->SetAsBoundaryNode(false);
3377
3378 c_vector<double, SPACE_DIM> new_node_1_location;
3379 new_node_1_location = intersection - mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3380 c_vector<double, SPACE_DIM> new_node_2_location;
3381 new_node_2_location = intersection + mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3382
3383 // Add new nodes which will always be boundary nodes
3384 unsigned new_node_1_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_1_location[0], new_node_1_location[1]));
3385 unsigned new_node_2_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_2_location[0], new_node_2_location[1]));
3386
3387 // Add the moved and new nodes to the element (this also updates the node)
3388 const double a_to_b_length = norm_2(vector_a_to_b);
3389 const c_vector<double, SPACE_DIM> a_to_node2_vector
3390 = this->GetVectorFromAtoB(new_node_2_location, vertexA);
3391 const double a_to_node2_length = norm_2(a_to_node2_vector);
3392 const c_vector<double, SPACE_DIM> a_to_nodeP_vector
3393 = this->GetVectorFromAtoB(intersection, vertexA);
3394 const double a_to_nodeP_length = norm_2(a_to_nodeP_vector);
3395 const c_vector<double, SPACE_DIM> a_to_node1_vector
3396 = this->GetVectorFromAtoB(new_node_1_location, vertexA);
3397 const double a_to_node1_length = norm_2(a_to_node1_vector);
3398
3399 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_2_global_index], node_A_local_index);
3400 if (mTrackMeshOperations)
3401 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, a_to_node2_length / a_to_b_length);
3402 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3403 if (mTrackMeshOperations)
3404 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, a_to_nodeP_length / a_to_node2_length);
3405 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_1_global_index], node_A_local_index);
3406 if (mTrackMeshOperations)
3407 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, a_to_node1_length / a_to_nodeP_length);
3408 rebuilt_elements.insert(elementIndex);
3409 // Add the new nodes to the original elements containing pNode (this also updates the node)
3410 if (next_node_1 == previous_node_2)
3411 {
3412 const unsigned inserted_local_index_1 = (local_index_1 + p_element_1->GetNumNodes() - 1) % (p_element_1->GetNumNodes());
3413 p_element_1->AddNode(this->mNodes[new_node_2_global_index], inserted_local_index_1);
3414 if (mTrackMeshOperations)
3415 {
3416 mOperationRecorder.RecordNewEdgeOperation(p_element_1, inserted_local_index_1);
3417 }
3418 p_element_2->AddNode(this->mNodes[new_node_1_global_index], local_index_2);
3419 if (mTrackMeshOperations)
3420 {
3421 mOperationRecorder.RecordNewEdgeOperation(p_element_2, local_index_2);
3422 }
3423 }
3424 else
3425 {
3426 assert(next_node_2 == previous_node_1);
3427 p_element_1->AddNode(this->mNodes[new_node_1_global_index], local_index_1);
3428 if (mTrackMeshOperations)
3429 {
3430 mOperationRecorder.RecordNewEdgeOperation(p_element_1, local_index_1);
3431 }
3432 const unsigned inserted_local_index_2 = (local_index_2 + p_element_2->GetNumNodes() - 1) % (p_element_2->GetNumNodes());
3433 p_element_2->AddNode(this->mNodes[new_node_2_global_index], inserted_local_index_2);
3434 if (mTrackMeshOperations)
3435 {
3436 mOperationRecorder.RecordNewEdgeOperation(p_element_2, local_index_2);
3437 }
3438 }
3439 rebuilt_elements.insert(p_element_1->GetIndex());
3440 rebuilt_elements.insert(p_element_2->GetIndex());
3441
3442 // Check the nodes are updated correctly
3443 assert(pNode->GetNumContainingElements() == 3);
3444 assert(this->mNodes[new_node_1_global_index]->GetNumContainingElements() == 2);
3445 assert(this->mNodes[new_node_2_global_index]->GetNumContainingElements() == 2);
3446 }
3447 }
3448 }
3449 else
3450 {
3451 EXCEPTION("Trying to merge a node, contained in more than 2 elements, into another element, this is not possible with the vertex mesh.");
3452 }
3453 for (unsigned i : rebuilt_elements)
3454 {
3455 this->GetElement(i)->RebuildEdges();
3456 }
3457 }
3458 else
3459 {
3460 NEVER_REACHED;
3461 }
3462}
3463
3464template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3465void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformVoidRemoval(
3466 [[maybe_unused]] Node<SPACE_DIM>* pNodeA,
3467 [[maybe_unused]] Node<SPACE_DIM>* pNodeB,
3468 [[maybe_unused]] Node<SPACE_DIM>* pNodeC)
3469{
3470 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
3471 {
3472 // Calculate void centroid
3473 c_vector<double, SPACE_DIM> nodes_midpoint = pNodeA->rGetLocation()
3474 + this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation()) / 3.0
3475 + this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeC->rGetLocation()) / 3.0;
3476 /*
3477 * In two steps, merge nodes A, B and C into a single node. This is implemented in such a way that
3478 * the ordering of their indices does not matter.
3479 */
3480 PerformNodeMerge(pNodeA, pNodeB);
3481 Node<SPACE_DIM>* p_merged_node = pNodeB;
3482
3483 if (pNodeB->IsDeleted())
3484 {
3485 p_merged_node = pNodeA;
3486 }
3487
3488 PerformNodeMerge(pNodeC, p_merged_node);
3489 if (p_merged_node->IsDeleted())
3490 {
3491 p_merged_node = pNodeC;
3492 }
3493
3494 p_merged_node->rGetModifiableLocation() = nodes_midpoint;
3495
3496 // Tag remaining node as non-boundary
3497 p_merged_node->SetAsBoundaryNode(false);
3498
3499 // Remove the deleted nodes and re-index
3500 RemoveDeletedNodes();
3501 }
3502 else
3503 {
3504 NEVER_REACHED;
3505 }
3506}
3507
3508template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3509void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::HandleHighOrderJunctions(Node<SPACE_DIM>* pNodeA, Node<SPACE_DIM>* pNodeB)
3510{
3511 unsigned node_a_rank = pNodeA->rGetContainingElementIndices().size();
3512 unsigned node_b_rank = pNodeB->rGetContainingElementIndices().size();
3513
3514 if ((node_a_rank > 3) && (node_b_rank > 3))
3515 {
3516 // The code can't handle this case
3517 EXCEPTION("Both nodes involved in a swap event are contained in more than three elements");
3518 }
3519 else // the rosette degree should increase in this case
3520 {
3521 assert(node_a_rank > 3 || node_b_rank > 3);
3522 this->PerformRosetteRankIncrease(pNodeA, pNodeB);
3523 this->RemoveDeletedNodes();
3524 }
3525}
3526
3527template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3528void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformRosetteRankIncrease(
3529 [[maybe_unused]] Node<SPACE_DIM>* pNodeA, [[maybe_unused]] Node<SPACE_DIM>* pNodeB)
3530{
3531 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
3532 {
3533 /*
3534 * One of the nodes will have 3 containing element indices, the other
3535 * will have at least four. We first identify which node is which.
3536 */
3537
3538 unsigned node_a_index = pNodeA->GetIndex();
3539 unsigned node_b_index = pNodeB->GetIndex();
3540
3541 unsigned node_a_rank = pNodeA->rGetContainingElementIndices().size();
3542 unsigned node_b_rank = pNodeB->rGetContainingElementIndices().size();
3543
3544 unsigned lo_rank_index = (node_a_rank < node_b_rank) ? node_a_index : node_b_index;
3545 unsigned hi_rank_index = (node_a_rank < node_b_rank) ? node_b_index : node_a_index;
3546
3547 // Get pointers to the nodes, sorted by index
3548 Node<SPACE_DIM>* p_lo_rank_node = this->GetNode(lo_rank_index);
3549 Node<SPACE_DIM>* p_hi_rank_node = this->GetNode(hi_rank_index);
3550
3551 // Find the sets of elements containing each of the nodes, sorted by index
3552 std::set<unsigned> lo_rank_elem_indices = p_lo_rank_node->rGetContainingElementIndices();
3553 std::set<unsigned> hi_rank_elem_indices = p_hi_rank_node->rGetContainingElementIndices();
3554
3581 for (std::set<unsigned>::const_iterator it = lo_rank_elem_indices.begin();
3582 it != lo_rank_elem_indices.end();
3583 ++it)
3584 {
3585 // Find the local index of lo_rank_node in this element
3586 unsigned lo_rank_local_index = this->mElements[*it]->GetNodeLocalIndex(lo_rank_index);
3587 assert(lo_rank_local_index < UINT_MAX); // double check this element contains lo_rank_node
3588
3589 /*
3590 * If this element already contains the hi_rank_node, we are in the situation of elements
3591 * C and D above, so we just remove lo_rank_node.
3592 *
3593 * Otherwise, we are in element E, so we must replace lo_rank_node with high-rank node,
3594 * and remove it from mNodes.
3595 *
3596 * We can check whether hi_rank_node is in this element using the set::count() method.
3597 */
3598
3599 if (hi_rank_elem_indices.count(*it) > 0)
3600 {
3601 // For node merge recording
3602 std::vector<unsigned> old_ids(this->mElements[*it]->GetNumEdges(),0);
3603
3604 for (unsigned i = 0; i < old_ids.size(); ++i)
3605 {
3606 old_ids[i] = this->mElements[*it]->GetEdge(i)->GetIndex();
3607 }
3608 // Delete lo_rank_node from current element
3609 this->mElements[*it]->DeleteNode(lo_rank_local_index);
3610
3611 // Record edge shrinkage
3612 const unsigned hi_rank_local_index = this->mElements[*it]->GetNodeLocalIndex(hi_rank_index);
3613 if (mTrackMeshOperations)
3614 {
3615 mOperationRecorder.RecordNodeMergeOperation(old_ids,
3616 this->mElements[*it],
3617 std::pair<unsigned, unsigned>(hi_rank_local_index, lo_rank_local_index));
3618 }
3619 }
3620 else
3621 {
3622 // Update lo_rank_node with all information (including index and location) of hi_rank_node
3623 this->mElements[*it]->UpdateNode(lo_rank_local_index, p_hi_rank_node);
3624 }
3625 this->mElements[*it]->RebuildEdges();
3626 }
3627
3628 // Tidy up the mesh by ensuring the global instance of lo_rank_node is deleted
3629 assert(!(this->mNodes[lo_rank_index]->IsDeleted()));
3630 this->mNodes[lo_rank_index]->MarkAsDeleted();
3631 this->mDeletedNodeIndices.push_back(lo_rank_index);
3632
3633 for (unsigned i : lo_rank_elem_indices)
3634 {
3635 this->GetElement(i)->RebuildEdges();
3636 }
3637 for (unsigned i : hi_rank_elem_indices)
3638 {
3639 this->GetElement(i)->RebuildEdges();
3640 }
3641 }
3642 else
3643 {
3644 NEVER_REACHED;
3645 }
3646}
3647
3648template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3649void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformProtorosetteResolution(
3650 [[maybe_unused]] Node<SPACE_DIM>* pProtorosetteNode)
3651{
3652 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
3653 {
3654 // Double check we are dealing with a protorosette
3655 assert(pProtorosetteNode->rGetContainingElementIndices().size() == 4);
3656
3657 // Get random number (0, 1, 2 or 3), as the resolution axis is assumed to be random
3658 unsigned random_elem_increment = RandomNumberGenerator::Instance()->randMod(4);
3659
3660 // Find global indices of elements around the protorosette node
3661 std::set<unsigned> protorosette_node_containing_elem_indices = pProtorosetteNode->rGetContainingElementIndices();
3662
3663 // Select random element by advancing iterator a random number times
3664 std::set<unsigned>::const_iterator elem_index_iter(protorosette_node_containing_elem_indices.begin());
3665 advance(elem_index_iter, random_elem_increment);
3666
3694 /*
3695 * We need to find the global indices of elements B, C and D. We do this with set intersections.
3696 */
3697
3698 unsigned elem_a_idx = *elem_index_iter;
3699 unsigned elem_b_idx = UINT_MAX;
3700 unsigned elem_c_idx = UINT_MAX;
3701 unsigned elem_d_idx = UINT_MAX;
3702
3703 // Get pointer to element we've chosen at random (element A)
3704 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_a = this->GetElement(elem_a_idx);
3705
3706 // Get all necessary info about element A and the protorosette node
3707 unsigned num_nodes_elem_a = p_elem_a->GetNumNodes();
3708 unsigned protorosette_node_global_idx = pProtorosetteNode->GetIndex();
3709 unsigned protorosette_node_local_idx = p_elem_a->GetNodeLocalIndex(protorosette_node_global_idx);
3710
3711 // Find global indices of previous (cw) and next (ccw) nodes, locally, from the protorosette node, in element A
3712 unsigned prev_node_global_idx = p_elem_a->GetNodeGlobalIndex((protorosette_node_local_idx + num_nodes_elem_a - 1) % num_nodes_elem_a);
3713 unsigned next_node_global_idx = p_elem_a->GetNodeGlobalIndex((protorosette_node_local_idx + 1) % num_nodes_elem_a);
3714
3715 // Get the set of elements the previous and next nodes are contained in
3716 Node<SPACE_DIM>* p_prev_node = this->GetNode(prev_node_global_idx);
3717 Node<SPACE_DIM>* p_next_node = this->GetNode(next_node_global_idx);
3718 std::set<unsigned> prev_node_elem_indices = p_prev_node->rGetContainingElementIndices();
3719 std::set<unsigned> next_node_elem_indices = p_next_node->rGetContainingElementIndices();
3720
3721 // Perform set intersections with the set of element indices which the protorosette node is contained in
3722 std::set<unsigned> intersection_with_prev;
3723 std::set<unsigned> intersection_with_next;
3724
3725 // This intersection should contain just global indices for elements A and B
3726 std::set_intersection(protorosette_node_containing_elem_indices.begin(),
3727 protorosette_node_containing_elem_indices.end(),
3728 prev_node_elem_indices.begin(),
3729 prev_node_elem_indices.end(),
3730 std::inserter(intersection_with_prev, intersection_with_prev.begin()));
3731
3732 // This intersection should contain just global indices for elements A and D
3733 std::set_intersection(protorosette_node_containing_elem_indices.begin(),
3734 protorosette_node_containing_elem_indices.end(),
3735 next_node_elem_indices.begin(),
3736 next_node_elem_indices.end(),
3737 std::inserter(intersection_with_next, intersection_with_next.begin()));
3738
3739 assert(intersection_with_prev.size() == 2);
3740 assert(intersection_with_next.size() == 2);
3741
3742 // Get global index of element B
3743 if (*intersection_with_prev.begin() != elem_a_idx)
3744 {
3745 elem_b_idx = *(intersection_with_prev.begin());
3746 }
3747 else
3748 {
3749 elem_b_idx = *(++(intersection_with_prev.begin()));
3750 }
3751 assert(elem_b_idx < UINT_MAX);
3752
3753 // Get global index of element D
3754 if (*intersection_with_next.begin() != elem_a_idx)
3755 {
3756 elem_d_idx = *(intersection_with_next.begin());
3757 }
3758 else
3759 {
3760 elem_d_idx = *(++(intersection_with_next.begin()));
3761 }
3762 assert(elem_d_idx < UINT_MAX);
3763
3764 // By elimination, the remaining unassigned index in the original set must be global index of element C
3765 for (elem_index_iter = protorosette_node_containing_elem_indices.begin();
3766 elem_index_iter != protorosette_node_containing_elem_indices.end();
3767 ++elem_index_iter)
3768 {
3769 if ((*elem_index_iter != elem_a_idx) && (*elem_index_iter != elem_b_idx) && (*elem_index_iter != elem_d_idx))
3770 {
3771 elem_c_idx = *elem_index_iter;
3772 }
3773 }
3774 assert(elem_c_idx < UINT_MAX);
3775
3790 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_b = this->GetElement(elem_b_idx);
3791 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_c = this->GetElement(elem_c_idx);
3792 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_d = this->GetElement(elem_d_idx);
3793
3794 double swap_distance = (this->mCellRearrangementRatio) * (this->mCellRearrangementThreshold);
3795
3796 // Get normalized vectors to centre of elements A and B from protorosette node
3797 c_vector<double, SPACE_DIM> node_to_elem_a_centre = this->GetCentroidOfElement(elem_a_idx) - pProtorosetteNode->rGetLocation();
3798 node_to_elem_a_centre /= norm_2(node_to_elem_a_centre);
3799
3800 c_vector<double, SPACE_DIM> node_to_elem_c_centre = this->GetCentroidOfElement(elem_c_idx) - pProtorosetteNode->rGetLocation();
3801 node_to_elem_c_centre /= norm_2(node_to_elem_c_centre);
3802
3803 // Calculate new node locations
3804 c_vector<double, SPACE_DIM> new_location_of_protorosette_node = pProtorosetteNode->rGetLocation() + (0.5 * swap_distance) * node_to_elem_a_centre;
3805 c_vector<double, SPACE_DIM> location_of_new_node = pProtorosetteNode->rGetLocation() + (0.5 * swap_distance) * node_to_elem_c_centre;
3806
3807 // Move protorosette node to new location
3808 pProtorosetteNode->rGetModifiableLocation() = new_location_of_protorosette_node;
3809
3810 // Create new node in correct location
3811 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(this->GetNumNodes(), location_of_new_node, false));
3812 Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
3813
3824 unsigned local_idx_elem_b = p_elem_b->GetNodeLocalIndex(protorosette_node_global_idx);
3825 local_idx_elem_b = (local_idx_elem_b + p_elem_b->GetNumNodes() - 1) % p_elem_b->GetNumNodes();
3826 unsigned local_idx_elem_c = p_elem_c->GetNodeLocalIndex(protorosette_node_global_idx);
3827 unsigned local_idx_elem_d = p_elem_d->GetNodeLocalIndex(protorosette_node_global_idx);
3828
3829 p_elem_b->AddNode(p_new_node, local_idx_elem_b);
3830 if (mTrackMeshOperations)
3831 mOperationRecorder.RecordNewEdgeOperation(p_elem_b, local_idx_elem_b);
3832 p_elem_c->AddNode(p_new_node, local_idx_elem_c);
3833 p_elem_d->AddNode(p_new_node, local_idx_elem_d);
3834 if (mTrackMeshOperations)
3835 mOperationRecorder.RecordNewEdgeOperation(p_elem_d, local_idx_elem_d);
3836
3837 // All that is left is to remove the original protorosette node from element C
3838 p_elem_c->DeleteNode(p_elem_c->GetNodeLocalIndex(protorosette_node_global_idx));
3839 }
3840 else
3841 {
3842 NEVER_REACHED;
3843 }
3844}
3845
3846template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3847void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformRosetteRankDecrease(
3848 [[maybe_unused]] Node<SPACE_DIM>* pRosetteNode)
3849{
3850 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
3851 {
3852 unsigned rosette_rank = pRosetteNode->rGetContainingElementIndices().size();
3853
3854 // Double check we're dealing with a rosette
3855 assert(rosette_rank > 4);
3856
3857 // Get random number in [0, 1, ..., n) where n is rank of rosette, as the resolution axis is assumed to be random
3858 unsigned random_elem_increment = RandomNumberGenerator::Instance()->randMod(rosette_rank);
3859
3860 // Find global indices of elements around the protorosette node
3861 std::set<unsigned> rosette_node_containing_elem_indices = pRosetteNode->rGetContainingElementIndices();
3862
3863 // Select random element by advancing iterator a random number times
3864 std::set<unsigned>::const_iterator elem_index_iter(rosette_node_containing_elem_indices.begin());
3865 advance(elem_index_iter, random_elem_increment);
3866
3905 /*
3906 * We need to find the global indices of elements N and P. We do this with set intersections.
3907 */
3908
3909 // Get the vertex element S (which we randomly selected)
3910 unsigned elem_s_idx = *elem_index_iter;
3911 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_s = this->GetElement(elem_s_idx);
3912
3913 unsigned elem_n_idx = UINT_MAX;
3914 unsigned elem_p_idx = UINT_MAX;
3915
3916 // Get all necessary info about element S and the rosette node
3917 unsigned num_nodes_elem_s = p_elem_s->GetNumNodes();
3918 unsigned rosette_node_global_idx = pRosetteNode->GetIndex();
3919 unsigned rosette_node_local_idx = p_elem_s->GetNodeLocalIndex(rosette_node_global_idx);
3920
3921 // Find global indices of previous (cw) and next (ccw) nodes, locally, from the rosette node, in element S
3922 unsigned prev_node_global_idx = p_elem_s->GetNodeGlobalIndex((rosette_node_local_idx + num_nodes_elem_s - 1) % num_nodes_elem_s);
3923 unsigned next_node_global_idx = p_elem_s->GetNodeGlobalIndex((rosette_node_local_idx + 1) % num_nodes_elem_s);
3924
3925 // Get the set of elements that the previous and next nodes are contained in
3926 Node<SPACE_DIM>* p_prev_node = this->GetNode(prev_node_global_idx);
3927 Node<SPACE_DIM>* p_next_node = this->GetNode(next_node_global_idx);
3928 std::set<unsigned> prev_node_elem_indices = p_prev_node->rGetContainingElementIndices();
3929 std::set<unsigned> next_node_elem_indices = p_next_node->rGetContainingElementIndices();
3930
3931 // Perform set intersections with the set of element indices that the rosette node is contained in
3932 std::set<unsigned> intersection_with_prev;
3933 std::set<unsigned> intersection_with_next;
3934
3935 // This intersection should contain just global indices for elements S and N
3936 std::set_intersection(rosette_node_containing_elem_indices.begin(),
3937 rosette_node_containing_elem_indices.end(),
3938 prev_node_elem_indices.begin(),
3939 prev_node_elem_indices.end(),
3940 std::inserter(intersection_with_prev, intersection_with_prev.begin()));
3941
3942 // This intersection should contain just global indices for elements S and P
3943 std::set_intersection(rosette_node_containing_elem_indices.begin(),
3944 rosette_node_containing_elem_indices.end(),
3945 next_node_elem_indices.begin(),
3946 next_node_elem_indices.end(),
3947 std::inserter(intersection_with_next, intersection_with_next.begin()));
3948
3949 assert(intersection_with_prev.size() == 2);
3950 assert(intersection_with_next.size() == 2);
3951
3952 // Get global index of element N
3953 if (*intersection_with_prev.begin() != elem_s_idx)
3954 {
3955 elem_n_idx = *intersection_with_prev.begin();
3956 }
3957 else
3958 {
3959 elem_n_idx = *(++(intersection_with_prev.begin()));
3960 }
3961 assert(elem_n_idx < UINT_MAX);
3962
3963 // Get global index of element P
3964 if (*intersection_with_next.begin() != elem_s_idx)
3965 {
3966 elem_p_idx = *intersection_with_next.begin();
3967 }
3968 else
3969 {
3970 elem_p_idx = *(++(intersection_with_next.begin()));
3971 }
3972 assert(elem_p_idx < UINT_MAX);
3973
3984 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_n = this->GetElement(elem_p_idx);
3985 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_p = this->GetElement(elem_n_idx);
3986
3987 double swap_distance = (this->mCellRearrangementRatio) * (this->mCellRearrangementThreshold);
3988
3989 // Calculate location of new node
3990 c_vector<double, 2> node_to_selected_elem = this->GetCentroidOfElement(elem_s_idx) - pRosetteNode->rGetLocation();
3991 node_to_selected_elem /= norm_2(node_to_selected_elem);
3992 c_vector<double, 2> new_node_location = pRosetteNode->rGetLocation() + (swap_distance * node_to_selected_elem);
3993
3994 // Create new node in correct location
3995 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(this->GetNumNodes(), new_node_location, false));
3996 Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
3997
4009 // Add new node, and remove rosette node, from element S
4010 unsigned node_local_idx_in_elem_s = p_elem_s->GetNodeLocalIndex(rosette_node_global_idx);
4011 p_elem_s->AddNode(p_new_node, node_local_idx_in_elem_s);
4012 p_elem_s->DeleteNode(node_local_idx_in_elem_s);
4013
4014 // Add new node to element N
4015 unsigned node_local_idx_in_elem_n = p_elem_n->GetNodeLocalIndex(rosette_node_global_idx);
4016 node_local_idx_in_elem_n = (node_local_idx_in_elem_n + p_elem_n->GetNumNodes() - 1) % p_elem_n->GetNumNodes();
4017 p_elem_n->AddNode(p_new_node, node_local_idx_in_elem_n);
4018 if (mTrackMeshOperations)
4019 {
4020 mOperationRecorder.RecordNewEdgeOperation(p_elem_n, node_local_idx_in_elem_n);
4021 }
4022 // Add new node to element P
4023 unsigned node_local_idx_in_elem_p = p_elem_p->GetNodeLocalIndex(rosette_node_global_idx);
4024 p_elem_p->AddNode(p_new_node, node_local_idx_in_elem_p);
4025 if (mTrackMeshOperations)
4026 {
4027 mOperationRecorder.RecordNewEdgeOperation(p_elem_p, node_local_idx_in_elem_p);
4028 }
4029 }
4030}
4031
4032template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
4033void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::CheckForRosettes()
4034{
4035 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
4036 {
4045 // Vectors to store the nodes that need resolution events
4046 std::vector<Node<SPACE_DIM>*> protorosette_nodes;
4047 std::vector<Node<SPACE_DIM>*> rosette_nodes;
4048
4049 // First loop in which we populate these vectors
4050 unsigned num_nodes = this->GetNumAllNodes();
4051 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
4052 {
4053 Node<SPACE_DIM>* current_node = this->GetNode(node_idx);
4054 unsigned node_rank = current_node->rGetContainingElementIndices().size();
4055
4056 if (node_rank < 4)
4057 {
4058 // Nothing to do if the node is not high-rank
4059 continue;
4060 }
4061 else if (node_rank == 4)
4062 {
4063 // For protorosette nodes, we check against a random number to decide if resolution is necessary
4064 if (mProtorosetteResolutionProbabilityPerTimestep >= RandomNumberGenerator::Instance()->ranf())
4065 {
4066 protorosette_nodes.push_back(current_node);
4067 }
4068 }
4069 else // if (node_rank > 4)
4070 {
4071 // For rosette nodes, we check against a random number to decide if resolution is necessary
4072 if (mRosetteResolutionProbabilityPerTimestep >= RandomNumberGenerator::Instance()->ranf())
4073 {
4074 rosette_nodes.push_back(current_node);
4075 }
4076 }
4077 }
4078
4086 // First, resolve any protorosettes
4087 for (unsigned node_idx = 0; node_idx < protorosette_nodes.size(); node_idx++)
4088 {
4089 Node<SPACE_DIM>* current_node = protorosette_nodes[node_idx];
4090
4091 // Verify that node has not been marked for deletion, and that it is still contained in four elements
4092 assert(!(current_node->IsDeleted()));
4093 assert(current_node->rGetContainingElementIndices().size() == 4);
4094
4095 // Perform protorosette resolution
4096 this->PerformProtorosetteResolution(current_node);
4097 }
4098
4099 // Finally, resolve any rosettes
4100 for (unsigned node_idx = 0; node_idx < rosette_nodes.size(); node_idx++)
4101 {
4102 Node<SPACE_DIM>* current_node = rosette_nodes[node_idx];
4103
4104 // Verify that node has not been marked for deletion, and that it is still contained in at least four elements
4105 assert(!(current_node->IsDeleted()));
4106 assert(current_node->rGetContainingElementIndices().size() > 4);
4107
4108 // Perform protorosette resolution
4109 this->PerformRosetteRankDecrease(current_node);
4110 }
4111 }
4112 else
4113 {
4114 NEVER_REACHED;
4115 }
4116}
4117
4118template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
4119c_vector<double, 2> MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::WidenEdgeOrCorrectIntersectionLocationIfNecessary(
4120 unsigned indexA, unsigned indexB, c_vector<double, 2> intersection)
4121{
4131 c_vector<double, SPACE_DIM> vertexA = this->GetNode(indexA)->rGetLocation();
4132 c_vector<double, SPACE_DIM> vertexB = this->GetNode(indexB)->rGetLocation();
4133 c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
4134
4135 if (norm_2(vector_a_to_b) < 4.0*mCellRearrangementRatio*mCellRearrangementThreshold)
4136 {
4137 WARNING("Trying to merge a node onto an edge which is too small.");
4138
4139 c_vector<double, SPACE_DIM> centre_a_and_b = vertexA + 0.5*vector_a_to_b;
4140
4141 vertexA = centre_a_and_b - 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*vector_a_to_b/norm_2(vector_a_to_b);
4142 ChastePoint<SPACE_DIM> vertex_A_point(vertexA);
4143 SetNode(indexA, vertex_A_point);
4144
4145 vertexB = centre_a_and_b + 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*vector_a_to_b/norm_2(vector_a_to_b);
4146 ChastePoint<SPACE_DIM> vertex_B_point(vertexB);
4147 SetNode(indexB, vertex_B_point);
4148
4149 intersection = centre_a_and_b;
4150 }
4151
4152 // Reset distances
4153 vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
4154 c_vector<double, 2> edge_ab_unit_vector = vector_a_to_b/norm_2(vector_a_to_b);
4155
4171 if (norm_2(intersection - vertexA) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
4172 {
4173 intersection = vertexA + 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
4174 }
4175 if (norm_2(intersection - vertexB) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
4176 {
4177 intersection = vertexB - 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
4178 }
4179 return intersection;
4180}
4181
4182template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
4183void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetMeshOperationTracking(const bool track)
4184{
4185 mTrackMeshOperations = track;
4186 if (track)
4187 {
4188 mOperationRecorder.SetEdgeHelper(&(this->mEdgeHelper));
4189 }
4190}
4191
4192template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
4198
4199// Explicit instantiation
4200template class MutableVertexMesh<1, 1>;
4201template class MutableVertexMesh<1, 2>;
4202template class MutableVertexMesh<1, 3>;
4203template class MutableVertexMesh<2, 2>;
4204template class MutableVertexMesh<2, 3>;
4205template class MutableVertexMesh<3, 3>;
4206
4207// Serialization for Boost >= 1.36
4208#include "SerializationExportWrapperForCpp.hpp"
4209EXPORT_TEMPLATE_CLASS_ALL_DIMS(MutableVertexMesh)
EXCEPTION
#define EXCEPTION(message)
Definition Exception.hpp:156
NEVER_REACHED
#define NEVER_REACHED
Definition Exception.hpp:227
EXPORT_TEMPLATE_CLASS_ALL_DIMS
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
Definition SerializationExportWrapper.hpp:358
AbstractElement::ReplaceNode
void ReplaceNode(Node< SPACE_DIM > *pOldNode, Node< SPACE_DIM > *pNewNode)
Definition AbstractElement.cpp:74
AbstractElement::GetNode
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
Definition AbstractElement.cpp:116
AbstractElement::GetNodeLocation
double GetNodeLocation(unsigned localIndex, unsigned dimension) const
Definition AbstractElement.cpp:89
AbstractElement::GetNumNodes
unsigned GetNumNodes() const
Definition AbstractElement.cpp:123
AbstractElement::GetNodeGlobalIndex
unsigned GetNodeGlobalIndex(unsigned localIndex) const
Definition AbstractElement.cpp:109
AbstractElement::GetIndex
unsigned GetIndex() const
Definition AbstractElement.cpp:141
AbstractMesh::mMeshChangesDuringSimulation
bool mMeshChangesDuringSimulation
Definition AbstractMesh.hpp:125
AbstractMesh::mNodes
std::vector< Node< SPACE_DIM > * > mNodes
Definition AbstractMesh.hpp:96
MutableElement::DeleteNode
void DeleteNode(const unsigned &rIndex)
Definition MutableElement.cpp:136
MutableElement::GetEdge
Edge< SPACE_DIM > * GetEdge(unsigned localIndex) const
Definition MutableElement.cpp:312
MutableElement::SetEdgeHelper
void SetEdgeHelper(EdgeHelper< SPACE_DIM > *pEdgeHelper)
Definition MutableElement.cpp:270
MutableElement::GetNodeLocalIndex
unsigned GetNodeLocalIndex(unsigned globalIndex) const
Definition MutableElement.cpp:225
MutableElement::AddNode
void AddNode(Node< SPACE_DIM > *pNode, const unsigned &rIndex)
Definition MutableElement.cpp:173
MutableElement::GetNumEdges
unsigned GetNumEdges() const
Definition MutableElement.cpp:319
MutableVertexMesh::SetCheckForT3Swaps
void SetCheckForT3Swaps(bool checkForT3Swaps)
Definition MutableVertexMesh.cpp:284
MutableVertexMesh::RemoveDeletedNodesAndElements
void RemoveDeletedNodesAndElements(VertexElementMap &rElementMap)
Definition MutableVertexMesh.cpp:931
MutableVertexMesh::PerformNodeMerge
void PerformNodeMerge(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
Definition MutableVertexMesh.cpp:1769
MutableVertexMesh::GetRosetteResolutionProbabilityPerTimestep
double GetRosetteResolutionProbabilityPerTimestep() const
Definition MutableVertexMesh.cpp:179
MutableVertexMesh::CheckForSwapsFromShortEdges
virtual bool CheckForSwapsFromShortEdges()
Definition MutableVertexMesh.cpp:1058
MutableVertexMesh::GetLastT2SwapLocation
c_vector< double, SPACE_DIM > GetLastT2SwapLocation()
Definition MutableVertexMesh.cpp:323
MutableVertexMesh::GetCellRearrangementThreshold
double GetCellRearrangementThreshold() const
Definition MutableVertexMesh.cpp:149
MutableVertexMesh::IdentifySwapType
virtual void IdentifySwapType(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
Definition MutableVertexMesh.cpp:1292
MutableVertexMesh::GetProtorosetteFormationProbability
double GetProtorosetteFormationProbability() const
Definition MutableVertexMesh.cpp:167
MutableVertexMesh::DeleteElementPriorToReMesh
void DeleteElementPriorToReMesh(unsigned index)
Definition MutableVertexMesh.cpp:819
MutableVertexMesh::PerformRosetteRankIncrease
void PerformRosetteRankIncrease(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
Definition MutableVertexMesh.cpp:3528
MutableVertexMesh::PerformT1Swap
void PerformT1Swap(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB, std::set< unsigned > &rElementsContainingNodes)
Definition MutableVertexMesh.cpp:1825
MutableVertexMesh::PerformVoidRemoval
void PerformVoidRemoval(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB, Node< SPACE_DIM > *pNodeC)
Definition MutableVertexMesh.cpp:3465
MutableVertexMesh::AddNode
unsigned AddNode(Node< SPACE_DIM > *pNewNode)
Definition MutableVertexMesh.cpp:365
MutableVertexMesh::SetCheckForInternalIntersections
void SetCheckForInternalIntersections(bool checkForInternalIntersections)
Definition MutableVertexMesh.cpp:278
MutableVertexMesh::GetCellRearrangementRatio
double GetCellRearrangementRatio() const
Definition MutableVertexMesh.cpp:161
MutableVertexMesh::GetDistanceForT3SwapChecking
double GetDistanceForT3SwapChecking() const
Definition MutableVertexMesh.cpp:191
MutableVertexMesh::GetLocationsOfT3Swaps
std::vector< c_vector< double, SPACE_DIM > > GetLocationsOfT3Swaps()
Definition MutableVertexMesh.cpp:329
MutableVertexMesh::SetRosetteResolutionProbabilityPerTimestep
void SetRosetteResolutionProbabilityPerTimestep(double rosetteResolutionProbabilityPerTimestep)
Definition MutableVertexMesh.cpp:261
MutableVertexMesh::CheckForT2Swaps
bool CheckForT2Swaps(VertexElementMap &rElementMap)
Definition MutableVertexMesh.cpp:1126
MutableVertexMesh::SetProtorosetteFormationProbability
void SetProtorosetteFormationProbability(double protorosetteFormationProbability)
Definition MutableVertexMesh.cpp:227
MutableVertexMesh::SetDistanceForT3SwapChecking
void SetDistanceForT3SwapChecking(double distanceForT3SwapChecking)
Definition MutableVertexMesh.cpp:185
MutableVertexMesh::DivideElementAlongGivenAxis
unsigned DivideElementAlongGivenAxis(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, c_vector< double, SPACE_DIM > axisOfDivision, bool placeOriginalElementBelow=false)
Definition MutableVertexMesh.cpp:410
MutableVertexMesh::GetNumNodes
unsigned GetNumNodes() const
Definition MutableVertexMesh.cpp:299
MutableVertexMesh::SetProtorosetteResolutionProbabilityPerTimestep
void SetProtorosetteResolutionProbabilityPerTimestep(double protorosetteResolutionProbabilityPerTimestep)
Definition MutableVertexMesh.cpp:244
MutableVertexMesh::DeleteNodePriorToReMesh
void DeleteNodePriorToReMesh(unsigned index)
Definition MutableVertexMesh.cpp:848
MutableVertexMesh::WidenEdgeOrCorrectIntersectionLocationIfNecessary
c_vector< double, 2 > WidenEdgeOrCorrectIntersectionLocationIfNecessary(unsigned indexA, unsigned indexB, c_vector< double, 2 > intersection)
Definition MutableVertexMesh.cpp:4119
MutableVertexMesh::GetLocationsOfT1Swaps
std::vector< c_vector< double, SPACE_DIM > > GetLocationsOfT1Swaps()
Definition MutableVertexMesh.cpp:311
MutableVertexMesh::PerformT3Swap
void PerformT3Swap(Node< SPACE_DIM > *pNode, unsigned elementIndex)
Definition MutableVertexMesh.cpp:2514
MutableVertexMesh::GetNumElements
unsigned GetNumElements() const
Definition MutableVertexMesh.cpp:305
MutableVertexMesh::GetProtorosetteResolutionProbabilityPerTimestep
double GetProtorosetteResolutionProbabilityPerTimestep() const
Definition MutableVertexMesh.cpp:173
MutableVertexMesh::DivideElementAlongShortAxis
unsigned DivideElementAlongShortAxis(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, bool placeOriginalElementBelow=false)
Definition MutableVertexMesh.cpp:668
MutableVertexMesh::PerformT2Swap
void PerformT2Swap(VertexElement< ELEMENT_DIM, SPACE_DIM > &rElement)
Definition MutableVertexMesh.cpp:2394
MutableVertexMesh::PerformProtorosetteResolution
void PerformProtorosetteResolution(Node< SPACE_DIM > *pProtorosetteNode)
Definition MutableVertexMesh.cpp:3649
MutableVertexMesh::AddElement
unsigned AddElement(VertexElement< ELEMENT_DIM, SPACE_DIM > *pNewElement)
Definition MutableVertexMesh.cpp:384
MutableVertexMesh::GetLocationsOfIntersectionSwaps
std::vector< c_vector< double, SPACE_DIM > > GetLocationsOfIntersectionSwaps()
Definition MutableVertexMesh.cpp:341
MutableVertexMesh::DivideEdge
void DivideEdge(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
Definition MutableVertexMesh.cpp:855
MutableVertexMesh::SetMeshOperationTracking
void SetMeshOperationTracking(const bool track)
Definition MutableVertexMesh.cpp:4183
MutableVertexMesh::HandleHighOrderJunctions
virtual void HandleHighOrderJunctions(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
Definition MutableVertexMesh.cpp:3509
MutableVertexMesh::SetCellRearrangementRatio
void SetCellRearrangementRatio(double cellRearrangementRatio)
Definition MutableVertexMesh.cpp:221
MutableVertexMesh::PerformRosetteRankDecrease
void PerformRosetteRankDecrease(Node< SPACE_DIM > *pRosetteNode)
Definition MutableVertexMesh.cpp:3847
MutableVertexMesh::SetNode
virtual void SetNode(unsigned nodeIndex, ChastePoint< SPACE_DIM > point)
Definition MutableVertexMesh.cpp:404
MutableVertexMesh::GetOperationRecorder
VertexMeshOperationRecorder< ELEMENT_DIM, SPACE_DIM > * GetOperationRecorder()
Definition MutableVertexMesh.cpp:4193
MutableVertexMesh::PerformIntersectionSwap
void PerformIntersectionSwap(Node< SPACE_DIM > *pNode, unsigned elementIndex)
Definition MutableVertexMesh.cpp:2025
MutableVertexMesh::GetCheckForT3Swaps
bool GetCheckForT3Swaps() const
Definition MutableVertexMesh.cpp:203
MutableVertexMesh::GetT2Threshold
double GetT2Threshold() const
Definition MutableVertexMesh.cpp:155
MutableVertexMesh::DivideElement
unsigned DivideElement(VertexElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned nodeAIndex, unsigned nodeBIndex, bool placeOriginalElementBelow=false)
Definition MutableVertexMesh.cpp:686
MutableVertexMesh::GetCheckForInternalIntersections
bool GetCheckForInternalIntersections() const
Definition MutableVertexMesh.cpp:197
MutableVertexMesh::SetT2Threshold
void SetT2Threshold(double t2Threshold)
Definition MutableVertexMesh.cpp:215
MutableVertexMesh::SetCellRearrangementThreshold
void SetCellRearrangementThreshold(double cellRearrangementThreshold)
Definition MutableVertexMesh.cpp:209
Node
Definition Node.hpp:59
Node::SetPoint
void SetPoint(ChastePoint< SPACE_DIM > point)
Definition Node.cpp:115
Node::rGetContainingElementIndices
std::set< unsigned > & rGetContainingElementIndices()
Definition Node.cpp:300
Node::rGetModifiableLocation
c_vector< double, SPACE_DIM > & rGetModifiableLocation()
Definition Node.cpp:151
Node::SetIndex
void SetIndex(unsigned index)
Definition Node.cpp:121
Node::AddElement
void AddElement(unsigned index)
Definition Node.cpp:268
Node::IsDeleted
bool IsDeleted() const
Definition Node.cpp:412
Node::rGetLocation
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition Node.cpp:139
Node::IsBoundaryNode
bool IsBoundaryNode() const
Definition Node.cpp:164
Node::GetIndex
unsigned GetIndex() const
Definition Node.cpp:158
Node::SetAsBoundaryNode
void SetAsBoundaryNode(bool value=true)
Definition Node.cpp:127
RandomNumberGenerator::Instance
static RandomNumberGenerator * Instance()
Definition RandomNumberGenerator.cpp:59
RandomNumberGenerator::randMod
unsigned randMod(unsigned base)
Definition RandomNumberGenerator.cpp:77
VertexElementMap::Resize
void Resize(unsigned size)
Definition VertexElementMap.cpp:50
VertexElementMap::SetNewIndex
void SetNewIndex(unsigned oldIndex, unsigned newIndex)
Definition VertexElementMap.cpp:63
VertexElementMap::SetDeleted
void SetDeleted(unsigned index)
Definition VertexElementMap.cpp:68
VertexElement::GetFace
VertexElement< ELEMENT_DIM-1, SPACE_DIM > * GetFace(unsigned index) const
Definition VertexElement.cpp:147
VertexMesh::mFaces
std::vector< VertexElement< ELEMENT_DIM - 1, SPACE_DIM > * > mFaces
Definition VertexMesh.hpp:87
VertexMesh::GenerateEdgesFromElements
void GenerateEdgesFromElements(std::vector< VertexElement< ELEMENT_DIM, SPACE_DIM > * > &rElements)
Definition VertexMesh.cpp:569
VertexMesh::Clear
virtual void Clear()
Definition VertexMesh.cpp:788
VertexMesh::mElements
std::vector< VertexElement< ELEMENT_DIM, SPACE_DIM > * > mElements
Definition VertexMesh.hpp:84
CellDivisionInfo::mDaughterLongAxis1
c_vector< double, SPACE_DIM > mDaughterLongAxis1
Definition CellDivisionInfo.hpp:78
CellDivisionInfo::mDaughterLocation2
c_vector< double, SPACE_DIM > mDaughterLocation2
Definition CellDivisionInfo.hpp:81
CellDivisionInfo::mLocation
c_vector< double, SPACE_DIM > mLocation
Definition CellDivisionInfo.hpp:72
CellDivisionInfo::mDaughterLongAxis2
c_vector< double, SPACE_DIM > mDaughterLongAxis2
Definition CellDivisionInfo.hpp:84
CellDivisionInfo::mDaughterLocation1
c_vector< double, SPACE_DIM > mDaughterLocation1
Definition CellDivisionInfo.hpp:75
CellDivisionInfo::mDivisionAxis
c_vector< double, SPACE_DIM > mDivisionAxis
Definition CellDivisionInfo.hpp:87
T1SwapInfo::mLocation
c_vector< double, SPACE_DIM > mLocation
Definition T1SwapInfo.hpp:69
T1SwapInfo::mPreSwapEdge
c_vector< double, SPACE_DIM > mPreSwapEdge
Definition T1SwapInfo.hpp:75
T1SwapInfo::mPostSwapEdge
c_vector< double, SPACE_DIM > mPostSwapEdge
Definition T1SwapInfo.hpp:81
T2SwapInfo::mCellId
unsigned mCellId
Definition T2SwapInfo.hpp:68
T2SwapInfo::mLocation
c_vector< double, SPACE_DIM > mLocation
Definition T2SwapInfo.hpp:71
T3SwapInfo::mLocation
c_vector< double, SPACE_DIM > mLocation
Definition T3SwapInfo.hpp:67