Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
MutableMesh.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include <map>
37#include <cstring>
38
39#include "MutableMesh.hpp"
40#include "OutputFileHandler.hpp"
41
42//Jonathan Shewchuk's triangle and Hang Si's tetgen
43#define REAL double
44#define VOID void
45#include "triangle.h"
46#include "tetgen.h"
47#undef REAL
48#undef VOID
49
50
51template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
53 : mAddedNodes(false)
54{
56}
57
58template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
60{
61 this->mMeshChangesDuringSimulation = true;
62 Clear();
63 for (unsigned index=0; index<nodes.size(); index++)
64 {
65 Node<SPACE_DIM>* p_temp_node = nodes[index];
66 this->mNodes.push_back(p_temp_node);
67 }
68 mAddedNodes = true;
69 NodeMap node_map(nodes.size());
70 ReMesh(node_map);
71}
72
73template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
78
79template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
81{
82 if (mDeletedNodeIndices.empty())
83 {
84 pNewNode->SetIndex(this->mNodes.size());
85 this->mNodes.push_back(pNewNode);
86 }
87 else
88 {
89 unsigned index = mDeletedNodeIndices.back();
90 pNewNode->SetIndex(index);
91 mDeletedNodeIndices.pop_back();
92 delete this->mNodes[index];
93 this->mNodes[index] = pNewNode;
94 }
95 mAddedNodes = true;
96 return pNewNode->GetIndex();
97}
98
99template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
101{
102 unsigned new_elt_index;
103
104 if (mDeletedElementIndices.empty())
105 {
106 new_elt_index = this->mElements.size();
107 this->mElements.push_back(pNewElement);
108 pNewElement->ResetIndex(new_elt_index);
109 }
110 else
111 {
112 unsigned index = mDeletedElementIndices.back();
113 pNewElement->ResetIndex(index);
114 mDeletedElementIndices.pop_back();
115 delete this->mElements[index];
116 this->mElements[index] = pNewElement;
117 }
118
119 return pNewElement->GetIndex();
120}
121
122
123template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
125{
126 mDeletedElementIndices.clear();
127 mDeletedBoundaryElementIndices.clear();
128 mDeletedNodeIndices.clear();
129 mAddedNodes = false;
130
132}
133
134template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
136{
137 return this->mBoundaryElements.size() - mDeletedBoundaryElementIndices.size();
138}
139
140template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
142{
143 return this->mElements.size() - mDeletedElementIndices.size();
144}
145
146template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
148{
149 return this->mNodes.size() - mDeletedNodeIndices.size();
150}
151
158template<>
159void MutableMesh<1, 1>::RescaleMeshFromBoundaryNode(ChastePoint<1> updatedPoint, unsigned boundaryNodeIndex)
160{
161 assert(this->GetNode(boundaryNodeIndex)->IsBoundaryNode());
162 double scaleFactor = updatedPoint[0] / this->GetNode(boundaryNodeIndex)->GetPoint()[0];
163 double temp;
164 for (unsigned i=0; i < boundaryNodeIndex+1; i++)
165 {
166 temp = scaleFactor * this->mNodes[i]->GetPoint()[0];
167 ChastePoint<1> newPoint(temp);
168 this->mNodes[i]->SetPoint(newPoint);
169 }
170 this->RefreshMesh();
172
173template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
176 bool concreteMove)
177{
178 this->mNodes[index]->SetPoint(point);
179
180 if (concreteMove)
182 for (typename Node<SPACE_DIM>::ContainingBoundaryElementIterator it = this->mNodes[index]->ContainingBoundaryElementsBegin();
183 it != this->mNodes[index]->ContainingBoundaryElementsEnd();
184 ++it)
185 {
186 try
187 {
188 this->GetBoundaryElement(*it)->CalculateWeightedDirection(this->mBoundaryElementWeightedDirections[ (*it) ],
189 this->mBoundaryElementJacobianDeterminants[ (*it) ]);
190 }
191 catch (Exception&)
192 {
193 EXCEPTION("Moving node caused a boundary element to have a non-positive Jacobian determinant");
194 }
195 }
196 for (typename Node<SPACE_DIM>::ContainingElementIterator it = this->mNodes[index]->ContainingElementsBegin();
197 it != this->mNodes[index]->ContainingElementsEnd();
198 ++it)
199 {
200 if (ELEMENT_DIM == SPACE_DIM)
202 try
203 {
204 this->GetElement(*it)->CalculateInverseJacobian(this->mElementJacobians[ (*it) ],
205 this->mElementJacobianDeterminants[ (*it) ],
206 this->mElementInverseJacobians[ (*it) ]);
207 }
208 catch (Exception&)
209 {
210 EXCEPTION("Moving node caused an element to have a non-positive Jacobian determinant");
211 }
212 }
213 else
214 {
215 c_vector<double,SPACE_DIM> previous_direction = this->mElementWeightedDirections[ (*it) ];
216
217 this->GetElement(*it)->CalculateWeightedDirection(this->mElementWeightedDirections[ (*it) ],
218 this->mElementJacobianDeterminants[ (*it) ]);
219
220 if (inner_prod(previous_direction, this->mElementWeightedDirections[ (*it) ]) < 0)
221 {
222 EXCEPTION("Moving node caused an subspace element to change direction");
223 }
224
225 }
227 }
228}
229
230template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
232{
233 if (this->mNodes[index]->IsDeleted())
235 EXCEPTION("Trying to delete a deleted node");
236 }
237 unsigned target_index = (unsigned)(-1);
238 bool found_target=false;
239 for (typename Node<SPACE_DIM>::ContainingElementIterator it = this->mNodes[index]->ContainingElementsBegin();
240 !found_target && it != this->mNodes[index]->ContainingElementsEnd();
241 ++it)
242 {
243 Element <ELEMENT_DIM,SPACE_DIM>* p_element = this->GetElement(*it);
244 for (unsigned i=0; i<=ELEMENT_DIM && !found_target; i++)
245 {
246 target_index = p_element->GetNodeGlobalIndex(i);
247 try
248 {
249 MoveMergeNode(index, target_index, false);
250 found_target = true;
251 }
252 catch (Exception&)
253 {
254 // Just try the next node
256 }
257 }
258 if (!found_target)
259 {
260 EXCEPTION("Failure to delete node");
261 }
262
263 MoveMergeNode(index, target_index);
264}
265
266template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
268{
269 assert(!this->mElements[index]->IsDeleted());
270 this->mElements[index]->MarkAsDeleted();
271 mDeletedElementIndices.push_back(index);
272
273 // Delete any nodes that are no longer attached to mesh
274 for (unsigned node_index = 0; node_index < this->mElements[index]->GetNumNodes(); ++node_index)
275 {
276 if (this->mElements[index]->GetNode(node_index)->GetNumContainingElements() == 0u)
277 {
278 if (this->mElements[index]->GetNode(node_index)->GetNumBoundaryElements() == 0u)
279 {
280 this->mElements[index]->GetNode(node_index)->MarkAsDeleted();
281 mDeletedNodeIndices.push_back(this->mElements[index]->GetNode(node_index)->GetIndex());
282 }
283 else if (this->mElements[index]->GetNode(node_index)->GetNumBoundaryElements() == 1u && ELEMENT_DIM == 1)
284 {
285 std::set<unsigned> indices = this->mElements[index]->GetNode(node_index)->rGetContainingBoundaryElementIndices();
286 assert(indices.size() == 1u);
287 this->mBoundaryElements[*indices.begin()]->MarkAsDeleted();
288 mDeletedBoundaryElementIndices.push_back(*indices.begin());
290 this->mElements[index]->GetNode(node_index)->MarkAsDeleted();
291 mDeletedNodeIndices.push_back(this->mElements[index]->GetNode(node_index)->GetIndex());
292 }
293 }
294 }
295}
296
297template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
299{
300 this->mNodes[index]->MarkAsDeleted();
301 mDeletedNodeIndices.push_back(index);
302}
304template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
306 unsigned targetIndex,
307 bool concreteMove)
308{
309 if (this->mNodes[index]->IsDeleted())
310 {
311 EXCEPTION("Trying to move a deleted node");
312 }
313
314 if (index == targetIndex)
315 {
316 EXCEPTION("Trying to merge a node with itself");
317 }
318 if (this->mNodes[index]->IsBoundaryNode())
320 if (!this->mNodes[targetIndex]->IsBoundaryNode())
321 {
322 EXCEPTION("A boundary node can only be moved on to another boundary node");
323 }
324 }
325 std::set<unsigned> unshared_element_indices;
326 std::set_difference(this->mNodes[index]->rGetContainingElementIndices().begin(),
327 this->mNodes[index]->rGetContainingElementIndices().end(),
328 this->mNodes[targetIndex]->rGetContainingElementIndices().begin(),
329 this->mNodes[targetIndex]->rGetContainingElementIndices().end(),
330 std::inserter(unshared_element_indices, unshared_element_indices.begin()));
331
332
333 if (unshared_element_indices.size() == this->mNodes[index]->rGetContainingElementIndices().size())
335 EXCEPTION("These nodes cannot be merged since they are not neighbours");
336 }
337
338 std::set<unsigned> unshared_boundary_element_indices;
339 std::set_difference(this->mNodes[index]->rGetContainingBoundaryElementIndices().begin(),
340 this->mNodes[index]->rGetContainingBoundaryElementIndices().end(),
341 this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().begin(),
342 this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().end(),
343 std::inserter(unshared_boundary_element_indices, unshared_boundary_element_indices.begin()));
344
345
346 if (this->mNodes[index]->IsBoundaryNode())
347 {
348 if (unshared_boundary_element_indices.size()
349 == this->mNodes[index]->rGetContainingBoundaryElementIndices().size())
350 {
351 //May be redundant (only thrown in 1D tests)
352 EXCEPTION("These nodes cannot be merged since they are not neighbours on the boundary");
354 }
355
356 this->mNodes[index]->rGetModifiableLocation() = this->mNodes[targetIndex]->rGetLocation();
357
358 for (std::set<unsigned>::const_iterator element_iter=unshared_element_indices.begin();
359 element_iter != unshared_element_indices.end();
360 element_iter++)
361 {
362 try
363 {
364 if (SPACE_DIM == ELEMENT_DIM)
365 {
366 this->GetElement(*element_iter)->CalculateInverseJacobian(this->mElementJacobians[(*element_iter)],
367 this->mElementJacobianDeterminants[(*element_iter)],
368 this->mElementInverseJacobians[ (*element_iter) ]);
369 }
370 else
371 {
372 this->GetElement(*element_iter)->CalculateWeightedDirection(this->mElementWeightedDirections[(*element_iter)],
373 this->mElementJacobianDeterminants[(*element_iter)]);
374 }
375
376 if (concreteMove)
377 {
378 this->GetElement(*element_iter)->ReplaceNode(this->mNodes[index], this->mNodes[targetIndex]);
379 }
380
381 }
382 catch (Exception&)
383 {
384 EXCEPTION("Moving node caused an element to have a non-positive Jacobian determinant");
385 }
386 }
387
388 for (std::set<unsigned>::const_iterator boundary_element_iter=
389 unshared_boundary_element_indices.begin();
390 boundary_element_iter != unshared_boundary_element_indices.end();
391 boundary_element_iter++)
392 {
393
394 this->GetBoundaryElement(*boundary_element_iter)->CalculateWeightedDirection(this->mBoundaryElementWeightedDirections[(*boundary_element_iter)],
395 this->mBoundaryElementJacobianDeterminants[(*boundary_element_iter)]);
396
397 if (concreteMove)
398 {
399 this->GetBoundaryElement(*boundary_element_iter)->ReplaceNode(this->mNodes[index], this->mNodes[targetIndex]);
400 }
401 }
402
403 std::set<unsigned> shared_element_indices;
404 std::set_intersection(this->mNodes[index]->rGetContainingElementIndices().begin(),
405 this->mNodes[index]->rGetContainingElementIndices().end(),
406 this->mNodes[targetIndex]->rGetContainingElementIndices().begin(),
407 this->mNodes[targetIndex]->rGetContainingElementIndices().end(),
408 std::inserter(shared_element_indices, shared_element_indices.begin()));
409
410 for (std::set<unsigned>::const_iterator element_iter=shared_element_indices.begin();
411 element_iter != shared_element_indices.end();
412 element_iter++)
413 {
414 if (concreteMove)
415 {
416 this->GetElement(*element_iter)->MarkAsDeleted();
417 this->mElementJacobianDeterminants[ (*element_iter) ] = 0.0; //This used to be done in MarkAsDeleted
418 mDeletedElementIndices.push_back(*element_iter);
419 }
420 else
421 {
422 this->mElementJacobianDeterminants[ (*element_iter) ] = 0.0;
423 }
424 }
425
426
427 std::set<unsigned> shared_boundary_element_indices;
428 std::set_intersection(this->mNodes[index]->rGetContainingBoundaryElementIndices().begin(),
429 this->mNodes[index]->rGetContainingBoundaryElementIndices().end(),
430 this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().begin(),
431 this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().end(),
432 std::inserter(shared_boundary_element_indices, shared_boundary_element_indices.begin()));
433
434 for (std::set<unsigned>::const_iterator boundary_element_iter=shared_boundary_element_indices.begin();
435 boundary_element_iter != shared_boundary_element_indices.end();
436 boundary_element_iter++)
437 {
438 if (concreteMove)
439 {
440 this->GetBoundaryElement(*boundary_element_iter)->MarkAsDeleted();
441 this->mBoundaryElementJacobianDeterminants[ (*boundary_element_iter) ] = 0.0; //This used to be done in MarkAsDeleted
442 mDeletedBoundaryElementIndices.push_back(*boundary_element_iter);
443 }
444 else
445 {
446 this->mBoundaryElementJacobianDeterminants[ (*boundary_element_iter) ] = 0.0;
447 this->mBoundaryElementWeightedDirections[ (*boundary_element_iter) ] = zero_vector<double>(SPACE_DIM);
448 }
449 }
450
451 if (concreteMove)
452 {
453 this->mNodes[index]->MarkAsDeleted();
454 mDeletedNodeIndices.push_back(index);
455 }
456}
457
458template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
462{
463 //Check that the point is in the element
464 if (pElement->IncludesPoint(point, true) == false)
465 {
466 EXCEPTION("RefineElement could not be started (point is not in element)");
467 }
468
469 // Add a new node from the point that is passed to RefineElement
470 unsigned new_node_index = AddNode(new Node<SPACE_DIM>(0, point.rGetLocation()));
471 // Note: the first argument is the index of the node, which is going to be
472 // overridden by AddNode, so it can safely be ignored
473
474 // This loop constructs the extra elements which are going to fill the space
475 for (unsigned i = 0; i < ELEMENT_DIM; i++)
476 {
477
478 // First, make a copy of the current element making sure we update its index
479 unsigned new_elt_index;
480 if (mDeletedElementIndices.empty())
481 {
482 new_elt_index = this->mElements.size();
483 }
484 else
485 {
486 new_elt_index = mDeletedElementIndices.back();
487 mDeletedElementIndices.pop_back();
488 }
489
490 Element<ELEMENT_DIM,SPACE_DIM>* p_new_element=
491 new Element<ELEMENT_DIM,SPACE_DIM>(*pElement, new_elt_index);
492
493 // Second, update the node in the element with the new one
494 p_new_element->UpdateNode(ELEMENT_DIM-1-i, this->mNodes[new_node_index]);
495
496 // Third, add the new element to the set
497 if ((unsigned) new_elt_index == this->mElements.size())
498 {
499 this->mElements.push_back(p_new_element);
500 }
501 else
502 {
503 delete this->mElements[new_elt_index];
504 this->mElements[new_elt_index] = p_new_element;
505 }
506 }
507
508 // Lastly, update the last node in the element to be refined
509 pElement->UpdateNode(ELEMENT_DIM, this->mNodes[new_node_index]);
510
511 return new_node_index;
512}
513
514template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
516{
517 if (!this->mNodes[index]->IsBoundaryNode() )
518 {
519 EXCEPTION(" You may only delete a boundary node ");
520 }
521
522 this->mNodes[index]->MarkAsDeleted();
523 mDeletedNodeIndices.push_back(index);
524
525 // Update the boundary node vector
526 typename std::vector<Node<SPACE_DIM>*>::iterator b_node_iter
527 = std::find(this->mBoundaryNodes.begin(), this->mBoundaryNodes.end(), this->mNodes[index]);
528 this->mBoundaryNodes.erase(b_node_iter);
529
530 // Remove boundary elements containing this node
531 std::set<unsigned> boundary_element_indices = this->mNodes[index]->rGetContainingBoundaryElementIndices();
532 std::set<unsigned>::const_iterator boundary_element_indices_iterator = boundary_element_indices.begin();
533 while (boundary_element_indices_iterator != boundary_element_indices.end())
534 {
535 BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_boundary_element = this->GetBoundaryElement(*boundary_element_indices_iterator);
536 p_boundary_element->MarkAsDeleted();
537 mDeletedBoundaryElementIndices.push_back(*boundary_element_indices_iterator);
538 boundary_element_indices_iterator++;
539 }
540
541 // Remove elements containing this node
542 std::set<unsigned> element_indices = this->mNodes[index]->rGetContainingElementIndices();
543 std::set<unsigned>::const_iterator element_indices_iterator = element_indices.begin();
544 while (element_indices_iterator != element_indices.end())
545 {
546 Element<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(*element_indices_iterator);
547 for (unsigned i=0; i<p_element->GetNumNodes(); i++)
548 {
549 Node<SPACE_DIM>* p_node = p_element->GetNode(i);
550 if (!p_node->IsDeleted())
551 {
552 p_node->SetAsBoundaryNode();
553 // Update the boundary node vector
554 this->mBoundaryNodes.push_back(p_node);
555 }
556 }
557 p_element->MarkAsDeleted();
558 mDeletedElementIndices.push_back(p_element->GetIndex());
559 element_indices_iterator++;
560 }
561}
562
563template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
565{
566 assert(!mAddedNodes);
567 map.Resize(this->GetNumAllNodes());
568
569 std::vector<Element<ELEMENT_DIM, SPACE_DIM> *> live_elements;
570
571 for (unsigned i=0; i<this->mElements.size(); i++)
572 {
573 assert(i==this->mElements[i]->GetIndex()); // We need this to be true to be able to reindex the Jacobian cache
574 if (this->mElements[i]->IsDeleted())
575 {
576 delete this->mElements[i];
577 }
578 else
579 {
580 live_elements.push_back(this->mElements[i]);
581
582 unsigned this_element_index = this->mElements[i]->GetIndex();
583 if (SPACE_DIM == ELEMENT_DIM)
584 {
585 this->mElementJacobians[live_elements.size()-1] = this->mElementJacobians[this_element_index];
586 this->mElementInverseJacobians[live_elements.size()-1] = this->mElementInverseJacobians[this_element_index];
587 }
588 else
589 {
590 this->mElementWeightedDirections[live_elements.size()-1] = this->mElementWeightedDirections[this_element_index];
591 }
592 this->mElementJacobianDeterminants[live_elements.size()-1] = this->mElementJacobianDeterminants[this_element_index];
593 }
594 }
595
596 assert(mDeletedElementIndices.size() == this->mElements.size()-live_elements.size());
597 mDeletedElementIndices.clear();
598 this->mElements = live_elements;
599 unsigned num_elements = this->mElements.size();
600
601 if (SPACE_DIM == ELEMENT_DIM)
602 {
603 this->mElementJacobians.resize(num_elements);
604 this->mElementInverseJacobians.resize(num_elements);
605 }
606 else
607 {
608 this->mElementWeightedDirections.resize(num_elements);
609 }
610 this->mElementJacobianDeterminants.resize(num_elements);
611
612 std::vector<Node<SPACE_DIM> *> live_nodes;
613 for (unsigned i=0; i<this->mNodes.size(); i++)
614 {
615 if (this->mNodes[i]->IsDeleted())
616 {
617 delete this->mNodes[i];
618 map.SetDeleted(i);
619 }
620 else
621 {
622 live_nodes.push_back(this->mNodes[i]);
623 // the nodes will have their index set to be the index into the live_nodes
624 // vector further down
625 map.SetNewIndex(i, (unsigned)(live_nodes.size()-1));
626 }
627 }
628
629 assert(mDeletedNodeIndices.size() == this->mNodes.size()-live_nodes.size());
630 this->mNodes = live_nodes;
631 mDeletedNodeIndices.clear();
632
633 std::vector<BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *> live_boundary_elements;
634 for (unsigned i=0; i<this->mBoundaryElements.size(); i++)
635 {
636 if (this->mBoundaryElements[i]->IsDeleted())
637 {
638 delete this->mBoundaryElements[i];
639 }
640 else
641 {
642 live_boundary_elements.push_back(this->mBoundaryElements[i]);
643
644 this->mBoundaryElementWeightedDirections[live_boundary_elements.size()-1] = this->mBoundaryElementWeightedDirections[this->mBoundaryElements[i]->GetIndex()];
645 this->mBoundaryElementJacobianDeterminants[live_boundary_elements.size()-1] = this->mBoundaryElementJacobianDeterminants[this->mBoundaryElements[i]->GetIndex()];
646 }
647 }
648
649 assert(mDeletedBoundaryElementIndices.size() == this->mBoundaryElements.size()-live_boundary_elements.size());
650 this->mBoundaryElements = live_boundary_elements;
651 mDeletedBoundaryElementIndices.clear();
652
653 unsigned num_boundary_elements = this->mBoundaryElements.size();
654
655 this->mBoundaryElementWeightedDirections.resize(num_boundary_elements);
656 this->mBoundaryElementJacobianDeterminants.resize(num_boundary_elements);
657
658 for (unsigned i=0; i<this->mNodes.size(); i++)
659 {
660 this->mNodes[i]->SetIndex(i);
661 }
662
663 for (unsigned i=0; i<this->mElements.size(); i++)
664 {
665 this->mElements[i]->ResetIndex(i);
666 }
667
668 for (unsigned i=0; i<this->mBoundaryElements.size(); i++)
669 {
670 this->mBoundaryElements[i]->ResetIndex(i);
671 }
672}
673
674template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
676{
677 // Make sure that we are in the correct dimension - this code will be eliminated at compile time
678 assert( ELEMENT_DIM == SPACE_DIM ); // LCOV_EXCL_LINE
679
680 // Avoid some triangle/tetgen errors: need at least four
681 // nodes for tetgen, and at least three for triangle
682 if (GetNumNodes() <= SPACE_DIM)
683 {
684 EXCEPTION("The number of nodes must exceed the spatial dimension.");
685 }
686
687 // Make sure the map is big enough
688 map.Resize(this->GetNumAllNodes());
689 if (mAddedNodes || !mDeletedNodeIndices.empty())
690 {
691 // Size of mesh is about to change
692 if (this->mpDistributedVectorFactory)
693 {
694 delete this->mpDistributedVectorFactory;
695 this->mpDistributedVectorFactory = new DistributedVectorFactory(this->GetNumNodes());
696 }
697 }
698 if (SPACE_DIM == 1)
699 {
700 // Store the node locations
701 std::vector<c_vector<double, SPACE_DIM> > old_node_locations;
702 unsigned new_index = 0;
703 for (unsigned i=0; i<this->GetNumAllNodes(); i++)
704 {
705 if (this->mNodes[i]->IsDeleted())
706 {
707 map.SetDeleted(i);
708 }
709 else
710 {
711 map.SetNewIndex(i, new_index);
712 old_node_locations.push_back(this->mNodes[i]->rGetLocation());
713 new_index++;
714 }
715 }
716
717 // Remove current data
718 Clear();
719
720 // Construct the nodes and boundary nodes
721 for (unsigned node_index=0; node_index<old_node_locations.size(); node_index++)
722 {
723 // As we're in 1D, the boundary nodes are simply at either end of the mesh
724 bool is_boundary_node = (node_index==0 || node_index==old_node_locations.size()-1);
725
726 Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index, old_node_locations[node_index], is_boundary_node);
727 this->mNodes.push_back(p_node);
728
729 if (is_boundary_node)
730 {
731 this->mBoundaryNodes.push_back(p_node);
732 }
733 }
734
735 // Create a map between node indices and node locations
736 std::map<double, unsigned> location_index_map;
737 for (unsigned i=0; i<this->mNodes.size(); i++)
738 {
739 location_index_map[this->mNodes[i]->rGetLocation()[0]] = this->mNodes[i]->GetIndex();
740 }
741
742 // Use this map to generate a vector of node indices that are ordered spatially
743 std::vector<unsigned> node_indices_ordered_spatially;
744 for (std::map<double, unsigned>::iterator iter = location_index_map.begin();
745 iter != location_index_map.end();
746 ++iter)
747 {
748 node_indices_ordered_spatially.push_back(iter->second);
749 }
750
751 // Construct the elements
752 this->mElements.reserve(old_node_locations.size()-1);
753 for (unsigned element_index=0; element_index<old_node_locations.size()-1; element_index++)
754 {
755 std::vector<Node<SPACE_DIM>*> nodes;
756 for (unsigned j=0; j<2; j++)
757 {
758 unsigned global_node_index = node_indices_ordered_spatially[element_index + j];
759 assert(global_node_index < this->mNodes.size());
760 nodes.push_back(this->mNodes[global_node_index]);
761 }
762 this->mElements.push_back(new Element<ELEMENT_DIM, SPACE_DIM>(element_index, nodes));
763 }
764
765 // Construct the two boundary elements - as we're in 1D, these are simply at either end of the mesh
766 std::vector<Node<SPACE_DIM>*> nodes;
767 nodes.push_back(this->mNodes[0]);
768 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(0, nodes));
769
770 nodes.clear();
771 nodes.push_back(this->mNodes[old_node_locations.size()-1]);
772 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(1, nodes));
773
774 this->RefreshJacobianCachedData();
775 }
776 else if (SPACE_DIM==2) // In 2D, remesh using triangle via library calls
777 {
778 struct triangulateio mesher_input, mesher_output;
779 this->InitialiseTriangulateIo(mesher_input);
780 this->InitialiseTriangulateIo(mesher_output);
781
782 this->ExportToMesher(map, mesher_input);
783
784 // Library call
785 triangulate((char*)"Qze", &mesher_input, &mesher_output, nullptr);
786
787 this->ImportFromMesher(mesher_output, mesher_output.numberoftriangles, mesher_output.trianglelist, mesher_output.numberofedges, mesher_output.edgelist, mesher_output.edgemarkerlist);
788
789 //Tidy up triangle
790 this->FreeTriangulateIo(mesher_input);
791 this->FreeTriangulateIo(mesher_output);
792 }
793 else // in 3D, remesh using tetgen
794 {
795
796 class tetgen::tetgenio mesher_input, mesher_output;
797
798 this->ExportToMesher(map, mesher_input);
799
800 // Library call
801 tetgen::tetrahedralize((char*)"Qz", &mesher_input, &mesher_output);
802
803 this->ImportFromMesher(mesher_output, mesher_output.numberoftetrahedra, mesher_output.tetrahedronlist, mesher_output.numberoftrifaces, mesher_output.trifacelist, nullptr);
804 }
805}
806
807template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
809{
810 NodeMap map(GetNumNodes());
811 ReMesh(map);
812}
813
814template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
815std::vector<c_vector<unsigned, 5> > MutableMesh<ELEMENT_DIM, SPACE_DIM>::SplitLongEdges(double cutoffLength)
816{
817 assert(ELEMENT_DIM == 2); // LCOV_EXCL_LINE
818 assert(SPACE_DIM == 3); // LCOV_EXCL_LINE
819
820 std::vector<c_vector<unsigned, 5> > history;
821
822 bool long_edge_exists = true;
823
824 while (long_edge_exists)
825 {
826 std::set<std::pair<unsigned, unsigned> > long_edges;
827
828 // Loop over elements to check for Long edges
829 for (typename AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::ElementIterator elem_iter = this->GetElementIteratorBegin();
830 elem_iter != this->GetElementIteratorEnd();
831 ++elem_iter)
832 {
833 unsigned num_nodes = ELEMENT_DIM+1;
834
835 // Loop over element vertices
836 for (unsigned local_index=0; local_index<num_nodes; local_index++)
837 {
838 // Find locations of current node (node a) and anticlockwise node (node b)
839 Node<SPACE_DIM>* p_node_a = elem_iter->GetNode(local_index);
840 unsigned local_index_plus_one = (local_index+1)%num_nodes;
841 Node<SPACE_DIM>* p_node_b = elem_iter->GetNode(local_index_plus_one);
842
843 // Find distance between nodes
844 double distance_between_nodes = this->GetDistanceBetweenNodes(p_node_a->GetIndex(), p_node_b->GetIndex());
845
846 if (distance_between_nodes > cutoffLength)
847 {
848 if (p_node_a->GetIndex() < p_node_b->GetIndex())
849 {
850 std::pair<unsigned, unsigned> long_edge(p_node_a->GetIndex(),p_node_b->GetIndex());
851 long_edges.insert(long_edge);
852 }
853 else
854 {
855 std::pair<unsigned, unsigned> long_edge(p_node_b->GetIndex(),p_node_a->GetIndex());
856 long_edges.insert(long_edge);
857 }
858 }
859 }
860 }
861
862 if (long_edges.size() > 0) //Split the edges in decreasing order.
863 {
864 while (long_edges.size() > 0)
865 {
866 double longest_edge = 0.0;
867 std::set<std::pair<unsigned, unsigned> >::iterator longest_edge_iter;
868
869 //Find the longest edge in the set and split it
870 for (std::set<std::pair<unsigned, unsigned> >::iterator edge_iter = long_edges.begin();
871 edge_iter != long_edges.end();
872 ++edge_iter)
873 {
874 unsigned node_a_global_index = edge_iter->first;
875 unsigned node_b_global_index = edge_iter->second;
876
877 double distance_between_nodes = this->GetDistanceBetweenNodes(node_a_global_index, node_b_global_index);
878
879 if (distance_between_nodes > longest_edge)
880 {
881 longest_edge = distance_between_nodes;
882 longest_edge_iter = edge_iter;
883 }
884 }
885 assert(longest_edge >0);
886
887 c_vector<unsigned, 3> new_node_index = SplitEdge(this->GetNode(longest_edge_iter->first), this->GetNode(longest_edge_iter->second));
888
889 c_vector<unsigned, 5> node_set;
890 node_set(0) = new_node_index[0];
891 node_set(1) = longest_edge_iter->first;
892 node_set(2) = longest_edge_iter->second;
893 node_set(3) = new_node_index[1];
894 node_set(4) = new_node_index[2];
895 history.push_back(node_set);
896
897 // Delete pair from set
898 long_edges.erase(*longest_edge_iter);
899 }
900 }
901 else
902 {
903 long_edge_exists = false;
904 }
905 }
906
907 return history;
908}
909
910template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
912{
913 c_vector<unsigned, 3> new_node_index_vector;
914
915 std::set<unsigned> elements_of_node_a = pNodeA->rGetContainingElementIndices();
916 std::set<unsigned> elements_of_node_b = pNodeB->rGetContainingElementIndices();
917
918 std::set<unsigned> intersection_elements;
919 std::set_intersection(elements_of_node_a.begin(), elements_of_node_a.end(),
920 elements_of_node_b.begin(), elements_of_node_b.end(),
921 std::inserter(intersection_elements, intersection_elements.begin()));
922
923 // Create the new node
924 c_vector<double, SPACE_DIM> new_node_location = pNodeA->rGetLocation() + 0.5*this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation());
925
926 bool is_boundary_node = intersection_elements.size() == 1; // If only in one element then its on the boundary
927
928 Node<SPACE_DIM>* p_new_node = new Node<SPACE_DIM>(0u,new_node_location,is_boundary_node); // Index is rewritten once added to mesh
929
930 unsigned new_node_index = this->AddNode(p_new_node);
931
932 new_node_index_vector[0] = new_node_index;
933
934 unsigned counter = 1;
935
936 for (std::set<unsigned>::const_iterator it = intersection_elements.begin(); it != intersection_elements.end(); ++it)
937 {
938 unsigned elementIndex = *it;
939
940 Element<ELEMENT_DIM,SPACE_DIM>* p_original_element = this->GetElement(elementIndex);
941
942 // First, make a copy of the current element and assign an unused index
943 Element<ELEMENT_DIM,SPACE_DIM>* p_new_element = new Element<ELEMENT_DIM,SPACE_DIM>(*p_original_element, UINT_MAX);
944
945 // Second, add the new element to the set of existing elements. This method will assign a proper index to the element.
946 AddElement(p_new_element);
947
948 // Third, update node a in the element with the new one
949 p_new_element->ReplaceNode(pNodeA, this->mNodes[new_node_index]);
950
951 // Last, update node b in the original element with the new one
952 p_original_element->ReplaceNode(pNodeB, this->mNodes[new_node_index]);
953
954 //Add node in both of these elements to new_node_index_vector (this enables us to add a new spring in the MeshBasedCellPopulation
955 unsigned other_node_index = UNSIGNED_UNSET;
956
957 if ((p_original_element->GetNodeGlobalIndex(0) != new_node_index) &&
958 (p_original_element->GetNodeGlobalIndex(0) != pNodeA->GetIndex()))
959 {
960 other_node_index = p_original_element->GetNodeGlobalIndex(0);
961 }
962 else if ((p_original_element->GetNodeGlobalIndex(1) != new_node_index) &&
963 (p_original_element->GetNodeGlobalIndex(1) != pNodeA->GetIndex()))
964 {
965 other_node_index = p_original_element->GetNodeGlobalIndex(1);
966 }
967 else if ((p_original_element->GetNodeGlobalIndex(2) != new_node_index) &&
968 (p_original_element->GetNodeGlobalIndex(2) != pNodeA->GetIndex()))
969 {
970 other_node_index = p_original_element->GetNodeGlobalIndex(2);
971 }
972 else
973 {
975 }
976 new_node_index_vector[counter] = other_node_index;
977 counter++;
978 }
979
980 assert(counter < 4);
981 assert(counter > 1);// need to be in at least one element
982
983 if (counter == 2) // only one new element
984 {
985 new_node_index_vector[2] = UNSIGNED_UNSET;
986 }
987
988 return new_node_index_vector;
989}
990
991template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
993{
994 assert(ELEMENT_DIM == SPACE_DIM); // LCOV_EXCL_LINE
995 unsigned num_nodes = pElement->GetNumNodes();
996 std::set<unsigned> neighbouring_elements_indices;
997 std::set< Element<ELEMENT_DIM,SPACE_DIM> *> neighbouring_elements;
998 std::set<unsigned> neighbouring_nodes_indices;
999
1000 // Form a set of neighbouring elements via the nodes
1001 for (unsigned i=0; i<num_nodes; i++)
1002 {
1003 Node<SPACE_DIM>* p_node = pElement->GetNode(i);
1004 neighbouring_elements_indices = p_node->rGetContainingElementIndices();
1005 for (std::set<unsigned>::const_iterator it = neighbouring_elements_indices.begin();
1006 it != neighbouring_elements_indices.end();
1007 ++it)
1008 {
1009 neighbouring_elements.insert(this->GetElement(*it));
1010 }
1011 }
1012 neighbouring_elements.erase(pElement);
1013
1014 // For each neighbouring element find the supporting nodes
1015 typedef typename std::set<Element<ELEMENT_DIM,SPACE_DIM> *>::const_iterator ElementIterator;
1016
1017 for (ElementIterator it = neighbouring_elements.begin();
1018 it != neighbouring_elements.end();
1019 ++it)
1020 {
1021 for (unsigned i=0; i<num_nodes; i++)
1022 {
1023 neighbouring_nodes_indices.insert((*it)->GetNodeGlobalIndex(i));
1024 }
1025 }
1026
1027 // Remove the nodes that support this element
1028 for (unsigned i = 0; i < num_nodes; i++)
1029 {
1030 neighbouring_nodes_indices.erase(pElement->GetNodeGlobalIndex(i));
1031 }
1032
1033 // Get the circumsphere information
1034 c_vector<double, SPACE_DIM+1> this_circum_centre = zero_vector<double>(SPACE_DIM+1);
1035
1036 this_circum_centre = pElement->CalculateCircumsphere(this->mElementJacobians[pElement->GetIndex()], this->mElementInverseJacobians[pElement->GetIndex()]);
1037
1038 // Copy the actually circumcentre into a smaller vector
1039 c_vector<double, ELEMENT_DIM> circum_centre = zero_vector<double>(ELEMENT_DIM);
1040 for (unsigned i=0; i<ELEMENT_DIM; i++)
1041 {
1042 circum_centre[i] = this_circum_centre[i];
1043 }
1044
1045 for (std::set<unsigned>::const_iterator it = neighbouring_nodes_indices.begin();
1046 it != neighbouring_nodes_indices.end();
1047 ++it)
1048 {
1049 c_vector<double, ELEMENT_DIM> node_location = this->GetNode(*it)->rGetLocation();
1050
1051 // Calculate vector from circumcenter to node
1052 node_location -= circum_centre;
1053
1054 // This is to calculate the squared distance between them
1055 double squared_distance = inner_prod(node_location, node_location);
1056
1057 // If the squared distance is less than the elements circum-radius(squared),
1058 // then the Voronoi property is violated.
1059 if (squared_distance < this_circum_centre[ELEMENT_DIM])
1060 {
1061 // We know the node is inside the circumsphere, but we don't know how far
1062 double radius = sqrt(this_circum_centre[ELEMENT_DIM]);
1063 double distance = radius - sqrt(squared_distance);
1064
1065 // If the node penetration is greater than supplied maximum penetration factor
1066 if (distance/radius > maxPenetration)
1067 {
1068 return false;
1069 }
1070 }
1071 }
1072 return true;
1073}
1074
1075template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1077{
1078 // Looping through all the elements in the mesh
1080 for (unsigned i=0; i<this->mElements.size(); i++)
1081 {
1082 // Check if the element is not deleted
1083 if (!this->mElements[i]->IsDeleted())
1084 {
1085 // Checking the Voronoi of the Element
1086 if (CheckIsVoronoi(this->mElements[i], maxPenetration) == false)
1087 {
1088 return false;
1089 }
1090 }
1091 }
1092 return true;
1093}
1094
1095// Explicit instantiation
1096template class MutableMesh<1,1>;
1097template class MutableMesh<1,2>;
1098template class MutableMesh<1,3>;
1099template class MutableMesh<2,2>;
1100template class MutableMesh<2,3>;
1101template class MutableMesh<3,3>;
1102
1103// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define NEVER_REACHED
const unsigned UNSIGNED_UNSET
Definition Exception.hpp:53
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
void ReplaceNode(Node< SPACE_DIM > *pOldNode, Node< SPACE_DIM > *pNewNode)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
unsigned GetNumNodes() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
unsigned GetIndex() const
bool mMeshChangesDuringSimulation
c_vector< double, DIM > & rGetLocation()
bool IncludesPoint(const ChastePoint< SPACE_DIM > &rTestPoint, bool strict=false)
Definition Element.cpp:310
void UpdateNode(const unsigned &rIndex, Node< SPACE_DIM > *pNode)
Definition Element.cpp:85
void ResetIndex(unsigned index)
Definition Element.cpp:100
void MarkAsDeleted()
Definition Element.cpp:74
c_vector< double, SPACE_DIM+1 > CalculateCircumsphere(c_matrix< double, SPACE_DIM, ELEMENT_DIM > &rJacobian, c_matrix< double, ELEMENT_DIM, SPACE_DIM > &rInverseJacobian)
Definition Element.cpp:114
unsigned GetNumElements() const
bool CheckIsVoronoi(Element< ELEMENT_DIM, SPACE_DIM > *pElement, double maxPenetration)
c_vector< unsigned, 3 > SplitEdge(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
void ReIndex(NodeMap &map)
unsigned AddElement(Element< ELEMENT_DIM, SPACE_DIM > *pNewElement)
unsigned GetNumNodes() const
void DeleteNodePriorToReMesh(unsigned index)
virtual void SetNode(unsigned index, ChastePoint< SPACE_DIM > point, bool concreteMove=true)
virtual void DeleteNode(unsigned index)
void RescaleMeshFromBoundaryNode(ChastePoint< 1 > updatedPoint, unsigned boundaryNodeIndex)
void DeleteBoundaryNodeAt(unsigned index)
std::vector< c_vector< unsigned, 5 > > SplitLongEdges(double cutoffLength)
virtual ~MutableMesh()
unsigned RefineElement(Element< ELEMENT_DIM, SPACE_DIM > *pElement, ChastePoint< SPACE_DIM > point)
unsigned GetNumBoundaryElements() const
virtual unsigned AddNode(Node< SPACE_DIM > *pNewNode)
void MoveMergeNode(unsigned index, unsigned targetIndex, bool concreteMove=true)
virtual void DeleteElement(unsigned index)
void SetDeleted(unsigned index)
Definition NodeMap.cpp:76
void SetNewIndex(unsigned oldIndex, unsigned newIndex)
Definition NodeMap.cpp:66
void Resize(unsigned size)
Definition NodeMap.cpp:52
Definition Node.hpp:59
std::set< unsigned > & rGetContainingElementIndices()
Definition Node.cpp:300
void SetIndex(unsigned index)
Definition Node.cpp:121
bool IsDeleted() const
Definition Node.cpp:412
const c_vector< double, SPACE_DIM > & rGetLocation() const
Definition Node.cpp:139
unsigned GetIndex() const
Definition Node.cpp:158
void SetAsBoundaryNode(bool value=true)
Definition Node.cpp:127
virtual void Clear()