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