42 unsigned numNodesPerCell,
43 double ellipseExponent,
44 double cellAspectRatio,
49 unsigned numFluidMeshPoints,
54 assert(numCellsWide > 0);
55 assert(numNodesPerCell > 3);
56 assert(ellipseExponent > 0.0);
57 assert(cellAspectRatio > 0.0);
58 assert(fabs(randomYMult) < 2.0);
59 assert( (leakyLaminas && numFluidMeshPoints < UINT_MAX) || (!leakyLaminas) );
62 bool override_nodes_per_cell = numFluidMeshPoints < UINT_MAX;
65 unit_vector<double> x_unit(2,0);
66 unit_vector<double> y_unit(2,1);
69 double cell_width = 1.0 /
double(numCellsWide);
70 double cell_height = cellAspectRatio * cell_width;
72 if (cell_height > 0.8)
75 cell_width = cell_height / cellAspectRatio;
79 if (override_nodes_per_cell)
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);
87 const double width_scale_fac = absoluteGap ==
DOUBLE_UNSET ? 0.95 : (cell_width - absoluteGap) / cell_width;
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();
104 double apical_cutoff = p_gen->GetHeightOfTopSurface();
105 double basal_cutoff = cell_height - apical_cutoff;
106 double periapical_cutoff = apical_cutoff - basal_cutoff;
108 if (apical_cutoff < 0.5 * cell_height || apical_cutoff > cell_height || periapical_cutoff < basal_cutoff)
110 EXCEPTION(
"Something went wrong calculating the height of top surface of the cell.");
114 std::vector<unsigned> corner_indices(4, UINT_MAX);
119 for (
unsigned location = 0; location < locations.size(); location++)
121 if (locations[location][1] > apical_cutoff)
123 corner_indices[RIGHT_APICAL_CORNER] = location;
129 for (
unsigned location =
unsigned(0.25 * numNodesPerCell); location < locations.size(); location++)
131 if (locations[location][1] < apical_cutoff)
133 corner_indices[LEFT_APICAL_CORNER] = location - 1;
139 for (
unsigned location =
unsigned(0.5 * numNodesPerCell); location < locations.size(); location++)
141 if (locations[location][1] < basal_cutoff)
143 corner_indices[LEFT_BASAL_CORNER] = location;
149 for (
unsigned location =
unsigned(0.75 * numNodesPerCell); location < locations.size(); location++)
151 if (locations[location][1] > basal_cutoff)
153 corner_indices[RIGHT_BASAL_CORNER] = location - 1;
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) )
163 EXCEPTION(
"At least one corner not tagged");
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]) )
171 EXCEPTION(
"Something went wrong when tagging corner locations");
174 if ( corner_indices[LEFT_APICAL_CORNER] - corner_indices[RIGHT_APICAL_CORNER] !=
175 corner_indices[RIGHT_BASAL_CORNER] - corner_indices[LEFT_BASAL_CORNER] )
177 EXCEPTION(
"Apical and basal surfaces are different sizes");
181 std::vector<ImmersedBoundaryElement<2,2>*> ib_elements;
182 std::vector<ImmersedBoundaryElement<1,2>*> ib_laminas;
183 std::vector<Node<2>*> nodes;
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;
190 double lamina_node_spacing = 2.0 * (cell_height + cell_width) / numNodesPerCell;
192 if (leakyLaminas && override_nodes_per_cell)
195 lamina_node_spacing = 2.0 /
static_cast<double>(numFluidMeshPoints);
202 double lam_height = absoluteGap ==
DOUBLE_UNSET ? y_offset[1] - 0.05 * cell_height : y_offset[1] - absoluteGap;
204 unsigned num_lamina_nodes =
static_cast<unsigned>(floor(1.0 / lamina_node_spacing));
206 std::vector<Node<2>*> nodes_this_elem;
208 for (
unsigned node_idx = 0; node_idx < num_lamina_nodes; node_idx++)
211 c_vector<double, 2> location = lam_height * y_unit + ( 0.5 / num_lamina_nodes +
double(node_idx) / num_lamina_nodes ) * x_unit;
215 p_node->SetRegion(LAMINA_REGION);
217 nodes.push_back(p_node);
218 nodes_this_elem.push_back(p_node);
228 if (randomYMult != 0.0)
230 for (
auto& p_node : nodes)
234 for (
auto& p_lam : ib_laminas)
238 EXCEPTION(
"Currently no random y variation allowed with an apical lamina");
242 double lam_height = y_offset[1] + cell_height;
244 unsigned num_lamina_nodes =
static_cast<unsigned>(floor(1.0 / lamina_node_spacing));
246 std::vector<Node<2>*> nodes_this_elem;
248 for (
unsigned node_idx = 0; node_idx < num_lamina_nodes; node_idx++)
251 c_vector<double, 2> location = lam_height * y_unit + ( 0.5 / num_lamina_nodes +
double(node_idx) / num_lamina_nodes ) * x_unit;
255 p_node->SetRegion(LAMINA_REGION);
257 nodes.push_back(p_node);
258 nodes_this_elem.push_back(p_node);
267 for (
unsigned elem_idx = 0; elem_idx < numCellsWide; elem_idx++)
269 std::vector<Node<2>*> nodes_this_elem;
270 double random_variation = p_rand_gen->
ranf() * randomYMult;
271 c_vector<double, 2> scaled_location;
273 for (
unsigned location = 0; location < locations.size(); location++)
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];
285 if (locations[location][1] <= basal_cutoff)
287 if (locations[location][0] < 0.5 * cell_width)
296 else if (locations[location][1] <= periapical_cutoff)
298 if (locations[location][0] < 0.5 * cell_width)
307 else if (locations[location][1] < apical_cutoff)
309 if (locations[location][0] < 0.5 * cell_width)
311 p_node->
SetRegion(LEFT_PERIAPICAL_REGION);
315 p_node->
SetRegion(RIGHT_PERIAPICAL_REGION);
320 if (locations[location][0] < 0.5 * cell_width)
330 nodes.push_back(p_node);
331 nodes_this_elem.push_back(p_node);
337 std::vector<Node<2>*>& r_elem_corners = ib_elements.back()->rGetCornerNodes();
338 for (
unsigned corner = 0; corner < 4; corner++)
340 r_elem_corners.push_back(ib_elements.back()->GetNode(corner_indices[corner]));
344 if (override_nodes_per_cell)