00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "MutableVertexMesh.hpp"
00030 #include "RandomNumberGenerator.hpp"
00031 #include "UblasCustomFunctions.hpp"
00032 #include "Warnings.hpp"
00033 #include "LogFile.hpp"
00034
00035 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00036 MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::MutableVertexMesh(std::vector<Node<SPACE_DIM>*> nodes,
00037 std::vector<VertexElement<ELEMENT_DIM,SPACE_DIM>*> vertexElements,
00038 double cellRearrangementThreshold,
00039 double t2Threshold,
00040 double cellRearrangementRatio)
00041 : mCellRearrangementThreshold(cellRearrangementThreshold),
00042 mCellRearrangementRatio(cellRearrangementRatio),
00043 mT2Threshold(t2Threshold)
00044 {
00045 assert(cellRearrangementThreshold > 0.0);
00046 assert(t2Threshold > 0.0);
00047
00048
00049 Clear();
00050
00051
00052 for (unsigned node_index=0; node_index<nodes.size(); node_index++)
00053 {
00054 Node<SPACE_DIM>* p_temp_node = nodes[node_index];
00055 this->mNodes.push_back(p_temp_node);
00056 }
00057 for (unsigned elem_index=0; elem_index<vertexElements.size(); elem_index++)
00058 {
00059 VertexElement<ELEMENT_DIM,SPACE_DIM>* p_temp_vertex_element = vertexElements[elem_index];
00060 this->mElements.push_back(p_temp_vertex_element);
00061 }
00062
00063 if (SPACE_DIM == 3)
00064 {
00065
00066 std::set<unsigned> faces_counted;
00067
00068
00069 for (unsigned elem_index=0; elem_index<this->mElements.size(); elem_index++)
00070 {
00071
00072 for (unsigned face_index=0; face_index<this->mElements[elem_index]->GetNumFaces(); face_index++)
00073 {
00074 VertexElement<ELEMENT_DIM-1, SPACE_DIM>* p_face = this->mElements[elem_index]->GetFace(face_index);
00075
00076
00077 if (faces_counted.find(p_face->GetIndex()) == faces_counted.end())
00078 {
00079 this->mFaces.push_back(p_face);
00080 faces_counted.insert(p_face->GetIndex());
00081 }
00082 }
00083 }
00084 }
00085
00086
00087 for (unsigned index=0; index<this->mElements.size(); index++)
00088 {
00089 VertexElement<ELEMENT_DIM,SPACE_DIM>* p_temp_vertex_element = this->mElements[index];
00090 for (unsigned node_index=0; node_index<p_temp_vertex_element->GetNumNodes(); node_index++)
00091 {
00092 Node<SPACE_DIM>* p_temp_node = p_temp_vertex_element->GetNode(node_index);
00093 p_temp_node->AddElement(p_temp_vertex_element->GetIndex());
00094 }
00095 }
00096
00097 this->mMeshChangesDuringSimulation = true;
00098 }
00099
00100
00101 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00102 MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::MutableVertexMesh()
00103 : mCellRearrangementThreshold(0.01),
00104 mCellRearrangementRatio(1.5),
00105 mT2Threshold(0.001)
00106 {
00107 this->mMeshChangesDuringSimulation = true;
00108 Clear();
00109 }
00110
00111
00112 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00113 MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::~MutableVertexMesh()
00114 {
00115 Clear();
00116 }
00117
00118 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00119 double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetCellRearrangementThreshold() const
00120 {
00121 return mCellRearrangementThreshold;
00122 }
00123
00124
00125 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00126 double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetT2Threshold() const
00127 {
00128 return mT2Threshold;
00129 }
00130
00131
00132 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00133 double MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetCellRearrangementRatio() const
00134 {
00135 return mCellRearrangementRatio;
00136 }
00137
00138
00139 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00140 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetCellRearrangementThreshold(double cellRearrangementThreshold)
00141 {
00142 mCellRearrangementThreshold = cellRearrangementThreshold;
00143 }
00144
00145
00146 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00147 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetT2Threshold(double t2Threshold)
00148 {
00149 mT2Threshold = t2Threshold;
00150 }
00151
00152
00153 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00154 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetCellRearrangementRatio(double cellRearrangementRatio)
00155 {
00156 mCellRearrangementRatio = cellRearrangementRatio;
00157 }
00158
00159
00160 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00161 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::Clear()
00162 {
00163 mDeletedNodeIndices.clear();
00164 mDeletedElementIndices.clear();
00165
00166 VertexMesh<ELEMENT_DIM, SPACE_DIM>::Clear();
00167 }
00168
00169
00170 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00171 unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNumNodes() const
00172 {
00173 return this->mNodes.size() - mDeletedNodeIndices.size();
00174 }
00175
00176
00177 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00178 unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNumElements() const
00179 {
00180 return this->mElements.size() - mDeletedElementIndices.size();
00181 }
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00192 std::vector< c_vector<double, SPACE_DIM> > MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLocationsOfT1Swaps()
00193 {
00194 return mLocationsOfT1Swaps;
00195 }
00196
00197 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00198 std::vector< c_vector<double, SPACE_DIM> > MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLocationsOfT3Swaps()
00199 {
00200 return mLocationsOfT3Swaps;
00201 }
00202
00203 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00204 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ClearLocationsOfT1Swaps()
00205 {
00206 mLocationsOfT1Swaps.clear();
00207 }
00208
00209 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00210 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ClearLocationsOfT3Swaps()
00211 {
00212 mLocationsOfT3Swaps.clear();
00213 }
00214
00215
00216 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00217 unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::AddNode(Node<SPACE_DIM>* pNewNode)
00218 {
00219 if (mDeletedNodeIndices.empty())
00220 {
00221 pNewNode->SetIndex(this->mNodes.size());
00222 this->mNodes.push_back(pNewNode);
00223 }
00224 else
00225 {
00226 unsigned index = mDeletedNodeIndices.back();
00227 pNewNode->SetIndex(index);
00228 mDeletedNodeIndices.pop_back();
00229 delete this->mNodes[index];
00230 this->mNodes[index] = pNewNode;
00231 }
00232 return pNewNode->GetIndex();
00233 }
00234
00235
00236 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00237 unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::AddElement(VertexElement<ELEMENT_DIM,SPACE_DIM>* pNewElement)
00238 {
00239 unsigned new_element_index = pNewElement->GetIndex();
00240
00241 if (new_element_index == this->mElements.size())
00242 {
00243 this->mElements.push_back(pNewElement);
00244 }
00245 else
00246 {
00247 this->mElements[new_element_index] = pNewElement;
00248 }
00249 pNewElement->RegisterWithNodes();
00250 return pNewElement->GetIndex();
00251 }
00252
00253
00254 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00255 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::SetNode(unsigned nodeIndex, ChastePoint<SPACE_DIM> point)
00256 {
00257 this->mNodes[nodeIndex]->SetPoint(point);
00258 }
00259
00260
00261 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00262 unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DivideElementAlongGivenAxis(VertexElement<ELEMENT_DIM,SPACE_DIM>* pElement, c_vector<double, SPACE_DIM> axisOfDivision,
00263 bool placeOriginalElementBelow)
00264 {
00265
00266
00267 assert(SPACE_DIM == 2);
00268 assert(ELEMENT_DIM == SPACE_DIM);
00269
00270
00271 c_vector<double, SPACE_DIM> centroid = GetCentroidOfElement(pElement->GetIndex());
00272
00273
00274 c_vector<double, SPACE_DIM> perp_axis;
00275 perp_axis(0) = -axisOfDivision(1);
00276 perp_axis(1) = axisOfDivision(0);
00277
00278
00279
00280
00281
00282
00283 unsigned num_nodes = pElement->GetNumNodes();
00284 std::vector<unsigned> intersecting_nodes;
00285 for (unsigned i=0; i<num_nodes; i++)
00286 {
00287 bool is_current_node_on_left = (inner_prod(GetVectorFromAtoB(pElement->GetNodeLocation(i), centroid), perp_axis) >= 0);
00288 bool is_next_node_on_left = (inner_prod(GetVectorFromAtoB(pElement->GetNodeLocation((i+1)%num_nodes), centroid), perp_axis) >= 0);
00289
00290 if (is_current_node_on_left != is_next_node_on_left)
00291 {
00292 intersecting_nodes.push_back(i);
00293 }
00294 }
00295
00296
00297 if (intersecting_nodes.size() != 2)
00298 {
00299 EXCEPTION("Cannot proceed with element division: the given axis of division does not cross two edges of the element");
00300 }
00301
00302 std::vector<unsigned> division_node_global_indices;
00303 unsigned nodes_added = 0;
00304
00305
00306 for (unsigned i=0; i<intersecting_nodes.size(); i++)
00307 {
00308
00309
00310
00311
00312
00313
00314
00315 Node<SPACE_DIM>* p_node_A = pElement->GetNode((intersecting_nodes[i]+nodes_added)%pElement->GetNumNodes());
00316 Node<SPACE_DIM>* p_node_B = pElement->GetNode((intersecting_nodes[i]+nodes_added+1)%pElement->GetNumNodes());
00317
00318
00319 std::set<unsigned> elems_containing_node_A = p_node_A->rGetContainingElementIndices();
00320 std::set<unsigned> elems_containing_node_B = p_node_B->rGetContainingElementIndices();
00321
00322
00323 c_vector<double, SPACE_DIM> position_a = p_node_A->rGetLocation();
00324 c_vector<double, SPACE_DIM> position_b = p_node_B->rGetLocation();
00325 c_vector<double, SPACE_DIM> a_to_b = GetVectorFromAtoB(position_a, position_b);
00326
00327
00328 c_vector<double, SPACE_DIM> intersection;
00329
00330 if (norm_2(a_to_b)< 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
00331 {
00332 WARNING("Edge is too small for normal division, putting node in the middle of a and b, there may be T1Swaps straight away.");
00334 intersection = position_a + 0.5*a_to_b;
00335 }
00336 else
00337 {
00338
00339 double determinant = a_to_b[0]*axisOfDivision[1] - a_to_b[1]*axisOfDivision[0];
00340
00341 c_vector<double, SPACE_DIM> moved_centroid = position_a + GetVectorFromAtoB(position_a, centroid);
00342
00343 double alpha = (moved_centroid[0]*a_to_b[1] - position_a[0]*a_to_b[1]
00344 -moved_centroid[1]*a_to_b[0] + position_a[1]*a_to_b[0])/determinant;
00345
00346 intersection = moved_centroid + alpha*axisOfDivision;
00347
00348
00349
00350
00351
00352 c_vector<double, SPACE_DIM> a_to_intersection = this->GetVectorFromAtoB(position_a, intersection);
00353 if (norm_2(a_to_intersection) < mCellRearrangementThreshold)
00354 {
00355 intersection = position_a + mCellRearrangementRatio*mCellRearrangementThreshold*a_to_b/norm_2(a_to_b);
00356 }
00357
00358 c_vector<double, SPACE_DIM> b_to_intersection = this->GetVectorFromAtoB(position_b, intersection);
00359 if (norm_2(b_to_intersection) < mCellRearrangementThreshold)
00360 {
00361 assert(norm_2(a_to_intersection) > mCellRearrangementThreshold);
00362
00363 intersection = position_b - mCellRearrangementRatio*mCellRearrangementThreshold*a_to_b/norm_2(a_to_b);
00364 }
00365 }
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375 bool is_boundary = false;
00376 if (p_node_A->IsBoundaryNode() && p_node_B->IsBoundaryNode())
00377 {
00378 if (elems_containing_node_A.size() !=2 ||
00379 elems_containing_node_B.size() !=2 ||
00380 elems_containing_node_A != elems_containing_node_B)
00381 {
00382 is_boundary = true;
00383 }
00384 }
00385
00386
00387 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, is_boundary, intersection[0], intersection[1]));
00388 nodes_added++;
00389
00390
00391
00392
00393 std::set<unsigned> shared_elements;
00394 std::set_intersection(elems_containing_node_A.begin(),
00395 elems_containing_node_A.end(),
00396 elems_containing_node_B.begin(),
00397 elems_containing_node_B.end(),
00398 std::inserter(shared_elements, shared_elements.begin()));
00399
00400
00401 for (std::set<unsigned>::iterator iter = shared_elements.begin();
00402 iter != shared_elements.end();
00403 ++iter)
00404 {
00405
00406 unsigned local_indexA = this->GetElement(*iter)->GetNodeLocalIndex(p_node_A->GetIndex());
00407 unsigned local_indexB = this->GetElement(*iter)->GetNodeLocalIndex(p_node_B->GetIndex());
00408
00409 unsigned index = local_indexB;
00410 if (local_indexB > local_indexA)
00411 {
00412 index = local_indexA;
00413 }
00414 if ((local_indexA == 0) && (local_indexB == this->GetElement(*iter)->GetNumNodes()-1))
00415 {
00416 index = local_indexB;
00417 }
00418 if ((local_indexB == 0) && (local_indexA == this->GetElement(*iter)->GetNumNodes()-1))
00419 {
00420 index = local_indexA;
00421 }
00422
00423 this->GetElement(*iter)->AddNode(index, this->GetNode(new_node_global_index));
00424 }
00425
00426 division_node_global_indices.push_back(new_node_global_index);
00427 }
00428
00429
00430 unsigned new_element_index = DivideElement(pElement,
00431 pElement->GetNodeLocalIndex(division_node_global_indices[0]),
00432 pElement->GetNodeLocalIndex(division_node_global_indices[1]),
00433 placeOriginalElementBelow);
00434 return new_element_index;
00435 }
00436
00437 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00438 unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DivideElementAlongShortAxis(VertexElement<ELEMENT_DIM,SPACE_DIM>* pElement,
00439 bool placeOriginalElementBelow)
00440 {
00441
00442 assert(SPACE_DIM == 2);
00443 assert(ELEMENT_DIM == SPACE_DIM);
00444
00445
00446 c_vector<double, SPACE_DIM> short_axis = GetShortAxisOfElement(pElement->GetIndex());
00447
00448 unsigned new_element_index = DivideElementAlongGivenAxis(pElement, short_axis, placeOriginalElementBelow);
00449 return new_element_index;
00450 }
00451
00452
00453 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00454 unsigned MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DivideElement(VertexElement<ELEMENT_DIM,SPACE_DIM>* pElement,
00455 unsigned nodeAIndex,
00456 unsigned nodeBIndex,
00457 bool placeOriginalElementBelow)
00458 {
00459
00460 assert(SPACE_DIM == 2);
00461 assert(ELEMENT_DIM == SPACE_DIM);
00462
00463
00464 assert(nodeBIndex != nodeAIndex);
00465
00466 unsigned node1_index = (nodeAIndex < nodeBIndex) ? nodeAIndex : nodeBIndex;
00467 unsigned node2_index = (nodeAIndex < nodeBIndex) ? nodeBIndex : nodeAIndex;
00468
00469
00470 unsigned num_nodes = pElement->GetNumNodes();
00471
00472
00473 std::vector<Node<SPACE_DIM>*> nodes_elem;
00474 for (unsigned i=0; i<num_nodes; i++)
00475 {
00476 nodes_elem.push_back(pElement->GetNode(i));
00477 }
00478
00479
00480 unsigned new_element_index;
00481 if (mDeletedElementIndices.empty())
00482 {
00483 new_element_index = this->mElements.size();
00484 }
00485 else
00486 {
00487 new_element_index = mDeletedElementIndices.back();
00488 mDeletedElementIndices.pop_back();
00489 delete this->mElements[new_element_index];
00490 }
00491
00492
00493 AddElement(new VertexElement<ELEMENT_DIM,SPACE_DIM>(new_element_index, nodes_elem));
00494
00501
00503 double height_midpoint_1 = 0.0;
00504 double height_midpoint_2 = 0.0;
00505 unsigned counter_1 = 0;
00506 unsigned counter_2 = 0;
00507
00508 for (unsigned i=0; i<num_nodes; i++)
00509 {
00510 if (i>=node1_index && i<=node2_index)
00511 {
00512 height_midpoint_1 += pElement->GetNode(i)->rGetLocation()[1];
00513 counter_1++;
00514 }
00515 if (i<=node1_index || i>=node2_index)
00516 {
00517 height_midpoint_2 += pElement->GetNode(i)->rGetLocation()[1];
00518 counter_2++;
00519 }
00520 }
00521 height_midpoint_1 /= (double)counter_1;
00522 height_midpoint_2 /= (double)counter_2;
00523
00524 for (unsigned i=num_nodes; i>0; i--)
00525 {
00526 if (i-1 < node1_index || i-1 > node2_index)
00527 {
00528 if (height_midpoint_1 < height_midpoint_2)
00529 {
00530 if (placeOriginalElementBelow)
00531 {
00532 pElement->DeleteNode(i-1);
00533 }
00534 else
00535 {
00536 this->mElements[new_element_index]->DeleteNode(i-1);
00537 }
00538 }
00539 else
00540 {
00541 if (placeOriginalElementBelow)
00542 {
00543 this->mElements[new_element_index]->DeleteNode(i-1);
00544 }
00545 else
00546 {
00547 pElement->DeleteNode(i-1);
00548 }
00549 }
00550 }
00551 else if (i-1 > node1_index && i-1 < node2_index)
00552 {
00553 if (height_midpoint_1 < height_midpoint_2)
00554 {
00555 if (placeOriginalElementBelow)
00556 {
00557 this->mElements[new_element_index]->DeleteNode(i-1);
00558 }
00559 else
00560 {
00561 pElement->DeleteNode(i-1);
00562 }
00563 }
00564 else
00565 {
00566 if (placeOriginalElementBelow)
00567 {
00568 pElement->DeleteNode(i-1);
00569 }
00570 else
00571 {
00572 this->mElements[new_element_index]->DeleteNode(i-1);
00573 }
00574 }
00575 }
00576 }
00577
00578 return new_element_index;
00579 }
00580
00581
00582 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00583 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DeleteElementPriorToReMesh(unsigned index)
00584 {
00585 assert(SPACE_DIM == 2);
00586
00587
00588 for (unsigned i=0; i<this->mElements[index]->GetNumNodes(); i++)
00589 {
00590 Node<SPACE_DIM>* p_node = this->mElements[index]->GetNode(i);
00591
00592 if (p_node->rGetContainingElementIndices().size()==1)
00593 {
00594 p_node->MarkAsDeleted();
00595 mDeletedNodeIndices.push_back(p_node->GetIndex());
00596 }
00597
00598
00599 p_node->SetAsBoundaryNode(true);
00600 }
00601
00602
00603 this->mElements[index]->MarkAsDeleted();
00604 mDeletedElementIndices.push_back(index);
00605 }
00606
00607
00608 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00609 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::DivideEdge(Node<SPACE_DIM>* pNodeA, Node<SPACE_DIM>* pNodeB)
00610 {
00611
00612 std::set<unsigned> elements_containing_nodeA = pNodeA->rGetContainingElementIndices();
00613 std::set<unsigned> elements_containing_nodeB = pNodeB->rGetContainingElementIndices();
00614
00615
00616 std::set<unsigned> shared_elements;
00617 std::set_intersection(elements_containing_nodeA.begin(),
00618 elements_containing_nodeA.end(),
00619 elements_containing_nodeB.begin(),
00620 elements_containing_nodeB.end(),
00621 std::inserter(shared_elements, shared_elements.begin()));
00622
00623
00624 assert(!shared_elements.empty());
00625 assert(!shared_elements.size()<=2);
00626
00627
00628 bool is_boundary_node = false;
00629 if (shared_elements.size()==1)
00630 {
00631
00632 assert((pNodeA->IsBoundaryNode())&&(pNodeB->IsBoundaryNode()));
00633 is_boundary_node = true;
00634 }
00635
00636
00637 Node<SPACE_DIM>* p_new_node = new Node<SPACE_DIM>(GetNumNodes(), is_boundary_node, 0.0, 0.0);
00638
00639
00640 c_vector<double, SPACE_DIM> new_node_position = pNodeA->rGetLocation() + 0.5*GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation());
00641 ChastePoint<SPACE_DIM> point(new_node_position);
00642 p_new_node->SetPoint(new_node_position);
00643
00644
00645 this->mNodes.push_back(p_new_node);
00646
00647
00648 for (std::set<unsigned>::iterator iter = shared_elements.begin();
00649 iter != shared_elements.end();
00650 ++iter)
00651 {
00652
00653 unsigned local_indexA = this->GetElement(*iter)->GetNodeLocalIndex(pNodeA->GetIndex());
00654 unsigned local_indexB = this->GetElement(*iter)->GetNodeLocalIndex(pNodeB->GetIndex());
00655
00656 unsigned index = local_indexB;
00657 if ( local_indexB > local_indexA )
00658 {
00659 index = local_indexA;
00660 }
00661 if ( (local_indexA == 0) && (local_indexB == this->GetElement(*iter)->GetNumNodes()-1))
00662 {
00663 index = local_indexB;
00664 }
00665 if ( (local_indexB == 0) && (local_indexA == this->GetElement(*iter)->GetNumNodes()-1))
00666 {
00667 index = local_indexA;
00668 }
00669
00670 this->GetElement(*iter)->AddNode(index, p_new_node);
00671 }
00672 }
00673
00674 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00675 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::RemoveDeletedNodesAndElements(VertexElementMap& rElementMap)
00676 {
00677
00678 assert(SPACE_DIM==2 || SPACE_DIM==3);
00679 assert(ELEMENT_DIM == SPACE_DIM);
00680
00681
00682 rElementMap.Resize(this->GetNumAllElements());
00683
00684
00685 std::vector<VertexElement<ELEMENT_DIM, SPACE_DIM>*> live_elements;
00686 for (unsigned i=0; i< this->mElements.size(); i++)
00687 {
00688 if (this->mElements[i]->IsDeleted())
00689 {
00690 delete this->mElements[i];
00691 rElementMap.SetDeleted(i);
00692 }
00693 else
00694 {
00695 live_elements.push_back(this->mElements[i]);
00696 rElementMap.SetNewIndex(i, (unsigned)(live_elements.size()-1));
00697 }
00698 }
00699
00700 assert(mDeletedElementIndices.size() == this->mElements.size() - live_elements.size());
00701 mDeletedElementIndices.clear();
00702 this->mElements = live_elements;
00703
00704 for (unsigned i=0; i<this->mElements.size(); i++)
00705 {
00706 this->mElements[i]->ResetIndex(i);
00707 }
00708
00709
00710 RemoveDeletedNodes();
00711 }
00712
00713 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00714 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::RemoveDeletedNodes()
00715 {
00716 std::vector<Node<SPACE_DIM>*> live_nodes;
00717 for (unsigned i=0; i<this->mNodes.size(); i++)
00718 {
00719 if (this->mNodes[i]->IsDeleted())
00720 {
00721 delete this->mNodes[i];
00722 }
00723 else
00724 {
00725 live_nodes.push_back(this->mNodes[i]);
00726 }
00727 }
00728
00729 assert(mDeletedNodeIndices.size() == this->mNodes.size() - live_nodes.size());
00730 this->mNodes = live_nodes;
00731 mDeletedNodeIndices.clear();
00732
00733 for (unsigned i=0; i<this->mNodes.size(); i++)
00734 {
00735 this->mNodes[i]->SetIndex(i);
00736 }
00737 }
00738
00739
00740 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00741 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ReMesh(VertexElementMap& rElementMap)
00742 {
00743
00744 assert(SPACE_DIM==2 || SPACE_DIM==3);
00745 assert(ELEMENT_DIM == SPACE_DIM);
00746
00747 if (SPACE_DIM==2)
00748 {
00749
00750
00751
00752
00753
00754
00755
00756 RemoveDeletedNodesAndElements(rElementMap);
00757
00758
00759 bool recheck_mesh = true;
00760 while (recheck_mesh == true)
00761 {
00762
00763 recheck_mesh = CheckForT2Swaps(rElementMap);
00764
00765 if (recheck_mesh == false)
00766 {
00767 recheck_mesh = CheckForT1Swaps(rElementMap);
00768 }
00769 }
00770
00771
00772 recheck_mesh = true;
00773 while (recheck_mesh == true)
00774 {
00775
00776 recheck_mesh = CheckForIntersections();
00777 }
00778
00779 RemoveDeletedNodes();
00780 }
00781 else
00782 {
00783 #define COVERAGE_IGNORE
00784 EXCEPTION("Remeshing has not been implemented in 3D (see #827 and #860)\n");
00785 #undef COVERAGE_IGNORE
00787 }
00788 }
00789
00790
00791 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00792 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::ReMesh()
00793 {
00794 VertexElementMap map(GetNumElements());
00795 ReMesh(map);
00796 }
00797
00798
00799 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00800 bool MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::CheckForT1Swaps(VertexElementMap& rElementMap)
00801 {
00802
00803
00804 for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
00805 elem_iter != this->GetElementIteratorEnd();
00806 ++elem_iter)
00807 {
00808 unsigned num_nodes = elem_iter->GetNumNodes();
00809 assert(num_nodes > 0);
00810
00811 unsigned new_num_nodes = num_nodes;
00812
00813
00814
00815
00816
00817
00818
00819
00820 for (unsigned local_index=0; local_index<num_nodes; local_index++)
00821 {
00822
00823 Node<SPACE_DIM>* p_current_node = elem_iter->GetNode(local_index);
00824 unsigned local_index_plus_one = (local_index+1)%new_num_nodes;
00825 Node<SPACE_DIM>* p_anticlockwise_node = elem_iter->GetNode(local_index_plus_one);
00826
00827
00828 double distance_between_nodes = this->GetDistanceBetweenNodes(p_current_node->GetIndex(), p_anticlockwise_node->GetIndex());
00829
00830 std::set<unsigned> elements_of_node_a = p_current_node->rGetContainingElementIndices();
00831 std::set<unsigned> elements_of_node_b = p_anticlockwise_node->rGetContainingElementIndices();
00832
00833 std::set<unsigned> all_elements;
00834 std::set_union(elements_of_node_a.begin(), elements_of_node_a.end(),
00835 elements_of_node_b.begin(), elements_of_node_b.end(),
00836 std::inserter(all_elements, all_elements.begin()));
00837
00838
00839 bool triangular_element = false;
00840
00841 for (std::set<unsigned>::const_iterator it = all_elements.begin();
00842 it != all_elements.end();
00843 ++it)
00844 {
00845 if (this->GetElement(*it)->GetNumNodes() <= 3)
00846 {
00847 triangular_element = true;
00848 }
00849 }
00850
00851
00852 if ((!triangular_element) && (distance_between_nodes < mCellRearrangementThreshold))
00853 {
00854
00855 IdentifySwapType(p_current_node, p_anticlockwise_node, rElementMap);
00856 return true;
00857 }
00858 }
00859 }
00860 return false;
00861 }
00862
00863 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00864 bool MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::CheckForT2Swaps(VertexElementMap& rElementMap)
00865 {
00866
00867
00868 for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
00869 elem_iter != this->GetElementIteratorEnd();
00870 ++elem_iter)
00871 {
00872 if (elem_iter->GetNumNodes() == 3u)
00873 {
00874
00875
00876
00877
00878 if (GetVolumeOfElement(elem_iter->GetIndex()) < GetT2Threshold())
00879 {
00880 PerformT2Swap(*elem_iter);
00881
00882
00883 RemoveDeletedNodesAndElements(rElementMap);
00884
00885 return true;
00886 }
00887 }
00888 }
00889 return false;
00890 }
00891
00892
00893 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00894 bool MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::CheckForIntersections()
00895 {
00896
00897
00898 for (typename AbstractMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator node_iter = this->GetNodeIteratorBegin();
00899 node_iter != this->GetNodeIteratorEnd();
00900 ++node_iter)
00901 {
00902 if (node_iter->IsBoundaryNode())
00903 {
00904 assert(!(node_iter->IsDeleted()));
00905
00906 for (typename VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexElementIterator elem_iter = this->GetElementIteratorBegin();
00907 elem_iter != this->GetElementIteratorEnd();
00908 ++elem_iter)
00909 {
00910 if (elem_iter->IsElementOnBoundary())
00911 {
00912
00913 unsigned elem_index = elem_iter->GetIndex();
00914
00915
00916 if (node_iter->rGetContainingElementIndices().count(elem_index) == 0)
00917 {
00918 if (ElementIncludesPoint(node_iter->rGetLocation(), elem_index))
00919 {
00920 PerformT3Swap(&(*node_iter), elem_index);
00921 return true;
00922 }
00923 }
00924 }
00925 }
00926 }
00927 }
00928
00929 return false;
00930 }
00931
00932 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00933 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::IdentifySwapType(Node<SPACE_DIM>* pNodeA, Node<SPACE_DIM>* pNodeB, VertexElementMap& rElementMap)
00934 {
00935
00936 assert(SPACE_DIM == 2);
00937 assert(ELEMENT_DIM == SPACE_DIM);
00938
00939
00940 std::set<unsigned> nodeA_elem_indices = pNodeA->rGetContainingElementIndices();
00941
00942 std::set<unsigned> nodeB_elem_indices = pNodeB->rGetContainingElementIndices();
00943
00944
00945 std::set<unsigned> all_indices, temp_set;
00946 std::set_union(nodeA_elem_indices.begin(), nodeA_elem_indices.end(),
00947 nodeB_elem_indices.begin(), nodeB_elem_indices.end(),
00948 std::inserter(temp_set, temp_set.begin()));
00949 all_indices.swap(temp_set);
00950
00951 if ((nodeA_elem_indices.size()>3) || (nodeB_elem_indices.size()>3))
00952 {
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963 EXCEPTION("A node is contained in more than three elements");
00964 }
00965 else
00966 {
00967 switch (all_indices.size())
00968 {
00969 case 1:
00970 {
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982 assert(pNodeA->IsBoundaryNode());
00983 assert(pNodeB->IsBoundaryNode());
00984
00985 PerformNodeMerge(pNodeA, pNodeB);
00986
00987
00988 RemoveDeletedNodes();
00989 break;
00990 }
00991 case 2:
00992 {
00993 if (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==2)
00994 {
00995 if (pNodeA->IsBoundaryNode() && pNodeB->IsBoundaryNode())
00996 {
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010 PerformT1Swap(pNodeA, pNodeB,all_indices);
01011 }
01012 else if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
01013 {
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030 EXCEPTION("There is a non boundary node contained only in 2 elements something has gone wrong.");
01031 }
01032 else
01033 {
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045 EXCEPTION("There are non-boundary nodes contained in only in 2 elements something has gone wrong.");
01046 }
01047 }
01048 else
01049 {
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061 PerformNodeMerge(pNodeA, pNodeB);
01062
01063
01064 RemoveDeletedNodes();
01065 }
01066 break;
01067 }
01068 case 3:
01069 {
01070 if (nodeA_elem_indices.size()==1 || nodeB_elem_indices.size()==1)
01071 {
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089 assert(pNodeA->IsBoundaryNode());
01090 assert(pNodeB->IsBoundaryNode());
01091
01092 EXCEPTION("There is a boundary node contained in three elements something has gone wrong.");
01093 }
01094 else if (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==2)
01095 {
01096
01097
01098 std::set<unsigned> element_A_not_B, temp_set;
01099 std::set_difference(all_indices.begin(), all_indices.end(),
01100 nodeB_elem_indices.begin(), nodeB_elem_indices.end(),
01101 std::inserter(temp_set, temp_set.begin()));
01102 element_A_not_B.swap(temp_set);
01103
01104
01105 assert(element_A_not_B.size()==1);
01106
01107 std::set<unsigned> element_B_not_A;
01108 std::set_difference(all_indices.begin(), all_indices.end(),
01109 nodeA_elem_indices.begin(), nodeA_elem_indices.end(),
01110 std::inserter(temp_set, temp_set.begin()));
01111 element_B_not_A.swap(temp_set);
01112
01113
01114 assert(element_B_not_A.size()==1);
01115
01116 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_A_not_B = this->mElements[*element_A_not_B.begin()];
01117 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_B_not_A = this->mElements[*element_B_not_A.begin()];
01118
01119 unsigned local_index_1 = p_element_A_not_B->GetNodeLocalIndex(pNodeA->GetIndex());
01120 unsigned next_node_1 = p_element_A_not_B->GetNodeGlobalIndex((local_index_1 + 1)%(p_element_A_not_B->GetNumNodes()));
01121 unsigned previous_node_1 = p_element_A_not_B->GetNodeGlobalIndex((local_index_1 + p_element_A_not_B->GetNumNodes() - 1)%(p_element_A_not_B->GetNumNodes()));
01122 unsigned local_index_2 = p_element_B_not_A->GetNodeLocalIndex(pNodeB->GetIndex());
01123 unsigned next_node_2 = p_element_B_not_A->GetNodeGlobalIndex((local_index_2 + 1)%(p_element_B_not_A->GetNumNodes()));
01124 unsigned previous_node_2 = p_element_B_not_A->GetNodeGlobalIndex((local_index_2 + p_element_B_not_A->GetNumNodes() - 1)%(p_element_B_not_A->GetNumNodes()));
01125
01126 if (next_node_1 == previous_node_2 || next_node_2 == previous_node_1)
01127 {
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143 assert(pNodeA->IsBoundaryNode());
01144 assert(pNodeB->IsBoundaryNode());
01145
01146
01147 unsigned nodeC_index;
01148 if (next_node_1 == previous_node_2 && next_node_2 != previous_node_1)
01149 {
01150 nodeC_index = next_node_1;
01151 }
01152 else if (next_node_2 == previous_node_1 && next_node_1 != previous_node_2)
01153 {
01154 nodeC_index = next_node_2;
01155 }
01156 else
01157 {
01158 assert(next_node_1 == previous_node_2 && next_node_2 == previous_node_1);
01159 EXCEPTION("Triangular element next to triangular void, not implemented yet.");
01160 }
01161
01162 PerformVoidRemoval(pNodeA, pNodeB, this->mNodes[nodeC_index]);
01163 }
01164 else
01165 {
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180 assert(pNodeA->IsBoundaryNode());
01181 assert(pNodeB->IsBoundaryNode());
01182
01183 PerformT1Swap(pNodeA, pNodeB, all_indices);
01184 }
01185 }
01186 else
01187 {
01188
01189
01190
01191
01192 assert ( (nodeA_elem_indices.size()==2 && nodeB_elem_indices.size()==3)
01193 || (nodeA_elem_indices.size()==3 && nodeB_elem_indices.size()==2) );
01194
01195
01196 assert(!(pNodeA->IsBoundaryNode() && pNodeB->IsBoundaryNode()));
01197
01198 if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
01199 {
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212 PerformT1Swap(pNodeA, pNodeB, all_indices);
01213 }
01214 else
01215 {
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230 EXCEPTION("There are non-boundary nodes contained only in 2 elements something has gone wrong.");
01231 }
01232 }
01233 break;
01234 }
01235 case 4:
01236 {
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248 PerformT1Swap(pNodeA, pNodeB, all_indices);
01249 break;
01250 }
01251 default:
01252
01253 NEVER_REACHED;
01254 }
01255 }
01256 }
01257
01258 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01259 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformNodeMerge(Node<SPACE_DIM>* pNodeA, Node<SPACE_DIM>* pNodeB)
01260 {
01261
01262 unsigned nodeA_index = pNodeA->GetIndex();
01263 unsigned nodeB_index = pNodeB->GetIndex();
01264
01265 unsigned lo_node_index = (nodeA_index < nodeB_index) ? nodeA_index : nodeB_index;
01266 unsigned hi_node_index = (nodeA_index < nodeB_index) ? nodeB_index : nodeA_index;
01267
01268
01269 Node<SPACE_DIM>* p_lo_node = this->GetNode(lo_node_index);
01270 Node<SPACE_DIM>* p_hi_node = this->GetNode(hi_node_index);
01271
01272
01273 std::set<unsigned> lo_node_elem_indices = p_lo_node->rGetContainingElementIndices();
01274 std::set<unsigned> hi_node_elem_indices = p_hi_node->rGetContainingElementIndices();
01275
01276
01277 c_vector<double, SPACE_DIM> node_midpoint = p_lo_node->rGetLocation() + 0.5*this->GetVectorFromAtoB(p_lo_node->rGetLocation(), p_hi_node->rGetLocation());
01278 c_vector<double, SPACE_DIM>& r_lo_node_location = p_lo_node->rGetModifiableLocation();
01279 r_lo_node_location = node_midpoint;
01280
01281
01282 for (std::set<unsigned>::const_iterator it = hi_node_elem_indices.begin();
01283 it != hi_node_elem_indices.end();
01284 ++it)
01285 {
01286
01287 unsigned hi_node_local_index = this->mElements[*it]->GetNodeLocalIndex(hi_node_index);
01288 assert(hi_node_local_index < UINT_MAX);
01289
01290
01291
01292
01293
01294 if (lo_node_elem_indices.count(*it) > 0)
01295 {
01296 this->mElements[*it]->DeleteNode(hi_node_local_index);
01297 }
01298 else
01299 {
01300
01301 this->mElements[*it]->UpdateNode(hi_node_local_index, p_lo_node);
01302 }
01303 }
01304 assert(!(this->mNodes[hi_node_index]->IsDeleted()));
01305 this->mNodes[hi_node_index]->MarkAsDeleted();
01306 mDeletedNodeIndices.push_back(hi_node_index);
01307 }
01308
01309 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01310 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformT1Swap(Node<SPACE_DIM>* pNodeA,
01311 Node<SPACE_DIM>* pNodeB,
01312 std::set<unsigned>& rElementsContainingNodes)
01313 {
01314
01315 assert(SPACE_DIM == 2);
01316 assert(ELEMENT_DIM == SPACE_DIM);
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352 std::set<unsigned> nodeA_elem_indices = pNodeA->rGetContainingElementIndices();
01353 std::set<unsigned> nodeB_elem_indices = pNodeB->rGetContainingElementIndices();
01354
01355 double distance_between_nodes_CD = mCellRearrangementRatio*mCellRearrangementThreshold;
01356
01357 c_vector<double, SPACE_DIM> nodeA_location = pNodeA->rGetLocation();
01358 c_vector<double, SPACE_DIM> nodeB_location = pNodeB->rGetLocation();
01359
01360 c_vector<double, SPACE_DIM> a_to_b = this->GetVectorFromAtoB(nodeA_location, nodeB_location);
01361
01362
01363 mLocationsOfT1Swaps.push_back(nodeA_location + 0.5*a_to_b);
01364
01365 c_vector<double, SPACE_DIM> perpendicular_vector;
01366 perpendicular_vector(0) = -a_to_b(1);
01367 perpendicular_vector(1) = a_to_b(0);
01368
01369 if (norm_2(a_to_b) < 1e-10)
01370 {
01371 EXCEPTION("Nodes are too close together, this shouldn't happen");
01372 }
01373
01374 c_vector<double, SPACE_DIM> c_to_d = distance_between_nodes_CD / norm_2(a_to_b) * perpendicular_vector;
01375 c_vector<double, SPACE_DIM> nodeC_location = nodeA_location + 0.5*a_to_b - 0.5*c_to_d;
01376 c_vector<double, SPACE_DIM> nodeD_location = nodeC_location + c_to_d;
01377
01378
01379
01380
01381
01382 c_vector<double, SPACE_DIM>& r_nodeA_location = pNodeA->rGetModifiableLocation();
01383 r_nodeA_location = nodeC_location;
01384
01385 c_vector<double, SPACE_DIM>& r_nodeB_location = pNodeB->rGetModifiableLocation();
01386 r_nodeB_location = nodeD_location;
01387
01388 for (std::set<unsigned>::const_iterator it = rElementsContainingNodes.begin();
01389 it != rElementsContainingNodes.end();
01390 ++it)
01391 {
01392 if (nodeA_elem_indices.find(*it) == nodeA_elem_indices.end())
01393 {
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404 unsigned nodeB_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeB->GetIndex());
01405 assert(nodeB_local_index < UINT_MAX);
01406
01407 this->mElements[*it]->AddNode(nodeB_local_index, pNodeA);
01408 }
01409 else if (nodeB_elem_indices.find(*it) == nodeB_elem_indices.end())
01410 {
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420 unsigned nodeA_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeA->GetIndex());
01421 assert(nodeA_local_index < UINT_MAX);
01422 this->mElements[*it]->AddNode(nodeA_local_index, pNodeB);
01423 }
01424 else
01425 {
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436 unsigned nodeA_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeA->GetIndex());
01437 assert(nodeA_local_index < UINT_MAX);
01438
01439 unsigned nodeB_local_index = this->mElements[*it]->GetNodeLocalIndex(pNodeB->GetIndex());
01440 assert(nodeB_local_index < UINT_MAX);
01441
01442 unsigned nodeB_local_index_plus_one = (nodeB_local_index + 1)%(this->mElements[*it]->GetNumNodes());
01443
01444 if (nodeA_local_index == nodeB_local_index_plus_one)
01445 {
01446
01447
01448
01449
01450 this->mElements[*it]->DeleteNode(nodeB_local_index);
01451 }
01452 else
01453 {
01454 assert(nodeB_local_index == (nodeA_local_index + 1)%(this->mElements[*it]->GetNumNodes()));
01455
01456
01457
01458
01459 this->mElements[*it]->DeleteNode(nodeA_local_index);
01460 }
01461 }
01462 }
01463
01464
01465 if (pNodeA->IsBoundaryNode() || pNodeB->IsBoundaryNode())
01466 {
01467 if (pNodeA->GetNumContainingElements()==3)
01468 {
01469 pNodeA->SetAsBoundaryNode(false);
01470 }
01471 else
01472 {
01473 pNodeA->SetAsBoundaryNode(true);
01474 }
01475 if (pNodeB->GetNumContainingElements()==3)
01476 {
01477 pNodeB->SetAsBoundaryNode(false);
01478 }
01479 else
01480 {
01481 pNodeB->SetAsBoundaryNode(true);
01482 }
01483 }
01484 }
01485
01486 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01487 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformT2Swap(VertexElement<ELEMENT_DIM,SPACE_DIM>& rElement)
01488 {
01489
01490 assert(SPACE_DIM == 2);
01491 assert(ELEMENT_DIM == SPACE_DIM);
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523 assert(rElement.GetNumNodes() == 3u);
01524
01525 bool is_node_on_boundary = false;
01526 for (unsigned i=0; i<3; i++)
01527 {
01528 if (rElement.GetNode(i)->IsBoundaryNode())
01529 {
01530 is_node_on_boundary= true;
01531 }
01532 }
01533
01534
01535 c_vector<double, SPACE_DIM> new_node_location = GetCentroidOfElement(rElement.GetIndex());
01536 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(GetNumNodes(), new_node_location, is_node_on_boundary));
01537 Node<SPACE_DIM>* p_new_node = this->GetNode(new_node_global_index);
01538
01539 std::set<unsigned> neighbouring_elements;
01540
01541 for (unsigned i=0; i<3; i++)
01542 {
01543 Node<SPACE_DIM>* p_node_a = rElement.GetNode((i+1)%3);
01544 Node<SPACE_DIM>* p_node_b = rElement.GetNode((i+2)%3);
01545
01546 std::set<unsigned> elements_of_node_a = p_node_a->rGetContainingElementIndices();
01547 std::set<unsigned> elements_of_node_b = p_node_b->rGetContainingElementIndices();
01548
01549 std::set<unsigned> common_elements;
01550
01551 std::set_intersection( elements_of_node_a.begin(), elements_of_node_a.end(),
01552 elements_of_node_b.begin(), elements_of_node_b.end(),
01553 std::inserter(common_elements, common_elements.begin()));
01554
01555 assert(common_elements.size() <= 2u);
01556 common_elements.erase(rElement.GetIndex());
01557
01558 if (common_elements.size() == 1u)
01559 {
01560 VertexElement<ELEMENT_DIM,SPACE_DIM>* p_neighbouring_element = this->GetElement(*(common_elements.begin()));
01561
01562 if (p_neighbouring_element->GetNumNodes() < 4u)
01563 {
01564 EXCEPTION("One of the neighbours of a small triangular element is also a triangle - dealing with this has not been implemented yet");
01565 }
01566
01567 p_neighbouring_element-> ReplaceNode(p_node_a, p_new_node);
01568 p_neighbouring_element->DeleteNode(p_neighbouring_element->GetNodeLocalIndex(p_node_b->GetIndex()));
01569 }
01570 else
01571 {
01572 assert( (p_node_a->IsBoundaryNode()) && (p_node_b->IsBoundaryNode()) );
01573 }
01574 }
01575
01576
01577 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(0));
01578 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(1));
01579 mDeletedNodeIndices.push_back(rElement.GetNodeGlobalIndex(2));
01580 rElement.GetNode(0)->MarkAsDeleted();
01581 rElement.GetNode(1)->MarkAsDeleted();
01582 rElement.GetNode(2)->MarkAsDeleted();
01583
01584 mDeletedElementIndices.push_back(rElement.GetIndex());
01585 rElement.MarkAsDeleted();
01586 }
01587
01588 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01589 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformT3Swap(Node<SPACE_DIM>* pNode, unsigned elementIndex)
01590 {
01591
01592 assert(SPACE_DIM == 2);
01593 assert(ELEMENT_DIM == SPACE_DIM);
01594
01595
01596 assert(pNode->IsBoundaryNode());
01597
01598
01599 std::set<unsigned> elements_containing_intersecting_node = pNode->rGetContainingElementIndices();
01600
01601
01602 unsigned node_A_local_index = GetLocalIndexForElementEdgeClosestToPoint(pNode->rGetLocation(), elementIndex);
01603
01604
01605 c_vector<double, SPACE_DIM> node_location = pNode->rGetModifiableLocation();
01606
01607
01608 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elementIndex);
01609 unsigned num_nodes = p_element->GetNumNodes();
01610
01611
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01647
01648
01649
01650 unsigned vertexA_index = p_element->GetNodeGlobalIndex(node_A_local_index);
01651 unsigned vertexB_index = p_element->GetNodeGlobalIndex((node_A_local_index+1)%num_nodes);
01652
01653
01654 if (!this->mNodes[vertexA_index]->IsBoundaryNode() || !this->mNodes[vertexB_index]->IsBoundaryNode())
01655 {
01656 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.");
01657 }
01658
01659
01660 c_vector<double, SPACE_DIM> vertexA = p_element->GetNodeLocation(node_A_local_index);
01661 c_vector<double, SPACE_DIM> vertexB = p_element->GetNodeLocation((node_A_local_index+1)%num_nodes);
01662 c_vector<double, SPACE_DIM> vector_a_to_point = this->GetVectorFromAtoB(vertexA, node_location);
01663
01664 c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
01665
01666 c_vector<double, SPACE_DIM> edge_ab_unit_vector = vector_a_to_b/norm_2(vector_a_to_b);
01667 c_vector<double, SPACE_DIM> intersection = vertexA + edge_ab_unit_vector*inner_prod(vector_a_to_point, edge_ab_unit_vector);
01668
01669
01670 mLocationsOfT3Swaps.push_back(intersection);
01671
01672
01673
01674
01675
01676
01677
01678
01679 if (norm_2(vector_a_to_b) < 4.0*mCellRearrangementRatio*mCellRearrangementThreshold)
01680 {
01681 WARNING("Trying to merge a node onto an edge which is too small.");
01682
01683 c_vector<double, SPACE_DIM> centre_a_and_b = vertexA + 0.5*vector_a_to_b;
01684
01685 vertexA = centre_a_and_b - 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*vector_a_to_b/norm_2(vector_a_to_b);
01686 ChastePoint<SPACE_DIM> vertex_A_point(vertexA);
01687 SetNode(p_element->GetNodeGlobalIndex(node_A_local_index), vertex_A_point);
01688
01689 vertexB = centre_a_and_b + 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*vector_a_to_b/norm_2(vector_a_to_b);
01690 ChastePoint<SPACE_DIM> vertex_B_point(vertexB);
01691 SetNode(p_element->GetNodeGlobalIndex((node_A_local_index+1)%num_nodes), vertex_B_point);
01692
01693
01694 vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
01695 edge_ab_unit_vector = vector_a_to_b/norm_2(vector_a_to_b);
01696
01697
01698 intersection = centre_a_and_b;
01699 }
01700
01701
01702
01703
01704
01705
01706
01707
01708 if (norm_2(intersection - vertexA) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
01709 {
01710 intersection = vertexA + 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
01711 }
01712 if (norm_2(intersection - vertexB) < 2.0*mCellRearrangementRatio*mCellRearrangementThreshold)
01713 {
01714 intersection = vertexB - 2.0*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
01715 }
01716
01717 if (pNode->GetNumContainingElements() == 1)
01718 {
01719
01720 unsigned intersecting_element_index = *elements_containing_intersecting_node.begin();
01721
01722
01723 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_intersecting_element = this->GetElement(intersecting_element_index);
01724
01725 unsigned local_index = p_intersecting_element->GetNodeLocalIndex(pNode->GetIndex());
01726 unsigned next_node = p_intersecting_element->GetNodeGlobalIndex((local_index + 1)%(p_intersecting_element->GetNumNodes()));
01727 unsigned previous_node = p_intersecting_element->GetNodeGlobalIndex((local_index + p_intersecting_element->GetNumNodes() - 1)%(p_intersecting_element->GetNumNodes()));
01728
01729
01730 if (next_node == vertexA_index || previous_node == vertexA_index || next_node == vertexB_index || previous_node == vertexB_index)
01731 {
01732 unsigned common_vertex_index;
01733
01734 if (next_node == vertexA_index || previous_node == vertexA_index)
01735 {
01736 common_vertex_index = vertexA_index;
01737 }
01738 else
01739 {
01740 common_vertex_index = vertexB_index;
01741 }
01742
01743 assert(this->mNodes[common_vertex_index]->GetNumContainingElements()>1);
01744
01745 std::set<unsigned> elements_containing_common_vertex = this->mNodes[common_vertex_index]->rGetContainingElementIndices();
01746 std::set<unsigned>::const_iterator it = elements_containing_common_vertex.begin();
01747 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*it);
01748 it++;
01749 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*it);
01750
01751
01752 unsigned num_common_vertices = 0;
01753 for (unsigned i=0; i<p_element_common_1->GetNumNodes(); i++)
01754 {
01755 for (unsigned j=0; j<p_element_common_2->GetNumNodes(); j++)
01756 {
01757 if (p_element_common_1->GetNodeGlobalIndex(i)==p_element_common_2->GetNodeGlobalIndex(j))
01758 {
01759 num_common_vertices++;
01760 }
01761 }
01762 }
01763
01764 if (num_common_vertices == 1 || this->mNodes[common_vertex_index]->GetNumContainingElements() > 2)
01765 {
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780 pNode->rGetModifiableLocation() = intersection;
01781
01782
01783 this->GetElement(elementIndex)->AddNode(node_A_local_index, pNode);
01784
01785
01786 assert(pNode->GetNumContainingElements() == 2);
01787 }
01788 else if (num_common_vertices == 2)
01789 {
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807 pNode->rGetModifiableLocation() = intersection;
01808
01809
01810 this->GetElement(elementIndex)->ReplaceNode(this->mNodes[common_vertex_index], pNode);
01811
01812
01813 this->GetElement(intersecting_element_index)->DeleteNode(this->GetElement(intersecting_element_index)->GetNodeLocalIndex(common_vertex_index));
01814 assert(this->mNodes[common_vertex_index]->GetNumContainingElements()==0);
01815
01816 this->mNodes[common_vertex_index]->MarkAsDeleted();
01817 mDeletedNodeIndices.push_back(common_vertex_index);
01818
01819
01820 assert(pNode->GetNumContainingElements() == 2);
01821 }
01822 else
01823 {
01824
01825 NEVER_REACHED;
01826 }
01827 }
01828 else
01829 {
01830
01831
01832
01833
01834
01835
01836
01837
01838
01839
01840
01841 pNode->rGetModifiableLocation() = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
01842
01843 c_vector<double, SPACE_DIM> new_node_location = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
01844
01845
01846 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
01847
01848
01849 this->GetElement(elementIndex)->AddNode(node_A_local_index, pNode);
01850 this->GetElement(elementIndex)->AddNode(node_A_local_index, this->mNodes[new_node_global_index]);
01851
01852
01853 this->GetElement(intersecting_element_index)->AddNode(this->GetElement(intersecting_element_index)->GetNodeLocalIndex(pNode->GetIndex()), this->mNodes[new_node_global_index]);
01854
01855
01856 assert(pNode->GetNumContainingElements() == 2);
01857 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
01858 }
01859 }
01860 else if (pNode->GetNumContainingElements() == 2)
01861 {
01862
01863 std::set<unsigned>::const_iterator it = elements_containing_intersecting_node.begin();
01864 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_intersection_1 = this->GetElement(*it);
01865 it++;
01866 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_intersection_2 = this->GetElement(*it);
01867
01868 unsigned local_index_1 = p_element_intersection_1->GetNodeLocalIndex(pNode->GetIndex());
01869 unsigned next_node_1 = p_element_intersection_1->GetNodeGlobalIndex((local_index_1 + 1)%(p_element_intersection_1->GetNumNodes()));
01870 unsigned previous_node_1 = p_element_intersection_1->GetNodeGlobalIndex((local_index_1 + p_element_intersection_1->GetNumNodes() - 1)%(p_element_intersection_1->GetNumNodes()));
01871 unsigned local_index_2 = p_element_intersection_2->GetNodeLocalIndex(pNode->GetIndex());
01872 unsigned next_node_2 = p_element_intersection_2->GetNodeGlobalIndex((local_index_2 + 1)%(p_element_intersection_2->GetNumNodes()));
01873 unsigned previous_node_2 = p_element_intersection_2->GetNodeGlobalIndex((local_index_2 + p_element_intersection_2->GetNumNodes() - 1)%(p_element_intersection_2->GetNumNodes()));
01874
01875
01876 if ((next_node_1 == vertexA_index || previous_node_1 == vertexA_index || next_node_2 == vertexA_index || previous_node_2 == vertexA_index) &&
01877 (next_node_1 == vertexB_index || previous_node_1 == vertexB_index || next_node_2 == vertexB_index || previous_node_2 == vertexB_index))
01878 {
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893 assert(pNode->IsBoundaryNode());
01894 assert(this->mNodes[vertexA_index]->IsBoundaryNode());
01895 assert(this->mNodes[vertexB_index]->IsBoundaryNode());
01896
01897
01898 pNode->rGetModifiableLocation() = intersection;
01899 pNode->SetAsBoundaryNode(false);
01900
01901
01902 this->GetElement(elementIndex)->AddNode(node_A_local_index, pNode);
01903
01904
01905 std::set<unsigned> elements_containing_vertex_A = this->mNodes[vertexA_index]->rGetContainingElementIndices();
01906 for (std::set<unsigned>::const_iterator iter = elements_containing_vertex_A.begin();
01907 iter != elements_containing_vertex_A.end();
01908 iter++)
01909 {
01910 this->GetElement(*iter)->DeleteNode(this->GetElement(*iter)->GetNodeLocalIndex(vertexA_index));
01911 }
01912
01913
01914 assert(this->mNodes[vertexA_index]->GetNumContainingElements()==0);
01915
01916 this->mNodes[vertexA_index]->MarkAsDeleted();
01917 mDeletedNodeIndices.push_back(vertexA_index);
01918
01919
01920 std::set<unsigned> elements_containing_vertex_B = this->mNodes[vertexB_index]->rGetContainingElementIndices();
01921 for (std::set<unsigned>::const_iterator iter = elements_containing_vertex_B.begin();
01922 iter != elements_containing_vertex_B.end();
01923 iter++)
01924 {
01925 this->GetElement(*iter)->DeleteNode(this->GetElement(*iter)->GetNodeLocalIndex(vertexB_index));
01926 }
01927
01928
01929 assert(this->mNodes[vertexB_index]->GetNumContainingElements()==0);
01930
01931 this->mNodes[vertexB_index]->MarkAsDeleted();
01932 mDeletedNodeIndices.push_back(vertexB_index);
01933 }
01934 else
01935 {
01936 if (next_node_1 == vertexA_index || previous_node_1 == vertexA_index || next_node_2 == vertexA_index || previous_node_2 == vertexA_index)
01937 {
01938
01939
01940 assert(this->mNodes[vertexA_index]->GetNumContainingElements()>1);
01941
01942 std::set<unsigned> elements_containing_vertex_A = this->mNodes[vertexA_index]->rGetContainingElementIndices();
01943 std::set<unsigned>::const_iterator iter = elements_containing_vertex_A.begin();
01944 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*iter);
01945 iter++;
01946 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*iter);
01947
01948
01949 unsigned num_common_vertices = 0;
01950 for (unsigned i=0; i<p_element_common_1->GetNumNodes(); i++)
01951 {
01952 for (unsigned j=0; j<p_element_common_2->GetNumNodes(); j++)
01953 {
01954 if (p_element_common_1->GetNodeGlobalIndex(i) == p_element_common_2->GetNodeGlobalIndex(j))
01955 {
01956 num_common_vertices++;
01957 }
01958 }
01959 }
01960
01961 if (num_common_vertices == 1 || this->mNodes[vertexA_index]->GetNumContainingElements() > 2)
01962 {
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977 pNode->rGetModifiableLocation() = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
01978 pNode->SetAsBoundaryNode(false);
01979
01980 c_vector<double, SPACE_DIM> new_node_location = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
01981
01982
01983 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
01984
01985
01986 this->GetElement(elementIndex)->AddNode(node_A_local_index, this->mNodes[new_node_global_index]);
01987 this->GetElement(elementIndex)->AddNode(node_A_local_index, pNode);
01988
01989
01990 if (next_node_1 == previous_node_2)
01991 {
01992 p_element_intersection_1->AddNode((local_index_1 + p_element_intersection_1->GetNumNodes() - 1)%(p_element_intersection_1->GetNumNodes()), this->mNodes[new_node_global_index]);
01993 }
01994 else
01995 {
01996 assert(next_node_2 == previous_node_1);
01997
01998 p_element_intersection_2->AddNode((local_index_2 + p_element_intersection_2->GetNumNodes() - 1)%(p_element_intersection_2->GetNumNodes()), this->mNodes[new_node_global_index]);
01999 }
02000
02001
02002 assert(pNode->GetNumContainingElements() == 3);
02003 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
02004 }
02005 else if (num_common_vertices == 2)
02006 {
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022 pNode->rGetModifiableLocation() = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
02023 pNode->SetAsBoundaryNode(false);
02024
02025 c_vector<double, SPACE_DIM> new_node_location = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
02026
02027
02028 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
02029
02030
02031 this->GetElement(elementIndex)->AddNode(node_A_local_index, this->mNodes[new_node_global_index]);
02032 this->GetElement(elementIndex)->AddNode(node_A_local_index, pNode);
02033
02034
02035 if (next_node_1 == previous_node_2)
02036 {
02037 p_element_intersection_1->AddNode((local_index_1 + p_element_intersection_1->GetNumNodes() - 1)%(p_element_intersection_1->GetNumNodes()), this->mNodes[new_node_global_index]);
02038 }
02039 else
02040 {
02041 assert(next_node_2 == previous_node_1);
02042
02043 p_element_intersection_2->AddNode((local_index_2 + p_element_intersection_2->GetNumNodes() - 1)%(p_element_intersection_2->GetNumNodes()), this->mNodes[new_node_global_index]);
02044 }
02045
02046
02047 p_element_common_1->DeleteNode(p_element_common_1->GetNodeLocalIndex(vertexA_index));
02048 p_element_common_2->DeleteNode(p_element_common_2->GetNodeLocalIndex(vertexA_index));
02049
02050 assert(this->mNodes[vertexA_index]->GetNumContainingElements()==0);
02051
02052 this->mNodes[vertexA_index]->MarkAsDeleted();
02053 mDeletedNodeIndices.push_back(vertexA_index);
02054
02055
02056 assert(pNode->GetNumContainingElements() == 3);
02057 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
02058 }
02059 else
02060 {
02061
02062 NEVER_REACHED;
02063 }
02064
02065 }
02066 else if (next_node_1 == vertexB_index || previous_node_1 == vertexB_index || next_node_2 == vertexB_index || previous_node_2 == vertexB_index)
02067 {
02068
02069
02070 assert(this->mNodes[vertexB_index]->GetNumContainingElements()>1);
02071
02072 std::set<unsigned> elements_containing_vertex_B = this->mNodes[vertexB_index]->rGetContainingElementIndices();
02073 std::set<unsigned>::const_iterator iter = elements_containing_vertex_B.begin();
02074 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_1 = this->GetElement(*iter);
02075 iter++;
02076 VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element_common_2 = this->GetElement(*iter);
02077
02078
02079 unsigned num_common_vertices = 0;
02080 for (unsigned i=0; i<p_element_common_1->GetNumNodes(); i++)
02081 {
02082 for (unsigned j=0; j<p_element_common_2->GetNumNodes(); j++)
02083 {
02084 if (p_element_common_1->GetNodeGlobalIndex(i)==p_element_common_2->GetNodeGlobalIndex(j))
02085 {
02086 num_common_vertices++;
02087 }
02088 }
02089 }
02090
02091 if (num_common_vertices == 1 || this->mNodes[vertexB_index]->GetNumContainingElements() > 2)
02092 {
02093
02094
02095
02096
02097
02098
02099
02100
02101
02102
02103
02104
02105
02106
02107 pNode->rGetModifiableLocation() = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
02108 pNode->SetAsBoundaryNode(false);
02109
02110 c_vector<double, SPACE_DIM> new_node_location = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
02111
02112
02113 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
02114
02115
02116 this->GetElement(elementIndex)->AddNode(node_A_local_index, pNode);
02117 this->GetElement(elementIndex)->AddNode(node_A_local_index, this->mNodes[new_node_global_index]);
02118
02119
02120 if (next_node_1 == previous_node_2)
02121 {
02122 p_element_intersection_2->AddNode(local_index_2, this->mNodes[new_node_global_index]);
02123 }
02124 else
02125 {
02126 assert(next_node_2 == previous_node_1);
02127
02128 p_element_intersection_1->AddNode(local_index_1, this->mNodes[new_node_global_index]);
02129 }
02130
02131
02132 assert(pNode->GetNumContainingElements() == 3);
02133 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
02134 }
02135 else if (num_common_vertices == 2)
02136 {
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152 pNode->rGetModifiableLocation() = intersection + 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
02153 pNode->SetAsBoundaryNode(false);
02154
02155 c_vector<double, SPACE_DIM> new_node_location = intersection - 0.5*mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
02156
02157
02158 unsigned new_node_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_location[0], new_node_location[1]));
02159
02160
02161 this->GetElement(elementIndex)->AddNode(node_A_local_index, pNode);
02162 this->GetElement(elementIndex)->AddNode(node_A_local_index, this->mNodes[new_node_global_index]);
02163
02164
02165 if (next_node_1 == previous_node_2)
02166 {
02167 p_element_intersection_2->AddNode(local_index_2, this->mNodes[new_node_global_index]);
02168 }
02169 else
02170 {
02171 assert(next_node_2 == previous_node_1);
02172
02173 p_element_intersection_1->AddNode(local_index_1, this->mNodes[new_node_global_index]);
02174 }
02175
02176
02177 p_element_common_1->DeleteNode(p_element_common_1->GetNodeLocalIndex(vertexB_index));
02178 p_element_common_2->DeleteNode(p_element_common_2->GetNodeLocalIndex(vertexB_index));
02179
02180 assert(this->mNodes[vertexB_index]->GetNumContainingElements()==0);
02181
02182 this->mNodes[vertexB_index]->MarkAsDeleted();
02183 mDeletedNodeIndices.push_back(vertexB_index);
02184
02185
02186 assert(pNode->GetNumContainingElements() == 3);
02187 assert(this->mNodes[new_node_global_index]->GetNumContainingElements() == 2);
02188 }
02189 else
02190 {
02191
02192 NEVER_REACHED;
02193 }
02194 }
02195 else
02196 {
02197
02198
02199
02200
02201
02202
02203
02204
02205
02206
02207
02208 pNode->rGetModifiableLocation() = intersection;
02209 pNode->SetAsBoundaryNode(false);
02210
02211 c_vector<double, SPACE_DIM> new_node_1_location = intersection - mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
02212 c_vector<double, SPACE_DIM> new_node_2_location = intersection + mCellRearrangementRatio*mCellRearrangementThreshold*edge_ab_unit_vector;
02213
02214
02215 unsigned new_node_1_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_1_location[0], new_node_1_location[1]));
02216 unsigned new_node_2_global_index = this->AddNode(new Node<SPACE_DIM>(0, true, new_node_2_location[0], new_node_2_location[1]));
02217
02218
02219 this->GetElement(elementIndex)->AddNode(node_A_local_index, this->mNodes[new_node_2_global_index]);
02220 this->GetElement(elementIndex)->AddNode(node_A_local_index, pNode);
02221 this->GetElement(elementIndex)->AddNode(node_A_local_index, this->mNodes[new_node_1_global_index]);
02222
02223
02224 if (next_node_1 == previous_node_2)
02225 {
02226 p_element_intersection_1->AddNode((local_index_1 + p_element_intersection_1->GetNumNodes() - 1)%(p_element_intersection_1->GetNumNodes()), this->mNodes[new_node_2_global_index]);
02227 p_element_intersection_2->AddNode(local_index_2, this->mNodes[new_node_1_global_index]);
02228 }
02229 else
02230 {
02231 assert(next_node_2 == previous_node_1);
02232
02233 p_element_intersection_1->AddNode(local_index_1, this->mNodes[new_node_1_global_index]);
02234 p_element_intersection_2->AddNode((local_index_2 + p_element_intersection_2->GetNumNodes() - 1)%(p_element_intersection_2->GetNumNodes()), this->mNodes[new_node_2_global_index]);
02235 }
02236
02237
02238 assert(pNode->GetNumContainingElements() == 3);
02239 assert(this->mNodes[new_node_1_global_index]->GetNumContainingElements() == 2);
02240 assert(this->mNodes[new_node_2_global_index]->GetNumContainingElements() == 2);
02241 }
02242 }
02243 }
02244 else
02245 {
02246 EXCEPTION("Trying to merge a node, contained in more than 2 elements, into another element, this is not possible with the vertex mesh.");
02247 }
02248 }
02249
02250 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
02251 void MutableVertexMesh<ELEMENT_DIM, SPACE_DIM>::PerformVoidRemoval(Node<SPACE_DIM>* pNodeA, Node<SPACE_DIM>* pNodeB, Node<SPACE_DIM>* pNodeC)
02252 {
02253 unsigned nodeA_index = pNodeA->GetIndex();
02254 unsigned nodeB_index = pNodeB->GetIndex();
02255 unsigned nodeC_index = pNodeC->GetIndex();
02256
02257 c_vector<double, SPACE_DIM> nodes_midpoint = pNodeA->rGetLocation() + this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeB->rGetLocation())/3.0
02258 + this->GetVectorFromAtoB(pNodeA->rGetLocation(), pNodeC->rGetLocation())/3.0;
02259
02260 Node<SPACE_DIM>* p_low_node_A_B = (nodeA_index < nodeB_index) ? pNodeA : pNodeB;
02261 Node<SPACE_DIM>* p_low_node = (p_low_node_A_B->GetIndex() < nodeC_index) ? p_low_node_A_B : pNodeC;
02262 PerformNodeMerge(pNodeA,pNodeB);
02263 PerformNodeMerge(p_low_node_A_B,pNodeC);
02264
02265 c_vector<double, SPACE_DIM>& r_low_node_location = p_low_node->rGetModifiableLocation();
02266
02267 r_low_node_location = nodes_midpoint;
02268
02269
02270 p_low_node->SetAsBoundaryNode(false);
02271
02272
02273 RemoveDeletedNodes();
02274 }
02275
02276
02278
02280
02281 template class MutableVertexMesh<1,1>;
02282 template class MutableVertexMesh<1,2>;
02283 template class MutableVertexMesh<1,3>;
02284 template class MutableVertexMesh<2,2>;
02285 template class MutableVertexMesh<2,3>;
02286 template class MutableVertexMesh<3,3>;
02287
02288
02289
02290 #include "SerializationExportWrapperForCpp.hpp"
02291 EXPORT_TEMPLATE_CLASS_ALL_DIMS(MutableVertexMesh);