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
00037 #define REAL double
00038 #define VOID void
00039 #include "triangle.h"
00040 #undef REAL
00041 #undef VOID
00042 
00043 
00044 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00045 MutableMesh<ELEMENT_DIM, SPACE_DIM>::MutableMesh()
00046     : mAddedNodes(false)
00047 {
00048     this->mMeshChangesDuringSimulation = true;
00049 }
00050 
00051 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00052 MutableMesh<ELEMENT_DIM, SPACE_DIM>::MutableMesh(std::vector<Node<SPACE_DIM> *> nodes)
00053 {
00054     this->mMeshChangesDuringSimulation = true;
00055     Clear();
00056     for (unsigned index=0; index<nodes.size(); index++)
00057     {
00058         Node<SPACE_DIM>* p_temp_node = nodes[index];
00059         this->mNodes.push_back(p_temp_node);
00060     }
00061     mAddedNodes = true;
00062     NodeMap node_map(nodes.size());
00063     ReMesh(node_map);
00064 }
00065 
00066 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00067 MutableMesh<ELEMENT_DIM, SPACE_DIM>::~MutableMesh()
00068 {
00069     Clear();
00070 }
00071 
00072 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00073 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::AddNode(Node<SPACE_DIM>* pNewNode)
00074 {
00075     if (mDeletedNodeIndices.empty())
00076     {
00077         pNewNode->SetIndex(this->mNodes.size());
00078         this->mNodes.push_back(pNewNode);
00079     }
00080     else
00081     {
00082         unsigned index = mDeletedNodeIndices.back();
00083         pNewNode->SetIndex(index);
00084         mDeletedNodeIndices.pop_back();
00085         delete this->mNodes[index];
00086         this->mNodes[index] = pNewNode;
00087     }
00088     mAddedNodes = true;
00089     return pNewNode->GetIndex();
00090 }
00091 
00092 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00093 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::Clear()
00094 {
00095     mDeletedElementIndices.clear();
00096     mDeletedBoundaryElementIndices.clear();
00097     mDeletedNodeIndices.clear();
00098     mAddedNodes = false;
00099 
00100     TetrahedralMesh<ELEMENT_DIM, SPACE_DIM>::Clear();
00101 }
00102 
00103 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00104 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::GetNumBoundaryElements() const
00105 {
00106     return this->mBoundaryElements.size() - mDeletedBoundaryElementIndices.size();
00107 }
00108 
00109 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00110 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::GetNumElements() const
00111 {
00112     return this->mElements.size() - mDeletedElementIndices.size();
00113 }
00114 
00115 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00116 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::GetNumNodes() const
00117 {
00118     return this->mNodes.size() - mDeletedNodeIndices.size();
00119 }
00120 
00127 template<>
00128 void MutableMesh<1, 1>::RescaleMeshFromBoundaryNode(ChastePoint<1> updatedPoint, unsigned boundaryNodeIndex)
00129 {
00130     assert(this->GetNode(boundaryNodeIndex)->IsBoundaryNode());
00131     double scaleFactor = updatedPoint[0] / this->GetNode(boundaryNodeIndex)->GetPoint()[0];
00132     double temp;
00133     for (unsigned i=0; i < boundaryNodeIndex+1; i++)
00134     {
00135         temp = scaleFactor * this->mNodes[i]->GetPoint()[0];
00136         ChastePoint<1> newPoint(temp);
00137         this->mNodes[i]->SetPoint(newPoint);
00138     }
00139     this->RefreshMesh();
00140 }
00141 
00142 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00143 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::SetNode(unsigned index,
00144         ChastePoint<SPACE_DIM> point,
00145         bool concreteMove)
00146 {
00147     this->mNodes[index]->SetPoint(point);
00148 
00149     if (concreteMove)
00150     {
00151         for (typename Node<SPACE_DIM>::ContainingElementIterator it = this->mNodes[index]->ContainingElementsBegin();
00152              it != this->mNodes[index]->ContainingElementsEnd();
00153              ++it)
00154         {
00155             if (ELEMENT_DIM == SPACE_DIM)
00156             {
00157                 try
00158                 {
00159                     GetElement(*it)->CalculateInverseJacobian(this->mElementJacobians[ (*it) ],
00160                                                               this->mElementJacobianDeterminants[ (*it) ],
00161                                                               this->mElementInverseJacobians[ (*it) ]);
00162                 }
00163                 catch (Exception e)
00164                 {
00165                         EXCEPTION("Moving node caused an element to have a non-positive Jacobian determinant");
00166                 }
00167             }
00168             else
00169             {
00170                 c_vector<double,SPACE_DIM> previous_direction = this->mElementWeightedDirections[ (*it) ];
00171 
00172                 GetElement(*it)->CalculateWeightedDirection(this->mElementWeightedDirections[ (*it) ],
00173                                                             this->mElementJacobianDeterminants[ (*it) ]);
00174 
00175                 if ( inner_prod(previous_direction, this->mElementWeightedDirections[ (*it) ]) < 0)
00176                 {
00177                     EXCEPTION("Moving node caused an subspace element to change direction");
00178                 }
00179 
00180             }
00181         }
00182         for (typename Node<SPACE_DIM>::ContainingBoundaryElementIterator it = this->mNodes[index]->ContainingBoundaryElementsBegin();
00183              it != this->mNodes[index]->ContainingBoundaryElementsEnd();
00184              ++it)
00185         {
00186             try
00187             {
00188                 GetBoundaryElement(*it)->CalculateWeightedDirection(this->mBoundaryElementWeightedDirections[ (*it) ],
00189                                                                     this->mBoundaryElementJacobianDeterminants[ (*it) ]);
00190             }
00191             catch (Exception e)
00192             {
00193                 EXCEPTION("Moving node caused a boundary element to have a non-positive Jacobian determinant");
00194             }
00195         }
00196     }
00197 }
00198 
00199 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00200 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::DeleteNode(unsigned index)
00201 {
00202     if (this->mNodes[index]->IsDeleted())
00203     {
00204         EXCEPTION("Trying to delete a deleted node");
00205     }
00206     unsigned target_index = (unsigned)(-1);
00207     bool found_target=false;
00208     for (typename Node<SPACE_DIM>::ContainingElementIterator it = this->mNodes[index]->ContainingElementsBegin();
00209          !found_target && it != this->mNodes[index]->ContainingElementsEnd();
00210          ++it)
00211     {
00212         Element <ELEMENT_DIM,SPACE_DIM>* p_element = GetElement(*it);
00213         for (unsigned i=0; i<=ELEMENT_DIM && !found_target; i++)
00214         {
00215             target_index = p_element->GetNodeGlobalIndex(i);
00216             try
00217             {
00218                 MoveMergeNode(index, target_index, false);
00219                 found_target = true;
00220             }
00221             catch (Exception e)
00222             {
00223                 // Just try the next node
00224             }
00225         }
00226     }
00227     if (!found_target)
00228     {
00229         EXCEPTION("Failure to delete node");
00230     }
00231 
00232     MoveMergeNode(index, target_index);
00233 }
00234 
00235 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00236 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::DeleteNodePriorToReMesh(unsigned index)
00237 {
00238     this->mNodes[index]->MarkAsDeleted();
00239     mDeletedNodeIndices.push_back(index);
00240 }
00241 
00242 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00243 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::MoveMergeNode(unsigned index,
00244         unsigned targetIndex,
00245         bool concreteMove)
00246 {
00247     if (this->mNodes[index]->IsDeleted())
00248     {
00249         EXCEPTION("Trying to move a deleted node");
00250     }
00251 
00252     if (index == targetIndex)
00253     {
00254         EXCEPTION("Trying to merge a node with itself");
00255     }
00256     if (this->mNodes[index]->IsBoundaryNode())
00257     {
00258         if (!this->mNodes[targetIndex]->IsBoundaryNode())
00259         {
00260             EXCEPTION("A boundary node can only be moved on to another boundary node");
00261         }
00262     }
00263     std::set<unsigned> unshared_element_indices;
00264     std::set_difference(this->mNodes[index]->rGetContainingElementIndices().begin(),
00265                         this->mNodes[index]->rGetContainingElementIndices().end(),
00266                         this->mNodes[targetIndex]->rGetContainingElementIndices().begin(),
00267                         this->mNodes[targetIndex]->rGetContainingElementIndices().end(),
00268                         std::inserter(unshared_element_indices, unshared_element_indices.begin()));
00269 
00270 
00271     if (unshared_element_indices.size() == this->mNodes[index]->rGetContainingElementIndices().size())
00272     {
00273         EXCEPTION("These nodes cannot be merged since they are not neighbours");
00274     }
00275 
00276     std::set<unsigned> unshared_boundary_element_indices;
00277     std::set_difference(this->mNodes[index]->rGetContainingBoundaryElementIndices().begin(),
00278                         this->mNodes[index]->rGetContainingBoundaryElementIndices().end(),
00279                         this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().begin(),
00280                         this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().end(),
00281                         std::inserter(unshared_boundary_element_indices, unshared_boundary_element_indices.begin()));
00282 
00283 
00284     if (this->mNodes[index]->IsBoundaryNode())
00285     {
00286         if (unshared_boundary_element_indices.size()
00287             == this->mNodes[index]->rGetContainingBoundaryElementIndices().size())
00288         {
00289             //May be redundant (only thrown in 1D tests)
00290             EXCEPTION("These nodes cannot be merged since they are not neighbours on the boundary");
00291         }
00292     }
00293 
00294     this->mNodes[index]->rGetModifiableLocation() = this->mNodes[targetIndex]->rGetLocation();
00295 
00296     for (std::set<unsigned>::const_iterator element_iter=unshared_element_indices.begin();
00297              element_iter != unshared_element_indices.end();
00298              element_iter++)
00299     {
00300         try
00301         {
00302             if (SPACE_DIM == ELEMENT_DIM)
00303             {
00304                 this->GetElement(*element_iter)->CalculateInverseJacobian(this->mElementJacobians[(*element_iter)],
00305                                                                           this->mElementJacobianDeterminants[(*element_iter)],
00306                                                                           this->mElementInverseJacobians[ (*element_iter) ]);
00307             }
00308             else
00309             {
00310                 this->GetElement(*element_iter)->CalculateWeightedDirection(this->mElementWeightedDirections[(*element_iter)],
00311                                                                             this->mElementJacobianDeterminants[(*element_iter)]);
00312             }
00313 
00314             if (concreteMove)
00315             {
00316                 this->GetElement(*element_iter)->ReplaceNode(this->mNodes[index], this->mNodes[targetIndex]);
00317             }
00318 
00319         }
00320         catch (Exception e)
00321         {
00322             EXCEPTION("Moving node caused an element to have a non-positive Jacobian determinant");
00323         }
00324     }
00325 
00326     for (std::set<unsigned>::const_iterator boundary_element_iter=
00327                  unshared_boundary_element_indices.begin();
00328              boundary_element_iter != unshared_boundary_element_indices.end();
00329              boundary_element_iter++)
00330     {
00331 
00332         this->GetBoundaryElement(*boundary_element_iter)->CalculateWeightedDirection(this->mBoundaryElementWeightedDirections[(*boundary_element_iter)],
00333                                                                                      this->mBoundaryElementJacobianDeterminants[(*boundary_element_iter)]);
00334 
00335         if (concreteMove)
00336         {
00337             this->GetBoundaryElement(*boundary_element_iter)->ReplaceNode(this->mNodes[index], this->mNodes[targetIndex]);
00338         }
00339     }
00340 
00341     std::set<unsigned> shared_element_indices;
00342     std::set_intersection(this->mNodes[index]->rGetContainingElementIndices().begin(),
00343                           this->mNodes[index]->rGetContainingElementIndices().end(),
00344                           this->mNodes[targetIndex]->rGetContainingElementIndices().begin(),
00345                           this->mNodes[targetIndex]->rGetContainingElementIndices().end(),
00346                           std::inserter(shared_element_indices, shared_element_indices.begin()));
00347 
00348     for (std::set<unsigned>::const_iterator element_iter=shared_element_indices.begin();
00349              element_iter != shared_element_indices.end();
00350              element_iter++)
00351     {
00352         if (concreteMove)
00353         {
00354             this->GetElement(*element_iter)->MarkAsDeleted();
00355             this->mElementJacobianDeterminants[ (*element_iter) ] = 0.0; //This used to be done in MarkAsDeleted
00356             mDeletedElementIndices.push_back(*element_iter);
00357         }
00358         else
00359         {
00360             this->mElementJacobianDeterminants[ (*element_iter) ] = 0.0;
00361         }
00362     }
00363 
00364 
00365     std::set<unsigned> shared_boundary_element_indices;
00366     std::set_intersection(this->mNodes[index]->rGetContainingBoundaryElementIndices().begin(),
00367                           this->mNodes[index]->rGetContainingBoundaryElementIndices().end(),
00368                           this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().begin(),
00369                           this->mNodes[targetIndex]->rGetContainingBoundaryElementIndices().end(),
00370                           std::inserter(shared_boundary_element_indices, shared_boundary_element_indices.begin()));
00371 
00372     for (std::set<unsigned>::const_iterator boundary_element_iter=shared_boundary_element_indices.begin();
00373              boundary_element_iter != shared_boundary_element_indices.end();
00374              boundary_element_iter++)
00375     {
00376         if (concreteMove)
00377         {
00378             this->GetBoundaryElement(*boundary_element_iter)->MarkAsDeleted();
00379             this->mBoundaryElementJacobianDeterminants[ (*boundary_element_iter) ] = 0.0; //This used to be done in MarkAsDeleted
00380             mDeletedBoundaryElementIndices.push_back(*boundary_element_iter);
00381         }
00382         else
00383         {
00384             this->mBoundaryElementJacobianDeterminants[ (*boundary_element_iter) ] = 0.0;
00385             this->mBoundaryElementWeightedDirections[ (*boundary_element_iter) ] = zero_vector<double>(SPACE_DIM);
00386         }
00387     }
00388 
00389     if (concreteMove)
00390     {
00391         this->mNodes[index]->MarkAsDeleted();
00392         mDeletedNodeIndices.push_back(index);
00393     }
00394 }
00395 
00396 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00397 unsigned MutableMesh<ELEMENT_DIM, SPACE_DIM>::RefineElement(
00398     Element<ELEMENT_DIM,SPACE_DIM>* pElement,
00399     ChastePoint<SPACE_DIM> point)
00400 {
00401     //Check that the point is in the element
00402     if (pElement->IncludesPoint(point, true) == false)
00403     {
00404         EXCEPTION("RefineElement could not be started (point is not in element)");
00405     }
00406 
00407     // Add a new node from the point that is passed to RefineElement
00408     unsigned new_node_index = AddNode(new Node<SPACE_DIM>(0, point.rGetLocation()));
00409     // Note: the first argument is the index of the node, which is going to be
00410     //       overriden by AddNode, so it can safely be ignored
00411 
00412     // This loop constructs the extra elements which are going to fill the space
00413     for (unsigned i = 0; i < ELEMENT_DIM; i++)
00414     {
00415 
00416         // First, make a copy of the current element making sure we update its index
00417         unsigned new_elt_index;
00418         if (mDeletedElementIndices.empty())
00419         {
00420             new_elt_index = this->mElements.size();
00421         }
00422         else
00423         {
00424             new_elt_index = mDeletedElementIndices.back();
00425             mDeletedElementIndices.pop_back();
00426         }
00427 
00428         Element<ELEMENT_DIM,SPACE_DIM>* p_new_element=
00429             new Element<ELEMENT_DIM,SPACE_DIM>(*pElement, new_elt_index);
00430 
00431         // Second, update the node in the element with the new one
00432         p_new_element->UpdateNode(ELEMENT_DIM-1-i, this->mNodes[new_node_index]);
00433 
00434         // Third, add the new element to the set
00435         if ((unsigned) new_elt_index == this->mElements.size())
00436         {
00437             this->mElements.push_back(p_new_element);
00438         }
00439         else
00440         {
00441             delete this->mElements[new_elt_index];
00442             this->mElements[new_elt_index] = p_new_element;
00443         }
00444     }
00445 
00446     // Lastly, update the last node in the element to be refined
00447     pElement->UpdateNode(ELEMENT_DIM, this->mNodes[new_node_index]);
00448 
00449     return new_node_index;
00450 }
00451 
00452 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00453 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::DeleteBoundaryNodeAt(unsigned index)
00454 {
00455     if (!this->mNodes[index]->IsBoundaryNode() )
00456     {
00457         EXCEPTION(" You may only delete a boundary node ");
00458     }
00459 
00460     this->mNodes[index]->MarkAsDeleted();
00461     mDeletedNodeIndices.push_back(index);
00462 
00463     // Update the boundary node vector
00464     typename std::vector<Node<SPACE_DIM>*>::iterator b_node_iter
00465     = std::find(this->mBoundaryNodes.begin(), this->mBoundaryNodes.end(), this->mNodes[index]);
00466     this->mBoundaryNodes.erase(b_node_iter);
00467 
00468     // Remove boundary elements containing this node
00469     std::set<unsigned> boundary_element_indices = this->mNodes[index]->rGetContainingBoundaryElementIndices();
00470     std::set<unsigned>::const_iterator boundary_element_indices_iterator = boundary_element_indices.begin();
00471     while (boundary_element_indices_iterator != boundary_element_indices.end())
00472     {
00473         BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>* p_boundary_element = this->GetBoundaryElement(*boundary_element_indices_iterator);
00474         p_boundary_element->MarkAsDeleted();
00475         mDeletedBoundaryElementIndices.push_back(*boundary_element_indices_iterator);
00476         boundary_element_indices_iterator++;
00477     }
00478 
00479     // Remove elements containing this node
00480     std::set<unsigned> element_indices = this->mNodes[index]->rGetContainingElementIndices();
00481     std::set<unsigned>::const_iterator element_indices_iterator = element_indices.begin();
00482     while (element_indices_iterator != element_indices.end())
00483     {
00484         Element<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(*element_indices_iterator);
00485         for (unsigned i=0; i<p_element->GetNumNodes(); i++)
00486         {
00487             Node<SPACE_DIM>* p_node = p_element->GetNode(i);
00488             if (!p_node->IsDeleted())
00489             {
00490                 p_node->SetAsBoundaryNode();
00491                 // Update the boundary node vector
00492                 this->mBoundaryNodes.push_back(p_node);
00493             }
00494         }
00495         p_element->MarkAsDeleted();
00496         mDeletedElementIndices.push_back(p_element->GetIndex());
00497         element_indices_iterator++;
00498     }
00499 }
00500 
00501 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00502 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::ReIndex(NodeMap& map)
00503 {
00504     assert(!mAddedNodes);
00505     map.Resize(this->GetNumAllNodes());
00506 
00507     std::vector<Element<ELEMENT_DIM, SPACE_DIM> *> live_elements;
00508 
00509     for (unsigned i=0; i<this->mElements.size(); i++)
00510     {
00511         assert(i==this->mElements[i]->GetIndex()); // We need this to be true to be able to reindex the jacobian cache
00512         if (this->mElements[i]->IsDeleted())
00513         {
00514             delete this->mElements[i];
00515         }
00516         else
00517         {
00518             live_elements.push_back(this->mElements[i]);
00519 
00520             if (SPACE_DIM == ELEMENT_DIM)
00521             {
00522                 this->mElementJacobians[live_elements.size()-1] = this->mElementJacobians[this->mElements[i]->GetIndex()];
00523                 this->mElementInverseJacobians[live_elements.size()-1] = this->mElementInverseJacobians[this->mElements[i]->GetIndex()];
00524             }
00525             else
00526             {
00527                 this->mElementWeightedDirections[live_elements.size()-1] = this->mElementWeightedDirections[this->mElements[i]->GetIndex()];
00528             }
00529             this->mElementJacobianDeterminants[live_elements.size()-1] = this->mElementJacobianDeterminants[this->mElements[i]->GetIndex()];
00530         }
00531     }
00532 
00533     assert(mDeletedElementIndices.size() == this->mElements.size()-live_elements.size());
00534     mDeletedElementIndices.clear();
00535     this->mElements = live_elements;
00536     unsigned num_elements = this->mElements.size();
00537 
00538     if (SPACE_DIM == ELEMENT_DIM)
00539     {
00540         this->mElementJacobians.resize(num_elements);
00541         this->mElementInverseJacobians.resize(num_elements);
00542     }
00543     else
00544     {
00545         this->mElementWeightedDirections.resize(num_elements);
00546     }
00547     this->mElementJacobianDeterminants.resize(num_elements);
00548 
00549     std::vector<Node<SPACE_DIM> *> live_nodes;
00550     for (unsigned i=0; i<this->mNodes.size(); i++)
00551     {
00552         if (this->mNodes[i]->IsDeleted())
00553         {
00554             delete this->mNodes[i];
00555             map.SetDeleted(i);
00556         }
00557         else
00558         {
00559             live_nodes.push_back(this->mNodes[i]);
00560             // the nodes will have their index set to be the index into the live_nodes
00561             // vector further down
00562             map.SetNewIndex(i, (unsigned)(live_nodes.size()-1));
00563         }
00564     }
00565 
00566     assert(mDeletedNodeIndices.size() == this->mNodes.size()-live_nodes.size());
00567     this->mNodes = live_nodes;
00568     mDeletedNodeIndices.clear();
00569 
00570     std::vector<BoundaryElement<ELEMENT_DIM-1, SPACE_DIM> *> live_boundary_elements;
00571     for (unsigned i=0; i<this->mBoundaryElements.size(); i++)
00572     {
00573         if (this->mBoundaryElements[i]->IsDeleted())
00574         {
00575             delete this->mBoundaryElements[i];
00576         }
00577         else
00578         {
00579             live_boundary_elements.push_back(this->mBoundaryElements[i]);
00580 
00581             this->mBoundaryElementWeightedDirections[live_boundary_elements.size()-1] = this->mBoundaryElementWeightedDirections[this->mBoundaryElements[i]->GetIndex()];
00582             this->mBoundaryElementJacobianDeterminants[live_boundary_elements.size()-1] = this->mBoundaryElementJacobianDeterminants[this->mBoundaryElements[i]->GetIndex()];
00583         }
00584     }
00585 
00586     assert(mDeletedBoundaryElementIndices.size() == this->mBoundaryElements.size()-live_boundary_elements.size());
00587     this->mBoundaryElements = live_boundary_elements;
00588     mDeletedBoundaryElementIndices.clear();
00589 
00590     unsigned num_boundary_elements = this->mBoundaryElements.size();
00591 
00592     this->mBoundaryElementWeightedDirections.resize(num_boundary_elements);
00593     this->mBoundaryElementJacobianDeterminants.resize(num_boundary_elements);
00594 
00595     for (unsigned i=0; i<this->mNodes.size(); i++)
00596     {
00597         this->mNodes[i]->SetIndex(i);
00598     }
00599     for (unsigned i=0; i<this->mElements.size(); i++)
00600     {
00601 
00602         this->mElements[i]->ResetIndex(i);
00603     }
00604     for (unsigned i=0; i<this->mBoundaryElements.size(); i++)
00605     {
00606         this->mBoundaryElements[i]->ResetIndex(i);
00607     }
00608 }
00609 
00610 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00611 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::ReMesh(NodeMap& map)
00612 {
00613     // Make sure that we are in the correct dimension - this code will be eliminated at compile time
00614     #define COVERAGE_IGNORE
00615     assert( ELEMENT_DIM == SPACE_DIM );
00616     #undef COVERAGE_IGNORE
00617 
00618     // Avoid some triangle/tetgen errors: need at least four
00619     // nodes for tetgen, and at least three for triangle
00620     if (GetNumNodes() <= SPACE_DIM)
00621     {
00622         EXCEPTION("The number of nodes must exceed the spatial dimension.");
00623     }
00624 
00625     // Make sure the map is big enough
00626     map.Resize(this->GetNumAllNodes());
00627     if (mAddedNodes || !mDeletedNodeIndices.empty())
00628     {
00629         //Size of mesh is about to change
00630         if (this->mpDistributedVectorFactory)
00631         {
00632             delete this->mpDistributedVectorFactory;
00633             this->mpDistributedVectorFactory =  new DistributedVectorFactory(this->GetNumNodes());
00634         }
00635     }
00636     if (SPACE_DIM==1)
00637     {
00638         // Store the node locations
00639         std::vector<c_vector<double, SPACE_DIM> > old_node_locations;
00640         unsigned new_index = 0;
00641         for (unsigned i=0; i<this->GetNumAllNodes(); i++)
00642         {
00643             if (this->mNodes[i]->IsDeleted())
00644             {
00645                 map.SetDeleted(i);
00646             }
00647             else
00648             {
00649                 map.SetNewIndex(i, new_index);
00650                 old_node_locations.push_back(this->mNodes[i]->rGetLocation());
00651                 new_index++;
00652             }
00653         }
00654 
00655         // Remove current data
00656         Clear();
00657 
00658         // Construct the nodes and boundary nodes
00659         for (unsigned node_index=0; node_index<old_node_locations.size(); node_index++)
00660         {
00661             Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index, old_node_locations[node_index], false);
00662             this->mNodes.push_back(p_node);
00663 
00664             // As we're in 1D, the boundary nodes are simply at either end of the mesh
00665             if ( node_index==0 || node_index==old_node_locations.size()-1 )
00666             {
00667                 this->mBoundaryNodes.push_back(p_node);
00668             }
00669         }
00670 
00671         // Create a map between node indices and node locations
00672         std::map<double, unsigned> location_index_map;
00673         for (unsigned i=0; i<this->mNodes.size(); i++)
00674         {
00675             location_index_map[this->mNodes[i]->rGetLocation()[0]] = this->mNodes[i]->GetIndex();
00676         }
00677 
00678         // Use this map to generate a vector of node indices that are ordered spatially
00679         std::vector<unsigned> node_indices_ordered_spatially;
00680         for (std::map<double, unsigned>::iterator iter = location_index_map.begin();
00681              iter != location_index_map.end();
00682              ++iter)
00683         {
00684             node_indices_ordered_spatially.push_back(iter->second);
00685         }
00686 
00687         // Construct the elements
00688         this->mElements.reserve(old_node_locations.size()-1);
00689         for (unsigned element_index=0; element_index<old_node_locations.size()-1; element_index++)
00690         {
00691             std::vector<Node<SPACE_DIM>*> nodes;
00692             for (unsigned j=0; j<2; j++)
00693             {
00694                 unsigned global_node_index = node_indices_ordered_spatially[element_index + j];
00695                 assert(global_node_index < this->mNodes.size());
00696                 nodes.push_back(this->mNodes[global_node_index]);
00697             }
00698             this->mElements.push_back(new Element<ELEMENT_DIM, SPACE_DIM>(element_index, nodes));
00699         }
00700 
00701         // Construct the two boundary elements - as we're in 1D, these are simply at either end of the mesh
00702         std::vector<Node<SPACE_DIM>*> nodes;
00703         nodes.push_back(this->mNodes[0]);
00704         this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(0, nodes));
00705 
00706         nodes.clear();
00707         nodes.push_back(this->mNodes[old_node_locations.size()-1]);
00708         this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(1, nodes));
00709 
00710         this->RefreshJacobianCachedData();
00711     }
00712     else if (SPACE_DIM==2)  // In 2D, remesh using triangle via library calls
00713     {
00714         struct triangulateio triangle_input;
00715         triangle_input.pointlist = (double *) malloc(GetNumNodes() * 2 * sizeof(double));
00716         triangle_input.numberofpoints = GetNumNodes();
00717         triangle_input.numberofpointattributes = 0;
00718         triangle_input.pointmarkerlist = NULL;
00719         triangle_input.numberofsegments = 0;
00720         triangle_input.numberofholes = 0;
00721         triangle_input.numberofregions = 0;
00722 
00723         unsigned new_index = 0;
00724         for (unsigned i=0; i<this->GetNumAllNodes(); i++)
00725         {
00726             if (this->mNodes[i]->IsDeleted())
00727             {
00728                 map.SetDeleted(i);
00729             }
00730             else
00731             {
00732                 map.SetNewIndex(i, new_index);
00733                 triangle_input.pointlist[2*new_index] = this->mNodes[i]->rGetLocation()[0];
00734                 triangle_input.pointlist[2*new_index + 1] = this->mNodes[i]->rGetLocation()[1];
00735                 new_index++;
00736             }
00737         }
00738 
00739         // Make structure for output
00740         struct triangulateio triangle_output;
00741         triangle_output.pointlist = NULL;
00742         triangle_output.pointattributelist = (double *) NULL;
00743         triangle_output.pointmarkerlist = (int *) NULL;
00744         triangle_output.trianglelist = (int *) NULL;
00745         triangle_output.triangleattributelist = (double *) NULL;
00746         triangle_output.edgelist = (int *) NULL;
00747         triangle_output.edgemarkerlist = (int *) NULL;
00748 
00749         // Library call
00750         triangulate((char*)"Qze", &triangle_input, &triangle_output, NULL);
00751 
00752         assert(triangle_output.numberofcorners == 3);
00753 
00754         // Remove current data
00755         Clear();
00756 
00757         // Construct the nodes
00758         for (unsigned node_index=0; node_index<(unsigned)triangle_output.numberofpoints; node_index++)
00759         {
00760             if (triangle_output.pointmarkerlist[node_index] == 1)
00761             {
00762                 // Boundary node
00763                 Node<SPACE_DIM>* p_node = new Node<SPACE_DIM>(node_index, true,
00764                   triangle_output.pointlist[node_index * 2],
00765                   triangle_output.pointlist[node_index * 2+1]);
00766                 this->mNodes.push_back(p_node);
00767                 this->mBoundaryNodes.push_back(p_node);
00768             }
00769             else
00770             {
00771                 this->mNodes.push_back(new Node<SPACE_DIM>(node_index, false,
00772                   triangle_output.pointlist[node_index * 2],
00773                   triangle_output.pointlist[node_index * 2+1]));
00774             }
00775         }
00776 
00777         // Construct the elements
00778         this->mElements.reserve(triangle_output.numberoftriangles);
00779         for (unsigned element_index=0; element_index<(unsigned)triangle_output.numberoftriangles; element_index++)
00780         {
00781             std::vector<Node<SPACE_DIM>*> nodes;
00782             for (unsigned j=0; j<3; j++)
00783             {
00784                 unsigned global_node_index = triangle_output.trianglelist[element_index*3 + j];
00785                 assert(global_node_index < this->mNodes.size());
00786                 nodes.push_back(this->mNodes[global_node_index]);
00787             }
00788             this->mElements.push_back(new Element<ELEMENT_DIM, SPACE_DIM>(element_index, nodes));
00789         }
00790 
00791         // Construct the edges
00792         // too big mBoundaryElements.reserve(triangle_output.numberoftriangles);
00793         unsigned next_boundary_element_index = 0;
00794         for (unsigned boundary_element_index=0; boundary_element_index<(unsigned)triangle_output.numberofedges; boundary_element_index++)
00795         {
00796             if (triangle_output.edgemarkerlist[boundary_element_index] == 1)
00797             {
00798                 std::vector<Node<SPACE_DIM>*> nodes;
00799                 for (unsigned j=0; j<2; j++)
00800                 {
00801                     unsigned global_node_index = triangle_output.edgelist[boundary_element_index*2 + j];
00802                     assert(global_node_index < this->mNodes.size());
00803                     nodes.push_back(this->mNodes[global_node_index]);
00804                 }
00805                 this->mBoundaryElements.push_back(new BoundaryElement<ELEMENT_DIM-1, SPACE_DIM>(next_boundary_element_index++, nodes));
00806             }
00807         }
00808 
00809         free(triangle_input.pointlist);
00810 
00811         free(triangle_output.pointlist);
00812         free(triangle_output.pointattributelist);
00813         free(triangle_output.pointmarkerlist);
00814         free(triangle_output.trianglelist);
00815         free(triangle_output.triangleattributelist);
00816         free(triangle_output.edgelist);
00817         free(triangle_output.edgemarkerlist);
00818 
00819         this->RefreshJacobianCachedData();
00820     }
00821     else // in 3D, remesh using tetgen
00822     {
00823         std::stringstream pid;
00824         pid << getpid();
00825 
00826         OutputFileHandler handler("");
00827         std::string full_name = handler.GetOutputDirectoryFullPath() + "temp_" + pid.str() + ".";
00828 
00829         // Only the master process should do IO and call the mesher
00830         if (PetscTools::AmMaster())
00831         {
00832             std::string node_file_name = "temp_" + pid.str() + ".node";
00833             {
00834                 out_stream node_file = handler.OpenOutputFile(node_file_name);
00835 
00836                 (*node_file) << GetNumNodes() << "\t" << SPACE_DIM << "\t0\t0\n";
00837 
00838                 unsigned new_index = 0;
00839 
00840                 for (unsigned i=0; i<this->GetNumAllNodes(); i++)
00841                 {
00842                     if (this->mNodes[i]->IsDeleted())
00843                     {
00844                         map.SetDeleted(i);
00845                     }
00846                     else
00847                     {
00848                         map.SetNewIndex(i, new_index);
00849                         new_index++;
00850                         const c_vector<double, SPACE_DIM> node_loc = this->mNodes[i]->rGetLocation();
00851                         (*node_file) << i << "\t" << node_loc[0] << "\t" << node_loc[1] << "\t" << node_loc[2] << "\n";
00852                     }
00853                 }
00854                 node_file->close();
00855             }
00856 
00857             std::string binary_name;
00858 
00859             binary_name="tetgen"; //Assume it's in the path
00860             std::string command = binary_name + " -Q " + full_name + "node";
00861 
00862             // Tetgen's quiet mode isn't as quiet as Triangle's
00863             command += " > /dev/null";
00864 
00865             int return_value = system(command.c_str());
00866             if (return_value != 0)
00867             {
00868 #define COVERAGE_IGNORE
00869                 EXCEPTION("The tetgen mesher did not succeed in remeshing.  This functionality relies on tetgen.  Do you have tetgen from http://tetgen.berlios.de/ in your path?");
00870 #undef COVERAGE_IGNORE
00871             }
00872         }
00873         // Wait for the new mesh to be available and communicate its name
00874 #ifndef SPECIAL_SERIAL
00875         if (!PetscTools::IsSequential())
00876         {
00877             char full_name_comm[200]; 
00878             strcpy(full_name_comm, full_name.c_str());
00879             MPI_Bcast(full_name_comm, 200, MPI_CHAR, 0, MPI_COMM_WORLD);
00880             full_name = full_name_comm;
00881         }
00882 #endif //SPECIAL_SERIAL
00883 
00884         // Clear all current data
00885         Clear();
00886 
00887         // Read the new mesh back from file
00888         TrianglesMeshReader<ELEMENT_DIM, SPACE_DIM> mesh_reader(full_name+"1");
00889         ConstructFromMeshReader(mesh_reader);
00890 
00891         // Make sure the file is not deleted before all the processors have read it
00892         PetscTools::Barrier("MutableMesh::ReMesh");
00893 
00894         if (PetscTools::AmMaster())
00895         {
00896             // delete the temporary files (one by one rather than using * to make it impossible
00897             // to ever accidentally end up with "rm -f *").
00898 
00900             std::string remove_command = "rm -f " + full_name + "node";
00901             system(remove_command.c_str());
00902             //std::cout << remove_command << "\n";
00903 
00904             remove_command = "rm -f " + full_name + "1.node";
00905             system(remove_command.c_str());
00906             //std::cout << remove_command << "\n";
00907 
00908             remove_command = "rm -f " + full_name + "1.ele";
00909             system(remove_command.c_str());
00910             //std::cout << remove_command << "\n";
00911 
00912             remove_command = "rm -f " + full_name + "1.face";
00913             system(remove_command.c_str());
00914             //std::cout << remove_command << "\n";
00915         }
00916     }
00917 }
00918 
00919 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00920 void MutableMesh<ELEMENT_DIM, SPACE_DIM>::ReMesh()
00921 {
00922     NodeMap map(GetNumNodes());
00923     ReMesh(map);
00924 }
00925 
00926 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00927 bool MutableMesh<ELEMENT_DIM, SPACE_DIM>::CheckIsVoronoi(Element<ELEMENT_DIM, SPACE_DIM>* pElement, double maxPenetration)
00928 {
00929     assert(ELEMENT_DIM == SPACE_DIM);
00930     unsigned num_nodes = pElement->GetNumNodes();
00931     std::set<unsigned> neighbouring_elements_indices;
00932     std::set< Element<ELEMENT_DIM,SPACE_DIM> *> neighbouring_elements;
00933     std::set<unsigned> neighbouring_nodes_indices;
00934 
00935     // Form a set of neighbouring elements via the nodes
00936     for (unsigned i=0; i<num_nodes; i++)
00937     {
00938         Node<SPACE_DIM>* p_node = pElement->GetNode(i);
00939         neighbouring_elements_indices = p_node->rGetContainingElementIndices();
00941         for (std::set<unsigned>::const_iterator it = neighbouring_elements_indices.begin();
00942              it != neighbouring_elements_indices.end();
00943              ++it)
00944         {
00945             neighbouring_elements.insert(this->GetElement(*it));
00946         }
00947     }
00948     neighbouring_elements.erase(pElement);
00949 
00950     // For each neighbouring element find the supporting nodes
00951     typedef typename std::set<Element<ELEMENT_DIM,SPACE_DIM> *>::const_iterator ElementIterator;
00952 
00953     for (ElementIterator it = neighbouring_elements.begin();
00954          it != neighbouring_elements.end();
00955          ++it)
00956     {
00957         for (unsigned i=0; i<num_nodes; i++)
00958         {
00959             neighbouring_nodes_indices.insert((*it)->GetNodeGlobalIndex(i));
00960         }
00961     }
00962 
00963     // Remove the nodes that support this element
00964     for (unsigned i = 0; i < num_nodes; i++)
00965     {
00966         neighbouring_nodes_indices.erase(pElement->GetNodeGlobalIndex(i));
00967     }
00968 
00969     // Get the circumsphere information
00970     c_vector<double, SPACE_DIM+1> this_circum_centre;
00971 
00972     this_circum_centre = pElement->CalculateCircumsphere(this->mElementJacobians[pElement->GetIndex()], this->mElementInverseJacobians[pElement->GetIndex()]);
00973 
00974     // Copy the actualy circumcentre into a smaller vector
00975     c_vector<double, ELEMENT_DIM> circum_centre;
00976     for (unsigned i=0; i<ELEMENT_DIM; i++)
00977     {
00978         circum_centre[i] = this_circum_centre[i];
00979     }
00980 
00981     for (std::set<unsigned>::const_iterator it = neighbouring_nodes_indices.begin();
00982          it != neighbouring_nodes_indices.end();
00983          ++it)
00984     {
00985         c_vector<double, ELEMENT_DIM> node_location = this->GetNode(*it)->rGetLocation();
00986 
00987         // Calculate vector from circumcenter to node
00988         node_location -= circum_centre;
00989 
00990         // This is to calculate the squared distance betweeen them
00991         double squared_distance = inner_prod(node_location, node_location);
00992 
00993         // If the squared idstance is less than the elements circum-radius(squared),
00994         // then the Voronoi property is violated.
00995         if (squared_distance < this_circum_centre[ELEMENT_DIM])
00996         {
00997             // We know the node is inside the circumsphere, but we don't know how far
00998             double radius = sqrt(this_circum_centre[ELEMENT_DIM]);
00999             double distance = radius - sqrt(squared_distance);
01000 
01001             // If the node penetration is greater than supplied maximum penetration factor
01002             if (distance/radius > maxPenetration)
01003             {
01004                 return false;
01005             }
01006         }
01007     }
01008     return true;
01009 }
01010 
01011 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01012 bool MutableMesh<ELEMENT_DIM, SPACE_DIM>::CheckIsVoronoi(double maxPenetration)
01013 {
01014     // Looping through all the elements in the mesh
01016     for (unsigned i=0; i<this->mElements.size(); i++)
01017     {
01018         // Check if the element is not deleted
01019         if (!this->mElements[i]->IsDeleted())
01020         {
01021             // Checking the Voronoi of the Element
01022             if (CheckIsVoronoi(this->mElements[i], maxPenetration) == false)
01023             {
01024                 return false;
01025             }
01026         }
01027     }
01028     return true;
01029 }
01030 
01031 
01033 // Explicit instantiation
01035 
01036 template class MutableMesh<1,1>;
01037 template class MutableMesh<1,2>;
01038 template class MutableMesh<1,3>;
01039 template class MutableMesh<2,2>;
01040 template class MutableMesh<2,3>;
01041 template class MutableMesh<3,3>;
01042 
01043 
01044 // Serialization for Boost >= 1.36
01045 #include "SerializationExportWrapperForCpp.hpp"
01046 EXPORT_TEMPLATE_CLASS_ALL_DIMS(MutableMesh);

Generated by  doxygen 1.6.2