Chaste Commit::6e4f5fe395bca70eb7641cf6e0e87f450383ca5a
MutableVertexMesh.cpp
1/*
2
3Copyright (c) 2005-2026, 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 = this->GetVectorFromAtoB(this->GetElement(elementIndex)->GetNode(intersected_edge)->rGetLocation(), this->mNodes[node_B_index]->rGetLocation());
2238 c_vector<double, SPACE_DIM> vector_A_to_B = this->GetVectorFromAtoB(this->mNodes[node_B_index]->rGetLocation(), this->mNodes[node_A_index]->rGetLocation());
2239 // Insertion of node A splits the intersected edge
2240 if (mTrackMeshOperations)
2241 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), intersected_edge,
2242 norm_2(vector_A_to_B) / norm_2(vector_B_to_node));
2243
2244 if (element_2_index != UINT_MAX)
2245 {
2246 assert(node_B_local_index_in_2 != UINT_MAX);
2247 this->mElements[element_2_index]->DeleteNode(node_B_local_index_in_2);
2248
2249 if (mTrackMeshOperations)
2250 mOperationRecorder.RecordEdgeMergeOperation(this->mElements[element_2_index], node_B_local_index_in_2);
2251 }
2252 if (element_4_index != UINT_MAX)
2253 {
2254 assert(node_A_local_index_in_4 != UINT_MAX);
2255 this->mElements[element_4_index]->DeleteNode(node_A_local_index_in_4);
2256 if (mTrackMeshOperations)
2257 mOperationRecorder.RecordEdgeMergeOperation(this->mElements[element_4_index], node_A_local_index_in_4);
2258 }
2259 }
2260 else
2261 {
2262 assert((intersected_edge + 1) % num_nodes == node_B_local_index_in_3);
2263
2264 // Add node B to element 1 before node A and add node A to element 3 before node B
2265 if (element_1_index != UINT_MAX)
2266 {
2267 assert(node_A_local_index_in_1 != UINT_MAX);
2268 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();
2269 this->mElements[element_1_index]->AddNode(this->mNodes[node_B_index], node_before_A_in_1);
2270 // Insertion of node B creates a new edge
2271 if (mTrackMeshOperations)
2272 mOperationRecorder.RecordNewEdgeOperation(this->mElements[element_1_index], node_before_A_in_1 + 1);
2273 }
2274
2275 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();
2276 this->mElements[element_3_index]->AddNode(this->mNodes[node_A_index], node_before_B_in_3);
2277 c_vector<double, SPACE_DIM> vector_B_to_node = this->GetVectorFromAtoB(this->GetElement(elementIndex)->GetNode(intersected_edge)->rGetLocation(), this->mNodes[node_B_index]->rGetLocation());
2278 c_vector<double, SPACE_DIM> vector_A_to_B = this->GetVectorFromAtoB(this->mNodes[node_B_index]->rGetLocation(), this->mNodes[node_A_index]->rGetLocation());
2279 // Insertion of node A splits the intersected edge
2280 if (mTrackMeshOperations)
2281 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), intersected_edge,
2282 norm_2(vector_A_to_B) / norm_2(vector_B_to_node));
2283
2284 // Remove node A from element 2 and remove node B from element 4
2285 if (element_2_index != UINT_MAX)
2286 {
2287 assert(node_A_local_index_in_2 != UINT_MAX);
2288 this->mElements[element_2_index]->DeleteNode(node_A_local_index_in_2);
2289 if (mTrackMeshOperations)
2290 mOperationRecorder.RecordEdgeMergeOperation(this->mElements[element_2_index], node_A_local_index_in_2);
2291 }
2292 if (element_4_index != UINT_MAX)
2293 {
2294 assert(node_B_local_index_in_4 != UINT_MAX);
2295 this->mElements[element_4_index]->DeleteNode(node_B_local_index_in_4);
2296 if (mTrackMeshOperations)
2297 mOperationRecorder.RecordEdgeMergeOperation(this->mElements[element_4_index], node_B_local_index_in_4);
2298 }
2299 }
2300 if (element_1_index != UINT_MAX)
2301 {
2302 this->mElements[element_1_index]->RebuildEdges();
2303 }
2304 if (element_2_index != UINT_MAX)
2305 {
2306 this->mElements[element_2_index]->RebuildEdges();
2307 }
2308 if (element_3_index != UINT_MAX)
2309 {
2310 this->mElements[element_3_index]->RebuildEdges();
2311 }
2312 if (element_4_index != UINT_MAX)
2313 {
2314 this->mElements[element_4_index]->RebuildEdges();
2315 }
2316 // Set as boundary nodes if intersecting node is a boundary node
2317 if (all_elements_containing_intersecting_node.size() == 2)
2318 {
2319 // Case 1: element 1 is missing
2320 if (elements_containing_intersecting_node.size() == 2)
2321 {
2322 assert(this->mNodes[node_A_index]->IsBoundaryNode() == true);
2323 assert(this->mNodes[node_B_index]->IsBoundaryNode() == false);
2324 this->mNodes[node_B_index]->SetAsBoundaryNode(true);
2325 }
2326 else if (elements_containing_intersecting_node.size() == 1)
2327 {
2328 assert(this->mNodes[node_A_index]->IsBoundaryNode() == true);
2329 assert(this->mNodes[node_B_index]->IsBoundaryNode() == true);
2330 // Case 2: element 2 is missing
2331 if (element_2_index == UINT_MAX)
2332 {
2333 if (intersected_edge == node_B_local_index_in_3)
2334 {
2335 this->mNodes[node_B_index]->SetAsBoundaryNode(false);
2336 }
2337 else
2338 {
2339 this->mNodes[node_A_index]->SetAsBoundaryNode(false);
2340 }
2341 }
2342 // Case 3: element 4 is missing
2343 if (element_4_index == UINT_MAX)
2344 {
2345 if (intersected_edge == node_B_local_index_in_3)
2346 {
2347 this->mNodes[node_A_index]->SetAsBoundaryNode(false);
2348 }
2349 else
2350 {
2351 this->mNodes[node_B_index]->SetAsBoundaryNode(false);
2352 }
2353 }
2354 }
2355 else
2356 {
2357 NEVER_REACHED;
2358 }
2359 }
2360#ifndef NDEBUG
2361 // Case 4: elements 1 and 2 are missing
2362 // Case 5: elements 1 and 4 are missing
2363 else if (all_elements_containing_intersecting_node.size() == 1)
2364 {
2365 if (elements_containing_intersecting_node.size() == 1)
2366 {
2367 assert(this->mNodes[node_A_index]->IsBoundaryNode() == true);
2368 assert(this->mNodes[node_B_index]->IsBoundaryNode() == true);
2369 }
2370 else
2371 {
2372 NEVER_REACHED;
2373 }
2374 }
2375 // Node A is not a boundary node
2376 else
2377 {
2378 assert(this->mNodes[node_A_index]->IsBoundaryNode() == false);
2379 assert(this->mNodes[node_B_index]->IsBoundaryNode() == false);
2380 }
2381#endif
2382 }
2383 else
2384 {
2385 NEVER_REACHED;
2386 }
2387}
2388
2389template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2390void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformT2Swap([[maybe_unused]] VertexElement<ELEMENT_DIM, SPACE_DIM>& rElement)
2391{
2392 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
2393 {
2394 // The given element must be triangular for us to be able to perform a T2 swap on it
2395 assert(rElement.GetNumNodes() == 3);
2396 // Note that we define this vector before setting it, as otherwise the profiling build will break (see #2367)
2397 c_vector<double, SPACE_DIM> new_node_location;
2398 new_node_location = this->GetCentroidOfElement(rElement.GetIndex());
2399 mLastT2SwapLocation = new_node_location;
2400
2401 T2SwapInfo<SPACE_DIM> swap_info;
2402 swap_info.mCellId = rElement.GetIndex();
2403 swap_info.mLocation = new_node_location;
2404 mOperationRecorder.RecordT2Swap(swap_info);
2405
2406 // If this element has no neighbours, delete the element, its nodes and exit function
2407 if (this->GetNeighbouringElementIndices(rElement.GetIndex()).size() == 0)
2408 {
2409 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(0));
2410 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(1));
2411 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(2));
2412
2413 rElement.GetNode(0)->MarkAsDeleted();
2414 rElement.GetNode(1)->MarkAsDeleted();
2415 rElement.GetNode(2)->MarkAsDeleted();
2416
2417 mDeletedElementIndices.push_back(rElement.GetIndex());
2418 rElement.MarkAsDeleted();
2419
2420 return;
2421 }
2422
2423 // Create a new node at the element's centroid; this will be a boundary node if any existing nodes were on the boundary
2424 bool is_node_on_boundary = false;
2425 for (unsigned i = 0; i < 3; i++)
2426 {
2427 if (rElement.GetNode(i)->IsBoundaryNode())
2428 {
2429 is_node_on_boundary = true;
2430 break;
2431 }
2432 }
2433 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(GetNumNodes(), new_node_location, is_node_on_boundary));
2434 Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
2435
2436 std::set<unsigned> neigh_indices;
2437
2438 // Loop over each of the three nodes contained in rElement
2439 for (unsigned i = 0; i < 3; i++)
2440 {
2441 // For each node, find the set of other elements containing it
2442 Node<SPACE_DIM>* p_node = rElement.GetNode(i);
2443
2444 std::set<unsigned> containing_elements = p_node->rGetContainingElementIndices();
2445 containing_elements.erase(rElement.GetIndex());
2446 // For each of these elements...
2447 for (std::set<unsigned>::iterator elem_iter = containing_elements.begin(); elem_iter != containing_elements.end(); ++elem_iter)
2448 {
2449 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_this_elem = this->GetElement(*elem_iter);
2450
2451 neigh_indices.insert(p_this_elem->GetIndex());
2452
2453 // ...throw an exception if the element is triangular...
2454 if (p_this_elem->GetNumNodes() < 4)
2455 {
2456 EXCEPTION("One of the neighbours of a small triangular element is also a triangle - dealing with this has not been implemented yet");
2457 }
2458
2459 // ...otherwise, replace p_node with p_new_node unless this has already happened (in which case, delete p_node from the element)
2460 if (p_this_elem->GetNodeLocalIndex(new_node_global_index) == UINT_MAX)
2461 {
2462 p_this_elem->ReplaceNode(p_node, p_new_node);
2463 }
2464 else
2465 {
2466 std::vector<unsigned> old_ids;
2467 std::pair<unsigned, unsigned> node_pair;
2468 if (mTrackMeshOperations)
2469 {
2470 for (unsigned k = 0; k < p_this_elem->GetNumEdges(); ++k)
2471 {
2472 old_ids.push_back(p_this_elem->GetEdge(k)->GetIndex());
2473 }
2474 node_pair.first = p_this_elem->GetNodeLocalIndex(new_node_global_index);
2475 node_pair.second = p_this_elem->GetNodeLocalIndex(p_node->GetIndex());
2476 }
2477 p_this_elem->DeleteNode(p_this_elem->GetNodeLocalIndex(p_node->GetIndex()));
2478 if (mTrackMeshOperations)
2479 {
2480 mOperationRecorder.RecordNodeMergeOperation(old_ids, p_this_elem, node_pair, true);
2481 }
2482 }
2483 }
2484 }
2485
2486 // We also have to mark pElement, pElement->GetNode(0), pElement->GetNode(1), and pElement->GetNode(2) as deleted
2487 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(0));
2488 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(1));
2489 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(2));
2490
2491 rElement.GetNode(0)->MarkAsDeleted();
2492 rElement.GetNode(1)->MarkAsDeleted();
2493 rElement.GetNode(2)->MarkAsDeleted();
2494
2495 mDeletedElementIndices.push_back(rElement.GetIndex());
2496 rElement.MarkAsDeleted();
2497
2498 for (unsigned i : neigh_indices)
2499 {
2500 this->GetElement(i)->RebuildEdges();
2501 }
2502 }
2503 else
2504 {
2505 NEVER_REACHED;
2506 }
2507}
2508
2509template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2510void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformT3Swap(
2511 [[maybe_unused]] Node<SPACE_DIM>* pNode, [[maybe_unused]] unsigned elementIndex)
2512{
2513 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
2514 {
2515
2516 assert(pNode->IsBoundaryNode());
2517
2518 // Get element
2519 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elementIndex);
2520 unsigned num_nodes = p_element->GetNumNodes();
2521
2522 // Store the index of the elements containing the intersecting node
2523 std::set<unsigned> elements_containing_intersecting_node = pNode->rGetContainingElementIndices();
2524
2525 // Get the local index of the node in the intersected element after which the new node is to be added
2526 unsigned node_A_local_index = this->GetLocalIndexForElementEdgeClosestToPoint(pNode->rGetLocation(), elementIndex);
2527
2528 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2529 c_vector<double, SPACE_DIM> node_location;
2530 node_location = pNode->rGetModifiableLocation();
2531
2532 // Get the nodes at either end of the edge to be divided
2533 unsigned vertexA_index = p_element->GetNodeGlobalIndex(node_A_local_index);
2534 unsigned vertexB_index = p_element->GetNodeGlobalIndex((node_A_local_index + 1) % num_nodes);
2535
2536 // Check these nodes are also boundary nodes if this fails then the elements have become concave and you need a smaller timestep
2537 if (!this->mNodes[vertexA_index]->IsBoundaryNode() || !this->mNodes[vertexB_index]->IsBoundaryNode())
2538 {
2539 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.");
2540 }
2541
2542 // Get the nodes at either end of the edge to be divided and calculate intersection
2543 c_vector<double, SPACE_DIM> vertexA = p_element->GetNodeLocation(node_A_local_index);
2544 c_vector<double, SPACE_DIM> vertexB = p_element->GetNodeLocation((node_A_local_index + 1) % num_nodes);
2545 c_vector<double, SPACE_DIM> vector_a_to_point = this->GetVectorFromAtoB(vertexA, node_location);
2546
2547 c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
2548
2549 c_vector<double, SPACE_DIM> edge_ab_unit_vector = vector_a_to_b / norm_2(vector_a_to_b);
2550 c_vector<double, SPACE_DIM> intersection = vertexA + edge_ab_unit_vector * inner_prod(vector_a_to_point, edge_ab_unit_vector);
2551
2552 // Store the location of the T3 swap, the location of the intersection with the edge
2554 // is called (see #2401) - we should correct this in these cases!
2555 T3SwapInfo<SPACE_DIM> swap_info;
2556 swap_info.mLocation = intersection;
2557 mOperationRecorder.RecordT3Swap(swap_info);
2558
2559 // For rebuilding edges
2560 std::set<unsigned> rebuilt_elements;
2561 if (pNode->GetNumContainingElements() == 1)
2562 {
2563 // Get the index of the element containing the intersecting node
2564 unsigned intersecting_element_index = *elements_containing_intersecting_node.begin();
2565
2566 // Get element
2567 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_intersecting_element = this->GetElement(intersecting_element_index);
2568
2569 //Edge map before swap
2570 const unsigned num_edges = p_intersecting_element->GetNumEdges();
2571 std::vector<unsigned> old_ids(num_edges);
2572 for (unsigned i=0; i<num_edges; ++i)
2573 {
2574 old_ids[i] = p_intersecting_element->GetEdge(i)->GetIndex();
2575 }
2576
2577 unsigned local_index = p_intersecting_element->GetNodeLocalIndex(pNode->GetIndex());
2578 unsigned next_node = p_intersecting_element->GetNodeGlobalIndex((local_index + 1) % (p_intersecting_element->GetNumNodes()));
2579 unsigned previous_node = p_intersecting_element->GetNodeGlobalIndex((local_index + p_intersecting_element->GetNumNodes() - 1) % (p_intersecting_element->GetNumNodes()));
2580
2581 // Check to see if the nodes adjacent to the intersecting node are contained in the intersected element between vertices A and B
2582 if (next_node == vertexA_index || previous_node == vertexA_index || next_node == vertexB_index || previous_node == vertexB_index)
2583 {
2584 unsigned common_vertex_index;
2585
2586 if (next_node == vertexA_index || previous_node == vertexA_index)
2587 {
2588 common_vertex_index = vertexA_index;
2589 }
2590 else
2591 {
2592 common_vertex_index = vertexB_index;
2593 }
2594
2595 assert(this->mNodes[common_vertex_index]->GetNumContainingElements() > 1);
2596
2597 std::set<unsigned> elements_containing_common_vertex = this->mNodes[common_vertex_index]->rGetContainingElementIndices();
2598 std::set<unsigned>::const_iterator it = elements_containing_common_vertex.begin();
2599 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*it);
2600 it++;
2601 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*it);
2602
2603 // Find the number and indices of common vertices between element_1 and element_2
2604 unsigned num_common_vertices = 0;
2605 std::vector<unsigned> common_vertex_indices;
2606 for (unsigned i = 0; i < p_element_common_1->GetNumNodes(); i++)
2607 {
2608 for (unsigned j = 0; j < p_element_common_2->GetNumNodes(); j++)
2609 {
2610 if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
2611 {
2612 num_common_vertices++;
2613 common_vertex_indices.push_back(p_element_common_1->GetNodeGlobalIndex(i));
2614 }
2615 }
2616 }
2617
2618 if (num_common_vertices == 1 || this->mNodes[common_vertex_index]->GetNumContainingElements() > 2)
2619 {
2620 /*
2621 * This is the situation here.
2622 *
2623 * From To
2624 * _ _
2625 * | <--- |
2626 * | /\ |\
2627 * | / \ | \
2628 * _|/____\ _|__\
2629 *
2630 * The edge goes from vertexA--vertexB to vertexA--pNode--vertexB
2631 */
2632
2633 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2634 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2635
2636 // Move original node
2637 pNode->rGetModifiableLocation() = intersection;
2638
2639 // Record the edge split data
2640 const double a_to_b_length = norm_2(vector_a_to_b);
2641 c_vector<double, SPACE_DIM> vector_a_to_node = this->GetVectorFromAtoB(vertexA, intersection);
2642 const double a_to_node_length = norm_2(vector_a_to_node);
2643
2644 // Add the moved nodes to the element (this also updates the node)
2645 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2646
2647 if (mTrackMeshOperations)
2648 mOperationRecorder.RecordEdgeSplitOperation(p_element, node_A_local_index, a_to_node_length / a_to_b_length);
2649 rebuilt_elements.insert(p_element->GetIndex());
2650 // Check the nodes are updated correctly
2651 assert(pNode->GetNumContainingElements() == 2);
2652 }
2653 else if (num_common_vertices == 2)
2654 {
2655 // The two elements must have an edge in common. Find whether the common edge is the same as the
2656 // edge that is merged onto.
2657
2658 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))
2659 {
2660 /*
2661 * Due to a previous T3 swap the situation looks like this.
2662 *
2663 * pNode
2664 * \ |\ /
2665 * \ | \ /
2666 * \_______|__\/
2667 * /A | B
2668 * / \
2669 *
2670 * A T3 Swap would merge pNode onto an edge of its own element.
2671 * We prevent this by just removing pNode. By doing this we also avoid the
2672 * intersecting element to be concave.
2673 */
2674
2675 // Obtain necessary data for event recording
2676 const unsigned downstream_index = (local_index + p_intersecting_element->GetNumNodes() - 1) % (p_intersecting_element->GetNumNodes());
2677
2678 // Delete pNode in the intersecting element
2679 p_intersecting_element->DeleteNode(local_index);
2680
2681 // We record node deletion as node merging.
2682 if (mTrackMeshOperations)
2683 {
2684 mOperationRecorder.RecordNodeMergeOperation(old_ids, p_intersecting_element, std::pair<unsigned, unsigned>(downstream_index, local_index));
2685 }
2686 rebuilt_elements.insert(p_intersecting_element->GetIndex());
2687
2688 // Mark all three nodes as deleted
2689 pNode->MarkAsDeleted();
2690 mDeletedNodeIndices.push_back(pNode->GetIndex());
2691 }
2692 else
2693 {
2694 /*
2695 * This is the situation here.
2696 *
2697 * C is common_vertex D is the other one.
2698 *
2699 * From To
2700 * _ D _
2701 * | <--- |
2702 * | /\ |\
2703 * C|/ \ | \
2704 * _|____\ _|__\
2705 *
2706 * The edge goes from vertexC--vertexB to vertexC--pNode--vertexD
2707 * then vertex B is removed as it is no longer needed.
2708 */
2709
2710 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2711 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2712
2713 // Move original node
2714 pNode->rGetModifiableLocation() = intersection;
2715
2716 // Replace common_vertex with the moved node (this also updates the nodes)
2717 this->GetElement(elementIndex)->ReplaceNode(this->mNodes[common_vertex_index], pNode);
2718
2719 // Remove common_vertex
2720 unsigned common_vertex_local_index = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(common_vertex_index);
2721 this->GetElement(intersecting_element_index)->DeleteNode(common_vertex_local_index);
2722 assert(this->mNodes[common_vertex_index]->GetNumContainingElements() == 0);
2723
2724 // Record edge merging in the intersecting element
2725 if (mTrackMeshOperations)
2726 mOperationRecorder.RecordEdgeMergeOperation(p_intersecting_element, common_vertex_local_index);
2727 rebuilt_elements.insert(p_intersecting_element->GetIndex());
2728
2729 this->mNodes[common_vertex_index]->MarkAsDeleted();
2730 mDeletedNodeIndices.push_back(common_vertex_index);
2731
2732 // Check the nodes are updated correctly
2733 assert(pNode->GetNumContainingElements() == 2);
2734 }
2735 }
2736 else if (num_common_vertices == 4)
2737 {
2738 /*
2739 * The two elements share edges CA and BD due to previous swaps but not the edge AB
2740 *
2741 * From To
2742 * D___ D___
2743 * | |
2744 * B|\ |
2745 * | \ |
2746 * | / |
2747 * A|/ |
2748 * C_|__ C_|__
2749 *
2750 * We just remove the intersecting node as well as vertices A and B.
2751 */
2752
2753 // Delete node A and B in the intersected element
2754 this->GetElement(elementIndex)->DeleteNode(node_A_local_index);
2755 if (mTrackMeshOperations)
2756 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(elementIndex), node_A_local_index);
2757 unsigned node_B_local_index = this->GetElement(elementIndex)->GetNodeLocalIndex(vertexB_index);
2758 this->GetElement(elementIndex)->DeleteNode(node_B_local_index);
2759 if (mTrackMeshOperations)
2760 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(elementIndex), node_B_local_index);
2761
2762 // Delete nodes A and B in the intersecting element
2763 unsigned node_A_local_index_intersecting_element = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(vertexA_index);
2764 this->GetElement(intersecting_element_index)->DeleteNode(node_A_local_index_intersecting_element);
2765 if (mTrackMeshOperations)
2766 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(intersecting_element_index), node_A_local_index_intersecting_element);
2767 unsigned node_B_local_index_intersecting_element = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(vertexB_index);
2768 this->GetElement(intersecting_element_index)->DeleteNode(node_B_local_index_intersecting_element);
2769 if (mTrackMeshOperations)
2770 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(intersecting_element_index), node_B_local_index_intersecting_element);
2771
2772 // Delete pNode in the intersecting element
2773 unsigned p_node_local_index = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(pNode->GetIndex());
2774 this->GetElement(intersecting_element_index)->DeleteNode(p_node_local_index);
2775 if (mTrackMeshOperations)
2776 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(intersecting_element_index), p_node_local_index);
2777
2778 rebuilt_elements.insert(elementIndex);
2779 rebuilt_elements.insert(intersecting_element_index);
2780 // Mark all three nodes as deleted
2781 pNode->MarkAsDeleted();
2782 mDeletedNodeIndices.push_back(pNode->GetIndex());
2783 this->mNodes[vertexA_index]->MarkAsDeleted();
2784 mDeletedNodeIndices.push_back(vertexA_index);
2785 this->mNodes[vertexB_index]->MarkAsDeleted();
2786 mDeletedNodeIndices.push_back(vertexB_index);
2787 }
2788 else
2789 {
2790 // This can't happen as nodes can't be on the internal edge of 2 elements.
2791 NEVER_REACHED;
2792 }
2793 }
2794 else
2795 {
2796 /*
2797 * From To
2798 * ____ _______
2799 * / \
2800 * /\ ^ / \
2801 * / \ |
2802 *
2803 * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
2804 */
2805
2806 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2807 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2808 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
2809
2810 // Move original node
2811 pNode->rGetModifiableLocation() = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
2812
2813 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
2814 c_vector<double, SPACE_DIM> new_node_location;
2815 new_node_location = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
2816
2817 // Add new node which will always be a boundary node
2818 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
2819
2820 // Add the moved node (this also updates the node)
2821 // and record the split of the edge vertexA--vertexB due to insertion of pNode
2822 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2823 // Recording edge split
2824 c_vector<double, SPACE_DIM> vector_a_to_node = this->GetVectorFromAtoB(pNode->rGetLocation(), vertexA);
2825 if (mTrackMeshOperations)
2826 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index,
2827 norm_2(vector_a_to_node) / norm_2(vector_a_to_b));
2828
2829 // Add the new node to the element (this also updates the node)
2830 // and record the split of the edge vertexA--pNode by insertion of the new node
2831 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
2832
2833 // Recording edge split
2834 c_vector<double, SPACE_DIM> vector_a_to_new_node = this->GetVectorFromAtoB(new_node_location, vertexA);
2835 if (mTrackMeshOperations)
2836 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index,
2837 norm_2(vector_a_to_new_node) / norm_2(vector_a_to_node));
2838 // New node must be between vertexA and pNode
2839 assert(norm_2(vector_a_to_new_node) / norm_2(vector_a_to_node) < 1);
2840
2841 // Add the new node to the original element containing pNode (this also updates the node)
2842 const unsigned new_node_local_index = this->GetElement(intersecting_element_index)->GetNodeLocalIndex(pNode->GetIndex());
2843 this->GetElement(intersecting_element_index)->AddNode(this->mNodes[new_node_global_index], new_node_local_index);
2844 if (mTrackMeshOperations)
2845 {
2846 mOperationRecorder.RecordNewEdgeOperation(this->GetElement(intersecting_element_index), new_node_local_index);
2847 }
2848
2849 // The nodes must have been updated correctly
2850 assert(pNode->GetNumContainingElements() == 2);
2851 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
2852
2853 rebuilt_elements.insert(elementIndex);
2854 rebuilt_elements.insert(intersecting_element_index);
2855 }
2856 }
2857 else if (pNode->GetNumContainingElements() == 2)
2858 {
2859 // Find the nodes contained in elements containing the intersecting node
2860 std::set<unsigned>::const_iterator it = elements_containing_intersecting_node.begin();
2861
2862 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_1 = this->GetElement(*it);
2863 unsigned num_nodes_elem_1 = p_element_1->GetNumNodes();
2864 it++;
2865
2866 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_2 = this->GetElement(*it);
2867 unsigned num_nodes_elem_2 = p_element_2->GetNumNodes();
2868
2869 unsigned node_global_index = pNode->GetIndex();
2870
2871 unsigned local_index_1 = p_element_1->GetNodeLocalIndex(node_global_index);
2872 unsigned next_node_1 = p_element_1->GetNodeGlobalIndex((local_index_1 + 1) % num_nodes_elem_1);
2873 unsigned previous_node_1 = p_element_1->GetNodeGlobalIndex((local_index_1 + num_nodes_elem_1 - 1) % num_nodes_elem_1);
2874
2875 unsigned local_index_2 = p_element_2->GetNodeLocalIndex(node_global_index);
2876 unsigned next_node_2 = p_element_2->GetNodeGlobalIndex((local_index_2 + 1) % num_nodes_elem_2);
2877 unsigned previous_node_2 = p_element_2->GetNodeGlobalIndex((local_index_2 + num_nodes_elem_2 - 1) % num_nodes_elem_2);
2878
2879 // Check to see if the nodes adjacent to the intersecting node are contained in the intersected element between vertices A and B
2880 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))
2881 {
2882 /*
2883 * Here we have
2884 * __
2885 * /| /
2886 * __ / | --> ___/
2887 * \ | \
2888 * \|__ \
2889 *
2890 * Where the node on the left has overlapped the edge A B
2891 *
2892 * Move p_node to the intersection on A B and merge AB and p_node
2893 */
2894
2895 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
2896 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
2897 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
2898
2899 // Check they are all boundary nodes
2900 assert(pNode->IsBoundaryNode());
2901 assert(this->mNodes[vertexA_index]->IsBoundaryNode());
2902 assert(this->mNodes[vertexB_index]->IsBoundaryNode());
2903
2904 // Move p_node to the intersection with the edge AB
2905 pNode->rGetModifiableLocation() = intersection;
2906 pNode->SetAsBoundaryNode(false);
2907
2908 // Add pNode to the intersected element
2909 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
2910
2911 // Record this as edge split
2912 const double a_to_b_length = norm_2(vector_a_to_b);
2913 c_vector<double, SPACE_DIM> vector_a_to_p = this->GetVectorFromAtoB(intersection, vertexA);
2914 const double split_ratio = norm_2(vector_a_to_p) / a_to_b_length;
2915 if (mTrackMeshOperations)
2916 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, split_ratio);
2917 rebuilt_elements.insert(elementIndex);
2918
2919 // Remove vertex A from elements and record this as edge merge operation
2920 std::set<unsigned> elements_containing_vertex_A = this->mNodes[vertexA_index]->rGetContainingElementIndices();
2921 for (std::set<unsigned>::const_iterator iter = elements_containing_vertex_A.begin();
2922 iter != elements_containing_vertex_A.end();
2923 iter++)
2924 {
2925 const unsigned this_vertexA_local_index = this->GetElement(*iter)->GetNodeLocalIndex(vertexA_index);
2926 this->GetElement(*iter)->DeleteNode(this_vertexA_local_index);
2927 if (mTrackMeshOperations)
2928 {
2929 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(*iter), this_vertexA_local_index);
2930 }
2931 rebuilt_elements.insert(this->GetElement(*iter)->GetIndex());
2932 }
2933
2934 // Remove vertex A from the mesh
2935 assert(this->mNodes[vertexA_index]->GetNumContainingElements() == 0);
2936 this->mNodes[vertexA_index]->MarkAsDeleted();
2937 mDeletedNodeIndices.push_back(vertexA_index);
2938
2939 // Remove vertex B from elements and record this as edge merge operation
2940 std::set<unsigned> elements_containing_vertex_B = this->mNodes[vertexB_index]->rGetContainingElementIndices();
2941 for (std::set<unsigned>::const_iterator iter = elements_containing_vertex_B.begin();
2942 iter != elements_containing_vertex_B.end();
2943 iter++)
2944 {
2945 const unsigned this_vertexB_local_index = this->GetElement(*iter)->GetNodeLocalIndex(vertexB_index);
2946 this->GetElement(*iter)->DeleteNode(this_vertexB_local_index);
2947 if (mTrackMeshOperations)
2948 {
2949 mOperationRecorder.RecordEdgeMergeOperation(this->GetElement(*iter), this_vertexB_local_index);
2950 }
2951 rebuilt_elements.insert(this->GetElement(*iter)->GetIndex());
2952 }
2953
2954 // Remove vertex B from the mesh
2955 assert(this->mNodes[vertexB_index]->GetNumContainingElements() == 0);
2956 this->mNodes[vertexB_index]->MarkAsDeleted();
2957 mDeletedNodeIndices.push_back(vertexB_index);
2958 }
2959 else
2960 {
2961 if (next_node_1 == vertexA_index || previous_node_1 == vertexA_index || next_node_2 == vertexA_index || previous_node_2 == vertexA_index)
2962 {
2963 // Get elements containing vertexA_index (the common vertex)
2964
2965 assert(this->mNodes[vertexA_index]->GetNumContainingElements() > 1);
2966
2967 std::set<unsigned> elements_containing_vertex_A = this->mNodes[vertexA_index]->rGetContainingElementIndices();
2968 std::set<unsigned>::const_iterator iter = elements_containing_vertex_A.begin();
2969 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*iter);
2970 iter++;
2971 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*iter);
2972
2973 // Calculate the number of common vertices between element_1 and element_2
2974 unsigned num_common_vertices = 0;
2975 for (unsigned i = 0; i < p_element_common_1->GetNumNodes(); i++)
2976 {
2977 for (unsigned j = 0; j < p_element_common_2->GetNumNodes(); j++)
2978 {
2979 if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
2980 {
2981 num_common_vertices++;
2982 }
2983 }
2984 }
2985
2986 if (num_common_vertices == 1 || this->mNodes[vertexA_index]->GetNumContainingElements() > 2)
2987 {
2988 /*
2989 * From To
2990 * _ B _ B
2991 * | <--- |
2992 * | /|\ |\
2993 * | / | \ | \
2994 * | / | \ |\ \
2995 * _|/___|___\ _|_\_\
2996 * A A
2997 *
2998 * The edge goes from vertexA--vertexB to vertexA--pNode--new_node--vertexB
2999 */
3000
3001 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3002 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3003 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3004
3005 // Move original node and change to non-boundary node
3006 pNode->rGetModifiableLocation() = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3007 pNode->SetAsBoundaryNode(false);
3008
3009 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
3010 c_vector<double, SPACE_DIM> new_node_location;
3011 new_node_location = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3012
3013 // Add new node, which will always be a boundary node
3014 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
3015
3016 // Add the moved nodes to the element (this also updates the node)
3017 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
3018
3019 const double a_to_b_length = norm_2(vector_a_to_b);
3020 c_vector<double, SPACE_DIM> a_to_new_node_vector = this->GetVectorFromAtoB(vertexA, new_node_location);
3021 const double a_to_new_length = norm_2(a_to_new_node_vector);
3022 const double ratio_new_node = a_to_new_length / a_to_b_length;
3023 if (mTrackMeshOperations)
3024 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_new_node);
3025 rebuilt_elements.insert(elementIndex);
3026
3027 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3028
3029 c_vector<double, SPACE_DIM> a_to_p_vector = this->GetVectorFromAtoB(vertexA, pNode->rGetLocation());
3030 const double ratio_p_node = norm_2(a_to_p_vector) / a_to_new_length;
3031 assert(ratio_p_node <= 1);
3032 if (mTrackMeshOperations)
3033 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_p_node);
3034
3035 // Add the new nodes to the original elements containing pNode (this also updates the node)
3036 if (next_node_1 == previous_node_2)
3037 {
3038 const unsigned insertion_local_index = (local_index_1 + p_element_1->GetNumNodes() - 1) % (p_element_1->GetNumNodes());
3039 p_element_1->AddNode(this->mNodes[new_node_global_index], insertion_local_index);
3040 if (mTrackMeshOperations)
3041 {
3042 mOperationRecorder.RecordNewEdgeOperation(p_element_1, insertion_local_index);
3043 }
3044 rebuilt_elements.insert(p_element_1->GetIndex());
3045 }
3046 else
3047 {
3048 assert(next_node_2 == previous_node_1);
3049 const unsigned insertion_local_index = (local_index_2 + p_element_2->GetNumNodes() - 1) % (p_element_2->GetNumNodes());
3050 p_element_2->AddNode(this->mNodes[new_node_global_index], insertion_local_index);
3051 if (mTrackMeshOperations)
3052 {
3053 mOperationRecorder.RecordNewEdgeOperation(p_element_2, insertion_local_index);
3054 }
3055 rebuilt_elements.insert(p_element_2->GetIndex());
3056 }
3057
3058 // Check the nodes are updated correctly
3059 assert(pNode->GetNumContainingElements() == 3);
3060 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
3061 }
3062 else if (num_common_vertices == 2)
3063 {
3064 /*
3065 * From To
3066 * _ B _ B
3067 * |<--- |
3068 * | /|\ |\
3069 * |/ | \ | \
3070 * | | \ |\ \
3071 * _|__|___\ _|_\_\
3072 * A A
3073 *
3074 * The edge goes from vertexA--vertexB to vertexA--pNode--new_node--vertexB
3075 * then vertexA is removed
3076 */
3077
3078 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3079 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3080 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3081
3082 // Move original node and change to non-boundary node
3083 pNode->rGetModifiableLocation() = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3084 pNode->SetAsBoundaryNode(false);
3085
3086 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
3087 c_vector<double, SPACE_DIM> new_node_location;
3088 new_node_location = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3089
3090 // Add new node, which will always be a boundary node
3091 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
3092
3093 // Add the moved nodes to the element (this also updates the node)
3094 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
3095
3096 const double a_to_b_length = norm_2(vector_a_to_b);
3097 c_vector<double, SPACE_DIM> a_to_new_node_vector = this->GetVectorFromAtoB(vertexA, new_node_location);
3098 const double a_to_new_length = norm_2(a_to_new_node_vector);
3099 const double ratio_new_node = a_to_new_length / a_to_b_length;
3100 if (mTrackMeshOperations)
3101 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_new_node);
3102
3103 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3104
3105 c_vector<double, SPACE_DIM> a_to_p_vector = this->GetVectorFromAtoB(vertexA, pNode->rGetLocation());
3106 const double ratio_p_node = norm_2(a_to_p_vector) / a_to_new_length;
3107 assert(ratio_p_node <= 1);
3108 if (mTrackMeshOperations)
3109 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_p_node);
3110
3111 rebuilt_elements.insert(elementIndex);
3112 // Add the new nodes to the original elements containing pNode (this also updates the node)
3113 if (next_node_1 == previous_node_2)
3114 {
3115 const unsigned insertion_local_index = (local_index_1 + p_element_1->GetNumNodes() - 1) % (p_element_1->GetNumNodes());
3116 p_element_1->AddNode(this->mNodes[new_node_global_index], insertion_local_index);
3117 if (mTrackMeshOperations)
3118 {
3119 mOperationRecorder.RecordNewEdgeOperation(p_element_1, insertion_local_index);
3120 }
3121 rebuilt_elements.insert(p_element_1->GetIndex());
3122 }
3123 else
3124 {
3125 assert(next_node_2 == previous_node_1);
3126 const unsigned insertion_local_index = (local_index_2 + p_element_2->GetNumNodes() - 1) % (p_element_2->GetNumNodes());
3127 p_element_2->AddNode(this->mNodes[new_node_global_index], insertion_local_index);
3128 if (mTrackMeshOperations)
3129 {
3130 mOperationRecorder.RecordNewEdgeOperation(p_element_2, insertion_local_index);
3131 }
3132 rebuilt_elements.insert(p_element_2->GetIndex());
3133 }
3134
3135 // Remove vertex A from the mesh
3136 const unsigned local_index_1 = p_element_common_1->GetNodeLocalIndex(vertexA_index);
3137 const unsigned local_index_2 = p_element_common_2->GetNodeLocalIndex(vertexA_index);
3138 p_element_common_1->DeleteNode(local_index_1);
3139 if (mTrackMeshOperations)
3140 {
3141 mOperationRecorder.RecordEdgeMergeOperation(p_element_common_1, local_index_1);
3142 }
3143 p_element_common_2->DeleteNode(local_index_2);
3144 if (mTrackMeshOperations)
3145 {
3146 mOperationRecorder.RecordEdgeMergeOperation(p_element_common_2, local_index_2);
3147 }
3148 assert(this->mNodes[vertexA_index]->GetNumContainingElements() == 0);
3149 rebuilt_elements.insert(p_element_common_1->GetIndex());
3150 rebuilt_elements.insert(p_element_common_2->GetIndex());
3151
3152 this->mNodes[vertexA_index]->MarkAsDeleted();
3153 mDeletedNodeIndices.push_back(vertexA_index);
3154
3155 // Check the nodes are updated correctly
3156 assert(pNode->GetNumContainingElements() == 3);
3157 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
3158 }
3159 else
3160 {
3161 // This can't happen as nodes can't be on the internal edge of two elements
3162 NEVER_REACHED;
3163 }
3164 }
3165 else if (next_node_1 == vertexB_index || previous_node_1 == vertexB_index || next_node_2 == vertexB_index || previous_node_2 == vertexB_index)
3166 {
3167 // Get elements containing vertexB_index (the common vertex)
3168
3169 assert(this->mNodes[vertexB_index]->GetNumContainingElements() > 1);
3170
3171 std::set<unsigned> elements_containing_vertex_B = this->mNodes[vertexB_index]->rGetContainingElementIndices();
3172 std::set<unsigned>::const_iterator iter = elements_containing_vertex_B.begin();
3173 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*iter);
3174 iter++;
3175 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*iter);
3176
3177 // Calculate the number of common vertices between element_1 and element_2
3178 unsigned num_common_vertices = 0;
3179 for (unsigned i = 0; i < p_element_common_1->GetNumNodes(); i++)
3180 {
3181 for (unsigned j = 0; j < p_element_common_2->GetNumNodes(); j++)
3182 {
3183 if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
3184 {
3185 num_common_vertices++;
3186 }
3187 }
3188 }
3189
3190 if (num_common_vertices == 1 || this->mNodes[vertexB_index]->GetNumContainingElements() > 2)
3191 {
3192 /*
3193 * From To
3194 * _B_________ _B____
3195 * |\ | / | / /
3196 * | \ | / |/ /
3197 * | \ | / | /
3198 * | \|/ |/
3199 * _| <--- _|
3200 * A
3201 *
3202 * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
3203 */
3204
3205 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3206 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3207 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3208
3209 // Move original node and change to non-boundary node
3210 pNode->rGetModifiableLocation() = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3211 pNode->SetAsBoundaryNode(false);
3212
3213 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
3214 c_vector<double, SPACE_DIM> new_node_location;
3215 new_node_location = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3216
3217 // Add new node which will always be a boundary node
3218 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
3219
3220 // Add the moved nodes to the element (this also updates the node)
3221 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3222 const double a_to_b_length = norm_2(vector_a_to_b);
3223 c_vector<double, SPACE_DIM> a_to_p_vector = this->GetVectorFromAtoB(vertexA, pNode->rGetLocation());
3224 const double ratio_p_node = norm_2(a_to_p_vector) / a_to_b_length;
3225 assert(ratio_p_node <= 1);
3226 if (mTrackMeshOperations)
3227 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_p_node);
3228
3229 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
3230 c_vector<double, SPACE_DIM> a_to_new_node_vector = this->GetVectorFromAtoB(vertexA, new_node_location);
3231 const double a_to_new_length = norm_2(a_to_new_node_vector);
3232 const double ratio_new_node = a_to_new_length / norm_2(a_to_p_vector);
3233 assert(ratio_new_node <= 1);
3234 if (mTrackMeshOperations)
3235 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_new_node);
3236 rebuilt_elements.insert(elementIndex);
3237 // Add the new nodes to the original elements containing pNode (this also updates the node)
3238 if (next_node_1 == previous_node_2)
3239 {
3240 p_element_2->AddNode(this->mNodes[new_node_global_index], local_index_2);
3241 if (mTrackMeshOperations)
3242 mOperationRecorder.RecordNewEdgeOperation(p_element_2, local_index_2);
3243 rebuilt_elements.insert(p_element_2->GetIndex());
3244 }
3245 else
3246 {
3247 assert(next_node_2 == previous_node_1);
3248 p_element_1->AddNode(this->mNodes[new_node_global_index], local_index_1);
3249 if (mTrackMeshOperations)
3250 mOperationRecorder.RecordNewEdgeOperation(p_element_1, local_index_1);
3251 rebuilt_elements.insert(p_element_1->GetIndex());
3252 }
3253
3254 // Check the nodes are updated correctly
3255 assert(pNode->GetNumContainingElements() == 3);
3256 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
3257 }
3258 else if (num_common_vertices == 2)
3259 {
3260 /*
3261 * From To
3262 * _B_______ _B____
3263 * | | / | / /
3264 * | | / |/ /
3265 * |\ | / | /
3266 * | \|/ |/
3267 * _| <--- _|
3268 * A
3269 *
3270 * The edge goes from vertexA--vertexB to vertexA--new_node--pNode--vertexB
3271 * then vertexB is removed
3272 */
3273
3274 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3275 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3276 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3277
3278 // Move original node and change to non-boundary node
3279 pNode->rGetModifiableLocation() = intersection + 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3280 pNode->SetAsBoundaryNode(false);
3281
3282 // Note that we define this vector before setting it as otherwise the profiling build will break (see #2367)
3283 c_vector<double, SPACE_DIM> new_node_location;
3284 new_node_location = intersection - 0.5 * mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3285
3286 // Add new node which will always be a boundary node
3287 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
3288
3289 // Add the moved nodes to the element (this also updates the node)
3290 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3291 const double a_to_b_length = norm_2(vector_a_to_b);
3292 c_vector<double, SPACE_DIM> a_to_p_vector = this->GetVectorFromAtoB(vertexA, pNode->rGetLocation());
3293 const double ratio_p_node = norm_2(a_to_p_vector) / a_to_b_length;
3294 assert(ratio_p_node <= 1);
3295 if (mTrackMeshOperations)
3296 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_p_node);
3297
3298 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_global_index], node_A_local_index);
3299 c_vector<double, SPACE_DIM> a_to_new_node_vector = this->GetVectorFromAtoB(vertexA, new_node_location);
3300 const double a_to_new_length = norm_2(a_to_new_node_vector);
3301 const double ratio_new_node = a_to_new_length / norm_2(a_to_p_vector);
3302 assert(ratio_new_node <= 1);
3303 if (mTrackMeshOperations)
3304 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, ratio_new_node);
3305 rebuilt_elements.insert(elementIndex);
3306 // Add the new nodes to the original elements containing pNode (this also updates the node)
3307 if (next_node_1 == previous_node_2)
3308 {
3309 p_element_2->AddNode(this->mNodes[new_node_global_index], local_index_2);
3310 if (mTrackMeshOperations)
3311 mOperationRecorder.RecordNewEdgeOperation(p_element_2, local_index_2);
3312 rebuilt_elements.insert(p_element_2->GetIndex());
3313 }
3314 else
3315 {
3316 assert(next_node_2 == previous_node_1);
3317 p_element_1->AddNode(this->mNodes[new_node_global_index], local_index_1);
3318 if (mTrackMeshOperations)
3319 mOperationRecorder.RecordNewEdgeOperation(p_element_1, local_index_1);
3320 rebuilt_elements.insert(p_element_1->GetIndex());
3321 }
3322
3323 // Remove vertex B from the mesh
3324 const unsigned local_index_1 = p_element_common_1->GetNodeLocalIndex(vertexB_index);
3325 const unsigned local_index_2 = p_element_common_2->GetNodeLocalIndex(vertexB_index);
3326 p_element_common_1->DeleteNode(local_index_1);
3327 if (mTrackMeshOperations)
3328 {
3329 mOperationRecorder.RecordEdgeMergeOperation(p_element_common_1, local_index_1);
3330 }
3331 p_element_common_2->DeleteNode(local_index_2);
3332 if (mTrackMeshOperations)
3333 {
3334 mOperationRecorder.RecordEdgeMergeOperation(p_element_common_2, local_index_2);
3335 }
3336 rebuilt_elements.insert(p_element_common_1->GetIndex());
3337 rebuilt_elements.insert(p_element_common_2->GetIndex());
3338
3339 assert(this->mNodes[vertexB_index]->GetNumContainingElements() == 0);
3340
3341 this->mNodes[vertexB_index]->MarkAsDeleted();
3342 mDeletedNodeIndices.push_back(vertexB_index);
3343
3344 // Check the nodes are updated correctly
3345 assert(pNode->GetNumContainingElements() == 3);
3346 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
3347 }
3348 else
3349 {
3350 // This can't happen as nodes can't be on the internal edge of two elements
3351 NEVER_REACHED;
3352 }
3353 }
3354 else
3355 {
3356 /*
3357 * From To
3358 * _____ _______
3359 * / | \
3360 * /|\ ^ / | \
3361 * / | \ |
3362 *
3363 * The edge goes from vertexA--vertexB to vertexA--new_node_1--pNode--new_node_2--vertexB
3364 */
3365
3366 // Check whether the intersection location fits into the edge and update distances and vertex positions afterwards.
3367 intersection = this->WidenEdgeOrCorrectIntersectionLocationIfNecessary(vertexA_index, vertexB_index, intersection);
3368 edge_ab_unit_vector = this->GetPreviousEdgeGradientOfElementAtNode(p_element, (node_A_local_index + 1) % num_nodes);
3369
3370 // Move original node and change to non-boundary node
3371 pNode->rGetModifiableLocation() = intersection;
3372 pNode->SetAsBoundaryNode(false);
3373
3374 c_vector<double, SPACE_DIM> new_node_1_location;
3375 new_node_1_location = intersection - mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3376 c_vector<double, SPACE_DIM> new_node_2_location;
3377 new_node_2_location = intersection + mCellRearrangementRatio * mCellRearrangementThreshold * edge_ab_unit_vector;
3378
3379 // Add new nodes which will always be boundary nodes
3380 unsigned new_node_1_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_1_location[0], new_node_1_location[1]));
3381 unsigned new_node_2_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_2_location[0], new_node_2_location[1]));
3382
3383 // Add the moved and new nodes to the element (this also updates the node)
3384 const double a_to_b_length = norm_2(vector_a_to_b);
3385 const c_vector<double, SPACE_DIM> a_to_node2_vector = this->GetVectorFromAtoB(new_node_2_location, vertexA);
3386 const double a_to_node2_length = norm_2(a_to_node2_vector);
3387 const c_vector<double, SPACE_DIM> a_to_nodeP_vector = this->GetVectorFromAtoB(intersection, vertexA);
3388 const double a_to_nodeP_length = norm_2(a_to_nodeP_vector);
3389 const c_vector<double, SPACE_DIM> a_to_node1_vector = this->GetVectorFromAtoB(new_node_1_location, vertexA);
3390 const double a_to_node1_length = norm_2(a_to_node1_vector);
3391
3392 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_2_global_index], node_A_local_index);
3393 if (mTrackMeshOperations)
3394 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, a_to_node2_length / a_to_b_length);
3395 this->GetElement(elementIndex)->AddNode(pNode, node_A_local_index);
3396 if (mTrackMeshOperations)
3397 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, a_to_nodeP_length / a_to_node2_length);
3398 this->GetElement(elementIndex)->AddNode(this->mNodes[new_node_1_global_index], node_A_local_index);
3399 if (mTrackMeshOperations)
3400 mOperationRecorder.RecordEdgeSplitOperation(this->GetElement(elementIndex), node_A_local_index, a_to_node1_length / a_to_nodeP_length);
3401 rebuilt_elements.insert(elementIndex);
3402 // Add the new nodes to the original elements containing pNode (this also updates the node)
3403 if (next_node_1 == previous_node_2)
3404 {
3405 const unsigned inserted_local_index_1 = (local_index_1 + p_element_1->GetNumNodes() - 1) % (p_element_1->GetNumNodes());
3406 p_element_1->AddNode(this->mNodes[new_node_2_global_index], inserted_local_index_1);
3407 if (mTrackMeshOperations)
3408 {
3409 mOperationRecorder.RecordNewEdgeOperation(p_element_1, inserted_local_index_1);
3410 }
3411 p_element_2->AddNode(this->mNodes[new_node_1_global_index], local_index_2);
3412 if (mTrackMeshOperations)
3413 {
3414 mOperationRecorder.RecordNewEdgeOperation(p_element_2, local_index_2);
3415 }
3416 }
3417 else
3418 {
3419 assert(next_node_2 == previous_node_1);
3420 p_element_1->AddNode(this->mNodes[new_node_1_global_index], local_index_1);
3421 if (mTrackMeshOperations)
3422 {
3423 mOperationRecorder.RecordNewEdgeOperation(p_element_1, local_index_1);
3424 }
3425 const unsigned inserted_local_index_2 = (local_index_2 + p_element_2->GetNumNodes() - 1) % (p_element_2->GetNumNodes());
3426 p_element_2->AddNode(this->mNodes[new_node_2_global_index], inserted_local_index_2);
3427 if (mTrackMeshOperations)
3428 {
3429 mOperationRecorder.RecordNewEdgeOperation(p_element_2, local_index_2);
3430 }
3431 }
3432 rebuilt_elements.insert(p_element_1->GetIndex());
3433 rebuilt_elements.insert(p_element_2->GetIndex());
3434
3435 // Check the nodes are updated correctly
3436 assert(pNode->GetNumContainingElements() == 3);
3437 assert(this->mNodes[new_node_1_global_index]->GetNumContainingElements() == 2);
3438 assert(this->mNodes[new_node_2_global_index]->GetNumContainingElements() == 2);
3439 }
3440 }
3441 }
3442 else
3443 {
3444 EXCEPTION("Trying to merge a node, contained in more than 2 elements, into another element, this is not possible with the vertex mesh.");
3445 }
3446 for (unsigned i : rebuilt_elements)
3447 {
3448 this->GetElement(i)->RebuildEdges();
3449 }
3450 }
3451 else
3452 {
3453 NEVER_REACHED;
3454 }
3455}
3456
3457template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3458void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformVoidRemoval(
3459 [[maybe_unused]] Node<SPACE_DIM>* pNodeA,
3460 [[maybe_unused]] Node<SPACE_DIM>* pNodeB,
3461 [[maybe_unused]] Node<SPACE_DIM>* pNodeC)
3462{
3463 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
3464 {
3465 // Calculate void centroid
3466 c_vector<double, SPACE_DIM> nodes_midpoint = pNodeA->rGetLocation()
3467 + this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation()) / 3.0
3468 + this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeC->rGetLocation()) / 3.0;
3469 /*
3470 * In two steps, merge nodes A, B and C into a single node. This is implemented in such a way that
3471 * the ordering of their indices does not matter.
3472 */
3473 PerformNodeMerge(pNodeA, pNodeB);
3474 Node<SPACE_DIM>* p_merged_node = pNodeB;
3475
3476 if (pNodeB->IsDeleted())
3477 {
3478 p_merged_node = pNodeA;
3479 }
3480
3481 PerformNodeMerge(pNodeC, p_merged_node);
3482 if (p_merged_node->IsDeleted())
3483 {
3484 p_merged_node = pNodeC;
3485 }
3486
3487 p_merged_node->rGetModifiableLocation() = nodes_midpoint;
3488
3489 // Tag remaining node as non-boundary
3490 p_merged_node->SetAsBoundaryNode(false);
3491
3492 // Remove the deleted nodes and re-index
3493 RemoveDeletedNodes();
3494 }
3495 else
3496 {
3497 NEVER_REACHED;
3498 }
3499}
3500
3501template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3502void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::HandleHighOrderJunctions(Node<SPACE_DIM>* pNodeA, Node<SPACE_DIM>* pNodeB)
3503{
3504 unsigned node_a_rank = pNodeA->rGetContainingElementIndices().size();
3505 unsigned node_b_rank = pNodeB->rGetContainingElementIndices().size();
3506
3507 if ((node_a_rank > 3) && (node_b_rank > 3))
3508 {
3509 // The code can't handle this case
3510 EXCEPTION("Both nodes involved in a swap event are contained in more than three elements");
3511 }
3512 else // the rosette degree should increase in this case
3513 {
3514 assert(node_a_rank > 3 || node_b_rank > 3);
3515 this->PerformRosetteRankIncrease(pNodeA, pNodeB);
3516 this->RemoveDeletedNodes();
3517 }
3518}
3519
3520template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3521void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformRosetteRankIncrease(
3522 [[maybe_unused]] Node<SPACE_DIM>* pNodeA, [[maybe_unused]] Node<SPACE_DIM>* pNodeB)
3523{
3524 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
3525 {
3526 /*
3527 * One of the nodes will have 3 containing element indices, the other
3528 * will have at least four. We first identify which node is which.
3529 */
3530
3531 unsigned node_a_index = pNodeA->GetIndex();
3532 unsigned node_b_index = pNodeB->GetIndex();
3533
3534 unsigned node_a_rank = pNodeA->rGetContainingElementIndices().size();
3535 unsigned node_b_rank = pNodeB->rGetContainingElementIndices().size();
3536
3537 unsigned lo_rank_index = (node_a_rank < node_b_rank) ? node_a_index : node_b_index;
3538 unsigned hi_rank_index = (node_a_rank < node_b_rank) ? node_b_index : node_a_index;
3539
3540 // Get pointers to the nodes, sorted by index
3541 Node<SPACE_DIM>* p_lo_rank_node = this->GetNode(lo_rank_index);
3542 Node<SPACE_DIM>* p_hi_rank_node = this->GetNode(hi_rank_index);
3543
3544 // Find the sets of elements containing each of the nodes, sorted by index
3545 std::set<unsigned> lo_rank_elem_indices = p_lo_rank_node->rGetContainingElementIndices();
3546 std::set<unsigned> hi_rank_elem_indices = p_hi_rank_node->rGetContainingElementIndices();
3547
3574 for (std::set<unsigned>::const_iterator it = lo_rank_elem_indices.begin();
3575 it != lo_rank_elem_indices.end();
3576 ++it)
3577 {
3578 // Find the local index of lo_rank_node in this element
3579 unsigned lo_rank_local_index = this->mElements[*it]->GetNodeLocalIndex(lo_rank_index);
3580 assert(lo_rank_local_index < UINT_MAX); // double check this element contains lo_rank_node
3581
3582 /*
3583 * If this element already contains the hi_rank_node, we are in the situation of elements
3584 * C and D above, so we just remove lo_rank_node.
3585 *
3586 * Otherwise, we are in element E, so we must replace lo_rank_node with high-rank node,
3587 * and remove it from mNodes.
3588 *
3589 * We can check whether hi_rank_node is in this element using the set::count() method.
3590 */
3591
3592 if (hi_rank_elem_indices.count(*it) > 0)
3593 {
3594 // For node merge recording
3595 std::vector<unsigned> old_ids(this->mElements[*it]->GetNumEdges(),0);
3596
3597 for (unsigned i = 0; i < old_ids.size(); ++i)
3598 {
3599 old_ids[i] = this->mElements[*it]->GetEdge(i)->GetIndex();
3600 }
3601 // Delete lo_rank_node from current element
3602 this->mElements[*it]->DeleteNode(lo_rank_local_index);
3603
3604 // Record edge shrinkage
3605 const unsigned hi_rank_local_index = this->mElements[*it]->GetNodeLocalIndex(hi_rank_index);
3606 if (mTrackMeshOperations)
3607 {
3608 mOperationRecorder.RecordNodeMergeOperation(old_ids,
3609 this->mElements[*it],
3610 std::pair<unsigned, unsigned>(hi_rank_local_index, lo_rank_local_index));
3611 }
3612 }
3613 else
3614 {
3615 // Update lo_rank_node with all information (including index and location) of hi_rank_node
3616 this->mElements[*it]->UpdateNode(lo_rank_local_index, p_hi_rank_node);
3617 }
3618 this->mElements[*it]->RebuildEdges();
3619 }
3620
3621 // Tidy up the mesh by ensuring the global instance of lo_rank_node is deleted
3622 assert(!(this->mNodes[lo_rank_index]->IsDeleted()));
3623 this->mNodes[lo_rank_index]->MarkAsDeleted();
3624 this->mDeletedNodeIndices.push_back(lo_rank_index);
3625
3626 for (unsigned i : lo_rank_elem_indices)
3627 {
3628 this->GetElement(i)->RebuildEdges();
3629 }
3630 for (unsigned i : hi_rank_elem_indices)
3631 {
3632 this->GetElement(i)->RebuildEdges();
3633 }
3634 }
3635 else
3636 {
3637 NEVER_REACHED;
3638 }
3639}
3640
3641template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3642void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformProtorosetteResolution(
3643 [[maybe_unused]] Node<SPACE_DIM>* pProtorosetteNode)
3644{
3645 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
3646 {
3647 // Double check we are dealing with a protorosette
3648 assert(pProtorosetteNode->rGetContainingElementIndices().size() == 4);
3649
3650 // Get random number (0, 1, 2 or 3), as the resolution axis is assumed to be random
3651 unsigned random_elem_increment = RandomNumberGenerator::Instance()->randMod(4);
3652
3653 // Find global indices of elements around the protorosette node
3654 std::set<unsigned> protorosette_node_containing_elem_indices = pProtorosetteNode->rGetContainingElementIndices();
3655
3656 // Select random element by advancing iterator a random number times
3657 std::set<unsigned>::const_iterator elem_index_iter(protorosette_node_containing_elem_indices.begin());
3658 advance(elem_index_iter, random_elem_increment);
3659
3687 /*
3688 * We need to find the global indices of elements B, C and D. We do this with set intersections.
3689 */
3690
3691 unsigned elem_a_idx = *elem_index_iter;
3692 unsigned elem_b_idx = UINT_MAX;
3693 unsigned elem_c_idx = UINT_MAX;
3694 unsigned elem_d_idx = UINT_MAX;
3695
3696 // Get pointer to element we've chosen at random (element A)
3697 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_a = this->GetElement(elem_a_idx);
3698
3699 // Get all necessary info about element A and the protorosette node
3700 unsigned num_nodes_elem_a = p_elem_a->GetNumNodes();
3701 unsigned protorosette_node_global_idx = pProtorosetteNode->GetIndex();
3702 unsigned protorosette_node_local_idx = p_elem_a->GetNodeLocalIndex(protorosette_node_global_idx);
3703
3704 // Find global indices of previous (cw) and next (ccw) nodes, locally, from the protorosette node, in element A
3705 unsigned prev_node_global_idx = p_elem_a->GetNodeGlobalIndex((protorosette_node_local_idx + num_nodes_elem_a - 1) % num_nodes_elem_a);
3706 unsigned next_node_global_idx = p_elem_a->GetNodeGlobalIndex((protorosette_node_local_idx + 1) % num_nodes_elem_a);
3707
3708 // Get the set of elements the previous and next nodes are contained in
3709 Node<SPACE_DIM>* p_prev_node = this->GetNode(prev_node_global_idx);
3710 Node<SPACE_DIM>* p_next_node = this->GetNode(next_node_global_idx);
3711 std::set<unsigned> prev_node_elem_indices = p_prev_node->rGetContainingElementIndices();
3712 std::set<unsigned> next_node_elem_indices = p_next_node->rGetContainingElementIndices();
3713
3714 // Perform set intersections with the set of element indices which the protorosette node is contained in
3715 std::set<unsigned> intersection_with_prev;
3716 std::set<unsigned> intersection_with_next;
3717
3718 // This intersection should contain just global indices for elements A and B
3719 std::set_intersection(protorosette_node_containing_elem_indices.begin(),
3720 protorosette_node_containing_elem_indices.end(),
3721 prev_node_elem_indices.begin(),
3722 prev_node_elem_indices.end(),
3723 std::inserter(intersection_with_prev, intersection_with_prev.begin()));
3724
3725 // This intersection should contain just global indices for elements A and D
3726 std::set_intersection(protorosette_node_containing_elem_indices.begin(),
3727 protorosette_node_containing_elem_indices.end(),
3728 next_node_elem_indices.begin(),
3729 next_node_elem_indices.end(),
3730 std::inserter(intersection_with_next, intersection_with_next.begin()));
3731
3732 assert(intersection_with_prev.size() == 2);
3733 assert(intersection_with_next.size() == 2);
3734
3735 // Get global index of element B
3736 if (*intersection_with_prev.begin() != elem_a_idx)
3737 {
3738 elem_b_idx = *(intersection_with_prev.begin());
3739 }
3740 else
3741 {
3742 elem_b_idx = *(++(intersection_with_prev.begin()));
3743 }
3744 assert(elem_b_idx < UINT_MAX);
3745
3746 // Get global index of element D
3747 if (*intersection_with_next.begin() != elem_a_idx)
3748 {
3749 elem_d_idx = *(intersection_with_next.begin());
3750 }
3751 else
3752 {
3753 elem_d_idx = *(++(intersection_with_next.begin()));
3754 }
3755 assert(elem_d_idx < UINT_MAX);
3756
3757 // By elimination, the remaining unassigned index in the original set must be global index of element C
3758 for (elem_index_iter = protorosette_node_containing_elem_indices.begin();
3759 elem_index_iter != protorosette_node_containing_elem_indices.end();
3760 ++elem_index_iter)
3761 {
3762 if ((*elem_index_iter != elem_a_idx) && (*elem_index_iter != elem_b_idx) && (*elem_index_iter != elem_d_idx))
3763 {
3764 elem_c_idx = *elem_index_iter;
3765 }
3766 }
3767 assert(elem_c_idx < UINT_MAX);
3768
3783 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_b = this->GetElement(elem_b_idx);
3784 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_c = this->GetElement(elem_c_idx);
3785 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_d = this->GetElement(elem_d_idx);
3786
3787 double swap_distance = (this->mCellRearrangementRatio) * (this->mCellRearrangementThreshold);
3788
3789 // Get normalized vectors to centre of elements A and B from protorosette node
3790 c_vector<double, SPACE_DIM> node_to_elem_a_centre = this->GetCentroidOfElement(elem_a_idx) - pProtorosetteNode->rGetLocation();
3791 node_to_elem_a_centre /= norm_2(node_to_elem_a_centre);
3792
3793 c_vector<double, SPACE_DIM> node_to_elem_c_centre = this->GetCentroidOfElement(elem_c_idx) - pProtorosetteNode->rGetLocation();
3794 node_to_elem_c_centre /= norm_2(node_to_elem_c_centre);
3795
3796 // Calculate new node locations
3797 c_vector<double, SPACE_DIM> new_location_of_protorosette_node = pProtorosetteNode->rGetLocation() + (0.5 * swap_distance) * node_to_elem_a_centre;
3798 c_vector<double, SPACE_DIM> location_of_new_node = pProtorosetteNode->rGetLocation() + (0.5 * swap_distance) * node_to_elem_c_centre;
3799
3800 // Move protorosette node to new location
3801 pProtorosetteNode->rGetModifiableLocation() = new_location_of_protorosette_node;
3802
3803 // Create new node in correct location
3804 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(this->GetNumNodes(), location_of_new_node, false));
3805 Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
3806
3817 unsigned local_idx_elem_b = p_elem_b->GetNodeLocalIndex(protorosette_node_global_idx);
3818 local_idx_elem_b = (local_idx_elem_b + p_elem_b->GetNumNodes() - 1) % p_elem_b->GetNumNodes();
3819 unsigned local_idx_elem_c = p_elem_c->GetNodeLocalIndex(protorosette_node_global_idx);
3820 unsigned local_idx_elem_d = p_elem_d->GetNodeLocalIndex(protorosette_node_global_idx);
3821
3822 p_elem_b->AddNode(p_new_node, local_idx_elem_b);
3823 if (mTrackMeshOperations)
3824 mOperationRecorder.RecordNewEdgeOperation(p_elem_b, local_idx_elem_b);
3825 p_elem_c->AddNode(p_new_node, local_idx_elem_c);
3826 p_elem_d->AddNode(p_new_node, local_idx_elem_d);
3827 if (mTrackMeshOperations)
3828 mOperationRecorder.RecordNewEdgeOperation(p_elem_d, local_idx_elem_d);
3829
3830 // All that is left is to remove the original protorosette node from element C
3831 p_elem_c->DeleteNode(p_elem_c->GetNodeLocalIndex(protorosette_node_global_idx));
3832 }
3833 else
3834 {
3835 NEVER_REACHED;
3836 }
3837}
3838
3839template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
3840void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformRosetteRankDecrease(
3841 [[maybe_unused]] Node<SPACE_DIM>* pRosetteNode)
3842{
3843 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
3844 {
3845 unsigned rosette_rank = pRosetteNode->rGetContainingElementIndices().size();
3846
3847 // Double check we're dealing with a rosette
3848 assert(rosette_rank > 4);
3849
3850 // Get random number in [0, 1, ..., n) where n is rank of rosette, as the resolution axis is assumed to be random
3851 unsigned random_elem_increment = RandomNumberGenerator::Instance()->randMod(rosette_rank);
3852
3853 // Find global indices of elements around the protorosette node
3854 std::set<unsigned> rosette_node_containing_elem_indices = pRosetteNode->rGetContainingElementIndices();
3855
3856 // Select random element by advancing iterator a random number times
3857 std::set<unsigned>::const_iterator elem_index_iter(rosette_node_containing_elem_indices.begin());
3858 advance(elem_index_iter, random_elem_increment);
3859
3898 /*
3899 * We need to find the global indices of elements N and P. We do this with set intersections.
3900 */
3901
3902 // Get the vertex element S (which we randomly selected)
3903 unsigned elem_s_idx = *elem_index_iter;
3904 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_s = this->GetElement(elem_s_idx);
3905
3906 unsigned elem_n_idx = UINT_MAX;
3907 unsigned elem_p_idx = UINT_MAX;
3908
3909 // Get all necessary info about element S and the rosette node
3910 unsigned num_nodes_elem_s = p_elem_s->GetNumNodes();
3911 unsigned rosette_node_global_idx = pRosetteNode->GetIndex();
3912 unsigned rosette_node_local_idx = p_elem_s->GetNodeLocalIndex(rosette_node_global_idx);
3913
3914 // Find global indices of previous (cw) and next (ccw) nodes, locally, from the rosette node, in element S
3915 unsigned prev_node_global_idx = p_elem_s->GetNodeGlobalIndex((rosette_node_local_idx + num_nodes_elem_s - 1) % num_nodes_elem_s);
3916 unsigned next_node_global_idx = p_elem_s->GetNodeGlobalIndex((rosette_node_local_idx + 1) % num_nodes_elem_s);
3917
3918 // Get the set of elements that the previous and next nodes are contained in
3919 Node<SPACE_DIM>* p_prev_node = this->GetNode(prev_node_global_idx);
3920 Node<SPACE_DIM>* p_next_node = this->GetNode(next_node_global_idx);
3921 std::set<unsigned> prev_node_elem_indices = p_prev_node->rGetContainingElementIndices();
3922 std::set<unsigned> next_node_elem_indices = p_next_node->rGetContainingElementIndices();
3923
3924 // Perform set intersections with the set of element indices that the rosette node is contained in
3925 std::set<unsigned> intersection_with_prev;
3926 std::set<unsigned> intersection_with_next;
3927
3928 // This intersection should contain just global indices for elements S and N
3929 std::set_intersection(rosette_node_containing_elem_indices.begin(),
3930 rosette_node_containing_elem_indices.end(),
3931 prev_node_elem_indices.begin(),
3932 prev_node_elem_indices.end(),
3933 std::inserter(intersection_with_prev, intersection_with_prev.begin()));
3934
3935 // This intersection should contain just global indices for elements S and P
3936 std::set_intersection(rosette_node_containing_elem_indices.begin(),
3937 rosette_node_containing_elem_indices.end(),
3938 next_node_elem_indices.begin(),
3939 next_node_elem_indices.end(),
3940 std::inserter(intersection_with_next, intersection_with_next.begin()));
3941
3942 assert(intersection_with_prev.size() == 2);
3943 assert(intersection_with_next.size() == 2);
3944
3945 // Get global index of element N
3946 if (*intersection_with_prev.begin() != elem_s_idx)
3947 {
3948 elem_n_idx = *intersection_with_prev.begin();
3949 }
3950 else
3951 {
3952 elem_n_idx = *(++(intersection_with_prev.begin()));
3953 }
3954 assert(elem_n_idx < UINT_MAX);
3955
3956 // Get global index of element P
3957 if (*intersection_with_next.begin() != elem_s_idx)
3958 {
3959 elem_p_idx = *intersection_with_next.begin();
3960 }
3961 else
3962 {
3963 elem_p_idx = *(++(intersection_with_next.begin()));
3964 }
3965 assert(elem_p_idx < UINT_MAX);
3966
3977 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_n = this->GetElement(elem_p_idx);
3978 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_elem_p = this->GetElement(elem_n_idx);
3979
3980 double swap_distance = (this->mCellRearrangementRatio) * (this->mCellRearrangementThreshold);
3981
3982 // Calculate location of new node
3983 c_vector<double, 2> node_to_selected_elem = this->GetCentroidOfElement(elem_s_idx) - pRosetteNode->rGetLocation();
3984 node_to_selected_elem /= norm_2(node_to_selected_elem);
3985 c_vector<double, 2> new_node_location = pRosetteNode->rGetLocation() + (swap_distance * node_to_selected_elem);
3986
3987 // Create new node in correct location
3988 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(this->GetNumNodes(), new_node_location, false));
3989 Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
3990
4002 // Add new node, and remove rosette node, from element S
4003 unsigned node_local_idx_in_elem_s = p_elem_s->GetNodeLocalIndex(rosette_node_global_idx);
4004 p_elem_s->AddNode(p_new_node, node_local_idx_in_elem_s);
4005 p_elem_s->DeleteNode(node_local_idx_in_elem_s);
4006
4007 // Add new node to element N
4008 unsigned node_local_idx_in_elem_n = p_elem_n->GetNodeLocalIndex(rosette_node_global_idx);
4009 node_local_idx_in_elem_n = (node_local_idx_in_elem_n + p_elem_n->GetNumNodes() - 1) % p_elem_n->GetNumNodes();
4010 p_elem_n->AddNode(p_new_node, node_local_idx_in_elem_n);
4011 if (mTrackMeshOperations)
4012 {
4013 mOperationRecorder.RecordNewEdgeOperation(p_elem_n, node_local_idx_in_elem_n);
4014 }
4015 // Add new node to element P
4016 unsigned node_local_idx_in_elem_p = p_elem_p->GetNodeLocalIndex(rosette_node_global_idx);
4017 p_elem_p->AddNode(p_new_node, node_local_idx_in_elem_p);
4018 if (mTrackMeshOperations)
4019 {
4020 mOperationRecorder.RecordNewEdgeOperation(p_elem_p, node_local_idx_in_elem_p);
4021 }
4022 }
4023}
4024
4025template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
4026void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::CheckForRosettes()
4027{
4028 if constexpr (ELEMENT_DIM == 2 && SPACE_DIM == 2)
4029 {
4038 // Vectors to store the nodes that need resolution events
4039 std::vector<Node<SPACE_DIM>*> protorosette_nodes;
4040 std::vector<Node<SPACE_DIM>*> rosette_nodes;
4041
4042 // First loop in which we populate these vectors
4043 unsigned num_nodes = this->GetNumAllNodes();
4044 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
4045 {
4046 Node<SPACE_DIM>* current_node = this->GetNode(node_idx);
4047 unsigned node_rank = current_node->rGetContainingElementIndices().size();
4048
4049 if (node_rank < 4)
4050 {
4051 // Nothing to do if the node is not high-rank
4052 continue;
4053 }
4054 else if (node_rank == 4)
4055 {
4056 // For protorosette nodes, we check against a random number to decide if resolution is necessary
4057 if (mProtorosetteResolutionProbabilityPerTimestep >= RandomNumberGenerator::Instance()->ranf())
4058 {
4059 protorosette_nodes.push_back(current_node);
4060 }
4061 }
4062 else // if (node_rank > 4)
4063 {
4064 // For rosette nodes, we check against a random number to decide if resolution is necessary
4065 if (mRosetteResolutionProbabilityPerTimestep >= RandomNumberGenerator::Instance()->ranf())
4066 {
4067 rosette_nodes.push_back(current_node);
4068 }
4069 }
4070 }
4071
4079 // First, resolve any protorosettes
4080 for (unsigned node_idx = 0; node_idx < protorosette_nodes.size(); node_idx++)
4081 {
4082 Node<SPACE_DIM>* current_node = protorosette_nodes[node_idx];
4083
4084 // Verify that node has not been marked for deletion, and that it is still contained in four elements
4085 assert(!(current_node->IsDeleted()));
4086 assert(current_node->rGetContainingElementIndices().size() == 4);
4087
4088 // Perform protorosette resolution
4089 this->PerformProtorosetteResolution(current_node);
4090 }
4091
4092 // Finally, resolve any rosettes
4093 for (unsigned node_idx = 0; node_idx < rosette_nodes.size(); node_idx++)
4094 {
4095 Node<SPACE_DIM>* current_node = rosette_nodes[node_idx];
4096
4097 // Verify that node has not been marked for deletion, and that it is still contained in at least four elements
4098 assert(!(current_node->IsDeleted()));
4099 assert(current_node->rGetContainingElementIndices().size() > 4);
4100
4101 // Perform protorosette resolution
4102 this->PerformRosetteRankDecrease(current_node);
4103 }
4104 }
4105 else
4106 {
4107 NEVER_REACHED;
4108 }
4109}
4110
4111template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
4112c_vector<double, 2> MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::WidenEdgeOrCorrectIntersectionLocationIfNecessary(
4113 unsigned indexA, unsigned indexB, c_vector<double, 2> intersection)
4114{
4124 c_vector<double, SPACE_DIM> vertexA = this->GetNode(indexA)->rGetLocation();
4125 c_vector<double, SPACE_DIM> vertexB = this->GetNode(indexB)->rGetLocation();
4126 c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
4127
4128 if (norm_2(vector_a_to_b) < 4.0*mCellRearrangementRatio*mCellRearrangementThreshold)
4129 {
4130 WARNING("Trying to merge a node onto an edge which is too small.");
4131
4132 c_vector<double, SPACE_DIM> centre_a_and_b = vertexA + 0.5*vector_a_to_b;
4133
4134 vertexA = centre_a_and_b - 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*vector_a_to_b/norm_2(vector_a_to_b);
4135 ChastePoint<SPACE_DIM> vertex_A_point(vertexA);
4136 SetNode(indexA, vertex_A_point);
4137
4138 vertexB = centre_a_and_b + 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*vector_a_to_b/norm_2(vector_a_to_b);
4139 ChastePoint<SPACE_DIM> vertex_B_point(vertexB);
4140 SetNode(indexB, vertex_B_point);
4141
4142 intersection = centre_a_and_b;
4143 }
4144
4145 // Reset distances
4146 vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
4147 c_vector<double, 2> edge_ab_unit_vector = vector_a_to_b/norm_2(vector_a_to_b);
4148
4164 if (norm_2(intersection - vertexA) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
4165 {
4166 intersection = vertexA + 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
4167 }
4168 if (norm_2(intersection - vertexB) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
4169 {
4170 intersection = vertexB - 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
4171 }
4172 return intersection;
4173}
4174
4175template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
4176void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetMeshOperationTracking(const bool track)
4177{
4178 mTrackMeshOperations = track;
4179 if (track)
4180 {
4181 mOperationRecorder.SetEdgeHelper(&(this->mEdgeHelper));
4182 }
4183}
4184
4185template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
4191
4192// Explicit instantiation
4193template class MutableVertexMesh<1, 1>;
4194template class MutableVertexMesh<1, 2>;
4195template class MutableVertexMesh<1, 3>;
4196template class MutableVertexMesh<2, 2>;
4197template class MutableVertexMesh<2, 3>;
4198template class MutableVertexMesh<3, 3>;
4199
4200// Serialization for Boost >= 1.36
4201#include "SerializationExportWrapperForCpp.hpp"
4202EXPORT_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:3521
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:3458
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:4112
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:2510
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:2390
MutableVertexMesh::PerformProtorosetteResolution
void PerformProtorosetteResolution(Node< SPACE_DIM > *pProtorosetteNode)
Definition MutableVertexMesh.cpp:3642
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:4176
MutableVertexMesh::HandleHighOrderJunctions
virtual void HandleHighOrderJunctions(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
Definition MutableVertexMesh.cpp:3502
MutableVertexMesh::SetCellRearrangementRatio
void SetCellRearrangementRatio(double cellRearrangementRatio)
Definition MutableVertexMesh.cpp:221
MutableVertexMesh::PerformRosetteRankDecrease
void PerformRosetteRankDecrease(Node< SPACE_DIM > *pRosetteNode)
Definition MutableVertexMesh.cpp:3840
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:4186
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:570
VertexMesh::Clear
virtual void Clear()
Definition VertexMesh.cpp:801
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