Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
CylindricalHoneycombMeshGenerator.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 "CylindricalHoneycombMeshGenerator.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
45CylindricalHoneycombMeshGenerator::CylindricalHoneycombMeshGenerator(unsigned numNodesAlongWidth, unsigned numNodesAlongLength, unsigned ghosts, double scaleFactor)
46{
47 mpMesh = nullptr;
48 mDomainWidth = numNodesAlongWidth*scaleFactor;
49 mNumCellWidth = numNodesAlongWidth; //*1 because cells are considered to be size one
50 mNumCellLength = numNodesAlongLength;
51 mMeshFilename = "mesh";
52 mGhostNodeIndices.clear();
53 // The code below won't work in parallel
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 mesh filename
60 std::stringstream pid;
61 pid << getpid();
62
63 OutputFileHandler output_file_handler("2D_temporary_honeycomb_mesh_"+ pid.str());
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_length += 2*ghosts;
75
76 unsigned num_nodes = num_nodes_along_width*num_nodes_along_length;
77 unsigned num_elem_along_width = num_nodes_along_width-1;
78 unsigned num_elem_along_length = num_nodes_along_length-1;
79 unsigned num_elem = 2*num_elem_along_width*num_elem_along_length;
80 unsigned num_edges = 3*num_elem_along_width*num_elem_along_length + num_elem_along_width + num_elem_along_length;
81
82 double x0 = 0;
83 double y0 = -vertical_spacing*ghosts;
84
85 mBottom = -vertical_spacing*ghosts;
86 mTop = mBottom + vertical_spacing*(num_nodes_along_length-1);
87
88 // Write node file
89 out_stream p_node_file = output_file_handler.OpenOutputFile(mMeshFilename+".node");
90 (*p_node_file) << std::scientific;
91 //(*p_node_file) << std::setprecision(20);
92 (*p_node_file) << num_nodes << "\t2\t0\t1" << std::endl;
93
94 unsigned node = 0;
95 for (unsigned i=0; i<num_nodes_along_length; i++)
96 {
97 for (unsigned j=0; j<num_nodes_along_width; j++)
98 {
99 if (i<ghosts || i>=(ghosts+mNumCellLength))
100 {
101 mGhostNodeIndices.insert(node);
102 }
103 unsigned boundary = 0;
104 if ((i==0) || (i==num_nodes_along_length-1))
105 {
106 boundary = 1;
107 }
108
109 const double adjustment = i % 2 != 0 ? 0.5 : 0.0;
110 double x = x0 + horizontal_spacing*((double)j + adjustment);
111 double y = y0 + vertical_spacing*(double)i;
112
113 // Avoid floating point errors which upset OffLatticeSimulation
114 if ((y<0.0) && (y>-1e-12))
115 {
116 // Difficult to cover - just corrects floating point errors that have occurred from time to time!
117 // LCOV_EXCL_START
118 y = 0.0;
119 // LCOV_EXCL_STOP
120 }
121
122 (*p_node_file) << node++ << "\t" << x << "\t" << y << "\t" << boundary << std::endl;
123 }
124 }
125 p_node_file->close();
126
127 // Write element file and edge file
128 out_stream p_elem_file = output_file_handler.OpenOutputFile(mMeshFilename+".ele");
129 (*p_elem_file) << std::scientific;
130
131 out_stream p_edge_file = output_file_handler.OpenOutputFile(mMeshFilename+".edge");
132 (*p_node_file) << std::scientific;
133
134 (*p_elem_file) << num_elem << "\t3\t0" << std::endl;
135 (*p_edge_file) << num_edges << "\t1" << std::endl;
136
137 unsigned elem = 0;
138 unsigned edge = 0;
139 for (unsigned i=0; i<num_elem_along_length; i++)
140 {
141 for (unsigned j=0; j < num_elem_along_width; j++)
142 {
143 unsigned node0 = i*num_nodes_along_width + j;
144 unsigned node1 = i*num_nodes_along_width + j+1;
145 unsigned node2 = (i+1)*num_nodes_along_width + j;
146
147 if (i%2 != 0)
148 {
149 node2 = node2 + 1;
150 }
151
152 (*p_elem_file) << elem++ << "\t" << node0 << "\t" << node1 << "\t" << node2 << std::endl;
153
154 unsigned horizontal_edge_is_boundary_edge = 0;
155 unsigned vertical_edge_is_boundary_edge = 0;
156 if (i==0)
157 {
158 horizontal_edge_is_boundary_edge = 1;
159 }
160
161 (*p_edge_file) << edge++ << "\t" << node0 << "\t" << node1 << "\t" << horizontal_edge_is_boundary_edge << std::endl;
162 (*p_edge_file) << edge++ << "\t" << node1 << "\t" << node2 << "\t" << 0 << std::endl;
163 (*p_edge_file) << edge++ << "\t" << node2 << "\t" << node0 << "\t" << vertical_edge_is_boundary_edge << std::endl;
164
165 node0 = i*num_nodes_along_width + j + 1;
166
167 if (i%2 != 0)
168 {
169 node0 = node0 - 1;
170 }
171 node1 = (i+1)*num_nodes_along_width + j+1;
172 node2 = (i+1)*num_nodes_along_width + j;
173
174 (*p_elem_file) << elem++ << "\t" << node0 << "\t" << node1 << "\t" << node2 << std::endl;
175 }
176 }
177
178 for (unsigned i=0; i<num_elem_along_length; i++)
179 {
180 unsigned node0, node1;
181
182 if (i%2==0)
183 {
184 node0 = (i+1)*num_nodes_along_width - 1;
185 node1 = (i+2)*num_nodes_along_width - 1;
186 }
187 else
188 {
189 node0 = (i+1)*num_nodes_along_width;
190 node1 = (i)*num_nodes_along_width;
191 }
192 (*p_edge_file) << edge++ << "\t" << node0 << "\t" << node1 << "\t" << 1 << std::endl;
193 }
194
195 for (unsigned j=0; j<num_elem_along_width; j++)
196 {
197 unsigned node0 = num_nodes_along_width*(num_nodes_along_length-1) + j;
198 unsigned node1 = num_nodes_along_width*(num_nodes_along_length-1) + j+1;
199 (*p_edge_file) << edge++ << "\t" << node1 << "\t" << node0 << "\t" << 1 << std::endl;
200 }
201
202 p_elem_file->close();
203 p_edge_file->close();
204
205 // Having written the mesh to file, now construct it using TrianglesMeshReader.
206 // Nested scope so the reader closes files before we delete them below.
207 {
208 TrianglesMeshReader<2,2> mesh_reader(output_file_handler.GetOutputDirectoryFullPath() + mMeshFilename);
209 mpMesh = boost::make_shared<Cylindrical2dMesh>(mDomainWidth);
210 mpMesh->ConstructFromMeshReader(mesh_reader);
211 }
212
213 // Make the mesh cylindrical (we use Triangle library mode inside this ReMesh call)
214 mpMesh->ReMesh();
215
216 // Delete the temporary files
217 output_file_handler.FindFile("").Remove();
218
219 // The original files have been deleted, it is better if the mesh object forgets about them
220 mpMesh->SetMeshHasChangedSinceLoading();
221}
222
223boost::shared_ptr<MutableMesh<2,2> > CylindricalHoneycombMeshGenerator::GetMesh()
224{
225 EXCEPTION("A cylindrical mesh was created but a normal mesh is being requested.");
226 return mpMesh; // Not really
227}
228
229boost::shared_ptr<Cylindrical2dMesh> CylindricalHoneycombMeshGenerator::GetCylindricalMesh()
230{
231 return boost::static_pointer_cast<Cylindrical2dMesh>(mpMesh);
232}
#define EXCEPTION(message)
boost::shared_ptr< Cylindrical2dMesh > GetCylindricalMesh()
boost::shared_ptr< MutableMesh< 2, 2 > > GetMesh()
CylindricalHoneycombMeshGenerator(unsigned numNodesAlongWidth, unsigned numNodesAlongLength, unsigned ghosts=3, double scaleFactor=1.0)
void Remove() const
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()