QuadraticMesh.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 "QuadraticMesh.hpp"
00030 #include "OutputFileHandler.hpp"
00031 #include "TrianglesMeshReader.hpp"
00032 
00033 //Jonathan Shewchuk's triangle and Hang Si's tetgen
00034 #define REAL double
00035 #define VOID void
00036 #include "triangle.h"
00037 #include "tetgen.h"
00038 #undef REAL
00039 #undef VOID
00040 
00041 template<unsigned DIM>
00042 void QuadraticMesh<DIM>::CountAndCheckVertices()
00043 {
00044     // count the number of vertices, and also check all vertices come before the
00045     // rest of the nodes (as this is assumed in
00046     // NonlinearElasticitySolver<DIM>::AllocateMatrixMemory() )
00047     //
00048     mNumVertices = 0;
00049     bool vertices_mode = true;
00050     for (unsigned i=0; i<this->GetNumNodes(); i++)
00051     {
00052         bool is_internal = this->GetNode(i)->IsInternal();
00053         if (is_internal==false)
00054         {
00055             mNumVertices++;
00056         }
00057         if ((vertices_mode == false)  && (is_internal==false ) )
00058         {
00059             //Covered in the 1D case by special mesh file data
00060             EXCEPTION("The quadratic mesh doesn't appear to have all vertices before the rest of the nodes");
00061         }
00062         if ( (vertices_mode == true)  && (is_internal==true) )
00063         {
00064             //This is the switch from vertices to internal nodes
00065             vertices_mode = false;
00066         }
00067     }
00068 }
00069 
00070 template<unsigned DIM>
00071 QuadraticMesh<DIM>::QuadraticMesh(double spaceStep, double width, double height, double depth)
00072 {
00073     this->ConstructRegularSlabMesh(spaceStep, width, height, depth);
00074 }
00075 
00076 template<unsigned DIM>
00077 void QuadraticMesh<DIM>::ConstructLinearMesh(unsigned numElemX)
00078 {
00079     this->mNodes.resize(2*numElemX+1);
00080     mNumVertices = numElemX+1;
00081     
00082     // create the left-most node
00083     Node<DIM>* p_edge_node = new Node<DIM>(0, true, 0.0);
00084     this->mNodes[0] = p_edge_node; // create nodes
00085     this->mBoundaryNodes.push_back(p_edge_node);
00086     this->mBoundaryElements.push_back(new BoundaryElement<DIM-1,DIM>(0, p_edge_node) );
00087     
00088     for (unsigned element_index=0; element_index<numElemX; element_index++)
00089     {
00090         unsigned right_node_index = element_index+1;
00091         unsigned mid_node_index = mNumVertices + element_index;
00092         
00093         double x_value_right_x = right_node_index;
00094         double x_value_mid_node = x_value_right_x-0.5;
00095         
00096         bool is_boundary = (element_index+1==numElemX);
00097         Node<DIM>* p_right_node = new Node<DIM>(right_node_index, is_boundary, x_value_right_x); 
00098         Node<DIM>* p_mid_node   = new Node<DIM>(mid_node_index, false, x_value_mid_node);
00099 
00100         this->mNodes[right_node_index] = p_right_node;
00101         this->mNodes[mid_node_index] = p_mid_node;
00102         
00103         if (element_index+1==numElemX) // right boundary
00104         {
00105             this->mBoundaryNodes.push_back(p_right_node);
00106             this->mBoundaryElements.push_back(new BoundaryElement<DIM-1,DIM>(1, p_right_node) );
00107         }
00108 
00109         std::vector<Node<DIM>*> nodes;
00110         nodes.push_back(this->mNodes[right_node_index-1]);
00111         nodes.push_back(this->mNodes[right_node_index]);
00112         nodes.push_back(this->mNodes[mid_node_index]);
00113         this->mElements.push_back(new Element<DIM,DIM>(element_index, nodes) );
00114     }
00115     
00116     this->RefreshMesh();
00117 }
00118 
00119 
00120 
00121 
00122 template<unsigned DIM>
00123 void QuadraticMesh<DIM>::ConstructRectangularMesh(unsigned numElemX, unsigned numElemY, bool unused)
00124 {
00125     assert(DIM==2);
00126 
00127     assert(numElemX>0);
00128     assert(numElemY>0);
00129     assert(unused);
00130     
00131     this->mMeshIsLinear=false;
00132     unsigned num_nodes=(numElemX+1)*(numElemY+1);
00133     struct triangulateio mesher_input;
00134 
00135     this->InitialiseTriangulateIo(mesher_input);
00136     mesher_input.pointlist = (double *) malloc( num_nodes * DIM * sizeof(double));
00137     mesher_input.numberofpoints = num_nodes;
00138 
00139     unsigned new_index = 0;
00140     for (unsigned j=0; j<=numElemY; j++)
00141     {
00142         double y = j;
00143         for (unsigned i=0; i<=numElemX; i++)
00144         {
00145             double x = i;
00146 
00147             mesher_input.pointlist[DIM*new_index] = x;
00148             mesher_input.pointlist[DIM*new_index + 1] = y;
00149             new_index++;
00150         }
00151     }
00152 
00153     // Make structure for output
00154     struct triangulateio mesher_output;
00155 
00156     this->InitialiseTriangulateIo(mesher_output);
00157 
00158     // Library call
00159     triangulate((char*)"Qzeo2", &mesher_input, &mesher_output, NULL);
00160   
00161     assert(mesher_output.numberofcorners == (DIM+1)*(DIM+2)/2);//Nodes per element (including internals, one per edge)
00162   
00163     this->ImportFromMesher(mesher_output, mesher_output.numberoftriangles, mesher_output.trianglelist, mesher_output.numberofedges, mesher_output.edgelist, mesher_output.edgemarkerlist);
00164     
00165     CountAndCheckVertices();
00166     AddNodesToBoundaryElements(NULL);
00167 
00168     this->FreeTriangulateIo(mesher_input);
00169     this->FreeTriangulateIo(mesher_output);
00170 }
00171 
00172 
00173 
00174 template<unsigned DIM>
00175 void QuadraticMesh<DIM>::ConstructCuboid(unsigned numElemX, unsigned numElemY, unsigned numElemZ)
00176 {
00177     assert(DIM==3);
00178 
00179 
00180     assert(numElemX>0);
00181     assert(numElemY>0);
00182     assert(numElemZ>0);
00183     this->mMeshIsLinear=false;
00184     unsigned num_nodes=(numElemX+1)*(numElemY+1)*(numElemZ+1);
00185     
00186     struct tetgen::tetgenio mesher_input;
00187     mesher_input.pointlist =  new double[num_nodes * DIM];
00188     mesher_input.numberofpoints = num_nodes;
00189     unsigned new_index = 0;
00190     for (unsigned k=0; k<=numElemZ; k++)
00191     {
00192         double z = k;
00193         for (unsigned j=0; j<=numElemY; j++)
00194         {
00195             double y = j;
00196             for (unsigned i=0; i<=numElemX; i++)
00197             {
00198                 double x = i;
00199                 mesher_input.pointlist[DIM*new_index] = x;
00200                 mesher_input.pointlist[DIM*new_index + 1] = y;
00201                 mesher_input.pointlist[DIM*new_index + 2] = z;
00202                 new_index++;
00203             }
00204         }
00205     }
00206 
00207     
00208     // Library call
00209     struct tetgen::tetgenio mesher_output;
00210     tetgen::tetrahedralize((char*)"Qo2", &mesher_input, &mesher_output, NULL);
00211     
00212     
00213     
00214     assert(mesher_output.numberofcorners == (DIM+1)*(DIM+2)/2);//Nodes per element (including internals, one per edge)
00215 
00216     this->ImportFromMesher(mesher_output, mesher_output.numberoftetrahedra, mesher_output.tetrahedronlist, mesher_output.numberoftrifaces, mesher_output.trifacelist, NULL);
00217     
00218     CountAndCheckVertices();
00219     AddNodesToBoundaryElements(NULL);
00220 }
00221 
00222 
00223 template<unsigned DIM>
00224 unsigned QuadraticMesh<DIM>::GetNumVertices()
00225 {
00226     return mNumVertices;
00227 }
00228 
00229 
00230 
00231 template<unsigned DIM>
00232 void QuadraticMesh<DIM>::ConstructFromLinearMeshReader(AbstractMeshReader<DIM, DIM>& rMeshReader)
00233 {
00234     assert(DIM != 1);
00235     
00236     //Make a linear mesh
00237     TetrahedralMesh<DIM,DIM>::ConstructFromMeshReader(rMeshReader);
00238     
00239     NodeMap unused_map(this->GetNumNodes());
00240     
00241     if (DIM==2)  // In 2D, remesh using triangle via library calls
00242     {
00243         struct triangulateio mesher_input, mesher_output;
00244         this->InitialiseTriangulateIo(mesher_input);
00245         this->InitialiseTriangulateIo(mesher_output);
00246 
00247         mesher_input.numberoftriangles = this->GetNumElements();
00248         mesher_input.trianglelist = (int *) malloc(this->GetNumElements() * (DIM+1) * sizeof(int));
00249         this->ExportToMesher(unused_map, mesher_input, mesher_input.trianglelist);
00250 
00251         // Library call
00252         triangulate((char*)"Qzero2", &mesher_input, &mesher_output, NULL);
00253         
00254         this->ImportFromMesher(mesher_output, mesher_output.numberoftriangles, mesher_output.trianglelist, mesher_output.numberofedges, mesher_output.edgelist, mesher_output.edgemarkerlist);
00255         CountAndCheckVertices();
00256         AddNodesToBoundaryElements(NULL);
00257 
00258         //Tidy up triangle
00259         this->FreeTriangulateIo(mesher_input);
00260         this->FreeTriangulateIo(mesher_output);
00261     }
00262     else // in 3D, remesh using tetgen
00263     {
00264 
00265         struct tetgen::tetgenio mesher_input, mesher_output;
00266 
00267         mesher_input.numberoftetrahedra = this->GetNumElements();
00268         mesher_input.tetrahedronlist = new int[this->GetNumElements() * (DIM+1)];
00269         this->ExportToMesher(unused_map, mesher_input, mesher_input.tetrahedronlist);
00270 
00271         // Library call
00272         tetgen::tetrahedralize((char*)"Qzro2", &mesher_input, &mesher_output);
00273         
00274         this->ImportFromMesher(mesher_output, mesher_output.numberoftetrahedra, mesher_output.tetrahedronlist, mesher_output.numberoftrifaces, mesher_output.trifacelist, NULL);
00275         CountAndCheckVertices();
00276         AddNodesToBoundaryElements(NULL);
00277     }    
00278 }
00279 
00280 
00281 template<unsigned DIM>
00282 void QuadraticMesh<DIM>::ConstructFromMeshReader(AbstractMeshReader<DIM, DIM>& rAbsMeshReader)
00283 {
00284     TrianglesMeshReader<DIM, DIM>* p_mesh_reader=dynamic_cast<TrianglesMeshReader<DIM, DIM>*>(&rAbsMeshReader);
00285     assert(p_mesh_reader != NULL);
00286 
00287 
00288     if (p_mesh_reader->GetOrderOfElements() == 1)
00289     {
00290         EXCEPTION("Supplied mesh reader is reading a linear mesh into quadratic mesh.  Consider using ConstructFromLinearMeshReader");
00291     }
00292 
00293     TetrahedralMesh<DIM,DIM>::ConstructFromMeshReader(*p_mesh_reader);
00294     assert(this->GetNumBoundaryElements()>0);
00295 
00296 
00297     p_mesh_reader->Reset();
00298 
00299 
00300     // add the extra nodes (1 extra node in 1D, 3 in 2D, 6 in 3D) to the element
00301     // data.
00302     for (unsigned i=0; i<this->GetNumElements(); i++)
00303     {
00304         std::vector<unsigned> nodes = p_mesh_reader->GetNextElementData().NodeIndices;
00305         assert(nodes.size()==(DIM+1)*(DIM+2)/2);
00306         assert(this->GetElement(i)->GetNumNodes()==DIM+1); // element is initially linear
00307         // add extra nodes to make it a quad element
00308         for (unsigned j=DIM+1; j<(DIM+1)*(DIM+2)/2; j++)
00309         {
00310             this->GetElement(i)->AddNode( this->GetNode(nodes[j]) );
00311             this->GetNode(nodes[j])->AddElement(this->GetElement(i)->GetIndex());
00312             this->GetNode(nodes[j])->MarkAsInternal();
00313         }
00314     }
00315 
00316     CountAndCheckVertices();
00317 
00318     if (DIM > 1)
00319     {
00320         // if  OrderOfBoundaryElements is 2 it can read in the extra nodes for each boundary element, other have to compute them.
00321         if (p_mesh_reader->GetOrderOfBoundaryElements() == 2u)
00322         {
00323             p_mesh_reader->Reset();
00324             for (typename TetrahedralMesh<DIM,DIM>::BoundaryElementIterator iter
00325                   = this->GetBoundaryElementIteratorBegin();
00326                  iter != this->GetBoundaryElementIteratorEnd();
00327                  ++iter)
00328             {
00329                 std::vector<unsigned> nodes = p_mesh_reader->GetNextFaceData().NodeIndices;
00330 
00331                 assert((*iter)->GetNumNodes()==DIM); // so far just the vertices
00332                 assert(nodes.size()==DIM*(DIM+1)/2); // the reader should have got 6 nodes (3d) for each face
00333 
00334                 for (unsigned j=DIM; j<DIM*(DIM+1)/2; j++)
00335                 {
00336                     Node <DIM> *p_internal_node = this->GetNode(nodes[j]);
00337                     (*iter)->AddNode( p_internal_node );
00338                     // add node to the boundary node list
00339                     if (!p_internal_node->IsBoundaryNode())
00340                     {
00341                         p_internal_node->SetAsBoundaryNode();
00342                         this->mBoundaryNodes.push_back(p_internal_node);
00343                     }
00344                 }
00345             }
00346         }
00347         else
00348         {
00349             AddNodesToBoundaryElements(p_mesh_reader);
00350         }
00351     }
00352 
00353     // Check each boundary element has a quadratic number of nodes
00354 #ifndef NDEBUG
00355     unsigned expected_num_nodes = DIM*(DIM+1)/2;
00356     for (typename TetrahedralMesh<DIM,DIM>::BoundaryElementIterator iter
00357           = this->GetBoundaryElementIteratorBegin();
00358           iter != this->GetBoundaryElementIteratorEnd();
00359           ++iter)
00360     {
00361         assert((*iter)->GetNumNodes()==expected_num_nodes);
00362     }
00363 #endif
00364 }
00365 
00366 template<unsigned DIM>
00367 void QuadraticMesh<DIM>::AddNodesToBoundaryElements(TrianglesMeshReader<DIM,DIM>* pMeshReader)
00368  {
00369     // Loop over all boundary elements, find the equivalent face from all
00370     // the elements, and add the extra nodes to the boundary element
00371     bool boundary_element_file_has_containing_element_info=false;
00372 
00373     if (pMeshReader)
00374     {
00375         boundary_element_file_has_containing_element_info=pMeshReader->GetReadContainingElementOfBoundaryElement();
00376     }
00377 
00378     if (DIM>1)
00379     {
00380         if(boundary_element_file_has_containing_element_info)
00381         {
00382             pMeshReader->Reset();
00383         }
00384 
00385         for (typename TetrahedralMesh<DIM,DIM>::BoundaryElementIterator iter
00386                = this->GetBoundaryElementIteratorBegin();
00387              iter != this->GetBoundaryElementIteratorEnd();
00388              ++iter)
00389         {
00390 
00391             // collect the nodes of this boundary element in a set
00392             std::set<unsigned> boundary_element_node_indices;
00393             for (unsigned i=0; i<DIM; i++)
00394             {
00395                 boundary_element_node_indices.insert( (*iter)->GetNodeGlobalIndex(i) );
00396             }
00397 
00398             bool found_this_boundary_element = false;
00399 
00400             // Loop over elements surrounding this face, then loop over each face of that element, and see if it matches
00401             // this boundary element.
00402             // Note, if we know what elem it should be in (boundary_element_file_has_containing_element_info==true)
00403             // we will reset elem_index immediately (below)
00404             Node<DIM>* p_representative_node = (*iter)->GetNode(0);
00405             for (typename Node<DIM>::ContainingElementIterator element_iter = p_representative_node->ContainingElementsBegin();
00406                  element_iter != p_representative_node->ContainingElementsEnd();
00407                  ++element_iter)
00408             {     
00409                 unsigned elem_index = *element_iter;
00410                 
00411                 // we know what elem it should be in
00412                 if(boundary_element_file_has_containing_element_info)
00413                 {
00414                     elem_index = pMeshReader->GetNextFaceData().ContainingElement;
00415                 }
00416 
00417                 Element<DIM,DIM>* p_element = this->GetElement(elem_index);
00418 
00419                 // for each element, loop over faces (the opposites to a node)
00420                 for (unsigned face=0; face<DIM+1; face++)
00421                 {
00422                     // collect the node indices for this face
00423                     std::set<unsigned> node_indices;
00424                     for (unsigned local_node_index=0; local_node_index<DIM+1; local_node_index++)
00425                     {
00426                         if (local_node_index!=face)
00427                         {
00428                             node_indices.insert( p_element->GetNodeGlobalIndex(local_node_index) );
00429                         }
00430                     }
00431 
00432                     assert(node_indices.size()==DIM);
00433 
00434                     // see if this face matches the boundary element,
00435                     // and call AddExtraBoundaryNodes() if so
00436                     if (node_indices==boundary_element_node_indices)
00437                     {
00438                         AddExtraBoundaryNodes(*iter, p_element, face);
00439 
00440                         found_this_boundary_element = true;
00441                         break;
00442                     }
00443                 }
00444 
00445                 // if the containing element info was given, we should certainly have found the
00446                 // face first time.
00447                 if(boundary_element_file_has_containing_element_info && !found_this_boundary_element)
00448                 {
00449                     #define COVERAGE_IGNORE
00450                     std::cout << (*iter)->GetIndex() << " " <<  pMeshReader->GetNextFaceData().ContainingElement << "\n";
00451                     std::stringstream ss;
00452                     ss << "Boundary element " << (*iter)->GetIndex()
00453                        << "wasn't found in the containing element given for it "
00454                        << elem_index;
00455 
00456                     EXCEPTION(ss.str());
00457                     #undef COVERAGE_IGNORE
00458                 }
00459 
00460                 if (found_this_boundary_element)
00461                 {
00462                     break;
00463                 }
00464             }
00465 
00466             if (!found_this_boundary_element)
00467             {
00468                 #define COVERAGE_IGNORE
00469                 EXCEPTION("Unable to find a face of an element which matches one of the boundary elements");
00470                 #undef COVERAGE_IGNORE
00471             }
00472         }
00473     }
00474 }
00475 
00476 
00477 template<unsigned DIM>
00478 void QuadraticMesh<DIM>::AddNodeToBoundaryElement(BoundaryElement<DIM-1,DIM>* pBoundaryElement,
00479                                                   Element<DIM,DIM>* pElement,
00480                                                   unsigned internalNode)
00481 {
00482     assert(DIM>1);
00483     assert(internalNode >= DIM+1);
00484     assert(internalNode < (DIM+1)*(DIM+2)/2);
00485     Node<DIM>* p_internal_node = pElement->GetNode(internalNode);
00486 
00487     // add node to the boundary node list
00488     if (!p_internal_node->IsBoundaryNode())
00489     {
00490         p_internal_node->SetAsBoundaryNode();
00491         this->mBoundaryNodes.push_back(p_internal_node);
00492     }
00493 
00494     pBoundaryElement->AddNode( p_internal_node );
00495 }
00496 
00497 
00498 template<unsigned DIM>
00499 void QuadraticMesh<DIM>::AddExtraBoundaryNodes(BoundaryElement<DIM-1,DIM>* pBoundaryElement,
00500                                                Element<DIM,DIM>* pElement,
00501                                                unsigned nodeIndexOppositeToFace)
00502 {
00503     assert(DIM!=1);
00504     if (DIM==2)
00505     {
00506         assert(nodeIndexOppositeToFace<3);
00507         // the single internal node of the elements face will be numbered 'face+3'
00508         AddNodeToBoundaryElement(pBoundaryElement, pElement, nodeIndexOppositeToFace+3);
00509     }
00510     else
00511     {
00512         assert(DIM==3);
00513 
00514         unsigned b_elem_n0 = pBoundaryElement->GetNodeGlobalIndex(0);
00515         unsigned b_elem_n1 = pBoundaryElement->GetNodeGlobalIndex(1);
00516 
00517         unsigned offset;
00518         bool reverse;
00519 
00520         if (nodeIndexOppositeToFace==0)
00521         {
00522             // face opposite to node 0 = {1,2,3}, with corresponding internals {9,8,5}
00523             HelperMethod1(b_elem_n0, b_elem_n1, pElement, 1, 2, 3, offset, reverse);
00524             HelperMethod2(pBoundaryElement, pElement, 9, 8, 5, offset, reverse);
00525         }
00526         else if (nodeIndexOppositeToFace==1)
00527         {
00528             // face opposite to node 1 = {2,0,3}, with corresponding internals {7,9,6}
00529             HelperMethod1(b_elem_n0, b_elem_n1, pElement, 2, 0, 3, offset, reverse);
00530             HelperMethod2(pBoundaryElement, pElement, 7, 9, 6, offset, reverse);
00531         }
00532         else if (nodeIndexOppositeToFace==2)
00533         {
00534             // face opposite to node 2 = {0,1,3}, with corresponding internals {8,7,4}
00535             HelperMethod1(b_elem_n0, b_elem_n1, pElement, 0, 1, 3, offset, reverse);
00536             HelperMethod2(pBoundaryElement, pElement, 8, 7, 4, offset, reverse);
00537         }
00538         else
00539         {
00540             assert(nodeIndexOppositeToFace==3);
00541             // face opposite to node 3 = {0,1,2}, with corresponding internals {5,6,4}
00542             HelperMethod1(b_elem_n0, b_elem_n1, pElement, 0, 1, 2, offset, reverse);
00543             HelperMethod2(pBoundaryElement, pElement, 5, 6, 4, offset, reverse);
00544         }
00545     }
00546 }
00547 
00548 template<unsigned DIM>
00549 void QuadraticMesh<DIM>::WriteBoundaryElementFile(std::string directory, std::string fileName)
00550 {
00551     OutputFileHandler handler(directory, false);
00552     out_stream p_file = handler.OpenOutputFile(fileName);
00553 
00554     unsigned expected_num_nodes;
00555     assert(DIM > 1);
00556     if (DIM == 2)
00557     {
00558         expected_num_nodes = 3;
00559     }
00560     else if (DIM == 3)
00561     {
00562         expected_num_nodes = 6;
00563     }
00564 
00565     unsigned num_elements = 0;
00566 
00567     for (typename TetrahedralMesh<DIM,DIM>::BoundaryElementIterator iter
00568           = this->GetBoundaryElementIteratorBegin();
00569           iter != this->GetBoundaryElementIteratorEnd();
00570           ++iter)
00571     {
00572         assert((*iter)->GetNumNodes()==expected_num_nodes);
00573         num_elements++;
00574     }
00575 
00576     *p_file << num_elements << " 0\n";
00577 
00578     unsigned counter = 0;
00579     for (typename TetrahedralMesh<DIM,DIM>::BoundaryElementIterator iter
00580           = this->GetBoundaryElementIteratorBegin();
00581           iter != this->GetBoundaryElementIteratorEnd();
00582           ++iter)
00583     {
00584         *p_file << counter++ << " ";
00585         for (unsigned i=0; i<(*iter)->GetNumNodes(); i++)
00586         {
00587             *p_file << (*iter)->GetNodeGlobalIndex(i) << " ";
00588         }
00589         *p_file << "\n";
00590     }
00591 
00592     p_file->close();
00593 }
00594 
00596 // two unpleasant helper methods for AddExtraBoundaryNodes()
00598 
00599 #define COVERAGE_IGNORE 
00600 template<unsigned DIM>
00601 void QuadraticMesh<DIM>::HelperMethod1(unsigned boundaryElemNode0, unsigned boundaryElemNode1,
00602                                        Element<DIM,DIM>* pElement,
00603                                        unsigned node0, unsigned node1, unsigned node2,
00604                                        unsigned& rOffset,
00605                                        bool& rReverse)
00606 {
00607     if (pElement->GetNodeGlobalIndex(node0)==boundaryElemNode0)
00608     {
00609         rOffset = 0;
00610         if (pElement->GetNodeGlobalIndex(node1)==boundaryElemNode1)
00611         {
00612             rReverse = false;
00613         }
00614         else
00615         {
00616             assert(pElement->GetNodeGlobalIndex(node2)==boundaryElemNode1);
00617             rReverse = true;
00618         }
00619     }
00620     else if (pElement->GetNodeGlobalIndex(node1)==boundaryElemNode0)
00621     {
00622         rOffset = 1;
00623         if (pElement->GetNodeGlobalIndex(node2)==boundaryElemNode1)
00624         {
00625             rReverse = false;
00626         }
00627         else
00628         {
00629             assert(pElement->GetNodeGlobalIndex(node0)==boundaryElemNode1);
00630             rReverse = true;
00631         }
00632     }
00633     else
00634     {
00635         assert(pElement->GetNodeGlobalIndex(node2)==boundaryElemNode0);
00636         rOffset = 2;
00637         if (pElement->GetNodeGlobalIndex(node0)==boundaryElemNode1)
00638         {
00639             rReverse = false;
00640         }
00641         else
00642         {
00643             assert(pElement->GetNodeGlobalIndex(node1)==boundaryElemNode1);
00644             rReverse = true;
00645         }
00646     }
00647 }
00648 #undef COVERAGE_IGNORE 
00649 
00650 
00651 #define COVERAGE_IGNORE 
00652 template<unsigned DIM>
00653 void QuadraticMesh<DIM>::HelperMethod2(BoundaryElement<DIM-1,DIM>* pBoundaryElement,
00654                                        Element<DIM,DIM>* pElement,
00655                                        unsigned internalNode0, unsigned internalNode1, unsigned internalNode2,
00656                                        unsigned offset,
00657                                        bool reverse)
00658 {
00659     if (offset==1)
00660     {
00661         unsigned temp = internalNode0;
00662         internalNode0 = internalNode1;
00663         internalNode1 = internalNode2;
00664         internalNode2 = temp;
00665     }
00666     else if (offset == 2)
00667     {
00668         unsigned temp = internalNode0;
00669         internalNode0 = internalNode2;
00670         internalNode2 = internalNode1;
00671         internalNode1 = temp;
00672     }
00673 
00674     if (reverse)
00675     {
00676         unsigned temp = internalNode1;
00677         internalNode1 = internalNode2;
00678         internalNode2 = temp;
00679     }
00680 
00681     AddNodeToBoundaryElement(pBoundaryElement, pElement, internalNode0);
00682     AddNodeToBoundaryElement(pBoundaryElement, pElement, internalNode1);
00683     AddNodeToBoundaryElement(pBoundaryElement, pElement, internalNode2);
00684 }
00685 #undef COVERAGE_IGNORE 
00686 
00687 
00689 // Explicit instantiation
00691 
00692 
00693 template class QuadraticMesh<1>;
00694 template class QuadraticMesh<2>;
00695 template class QuadraticMesh<3>;
00696 
00697 
00698 // Serialization for Boost >= 1.36
00699 #include "SerializationExportWrapperForCpp.hpp"
00700 EXPORT_TEMPLATE_CLASS_SAME_DIMS(QuadraticMesh)

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