HoneycombMeshGenerator.cpp

00001 
00002 /*
00003 
00004 Copyright (C) University of Oxford, 2005-2011
00005 
00006 University of Oxford means the Chancellor, Masters and Scholars of the
00007 University of Oxford, having an administrative office at Wellington
00008 Square, Oxford OX1 2JD, UK.
00009 
00010 This file is part of Chaste.
00011 
00012 Chaste is free software: you can redistribute it and/or modify it
00013 under the terms of the GNU Lesser General Public License as published
00014 by the Free Software Foundation, either version 2.1 of the License, or
00015 (at your option) any later version.
00016 
00017 Chaste is distributed in the hope that it will be useful, but WITHOUT
00018 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00019 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00020 License for more details. The offer of Chaste under the terms of the
00021 License is subject to the License being interpreted in accordance with
00022 English Law and subject to any action against the University of Oxford
00023 being under the jurisdiction of the English Courts.
00024 
00025 You should have received a copy of the GNU Lesser General Public License
00026 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00027 
00028 */
00029 
00030 #include "HoneycombMeshGenerator.hpp"
00031 #include "RandomNumberGenerator.hpp"
00032 #include "UblasCustomFunctions.hpp"
00033 
00034 HoneycombMeshGenerator::HoneycombMeshGenerator(unsigned numNodesAlongWidth, unsigned numNodesAlongLength, unsigned ghosts, double scaleFactor)
00035   : mpMesh(NULL),
00036     mDomainWidth(numNodesAlongWidth*scaleFactor),
00037     mNumCellWidth(numNodesAlongWidth), //*1 because cells are considered to be size one
00038     mNumCellLength(numNodesAlongLength)
00039 {
00040     // The getpid code below won't work in parallel
00041     assert(PetscTools::IsSequential());
00042 
00043     // Get a unique mesh filename
00044     std::stringstream pid;
00045     pid << getpid();
00046     mMeshFilename = "2D_temporary_honeycomb_mesh_" + pid.str();
00047 
00048     mGhostNodeIndices.empty();
00049 
00050     OutputFileHandler output_file_handler("");
00051     std::string output_dir = output_file_handler.GetOutputDirectoryFullPath();
00052 
00053     if (PetscTools::AmMaster())
00054     {
00055         unsigned numNodesAlongWidth = mNumCellWidth;
00056         unsigned numNodesAlongLength = mNumCellLength;
00057         double horizontal_spacing = mDomainWidth / (double)numNodesAlongWidth;
00058         double vertical_spacing = (sqrt(3)/2)*horizontal_spacing;
00059 
00060         // This line is needed to define ghost nodes later
00061         mDomainDepth = (double)(numNodesAlongLength) * vertical_spacing;
00062 
00063         // Take account of ghost nodes
00064         numNodesAlongWidth = numNodesAlongWidth + 2*ghosts;
00065         numNodesAlongLength = numNodesAlongLength + 2*ghosts;
00066 
00067         unsigned num_nodes            = numNodesAlongWidth*numNodesAlongLength;
00068         unsigned num_elem_along_width = numNodesAlongWidth-1;
00069         unsigned num_elem_along_length = numNodesAlongLength-1;
00070         unsigned num_elem             = 2*num_elem_along_width*num_elem_along_length;
00071         unsigned num_edges            = 3*num_elem_along_width*num_elem_along_length + num_elem_along_width + num_elem_along_length;
00072 
00073         double x0 = -horizontal_spacing*ghosts;
00074         double y0 = -vertical_spacing*ghosts;
00075 
00076         mBottom = -vertical_spacing*ghosts;
00077         mTop = mBottom + vertical_spacing*(numNodesAlongLength-1);
00078 
00079         // Write node file
00080         out_stream p_node_file = output_file_handler.OpenOutputFile(mMeshFilename+".node");
00081         (*p_node_file) << std::scientific;
00082         (*p_node_file) << num_nodes << "\t2\t0\t1" << std::endl;
00083 
00084         unsigned node = 0;
00085         for (unsigned i=0; i<numNodesAlongLength; i++)
00086         {
00087             for (unsigned j=0; j<numNodesAlongWidth; j++)
00088             {
00089                 if (i<ghosts || i>=(ghosts+mNumCellLength))
00090                 {
00091                     mGhostNodeIndices.insert(node);
00092                 }
00093                 else if (j < ghosts || j >= (ghosts+mNumCellWidth))
00094                 {
00095                     mGhostNodeIndices.insert(node);
00096                 }
00097                 unsigned boundary = 0;
00098                 if (i==0 || i==numNodesAlongLength-1 || j==0 || j==numNodesAlongWidth-1)
00099                 {
00100                     boundary = 1;
00101                 }
00102 
00103                 double x = x0 + horizontal_spacing*((double)j + 0.25*(1.0+ SmallPow(-1,i+1)));
00104                 double y = y0 + vertical_spacing*(double)i;
00105 
00106                 // Avoid floating point errors which upset CellBasedSimulation
00107                 if ( (y<0.0) && (y>-1e-12) )
00108                 {
00109                     // Difficult to cover - just corrects floating point errors that have occurred from time to time!
00110                     #define COVERAGE_IGNORE
00111                     y = 0.0;
00112                     #undef COVERAGE_IGNORE
00113                 }
00114 
00115                 (*p_node_file) << node++ << "\t" << x << "\t" << y << "\t" << boundary << std::endl;
00116             }
00117         }
00118         p_node_file->close();
00119 
00120         // Write element file and edge file
00121         out_stream p_elem_file = output_file_handler.OpenOutputFile(mMeshFilename+".ele");
00122         (*p_elem_file) << std::scientific;
00123 
00124         out_stream p_edge_file = output_file_handler.OpenOutputFile(mMeshFilename+".edge");
00125         (*p_node_file) << std::scientific;
00126 
00127         (*p_elem_file) << num_elem << "\t3\t0" << std::endl;
00128         (*p_edge_file) << num_edges << "\t1" << std::endl;
00129 
00130         unsigned elem = 0;
00131         unsigned edge = 0;
00132         for (unsigned i=0; i<num_elem_along_length; i++)
00133         {
00134             for (unsigned j=0; j < num_elem_along_width; j++)
00135             {
00136                 unsigned node0 =     i*numNodesAlongWidth + j;
00137                 unsigned node1 =     i*numNodesAlongWidth + j+1;
00138                 unsigned node2 = (i+1)*numNodesAlongWidth + j;
00139 
00140                 if (i%2 != 0)
00141                 {
00142                     node2 = node2 + 1;
00143                 }
00144 
00145                 (*p_elem_file) << elem++ << "\t" << node0 << "\t" << node1 << "\t" << node2 << std::endl;
00146 
00147                 unsigned horizontal_edge_is_boundary_edge = 0;
00148                 unsigned vertical_edge_is_boundary_edge = 0;
00149                 if (i==0)
00150                 {
00151                     horizontal_edge_is_boundary_edge = 1;
00152                 }
00153                 if (j==0 && i%2==0)
00154                 {
00155                     vertical_edge_is_boundary_edge = 1;
00156                 }
00157 
00158                 (*p_edge_file) << edge++ << "\t" << node0 << "\t" << node1 << "\t" << horizontal_edge_is_boundary_edge << std::endl;
00159                 (*p_edge_file) << edge++ << "\t" << node1 << "\t" << node2 << "\t" << 0 << std::endl;
00160                 (*p_edge_file) << edge++ << "\t" << node2 << "\t" << node0 << "\t" << vertical_edge_is_boundary_edge << std::endl;
00161 
00162                 node0 = i*numNodesAlongWidth + j + 1;
00163 
00164                 if (i%2 != 0)
00165                 {
00166                     node0 = node0 - 1;
00167                 }
00168                 node1 = (i+1)*numNodesAlongWidth + j+1;
00169                 node2 = (i+1)*numNodesAlongWidth + j;
00170 
00171                 (*p_elem_file) << elem++ << "\t" << node0 << "\t" << node1 << "\t" << node2 << std::endl;
00172             }
00173         }
00174 
00175         for (unsigned i=0; i<num_elem_along_length; i++)
00176         {
00177             unsigned node0, node1;
00178 
00179             if (i%2==0)
00180             {
00181                  node0 = (i+1)*numNodesAlongWidth - 1;
00182                  node1 = (i+2)*numNodesAlongWidth - 1;
00183             }
00184             else
00185             {
00186                 node0 = (i+1)*numNodesAlongWidth;
00187                 node1 = (i)*numNodesAlongWidth;
00188             }
00189             (*p_edge_file) << edge++ << "\t" << node0 << "\t" << node1 << "\t" << 1 << std::endl;
00190         }
00191 
00192         for (unsigned j=0; j<num_elem_along_width; j++)
00193         {
00194             unsigned node0 = numNodesAlongWidth*(numNodesAlongLength-1) + j;
00195             unsigned node1 = numNodesAlongWidth*(numNodesAlongLength-1) + j+1;
00196             (*p_edge_file) << edge++ << "\t" << node1 << "\t" << node0 << "\t" << 1 << std::endl;
00197         }
00198 
00199         p_elem_file->close();
00200         p_edge_file->close();
00201     }
00202 
00203     // Wait for the new mesh to be available
00204     PetscTools::Barrier();
00205 
00206     // Having written the mesh to file, now construct it using TrianglesMeshReader
00207     TrianglesMeshReader<2,2> mesh_reader(output_dir + mMeshFilename);
00208     mpMesh = new MutableMesh<2,2>;
00209     mpMesh->ConstructFromMeshReader(mesh_reader);
00210 
00211     // Delete the temporary files
00212     std::string command = "rm " + output_dir + mMeshFilename + ".*";
00213     int return_value = system(command.c_str());
00214     if (return_value != 0)
00215     {
00216         // Can't figure out how to make this throw but seems as if it should be here?
00217         #define COVERAGE_IGNORE
00218         EXCEPTION("HoneycombMeshGenerator cannot delete temporary files\n");
00219         #undef COVERAGE_IGNORE
00220     }
00221 
00222     // The original files have been deleted, it is better if the mesh object forgets about them
00223     mpMesh->SetMeshHasChangedSinceLoading();
00224 }
00225 
00226 HoneycombMeshGenerator::~HoneycombMeshGenerator()
00227 {
00228     delete mpMesh;
00229 }
00230 
00231 MutableMesh<2,2>* HoneycombMeshGenerator::GetMesh()
00232 {
00233     return mpMesh;
00234 }
00235 
00236 std::vector<unsigned> HoneycombMeshGenerator::GetCellLocationIndices()
00237 {
00238     std::vector<unsigned> location_indices;
00239 
00240     for (unsigned i=0; i<mpMesh->GetNumNodes(); i++)
00241     {
00242         if (mGhostNodeIndices.find(i)==mGhostNodeIndices.end())
00243         {
00244             location_indices.push_back(i);
00245         }
00246     }
00247     return location_indices;
00248 }
00249 
00250 MutableMesh<2,2>* HoneycombMeshGenerator::GetCircularMesh(double radius)
00251 {
00252     // Centre the mesh at (0,0)
00253     c_vector<double,2> centre = zero_vector<double>(2);
00254     for (unsigned i=0; i<mpMesh->GetNumNodes(); i++)
00255     {
00256         centre += mpMesh->GetNode(i)->rGetLocation();
00257     }
00258     centre /= (double)mpMesh->GetNumNodes();
00259 
00260     mpMesh->Translate(-centre[0], -centre[1]);
00261 
00262     // Iterate over nodes, deleting any that lie more
00263     // than the specified radius from (0,0)
00264     for (unsigned i=0; i<mpMesh->GetNumAllNodes(); i++)
00265     {
00266         if ( norm_2(mpMesh->GetNode(i)->rGetLocation()) >= radius)
00267         {
00268             mpMesh->DeleteNodePriorToReMesh(i);
00269         }
00270         else
00271         {
00272             // Jiggle the data
00273             c_vector<double,2>& r_location = mpMesh->GetNode(i)->rGetModifiableLocation();
00274             c_vector<double,2> shift;
00275             RandomNumberGenerator* p_gen = RandomNumberGenerator::Instance();
00276             double max_jiggle = radius*5e-6;
00277             shift[0] = max_jiggle*(p_gen->ranf()-0.5);
00278             shift[1] = max_jiggle*(p_gen->ranf()-0.5);
00279             r_location += shift;
00280         }
00281     }
00282 
00283     // Remesh
00284     NodeMap map(mpMesh->GetNumNodes());
00285     mpMesh->ReMesh(map);
00286 
00287     return mpMesh;
00288 }
00289 
00290 double HoneycombMeshGenerator::GetDomainDepth()
00291 {
00292     return mDomainDepth;
00293 }
00294 
00295 double HoneycombMeshGenerator::GetDomainWidth()
00296 {
00297     return mDomainWidth;
00298 }

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