Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ToroidalHoneycombMeshGenerator.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include "ToroidalHoneycombMeshGenerator.hpp"
37
38#include <boost/foreach.hpp>
39#include <boost/make_shared.hpp>
40#include <boost/shared_ptr.hpp>
41#include "RandomNumberGenerator.hpp"
43#include "ChasteSyscalls.hpp"
44
45ToroidalHoneycombMeshGenerator::ToroidalHoneycombMeshGenerator(unsigned numNodesAlongWidth, unsigned numNodesAlongDepth, double widthScaleFactor, double depthScaleFactor)
46{
47 mpMesh = nullptr;
48 mDomainWidth = numNodesAlongWidth*widthScaleFactor;
49 mDomainDepth = 0.5*sqrt(3)*numNodesAlongDepth*depthScaleFactor;
50 mNumCellWidth = numNodesAlongWidth; //*1 because cells are considered to be size one
51 mNumCellLength = numNodesAlongDepth;
52 mMeshFilename = "mesh";
53
54 // The code below won't work in parallel
56
57 // An older version of the constructor might allow the wrong argument through to the scale factor
58 assert(widthScaleFactor > 0.0);
59 assert(depthScaleFactor > 0.0);
60
61 // Get a unique mesh filename
62 std::stringstream pid;
63 pid << getpid();
64
65 OutputFileHandler output_file_handler("2D_temporary_honeycomb_mesh_"+ pid.str());
66
67 unsigned num_nodes_along_width = mNumCellWidth;
68 unsigned num_nodes_along_depth = mNumCellLength;
69 double horizontal_spacing = mDomainWidth / (double)num_nodes_along_width;
70 double vertical_spacing = mDomainDepth / (double)num_nodes_along_depth;
71
72 unsigned num_nodes = num_nodes_along_width*num_nodes_along_depth;
73 unsigned num_elem_along_width = num_nodes_along_width-1;
74 unsigned num_elem_along_depth = num_nodes_along_depth-1;
75 unsigned num_elem = 2*num_elem_along_width*num_elem_along_depth;
76 unsigned num_edges = 3*num_elem_along_width*num_elem_along_depth + num_elem_along_width + num_elem_along_depth;
77
78 double x0 = 0.0;
79 double y0 = 0.0;
80
81 // Write node file
82 out_stream p_node_file = output_file_handler.OpenOutputFile(mMeshFilename+".node");
83 (*p_node_file) << std::scientific;
84 //(*p_node_file) << std::setprecision(20);
85 (*p_node_file) << num_nodes << "\t2\t0\t1" << std::endl;
86
87 unsigned node = 0;
88 for (unsigned i=0; i<num_nodes_along_depth; i++)
89 {
90 for (unsigned j=0; j<num_nodes_along_width; j++)
91 {
92 unsigned boundary = 0;
93 if ((i==0) || (i==num_nodes_along_depth-1))
94 {
95 boundary = 1;
96 }
97
98 const double adjustment = i % 2 != 0 ? 0.5 : 0.0;
99 double x = x0 + horizontal_spacing*((double)j + adjustment);
100 double y = y0 + vertical_spacing*(double)i;
101
102 // Avoid floating point errors which upset OffLatticeSimulation
103 if ((y<0.0) && (y>-1e-12))
104 {
105 // Difficult to cover - just corrects floating point errors that have occurred from time to time!
106 // LCOV_EXCL_START
107 y = 0.0;
108 // LCOV_EXCL_STOP
109 }
110
111 (*p_node_file) << node++ << "\t" << x << "\t" << y << "\t" << boundary << std::endl;
112 }
113 }
114 p_node_file->close();
115
116 // Write element file and edge file
117 out_stream p_elem_file = output_file_handler.OpenOutputFile(mMeshFilename+".ele");
118 (*p_elem_file) << std::scientific;
119
120 out_stream p_edge_file = output_file_handler.OpenOutputFile(mMeshFilename+".edge");
121 (*p_node_file) << std::scientific;
122
123 (*p_elem_file) << num_elem << "\t3\t0" << std::endl;
124 (*p_edge_file) << num_edges << "\t1" << std::endl;
125
126 unsigned elem = 0;
127 unsigned edge = 0;
128 for (unsigned i=0; i<num_elem_along_depth; i++)
129 {
130 for (unsigned j=0; j < num_elem_along_width; j++)
131 {
132 unsigned node0 = i*num_nodes_along_width + j;
133 unsigned node1 = i*num_nodes_along_width + j+1;
134 unsigned node2 = (i+1)*num_nodes_along_width + j;
135
136 if (i%2 != 0)
137 {
138 node2 = node2 + 1;
139 }
140
141 (*p_elem_file) << elem++ << "\t" << node0 << "\t" << node1 << "\t" << node2 << std::endl;
142
143 unsigned horizontal_edge_is_boundary_edge = 0;
144 unsigned vertical_edge_is_boundary_edge = 0;
145 if (i==0)
146 {
147 horizontal_edge_is_boundary_edge = 1;
148 }
149
150 (*p_edge_file) << edge++ << "\t" << node0 << "\t" << node1 << "\t" << horizontal_edge_is_boundary_edge << std::endl;
151 (*p_edge_file) << edge++ << "\t" << node1 << "\t" << node2 << "\t" << 0 << std::endl;
152 (*p_edge_file) << edge++ << "\t" << node2 << "\t" << node0 << "\t" << vertical_edge_is_boundary_edge << std::endl;
153
154 node0 = i*num_nodes_along_width + j + 1;
155
156 if (i%2 != 0)
157 {
158 node0 = node0 - 1;
159 }
160 node1 = (i+1)*num_nodes_along_width + j+1;
161 node2 = (i+1)*num_nodes_along_width + j;
162
163 (*p_elem_file) << elem++ << "\t" << node0 << "\t" << node1 << "\t" << node2 << std::endl;
164 }
165 }
166
167 for (unsigned i=0; i<num_elem_along_depth; i++)
168 {
169 unsigned node0, node1;
170
171 if (i%2==0)
172 {
173 node0 = (i+1)*num_nodes_along_width - 1;
174 node1 = (i+2)*num_nodes_along_width - 1;
175 }
176 else
177 {
178 node0 = (i+1)*num_nodes_along_width;
179 node1 = (i)*num_nodes_along_width;
180 }
181 (*p_edge_file) << edge++ << "\t" << node0 << "\t" << node1 << "\t" << 1 << std::endl;
182 }
183
184 for (unsigned j=0; j<num_elem_along_width; j++)
185 {
186 unsigned node0 = num_nodes_along_width*(num_nodes_along_depth-1) + j;
187 unsigned node1 = num_nodes_along_width*(num_nodes_along_depth-1) + j+1;
188 (*p_edge_file) << edge++ << "\t" << node1 << "\t" << node0 << "\t" << 1 << std::endl;
189 }
190
191 p_elem_file->close();
192 p_edge_file->close();
193
194 // Having written the mesh to file, now construct it using TrianglesMeshReader.
195 // Nested scope so the reader closes files before we delete them below.
196 {
197 TrianglesMeshReader<2,2> mesh_reader(output_file_handler.GetOutputDirectoryFullPath() + mMeshFilename);
198 mpMesh = boost::make_shared<Toroidal2dMesh >(mDomainWidth,mDomainDepth);
199 mpMesh->ConstructFromMeshReader(mesh_reader);
200 }
201
202 // Make the mesh Toroidal (we use Triangle library mode inside this ReMesh call)
203 mpMesh->ReMesh();
204
205 // Delete the temporary files
206 output_file_handler.FindFile("").Remove();
207
208 // The original files have been deleted, it is better if the mesh object forgets about them
209 mpMesh->SetMeshHasChangedSinceLoading();
210}
211
212boost::shared_ptr<MutableMesh<2,2> > ToroidalHoneycombMeshGenerator::GetMesh()
213{
214 EXCEPTION("A Toroidal mesh was created but a normal mesh is being requested.");
215 return mpMesh; // Not really
216}
217
218boost::shared_ptr<Toroidal2dMesh> ToroidalHoneycombMeshGenerator::GetToroidalMesh()
219{
220 return boost::static_pointer_cast<Toroidal2dMesh>(mpMesh);
221}
#define EXCEPTION(message)
void Remove() const
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()
boost::shared_ptr< Toroidal2dMesh > GetToroidalMesh()
ToroidalHoneycombMeshGenerator(unsigned numNodesAlongWidth, unsigned numNodesAlongDepth, double widthScaleFactor=1.0, double depthScaleFactor=1.0)
boost::shared_ptr< MutableMesh< 2, 2 > > GetMesh()