Cylindrical2dMesh.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 <map>
00030 
00031 /* These lines are very useful for debugging (visualize with 'showme').
00032 #include "TrianglesMeshWriter.hpp"
00033 TrianglesMeshWriter<2,2> mesh_writer("Cylindrical2dMeshDebug", "mesh", false);
00034 mesh_writer.WriteFilesUsingMesh(*this);
00035 */
00036 #include "Cylindrical2dMesh.hpp"
00037 
00038 Cylindrical2dMesh::Cylindrical2dMesh(double width)
00039   : MutableMesh<2,2>(),
00040     mWidth(width)
00041 {
00042     assert(width > 0.0);
00043 }
00044 
00045 Cylindrical2dMesh::~Cylindrical2dMesh()
00046 {
00047 }
00048 
00049 Cylindrical2dMesh::Cylindrical2dMesh(double width, std::vector<Node<2>* > nodes)
00050   : MutableMesh<2,2>(),
00051     mWidth(width)
00052 {
00053     assert(width > 0.0);
00054     for (unsigned index=0; index<nodes.size(); index++)
00055     {
00056         Node<2>* p_temp_node = nodes[index];
00057         double x = p_temp_node->rGetLocation()[0];
00058         x = x; // Fix optimised build
00059         assert( 0 <= x && x < width);
00060         mNodes.push_back(p_temp_node);
00061     }
00062 
00063     NodeMap node_map(nodes.size());
00064     ReMesh(node_map);
00065 }
00066 
00067 void Cylindrical2dMesh::UpdateTopAndBottom()
00068 {
00069     ChasteCuboid<2> bounding_box = CalculateBoundingBox();
00070     mBottom = bounding_box.rGetLowerCorner()[1];
00071     mTop = bounding_box.rGetUpperCorner()[1];
00072 }
00073 
00074 void Cylindrical2dMesh::CreateMirrorNodes()
00075 {
00076     double half_way = 0.5*mWidth;
00077 
00078     mLeftOriginals.clear();
00079     mLeftImages.clear();
00080     mImageToLeftOriginalNodeMap.clear();
00081     mRightOriginals.clear();
00082     mRightImages.clear();
00083     mImageToRightOriginalNodeMap.clear();
00084     mLeftPeriodicBoundaryElementIndices.clear();
00085     mRightPeriodicBoundaryElementIndices.clear();
00086 
00087     for (AbstractMesh<2,2>::NodeIterator node_iter = GetNodeIteratorBegin();
00088          node_iter != GetNodeIteratorEnd();
00089          ++node_iter)
00090     {
00091         c_vector<double, 2> location = node_iter->rGetLocation();
00092         unsigned this_node_index = node_iter->GetIndex();
00093         double this_node_x_location = location[0];
00094 
00095         // Check the mesh currently conforms to the dimensions given
00096         assert(0.0 <= location[0]);
00097         assert(location[0] <= mWidth);
00098 
00099         // Put the nodes which are to be mirrored in the relevant vectors
00100         if (this_node_x_location < half_way)
00101         {
00102             mLeftOriginals.push_back(this_node_index);
00103         }
00104         else
00105         {
00106             mRightOriginals.push_back(this_node_index);
00107         }
00108     }
00109 
00110     // For each left original node, create an image node and record its new index
00111     for (unsigned i=0; i<mLeftOriginals.size(); i++)
00112     {
00113         c_vector<double, 2> location = mNodes[mLeftOriginals[i]]->rGetLocation();
00114         location[0] = location[0] + mWidth;
00115 
00116         unsigned new_node_index = MutableMesh<2,2>::AddNode(new Node<2>(0, location));
00117         mLeftImages.push_back(new_node_index);
00118         mImageToLeftOriginalNodeMap[new_node_index] = mLeftOriginals[i];
00119     }
00120 
00121     // For each right original node, create an image node and record its new index
00122     for (unsigned i=0; i<mRightOriginals.size(); i++)
00123     {
00124         c_vector<double, 2> location = mNodes[mRightOriginals[i]]->rGetLocation();
00125         location[0] = location[0] - mWidth;
00126 
00127         unsigned new_node_index = MutableMesh<2,2>::AddNode(new Node<2>(0, location));
00128         mRightImages.push_back(new_node_index);
00129         mImageToRightOriginalNodeMap[new_node_index] = mRightOriginals[i];
00130     }
00131 
00132     assert(mRightOriginals.size() == mRightImages.size());
00133     assert(mLeftOriginals.size() == mLeftImages.size());
00134     assert(mImageToLeftOriginalNodeMap.size() == mLeftOriginals.size());
00135     assert(mImageToRightOriginalNodeMap.size() == mRightOriginals.size());
00136 }
00137 
00138 void Cylindrical2dMesh::CreateHaloNodes()
00139 {
00140     UpdateTopAndBottom();
00141 
00142     mTopHaloNodes.clear();
00143     mBottomHaloNodes.clear();
00144 
00145     unsigned num_halo_nodes = (unsigned)(floor(mWidth*2.0));
00146     double halo_node_separation = mWidth/((double)(num_halo_nodes));
00147     double y_top_coordinate = mTop + halo_node_separation;
00148     double y_bottom_coordinate = mBottom - halo_node_separation;
00149 
00150     c_vector<double, 2> location;
00151     for (unsigned i=0; i<num_halo_nodes; i++)
00152     {
00153        double x_coordinate = 0.5*halo_node_separation + (double)(i)*halo_node_separation;
00154 
00155        // Inserting top halo node in mesh
00156        location[0] = x_coordinate;
00157        location[1] = y_top_coordinate;
00158        unsigned new_node_index = MutableMesh<2,2>::AddNode(new Node<2>(0, location));
00159        mTopHaloNodes.push_back(new_node_index);
00160 
00161        location[1] = y_bottom_coordinate;
00162        new_node_index = MutableMesh<2,2>::AddNode(new Node<2>(0, location));
00163        mBottomHaloNodes.push_back(new_node_index);
00164     }
00165 }
00166 
00167 void Cylindrical2dMesh::ReMesh(NodeMap& rMap)
00168 {
00169     unsigned old_num_all_nodes = GetNumAllNodes();
00170 
00171     rMap.Resize(old_num_all_nodes);
00172     rMap.ResetToIdentity();
00173 
00174     // Flag the deleted nodes as deleted in the map
00175     for (unsigned i=0; i<old_num_all_nodes; i++)
00176     {
00177         if (mNodes[i]->IsDeleted())
00178         {
00179             rMap.SetDeleted(i);
00180         }
00181     }
00182 
00183     CreateHaloNodes();
00184 
00185     // Create mirrored nodes for the normal remesher to work with
00186     CreateMirrorNodes();
00187 
00188     /*
00189      * The mesh now has messed up boundary elements, but this
00190      * doesn't matter as the ReMesh() below doesn't read them in
00191      * and reconstructs the boundary elements.
00192      * 
00193      * Call ReMesh() on the parent class. Note that the mesh now has lots
00194      * of extra nodes which will be deleted, hence the name 'big_map'.
00195      */
00196     NodeMap big_map(GetNumAllNodes());
00197     MutableMesh<2,2>::ReMesh(big_map);
00198 
00199     /*
00200      * If the big_map isn't the identity map, the little map ('map') needs to be
00201      * altered accordingly before being passed to the user. Not sure how this all
00202      * works, so deal with this bridge when we get to it.
00203      */
00204     assert(big_map.IsIdentityMap());
00205 
00206     // Re-index the vectors according to the big nodemap, and set up the maps.
00207     mImageToLeftOriginalNodeMap.clear();
00208     mImageToRightOriginalNodeMap.clear();
00209 
00210     assert(mLeftOriginals.size() == mLeftImages.size());
00211     assert(mRightOriginals.size() == mRightImages.size());
00212 
00213     for (unsigned i=0; i<mLeftOriginals.size(); i++)
00214     {
00215         mLeftOriginals[i] = big_map.GetNewIndex(mLeftOriginals[i]);
00216         mLeftImages[i] = big_map.GetNewIndex(mLeftImages[i]);
00217         mImageToLeftOriginalNodeMap[mLeftImages[i]] = mLeftOriginals[i];
00218     }
00219 
00220     for (unsigned i=0; i<mRightOriginals.size(); i++)
00221     {
00222         mRightOriginals[i] = big_map.GetNewIndex(mRightOriginals[i]);
00223         mRightImages[i] = big_map.GetNewIndex(mRightImages[i]);
00224         mImageToRightOriginalNodeMap[mRightImages[i]] = mRightOriginals[i];
00225     }
00226 
00227     for (unsigned i=0; i<mTopHaloNodes.size(); i++)
00228     {
00229         mTopHaloNodes[i] = big_map.GetNewIndex(mTopHaloNodes[i]);
00230         mBottomHaloNodes[i] = big_map.GetNewIndex(mBottomHaloNodes[i]);
00231     }
00232 
00233     /*
00234      * Check that elements crossing the periodic boundary have been meshed
00235      * in the same way at each side.
00236      */
00237     CorrectNonPeriodicMesh();
00238 
00239     /*
00240      * Take the double-sized mesh, with its new boundary elements, and
00241      * remove the relevant nodes, elements and boundary elements to leave
00242      * a proper periodic mesh.
00243      */
00244     ReconstructCylindricalMesh();
00245 
00246     DeleteHaloNodes();
00247 
00248     /*
00249      * Create a random boundary element between two nodes of the first
00250      * element if it is not deleted. This is a temporary measure to get
00251      * around re-index crashing when there are no boundary elements.
00252      */
00253     unsigned num_elements = GetNumAllElements();
00254     bool boundary_element_made = false;
00255     unsigned elem_index = 0;
00256 
00257     while (elem_index<num_elements && !boundary_element_made)
00258     {
00259         Element<2,2>* p_element = GetElement(elem_index);
00260         if (!p_element->IsDeleted())
00261         {
00262             boundary_element_made = true;
00263             std::vector<Node<2>*> nodes;
00264             nodes.push_back(p_element->GetNode(0));
00265             nodes.push_back(p_element->GetNode(1));
00266             BoundaryElement<1,2>* p_boundary_element = new BoundaryElement<1,2>(0, nodes);
00267             p_boundary_element->RegisterWithNodes();
00268             mBoundaryElements.push_back(p_boundary_element);
00269             this->mBoundaryElementWeightedDirections.push_back(zero_vector<double>(2));
00270             this->mBoundaryElementJacobianDeterminants.push_back(0.0);
00271         }
00272         elem_index++;
00273     }
00274 
00275     // Now call ReIndex() to remove the temporary nodes which are marked as deleted
00276     NodeMap reindex_map(GetNumAllNodes());
00277     ReIndex(reindex_map);
00278     assert(!reindex_map.IsIdentityMap());  // maybe don't need this
00279 
00280     /*
00281      * Go through the reindex map and use it to populate the original NodeMap
00282      * (the one that is returned to the user).
00283      */
00284     for (unsigned i=0; i<rMap.Size(); i++) // only going up to be size of map, not size of reindex_map
00285     {
00286         if (reindex_map.IsDeleted(i))
00287         {
00288             /*
00289              * i < num_original_nodes and node is deleted, this should correspond to
00290              * a node that was labelled as before the remeshing, so should have already
00291              * been set as deleted in the map above.
00292              */
00293             assert(rMap.IsDeleted(i));
00294         }
00295         else
00296         {
00297             rMap.SetNewIndex(i, reindex_map.GetNewIndex(i));
00298         }
00299     }
00300 
00301     // We can now clear the index vectors and maps; they are only used for remeshing
00302     mLeftOriginals.clear();
00303     mLeftImages.clear();
00304     mImageToLeftOriginalNodeMap.clear();
00305     mRightOriginals.clear();
00306     mRightImages.clear();
00307     mImageToRightOriginalNodeMap.clear();
00308     mLeftPeriodicBoundaryElementIndices.clear();
00309     mRightPeriodicBoundaryElementIndices.clear();
00310 }
00311 
00312 void Cylindrical2dMesh::ReconstructCylindricalMesh()
00313 {
00314     /*
00315      * Figure out which elements have real nodes and image nodes in them
00316      * and replace image nodes with corresponding real ones.
00317      */
00318     for (MutableMesh<2,2>::ElementIterator elem_iter = GetElementIteratorBegin();
00319          elem_iter != GetElementIteratorEnd();
00320          ++elem_iter)
00321     {
00322         // Left images are on the right of the mesh
00323         unsigned number_of_left_image_nodes = 0;
00324         unsigned number_of_right_image_nodes = 0;
00325 
00326         for (unsigned i=0; i<3; i++)
00327         {
00328             unsigned this_node_index = elem_iter->GetNodeGlobalIndex(i);
00329 
00330             if (mImageToLeftOriginalNodeMap.find(this_node_index) != mImageToLeftOriginalNodeMap.end())
00331             {
00332                 number_of_left_image_nodes++;
00333             }
00334             else if (mImageToRightOriginalNodeMap.find(this_node_index) != mImageToRightOriginalNodeMap.end())
00335             {
00336                 number_of_right_image_nodes++;
00337             }
00338         }
00339 
00340         // Delete all the elements on the left hand side (images of right)...
00341         if (number_of_right_image_nodes >= 1)
00342         {
00343             elem_iter->MarkAsDeleted();
00344             mDeletedElementIndices.push_back(elem_iter->GetIndex());
00345         }
00346 
00347         // Delete only purely imaginary elements on the right (images of left nodes)
00348         if (number_of_left_image_nodes == 3)
00349         {
00350             elem_iter->MarkAsDeleted();
00351             mDeletedElementIndices.push_back(elem_iter->GetIndex());
00352         }
00353 
00354         /*
00355          * If some are images then replace them with the real nodes. There
00356          * can be elements with either two image nodes on the right (and one
00357          * real) or one image node on the right (and two real).
00358          */
00359         if (number_of_left_image_nodes == 1 || number_of_left_image_nodes == 2)
00360         {
00361             for (unsigned i=0; i<3; i++)
00362             {
00363                 unsigned this_node_index = elem_iter->GetNodeGlobalIndex(i);
00364                 std::map<unsigned, unsigned>::iterator it = mImageToLeftOriginalNodeMap.find(this_node_index);
00365                 if (it != mImageToLeftOriginalNodeMap.end())
00366                 {
00367                     elem_iter->ReplaceNode(mNodes[this_node_index], mNodes[it->second]);
00368                 }
00369             }
00370         }
00371     }
00372 
00373     /*
00374      * Figure out which boundary elements have real nodes and image nodes in them
00375      * and replace image nodes with corresponding real ones.
00376      */
00377     for (unsigned elem_index = 0; elem_index<GetNumAllBoundaryElements(); elem_index++)
00378     {
00379         BoundaryElement<1,2>* p_boundary_element = GetBoundaryElement(elem_index);
00380         if (!p_boundary_element->IsDeleted())
00381         {
00382             unsigned number_of_image_nodes = 0;
00383             for (unsigned i=0; i<2; i++)
00384             {
00385                 unsigned this_node_index = p_boundary_element->GetNodeGlobalIndex(i);
00386 
00387                 if (mImageToLeftOriginalNodeMap.find(this_node_index) != mImageToLeftOriginalNodeMap.end())
00388                 {
00389                     number_of_image_nodes++;
00390                 }
00391                 else if (mImageToRightOriginalNodeMap.find(this_node_index) != mImageToRightOriginalNodeMap.end())
00392                 {
00393                     number_of_image_nodes++;
00394                 }
00395             }
00396 
00397             if (number_of_image_nodes == 2)
00398             {
00399                 p_boundary_element->MarkAsDeleted();
00400                 mDeletedBoundaryElementIndices.push_back(p_boundary_element->GetIndex());
00401             }
00402 
00403             /*
00404              * To avoid having two copies of the boundary elements on the periodic
00405              * boundaries we only deal with the elements on the left image and
00406              * delete the ones on the right image.
00407              */
00408             if (number_of_image_nodes == 1)
00409             {
00410                 for (unsigned i=0; i<2; i++)
00411                 {
00412                     unsigned this_node_index = p_boundary_element->GetNodeGlobalIndex(i);
00413                     std::map<unsigned, unsigned>::iterator it = mImageToLeftOriginalNodeMap.find(this_node_index);
00414                     if (it != mImageToLeftOriginalNodeMap.end())
00415                     {
00416                         p_boundary_element->ReplaceNode(mNodes[this_node_index], mNodes[it->second]);
00417                     }
00418                     else
00419                     {
00420                         it = mImageToRightOriginalNodeMap.find(this_node_index);
00421                         if (it != mImageToRightOriginalNodeMap.end())
00422                         {
00423                             p_boundary_element->MarkAsDeleted();
00424                             mDeletedBoundaryElementIndices.push_back(p_boundary_element->GetIndex());
00425                         }
00426                     }
00427                 }
00428             }
00429         }
00430     }
00431 
00432     // Delete all image nodes unless they have already gone (halo nodes)
00433     for (unsigned i=0; i<mLeftImages.size(); i++)
00434     {
00435         mNodes[mLeftImages[i]]->MarkAsDeleted();
00436         mDeletedNodeIndices.push_back(mLeftImages[i]);
00437     }
00438 
00439     for (unsigned i=0; i<mRightImages.size(); i++)
00440     {
00441         mNodes[mRightImages[i]]->MarkAsDeleted();
00442         mDeletedNodeIndices.push_back(mRightImages[i]);
00443     }
00444 }
00445 
00446 void Cylindrical2dMesh::DeleteHaloNodes()
00447 {
00448     assert(mTopHaloNodes.size() == mBottomHaloNodes.size());
00449     for (unsigned i=0; i<mTopHaloNodes.size(); i++)
00450     {
00451         DeleteBoundaryNodeAt(mTopHaloNodes[i]);
00452         DeleteBoundaryNodeAt(mBottomHaloNodes[i]);
00453     }
00454 }
00455 
00456 c_vector<double, 2> Cylindrical2dMesh::GetVectorFromAtoB(const c_vector<double, 2>& rLocation1, const c_vector<double, 2>& rLocation2)
00457 {
00458     assert(mWidth > 0.0);
00459 
00460     c_vector<double, 2> location1 = rLocation1;
00461     c_vector<double, 2> location2 = rLocation2;
00462 
00463     location1[0] = fmod(location1[0], mWidth);
00464     location2[0] = fmod(location2[0], mWidth);
00465 
00466     c_vector<double, 2> vector = location2 - location1;
00467 
00468     /*
00469      * Handle the cylindrical condition here: if the points are more
00470      * than halfway around the cylinder apart, measure the other way.
00471      */
00472     if (vector[0] > 0.5*mWidth)
00473     {
00474         vector[0] -= mWidth;
00475     }
00476     else if (vector[0] < -0.5*mWidth)
00477     {
00478         vector[0] += mWidth;
00479     }
00480     return vector;
00481 }
00482 
00483 void Cylindrical2dMesh::SetNode(unsigned index, ChastePoint<2> point, bool concreteMove)
00484 {
00485     // Perform a periodic movement if necessary
00486     if (point.rGetLocation()[0] >= mWidth)
00487     {
00488         // Move point to the left
00489         point.SetCoordinate(0, point.rGetLocation()[0] - mWidth);
00490     }
00491     else if (point.rGetLocation()[0] < 0.0)
00492     {
00493         // Move point to the right
00494         point.SetCoordinate(0, point.rGetLocation()[0] + mWidth);
00495     }
00496 
00497     // Update the node's location
00498     MutableMesh<2,2>::SetNode(index, point, concreteMove);
00499 }
00500 
00501 double Cylindrical2dMesh::GetWidth(const unsigned& rDimension) const
00502 {
00503     double width = 0.0;
00504     assert(rDimension==0 || rDimension==1);
00505     if (rDimension==0)
00506     {
00507         width = mWidth;
00508     }
00509     else
00510     {
00511         width = MutableMesh<2,2>::GetWidth(rDimension);
00512     }
00513     return width;
00514 }
00515 
00516 unsigned Cylindrical2dMesh::AddNode(Node<2>* pNewNode)
00517 {
00518     unsigned node_index = MutableMesh<2,2>::AddNode(pNewNode);
00519 
00520     // If necessary move it to be back on the cylinder
00521     ChastePoint<2> new_node_point = pNewNode->GetPoint();
00522     SetNode(node_index, new_node_point, false);
00523 
00524     return node_index;
00525 }
00526 
00527 void Cylindrical2dMesh::CorrectNonPeriodicMesh()
00528 {
00529     GenerateVectorsOfElementsStraddlingPeriodicBoundaries();
00530 
00531     /*
00532      * Copy the member variables into new vectors, which we modify
00533      * by knocking out elements which pair up on each side.
00534      */
00535     std::set<unsigned> temp_left_hand_side_elements = mLeftPeriodicBoundaryElementIndices;
00536     std::set<unsigned> temp_right_hand_side_elements = mRightPeriodicBoundaryElementIndices;
00537     assert(mLeftPeriodicBoundaryElementIndices.size()==mRightPeriodicBoundaryElementIndices.size());
00538 
00539     // Go through all of the elements on the left periodic boundary
00540     for (std::set<unsigned>::iterator left_iter = mLeftPeriodicBoundaryElementIndices.begin();
00541          left_iter != mLeftPeriodicBoundaryElementIndices.end();
00542          ++left_iter)
00543     {
00544         unsigned elem_index = *left_iter;
00545         Element<2,2>* p_element = GetElement(elem_index);
00546 
00547         /*
00548          * Make lists of the nodes which the elements on the left contain and
00549          * the nodes which should be in a corresponding element on the right.
00550          */
00551         c_vector<unsigned,3> original_element_node_indices;
00552         c_vector<unsigned,3> corresponding_element_node_indices;
00553         for (unsigned i=0; i<3; i++)
00554         {
00555             original_element_node_indices[i] = p_element->GetNodeGlobalIndex(i);
00556             corresponding_element_node_indices[i] = GetCorrespondingNodeIndex(original_element_node_indices[i]);
00557         }
00558 
00559         // Search the right hand side elements for the corresponding element
00560         for (std::set<unsigned>::iterator right_iter = mRightPeriodicBoundaryElementIndices.begin();
00561              right_iter != mRightPeriodicBoundaryElementIndices.end();
00562              ++right_iter)
00563         {
00564             unsigned corresponding_elem_index = *right_iter;
00565             Element<2,2>* p_corresponding_element = GetElement(corresponding_elem_index);
00566 
00567             bool is_corresponding_node = true;
00568 
00569             for (unsigned i=0; i<3; i++)
00570             {
00571                 if ( (corresponding_element_node_indices[i] != p_corresponding_element->GetNodeGlobalIndex(0)) &&
00572                      (corresponding_element_node_indices[i] != p_corresponding_element->GetNodeGlobalIndex(1)) &&
00573                      (corresponding_element_node_indices[i] != p_corresponding_element->GetNodeGlobalIndex(2)) )
00574                 {
00575                     is_corresponding_node = false;
00576                     break;
00577                 }
00578             }
00579 
00580             if (is_corresponding_node)
00581             {
00582                 // Remove original and corresponding element from sets
00583                 temp_left_hand_side_elements.erase(elem_index);
00584                 temp_right_hand_side_elements.erase(corresponding_elem_index);
00585             }
00586         }
00587     }
00588 
00589     /*
00590      * If either of these ever throw you have more than one situation where the mesher has an option
00591      * of how to mesh. If it does ever throw you need to be cleverer and match up the
00592      * elements into as many pairs as possible on the left hand and right hand sides.
00593      */
00594     assert(temp_left_hand_side_elements.size() <= 2);
00595     assert(temp_right_hand_side_elements.size() <= 2);
00596 
00597     /*
00598      * Now we just have to use the first pair of elements and copy their info over to the other side.
00599      * First we need to get hold of both elements on either side.
00600      */
00601     if (temp_left_hand_side_elements.empty() || temp_right_hand_side_elements.empty())
00602     {
00603         assert(temp_right_hand_side_elements.empty());
00604         assert(temp_left_hand_side_elements.empty());
00605     }
00606     else
00607     {
00608         assert(temp_right_hand_side_elements.size() == 2 && temp_left_hand_side_elements.size() == 2);
00609         if (temp_right_hand_side_elements.size() == 2)
00610         {
00611             // Use the right hand side meshing and map to left
00612             UseTheseElementsToDecideMeshing(temp_right_hand_side_elements);
00613         }
00614         else
00615         {
00616             /*
00617              * If you get here there are more than two mixed up elements on the periodic edge.
00618              * We need to knock the pair out and then rerun this function. This shouldn't be
00619              * too hard to do but is as yet unnecessary.
00620              */
00621             NEVER_REACHED;
00622         }
00623     }
00624 }
00625 
00626 void Cylindrical2dMesh::UseTheseElementsToDecideMeshing(std::set<unsigned>& rMainSideElements)
00627 {
00628     assert(rMainSideElements.size() == 2);
00629 
00630     // We find the four nodes surrounding the dodgy meshing, on each side.
00631     std::set<unsigned> main_four_nodes;
00632     for (std::set<unsigned>::iterator left_iter = rMainSideElements.begin();
00633          left_iter != rMainSideElements.end();
00634          ++left_iter)
00635     {
00636         unsigned elem_index = *left_iter;
00637         Element<2,2>* p_element = GetElement(elem_index);
00638         for (unsigned i=0; i<3; i++)
00639         {
00640             unsigned index = p_element->GetNodeGlobalIndex(i);
00641             main_four_nodes.insert(index);
00642         }
00643     }
00644     assert(main_four_nodes.size() == 4);
00645 
00646     std::set<unsigned> other_four_nodes;
00647     for (std::set<unsigned>::iterator iter = main_four_nodes.begin();
00648          iter != main_four_nodes.end();
00649          ++iter)
00650     {
00651         other_four_nodes.insert(GetCorrespondingNodeIndex(*iter));
00652     }
00653     assert(other_four_nodes.size() == 4);
00654 
00655     /*
00656      * Find the elements surrounded by the nodes on the right
00657      * and change them to match the elements on the left.
00658      */
00659     std::vector<unsigned> corresponding_elements;
00660 
00661     // Loop over all elements
00662     for (MutableMesh<2,2>::ElementIterator elem_iter = GetElementIteratorBegin();
00663          elem_iter != GetElementIteratorEnd();
00664          ++elem_iter)
00665     {
00666         // Loop over the nodes of the element
00667         if (!(other_four_nodes.find(elem_iter->GetNodeGlobalIndex(0))==other_four_nodes.end()) &&
00668             !(other_four_nodes.find(elem_iter->GetNodeGlobalIndex(1))==other_four_nodes.end()) &&
00669             !(other_four_nodes.find(elem_iter->GetNodeGlobalIndex(2))==other_four_nodes.end()) )
00670         {
00671             corresponding_elements.push_back(elem_iter->GetIndex());
00672             elem_iter->MarkAsDeleted();
00673             mDeletedElementIndices.push_back(elem_iter->GetIndex());
00674         }
00675     }
00676     assert(corresponding_elements.size() == 2);
00677 
00678     // Now corresponding_elements contains the two elements which are going to be replaced by rMainSideElements
00679     unsigned num_elements = GetNumAllElements();
00680     for (std::set<unsigned>::iterator iter = rMainSideElements.begin();
00681          iter != rMainSideElements.end();
00682          ++iter)
00683     {
00684         Element<2,2>* p_main_element = GetElement(*iter);
00685         std::vector<Node<2>*> nodes;
00686 
00687         // Put corresponding nodes into a std::vector
00688         for (unsigned i=0; i<3; i++)
00689         {
00690             unsigned main_node = p_main_element->GetNodeGlobalIndex(i);
00691             nodes.push_back(this->GetNode(GetCorrespondingNodeIndex(main_node)));
00692         }
00693 
00694         // Make a new element
00695         Element<2,2>* p_new_element = new Element<2,2>(num_elements, nodes);
00696         this->mElements.push_back(p_new_element);
00697         this->mElementJacobians.push_back(zero_matrix<double>(2,2));
00698         this->mElementInverseJacobians.push_back(zero_matrix<double>(2,2));
00699         this->mElementJacobianDeterminants.push_back(0.0);
00700         num_elements++;
00701     }
00702 
00703     // Reindex to get rid of extra elements indices
00704     NodeMap map(GetNumAllNodes());
00705     this->ReIndex(map);
00706 }
00707 
00708 void Cylindrical2dMesh::GenerateVectorsOfElementsStraddlingPeriodicBoundaries()
00709 {
00710     mLeftPeriodicBoundaryElementIndices.clear();
00711     mRightPeriodicBoundaryElementIndices.clear();
00712 
00713     for (MutableMesh<2,2>::ElementIterator elem_iter = GetElementIteratorBegin();
00714          elem_iter != GetElementIteratorEnd();
00715          ++elem_iter)
00716     {
00717         // Left images are on the right of the mesh
00718         unsigned number_of_left_image_nodes = 0;
00719         unsigned number_of_right_image_nodes = 0;
00720         for (unsigned i=0; i<3; i++)
00721         {
00722             unsigned this_node_index = elem_iter->GetNodeGlobalIndex(i);
00723             if (mImageToLeftOriginalNodeMap.find(this_node_index) != mImageToLeftOriginalNodeMap.end())
00724             {
00725                 number_of_left_image_nodes++;
00726             }
00727             else if (mImageToRightOriginalNodeMap.find(this_node_index) != mImageToRightOriginalNodeMap.end())
00728             {
00729                 number_of_right_image_nodes++;
00730             }
00731         }
00732 
00733         // Elements on the left hand side (images of right nodes)
00734         if (number_of_right_image_nodes == 1 || number_of_right_image_nodes == 2)
00735         {
00736             mLeftPeriodicBoundaryElementIndices.insert(elem_iter->GetIndex());
00737         }
00738 
00739         // Elements on the right (images of left nodes)
00740         if (number_of_left_image_nodes == 1 || number_of_left_image_nodes == 2)
00741         {
00742             mRightPeriodicBoundaryElementIndices.insert(elem_iter->GetIndex());
00743         }
00744     }
00745 
00746     // Every boundary element on the left must have a corresponding element on the right
00747     assert(mLeftPeriodicBoundaryElementIndices.size() == mRightPeriodicBoundaryElementIndices.size());
00748 }
00749 
00750 unsigned Cylindrical2dMesh::GetCorrespondingNodeIndex(unsigned nodeIndex)
00751 {
00752     unsigned corresponding_node_index = UINT_MAX;
00753 
00754     // If nodeIndex is a member of mRightOriginals, then find the corresponding node index in mRightImages
00755     std::vector<unsigned>::iterator right_orig_iter = std::find(mRightOriginals.begin(), mRightOriginals.end(), nodeIndex);
00756     if (right_orig_iter != mRightOriginals.end())
00757     {
00758         corresponding_node_index = mRightImages[right_orig_iter - mRightOriginals.begin()];
00759     }
00760     else
00761     {
00762         // If nodeIndex is a member of mRightImages, then find the corresponding node index in mRightOriginals
00763         std::vector<unsigned>::iterator right_im_iter = std::find(mRightImages.begin(), mRightImages.end(), nodeIndex);
00764         if (right_im_iter != mRightImages.end())
00765         {
00766             corresponding_node_index = mRightOriginals[right_im_iter - mRightImages.begin()];
00767         }
00768         else
00769         {
00770             // If nodeIndex is a member of mLeftOriginals, then find the corresponding node index in mLeftImages
00771             std::vector<unsigned>::iterator left_orig_iter = std::find(mLeftOriginals.begin(), mLeftOriginals.end(), nodeIndex);
00772             if (left_orig_iter != mLeftOriginals.end())
00773             {
00774                 corresponding_node_index = mLeftImages[left_orig_iter - mLeftOriginals.begin()];
00775             }
00776             else
00777             {
00778                 // If nodeIndex is a member of mLeftImages, then find the corresponding node index in mLeftOriginals
00779                 std::vector<unsigned>::iterator left_im_iter = std::find(mLeftImages.begin(), mLeftImages.end(), nodeIndex);
00780                 if (left_im_iter != mLeftImages.end())
00781                 {
00782                     corresponding_node_index = mLeftOriginals[left_im_iter - mLeftImages.begin()];
00783                 }
00784             }
00785         }
00786     }
00787 
00788     // We must have found the corresponding node index
00789     assert(corresponding_node_index != UINT_MAX);
00790     return corresponding_node_index;
00791 }
00792 
00793 // Serialization for Boost >= 1.36
00794 #include "SerializationExportWrapperForCpp.hpp"
00795 CHASTE_CLASS_EXPORT(Cylindrical2dMesh)

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