VertexMesh.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 "VertexMesh.hpp"
00030 #include "RandomNumberGenerator.hpp"
00031 #include "UblasCustomFunctions.hpp"
00032 #include <list>
00033 
00038 bool IndexAngleComparison(const std::pair<unsigned, double> lhs, const std::pair<unsigned, double> rhs)
00039 {
00040     return lhs.second < rhs.second;
00041 }
00042 
00043 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00044 VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexMesh(std::vector<Node<SPACE_DIM>*> nodes,
00045                                                std::vector<VertexElement<ELEMENT_DIM,SPACE_DIM>*> vertexElements)
00046     : mpDelaunayMesh(NULL)
00047 {
00048 
00049     // Reset member variables and clear mNodes and mElements
00050     Clear();
00051 
00052     // Populate mNodes and mElements
00053     for (unsigned node_index=0; node_index<nodes.size(); node_index++)
00054     {
00055         Node<SPACE_DIM>* p_temp_node = nodes[node_index];
00056         this->mNodes.push_back(p_temp_node);
00057     }
00058     for (unsigned elem_index=0; elem_index<vertexElements.size(); elem_index++)
00059     {
00060         VertexElement<ELEMENT_DIM, SPACE_DIM>* p_temp_vertex_element = vertexElements[elem_index];
00061         mElements.push_back(p_temp_vertex_element);
00062     }
00063 
00064     // In 3D, populate mFaces
00065     if (SPACE_DIM == 3)
00066     {
00067         // Use a std::set to keep track of which faces have been added to mFaces
00068         std::set<unsigned> faces_counted;
00069 
00070         // Loop over mElements
00071         for (unsigned elem_index=0; elem_index<mElements.size(); elem_index++)
00072         {
00073             // Loop over faces of this element
00074             for (unsigned face_index=0; face_index<mElements[elem_index]->GetNumFaces(); face_index++)
00075             {
00076                 VertexElement<ELEMENT_DIM-1, SPACE_DIM>* p_face = mElements[elem_index]->GetFace(face_index);
00077                 unsigned global_index = p_face->GetIndex();
00078 
00079                 // If this face is not already contained in mFaces, add it, and update faces_counted
00080                 if (faces_counted.find(global_index) == faces_counted.end())
00081                 {
00082                     mFaces.push_back(p_face);
00083                     faces_counted.insert(global_index);
00084                 }
00085             }
00086         }
00087     }
00088 
00089     // Register elements with nodes
00090     for (unsigned index=0; index<mElements.size(); index++)
00091     {
00092         VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = mElements[index];
00093 
00094         unsigned element_index = p_element->GetIndex();
00095         unsigned num_nodes_in_element = p_element->GetNumNodes();
00096 
00097         for (unsigned node_index=0; node_index<num_nodes_in_element; node_index++)
00098         {
00099             p_element->GetNode(node_index)->AddElement(element_index);
00100         }
00101     }
00102 
00103     this->mMeshChangesDuringSimulation = false;
00104 }
00105 
00106 
00107 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00108 VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexMesh(std::vector<Node<SPACE_DIM>*> nodes,
00109                            std::vector<VertexElement<ELEMENT_DIM-1, SPACE_DIM>*> faces,
00110                            std::vector<VertexElement<ELEMENT_DIM, SPACE_DIM>*> vertexElements)
00111     : mpDelaunayMesh(NULL)
00112 {
00113     // Reset member variables and clear mNodes, mFaces and mElements
00114     Clear();
00115 
00116     // Populate mNodes mFaces and mElements
00117     for (unsigned node_index=0; node_index<nodes.size(); node_index++)
00118     {
00119         Node<SPACE_DIM>* p_temp_node = nodes[node_index];
00120         this->mNodes.push_back(p_temp_node);
00121     }
00122 
00123     for (unsigned face_index=0; face_index<faces.size(); face_index++)
00124     {
00125         VertexElement<ELEMENT_DIM-1, SPACE_DIM>* p_temp_face = faces[face_index];
00126         mFaces.push_back(p_temp_face);
00127     }
00128 
00129     for (unsigned elem_index=0; elem_index<vertexElements.size(); elem_index++)
00130     {
00131         VertexElement<ELEMENT_DIM, SPACE_DIM>* p_temp_vertex_element = vertexElements[elem_index];
00132         mElements.push_back(p_temp_vertex_element);
00133     }
00134 
00135     // Register elements with nodes
00136     for (unsigned index=0; index<mElements.size(); index++)
00137     {
00138         VertexElement<ELEMENT_DIM, SPACE_DIM>* p_temp_vertex_element = mElements[index];
00139         for (unsigned node_index=0; node_index<p_temp_vertex_element->GetNumNodes(); node_index++)
00140         {
00141             Node<SPACE_DIM>* p_temp_node = p_temp_vertex_element->GetNode(node_index);
00142             p_temp_node->AddElement(p_temp_vertex_element->GetIndex());
00143         }
00144     }
00145 
00146     this->mMeshChangesDuringSimulation = false;
00147 }
00148 
00149 
00156 template<>
00157 VertexMesh<2,2>::VertexMesh(TetrahedralMesh<2,2>& rMesh, bool isPeriodic)
00158     : mpDelaunayMesh(&rMesh)
00159 {
00160     if (!isPeriodic)        // If the mesh is non-periodic
00161     {
00162         // Reset member variables and clear mNodes, mFaces and mElements
00163         Clear();
00164     
00165         unsigned num_elements = mpDelaunayMesh->GetNumAllNodes();
00166         unsigned num_nodes = mpDelaunayMesh->GetNumAllElements();
00167     
00168         // Allocate memory for mNodes and mElements
00169         this->mNodes.reserve(num_nodes);
00170     
00171         // Create as many elements as there are nodes in the mesh
00172         mElements.reserve(num_elements);
00173         for (unsigned elem_index=0; elem_index<num_elements; elem_index++)
00174         {
00175             VertexElement<2,2>* p_element = new VertexElement<2,2>(elem_index);
00176             mElements.push_back(p_element);
00177         }
00178     
00179         // Populate mNodes
00180         GenerateVerticesFromElementCircumcentres(rMesh);
00181     
00182         // Loop over elements of the Delaunay mesh
00183         for (unsigned i=0; i<num_nodes; i++)
00184         {
00185             // Loop over nodes owned by this element in the Delaunay mesh
00186             for (unsigned local_index=0; local_index<3; local_index++)
00187             {
00188                 unsigned elem_index = mpDelaunayMesh->GetElement(i)->GetNodeGlobalIndex(local_index);
00189                 unsigned num_nodes_in_elem = mElements[elem_index]->GetNumNodes();
00190                 unsigned end_index = num_nodes_in_elem>0 ? num_nodes_in_elem-1 : 0;
00191     
00192                 mElements[elem_index]->AddNode(end_index, this->mNodes[i]);
00193             }
00194         }
00195     
00196         // Reorder mNodes anticlockwise
00197         for (unsigned elem_index=0; elem_index<mElements.size(); elem_index++)
00198         {
00204             std::list<std::pair<unsigned, double> > index_angle_list;
00205             for (unsigned local_index=0; local_index<mElements[elem_index]->GetNumNodes(); local_index++)
00206             {
00207                 c_vector<double, 2> vectorA = mpDelaunayMesh->GetNode(elem_index)->rGetLocation();
00208                 c_vector<double, 2> vectorB = mElements[elem_index]->GetNodeLocation(local_index);
00209                 c_vector<double, 2> centre_to_vertex = mpDelaunayMesh->GetVectorFromAtoB(vectorA, vectorB);
00210     
00211                 double angle = atan2(centre_to_vertex(1), centre_to_vertex(0));
00212                 unsigned global_index = mElements[elem_index]->GetNodeGlobalIndex(local_index);
00213     
00214                 std::pair<unsigned, double> pair(global_index, angle);
00215                 index_angle_list.push_back(pair);
00216             }
00217     
00218             // Sort the list in order of increasing angle
00219             index_angle_list.sort(IndexAngleComparison);
00220     
00221             // Create a new Voronoi element and pass in the appropriate Nodes, ordered anticlockwise
00222             VertexElement<2,2>* p_new_element = new VertexElement<2,2>(elem_index);
00223             unsigned count = 0;
00224             for (std::list<std::pair<unsigned, double> >::iterator list_iter = index_angle_list.begin();
00225                  list_iter != index_angle_list.end();
00226                  ++list_iter)
00227             {
00228                 unsigned local_index = count>1 ? count-1 : 0;
00229                 p_new_element->AddNode(local_index, mNodes[list_iter->first]);
00230                 count++;
00231             }
00232     
00233             // Replace the relevant member of mElements with this Voronoi element
00234             delete mElements[elem_index];
00235             mElements[elem_index] = p_new_element;
00236         }
00237     
00238         this->mMeshChangesDuringSimulation = false;
00239     
00240     }
00241     
00242     else    // For a periodic crypt
00243     {       
00244         // Reset member variables and clear mNodes, mFaces and mElements
00245         Clear();
00246     
00247         unsigned num_elements = mpDelaunayMesh->GetNumAllNodes();
00248         unsigned num_nodes = mpDelaunayMesh->GetNumAllElements();
00249     
00250         // Allocate memory for mNodes and mElements
00251         this->mNodes.reserve(num_nodes);
00252     
00253         // Create as many elements as there are nodes in the mesh
00254         mElements.reserve(num_elements);
00255         for (unsigned elem_index=0; elem_index<num_elements; elem_index++)
00256         {
00257             VertexElement<2,2>* p_element = new VertexElement<2,2>(elem_index);
00258             mElements.push_back(p_element);
00259         }
00260     
00261         // Populate mNodes
00262         GenerateVerticesFromElementCircumcentres(rMesh);
00263     
00264         // Loop over elements of the Delaunay mesh
00265         for (unsigned i=0; i<num_nodes; i++)
00266         {
00267             // Loop over nodes owned by this element in the Delaunay mesh
00268             for (unsigned local_index=0; local_index<3; local_index++)
00269             {
00270                 unsigned elem_index = mpDelaunayMesh->GetElement(i)->GetNodeGlobalIndex(local_index);
00271                 unsigned num_nodes_in_elem = mElements[elem_index]->GetNumNodes();
00272                 unsigned end_index = num_nodes_in_elem>0 ? num_nodes_in_elem-1 : 0;
00273     
00274                 mElements[elem_index]->AddNode(end_index, this->mNodes[i]);
00275             }
00276         }
00277     
00278         // Reorder mNodes anticlockwise
00279         for (unsigned elem_index=0; elem_index<mElements.size(); elem_index++)
00280         {
00286             std::list<std::pair<unsigned, double> > index_angle_list;
00287             for (unsigned local_index=0; local_index<mElements[elem_index]->GetNumNodes(); local_index++)
00288             {
00289                 c_vector<double, 2> vectorA = mpDelaunayMesh->GetNode(elem_index)->rGetLocation();
00290                 c_vector<double, 2> vectorB = mElements[elem_index]->GetNodeLocation(local_index);
00291                 c_vector<double, 2> centre_to_vertex = mpDelaunayMesh->GetVectorFromAtoB(vectorA, vectorB);
00292     
00293                 double angle = atan2(centre_to_vertex(1), centre_to_vertex(0));
00294                 unsigned global_index = mElements[elem_index]->GetNodeGlobalIndex(local_index);
00295     
00296                 std::pair<unsigned, double> pair(global_index, angle);
00297                 index_angle_list.push_back(pair);
00298             }
00299     
00300             // Sort the list in order of increasing angle
00301             index_angle_list.sort(IndexAngleComparison);
00302     
00303             // Create a new Voronoi element and pass in the appropriate Nodes, ordered anticlockwise
00304             VertexElement<2,2>* p_new_element = new VertexElement<2,2>(elem_index);
00305             unsigned count = 0;
00306             for (std::list<std::pair<unsigned, double> >::iterator list_iter = index_angle_list.begin();
00307                  list_iter != index_angle_list.end();
00308                  ++list_iter)
00309             {
00310                 unsigned local_index = count>1 ? count-1 : 0;
00311                 p_new_element->AddNode(local_index, mNodes[list_iter->first]);
00312                 count++;
00313             }
00314     
00315             // Replace the relevant member of mElements with this Voronoi element
00316             delete mElements[elem_index];
00317             mElements[elem_index] = p_new_element;
00318         }
00319     
00320         this->mMeshChangesDuringSimulation = false;     
00321     }
00322     
00323 }
00324 
00325 
00331 template<>
00332 VertexMesh<3,3>::VertexMesh(TetrahedralMesh<3,3>& rMesh)
00333     : mpDelaunayMesh(&rMesh)
00334 {
00335     // Reset member variables and clear mNodes, mFaces and mElements
00336     Clear();
00337 
00338     unsigned num_nodes = mpDelaunayMesh->GetNumAllElements();
00339 
00340     // Allocate memory for mNodes
00341     this->mNodes.reserve(num_nodes);
00342 
00343     // Populate mNodes
00344     GenerateVerticesFromElementCircumcentres(rMesh);
00345 
00346     std::map<unsigned, VertexElement<3,3>*> index_element_map;
00347     unsigned face_index = 0;
00348     unsigned element_index = 0;
00349 
00350     // Loop over each edge of the Delaunay mesh and populate mFaces and mElements
00351     for (TetrahedralMesh<3,3>::EdgeIterator edge_iterator = mpDelaunayMesh->EdgesBegin();
00352          edge_iterator != mpDelaunayMesh->EdgesEnd();
00353          ++edge_iterator)
00354     {
00355         Node<3>* p_node_a = edge_iterator.GetNodeA();
00356         Node<3>* p_node_b = edge_iterator.GetNodeB();
00357 
00358         if ( !(p_node_a->IsBoundaryNode() && p_node_b->IsBoundaryNode()) )
00359         {
00360             std::set<unsigned>& node_a_element_indices = p_node_a->rGetContainingElementIndices();
00361             std::set<unsigned>& node_b_element_indices = p_node_b->rGetContainingElementIndices();
00362             std::set<unsigned> edge_element_indices;
00363 
00364             std::set_intersection(node_a_element_indices.begin(),
00365                                   node_a_element_indices.end(),
00366                                   node_b_element_indices.begin(),
00367                                   node_b_element_indices.end(),
00368                                   std::inserter(edge_element_indices, edge_element_indices.begin()));
00369 
00370             c_vector<double,3> edge_vector = p_node_b->rGetLocation() - p_node_a->rGetLocation();
00371             c_vector<double,3> mid_edge = edge_vector*0.5 + p_node_a->rGetLocation();
00372 
00373             unsigned element0_index = *(edge_element_indices.begin());
00374 
00375             c_vector<double,3> basis_vector1 = mNodes[element0_index]->rGetLocation() - mid_edge;
00376 
00377             c_vector<double,3> basis_vector2;
00378             basis_vector2[0] = edge_vector[1]*basis_vector1[2] - edge_vector[2]*basis_vector1[1];
00379             basis_vector2[1] = edge_vector[2]*basis_vector1[0] - edge_vector[0]*basis_vector1[2];
00380             basis_vector2[2] = edge_vector[0]*basis_vector1[1] - edge_vector[1]*basis_vector1[0];
00381 
00387             std::list<std::pair<unsigned, double> > index_angle_list;
00388 
00389             // Loop over each element containing this edge (i.e. those containing both nodes of the edge)
00390             for (std::set<unsigned>::iterator index_iter = edge_element_indices.begin();
00391                  index_iter != edge_element_indices.end();
00392                  ++index_iter)
00393             {
00394                 // Calculate angle
00395                 c_vector<double, 3> vertex_vector =  mNodes[*index_iter]->rGetLocation() - mid_edge;
00396 
00397                 double local_vertex_dot_basis_vector1 = inner_prod(vertex_vector, basis_vector1);
00398                 double local_vertex_dot_basis_vector2 = inner_prod(vertex_vector, basis_vector2);
00399 
00400                 double angle = atan2(local_vertex_dot_basis_vector2, local_vertex_dot_basis_vector1);
00401 
00402                 std::pair<unsigned, double> pair(*index_iter, angle);
00403                 index_angle_list.push_back(pair);
00404             }
00405 
00406             // Sort the list in order of increasing angle
00407             index_angle_list.sort(IndexAngleComparison);
00408 
00409             // Create face
00410             VertexElement<2,3>* p_face = new VertexElement<2,3>(face_index);
00411             face_index++;
00412             unsigned count = 0;
00413             for (std::list<std::pair<unsigned, double> >::iterator list_iter = index_angle_list.begin();
00414                  list_iter != index_angle_list.end();
00415                  ++list_iter)
00416             {
00417                 unsigned local_index = count>1 ? count-1 : 0;
00418                 p_face->AddNode(local_index, mNodes[list_iter->first]);
00419                 count++;
00420             }
00421 
00422             // Add face to list of faces
00423             mFaces.push_back(p_face);
00424 
00425             // Add face to appropriate elements
00426             if (!p_node_a->IsBoundaryNode())
00427             {
00428                 unsigned node_a_index = p_node_a->GetIndex();
00429                 if (index_element_map[node_a_index])
00430                 {
00431                     // If there is already an element with this index, add the face to it...
00432                     index_element_map[node_a_index]->AddFace(p_face);
00433                 }
00434                 else
00435                 {
00436                     // ...otherwise create an element, add the face to it, and add to the map
00437                     mVoronoiElementIndexMap[node_a_index] = element_index;
00438                     VertexElement<3,3>* p_element = new VertexElement<3,3>(element_index);
00439                     element_index++;
00440                     p_element->AddFace(p_face);
00441                     index_element_map[node_a_index] = p_element;
00442                 }
00443             }
00444             if (!p_node_b->IsBoundaryNode())
00445             {
00446                 unsigned node_b_index = p_node_b->GetIndex();
00447                 if (index_element_map[node_b_index])
00448                 {
00449                     // If there is already an element with this index, add the face to it...
00450                     index_element_map[node_b_index]->AddFace(p_face);
00451                 }
00452                 else
00453                 {
00454                     // ...otherwise create an element, add the face to it, and add to the map
00455                     mVoronoiElementIndexMap[node_b_index] = element_index;
00456                     VertexElement<3,3>* p_element = new VertexElement<3,3>(element_index);
00457                     element_index++;
00458                     p_element->AddFace(p_face);
00459                     index_element_map[node_b_index] = p_element;
00460                 }
00461             }
00462         }
00463     }
00464 
00465     // Populate mElements
00466     unsigned elem_count = 0;
00467     for (std::map<unsigned, VertexElement<3,3>*>::iterator element_iter = index_element_map.begin();
00468          element_iter != index_element_map.end();
00469          ++element_iter)
00470     {
00471         mElements.push_back(element_iter->second);
00472         mVoronoiElementIndexMap[element_iter->first] = elem_count;
00473         elem_count++;
00474     }
00475 
00476     this->mMeshChangesDuringSimulation = false;
00477 }
00478 
00479 
00480 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00481 void VertexMesh<ELEMENT_DIM, SPACE_DIM>::GenerateVerticesFromElementCircumcentres(TetrahedralMesh<ELEMENT_DIM, SPACE_DIM>& rMesh)
00482 {
00483     c_matrix<double, SPACE_DIM, ELEMENT_DIM> jacobian;
00484     c_matrix<double, ELEMENT_DIM, SPACE_DIM> inverse_jacobian;
00485     double jacobian_det;
00486 
00487     // Loop over elements of the Delaunay mesh and populate mNodes
00488     for (unsigned i=0; i<rMesh.GetNumElements(); i++)
00489     {
00490         // Calculate the circumcentre of this element in the Delaunay mesh
00491         rMesh.GetInverseJacobianForElement(i, jacobian, jacobian_det, inverse_jacobian);
00492         c_vector<double, SPACE_DIM+1> circumsphere = rMesh.GetElement(i)->CalculateCircumsphere(jacobian, inverse_jacobian);
00493 
00494         c_vector<double, SPACE_DIM> circumcentre;
00495         for (unsigned j=0; j<SPACE_DIM; j++)
00496         {
00497             circumcentre(j) = circumsphere(j);
00498         }
00499 
00500         // Create a node in the Voronoi mesh at the location of this circumcentre
00501         this->mNodes.push_back(new Node<SPACE_DIM>(i, circumcentre));
00502     }
00503 }
00504 
00505 
00506 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00507 double VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetEdgeLength(unsigned elementIndex1, unsigned elementIndex2)
00508 {
00509     assert(SPACE_DIM==2);
00510 
00511     std::set<unsigned> node_indices_1;
00512     for (unsigned i=0; i<mElements[elementIndex1]->GetNumNodes(); i++)
00513     {
00514         node_indices_1.insert(mElements[elementIndex1]->GetNodeGlobalIndex(i));
00515     }
00516     std::set<unsigned> node_indices_2;
00517     for (unsigned i=0; i<mElements[elementIndex2]->GetNumNodes(); i++)
00518     {
00519         node_indices_2.insert(mElements[elementIndex2]->GetNodeGlobalIndex(i));
00520     }
00521 
00522     std::set<unsigned> shared_nodes;
00523     std::set_intersection(node_indices_1.begin(), node_indices_1.end(),
00524                           node_indices_2.begin(), node_indices_2.end(),
00525                           std::inserter(shared_nodes, shared_nodes.begin()));
00526 
00527     assert(shared_nodes.size() == 2);
00528 
00529     unsigned index1 = *(shared_nodes.begin());
00530     unsigned index2 = *(++(shared_nodes.begin()));
00531 
00532     c_vector<double, SPACE_DIM> node1_location = this->mNodes[index1]->rGetLocation();
00533     c_vector<double, SPACE_DIM> node2_location = this->mNodes[index2]->rGetLocation();
00534 
00535     double edge_length = norm_2(GetVectorFromAtoB(node1_location, node2_location));
00536     return edge_length;
00537 }
00538 
00539 
00540 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00541 VertexMesh<ELEMENT_DIM, SPACE_DIM>::VertexMesh()
00542 {
00543     mpDelaunayMesh = NULL;
00544     this->mMeshChangesDuringSimulation = false;
00545     Clear();
00546 }
00547 
00548 
00549 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00550 VertexMesh<ELEMENT_DIM, SPACE_DIM>::~VertexMesh()
00551 {
00552     Clear();
00553 }
00554 
00555 
00556 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00557 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::SolveNodeMapping(unsigned index) const
00558 {
00559     assert(index < this->mNodes.size());
00560     return index;
00561 }
00562 
00563 
00564 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00565 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::SolveElementMapping(unsigned index) const
00566 {
00567     assert(index < this->mElements.size());
00568     return index;
00569 }
00570 
00571 
00572 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00573 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::SolveBoundaryElementMapping(unsigned index) const
00574 {
00576 //    assert(index < this->mBoundaryElements.size() );
00577     return index;
00578 }
00579 
00580 
00581 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00582 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetDelaunayNodeIndexCorrespondingToVoronoiElementIndex(unsigned elementIndex)
00583 {
00584     unsigned node_index = UNSIGNED_UNSET;
00585 
00586     if (mVoronoiElementIndexMap.empty())
00587     {
00588         node_index = elementIndex;
00589     }
00590     else
00591     {
00592         for (std::map<unsigned, unsigned>::iterator iter = mVoronoiElementIndexMap.begin();
00593              iter != mVoronoiElementIndexMap.end();
00594              ++iter)
00595         {
00596             if (iter->second == elementIndex)
00597             {
00598                 node_index = iter->first;
00599                 break;
00600             }
00601         }
00602     }
00603     assert(node_index != UNSIGNED_UNSET);
00604     return node_index;
00605 }
00606 
00607 
00608 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00609 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetVoronoiElementIndexCorrespondingToDelaunayNodeIndex(unsigned nodeIndex)
00610 {
00611     unsigned element_index = UNSIGNED_UNSET;
00612 
00613     if (mVoronoiElementIndexMap.empty())
00614     {
00615         element_index = nodeIndex;
00616     }
00617     else
00618     {
00619         std::map<unsigned, unsigned>::iterator iter = mVoronoiElementIndexMap.find(nodeIndex);
00620 
00621         if (iter == mVoronoiElementIndexMap.end())
00622         {
00623             EXCEPTION("This index does not correspond to a VertexElement");
00624         }
00625         else
00626         {
00627             element_index = iter->second;
00628         }
00629     }
00630     assert(element_index != UNSIGNED_UNSET);
00631     return element_index;
00632 }
00633 
00634 
00635 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00636 void VertexMesh<ELEMENT_DIM, SPACE_DIM>::Clear()
00637 {
00638     // Delete elements
00639     for (unsigned i=0; i<mElements.size(); i++)
00640     {
00641         delete mElements[i];
00642     }
00643     mElements.clear();
00644 
00645     // Delete faces
00646     for (unsigned i=0; i<mFaces.size(); i++)
00647     {
00648         delete mFaces[i];
00649     }
00650     mFaces.clear();
00651 
00652     // Delete nodes
00653     for (unsigned i=0; i<this->mNodes.size(); i++)
00654     {
00655         delete this->mNodes[i];
00656     }
00657     this->mNodes.clear();
00658 }
00659 
00660 
00661 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00662 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNumNodes() const
00663 {
00664     return this->mNodes.size();
00665 }
00666 
00667 
00668 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00669 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNumElements() const
00670 {
00671     return mElements.size();
00672 }
00673 
00674 
00675 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00676 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNumAllElements() const
00677 {
00678     return mElements.size();
00679 }
00680 
00681 
00682 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00683 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNumFaces() const
00684 {
00685     return mFaces.size();
00686 }
00687 
00688 
00689 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00690 VertexElement<ELEMENT_DIM, SPACE_DIM>* VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetElement(unsigned index) const
00691 {
00692     assert(index < mElements.size());
00693     return mElements[index];
00694 }
00695 
00696 
00697 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00698 VertexElement<ELEMENT_DIM-1, SPACE_DIM>* VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetFace(unsigned index) const
00699 {
00700     assert(index < mFaces.size());
00701     return mFaces[index];
00702 }
00703 
00704 
00705 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00706 c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetCentroidOfElement(unsigned index)
00707 {
00708     VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
00709     unsigned num_nodes_in_element = p_element->GetNumNodes();
00710 
00711     c_vector<double, SPACE_DIM> centroid = zero_vector<double>(SPACE_DIM);
00712 
00713     switch (SPACE_DIM)
00714     {
00715         case 1:
00716         {
00717             centroid = 0.5*(p_element->GetNodeLocation(0) + p_element->GetNodeLocation(1));
00718         }
00719         break;
00720         case 2: 
00721         {
00722             c_vector<double, SPACE_DIM> current_node;
00723             c_vector<double, SPACE_DIM> anticlockwise_node;
00724 
00725             double temp_centroid_x = 0;
00726             double temp_centroid_y = 0;
00727 
00728             for (unsigned local_index=0; local_index<num_nodes_in_element; local_index++)
00729             {
00730                 // Find locations of current node and anticlockwise node
00731                 current_node = p_element->GetNodeLocation(local_index);
00732                 anticlockwise_node = p_element->GetNodeLocation((local_index+1)%num_nodes_in_element);
00733 
00734                 temp_centroid_x += (current_node[0]+anticlockwise_node[0])*(current_node[0]*anticlockwise_node[1]-current_node[1]*anticlockwise_node[0]);
00735                 temp_centroid_y += (current_node[1]+anticlockwise_node[1])*(current_node[0]*anticlockwise_node[1]-current_node[1]*anticlockwise_node[0]);
00736             }
00737 
00738             double vertex_area = GetVolumeOfElement(index);
00739             double centroid_coefficient = 1.0/(6.0*vertex_area);
00740 
00741             centroid(0) = centroid_coefficient*temp_centroid_x;
00742             centroid(1) = centroid_coefficient*temp_centroid_y;
00743         }
00744         break;
00745         case 3:
00746         {
00747             for (unsigned local_index=0; local_index<num_nodes_in_element; local_index++)
00748             {
00749                 centroid += p_element->GetNodeLocation(local_index);
00750             }
00751             centroid /= ((double) num_nodes_in_element);
00752         }
00753         break;
00754         default:
00755             NEVER_REACHED;
00756     }
00757     return centroid;
00758 }
00759 
00760 
00761 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00762 std::set<unsigned> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNeighbouringNodeIndices(unsigned nodeIndex)
00763 {
00764     // Create a set of neighbouring node indices
00765     std::set<unsigned> neighbouring_node_indices;
00766 
00767     // Find the indices of the elements owned by this node
00768     std::set<unsigned> containing_elem_indices = this->GetNode(nodeIndex)->rGetContainingElementIndices();
00769 
00770     // Iterate over these elements
00771     for (std::set<unsigned>::iterator elem_iter = containing_elem_indices.begin();
00772          elem_iter != containing_elem_indices.end();
00773          ++elem_iter)
00774     {
00775         // Find the local index of this node in this element
00776         unsigned local_index = GetElement(*elem_iter)->GetNodeLocalIndex(nodeIndex);
00777 
00778         // Find the global indices of the preceding and successive nodes in this element
00779         unsigned num_nodes = GetElement(*elem_iter)->GetNumNodes();
00780         unsigned previous_local_index = (local_index + num_nodes - 1)%num_nodes;
00781         unsigned next_local_index = (local_index + 1)%num_nodes;
00782 
00783         // Add the global indices of these two nodes to the set of neighbouring node indices
00784         neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(previous_local_index));
00785         neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(next_local_index));
00786     }
00787 
00788     return neighbouring_node_indices;
00789 }
00790 
00791 
00792 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00793 std::set<unsigned> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNeighbouringNodeNotAlsoInElement(unsigned nodeIndex, unsigned elemIndex)
00794 {
00795     // Get a pointer to this element
00796     VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elemIndex);
00797 
00798     // Get the indices of the nodes contained in this element
00799     std::set<unsigned> node_indices_in_this_element;
00800     for (unsigned local_index=0; local_index<p_element->GetNumNodes(); local_index++)
00801     {
00802         unsigned global_index = p_element->GetNodeGlobalIndex(local_index);
00803         node_indices_in_this_element.insert(global_index);
00804     }
00805 
00806     // Check that the node is in fact contained in the element
00807     if (node_indices_in_this_element.find(nodeIndex) == node_indices_in_this_element.end())
00808     {
00809         EXCEPTION("The given node is not contained in the given element.");
00810     }
00811 
00812     // Create a set of neighbouring node indices
00813     std::set<unsigned> neighbouring_node_indices_not_in_this_element;
00814 
00815     // Get the indices of this node's neighbours
00816     std::set<unsigned> node_neighbours = GetNeighbouringNodeIndices(nodeIndex);
00817 
00818     // Check if each neighbour is also in this element; if not, add it to the set
00819     for (std::set<unsigned>::iterator iter = node_neighbours.begin();
00820          iter != node_neighbours.end();
00821          ++iter)
00822     {
00823         if (node_indices_in_this_element.find(*iter) == node_indices_in_this_element.end())
00824         {
00825             neighbouring_node_indices_not_in_this_element.insert(*iter);
00826         }
00827     }
00828 
00829     return neighbouring_node_indices_not_in_this_element;
00830 }
00831 
00832 
00833 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00834 std::set<unsigned> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNeighbouringElementIndices(unsigned elementIndex)
00835 {
00836     // Get a pointer to this element
00837     VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = this->GetElement(elementIndex);
00838 
00839     // Create a set of neighbouring element indices
00840     std::set<unsigned> neighbouring_element_indices;
00841 
00842     // Loop over nodes owned by this element
00843     for (unsigned local_index=0; local_index<p_element->GetNumNodes(); local_index++)
00844     {
00845         // Get a pointer to this node
00846         Node<SPACE_DIM>* p_node = p_element->GetNode(local_index);
00847 
00848         // Find the indices of the elements owned by this node
00849         std::set<unsigned> containing_elem_indices = p_node->rGetContainingElementIndices();
00850 
00851         // Form the union of this set with the current set of neighbouring element indices
00852         std::set<unsigned> all_elements;
00853         std::set_union(neighbouring_element_indices.begin(), neighbouring_element_indices.end(),
00854                        containing_elem_indices.begin(), containing_elem_indices.end(),
00855                        std::inserter(all_elements, all_elements.begin()));
00856 
00857         // Update the set of neighbouring element indices
00858         neighbouring_element_indices = all_elements;
00859     }
00860 
00861     // Lastly remove this element's index from the set of neighbouring element indices
00862     neighbouring_element_indices.erase(elementIndex);
00863 
00864     return neighbouring_element_indices;
00865 }
00866 
00867 
00869 template<>
00870 void VertexMesh<1,1>::ConstructFromMeshReader(AbstractMeshReader<1,1>& rMeshReader)
00872 {
00873 }
00874 
00876 template<>
00877 void VertexMesh<1,2>::ConstructFromMeshReader(AbstractMeshReader<1,2>& rMeshReader)
00879 {
00880 }
00881 
00883 template<>
00884 void VertexMesh<1,3>::ConstructFromMeshReader(AbstractMeshReader<1,3>& rMeshReader)
00886 {
00887 }
00888 
00890 template<>
00891 void VertexMesh<2,3>::ConstructFromMeshReader(AbstractMeshReader<2,3>& rMeshReader)
00893 {
00894 }
00895 
00897 template<>
00898 void VertexMesh<2,2>::ConstructFromMeshReader(AbstractMeshReader<2,2>& rMeshReader)
00900 {
00901     // Store numbers of nodes and elements
00902     unsigned num_nodes = rMeshReader.GetNumNodes();
00903     unsigned num_elements = rMeshReader.GetNumElements();
00904 
00905     // Reserve memory for nodes
00906     this->mNodes.reserve(num_nodes);
00907 
00908     rMeshReader.Reset();
00909 
00910     // Add nodes
00911     std::vector<double> node_data;
00912     for (unsigned i=0; i<num_nodes; i++)
00913     {
00914         node_data = rMeshReader.GetNextNode();
00915         unsigned is_boundary_node = (unsigned) node_data[2];
00916         node_data.pop_back();
00917         this->mNodes.push_back(new Node<2>(i, node_data, is_boundary_node));
00918     }
00919 
00920     rMeshReader.Reset();
00921 
00922     // Reserve memory for nodes
00923     mElements.reserve(rMeshReader.GetNumElements());
00924 
00925     // Add elements
00926     for (unsigned elem_index=0; elem_index<num_elements; elem_index++)
00927     {
00928         // Get the data for this element
00929         ElementData element_data = rMeshReader.GetNextElementData();
00930 
00931         // Get the nodes owned by this element
00932         std::vector<Node<2>*> nodes;
00933         unsigned num_nodes_in_element = element_data.NodeIndices.size();
00934         for (unsigned j=0; j<num_nodes_in_element; j++)
00935         {
00936             assert(element_data.NodeIndices[j] < this->mNodes.size());
00937             nodes.push_back(this->mNodes[element_data.NodeIndices[j]]);
00938         }
00939 
00940         // Use nodes and index to construct this element
00941         VertexElement<2,2>* p_element = new VertexElement<2,2>(elem_index, nodes);
00942         mElements.push_back(p_element);
00943 
00944         if (rMeshReader.GetNumElementAttributes() > 0)
00945         {
00946             assert(rMeshReader.GetNumElementAttributes() == 1);
00947             unsigned attribute_value = element_data.AttributeValue;
00948             p_element->SetRegion(attribute_value);
00949         }
00950     }
00951 }
00952 
00954 template<>
00955 void VertexMesh<3,3>::ConstructFromMeshReader(AbstractMeshReader<3,3>& rMeshReader)
00957 {
00958     // Store numbers of nodes and elements
00959     unsigned num_nodes = rMeshReader.GetNumNodes();
00960     unsigned num_elements = rMeshReader.GetNumElements();
00961 
00962     // Reserve memory for nodes
00963     this->mNodes.reserve(num_nodes);
00964 
00965     rMeshReader.Reset();
00966 
00967     // Add nodes
00968     std::vector<double> node_data;
00969     for (unsigned i=0; i<num_nodes; i++)
00970     {
00971         node_data = rMeshReader.GetNextNode();
00972         unsigned is_boundary_node = (unsigned) node_data[3];
00973         node_data.pop_back();
00974         this->mNodes.push_back(new Node<3>(i, node_data, is_boundary_node));
00975     }
00976 
00977     rMeshReader.Reset();
00978 
00979     // Reserve memory for nodes
00980     mElements.reserve(rMeshReader.GetNumElements());
00981 
00982     // Use a std::set to keep track of which faces have been added to mFaces
00983     std::set<unsigned> faces_counted;
00984 
00985     // Add elements
00986     for (unsigned elem_index=0; elem_index<num_elements; elem_index++)
00987     {
00989         typedef VertexMeshReader<3,3> VERTEX_MESH_READER;
00990         assert(dynamic_cast<VERTEX_MESH_READER*>(&rMeshReader) != NULL);
00991 
00992         // Get the data for this element
00993         VertexElementData element_data = static_cast<VERTEX_MESH_READER*>(&rMeshReader)->GetNextElementDataWithFaces();
00994 
00995         // Get the nodes owned by this element
00996         std::vector<Node<3>*> nodes;
00997         unsigned num_nodes_in_element = element_data.NodeIndices.size();
00998         for (unsigned j=0; j<num_nodes_in_element; j++)
00999         {
01000             assert(element_data.NodeIndices[j] < this->mNodes.size());
01001             nodes.push_back(this->mNodes[element_data.NodeIndices[j]]);
01002         }
01003 
01004         // Get the faces owned by this element
01005         std::vector<VertexElement<2,3>*> faces;
01006         unsigned num_faces_in_element = element_data.Faces.size();
01007         for (unsigned i=0; i<num_faces_in_element; i++)
01008         {
01009             // Get the data for this face
01010             ElementData face_data = element_data.Faces[i];
01011 
01012             // Get the face index
01013             unsigned face_index = face_data.AttributeValue;
01014 
01015             // Get the nodes owned by this face
01016             std::vector<Node<3>*> nodes_in_face;
01017             unsigned num_nodes_in_face = face_data.NodeIndices.size();
01018             for (unsigned j=0; j<num_nodes_in_face; j++)
01019             {
01020                 assert(face_data.NodeIndices[j] < this->mNodes.size());
01021                 nodes_in_face.push_back(this->mNodes[face_data.NodeIndices[j]]);
01022             }
01023 
01024             // If this face index is not already encountered, create a new face and update faces_counted...
01025             if (faces_counted.find(face_index) == faces_counted.end())
01026             {
01027                 // Use nodes and index to construct this face
01028                 VertexElement<2,3>* p_face = new VertexElement<2,3>(face_index, nodes_in_face);
01029                 mFaces.push_back(p_face);
01030                 faces_counted.insert(face_index);
01031                 faces.push_back(p_face);
01032             }
01033             else
01034             {
01035                 // ... otherwise use the member of mFaces with this index
01036                 bool face_added = false;
01037                 for (unsigned k=0; k<mFaces.size(); k++)
01038                 {
01039                     if (mFaces[k]->GetIndex() == face_index)
01040                     {
01041                         faces.push_back(mFaces[k]);
01042                         face_added = true;
01043                         break;
01044                     }
01045                 }
01046                 assert(face_added == true);
01047             }
01048         }
01049 
01051         std::vector<bool> orientations = std::vector<bool>(num_faces_in_element, true);
01052 
01053         // Use faces and index to construct this element
01054         VertexElement<3,3>* p_element = new VertexElement<3,3>(elem_index, faces, orientations, nodes);
01055         mElements.push_back(p_element);
01056 
01057         if (rMeshReader.GetNumElementAttributes() > 0)
01058         {
01059             assert(rMeshReader.GetNumElementAttributes() == 1);
01060             unsigned attribute_value = element_data.AttributeValue;
01061             p_element->SetRegion(attribute_value);
01062         }
01063     }
01064 }
01065 
01066 
01067 template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01068 c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetVectorFromAtoB(
01069     const c_vector<double, SPACE_DIM>& rLocationA, const c_vector<double, SPACE_DIM>& rLocationB)
01070 {
01071     c_vector<double, SPACE_DIM> vector;
01072     if (mpDelaunayMesh)
01073     {
01074         vector = mpDelaunayMesh->GetVectorFromAtoB(rLocationA, rLocationB);
01075     }
01076     else
01077     {
01078         vector = AbstractMesh<ELEMENT_DIM, SPACE_DIM>::GetVectorFromAtoB(rLocationA, rLocationB);
01079     }
01080     return vector;
01081 }
01082 
01083 
01084 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01085 double VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetVolumeOfElement(unsigned index)
01086 {
01087     assert(SPACE_DIM == 2 || SPACE_DIM == 3);
01088 
01089     // Get pointer to this element
01090     VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
01091 
01092     double element_volume = 0.0;
01093     if (SPACE_DIM == 2)
01094     {
01095         c_vector<double, SPACE_DIM> first_node = p_element->GetNodeLocation(0);
01096 
01097         unsigned num_nodes_in_element = p_element->GetNumNodes();
01098 
01099         for (unsigned local_index=0; local_index<num_nodes_in_element; local_index++)
01100         {
01101             // Find locations of current node and anticlockwise node
01102             c_vector<double, SPACE_DIM> current_node = p_element->GetNodeLocation(local_index);
01103             c_vector<double, SPACE_DIM> anticlockwise_node = p_element->GetNodeLocation((local_index+1)%num_nodes_in_element);
01104 
01105             /*
01106              * In order to calculate the area we map the origin to (x[0],y[0])
01107              * then use GetVectorFromAtoB() to get node cooordiantes
01108              */
01109             c_vector<double, SPACE_DIM> transformed_current_node = GetVectorFromAtoB(first_node, current_node);
01110             c_vector<double, SPACE_DIM> transformed_anticlockwise_node = GetVectorFromAtoB(first_node, anticlockwise_node);
01111 
01112             element_volume += 0.5*(transformed_current_node[0]*transformed_anticlockwise_node[1]
01113                                    - transformed_anticlockwise_node[0]*transformed_current_node[1]);
01114         }
01115     }
01116     else
01117     {
01118         // Loop over faces and add up pyramid volumes
01119         c_vector<double, SPACE_DIM> pyramid_apex = p_element->GetNodeLocation(0);
01120         for (unsigned face_index=0; face_index<p_element->GetNumFaces(); face_index++)
01121         {
01122             // Get pointer to face
01123             VertexElement<ELEMENT_DIM-1, SPACE_DIM>* p_face = p_element->GetFace(face_index);
01124 
01125             // Get unit normal to this face
01126             c_vector<double, SPACE_DIM> unit_normal = GetUnitNormalToFace(p_face);
01127 
01128             // Calculate the perpendicular distance from the plane of the face to the chosen apex
01129             c_vector<double, SPACE_DIM> base_to_apex = GetVectorFromAtoB(p_face->GetNodeLocation(0), pyramid_apex);
01130             double perpendicular_distance = inner_prod(base_to_apex, unit_normal);
01131 
01132             // Calculate the area of the face
01133             double face_area = GetAreaOfFace(p_face);
01134 
01135             // Use these to calculate the volume of the pyramid formed by the face and the point pyramid_apex
01136             element_volume += face_area * perpendicular_distance / 3;
01137         }
01138     }
01139     return fabs(element_volume);
01140 }
01141 
01142 
01143 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01144 double VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetSurfaceAreaOfElement(unsigned index)
01145 {
01146     assert(SPACE_DIM == 2 || SPACE_DIM == 3);
01147 
01148     // Get pointer to this element
01149     VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
01150 
01151     double surface_area = 0.0;
01152     if (SPACE_DIM == 2)
01153     {
01154         unsigned num_nodes_in_element = p_element->GetNumNodes();
01155         for (unsigned local_index=0; local_index<num_nodes_in_element; local_index++)
01156         {
01157             // Find locations of current node and anticlockwise node
01158             unsigned current_node_index = p_element->GetNodeGlobalIndex(local_index);
01159             unsigned anticlockwise_node_index = p_element->GetNodeGlobalIndex((local_index+1)%num_nodes_in_element);
01160 
01161             surface_area += this->GetDistanceBetweenNodes(current_node_index, anticlockwise_node_index);
01162         }
01163     }
01164     else
01165     {
01166         // Loop over faces and add up areas
01167         for (unsigned face_index=0; face_index<p_element->GetNumFaces(); face_index++)
01168         {
01169             surface_area += GetAreaOfFace(p_element->GetFace(face_index));
01170         }
01171     }
01172     return surface_area;
01173 }
01174 
01175 
01177 //                        2D-specific methods                       //
01179 
01180 
01181 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01182 bool VertexMesh<ELEMENT_DIM, SPACE_DIM>::ElementIncludesPoint(const c_vector<double, SPACE_DIM>& rTestPoint, unsigned elementIndex)
01183 {
01184     // Make sure that we are in the correct dimension - this code will be eliminated at compile time
01185     assert(SPACE_DIM == 2); // only works in 2D at present
01186     assert(ELEMENT_DIM == SPACE_DIM);
01187 
01188     // Initialise boolean
01189     bool element_includes_point = false;
01190 
01191     // Get the element
01192     VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(elementIndex);
01193     unsigned num_nodes = p_element->GetNumNodes();
01194 
01195     // Remap the origin to the first vertex to allow alternative distance metrics to be used in subclasses
01196     c_vector<double, SPACE_DIM> first_vertex = p_element->GetNodeLocation(0);
01197 
01198     c_vector<double, SPACE_DIM> test_point =  GetVectorFromAtoB(first_vertex,rTestPoint);
01199 
01200     // Loop over edges of the element
01201     for (unsigned local_index=0; local_index<num_nodes; local_index++)
01202     {
01203         // Get the end points of this edge
01204         // Remap to the origin to allow alternative distance metrics to be used in subclasses
01205         c_vector<double, SPACE_DIM> vertexA = GetVectorFromAtoB(first_vertex, p_element->GetNodeLocation(local_index));
01206         c_vector<double, SPACE_DIM> vertexB = GetVectorFromAtoB(first_vertex, p_element->GetNodeLocation((local_index+1)%num_nodes));
01207 
01208 
01209         // Check if this edge crosses the ray running out horizontally (increasing x, fixed y) from the test point
01210         c_vector<double, SPACE_DIM> vector_a_to_point = GetVectorFromAtoB(vertexA, test_point);
01211         c_vector<double, SPACE_DIM> vector_b_to_point = GetVectorFromAtoB(vertexB, test_point);
01212         c_vector<double, SPACE_DIM> vector_a_to_b = GetVectorFromAtoB(vertexA, vertexB);
01213 
01214         // Pathological case - test point coincides with vertexA or vertexB
01215         if (    (norm_2(vector_a_to_point) < DBL_EPSILON)
01216              || (norm_2(vector_b_to_point) < DBL_EPSILON) )
01217         {
01218             return false;
01219         }
01220 
01221         // Pathological case - ray coincides with horizontal edge
01222         if ( (fabs(vector_a_to_b[1]) < DBL_EPSILON) &&
01223              (fabs(vector_a_to_point[1]) < DBL_EPSILON) &&
01224              (fabs(vector_b_to_point[1]) < DBL_EPSILON) )
01225         {
01226             if ( (vector_a_to_point[0]>0) != (vector_b_to_point[0]>0) )
01227             {
01228                 return false;
01229             }
01230         }
01231 
01232         // Non-pathological case
01233         // A and B on different sides of the line y = test_point[1]
01234         if ( (vertexA[1] > test_point[1]) != (vertexB[1] > test_point[1]) )
01235         {
01236             // intersection of y=test_point[1] and vector_a_to_b is on the Right of test_point
01237             if (test_point[0] < vertexA[0] + vector_a_to_b[0]*vector_a_to_point[1]/vector_a_to_b[1])
01238             {
01239                 element_includes_point = !element_includes_point;
01240             }
01241         }
01242     }
01243     return element_includes_point;
01244 }
01245 
01246 
01247 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01248 unsigned VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetLocalIndexForElementEdgeClosestToPoint(const c_vector<double, SPACE_DIM>& rTestPoint, unsigned elementIndex)
01249 {
01250     // Make sure that we are in the correct dimension - this code will be eliminated at compile time
01251     assert(SPACE_DIM == 2); // only works in 2D at present
01252     assert(ELEMENT_DIM == SPACE_DIM);
01253 
01254     // Get the element
01255     VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(elementIndex);
01256     unsigned num_nodes = p_element->GetNumNodes();
01257 
01258     double min_squared_normal_distance = DBL_MAX;
01259     unsigned min_distance_edge_index = UINT_MAX;
01260 
01261     // Loop over edges of the element
01262     for (unsigned local_index=0; local_index<num_nodes; local_index++)
01263     {
01264         // Get the end points of this edge
01265         c_vector<double, SPACE_DIM> vertexA = p_element->GetNodeLocation(local_index);
01266         c_vector<double, SPACE_DIM> vertexB = p_element->GetNodeLocation((local_index+1)%num_nodes);
01267 
01268         c_vector<double, SPACE_DIM> vector_a_to_point = this->GetVectorFromAtoB(vertexA, rTestPoint);
01269         c_vector<double, SPACE_DIM> vector_a_to_b = this->GetVectorFromAtoB(vertexA, vertexB);
01270 
01271         c_vector<double, SPACE_DIM> edge_ab_unit_vector = vector_a_to_b/norm_2(vector_a_to_b);
01272         double distance_parallel_to_edge = inner_prod(vector_a_to_point, edge_ab_unit_vector);
01273 
01274         double squared_distance_normal_to_edge = SmallPow(norm_2(vector_a_to_point), 2) - SmallPow(distance_parallel_to_edge, 2);
01275 
01276         //If the point lies almost bang on the supporting line of the edge, then snap to the line
01277         //This allows us to do floating point tie-breaks when line is exactly at a node
01278         if (squared_distance_normal_to_edge < DBL_EPSILON)
01279         {
01280             squared_distance_normal_to_edge = 0.0;
01281         }
01282         // Make sure node is within the confines of the edge and is the nearest edge to the node \this breks for convex elements
01283         if (squared_distance_normal_to_edge < min_squared_normal_distance &&
01284             distance_parallel_to_edge >=0 &&
01285             distance_parallel_to_edge <= norm_2(vector_a_to_b))
01286         {
01287             min_squared_normal_distance = squared_distance_normal_to_edge;
01288             min_distance_edge_index = local_index;
01289         }
01290     }
01291     return min_distance_edge_index;
01292 }
01293 
01294 
01295 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01296 c_vector<double, 3> VertexMesh<ELEMENT_DIM, SPACE_DIM>::CalculateMomentsOfElement(unsigned index)
01297 {
01298     assert(SPACE_DIM == 2);
01299 
01300     VertexElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
01301     unsigned num_nodes_in_element = p_element->GetNumNodes();
01302     c_vector<double, 2> centroid = GetCentroidOfElement(index);
01303 
01304     c_vector<double, 3> moments = zero_vector<double>(3);
01305 
01306     unsigned node_1;
01307     unsigned node_2;
01308 
01309     for (unsigned local_index=0; local_index<num_nodes_in_element; local_index++)
01310     {
01311         node_1 = local_index;
01312         node_2 = (local_index+1)%num_nodes_in_element;
01313 
01314         // Original position of nodes
01315         c_vector<double, 2> original_pos_1 = p_element->GetNodeLocation(node_1);
01316         c_vector<double, 2> original_pos_2 = p_element->GetNodeLocation(node_2);
01317 
01318         // Node position so centerd on origin
01319         c_vector<double, 2> pos_1 = this->GetVectorFromAtoB(centroid, original_pos_1);
01320         c_vector<double, 2> pos_2 = this->GetVectorFromAtoB(centroid, original_pos_2);
01321 
01322         /*
01323          * Note these formulae require the polygon to be centered on the origin
01324          */
01325         double a = pos_1(0)*pos_2(1)-pos_2(0)*pos_1(1);
01326 
01327         // Ixx
01328         moments(0) += (  pos_1(1)*pos_1(1)
01329                        + pos_1(1)*pos_2(1)
01330                        + pos_2(1)*pos_2(1) ) * a;
01331 
01332         // Iyy
01333         moments(1) += (  pos_1(0)*pos_1(0)
01334                        + pos_1(0)*pos_2(0)
01335                        + pos_2(0)*pos_2(0) ) * a;
01336 
01337         // Ixy
01338         moments(2) += (  pos_1(0)*pos_2(1)
01339                        + 2*pos_1(0)*pos_1(1)
01340                        + 2*pos_2(0)*pos_2(1)
01341                        + pos_2(0)*pos_1(1) ) * a;
01342     }
01343 
01344     moments(0) /= 12;
01345     moments(1) /= 12;
01346     moments(2) /= 24;
01347 
01348     return moments;
01349 }
01350 
01351 
01352 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01353 c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetShortAxisOfElement(unsigned index)
01354 {
01355     assert(SPACE_DIM == 2);
01356 
01357     c_vector<double, SPACE_DIM> short_axis = zero_vector<double>(SPACE_DIM);
01358     c_vector<double, 3> moments = CalculateMomentsOfElement(index);
01359 
01360     double largest_eigenvalue, discriminant;
01361 
01362     discriminant = sqrt((moments(0) - moments(1))*(moments(0) - moments(1)) + 4.0*moments(2)*moments(2));
01363     // This is always the largest eigenvalue as both eigenvalues are real as it is a
01364     // symmetric matrix
01365     largest_eigenvalue = (moments(0) + moments(1) + discriminant)*0.5;
01366     if (fabs(discriminant) < 1e-10)
01367     {
01368         // Return a random unit vector
01369         short_axis(0) = RandomNumberGenerator::Instance()->ranf();
01370         short_axis(1) = sqrt(1.0 - short_axis(0)*short_axis(0));
01371     }
01372     else
01373     {
01374         if (moments(2) == 0.0)
01375         {
01376             short_axis(0) = (moments(0) < moments(1)) ? 0.0 : 1.0;
01377             short_axis(1) = (moments(0) < moments(1)) ? 1.0 : 0.0;
01378         }
01379         else
01380         {
01381             short_axis(0) = 1.0;
01382             short_axis(1) = (moments(0) - largest_eigenvalue)/moments(2);
01383 
01384             double length_short_axis = norm_2(short_axis);
01385 
01386             short_axis /= length_short_axis;
01387         }
01388     }
01389     return short_axis;
01390 }
01391 
01392 
01393 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01394 c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetAreaGradientOfElementAtNode(VertexElement<ELEMENT_DIM,SPACE_DIM>* pElement, unsigned localIndex)
01395 {
01396     assert(SPACE_DIM==2);
01397 
01398     unsigned num_nodes_in_element = pElement->GetNumNodes();
01399     unsigned next_local_index = (localIndex+1)%num_nodes_in_element;
01400 
01401     // We add an extra num_nodes_in_element in the line below as otherwise this term can be negative, which breaks the % operator
01402     unsigned previous_local_index = (num_nodes_in_element+localIndex-1)%num_nodes_in_element;
01403 
01404     c_vector<double, SPACE_DIM> previous_node_location = pElement->GetNodeLocation(previous_local_index);
01405     c_vector<double, SPACE_DIM> next_node_location = pElement->GetNodeLocation(next_local_index);
01406     c_vector<double, SPACE_DIM> difference_vector = this->GetVectorFromAtoB(previous_node_location, next_node_location);
01407 
01408     c_vector<double, SPACE_DIM> area_gradient;
01409 
01410     area_gradient[0] = 0.5*difference_vector[1];
01411     area_gradient[1] = -0.5*difference_vector[0];
01412 
01413     return area_gradient;
01414 }
01415 
01416 
01417 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01418 c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetPreviousEdgeGradientOfElementAtNode(VertexElement<ELEMENT_DIM,SPACE_DIM>* pElement, unsigned localIndex)
01419 {
01420     assert(SPACE_DIM==2);
01421 
01422     unsigned num_nodes_in_element = pElement->GetNumNodes();
01423 
01424     // We add an extra localIndex-1 in the line below as otherwise this term can be negative, which breaks the % operator
01425     unsigned previous_local_index = (num_nodes_in_element+localIndex-1)%num_nodes_in_element;
01426 
01427     unsigned current_global_index = pElement->GetNodeGlobalIndex(localIndex);
01428     unsigned previous_global_index = pElement->GetNodeGlobalIndex(previous_local_index);
01429 
01430     double previous_edge_length = this->GetDistanceBetweenNodes(current_global_index, previous_global_index);
01431     assert(previous_edge_length > DBL_EPSILON);
01432 
01433     c_vector<double, SPACE_DIM> previous_edge_gradient = this->GetVectorFromAtoB(pElement->GetNodeLocation(previous_local_index), pElement->GetNodeLocation(localIndex))/previous_edge_length;
01434 
01435     return previous_edge_gradient;
01436 }
01437 
01438 
01439 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01440 c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetNextEdgeGradientOfElementAtNode(VertexElement<ELEMENT_DIM,SPACE_DIM>* pElement, unsigned localIndex)
01441 {
01442     assert(SPACE_DIM==2);
01443 
01444     unsigned next_local_index = (localIndex+1)%(pElement->GetNumNodes());
01445 
01446     unsigned current_global_index = pElement->GetNodeGlobalIndex(localIndex);
01447     unsigned next_global_index = pElement->GetNodeGlobalIndex(next_local_index);
01448 
01449     double next_edge_length = this->GetDistanceBetweenNodes(current_global_index, next_global_index);
01450     assert(next_edge_length > DBL_EPSILON);
01451 
01452     c_vector<double, SPACE_DIM> next_edge_gradient = this->GetVectorFromAtoB(pElement->GetNodeLocation(next_local_index), pElement->GetNodeLocation(localIndex))/next_edge_length;
01453 
01454     return next_edge_gradient;
01455 }
01456 
01457 
01458 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01459 c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetPerimeterGradientOfElementAtNode(VertexElement<ELEMENT_DIM,SPACE_DIM>* pElement, unsigned localIndex)
01460 {
01461     assert(SPACE_DIM==2);
01462 
01463     c_vector<double, SPACE_DIM> previous_edge_gradient = GetPreviousEdgeGradientOfElementAtNode(pElement, localIndex);
01464     c_vector<double, SPACE_DIM> next_edge_gradient = GetNextEdgeGradientOfElementAtNode(pElement, localIndex);
01465 
01466     return previous_edge_gradient + next_edge_gradient;
01467 }
01468 
01469 
01471 //                        3D-specific methods                       //
01473 
01474 
01475 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01476 c_vector<double, SPACE_DIM> VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetUnitNormalToFace(VertexElement<ELEMENT_DIM-1, SPACE_DIM>* pFace)
01477 {
01478     assert(SPACE_DIM == 3);
01479 
01480     // As we are in 3D, the face must have at least three vertices, so use its first three vertices
01481     c_vector<double, SPACE_DIM> v0 = pFace->GetNode(0)->rGetLocation();
01482     c_vector<double, SPACE_DIM> v1 = pFace->GetNode(1)->rGetLocation();
01483     c_vector<double, SPACE_DIM> v2 = pFace->GetNode(2)->rGetLocation();
01484 
01485     c_vector<double, SPACE_DIM> v1_minus_v0 = this->GetVectorFromAtoB(v0, v1);
01486     c_vector<double, SPACE_DIM> v2_minus_v0 = this->GetVectorFromAtoB(v0, v2);
01487 
01488     c_vector<double, SPACE_DIM> unit_normal = zero_vector<double>(SPACE_DIM);
01489     unit_normal(0) = v1_minus_v0(1)*v2_minus_v0(2) - v1_minus_v0(2)*v2_minus_v0(1);
01490     unit_normal(1) = v1_minus_v0(2)*v2_minus_v0(0) - v1_minus_v0(0)*v2_minus_v0(2);
01491     unit_normal(2) = v1_minus_v0(0)*v2_minus_v0(1) - v1_minus_v0(1)*v2_minus_v0(0);
01492 
01493     // Normalize the normal vector
01494     unit_normal /= norm_2(unit_normal);
01495 
01496     return unit_normal;
01497 }
01498 
01499 
01500 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
01501 double VertexMesh<ELEMENT_DIM, SPACE_DIM>::GetAreaOfFace(VertexElement<ELEMENT_DIM-1, SPACE_DIM>* pFace)
01502 {
01503     assert(SPACE_DIM == 3);
01504 
01505     // Get the unit normal to the plane of this face
01506     c_vector<double, SPACE_DIM> unit_normal = GetUnitNormalToFace(pFace);
01507 
01508     // Select the largest absolute coordinate to ignore for planar projection
01509     double abs_x = unit_normal[0]>0 ? unit_normal[0]>0 : -unit_normal[0];
01510     double abs_y = unit_normal[1]>0 ? unit_normal[1]>0 : -unit_normal[1];
01511     double abs_z = unit_normal[2]>0 ? unit_normal[2]>0 : -unit_normal[2];
01512 
01513     unsigned dim_to_ignore = 2; // ignore z coordinate
01514     double abs = abs_z;
01515 
01516     if (abs_x > abs_y)
01517     {
01518         if (abs_x > abs_z)
01519         {
01520             dim_to_ignore = 0; // ignore x coordinate
01521             abs = abs_x;
01522         }
01523     }
01524     else if (abs_y > abs_z)
01525     {
01526         dim_to_ignore = 1; // ignore y coordinate
01527         abs = abs_y;
01528     }
01529 
01530     // Compute area of the 2D projection
01531     c_vector<double, SPACE_DIM-1> current_vertex;
01532     c_vector<double, SPACE_DIM-1> anticlockwise_vertex;
01533 
01534     unsigned num_nodes_in_face = pFace->GetNumNodes();
01535 
01536     unsigned dim1 = dim_to_ignore==0 ? 1 : 0;
01537     unsigned dim2 = dim_to_ignore==2 ? 1 : 2;
01538 
01539     double face_area = 0.0;
01540     for (unsigned local_index=0; local_index<num_nodes_in_face; local_index++)
01541     {
01542         // Find locations of current vertex and anticlockwise vertex
01543         current_vertex[0] = pFace->GetNodeLocation(local_index, dim1);
01544         current_vertex[1] = pFace->GetNodeLocation(local_index, dim2);
01545         anticlockwise_vertex[0] = pFace->GetNodeLocation((local_index+1)%num_nodes_in_face, dim1);
01546         anticlockwise_vertex[1] = pFace->GetNodeLocation((local_index+1)%num_nodes_in_face, dim2);
01547 
01548         // It doesn't matter if the face is oriented clockwise or not, since we area only interested in the area
01549         face_area += 0.5*(current_vertex[0]*anticlockwise_vertex[1] - anticlockwise_vertex[0]*current_vertex[1]);
01550     }
01551 
01552     // Scale to get area before projection
01553     face_area /= abs;
01554     return fabs(face_area);
01555 }
01557 #if defined(__xlC__)
01558 template<>
01559 double VertexMesh<1,1>::GetAreaOfFace(VertexElement<0,1>* pFace)
01560 {
01561     NEVER_REACHED;
01562 }
01563 #endif
01564 
01565 
01567 // Explicit instantiation
01569 
01570 template class VertexMesh<1,1>;
01571 template class VertexMesh<1,2>;
01572 template class VertexMesh<1,3>;
01573 template class VertexMesh<2,2>;
01574 template class VertexMesh<2,3>;
01575 template class VertexMesh<3,3>;
01576 
01577 
01578 // Serialization for Boost >= 1.36
01579 #include "SerializationExportWrapperForCpp.hpp"
01580 EXPORT_TEMPLATE_CLASS_ALL_DIMS(VertexMesh)

Generated on Tue May 31 14:31:48 2011 for Chaste by  doxygen 1.5.5