MutableMesh.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include <map>
00030 #include <cstring>
00031 
00032 #include "MutableMesh.hpp"
00033 #include "OutputFileHandler.hpp"
00034 #include "TrianglesMeshReader.hpp"
00035 
00036 //Jonathan Shewchuk's triangle and Hang Si's tetgen
00037 #define REAL double
00038 #define VOID void
00039 #include "triangle.h"
00040 #include "tetgen.h"
00041 #undef REAL
00042 #undef VOID
00043 
00044 
00045 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00046 MutableMesh<ELEMENT_DIM, SPACE_DIM>::MutableMesh()
00047     : mAddedNodes(false)
00048 {
00049     this->mMeshChangesDuringSimulation = true;
00050 }
00051 
00052 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00053 MutableMesh<ELEMENT_DIM, SPACE_DIM>::MutableMesh(std::vector<Node<SPACE_DIM> *> nodes)
00054 {
00055     this->mMeshChangesDuringSimulation = true;
00056     Clear();
00057     for (unsigned index=0; index<nodes.size(); index++)
00058     {
00059         Node<SPACE_DIM>* p_temp_node = nodes[index];
00060         this->mNodes.push_back(p_temp_node);
00061     }
00062     mAddedNodes = true;
00063     NodeMap node_map(nodes.size());
00064     ReMesh(node_map);
00065 }
00066 
00067 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00068 MutableMesh<ELEMENT_DIM, SPACE_DIM>::~MutableMesh()
00069 {
00070     Clear();
00071 }
00072 
00073 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00074 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::AddNode(Node<SPACE_DIM>* pNewNode)
00075 {
00076     if (mDeletedNodeIndices.empty())
00077     {
00078         pNewNode->SetIndex(this->mNodes.size());
00079         this->mNodes.push_back(pNewNode);
00080     }
00081     else
00082     {
00083         unsigned index = mDeletedNodeIndices.back();
00084         pNewNode->SetIndex(index);
00085         mDeletedNodeIndices.pop_back();
00086         delete this->mNodes[index];
00087         this->mNodes[index] = pNewNode;
00088     }
00089     mAddedNodes = true;
00090     return pNewNode->GetIndex();
00091 }
00092 
00093 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00094 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::Clear()
00095 {
00096     mDeletedElementIndices.clear();
00097     mDeletedBoundaryElementIndices.clear();
00098     mDeletedNodeIndices.clear();
00099     mAddedNodes = false;
00100 
00101     TetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::Clear();
00102 }
00103 
00104 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00105 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::GetNumBoundaryElements() const
00106 {
00107     return this->mBoundaryElements.size() - mDeletedBoundaryElementIndices.size();
00108 }
00109 
00110 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00111 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::GetNumElements() const
00112 {
00113     return this->mElements.size() - mDeletedElementIndices.size();
00114 }
00115 
00116 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00117 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::GetNumNodes() const
00118 {
00119     return this->mNodes.size() - mDeletedNodeIndices.size();
00120 }
00121 
00128 template<>
00129 void MutableMesh<1, 1>::RescaleMeshFromBoundaryNode(ChastePoint<1> updatedPoint, unsigned boundaryNodeIndex)
00130 {
00131     assert(this->GetNode(boundaryNodeIndex)->IsBoundaryNode());
00132     double scaleFactor = updatedPoint[0] / this->GetNode(boundaryNodeIndex)->GetPoint()[0];
00133     double temp;
00134     for (unsigned i=0; i < boundaryNodeIndex+1; i++)
00135     {
00136         temp = scaleFactor * this->mNodes[i]->GetPoint()[0];
00137         ChastePoint<1> newPoint(temp);
00138         this->mNodes[i]->SetPoint(newPoint);
00139     }
00140     this->RefreshMesh();
00141 }
00142 
00143 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00144 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::SetNode(unsigned index,
00145         ChastePoint<SPACE_DIM> point,
00146         bool concreteMove)
00147 {
00148     this->mNodes[index]->SetPoint(point);
00149 
00150     if (concreteMove)
00151     {
00152         for (typename Node<SPACE_DIM>::ContainingElementIterator it = this->mNodes[index]->ContainingElementsBegin();
00153              it != this->mNodes[index]->ContainingElementsEnd();
00154              ++it)
00155         {
00156             if (ELEMENT_DIM == SPACE_DIM)
00157             {
00158                 try
00159                 {
00160                     GetElement(*it)->CalculateInverseJacobian(this->mElementJacobians[ (*it) ],
00161                                                               this->mElementJacobianDeterminants[ (*it) ],
00162                                                               this->mElementInverseJacobians[ (*it) ]);
00163                 }
00164                 catch (Exception e)
00165                 {
00166                         EXCEPTION("Moving node caused an element to have a non-positive Jacobian determinant");
00167                 }
00168             }
00169             else
00170             {
00171                 c_vector<double,SPACE_DIM> previous_direction = this->mElementWeightedDirections[ (*it) ];
00172 
00173                 GetElement(*it)->CalculateWeightedDirection(this->mElementWeightedDirections[ (*it) ],
00174                                                             this->mElementJacobianDeterminants[ (*it) ]);
00175 
00176                 if ( inner_prod(previous_direction, this->mElementWeightedDirections[ (*it) ]) < 0)
00177                 {
00178                     EXCEPTION("Moving node caused an subspace element to change direction");
00179                 }
00180 
00181             }
00182         }
00183         for (typename Node<SPACE_DIM>::ContainingBoundaryElementIterator it = this->mNodes[index]->ContainingBoundaryElementsBegin();
00184              it != this->mNodes[index]->ContainingBoundaryElementsEnd();
00185              ++it)
00186         {
00187             try
00188             {
00189                 GetBoundaryElement(*it)->CalculateWeightedDirection(this->mBoundaryElementWeightedDirections[ (*it) ],
00190                                                                     this->mBoundaryElementJacobianDeterminants[ (*it) ]);
00191             }
00192             catch (Exception e)
00193             {
00194                 EXCEPTION("Moving node caused a boundary element to have a non-positive Jacobian determinant");
00195             }
00196         }
00197     }
00198 }
00199 
00200 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00201 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::DeleteNode(unsigned index)
00202 {
00203     if (this->mNodes[index]->IsDeleted())
00204     {
00205         EXCEPTION("Trying to delete a deleted node");
00206     }
00207     unsigned target_index = (unsigned)(-1);
00208     bool found_target=false;
00209     for (typename Node<SPACE_DIM>::ContainingElementIterator it = this->mNodes[index]->ContainingElementsBegin();
00210          !found_target && it != this->mNodes[index]->ContainingElementsEnd();
00211          ++it)
00212     {
00213         Element <ELEMENT_DIM,SPACE_DIM>* p_element = GetElement(*it);
00214         for (unsigned i=0; i<=ELEMENT_DIM && !found_target; i++)
00215         {
00216             target_index = p_element->GetNodeGlobalIndex(i);
00217             try
00218             {
00219                 MoveMergeNode(index, target_index, false);
00220                 found_target = true;
00221             }
00222             catch (Exception e)
00223             {
00224                 // Just try the next node
00225             }
00226         }
00227     }
00228     if (!found_target)
00229     {
00230         EXCEPTION("Failure to delete node");
00231     }
00232 
00233     MoveMergeNode(index, target_index);
00234 }
00235 
00236 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00237 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::DeleteNodePriorToReMesh(unsigned index)
00238 {
00239     this->mNodes[index]->MarkAsDeleted();
00240     mDeletedNodeIndices.push_back(index);
00241 }
00242 
00243 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00244 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::MoveMergeNode(unsigned index,
00245         unsigned targetIndex,
00246         bool concreteMove)
00247 {
00248     if (this->mNodes[index]->IsDeleted())
00249     {
00250         EXCEPTION("Trying to move a deleted node");
00251     }
00252 
00253     if (index == targetIndex)
00254     {
00255         EXCEPTION("Trying to merge a node with itself");
00256     }
00257     if (this->mNodes[index]->IsBoundaryNode())
00258     {
00259         if (!this->mNodes[targetIndex]->IsBoundaryNode())
00260         {
00261             EXCEPTION("A boundary node can only be moved on to another boundary node");
00262         }
00263     }
00264     std::set<unsigned> unshared_element_indices;
00265     std::set_difference(this->mNodes[index]->rGetContainingElementIndices().begin(),
00266                         this->mNodes[index]->rGetContainingElementIndices().end(),
00267                         this->mNodes[targetIndex]->rGetContainingElementIndices().begin(),
00268                         this->mNodes[targetIndex]->rGetContainingElementIndices().end(),
00269                         std::inserter(unshared_element_indices, unshared_element_indices.begin()));
00270 
00271 
00272     if (unshared_element_indices.size() == this->mNodes[index]->rGetContainingElementIndices().size())
00273     {
00274         EXCEPTION("These nodes cannot be merged since they are not neighbours");
00275     }
00276 
00277     std::set<unsigned> unshared_boundary_element_indices;
00278     std::set_difference(this->mNodes[index]->rGetContainingBoundaryElementIndices().begin(),
00279                         this->mNodes[index]->rGetContainingBoundaryElementIndices().end(),
00280                         this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().begin(),
00281                         this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().end(),
00282                         std::inserter(unshared_boundary_element_indices, unshared_boundary_element_indices.begin()));
00283 
00284 
00285     if (this->mNodes[index]->IsBoundaryNode())
00286     {
00287         if (unshared_boundary_element_indices.size()
00288             == this->mNodes[index]->rGetContainingBoundaryElementIndices().size())
00289         {
00290             //May be redundant (only thrown in 1D tests)
00291             EXCEPTION("These nodes cannot be merged since they are not neighbours on the boundary");
00292         }
00293     }
00294 
00295     this->mNodes[index]->rGetModifiableLocation() = this->mNodes[targetIndex]->rGetLocation();
00296 
00297     for (std::set<unsigned>::const_iterator element_iter=unshared_element_indices.begin();
00298              element_iter != unshared_element_indices.end();
00299              element_iter++)
00300     {
00301         try
00302         {
00303             if (SPACE_DIM == ELEMENT_DIM)
00304             {
00305                 this->GetElement(*element_iter)->CalculateInverseJacobian(this->mElementJacobians[(*element_iter)],
00306                                                                           this->mElementJacobianDeterminants[(*element_iter)],
00307                                                                           this->mElementInverseJacobians[ (*element_iter) ]);
00308             }
00309             else
00310             {
00311                 this->GetElement(*element_iter)->CalculateWeightedDirection(this->mElementWeightedDirections[(*element_iter)],
00312                                                                             this->mElementJacobianDeterminants[(*element_iter)]);
00313             }
00314 
00315             if (concreteMove)
00316             {
00317                 this->GetElement(*element_iter)->ReplaceNode(this->mNodes[index], this->mNodes[targetIndex]);
00318             }
00319 
00320         }
00321         catch (Exception e)
00322         {
00323             EXCEPTION("Moving node caused an element to have a non-positive Jacobian determinant");
00324         }
00325     }
00326 
00327     for (std::set<unsigned>::const_iterator boundary_element_iter=
00328                  unshared_boundary_element_indices.begin();
00329              boundary_element_iter != unshared_boundary_element_indices.end();
00330              boundary_element_iter++)
00331     {
00332 
00333         this->GetBoundaryElement(*boundary_element_iter)->CalculateWeightedDirection(this->mBoundaryElementWeightedDirections[(*boundary_element_iter)],
00334                                                                                      this->mBoundaryElementJacobianDeterminants[(*boundary_element_iter)]);
00335 
00336         if (concreteMove)
00337         {
00338             this->GetBoundaryElement(*boundary_element_iter)->ReplaceNode(this->mNodes[index], this->mNodes[targetIndex]);
00339         }
00340     }
00341 
00342     std::set<unsigned> shared_element_indices;
00343     std::set_intersection(this->mNodes[index]->rGetContainingElementIndices().begin(),
00344                           this->mNodes[index]->rGetContainingElementIndices().end(),
00345                           this->mNodes[targetIndex]->rGetContainingElementIndices().begin(),
00346                           this->mNodes[targetIndex]->rGetContainingElementIndices().end(),
00347                           std::inserter(shared_element_indices, shared_element_indices.begin()));
00348 
00349     for (std::set<unsigned>::const_iterator element_iter=shared_element_indices.begin();
00350              element_iter != shared_element_indices.end();
00351              element_iter++)
00352     {
00353         if (concreteMove)
00354         {
00355             this->GetElement(*element_iter)->MarkAsDeleted();
00356             this->mElementJacobianDeterminants[ (*element_iter) ] = 0.0; //This used to be done in MarkAsDeleted
00357             mDeletedElementIndices.push_back(*element_iter);
00358         }
00359         else
00360         {
00361             this->mElementJacobianDeterminants[ (*element_iter) ] = 0.0;
00362         }
00363     }
00364 
00365 
00366     std::set<unsigned> shared_boundary_element_indices;
00367     std::set_intersection(this->mNodes[index]->rGetContainingBoundaryElementIndices().begin(),
00368                           this->mNodes[index]->rGetContainingBoundaryElementIndices().end(),
00369                           this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().begin(),
00370                           this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().end(),
00371                           std::inserter(shared_boundary_element_indices, shared_boundary_element_indices.begin()));
00372 
00373     for (std::set<unsigned>::const_iterator boundary_element_iter=shared_boundary_element_indices.begin();
00374              boundary_element_iter != shared_boundary_element_indices.end();
00375              boundary_element_iter++)
00376     {
00377         if (concreteMove)
00378         {
00379             this->GetBoundaryElement(*boundary_element_iter)->MarkAsDeleted();
00380             this->mBoundaryElementJacobianDeterminants[ (*boundary_element_iter) ] = 0.0; //This used to be done in MarkAsDeleted
00381             mDeletedBoundaryElementIndices.push_back(*boundary_element_iter);
00382         }
00383         else
00384         {
00385             this->mBoundaryElementJacobianDeterminants[ (*boundary_element_iter) ] = 0.0;
00386             this->mBoundaryElementWeightedDirections[ (*boundary_element_iter) ] = zero_vector<double>(SPACE_DIM);
00387         }
00388     }
00389 
00390     if (concreteMove)
00391     {
00392         this->mNodes[index]->MarkAsDeleted();
00393         mDeletedNodeIndices.push_back(index);
00394     }
00395 }
00396 
00397 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00398 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::RefineElement(
00399     Element<ELEMENT_DIM,SPACE_DIM>* pElement,
00400     ChastePoint<SPACE_DIM> point)
00401 {
00402     //Check that the point is in the element
00403     if (pElement->IncludesPoint(point, true) == false)
00404     {
00405         EXCEPTION("RefineElement could not be started (point is not in element)");
00406     }
00407 
00408     // Add a new node from the point that is passed to RefineElement
00409     unsigned new_node_index = AddNode(new Node<SPACE_DIM>(0, point.rGetLocation()));
00410     // Note: the first argument is the index of the node, which is going to be
00411     //       overriden by AddNode, so it can safely be ignored
00412 
00413     // This loop constructs the extra elements which are going to fill the space
00414     for (unsigned i = 0; i < ELEMENT_DIM; i++)
00415     {
00416 
00417         // First, make a copy of the current element making sure we update its index
00418         unsigned new_elt_index;
00419         if (mDeletedElementIndices.empty())
00420         {
00421             new_elt_index = this->mElements.size();
00422         }
00423         else
00424         {
00425             new_elt_index = mDeletedElementIndices.back();
00426             mDeletedElementIndices.pop_back();
00427         }
00428 
00429         Element<ELEMENT_DIM,SPACE_DIM>* p_new_element=
00430             new Element<ELEMENT_DIM,SPACE_DIM>(*pElement, new_elt_index);
00431 
00432         // Second, update the node in the element with the new one
00433         p_new_element->UpdateNode(ELEMENT_DIM-1-i, this->mNodes[new_node_index]);
00434 
00435         // Third, add the new element to the set
00436         if ((unsigned) new_elt_index == this->mElements.size())
00437         {
00438             this->mElements.push_back(p_new_element);
00439         }
00440         else
00441         {
00442             delete this->mElements[new_elt_index];
00443             this->mElements[new_elt_index] = p_new_element;
00444         }
00445     }
00446 
00447     // Lastly, update the last node in the element to be refined
00448     pElement->UpdateNode(ELEMENT_DIM, this->mNodes[new_node_index]);
00449 
00450     return new_node_index;
00451 }
00452 
00453 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00454 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::DeleteBoundaryNodeAt(unsigned index)
00455 {
00456     if (!this->mNodes[index]->IsBoundaryNode() )
00457     {
00458         EXCEPTION(" You may only delete a boundary node ");
00459     }
00460 
00461     this->mNodes[index]->MarkAsDeleted();
00462     mDeletedNodeIndices.push_back(index);
00463 
00464     // Update the boundary node vector
00465     typename std::vector<Node<SPACE_DIM>*>::iterator b_node_iter
00466     = std::find(this->mBoundaryNodes.begin(), this->mBoundaryNodes.end(), this->mNodes[index]);
00467     this->mBoundaryNodes.erase(b_node_iter);
00468 
00469     // Remove boundary elements containing this node
00470     std::set<unsigned> boundary_element_indices = this->mNodes[index]->rGetContainingBoundaryElementIndices();
00471     std::set<unsigned>::const_iterator boundary_element_indices_iterator = boundary_element_indices.begin();
00472     while (boundary_element_indices_iterator != boundary_element_indices.end())
00473     {
00474         BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_boundary_element = this->GetBoundaryElement(*boundary_element_indices_iterator);
00475         p_boundary_element->MarkAsDeleted();
00476         mDeletedBoundaryElementIndices.push_back(*boundary_element_indices_iterator);
00477         boundary_element_indices_iterator++;
00478     }
00479 
00480     // Remove elements containing this node
00481     std::set<unsigned> element_indices = this->mNodes[index]->rGetContainingElementIndices();
00482     std::set<unsigned>::const_iterator element_indices_iterator = element_indices.begin();
00483     while (element_indices_iterator != element_indices.end())
00484     {
00485         Element<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(*element_indices_iterator);
00486         for (unsigned i=0; i<p_element->GetNumNodes(); i++)
00487         {
00488             Node<SPACE_DIM>* p_node = p_element->GetNode(i);
00489             if (!p_node->IsDeleted())
00490             {
00491                 p_node->SetAsBoundaryNode();
00492                 // Update the boundary node vector
00493                 this->mBoundaryNodes.push_back(p_node);
00494             }
00495         }
00496         p_element->MarkAsDeleted();
00497         mDeletedElementIndices.push_back(p_element->GetIndex());
00498         element_indices_iterator++;
00499     }
00500 }
00501 
00502 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00503 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::ReIndex(NodeMap& map)
00504 {
00505     assert(!mAddedNodes);
00506     map.Resize(this->GetNumAllNodes());
00507 
00508     std::vector<Element<ELEMENT_DIM, SPACE_DIM> *> live_elements;
00509 
00510     for (unsigned i=0; i<this->mElements.size(); i++)
00511     {
00512         assert(i==this->mElements[i]->GetIndex()); // We need this to be true to be able to reindex the jacobian cache
00513         if (this->mElements[i]->IsDeleted())
00514         {
00515             delete this->mElements[i];
00516         }
00517         else
00518         {
00519             live_elements.push_back(this->mElements[i]);
00520 
00521             unsigned this_element_index = this->mElements[i]->GetIndex();
00522             if (SPACE_DIM == ELEMENT_DIM)
00523             {
00524                 this->mElementJacobians[live_elements.size()-1] = this->mElementJacobians[this_element_index];
00525                 this->mElementInverseJacobians[live_elements.size()-1] = this->mElementInverseJacobians[this_element_index];
00526             }
00527             else
00528             {
00529                 this->mElementWeightedDirections[live_elements.size()-1] = this->mElementWeightedDirections[this_element_index];
00530             }
00531             this->mElementJacobianDeterminants[live_elements.size()-1] = this->mElementJacobianDeterminants[this_element_index];
00532         }
00533     }
00534 
00535     assert(mDeletedElementIndices.size() == this->mElements.size()-live_elements.size());
00536     mDeletedElementIndices.clear();
00537     this->mElements = live_elements;
00538     unsigned num_elements = this->mElements.size();
00539 
00540     if (SPACE_DIM == ELEMENT_DIM)
00541     {
00542         this->mElementJacobians.resize(num_elements);
00543         this->mElementInverseJacobians.resize(num_elements);
00544     }
00545     else
00546     {
00547         this->mElementWeightedDirections.resize(num_elements);
00548     }
00549     this->mElementJacobianDeterminants.resize(num_elements);
00550 
00551     std::vector<Node<SPACE_DIM> *> live_nodes;
00552     for (unsigned i=0; i<this->mNodes.size(); i++)
00553     {
00554         if (this->mNodes[i]->IsDeleted())
00555         {
00556             delete this->mNodes[i];
00557             map.SetDeleted(i);
00558         }
00559         else
00560         {
00561             live_nodes.push_back(this->mNodes[i]);
00562             // the nodes will have their index set to be the index into the live_nodes
00563             // vector further down
00564             map.SetNewIndex(i, (unsigned)(live_nodes.size()-1));
00565         }
00566     }
00567 
00568     assert(mDeletedNodeIndices.size() == this->mNodes.size()-live_nodes.size());
00569     this->mNodes = live_nodes;
00570     mDeletedNodeIndices.clear();
00571 
00572     std::vector<BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *> live_boundary_elements;
00573     for (unsigned i=0; i<this->mBoundaryElements.size(); i++)
00574     {
00575         if (this->mBoundaryElements[i]->IsDeleted())
00576         {
00577             delete this->mBoundaryElements[i];
00578         }
00579         else
00580         {
00581             live_boundary_elements.push_back(this->mBoundaryElements[i]);
00582 
00583             this->mBoundaryElementWeightedDirections[live_boundary_elements.size()-1] = this->mBoundaryElementWeightedDirections[this->mBoundaryElements[i]->GetIndex()];
00584             this->mBoundaryElementJacobianDeterminants[live_boundary_elements.size()-1] = this->mBoundaryElementJacobianDeterminants[this->mBoundaryElements[i]->GetIndex()];
00585         }
00586     }
00587 
00588     assert(mDeletedBoundaryElementIndices.size() == this->mBoundaryElements.size()-live_boundary_elements.size());
00589     this->mBoundaryElements = live_boundary_elements;
00590     mDeletedBoundaryElementIndices.clear();
00591 
00592     unsigned num_boundary_elements = this->mBoundaryElements.size();
00593 
00594     this->mBoundaryElementWeightedDirections.resize(num_boundary_elements);
00595     this->mBoundaryElementJacobianDeterminants.resize(num_boundary_elements);
00596 
00597     for (unsigned i=0; i<this->mNodes.size(); i++)
00598     {
00599         this->mNodes[i]->SetIndex(i);
00600     }
00601     for (unsigned i=0; i<this->mElements.size(); i++)
00602     {
00603 
00604         this->mElements[i]->ResetIndex(i);
00605     }
00606     for (unsigned i=0; i<this->mBoundaryElements.size(); i++)
00607     {
00608         this->mBoundaryElements[i]->ResetIndex(i);
00609     }
00610 }
00611 
00612 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00613 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::ReMesh(NodeMap& map)
00614 {
00615     // Make sure that we are in the correct dimension - this code will be eliminated at compile time
00616     #define COVERAGE_IGNORE
00617     assert( ELEMENT_DIM == SPACE_DIM );
00618     #undef COVERAGE_IGNORE
00619 
00620     // Avoid some triangle/tetgen errors: need at least four
00621     // nodes for tetgen, and at least three for triangle
00622     if (GetNumNodes() <= SPACE_DIM)
00623     {
00624         EXCEPTION("The number of nodes must exceed the spatial dimension.");
00625     }
00626 
00627     // Make sure the map is big enough
00628     map.Resize(this->GetNumAllNodes());
00629     if (mAddedNodes || !mDeletedNodeIndices.empty())
00630     {
00631         //Size of mesh is about to change
00632         if (this->mpDistributedVectorFactory)
00633         {
00634             delete this->mpDistributedVectorFactory;
00635             this->mpDistributedVectorFactory = new DistributedVectorFactory(this->GetNumNodes());
00636         }
00637     }
00638     if (SPACE_DIM==1)
00639     {
00640         // Store the node locations
00641         std::vector<c_vector<double, SPACE_DIM> > old_node_locations;
00642         unsigned new_index = 0;
00643         for (unsigned i=0; i<this->GetNumAllNodes(); i++)
00644         {
00645             if (this->mNodes[i]->IsDeleted())
00646             {
00647                 map.SetDeleted(i);
00648             }
00649             else
00650             {
00651                 map.SetNewIndex(i, new_index);
00652                 old_node_locations.push_back(this->mNodes[i]->rGetLocation());
00653                 new_index++;
00654             }
00655         }
00656 
00657         // Remove current data
00658         Clear();
00659 
00660         // Construct the nodes and boundary nodes
00661         for (unsigned node_index=0; node_index<old_node_locations.size(); node_index++)
00662         {
00663             Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index, old_node_locations[node_index], false);
00664             this->mNodes.push_back(p_node);
00665 
00666             // As we're in 1D, the boundary nodes are simply at either end of the mesh
00667             if ( node_index==0 || node_index==old_node_locations.size()-1 )
00668             {
00669                 this->mBoundaryNodes.push_back(p_node);
00670             }
00671         }
00672 
00673         // Create a map between node indices and node locations
00674         std::map<double, unsigned> location_index_map;
00675         for (unsigned i=0; i<this->mNodes.size(); i++)
00676         {
00677             location_index_map[this->mNodes[i]->rGetLocation()[0]] = this->mNodes[i]->GetIndex();
00678         }
00679 
00680         // Use this map to generate a vector of node indices that are ordered spatially
00681         std::vector<unsigned> node_indices_ordered_spatially;
00682         for (std::map<double, unsigned>::iterator iter = location_index_map.begin();
00683              iter != location_index_map.end();
00684              ++iter)
00685         {
00686             node_indices_ordered_spatially.push_back(iter->second);
00687         }
00688 
00689         // Construct the elements
00690         this->mElements.reserve(old_node_locations.size()-1);
00691         for (unsigned element_index=0; element_index<old_node_locations.size()-1; element_index++)
00692         {
00693             std::vector<Node<SPACE_DIM>*> nodes;
00694             for (unsigned j=0; j<2; j++)
00695             {
00696                 unsigned global_node_index = node_indices_ordered_spatially[element_index + j];
00697                 assert(global_node_index < this->mNodes.size());
00698                 nodes.push_back(this->mNodes[global_node_index]);
00699             }
00700             this->mElements.push_back(new Element<ELEMENT_DIM, SPACE_DIM>(element_index, nodes));
00701         }
00702 
00703         // Construct the two boundary elements - as we're in 1D, these are simply at either end of the mesh
00704         std::vector<Node<SPACE_DIM>*> nodes;
00705         nodes.push_back(this->mNodes[0]);
00706         this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(0, nodes));
00707 
00708         nodes.clear();
00709         nodes.push_back(this->mNodes[old_node_locations.size()-1]);
00710         this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(1, nodes));
00711 
00712         this->RefreshJacobianCachedData();
00713     }
00714     else if (SPACE_DIM==2)  // In 2D, remesh using triangle via library calls
00715     {
00716         struct triangulateio mesher_input, mesher_output;
00717         this->InitialiseTriangulateIo(mesher_input);
00718         this->InitialiseTriangulateIo(mesher_output);
00719 
00720         this->ExportToMesher(map, mesher_input);
00721 
00722         // Library call
00723         triangulate((char*)"Qze", &mesher_input, &mesher_output, NULL);
00724         
00725         this->ImportFromMesher(mesher_output, mesher_output.numberoftriangles, mesher_output.trianglelist, mesher_output.numberofedges, mesher_output.edgelist, mesher_output.edgemarkerlist);
00726 
00727         //Tidy up triangle
00728         this->FreeTriangulateIo(mesher_input);
00729         this->FreeTriangulateIo(mesher_output);
00730     }
00731     else // in 3D, remesh using tetgen
00732     {
00733 
00734         struct tetgen::tetgenio mesher_input, mesher_output;
00735 
00736         this->ExportToMesher(map, mesher_input);
00737 
00738         // Library call
00739         tetgen::tetrahedralize((char*)"Qz", &mesher_input, &mesher_output);
00740         
00741         this->ImportFromMesher(mesher_output, mesher_output.numberoftetrahedra, mesher_output.tetrahedronlist, mesher_output.numberoftrifaces, mesher_output.trifacelist, NULL);
00742     }
00743 }
00744 
00745 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00746 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::ReMesh()
00747 {
00748     NodeMap map(GetNumNodes());
00749     ReMesh(map);
00750 }
00751 
00752 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00753 bool MutableMesh<ELEMENT_DIM, SPACE_DIM>::CheckIsVoronoi(Element<ELEMENT_DIM, SPACE_DIM>* pElement, double maxPenetration)
00754 {
00755     assert(ELEMENT_DIM == SPACE_DIM);
00756     unsigned num_nodes = pElement->GetNumNodes();
00757     std::set<unsigned> neighbouring_elements_indices;
00758     std::set< Element<ELEMENT_DIM,SPACE_DIM> *> neighbouring_elements;
00759     std::set<unsigned> neighbouring_nodes_indices;
00760 
00761     // Form a set of neighbouring elements via the nodes
00762     for (unsigned i=0; i<num_nodes; i++)
00763     {
00764         Node<SPACE_DIM>* p_node = pElement->GetNode(i);
00765         neighbouring_elements_indices = p_node->rGetContainingElementIndices();
00767         for (std::set<unsigned>::const_iterator it = neighbouring_elements_indices.begin();
00768              it != neighbouring_elements_indices.end();
00769              ++it)
00770         {
00771             neighbouring_elements.insert(this->GetElement(*it));
00772         }
00773     }
00774     neighbouring_elements.erase(pElement);
00775 
00776     // For each neighbouring element find the supporting nodes
00777     typedef typename std::set<Element<ELEMENT_DIM,SPACE_DIM> *>::const_iterator ElementIterator;
00778 
00779     for (ElementIterator it = neighbouring_elements.begin();
00780          it != neighbouring_elements.end();
00781          ++it)
00782     {
00783         for (unsigned i=0; i<num_nodes; i++)
00784         {
00785             neighbouring_nodes_indices.insert((*it)->GetNodeGlobalIndex(i));
00786         }
00787     }
00788 
00789     // Remove the nodes that support this element
00790     for (unsigned i = 0; i < num_nodes; i++)
00791     {
00792         neighbouring_nodes_indices.erase(pElement->GetNodeGlobalIndex(i));
00793     }
00794 
00795     // Get the circumsphere information
00796     c_vector<double, SPACE_DIM+1> this_circum_centre;
00797 
00798     this_circum_centre = pElement->CalculateCircumsphere(this->mElementJacobians[pElement->GetIndex()], this->mElementInverseJacobians[pElement->GetIndex()]);
00799 
00800     // Copy the actualy circumcentre into a smaller vector
00801     c_vector<double, ELEMENT_DIM> circum_centre;
00802     for (unsigned i=0; i<ELEMENT_DIM; i++)
00803     {
00804         circum_centre[i] = this_circum_centre[i];
00805     }
00806 
00807     for (std::set<unsigned>::const_iterator it = neighbouring_nodes_indices.begin();
00808          it != neighbouring_nodes_indices.end();
00809          ++it)
00810     {
00811         c_vector<double, ELEMENT_DIM> node_location = this->GetNode(*it)->rGetLocation();
00812 
00813         // Calculate vector from circumcenter to node
00814         node_location -= circum_centre;
00815 
00816         // This is to calculate the squared distance betweeen them
00817         double squared_distance = inner_prod(node_location, node_location);
00818 
00819         // If the squared idstance is less than the elements circum-radius(squared),
00820         // then the Voronoi property is violated.
00821         if (squared_distance < this_circum_centre[ELEMENT_DIM])
00822         {
00823             // We know the node is inside the circumsphere, but we don't know how far
00824             double radius = sqrt(this_circum_centre[ELEMENT_DIM]);
00825             double distance = radius - sqrt(squared_distance);
00826 
00827             // If the node penetration is greater than supplied maximum penetration factor
00828             if (distance/radius > maxPenetration)
00829             {
00830                 return false;
00831             }
00832         }
00833     }
00834     return true;
00835 }
00836 
00837 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00838 bool MutableMesh<ELEMENT_DIM, SPACE_DIM>::CheckIsVoronoi(double maxPenetration)
00839 {
00840     // Looping through all the elements in the mesh
00842     for (unsigned i=0; i<this->mElements.size(); i++)
00843     {
00844         // Check if the element is not deleted
00845         if (!this->mElements[i]->IsDeleted())
00846         {
00847             // Checking the Voronoi of the Element
00848             if (CheckIsVoronoi(this->mElements[i], maxPenetration) == false)
00849             {
00850                 return false;
00851             }
00852         }
00853     }
00854     return true;
00855 }
00856 
00857 
00859 // Explicit instantiation
00861 
00862 template class MutableMesh<1,1>;
00863 template class MutableMesh<1,2>;
00864 template class MutableMesh<1,3>;
00865 template class MutableMesh<2,2>;
00866 template class MutableMesh<2,3>;
00867 template class MutableMesh<3,3>;
00868 
00869 
00870 // Serialization for Boost >= 1.36
00871 #include "SerializationExportWrapperForCpp.hpp"
00872 EXPORT_TEMPLATE_CLASS_ALL_DIMS(MutableMesh)

Generated on Mon Nov 1 12:35:23 2010 for Chaste by  doxygen 1.5.5