Chaste  Release::3.4
HoneycombMeshGenerator.cpp
1 
2 /*
3 
4 Copyright (c) 2005-2016, University of Oxford.
5 All rights reserved.
6 
7 University of Oxford means the Chancellor, Masters and Scholars of the
8 University of Oxford, having an administrative office at Wellington
9 Square, Oxford OX1 2JD, UK.
10 
11 This file is part of Chaste.
12 
13 Redistribution and use in source and binary forms, with or without
14 modification, are permitted provided that the following conditions are met:
15  * Redistributions of source code must retain the above copyright notice,
16  this list of conditions and the following disclaimer.
17  * Redistributions in binary form must reproduce the above copyright notice,
18  this list of conditions and the following disclaimer in the documentation
19  and/or other materials provided with the distribution.
20  * Neither the name of the University of Oxford nor the names of its
21  contributors may be used to endorse or promote products derived from this
22  software without specific prior written permission.
23 
24 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
25 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
26 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
27 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
28 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
29 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
30 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
31 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
32 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
33 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34 
35 */
36 
37 #include "HoneycombMeshGenerator.hpp"
38 
39 #include <boost/foreach.hpp>
40 #include "TrianglesMeshReader.hpp"
41 #include "OutputFileHandler.hpp"
42 #include "RandomNumberGenerator.hpp"
43 #include "MathsCustomFunctions.hpp"
44 #include "ChasteSyscalls.hpp"
45 
46 HoneycombMeshGenerator::HoneycombMeshGenerator(unsigned numNodesAlongWidth, unsigned numNodesAlongLength, unsigned ghosts, double scaleFactor)
47  : mpMesh(NULL),
48  mMeshFilename("mesh"),
49  mDomainWidth(numNodesAlongWidth*scaleFactor),
50  mNumCellWidth(numNodesAlongWidth), //*1 because cells are considered to be size one
51  mNumCellLength(numNodesAlongLength)
52 {
53  // The code below won't work in parallel
54  assert(PetscTools::IsSequential());
55 
56  // An older version of the constructor might allow the wrong argument through to the scale factor
57  assert(scaleFactor > 0.0);
58 
59  // Get a unique temporary foldername
60  std::stringstream pid;
61  pid << getpid();
62  OutputFileHandler output_file_handler("2D_temporary_honeycomb_mesh_" + pid.str());
63  std::string output_dir = output_file_handler.GetOutputDirectoryFullPath();
64 
65  unsigned num_nodes_along_width = mNumCellWidth;
66  unsigned num_nodes_along_length = mNumCellLength;
67  double horizontal_spacing = mDomainWidth / (double)num_nodes_along_width;
68  double vertical_spacing = (sqrt(3.0)/2)*horizontal_spacing;
69 
70  // This line is needed to define ghost nodes later
71  mDomainDepth = (double)(num_nodes_along_length) * vertical_spacing;
72 
73  // Take account of ghost nodes
74  num_nodes_along_width += 2*ghosts;
75  num_nodes_along_length += 2*ghosts;
76 
77  unsigned num_nodes = num_nodes_along_width*num_nodes_along_length;
78  unsigned num_elem_along_width = num_nodes_along_width-1;
79  unsigned num_elem_along_length = num_nodes_along_length-1;
80  unsigned num_elem = 2*num_elem_along_width*num_elem_along_length;
81  unsigned num_edges = 3*num_elem_along_width*num_elem_along_length + num_elem_along_width + num_elem_along_length;
82 
83  double x0 = -horizontal_spacing*ghosts;
84  double y0 = -vertical_spacing*ghosts;
85 
86  mBottom = -vertical_spacing*ghosts;
87  mTop = mBottom + vertical_spacing*(num_nodes_along_length-1);
88 
89  // Write node file
90  out_stream p_node_file = output_file_handler.OpenOutputFile(mMeshFilename+".node");
91  (*p_node_file) << std::scientific;
92  (*p_node_file) << std::setprecision(20);
93  (*p_node_file) << num_nodes << "\t2\t0\t1" << std::endl;
94 
95  unsigned node = 0;
96  for (unsigned i=0; i<num_nodes_along_length; i++)
97  {
98  for (unsigned j=0; j<num_nodes_along_width; j++)
99  {
100  if (i<ghosts || i>=(ghosts+mNumCellLength))
101  {
102  mGhostNodeIndices.insert(node);
103  }
104  else if (j < ghosts || j >= (ghosts+mNumCellWidth))
105  {
106  mGhostNodeIndices.insert(node);
107  }
108  unsigned boundary = 0;
109  if (i==0 || i==num_nodes_along_length-1 || j==0 || j==num_nodes_along_width-1)
110  {
111  boundary = 1;
112  }
113 
114  double x = x0 + horizontal_spacing*((double)j + 0.25*(1.0+ SmallPow(-1.0,i+1)));
115  double y = y0 + vertical_spacing*(double)i;
116 
117  // Avoid floating point errors which upset OffLatticeSimulation
118  if ( (y<0.0) && (y>-1e-12) )
119  {
120  // Difficult to cover - just corrects floating point errors that have occurred from time to time!
121  #define COVERAGE_IGNORE
122  y = 0.0;
123  #undef COVERAGE_IGNORE
124  }
125 
126  (*p_node_file) << node++ << "\t" << x << "\t" << y << "\t" << boundary << std::endl;
127  }
128  }
129  p_node_file->close();
130 
131  // Write element file and edge file
132  out_stream p_elem_file = output_file_handler.OpenOutputFile(mMeshFilename+".ele");
133  (*p_elem_file) << std::scientific;
134 
135  out_stream p_edge_file = output_file_handler.OpenOutputFile(mMeshFilename+".edge");
136  (*p_node_file) << std::scientific;
137 
138  (*p_elem_file) << num_elem << "\t3\t0" << std::endl;
139  (*p_edge_file) << num_edges << "\t1" << std::endl;
140 
141  unsigned elem = 0;
142  unsigned edge = 0;
143  for (unsigned i=0; i<num_elem_along_length; i++)
144  {
145  for (unsigned j=0; j < num_elem_along_width; j++)
146  {
147  unsigned node0 = i*num_nodes_along_width + j;
148  unsigned node1 = i*num_nodes_along_width + j+1;
149  unsigned node2 = (i+1)*num_nodes_along_width + j;
150 
151  if (i%2 != 0)
152  {
153  node2 = node2 + 1;
154  }
155 
156  (*p_elem_file) << elem++ << "\t" << node0 << "\t" << node1 << "\t" << node2 << std::endl;
157 
158  unsigned horizontal_edge_is_boundary_edge = 0;
159  unsigned vertical_edge_is_boundary_edge = 0;
160  if (i==0)
161  {
162  horizontal_edge_is_boundary_edge = 1;
163  }
164  if (j==0 && i%2==0)
165  {
166  vertical_edge_is_boundary_edge = 1;
167  }
168 
169  (*p_edge_file) << edge++ << "\t" << node0 << "\t" << node1 << "\t" << horizontal_edge_is_boundary_edge << std::endl;
170  (*p_edge_file) << edge++ << "\t" << node1 << "\t" << node2 << "\t" << 0 << std::endl;
171  (*p_edge_file) << edge++ << "\t" << node2 << "\t" << node0 << "\t" << vertical_edge_is_boundary_edge << std::endl;
172 
173  node0 = i*num_nodes_along_width + j + 1;
174 
175  if (i%2 != 0)
176  {
177  node0 = node0 - 1;
178  }
179  node1 = (i+1)*num_nodes_along_width + j+1;
180  node2 = (i+1)*num_nodes_along_width + j;
181 
182  (*p_elem_file) << elem++ << "\t" << node0 << "\t" << node1 << "\t" << node2 << std::endl;
183  }
184  }
185 
186  for (unsigned i=0; i<num_elem_along_length; i++)
187  {
188  unsigned node0, node1;
189 
190  if (i%2==0)
191  {
192  node0 = (i+1)*num_nodes_along_width - 1;
193  node1 = (i+2)*num_nodes_along_width - 1;
194  }
195  else
196  {
197  node0 = (i+1)*num_nodes_along_width;
198  node1 = (i)*num_nodes_along_width;
199  }
200  (*p_edge_file) << edge++ << "\t" << node0 << "\t" << node1 << "\t" << 1 << std::endl;
201  }
202 
203  for (unsigned j=0; j<num_elem_along_width; j++)
204  {
205  unsigned node0 = num_nodes_along_width*(num_nodes_along_length-1) + j;
206  unsigned node1 = num_nodes_along_width*(num_nodes_along_length-1) + j+1;
207  (*p_edge_file) << edge++ << "\t" << node1 << "\t" << node0 << "\t" << 1 << std::endl;
208  }
209 
210  p_elem_file->close();
211  p_edge_file->close();
212 
213  // Having written the mesh to file, now construct it using TrianglesMeshReader.
214  // Nested scope so the reader closes files before we try deleting them below the scope.
215  {
216  TrianglesMeshReader<2,2> mesh_reader(output_dir + mMeshFilename);
217  mpMesh = new MutableMesh<2,2>;
218  mpMesh->ConstructFromMeshReader(mesh_reader);
219  }
220 
221  // Delete the temporary folder
222  output_file_handler.FindFile("").Remove();
223 
224  // The original files have been deleted, it is better if the mesh object forgets about them
226 }
227 
229 {
230  delete mpMesh;
231 }
232 
234 {
235  return mpMesh;
236 }
237 
239 {
240  std::vector<unsigned> location_indices;
241 
242  for (unsigned i=0; i<mpMesh->GetNumNodes(); i++)
243  {
244  if (mGhostNodeIndices.find(i)==mGhostNodeIndices.end())
245  {
246  location_indices.push_back(i);
247  }
248  }
249  return location_indices;
250 }
251 
253 {
254  if (!mGhostNodeIndices.empty())
255  {
256  EXCEPTION("Cannot call GetCircularMesh on a HoneycombMesh with ghost nodes");
257  }
258 
259  // Centre the mesh at (0,0)
260  c_vector<double,2> centre = zero_vector<double>(2);
261  for (unsigned i=0; i<mpMesh->GetNumNodes(); i++)
262  {
263  centre += mpMesh->GetNode(i)->rGetLocation();
264  }
265  centre /= (double)mpMesh->GetNumNodes();
266 
267  mpMesh->Translate(-centre[0], -centre[1]);
268 
269  // Iterate over nodes, deleting any that lie more than the specified radius from (0,0)
270  for (unsigned i=0; i<mpMesh->GetNumAllNodes(); i++)
271  {
272  if (norm_2(mpMesh->GetNode(i)->rGetLocation()) >= radius)
273  {
275  }
276  else
277  {
278  // Jiggle the data
279  c_vector<double,2>& r_location = mpMesh->GetNode(i)->rGetModifiableLocation();
280  c_vector<double,2> shift;
282  double max_jiggle = radius*5e-6;
283  shift[0] = max_jiggle*(p_gen->ranf()-0.5);
284  shift[1] = max_jiggle*(p_gen->ranf()-0.5);
285  r_location += shift;
286  }
287  }
288 
289  // Remesh
290  NodeMap map(mpMesh->GetNumNodes());
291  mpMesh->ReMesh(map);
292 
293  return mpMesh;
294 }
295 
297 {
298  return mDomainDepth;
299 }
300 
302 {
303  return mDomainWidth;
304 }
double SmallPow(double x, unsigned exponent)
void DeleteNodePriorToReMesh(unsigned index)
MutableMesh< 2, 2 > * GetCircularMesh(double radius)
virtual void Translate(const c_vector< double, SPACE_DIM > &rDisplacement)
Node< SPACE_DIM > * GetNode(unsigned index) const
#define EXCEPTION(message)
Definition: Exception.hpp:143
void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)
std::string GetOutputDirectoryFullPath() const
virtual unsigned GetNumAllNodes() const
unsigned GetNumNodes() const
static bool IsSequential()
Definition: PetscTools.cpp:91
std::vector< unsigned > GetCellLocationIndices()
std::set< unsigned > mGhostNodeIndices
virtual void ReMesh(NodeMap &map)
static RandomNumberGenerator * Instance()
virtual MutableMesh< 2, 2 > * GetMesh()
void SetMeshHasChangedSinceLoading()
MutableMesh< 2, 2 > * mpMesh