Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ImmersedBoundaryPalisadeMeshGenerator.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 "ImmersedBoundaryPalisadeMeshGenerator.hpp"
38#include "ImmersedBoundaryEnumerations.hpp"
39#include "RandomNumberGenerator.hpp"
40
42 unsigned numNodesPerCell,
43 double ellipseExponent,
44 double cellAspectRatio,
45 double randomYMult,
46 bool basalLamina,
47 bool apicalLamina,
48 bool leakyLaminas,
49 unsigned numFluidMeshPoints,
50 double absoluteGap)
51 : mpMesh(NULL)
52{
53 // Check for sensible input
54 assert(numCellsWide > 0);
55 assert(numNodesPerCell > 3);
56 assert(ellipseExponent > 0.0);
57 assert(cellAspectRatio > 0.0); // aspect ratio is cell height / cell width
58 assert(fabs(randomYMult) < 2.0);
59 assert( (leakyLaminas && numFluidMeshPoints < UINT_MAX) || (!leakyLaminas) );
60
61 // If the number of fluid mesh points is specified, calculate an appropriate number of nodes per cell automatically
62 bool override_nodes_per_cell = numFluidMeshPoints < UINT_MAX;
63
64 // Helper vectors
65 unit_vector<double> x_unit(2,0);
66 unit_vector<double> y_unit(2,1);
67
68 // Calculate cell sizes and rescale if necessary
69 double cell_width = 1.0 / double(numCellsWide);
70 double cell_height = cellAspectRatio * cell_width;
71
72 if (cell_height > 0.8)
73 {
74 cell_height = 0.8;
75 cell_width = cell_height / cellAspectRatio;
76 }
77
78 // If we are overriding the number of nodes per cell, alter it so that the spacing ratio is approx 0.5
79 if (override_nodes_per_cell)
80 {
81 double cell_perimeter = 2.0 * (cell_height + cell_width);
82 double fluid_mesh_spacing = 1.0 / static_cast<double>(numFluidMeshPoints);
83 double target_node_spacing = 0.5 * fluid_mesh_spacing;
84 numNodesPerCell = static_cast<unsigned>(cell_perimeter / target_node_spacing);
85 }
86
87 const double width_scale_fac = absoluteGap == DOUBLE_UNSET ? 0.95 : (cell_width - absoluteGap) / cell_width;
88
89 // Generate a reference superellipse
90 std::unique_ptr<SuperellipseGenerator> p_gen(new SuperellipseGenerator(numNodesPerCell, ellipseExponent, cell_width, cell_height, 0.0, 0.0));
91 std::vector<c_vector<double, 2> > locations = p_gen->GetPointsAsVectors();
92
93 /*
94 * The top and bottom heights are the heights at which there is maximum curvature in the superellipse.
95 * These define the apical and basal domains.
96 * _____
97 * __/_____\__ Apical cutoff
98 * __|_____|__ Periapical cutoff
99 * | |
100 * | |
101 * __|_____|__ Basal cutoff
102 * \_____/__ y = 0
103 */
104 double apical_cutoff = p_gen->GetHeightOfTopSurface();
105 double basal_cutoff = cell_height - apical_cutoff;
106 double periapical_cutoff = apical_cutoff - basal_cutoff;
107
108 if (apical_cutoff < 0.5 * cell_height || apical_cutoff > cell_height || periapical_cutoff < basal_cutoff)
109 {
110 EXCEPTION("Something went wrong calculating the height of top surface of the cell."); //LCOV_EXCL_LINE
111 }
112
113 // Corner 0 is left-apical, 1 is right-apical, 2 is right-basal, 3 is left-basal
114 std::vector<unsigned> corner_indices(4, UINT_MAX);
115
116 // Node 0 is middle of right lateral side
117
118 // First anticlockwise is RIGHT_APICAL_CORNER
119 for (unsigned location = 0; location < locations.size(); location++)
120 {
121 if (locations[location][1] > apical_cutoff)
122 {
123 corner_indices[RIGHT_APICAL_CORNER] = location;
124 break;
125 }
126 }
127
128 // Second anticlockwise is LEFT_APICAL_CORNER
129 for (unsigned location = unsigned(0.25 * numNodesPerCell); location < locations.size(); location++)
130 {
131 if (locations[location][1] < apical_cutoff)
132 {
133 corner_indices[LEFT_APICAL_CORNER] = location - 1;
134 break;
135 }
136 }
137
138 // Third anticlockwise is LEFT_BASAL_CORNER
139 for (unsigned location = unsigned(0.5 * numNodesPerCell); location < locations.size(); location++)
140 {
141 if (locations[location][1] < basal_cutoff)
142 {
143 corner_indices[LEFT_BASAL_CORNER] = location;
144 break;
145 }
146 }
147
148 // Fourth anticlockwise is RIGHT_BASAL_CORNER
149 for (unsigned location = unsigned(0.75 * numNodesPerCell); location < locations.size(); location++)
150 {
151 if (locations[location][1] > basal_cutoff)
152 {
153 corner_indices[RIGHT_BASAL_CORNER] = location - 1;
154 break;
155 }
156 }
157
158 if ( (corner_indices[0] == UINT_MAX) ||
159 (corner_indices[1] == UINT_MAX) ||
160 (corner_indices[2] == UINT_MAX) ||
161 (corner_indices[3] == UINT_MAX) )
162 {
163 EXCEPTION("At least one corner not tagged"); //LCOV_EXCL_LINE
164 }
165
166 // Should have RIGHT_APICAL_CORNER < LEFT_APICAL_CORNER < LEFT_BASAL_CORNER < RIGHT_BASAL_CORNER
167 if ( (corner_indices[RIGHT_APICAL_CORNER] > corner_indices[LEFT_APICAL_CORNER]) ||
168 (corner_indices[LEFT_APICAL_CORNER] > corner_indices[LEFT_BASAL_CORNER]) ||
169 (corner_indices[LEFT_BASAL_CORNER] > corner_indices[RIGHT_BASAL_CORNER]) )
170 {
171 EXCEPTION("Something went wrong when tagging corner locations"); //LCOV_EXCL_LINE
172 }
173
174 if ( corner_indices[LEFT_APICAL_CORNER] - corner_indices[RIGHT_APICAL_CORNER] !=
175 corner_indices[RIGHT_BASAL_CORNER] - corner_indices[LEFT_BASAL_CORNER] )
176 {
177 EXCEPTION("Apical and basal surfaces are different sizes"); //LCOV_EXCL_LINE
178 }
179
180 // Create vectors of immersed boundary elements, laminas, nodes
181 std::vector<ImmersedBoundaryElement<2,2>*> ib_elements;
182 std::vector<ImmersedBoundaryElement<1,2>*> ib_laminas;
183 std::vector<Node<2>*> nodes;
184
185 // Helper c_vector for offsets in x and y
186 c_vector<double, 2> x_offset = x_unit * cell_width;
187 c_vector<double, 2> y_offset = y_unit * (1.0 - cell_height) / 2.0;
188
189 // Aim to give the laminas roughly the same node-spacing as the other cells...
190 double lamina_node_spacing = 2.0 * (cell_height + cell_width) / numNodesPerCell;
191 //... unless the laminas are to be leaky, in which case the spacing should be bigger
192 if (leakyLaminas && override_nodes_per_cell)
193 {
194 // If laminas are leaky, increase the spacing so fluid can flow through
195 lamina_node_spacing = 2.0 / static_cast<double>(numFluidMeshPoints);
196 }
197
198 // Add the basal lamina, if there is one
199 if (basalLamina)
200 {
201 // The height of the lamina is offset by a proportion of the cell height
202 double lam_height = absoluteGap == DOUBLE_UNSET ? y_offset[1] - 0.05 * cell_height : y_offset[1] - absoluteGap;
203
204 unsigned num_lamina_nodes = static_cast<unsigned>(floor(1.0 / lamina_node_spacing));
205
206 std::vector<Node<2>*> nodes_this_elem;
207
208 for (unsigned node_idx = 0; node_idx < num_lamina_nodes; node_idx++)
209 {
210 // Calculate location of node
211 c_vector<double, 2> location = lam_height * y_unit + ( 0.5 / num_lamina_nodes + double(node_idx) / num_lamina_nodes ) * x_unit;
212
213 // Create the new node
214 Node<2>* p_node = new Node<2>(nodes.size(), location, true);
215 p_node->SetRegion(LAMINA_REGION);
216
217 nodes.push_back(p_node);
218 nodes_this_elem.push_back(p_node);
219 }
220
221 // Create the lamina element
222 ib_laminas.push_back(new ImmersedBoundaryElement<1,2>(0, nodes_this_elem));
223 }
224
225 // Add the apical lamina, if there is one
226 if (apicalLamina)
227 {
228 if (randomYMult != 0.0)
229 {
230 for (auto& p_node : nodes)
231 {
232 delete p_node;
233 }
234 for (auto& p_lam : ib_laminas)
235 {
236 delete p_lam;
237 }
238 EXCEPTION("Currently no random y variation allowed with an apical lamina");
239 }
240
241 // The height of the lamina is offset by a proportion of the cell height
242 double lam_height = y_offset[1] + cell_height;
243
244 unsigned num_lamina_nodes = static_cast<unsigned>(floor(1.0 / lamina_node_spacing));
245
246 std::vector<Node<2>*> nodes_this_elem;
247
248 for (unsigned node_idx = 0; node_idx < num_lamina_nodes; node_idx++)
249 {
250 // Calculate location of node
251 c_vector<double, 2> location = lam_height * y_unit + ( 0.5 / num_lamina_nodes + double(node_idx) / num_lamina_nodes ) * x_unit;
252
253 // Create the new node
254 Node<2>* p_node = new Node<2>(nodes.size(), location, true);
255 p_node->SetRegion(LAMINA_REGION);
256
257 nodes.push_back(p_node);
258 nodes_this_elem.push_back(p_node);
259 }
260
261 // Create the lamina element
262 ib_laminas.push_back(new ImmersedBoundaryElement<1,2>(1, nodes_this_elem));
263 }
264
266
267 for (unsigned elem_idx = 0; elem_idx < numCellsWide; elem_idx++)
268 {
269 std::vector<Node<2>*> nodes_this_elem;
270 double random_variation = p_rand_gen->ranf() * randomYMult;
271 c_vector<double, 2> scaled_location;
272
273 for (unsigned location = 0; location < locations.size(); location++)
274 {
275 unsigned node_index = nodes.size();
276 scaled_location = width_scale_fac * locations[location] + x_offset * (elem_idx + 0.5);
277 scaled_location[0] = fmod(scaled_location[0], 1.0);
278 scaled_location[1] *= (1.0 + random_variation);
279 scaled_location[1] += y_offset[1];
280
281 // Create the new node
282 Node<2>* p_node = new Node<2>(node_index, scaled_location, true);
283
284 // Tag the node region
285 if (locations[location][1] <= basal_cutoff)
286 {
287 if (locations[location][0] < 0.5 * cell_width)
288 {
289 p_node->SetRegion(LEFT_BASAL_REGION);
290 }
291 else
292 {
293 p_node->SetRegion(RIGHT_BASAL_REGION);
294 }
295 }
296 else if (locations[location][1] <= periapical_cutoff)
297 {
298 if (locations[location][0] < 0.5 * cell_width)
299 {
300 p_node->SetRegion(LEFT_LATERAL_REGION);
301 }
302 else
303 {
304 p_node->SetRegion(RIGHT_LATERAL_REGION);
305 }
306 }
307 else if (locations[location][1] < apical_cutoff)
308 {
309 if (locations[location][0] < 0.5 * cell_width)
310 {
311 p_node->SetRegion(LEFT_PERIAPICAL_REGION);
312 }
313 else
314 {
315 p_node->SetRegion(RIGHT_PERIAPICAL_REGION);
316 }
317 }
318 else
319 {
320 if (locations[location][0] < 0.5 * cell_width)
321 {
322 p_node->SetRegion(LEFT_APICAL_REGION);
323 }
324 else
325 {
326 p_node->SetRegion(RIGHT_APICAL_REGION);
327 }
328 }
329
330 nodes.push_back(p_node);
331 nodes_this_elem.push_back(p_node);
332 }
333
334 ib_elements.push_back(new ImmersedBoundaryElement<2,2>(elem_idx, nodes_this_elem));
335
336 // Pass in the corners calculated above
337 std::vector<Node<2>*>& r_elem_corners = ib_elements.back()->rGetCornerNodes();
338 for (unsigned corner = 0; corner < 4; corner++)
339 {
340 r_elem_corners.push_back(ib_elements.back()->GetNode(corner_indices[corner]));
341 }
342 }
343
344 if (override_nodes_per_cell)
345 {
346 mpMesh = new ImmersedBoundaryMesh<2, 2>(nodes, ib_elements, ib_laminas, numFluidMeshPoints, numFluidMeshPoints);
347 }
348 else
349 {
350 mpMesh = new ImmersedBoundaryMesh<2, 2>(nodes, ib_elements, ib_laminas);
351 }
352}
353
358
const double DOUBLE_UNSET
Definition Exception.hpp:57
#define EXCEPTION(message)
Definition Node.hpp:59
void SetRegion(unsigned region)
Definition Node.cpp:430
static RandomNumberGenerator * Instance()