Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
ImmersedBoundaryMesh.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 "ImmersedBoundaryMesh.hpp"
37
38#include <algorithm>
39#include <cmath>
40
41#include "ImmersedBoundaryEnumerations.hpp"
43#include "RandomNumberGenerator.hpp"
45#include "Warnings.hpp"
46
47
48template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
52 unsigned numGridPtsX,
53 unsigned numGridPtsY)
54 : mNumGridPtsX(numGridPtsX),
55 mNumGridPtsY(numGridPtsY),
56 mCharacteristicNodeSpacing(DOUBLE_UNSET),
57 mElementDivisionSpacing(DOUBLE_UNSET),
58 mNeighbourDist(0.1),
59 mSummaryOfNodeLocations(0.0),
60 mCellRearrangementThreshold(0.05)
61{
62 // Clear mNodes and mElements
63 Clear();
64
65 switch (SPACE_DIM)
66 {
67 case 2:
68 m2dVelocityGrids.resize(extents[2][mNumGridPtsX][mNumGridPtsY]);
69 break;
70
71 case 3:
72 EXCEPTION("Not implemented yet in 3D");
73 break;
74
75 //LCOV_EXCL_START
76 default:
78 break;
79 //LCOV_EXCL_STOP
80 }
81
82 // Populate mNodes, mElements, and mLaminas
83 for (unsigned node_it = 0; node_it < nodes.size(); node_it++)
84 {
85 Node<SPACE_DIM>* p_temp_node = nodes[node_it];
86 this->mNodes.push_back(p_temp_node);
87 }
88 for (unsigned elem_it = 0; elem_it < elements.size(); elem_it++)
89 {
90 ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>* p_temp_element = elements[elem_it];
91 mElements.push_back(p_temp_element);
92 }
93 for (unsigned lam_it = 0; lam_it < laminas.size(); lam_it++)
94 {
95 ImmersedBoundaryElement<ELEMENT_DIM - 1, SPACE_DIM>* p_temp_lamina = laminas[lam_it];
96 mLaminas.push_back(p_temp_lamina);
97 }
98
99 // Register elements with nodes
100 for (unsigned elem_it = 0; elem_it < mElements.size(); elem_it++)
101 {
103
104 unsigned element_index = p_element->GetIndex();
105 unsigned num_nodes_in_element = p_element->GetNumNodes();
106
107 for (unsigned node_idx = 0; node_idx < num_nodes_in_element; node_idx++)
108 {
109 p_element->GetNode(node_idx)->AddElement(element_index);
110 }
111 }
112
113 // Register laminas with nodes
114 //\todo is there a way we can register laminas with nodes?
115
116 // Set characteristic node spacing to the average distance between nodes in elements
117 double total_perimeter = 0.0;
118 unsigned total_nodes = 0;
119 for (unsigned elem_it = 0; elem_it < mElements.size(); elem_it++)
120 {
121 total_perimeter += this->GetSurfaceAreaOfElement(elem_it);
122 total_nodes += mElements[elem_it]->GetNumNodes();
123 }
124 mCharacteristicNodeSpacing = total_perimeter / double(total_nodes);
125
126 // Position fluid sources at the centroid of each element, and set strength to zero
127 for (unsigned elem_it = 0; elem_it < elements.size(); elem_it++)
128 {
129 unsigned elem_idx = mElements[elem_it]->GetIndex();
130
131 // Create a new fluid source at the correct location for each element
132 unsigned source_idx = static_cast<unsigned>(mElementFluidSources.size());
133 c_vector<double, SPACE_DIM> source_location = this->GetCentroidOfElement(elem_idx);
134 mElementFluidSources.push_back(std::make_shared<FluidSource<SPACE_DIM>>(source_idx, source_location));
135
136 // Set source parameters
137 mElementFluidSources.back()->SetAssociatedElementIndex(elem_idx);
138 mElementFluidSources.back()->SetStrength(0.0);
139
140 // Associate source with element
141 mElements[elem_it]->SetFluidSource(mElementFluidSources.back());
142 }
143
144 //Set up a number of sources to balance any active sources associated with elements
145 double balancing_source_spacing = 2.0 * mCharacteristicNodeSpacing;
147 // We start at the characteristic spacing in from the left-hand end, and place a source every 2 spacings
148 double current_location = mCharacteristicNodeSpacing;
149
150 while (current_location < 1.0)
151 {
152 // Create a new fluid source at the current x-location and zero y-location
153 unsigned source_idx = static_cast<unsigned>(mBalancingFluidSources.size());
154 mBalancingFluidSources.push_back(std::make_shared<FluidSource<SPACE_DIM>>(source_idx, current_location));
156 mBalancingFluidSources.back()->SetStrength(0.0);
157
158 // Increment the current location
159 current_location += balancing_source_spacing;
160 }
161
162 // Calculate a default neighbour dist, as half the root of the average element volume
163 constexpr double power = 1.0 / SPACE_DIM;
164 const double total_volume_of_elems = std::accumulate(mElements.begin(), mElements.end(), 0.0,
166 {
167 return d + this->GetVolumeOfElement(a->GetIndex());
168 });
169 mNeighbourDist = 0.5 * std::pow(total_volume_of_elems / mElements.size(), power);
170
171 this->mMeshChangesDuringSimulation = true;
172}
173
174template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
176 [[maybe_unused]] unsigned index) // [[maybe_unused]] due to unused-but-set-parameter warning in GCC 7,8,9
177{
178 if constexpr (SPACE_DIM == 2)
179 {
180 c_vector<double, 3> moments = CalculateMomentsOfElement(index);
181
182 double discriminant = sqrt((moments(0) - moments(1)) * (moments(0) - moments(1)) + 4.0 * moments(2) * moments(2));
183
184 // Note that as the matrix of second moments of area is symmetric, both its eigenvalues are real
185 double largest_eigenvalue = (moments(0) + moments(1) + discriminant) * 0.5;
186 double smallest_eigenvalue = (moments(0) + moments(1) - discriminant) * 0.5;
188 double elongation_shape_factor = sqrt(largest_eigenvalue / smallest_eigenvalue);
189 return elongation_shape_factor;
190 }
191 else
192 {
194 }
195}
196
197bool CustomComparisonForVectorX(c_vector<double, 2> vecA, c_vector<double, 2> vecB)
198{
199 return vecA[0] < vecB[0];
200}
201
202template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
204{
205 if constexpr (SPACE_DIM == 2)
206 {
207 double total_length = 0.0;
208
209 // Get the current elements
210 std::vector<c_vector<double, 2> > centroids(mElements.size());
211 for (unsigned elem_it = 0; elem_it < mElements.size(); elem_it++)
212 {
213 centroids[elem_it] = this->GetCentroidOfElement(mElements[elem_it]->GetIndex());
214 }
215
216 // Sort centroids by X
217 std::sort(centroids.begin(), centroids.end(), CustomComparisonForVectorX);
218
219 // Calculate piecewise linear length connecting centroids
220 for (unsigned cent_it = 1; cent_it < centroids.size(); cent_it++)
221 {
222 total_length += norm_2(centroids[cent_it - 1] - centroids[cent_it]);
223 }
224
225 double straight_line_length = norm_2(centroids[0] - centroids[centroids.size() - 1]);
226
227 // Avoid division by zero
228 if (straight_line_length < 0.00001)
229 {
230 return 0;
231 }
232
233 return total_length / straight_line_length;
234 }
235 else
236 {
238 }
239}
240
241bool CustomComparisonForSkewnessMeasure(std::pair<unsigned, c_vector<double, 2> > pairA, std::pair<unsigned, c_vector<double, 2> > pairB)
242{
243 return pairA.second[0] < pairB.second[0];
244}
245
246template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
248 [[maybe_unused]] unsigned elemIndex, // [[maybe_unused]] due to unused-but-set-parameter warning in GCC 7,8,9
249 c_vector<double, SPACE_DIM> axis)
250{
251 /*
252 * Method outline:
253 *
254 * Given an arbitrary axis and a closed polygon, we calculate the skewness of the mass distribution of the polygon
255 * perpendicular to the axis. This is used as a measure of asymmetry.
256 *
257 * To simplify calculating the mass distribution, we translate the centroid of the element to the origin and rotate
258 * about the centroid so the axis is vertical; then we sort all the nodes in ascending order of their x-coordinate.
259 *
260 * For each node in order, we need to know the length of the intersection through the node of the vertical line with
261 * the polygon. Once calculated, we have a piecewise-linear PDF for the mass distribution, which can be normalised
262 * by the surface area of the polygon.
263 *
264 * By integrating the pdf directly, we can calculate exactly the necessary moments of the distribution needed for
265 * the skewness.
266 *
267 * This method does not work for concave shapes, and emits a warning but continues anyway. Results will be close
268 * to correct for mildly concave shapes or may be correct depending on the alignment of the element.
269 */
270
271 // This method only works in 2D
272 if constexpr (SPACE_DIM == 2 && ELEMENT_DIM == 2)
273 {
274 // Get relevant info about the element
275 ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>* p_elem = this->GetElement(elemIndex);
276 unsigned num_nodes = p_elem->GetNumNodes();
277 double area_of_elem = this->GetVolumeOfElement(elemIndex);
278 c_vector<double, SPACE_DIM> centroid = this->GetCentroidOfElement(elemIndex);
279
280 // Get the unit axis and trig terms for rotation
281 c_vector<double, SPACE_DIM> unit_axis = axis / norm_2(axis);
282 double sin_theta = unit_axis[0];
283 double cos_theta = unit_axis[1];
284
285 // We need the (rotated) node locations in two orders - original and ordered left-to-right.
286 // For the latter we need to keep track of index, so we store that as part of a pair.
287 std::vector<c_vector<double, SPACE_DIM> > node_locations_original_order;
288 std::vector<std::pair<unsigned, c_vector<double, SPACE_DIM> > > ordered_locations;
289
290 // Get the node locations of the current element relative to its centroid, and rotate them
291 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
292 {
293 const c_vector<double, SPACE_DIM>& node_location = p_elem->GetNode(node_idx)->rGetLocation();
294
295 c_vector<double, SPACE_DIM> displacement = this->GetVectorFromAtoB(centroid, node_location);
296
297 c_vector<double, SPACE_DIM> rotated_location;
298 rotated_location[0] = cos_theta * displacement[0] - sin_theta * displacement[1];
299 rotated_location[1] = sin_theta * displacement[0] + cos_theta * displacement[1];
301 node_locations_original_order.push_back(rotated_location);
302 }
303
304 // Fill up a vector of identical points, and sort it so nodes are ordered in ascending x value
305 for (unsigned i = 0; i < node_locations_original_order.size(); i++)
306 {
307 ordered_locations.push_back(std::pair<unsigned, c_vector<double, SPACE_DIM> >(i, node_locations_original_order[i]));
308 }
309
310 std::sort(ordered_locations.begin(), ordered_locations.end(), CustomComparisonForSkewnessMeasure);
311
312 /*
313 * For each node, we must find every place where the axis (now rotated to be vertical) intersects the polygon:
314 *
315 * |
316 * __|______
317 * / | \
318 * / | \
319 * /____|___ |
320 * | | |
321 * _____|__| |
322 * \ | |
323 * \ | /
324 * \__|_____/
325 * |
326 * |
327 * ^
328 * For instance, the number of times the vertical intersects the polygon above is 4 and, for each node, we need to
329 * find all such intersections. We can do this by checking where the dot product of the the vector a with the unit
330 * x direction changes sign as we iterate over the original node locations, where a is the vector from the current
331 * node to the test node.
332 */
333
334 // For each node, we keep track of all the y-locations where the vertical through the node meets the polygon
335 std::vector<std::vector<double> > knots(num_nodes);
336
337 // Iterate over ordered locations from left to right
338 for (unsigned location = 0; location < num_nodes; location++)
339 {
340 // Get the two parts of the pair
341 unsigned this_idx = ordered_locations[location].first;
342 // this_location is the location of the node relative to centroid
343 c_vector<double, SPACE_DIM> this_location = ordered_locations[location].second;
344
345 // The y-coordinate of the current location is always a knot
346 // because we are passing the vertical through this node
347 knots[location].push_back(this_location[1]);
348
349 // To calculate all the intersection points, we need to iterate over every other location and see, sequentially,
350 // if the x-coordinate of location i+1 and i+2 crosses the x-coordinate of the current location.
351 // i.e. check whether each other edge crosses the vertical through this_location/current node
352 unsigned next_idx = (this_idx + 1) % num_nodes;
353 c_vector<double, SPACE_DIM> to_previous = node_locations_original_order[next_idx] - this_location;
354 c_vector<double, SPACE_DIM> to_next;
355
356 for (unsigned node_idx = this_idx + 2; node_idx < this_idx + num_nodes; node_idx++)
358 unsigned idx = node_idx % num_nodes;
359
360 to_next = node_locations_original_order[idx] - this_location;
361
362 // If the segment between to_previous and to_next intersects the vertical through this_location, the clause
363 // in the if statement below will be triggered
364 if (to_previous[0] * to_next[0] <= 0.0)
366 // Find how far between to_previous and to_next the point of intersection is
367 double interp = 0.5;
368 if (to_previous[0] - to_next[0] != 0.0)
369 {
370 interp = to_previous[0] / (to_previous[0] - to_next[0]);
371 }
372
373 assert(interp >= 0.0 && interp <= 1.0);
375 // Record the y-value of the intersection point
376 double new_intersection = this_location[1] + to_previous[1] + interp * (to_next[1] - to_previous[1]);
377 knots[location].push_back(new_intersection);
378 }
380 to_previous = to_next;
381 }
382
383 if (knots[location].size() > 2)
384 {
385 WARN_ONCE_ONLY("Axis intersects polygon more than 2 times (concavity) - check element is fairly convex.");
386 }
387 }
388
389 // For ease, construct a vector of the x-locations of all the nodes, in order
390 std::vector<double> ordered_x(num_nodes);
391 for (unsigned location = 0; location < num_nodes; location++)
392 {
393 ordered_x[location] = ordered_locations[location].second[0];
395
396 // Calculate the mass contributions at each x-location - this is the length of the intersection of the vertical
397 // through each location
398 std::vector<double> mass_contributions(num_nodes);
399 for (unsigned i = 0; i < num_nodes; i++)
400 {
401 std::sort(knots[i].begin(), knots[i].end());
402
403 switch (knots[i].size())
405 case 1:
406 mass_contributions[i] = 0.0;
407 break;
408
409 case 2:
410 mass_contributions[i] = knots[i][1] - knots[i][0];
411 break;
412
413 default:
414 mass_contributions[i] += knots[i][knots[i].size()-1] - knots[i][0];
415 }
417 // Normalise, so that these lengths define a pdf
418 mass_contributions[i] /= area_of_elem;
419 }
420
421 // Calculate moments. Because we just have a bunch of linear segments, we can integrate the pdf exactly
422 double e_x0 = 0.0;
423 double e_x1 = 0.0;
424 double e_x2 = 0.0;
425 double e_x3 = 0.0;
427 for (unsigned i = 1; i < num_nodes; i++)
428 {
429 double x0 = ordered_x[i - 1];
430 double x1 = ordered_x[i];
431
432 double fx0 = mass_contributions[i - 1];
433 double fx1 = mass_contributions[i];
435 // We need squared, cubed, ..., order 5 for each x
436 double x0_2 = x0 * x0;
437 double x0_3 = x0_2 * x0;
438 double x0_4 = x0_3 * x0;
439 double x0_5 = x0_4 * x0;
440
441 double x1_2 = x1 * x1;
442 double x1_3 = x1_2 * x1;
443 double x1_4 = x1_3 * x1;
444 double x1_5 = x1_4 * x1;
445
446 if (x1 - x0 > 0)
447 {
448 // Calculate y = mx + c for this section of the pdf
449 double m = (fx1 - fx0) / (x1 - x0);
450 double c = fx0 - m * x0;
451
452 e_x0 += m * (x1_2 - x0_2) / 2.0 + c * (x1 - x0);
453 e_x1 += m * (x1_3 - x0_3) / 3.0 + c * (x1_2 - x0_2) / 2.0;
454 e_x2 += m * (x1_4 - x0_4) / 4.0 + c * (x1_3 - x0_3) / 3.0;
455 e_x3 += m * (x1_5 - x0_5) / 5.0 + c * (x1_4 - x0_4) / 4.0;
456 }
457 }
458
459 // Check that we have correctly defined a pdf
460 if (fabs(e_x0 - 1.0) < 1e-6)
461 {
462 WARN_ONCE_ONLY("Mass distribution of element calculated incorrectly due to element concavity. Skewness may not be correct!");
463 }
464
465 // Calculate the standard deviation, and return the skewness
466 double sd = sqrt(e_x2 - e_x1 * e_x1);
467 return (e_x3 - 3.0 * e_x1 * sd * sd - e_x1 * e_x1 * e_x1) / (sd * sd * sd);
468 }
469 else
470 {
472 }
473}
474
475template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
478 ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>* p_elem = this->GetElement(index);
479
480 // Get the location of node zero as a reference point
481 c_vector<double, SPACE_DIM> ref_point = p_elem->GetNode(0)->rGetLocation();
482
483 // Vector to represent the n-dimensional 'bottom left'-most node
484 c_vector<double, SPACE_DIM> bottom_left = zero_vector<double>(SPACE_DIM);
485
486 // Vector to represent the n-dimensional 'top right'-most node
487 c_vector<double, SPACE_DIM> top_right = zero_vector<double>(SPACE_DIM);
489 // Loop over all nodes in the element and update bottom_left and top_right, relative to node zero to account for periodicity
490 for (unsigned node_idx = 0; node_idx < p_elem->GetNumNodes(); node_idx++)
491 {
492 c_vector<double, SPACE_DIM> vec_to_node = this->GetVectorFromAtoB(ref_point, p_elem->GetNode(node_idx)->rGetLocation());
493
494 for (unsigned dim = 0; dim < SPACE_DIM; dim++)
495 {
496 if (vec_to_node[dim] < bottom_left[dim])
497 {
498 bottom_left[dim] = vec_to_node[dim];
500 else if (vec_to_node[dim] > top_right[dim])
501 {
502 top_right[dim] = vec_to_node[dim];
503 }
504 }
505 }
506
507 // Create Chaste points, rescaled by the location of node zero
508 ChastePoint<SPACE_DIM> min(bottom_left + ref_point);
509 ChastePoint<SPACE_DIM> max(top_right + ref_point);
510
511 return ChasteCuboid<SPACE_DIM>(min, max);
512}
513
514template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
516{
517 this->mMeshChangesDuringSimulation = false;
518 Clear();
519}
520
521template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
526
527template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
529{
530 // Create a set of neighbouring node indices
531 std::set<unsigned> neighbouring_node_indices;
532
533 // Find the indices of the elements owned by this node
534 std::set<unsigned> containing_elem_indices = this->GetNode(nodeIndex)->rGetContainingElementIndices();
535
536 // Iterate over these elements
537 for (std::set<unsigned>::iterator elem_iter = containing_elem_indices.begin();
538 elem_iter != containing_elem_indices.end();
539 ++elem_iter)
540 {
541 // Find the local index of this node in this element
542 unsigned local_index = GetElement(*elem_iter)->GetNodeLocalIndex(nodeIndex);
543
544 // Find the global indices of the preceding and successive nodes in this element
545 unsigned num_nodes = GetElement(*elem_iter)->GetNumNodes();
546 unsigned previous_local_index = (local_index + num_nodes - 1) % num_nodes;
547 unsigned next_local_index = (local_index + 1) % num_nodes;
548
549 // Add the global indices of these two nodes to the set of neighbouring node indices
550 neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(previous_local_index));
551 neighbouring_node_indices.insert(GetElement(*elem_iter)->GetNodeGlobalIndex(next_local_index));
552 }
553
554 return neighbouring_node_indices;
555}
556
557template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
559{
560 assert(index < this->mNodes.size());
561 return index;
562}
563
564template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
566{
567 assert(index < this->mElements.size());
568 return index;
569}
571template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
573{
574 return index;
575}
576
577template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
579{
580 // Delete elements
581 for (unsigned i = 0; i < mElements.size(); i++)
582 {
583 delete mElements[i];
584 }
585 mElements.clear();
586
587 // Delete laminas
588 for (auto lamina : mLaminas)
589 {
590 delete(lamina);
591 }
592 mLaminas.clear();
593
594 // Delete nodes
595 for (unsigned i = 0; i < this->mNodes.size(); i++)
596 {
597 delete this->mNodes[i];
599 this->mNodes.clear();
600
601 mBalancingFluidSources.clear();
602 mElementFluidSources.clear();
603}
604
605template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
607{
608 return mCharacteristicNodeSpacing;
609}
610
611template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
613{
614 return mCharacteristicNodeSpacing / (1.0 / double(mNumGridPtsX));
615}
616
617template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
618const std::vector<Node<SPACE_DIM>*>& ImmersedBoundaryMesh<ELEMENT_DIM, SPACE_DIM>::rGetNodes() const
619{
620 return this->mNodes;
621}
622
623template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
626 return mNumGridPtsX;
627}
628
629template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
631{
632 return mNumGridPtsY;
633}
634
635template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
637{
638 mNumGridPtsX = meshPointsX;
639 m2dVelocityGrids.resize(extents[2][mNumGridPtsX][mNumGridPtsY]);
641
642template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
644{
645 mNumGridPtsY = meshPointsY;
646 m2dVelocityGrids.resize(extents[2][mNumGridPtsX][mNumGridPtsY]);
647}
648
649template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
651{
652 mNumGridPtsX = numGridPts;
653 mNumGridPtsY = numGridPts;
654 m2dVelocityGrids.resize(extents[2][mNumGridPtsX][mNumGridPtsY]);
655}
656
657template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
659{
660 mCharacteristicNodeSpacing = nodeSpacing;
661}
662
663template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
664std::vector<std::shared_ptr<FluidSource<SPACE_DIM>>>& ImmersedBoundaryMesh<ELEMENT_DIM, SPACE_DIM>::rGetElementFluidSources()
665{
666 mElementFluidSources.clear();
667 for (const auto& element : mElements) {
668 if (element->GetFluidSource() != nullptr) {
669 mElementFluidSources.push_back(element->GetFluidSource());
670 }
671 }
672 return mElementFluidSources;
674
675template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
676std::vector<std::shared_ptr<FluidSource<SPACE_DIM>>>& ImmersedBoundaryMesh<ELEMENT_DIM, SPACE_DIM>::rGetBalancingFluidSources()
677{
678 return mBalancingFluidSources;
679}
680
681template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
684 return m2dVelocityGrids;
685}
686
687template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
690 return m2dVelocityGrids;
691}
692
693template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
695 Node<SPACE_DIM>* pNewNode)
697 pNewNode->SetIndex(this->mNodes.size());
698 this->mNodes.push_back(pNewNode);
699 return this->mNodes.size() - 1;
700}
701
702template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
705 mCellRearrangementThreshold = cellRearrangementThreshold;
706}
707
708template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
710{
711 return mCellRearrangementThreshold;
713
714template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
715c_vector<double, SPACE_DIM> ImmersedBoundaryMesh<ELEMENT_DIM, SPACE_DIM>::GetVectorFromAtoB(const c_vector<double, SPACE_DIM>& rLocation1, const c_vector<double, SPACE_DIM>& rLocation2)
716{
717 // This code currently assumes the grid is precisely [0,1)x[0,1)
718 c_vector<double, SPACE_DIM> vector = rLocation2 - rLocation1;
719
720 /*
721 * Handle the periodic condition here: if the points are more
722 * than 0.5 apart in any direction, choose -(1.0-dist).
723 */
724 for (unsigned dim = 0; dim < SPACE_DIM; dim++)
725 {
726 if (fabs(vector[dim]) > 0.5)
727 {
728 vector[dim] = copysign(fabs(vector[dim]) - 1.0, -vector[dim]);
730 }
731
732 return vector;
733}
734
735template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
738 this->mNodes[nodeIndex]->SetPoint(point);
739}
740
741template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
743{
744 return this->mNodes.size();
745}
746
747template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
749{
750 return mElements.size();
751}
752
753template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
755{
756 return mElements.size();
758
759template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
761{
762 return mLaminas.size();
763}
764
765template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
767{
768 assert(index < mElements.size());
769 return mElements[index];
770}
771
772template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
774{
775 assert(index < mLaminas.size());
776 return mLaminas[index];
777}
778
779template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
781{
782 return mNeighbourDist;
783}
784
785template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
787{
788 mNeighbourDist = neighbourDist;
789}
790
791template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
793 [[maybe_unused]] unsigned index) // [[maybe_unused]] due to unused-but-set-parameter warning in GCC 7,8,9
794{
795 // Only implemented in 2D
796 if constexpr (SPACE_DIM == 2)
797 {
798 ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
799
800 unsigned num_nodes = p_element->GetNumNodes();
801 c_vector<double, SPACE_DIM> centroid = zero_vector<double>(SPACE_DIM);
802
803 double centroid_x = 0;
804 double centroid_y = 0;
805
806 // Note that we cannot use GetVolumeOfElement() below as it returns the absolute, rather than signed, area
807 double element_signed_area = 0.0;
808
809 // Map the first vertex to the origin and employ GetVectorFromAtoB() to allow for periodicity
810 c_vector<double, SPACE_DIM> first_node_location = p_element->GetNodeLocation(0);
811 c_vector<double, SPACE_DIM> pos_1 = zero_vector<double>(SPACE_DIM);
812
813 // Loop over vertices
814 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
815 {
816 c_vector<double, SPACE_DIM> next_node_location = p_element->GetNodeLocation((local_index + 1) % num_nodes);
817 c_vector<double, SPACE_DIM> pos_2 = GetVectorFromAtoB(first_node_location, next_node_location);
818
819 double this_x = pos_1[0];
820 double this_y = pos_1[1];
821 double next_x = pos_2[0];
822 double next_y = pos_2[1];
823
824 double signed_area_term = this_x * next_y - this_y * next_x;
825
826 centroid_x += (this_x + next_x) * signed_area_term;
827 centroid_y += (this_y + next_y) * signed_area_term;
828 element_signed_area += 0.5 * signed_area_term;
829
830 pos_1 = pos_2;
831 }
832
833 assert(element_signed_area != 0.0);
834
835 // Finally, map back and employ GetVectorFromAtoB() to allow for periodicity
836 centroid = first_node_location;
837 centroid[0] += centroid_x / (6.0 * element_signed_area);
838 centroid[1] += centroid_y / (6.0 * element_signed_area);
839
840 centroid[0] = centroid[0] < 0 ? centroid[0] + 1.0 : fmod(centroid[0], 1.0);
841 centroid[1] = centroid[1] < 0 ? centroid[1] + 1.0 : fmod(centroid[1], 1.0);
842
843 return centroid;
844 }
845 else
846 {
848 }
849}
850
852template <>
855{
856 EXCEPTION("ImmersedBoundaryMesh not yet supported for the specified dimensions");
857}
858
860template <>
863{
864 EXCEPTION("ImmersedBoundaryMesh not yet supported for the specified dimensions");
865}
866
868template <>
871{
872 EXCEPTION("ImmersedBoundaryMesh not yet supported for the specified dimensions");
873}
874
876template <>
879{
880 EXCEPTION("ImmersedBoundaryMesh not yet supported for the specified dimensions");
881}
882
884template <>
887{
888 ImmersedBoundaryMeshReader<2, 2>& rIBMeshReader = dynamic_cast<ImmersedBoundaryMeshReader<2, 2>&>(rMeshReader);
889
890 assert(!rIBMeshReader.HasNodePermutation());
891
892 // Store numbers of nodes and elements
893 unsigned num_nodes = rIBMeshReader.GetNumNodes();
894 unsigned num_elements = rIBMeshReader.GetNumElements();
895 unsigned num_laminas = rIBMeshReader.GetNumLaminas();
896 this->mCharacteristicNodeSpacing = rIBMeshReader.GetCharacteristicNodeSpacing();
897
898 // Add nodes
899 rIBMeshReader.Reset();
900 mNodes.reserve(num_nodes);
901 std::vector<double> node_data;
902 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
903 {
904 node_data = rIBMeshReader.GetNextNode();
905 unsigned is_boundary_node = (bool)node_data[2];
906 node_data.pop_back();
907 this->mNodes.push_back(new Node<2>(node_idx, node_data, is_boundary_node));
908 }
909
910 // Add laminas
911 rIBMeshReader.Reset();
912 mLaminas.reserve(num_laminas);
913 for (unsigned lam_idx = 0; lam_idx < num_laminas; lam_idx++)
914 {
915 // Get the data for this element
917
918 // Get the nodes owned by this element
919 std::vector<Node<2>*> nodes;
920 unsigned num_nodes_in_lamina = lamina_data.NodeIndices.size();
921 for (unsigned node_idx = 0; node_idx < num_nodes_in_lamina; node_idx++)
922 {
923 assert(lamina_data.NodeIndices[node_idx] < this->mNodes.size());
924 nodes.push_back(this->mNodes[lamina_data.NodeIndices[node_idx]]);
925 }
926
927 // Use nodes and index to construct this element
928 ImmersedBoundaryElement<1, 2>* p_lamina = new ImmersedBoundaryElement<1, 2>(lam_idx, nodes);
929 mLaminas.push_back(p_lamina);
930
931 if (rIBMeshReader.GetNumLaminaAttributes() > 0)
932 {
933 assert(rIBMeshReader.GetNumLaminaAttributes() == 1);
934 unsigned attribute_value = lamina_data.AttributeValue;
935 p_lamina->SetAttribute(attribute_value);
936 }
937 }
938
939 // Add elements
940 rIBMeshReader.Reset();
941 mElements.reserve(num_elements);
942 for (unsigned elem_idx = 0; elem_idx < num_elements; elem_idx++)
943 {
944 // Get the data for this element
946
947 // Get the nodes owned by this element
948 std::vector<Node<2>*> nodes;
949 unsigned num_nodes_in_element = element_data.NodeIndices.size();
950 for (unsigned node_idx = 0; node_idx < num_nodes_in_element; node_idx++)
951 {
952 assert(element_data.NodeIndices[node_idx] < this->mNodes.size());
953 nodes.push_back(this->mNodes[element_data.NodeIndices[node_idx]]);
954 }
955
956 // Use nodes and index to construct this element
957 ImmersedBoundaryElement<2, 2>* p_element = new ImmersedBoundaryElement<2, 2>(elem_idx, nodes);
958 mElements.push_back(p_element);
959
960 if (rIBMeshReader.GetNumElementAttributes() > 0)
961 {
962 assert(rIBMeshReader.GetNumElementAttributes() == 1);
963 unsigned attribute_value = element_data.AttributeValue;
964 p_element->SetAttribute(attribute_value);
965 }
966 }
967
968 // Get grid dimensions from grid file and set up grids accordingly
969 this->mNumGridPtsX = rIBMeshReader.GetNumGridPtsX();
970 this->mNumGridPtsY = rIBMeshReader.GetNumGridPtsY();
971 m2dVelocityGrids.resize(extents[2][mNumGridPtsX][mNumGridPtsY]);
972
973 // Construct the velocity grids from mesh reader
974 for (unsigned dim = 0; dim < 2; dim++)
975 {
976 for (unsigned grid_row = 0; grid_row < mNumGridPtsY; grid_row++)
977 {
978 std::vector<double> next_row = rIBMeshReader.GetNextGridRow();
979 assert(next_row.size() == mNumGridPtsX);
980
981 for (unsigned i = 0; i < mNumGridPtsX; i++)
982 {
983 m2dVelocityGrids[dim][i][grid_row] = next_row[i];
984 }
985 }
986 }
987}
988
990template <>
993{
994 EXCEPTION("ImmersedBoundaryMesh not yet supported for the specified dimensions");
995}
996
997template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
999{
1000 if constexpr (SPACE_DIM == 2)
1001 {
1002 // Get pointer to this element
1003 ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
1004
1005 double element_volume = 0.0;
1006
1007 // We map the first vertex to the origin and employ GetVectorFromAtoB() to allow for periodicity
1008 c_vector<double, SPACE_DIM> first_node_location = p_element->GetNodeLocation(0);
1009 c_vector<double, SPACE_DIM> pos_1 = zero_vector<double>(SPACE_DIM);
1010
1011 unsigned num_nodes = p_element->GetNumNodes();
1012 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1013 {
1014 c_vector<double, SPACE_DIM> next_node_location = p_element->GetNodeLocation((local_index + 1) % num_nodes);
1015 c_vector<double, SPACE_DIM> pos_2 = GetVectorFromAtoB(first_node_location, next_node_location);
1016
1017 double this_x = pos_1[0];
1018 double this_y = pos_1[1];
1019 double next_x = pos_2[0];
1020 double next_y = pos_2[1];
1021
1022 element_volume += 0.5 * (this_x * next_y - next_x * this_y);
1023
1024 pos_1 = pos_2;
1025 }
1026
1027 // We take the absolute value just in case the nodes were really oriented clockwise
1028 return fabs(element_volume);
1029 }
1030 else
1031 {
1033 }
1034}
1035
1036template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1038{
1039 if constexpr (SPACE_DIM == 2)
1040 {
1041 // Get pointer to this element
1042 ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
1043
1044 double surface_area = 0.0;
1045 unsigned num_nodes = p_element->GetNumNodes();
1046 unsigned this_node_index = p_element->GetNodeGlobalIndex(0);
1047 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1048 {
1049 unsigned next_node_index = p_element->GetNodeGlobalIndex((local_index + 1) % num_nodes);
1050 surface_area += this->GetDistanceBetweenNodes(this_node_index, next_node_index);
1051 this_node_index = next_node_index;
1052 }
1053
1054 return surface_area;
1055 }
1056 else
1057 {
1059 }
1060}
1061
1062template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1064{
1065 if constexpr (SPACE_DIM == 2)
1066 {
1067 UpdateNodeLocationsVoronoiDiagramIfOutOfDate();
1068 ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(elemIdx);
1069 double surface_area = 0.0;
1070
1071 // Sum contributions from each node in the element
1072 for (unsigned local_idx = 0; local_idx < p_element->GetNumNodes(); local_idx++)
1073 {
1074 const unsigned global_node_idx = p_element->GetNodeGlobalIndex(local_idx);
1075
1076 // Get the voronoi cell that corresponds to this node
1077 const unsigned voronoi_cell_id = mVoronoiCellIdsIndexedByNodeIndex[global_node_idx];
1078 const auto& voronoi_cell = mNodeLocationsVoronoiDiagram.cells()[voronoi_cell_id];
1079
1080 // Iterate over the edges of this voronoi cell
1081 auto p_edge = voronoi_cell.incident_edge();
1082 do
1083 {
1084 // The global node index corresponding to a voronoi cell is encoded in its 'color' variable
1085 const unsigned twin_node_idx = p_edge->twin()->cell()->color();
1086 const unsigned twin_elem_idx = *this->GetNode(twin_node_idx)->ContainingElementsBegin();
1087
1088 // Check if the nodes are in different elements
1089 if (elemIdx != twin_elem_idx)
1090 {
1091 surface_area += this->CalculateLengthOfVoronoiEdge(*p_edge);
1092 }
1093
1094 p_edge = p_edge->next();
1095 } while (p_edge != voronoi_cell.incident_edge());
1096 }
1097
1098 return surface_area;
1099 }
1100 else
1101 {
1103 }
1104}
1105
1106// Excluded from coverage as it is impossible to construct a 1 dimensional element currently
1107//LCOV_EXCL_START
1113template <>
1115{
1116 return 0.0;
1117}
1118//LCOV_EXCL_STOP
1119
1120template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1122{
1123 if (recalculate || (this->GetElement(index)->GetAverageNodeSpacing() == DOUBLE_UNSET))
1124 {
1125 double average_node_spacing = this->GetSurfaceAreaOfElement(index) / this->GetElement(index)->GetNumNodes();
1126 this->GetElement(index)->SetAverageNodeSpacing(average_node_spacing);
1127
1128 return average_node_spacing;
1129 }
1130 else
1131 {
1132 return this->GetElement(index)->GetAverageNodeSpacing();
1133 }
1134}
1135
1136template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1138{
1139 if (recalculate || (this->GetLamina(index)->GetAverageNodeSpacing() == DOUBLE_UNSET))
1140 {
1141 ImmersedBoundaryElement<ELEMENT_DIM - 1, SPACE_DIM>* p_lam = this->GetLamina(index);
1142
1143 // Explicitly calculate the average node spacing
1144 double average_node_spacing = 0.0;
1145 for (unsigned node_it = 1; node_it < p_lam->GetNumNodes(); node_it++)
1146 {
1147 average_node_spacing += this->GetDistanceBetweenNodes(p_lam->GetNodeGlobalIndex(node_it),
1148 p_lam->GetNodeGlobalIndex(node_it - 1));
1149 }
1150
1151 average_node_spacing /= (p_lam->GetNumNodes() - 1);
1152
1153 // Set it for quick retrieval next time
1154 p_lam->SetAverageNodeSpacing(average_node_spacing);
1155 return average_node_spacing;
1156 }
1157 else
1158 {
1159 return this->GetLamina(index)->GetAverageNodeSpacing();
1160 }
1161}
1162
1163template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1165{
1166 return mElementDivisionSpacing;
1167}
1168
1169template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1171{
1172 mElementDivisionSpacing = elementDivisionSpacing;
1173}
1174
1176// 2D-specific methods //
1178
1179template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1181{
1182 if constexpr (SPACE_DIM == 2)
1183 {
1184 // Define helper variables
1185 ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>* p_element = GetElement(index);
1186 unsigned num_nodes = p_element->GetNumNodes();
1187 c_vector<double, 3> moments = zero_vector<double>(3);
1188
1189 // Since we compute I_xx, I_yy and I_xy about the centroid, we must shift each vertex accordingly
1190 c_vector<double, SPACE_DIM> centroid = GetCentroidOfElement(index);
1191
1192 c_vector<double, SPACE_DIM> this_node_location = p_element->GetNodeLocation(0);
1193 c_vector<double, SPACE_DIM> pos_1 = this->GetVectorFromAtoB(centroid, this_node_location);
1194
1195 for (unsigned local_index = 0; local_index < num_nodes; local_index++)
1196 {
1197 unsigned next_index = (local_index + 1) % num_nodes;
1198 c_vector<double, SPACE_DIM> next_node_location = p_element->GetNodeLocation(next_index);
1199 c_vector<double, SPACE_DIM> pos_2 = this->GetVectorFromAtoB(centroid, next_node_location);
1200
1201 double signed_area_term = pos_1(0) * pos_2(1) - pos_2(0) * pos_1(1);
1202 // Ixx
1203 moments(0) += (pos_1(1) * pos_1(1) + pos_1(1) * pos_2(1) + pos_2(1) * pos_2(1)) * signed_area_term;
1204
1205 // Iyy
1206 moments(1) += (pos_1(0) * pos_1(0) + pos_1(0) * pos_2(0) + pos_2(0) * pos_2(0)) * signed_area_term;
1207
1208 // Ixy
1209 moments(2) += (pos_1(0) * pos_2(1) + 2 * pos_1(0) * pos_1(1) + 2 * pos_2(0) * pos_2(1) + pos_2(0) * pos_1(1)) * signed_area_term;
1210
1211 pos_1 = pos_2;
1212 }
1213
1214 moments(0) /= 12;
1215 moments(1) /= 12;
1216 moments(2) /= 24;
1217
1218 /*
1219 * If the nodes owned by the element were supplied in a clockwise rather
1220 * than anticlockwise manner, or if this arose as a result of enforcing
1221 * periodicity, then our computed quantities will be the wrong sign, so
1222 * we need to fix this.
1223 */
1224 if (moments(0) < 0.0)
1225 {
1226 moments(0) = -moments(0);
1227 moments(1) = -moments(1);
1228 moments(2) = -moments(2);
1229 }
1230 return moments;
1231 }
1232 else
1233 {
1235 }
1236}
1237
1238template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1240{
1241 if constexpr (SPACE_DIM == 2)
1242 {
1243 c_vector<double, SPACE_DIM> short_axis = zero_vector<double>(SPACE_DIM);
1244
1245 // Calculate the moments of the element about its centroid (recall that I_xx and I_yy must be non-negative)
1246 c_vector<double, 3> moments = CalculateMomentsOfElement(index);
1247
1248 // Normalise the moments vector to remove problem of a very small discriminant (see #2874)
1249 moments /= norm_2(moments);
1250
1251 // If the principal moments are equal...
1252 double discriminant = (moments(0) - moments(1)) * (moments(0) - moments(1)) + 4.0 * moments(2) * moments(2);
1253 if (fabs(discriminant) < DBL_EPSILON)
1254 {
1255 // ...then every axis through the centroid is a principal axis, so return a random unit vector
1256 short_axis(0) = RandomNumberGenerator::Instance()->ranf();
1257 short_axis(1) = sqrt(1.0 - short_axis(0) * short_axis(0));
1258 }
1259 else
1260 {
1261 // If the product of inertia is zero, then the coordinate axes are the principal axes
1262 if (fabs(moments(2)) < DBL_EPSILON)
1263 {
1264 if (moments(0) < moments(1))
1265 {
1266 short_axis(0) = 0.0;
1267 short_axis(1) = 1.0;
1268 }
1269 else
1270 {
1271 short_axis(0) = 1.0;
1272 short_axis(1) = 0.0;
1273 }
1274 }
1275 else
1276 {
1277 // Otherwise we find the eigenvector of the inertia matrix corresponding to the largest eigenvalue
1278 double lambda = 0.5 * (moments(0) + moments(1) + sqrt(discriminant));
1279
1280 short_axis(0) = 1.0;
1281 short_axis(1) = (moments(0) - lambda) / moments(2);
1282
1283 // Normalise the short axis before returning it
1284 short_axis /= norm_2(short_axis);
1285 }
1286 }
1287
1288 return short_axis;
1289 }
1290 else
1291 {
1293 }
1294}
1295
1296template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1298 c_vector<double, SPACE_DIM> axisOfDivision,
1299 bool placeOriginalElementBelow)
1300{
1301 if constexpr (SPACE_DIM == 2 && ELEMENT_DIM == 2)
1302 {
1303 // Get the centroid of the element
1304 c_vector<double, SPACE_DIM> centroid = this->GetCentroidOfElement(pElement->GetIndex());
1305
1306 // Create a vector perpendicular to the axis of division
1307 c_vector<double, SPACE_DIM> perp_axis;
1308 perp_axis(0) = -axisOfDivision(1);
1309 perp_axis(1) = axisOfDivision(0);
1310
1311 /*
1312 * Find which edges the axis of division crosses by finding any node
1313 * that lies on the opposite side of the axis of division to its next
1314 * neighbour.
1315 */
1316
1317 unsigned num_nodes = pElement->GetNumNodes();
1318 std::vector<unsigned> intersecting_nodes;
1319 bool is_current_node_on_left = (inner_prod(this->GetVectorFromAtoB(pElement->GetNodeLocation(0), centroid), perp_axis) >= 0);
1320 for (unsigned i = 0; i < num_nodes; i++)
1321 {
1322 bool is_next_node_on_left = (inner_prod(this->GetVectorFromAtoB(pElement->GetNodeLocation((i + 1) % num_nodes), centroid), perp_axis) >= 0);
1323 if (is_current_node_on_left != is_next_node_on_left)
1324 {
1325 intersecting_nodes.push_back(i);
1326 }
1327 is_current_node_on_left = is_next_node_on_left;
1328 }
1329
1330 // If the axis of division does not cross two edges then we cannot proceed
1331 if (intersecting_nodes.size() != 2)
1332 {
1333 EXCEPTION("Cannot proceed with element division: the given axis of division does not cross two edges of the element"); // LCOV_EXCL_LINE
1334 }
1335
1336 // Now call DivideElement() to divide the element using the nodes found above
1337 //unsigned new_element_index = 0;
1338 unsigned new_element_index = DivideElement(pElement,
1339 intersecting_nodes[0],
1340 intersecting_nodes[1],
1341 centroid,
1342 axisOfDivision);
1343
1344 return new_element_index;
1345 }
1346 else
1347 {
1349 }
1350}
1351
1352template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1354 bool placeOriginalElementBelow)
1355{
1356 if constexpr (SPACE_DIM == 2 && ELEMENT_DIM == 2)
1357 {
1358 c_vector<double, SPACE_DIM> short_axis = this->GetShortAxisOfElement(pElement->GetIndex());
1359 unsigned new_element_index = DivideElementAlongGivenAxis(pElement, short_axis, placeOriginalElementBelow);
1360 return new_element_index;
1361 }
1362 else
1363 {
1365 }
1366}
1367
1368template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1370 unsigned nodeAIndex,
1371 unsigned nodeBIndex,
1372 c_vector<double, SPACE_DIM> centroid,
1373 c_vector<double, SPACE_DIM> axisOfDivision)
1374{
1375 if constexpr (SPACE_DIM == 2 && ELEMENT_DIM == 2)
1376 {
1377 if (mElementDivisionSpacing == DOUBLE_UNSET)
1378 {
1379 EXCEPTION("The value of mElementDivisionSpacing has not been set.");
1380 }
1381
1382 /*
1383 * Method outline:
1384 *
1385 * Each element needs to end up with the same number of nodes as the original element, and those nodes will be
1386 * equally spaced around the outline of each of the two daughter elements.
1387 *
1388 * The two elements need to be divided by a distance of mElementDivisionSpacing, where the distance is measured
1389 * perpendicular to the axis of division.
1390 *
1391 * To achieve this, we find four 'corner' locations, each of which has a perpendicular distance from the axis of
1392 * half the required spacing, and are found by using the locations from the existing element as a stencil.
1393 */
1394
1395 double half_spacing = 0.5 * mElementDivisionSpacing;
1396
1397 // Get unit vectors in the direction of the division axis, and the perpendicular
1398 c_vector<double, SPACE_DIM> unit_axis = axisOfDivision / norm_2(axisOfDivision);
1399 c_vector<double, SPACE_DIM> unit_perp;
1400 unit_perp[0] = -unit_axis[1];
1401 unit_perp[1] = unit_axis[0];
1402
1403 unsigned num_nodes = pElement->GetNumNodes();
1404
1405 /*
1406 * We first identify the start and end indices of the nodes which will form the location stencil for each daughter
1407 * cell. Our starting point is the node indices already identified.
1408 *
1409 * In order to ensure the resulting gap between the elements is the correct size, we remove as many nodes as
1410 * necessary until the perpendicular distance between the centroid and the node is at least half the required
1411 * spacing.
1412 *
1413 * Finally, we move the relevant node to be exactly half the required spacing.
1414 */
1415 unsigned start_a = (nodeAIndex + 1) % num_nodes;
1416 unsigned end_a = nodeBIndex;
1417
1418 unsigned start_b = (nodeBIndex + 1) % num_nodes;
1419 unsigned end_b = nodeAIndex;
1420
1421 // Find correct start_a
1422 bool no_node_satisfied_condition_1 = true;
1423 for (unsigned i = start_a; i != end_a;)
1424 {
1425 c_vector<double, SPACE_DIM> centroid_to_i = this->GetVectorFromAtoB(centroid, pElement->GetNode(i)->rGetLocation());
1426 double perpendicular_dist = inner_prod(centroid_to_i, unit_perp);
1427
1428 if (fabs(perpendicular_dist) >= half_spacing)
1429 {
1430 no_node_satisfied_condition_1 = false;
1431 start_a = i;
1432
1433 // Calculate position so it's exactly 0.5 * elem_spacing perpendicular distance from the centroid
1434 c_vector<double, SPACE_DIM> new_location = pElement->GetNode(i)->rGetLocation();
1435 new_location -= unit_perp * copysign(fabs(perpendicular_dist) - half_spacing, perpendicular_dist);
1436
1437 pElement->GetNode(i)->SetPoint(ChastePoint<SPACE_DIM>(new_location));
1438 break;
1439 }
1440
1441 // Go to the next node
1442 i = (i + 1) % num_nodes;
1443 }
1444
1445 // Find correct end_a
1446 bool no_node_satisfied_condition_2 = true;
1447 for (unsigned i = end_a; i != start_a;)
1448 {
1449 c_vector<double, SPACE_DIM> centroid_to_i = this->GetVectorFromAtoB(centroid, pElement->GetNode(i)->rGetLocation());
1450 double perpendicular_dist = inner_prod(centroid_to_i, unit_perp);
1451
1452 if (fabs(perpendicular_dist) >= half_spacing)
1453 {
1454 no_node_satisfied_condition_2 = false;
1455 end_a = i;
1456
1457 // Calculate position so it's exactly 0.5 * elem_spacing perpendicular distance from the centroid
1458 c_vector<double, SPACE_DIM> new_location = pElement->GetNode(i)->rGetLocation();
1459 new_location -= unit_perp * copysign(fabs(perpendicular_dist) - half_spacing, perpendicular_dist);
1460
1461 pElement->GetNode(i)->SetPoint(ChastePoint<SPACE_DIM>(new_location));
1462 break;
1463 }
1464
1465 // Go to the previous node
1466 i = (i + num_nodes - 1) % num_nodes; // LCOV_EXCL_LINE
1467 }
1468
1469 // Find correct start_b
1470 bool no_node_satisfied_condition_3 = true;
1471 for (unsigned i = start_b; i != end_b;)
1472 {
1473 c_vector<double, SPACE_DIM> centroid_to_i = this->GetVectorFromAtoB(centroid, pElement->GetNode(i)->rGetLocation());
1474 double perpendicular_dist = inner_prod(centroid_to_i, unit_perp);
1475
1476 if (fabs(perpendicular_dist) >= half_spacing)
1477 {
1478 no_node_satisfied_condition_3 = false;
1479 start_b = i;
1480
1481 // Calculate position so it's exactly 0.5 * elem_spacing perpendicular distance from the centroid
1482 c_vector<double, SPACE_DIM> new_location = pElement->GetNode(i)->rGetLocation();
1483 new_location -= unit_perp * copysign(fabs(perpendicular_dist) - half_spacing, perpendicular_dist);
1484
1485 pElement->GetNode(i)->SetPoint(ChastePoint<SPACE_DIM>(new_location));
1486 break;
1487 }
1488
1489 // Go to the next node
1490 i = (i + 1) % num_nodes;
1491 }
1492
1493 // Find correct end_b
1494 bool no_node_satisfied_condition_4 = true;
1495 for (unsigned i = end_b; i != start_b;)
1496 {
1497 c_vector<double, SPACE_DIM> centroid_to_i = this->GetVectorFromAtoB(centroid, pElement->GetNode(i)->rGetLocation());
1498 double perpendicular_dist = inner_prod(centroid_to_i, unit_perp);
1499
1500 if (fabs(perpendicular_dist) >= half_spacing)
1501 {
1502 no_node_satisfied_condition_4 = false;
1503 end_b = i;
1504
1505 // Calculate position so it's exactly 0.5 * elem_spacing perpendicular distance from the centroid
1506 c_vector<double, SPACE_DIM> new_location = pElement->GetNode(i)->rGetLocation();
1507 new_location -= unit_perp * copysign(fabs(perpendicular_dist) - half_spacing, perpendicular_dist);
1508
1509 pElement->GetNode(i)->SetPoint(ChastePoint<SPACE_DIM>(new_location));
1510 break;
1511 }
1512
1513 // Go to the previous node
1514 i = (i + num_nodes - 1) % num_nodes;
1515 }
1516
1517 if (no_node_satisfied_condition_1 || no_node_satisfied_condition_2 || no_node_satisfied_condition_3 || no_node_satisfied_condition_4)
1518 {
1519 EXCEPTION("Could not space elements far enough apart during cell division. Cannot currently handle this case");
1520 }
1521
1522 /*
1523 * Create location stencils for each of the daughter cells
1524 */
1525 std::vector<c_vector<double, SPACE_DIM> > daughter_a_location_stencil;
1526 for (unsigned node_idx = start_a; node_idx != (end_a + 1) % num_nodes;)
1527 {
1528 daughter_a_location_stencil.push_back(c_vector<double, SPACE_DIM>(pElement->GetNode(node_idx)->rGetLocation()));
1529
1530 // Go to next node
1531 node_idx = (node_idx + 1) % num_nodes;
1532 }
1533
1534 std::vector<c_vector<double, SPACE_DIM> > daughter_b_location_stencil;
1535 for (unsigned node_idx = start_b; node_idx != (end_b + 1) % num_nodes;)
1536 {
1537 daughter_b_location_stencil.push_back(c_vector<double, SPACE_DIM>(pElement->GetNode(node_idx)->rGetLocation()));
1538
1539 // Go to next node
1540 node_idx = (node_idx + 1) % num_nodes;
1541 }
1542
1543 assert(daughter_a_location_stencil.size() > 1);
1544 assert(daughter_b_location_stencil.size() > 1);
1545
1546 // To help calculating cumulative distances, add the first location on to the end
1547 daughter_a_location_stencil.push_back(daughter_a_location_stencil[0]);
1548 daughter_b_location_stencil.push_back(daughter_b_location_stencil[0]);
1549
1550 // Calculate the cumulative distances around the stencils
1551 std::vector<double> cumulative_distances_a;
1552 std::vector<double> cumulative_distances_b;
1553 cumulative_distances_a.push_back(0.0);
1554 cumulative_distances_b.push_back(0.0);
1555 for (unsigned loc_idx = 1; loc_idx < daughter_a_location_stencil.size(); loc_idx++)
1556 {
1557 cumulative_distances_a.push_back(cumulative_distances_a.back() + norm_2(this->GetVectorFromAtoB(daughter_a_location_stencil[loc_idx - 1], daughter_a_location_stencil[loc_idx])));
1558 }
1559 for (unsigned loc_idx = 1; loc_idx < daughter_b_location_stencil.size(); loc_idx++)
1560 {
1561 cumulative_distances_b.push_back(cumulative_distances_b.back() + norm_2(this->GetVectorFromAtoB(daughter_b_location_stencil[loc_idx - 1], daughter_b_location_stencil[loc_idx])));
1562 }
1563
1564 // Find the target node spacing for each of the daughter elements
1565 double target_spacing_a = cumulative_distances_a.back() / (double)num_nodes;
1566 double target_spacing_b = cumulative_distances_b.back() / (double)num_nodes;
1567
1568 // Move the existing nodes into position to become daughter-A nodes
1569 unsigned last_idx_used = 0;
1570 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
1571 {
1572 double location_along_arc = (double)node_idx * target_spacing_a;
1573
1574 while (location_along_arc > cumulative_distances_a[last_idx_used + 1])
1575 {
1576 last_idx_used++;
1577 }
1578
1579 // Interpolant is the extra distance past the last index used divided by the length of the next line segment
1580 double interpolant = (location_along_arc - cumulative_distances_a[last_idx_used]) / (cumulative_distances_a[last_idx_used + 1] - cumulative_distances_a[last_idx_used]);
1581
1582 c_vector<double, SPACE_DIM> this_to_next = this->GetVectorFromAtoB(daughter_a_location_stencil[last_idx_used],
1583 daughter_a_location_stencil[last_idx_used + 1]);
1584
1585 c_vector<double, SPACE_DIM> new_location_a = daughter_a_location_stencil[last_idx_used] + interpolant * this_to_next;
1586
1587 pElement->GetNode(node_idx)->SetPoint(ChastePoint<SPACE_DIM>(new_location_a));
1588 }
1589
1590 // Create new nodes at positions around the daughter-B stencil
1591 last_idx_used = 0;
1592 std::vector<Node<SPACE_DIM>*> new_nodes_vec;
1593 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
1594 {
1595 double location_along_arc = (double)node_idx * target_spacing_b;
1596
1597 while (location_along_arc > cumulative_distances_b[last_idx_used + 1])
1598 {
1599 last_idx_used++;
1600 }
1601
1602 // Interpolant is the extra distance past the last index used divided by the length of the next line segment
1603 double interpolant = (location_along_arc - cumulative_distances_b[last_idx_used]) / (cumulative_distances_b[last_idx_used + 1] - cumulative_distances_b[last_idx_used]);
1604
1605 c_vector<double, SPACE_DIM> this_to_next = this->GetVectorFromAtoB(daughter_b_location_stencil[last_idx_used],
1606 daughter_b_location_stencil[last_idx_used + 1]);
1607
1608 c_vector<double, SPACE_DIM> new_location_b = daughter_b_location_stencil[last_idx_used] + interpolant * this_to_next;
1609
1610 unsigned new_node_idx = this->mNodes.size();
1611 this->mNodes.push_back(new Node<SPACE_DIM>(new_node_idx, new_location_b, true));
1612 new_nodes_vec.push_back(this->mNodes.back());
1613 }
1614
1615 // Copy node attributes
1616 for (unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
1617 {
1618 new_nodes_vec[node_idx]->SetRegion(pElement->GetNode(node_idx)->GetRegion());
1619
1620 for (unsigned node_attribute = 0; node_attribute < pElement->GetNode(node_idx)->GetNumNodeAttributes(); node_attribute++)
1621 {
1622 new_nodes_vec[node_idx]->AddNodeAttribute(pElement->GetNode(node_idx)->rGetNodeAttributes()[node_attribute]);
1623 }
1624 }
1625
1626 // Create the new element
1627 unsigned new_elem_idx = this->mElements.size();
1628 this->mElements.push_back(new ImmersedBoundaryElement<ELEMENT_DIM, SPACE_DIM>(new_elem_idx, new_nodes_vec));
1629 this->mElements.back()->RegisterWithNodes();
1630
1631 // Copy any element attributes
1632 for (unsigned elem_attribute = 0; elem_attribute < pElement->GetNumElementAttributes(); elem_attribute++)
1633 {
1634 this->mElements.back()->AddElementAttribute(pElement->rGetElementAttributes()[elem_attribute]);
1635 }
1636
1637 // Add the necessary corners to keep consistency with the other daughter element
1638 for (unsigned corner = 0; corner < pElement->rGetCornerNodes().size(); corner++)
1639 {
1640 this->mElements.back()->rGetCornerNodes().push_back(pElement->rGetCornerNodes()[corner]);
1641 }
1642
1643 // Update fluid source location for the existing element
1644 pElement->GetFluidSource()->rGetModifiableLocation() = this->GetCentroidOfElement(pElement->GetIndex());
1645
1646 // Add a fluid source for the new element
1647 c_vector<double, SPACE_DIM> new_centroid = this->GetCentroidOfElement(new_elem_idx);
1648 mElementFluidSources.push_back(std::make_shared<FluidSource<SPACE_DIM>>(new_elem_idx, new_centroid));
1649
1650 // Set source parameters
1651 mElementFluidSources.back()->SetAssociatedElementIndex(new_elem_idx);
1652 mElementFluidSources.back()->SetStrength(0.0);
1653
1654 // Associate source with element
1655 mElements[new_elem_idx]->SetFluidSource(mElementFluidSources.back());
1656
1657 return new_elem_idx;
1658 }
1659 else
1660 {
1662 }
1663}
1664
1665template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1667{
1668 // Iterate over elements and remesh each one
1669 for (auto p_elem : mElements)
1670 {
1671 ReMeshElement(p_elem, randomOrder);
1672 }
1673
1674 // Iterate over laminas and remesh each one
1675 for (auto p_lam : mLaminas)
1676 {
1677 ReMeshLamina(p_lam, randomOrder);
1678 }
1679
1680 // Reposition fluid sources to the centroid of cells
1681 for (auto& p_source : mElementFluidSources)
1682 {
1683 p_source->rGetModifiableLocation() = this->GetCentroidOfElement(p_source->GetAssociatedElementIndex());
1684 }
1685}
1686
1687template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1689 bool randomOrder)
1690{
1691 if constexpr (SPACE_DIM == 2)
1692 {
1693 const unsigned num_nodes = pElement->GetNumNodes();
1694
1695 // Start at a random location in the vector of element nodes if requested
1696 const unsigned start_idx = randomOrder ? RandomNumberGenerator::Instance()->randMod(num_nodes) : 0u;
1697
1698 // Straighten out locations to a contiguous polygon rather than wrapped around the domain due to periodicity
1699 std::vector<c_vector<double, SPACE_DIM>> locations_straightened;
1700 locations_straightened.reserve(num_nodes);
1701
1702 locations_straightened.emplace_back(pElement->GetNodeLocation(start_idx));
1703
1704 for (unsigned node_idx = 1; node_idx < num_nodes; ++node_idx)
1705 {
1706 const unsigned prev_idx = AdvanceMod(start_idx, node_idx - 1, num_nodes);
1707 const unsigned this_idx = AdvanceMod(start_idx, node_idx, num_nodes);
1708
1709 const c_vector<double, SPACE_DIM>& r_last_location_added = locations_straightened.back();
1710
1711 const c_vector<double, SPACE_DIM>& r_prev_location = pElement->GetNodeLocation(prev_idx);
1712 const c_vector<double, SPACE_DIM>& r_this_location = pElement->GetNodeLocation(this_idx);
1713
1714 locations_straightened.emplace_back(r_last_location_added +
1715 this->GetVectorFromAtoB(r_prev_location, r_this_location));
1716 }
1717
1718 assert(locations_straightened.size() == num_nodes);
1719
1720 const bool closed_path = true;
1721 const bool permute_order = false;
1722 const std::size_t num_pts_to_place = num_nodes;
1723
1724 std::vector<c_vector<double, SPACE_DIM>> evenly_spaced_locations = EvenlySpaceAlongPath(
1725 locations_straightened,
1726 closed_path,
1727 permute_order,
1728 num_pts_to_place
1729 );
1730
1731 assert(evenly_spaced_locations.size() == num_nodes);
1732
1733 // Conform all locations to geometry
1734 for (c_vector<double, SPACE_DIM>& r_loc : evenly_spaced_locations)
1735 {
1736 ConformToGeometry(r_loc);
1737 }
1738
1739 // Update the node locations
1740 for (unsigned node_idx = 0; node_idx < num_nodes; ++node_idx)
1741 {
1742 const unsigned this_idx = AdvanceMod(node_idx, start_idx, num_nodes);
1743 pElement->GetNode(this_idx)->rGetModifiableLocation() = evenly_spaced_locations[node_idx];
1744 }
1745 }
1746 else
1747 {
1749 }
1750}
1751
1752template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1754 bool randomOrder)
1755{
1756 if constexpr (SPACE_DIM == 2)
1757 {
1758 const unsigned num_nodes = pLamina->GetNumNodes();
1759
1760 // Start at a random location in the vector of element nodes
1761 const unsigned start_idx = randomOrder ? RandomNumberGenerator::Instance()->randMod(num_nodes) : 0u;
1762
1763 // Straighten out locations to a contiguous line rather than wrapped around the domain due to periodicity
1764 std::vector<c_vector<double, SPACE_DIM>> locations_straightened;
1765 locations_straightened.reserve(1 + num_nodes);
1766
1767 locations_straightened.emplace_back(pLamina->GetNodeLocation(start_idx));
1768
1769 // We go to 1 + num_nodes so that we add on a node in a location congruent to the first location added
1770 for (unsigned node_idx = 1; node_idx < 1 + num_nodes; ++node_idx)
1771 {
1772 const unsigned prev_idx = AdvanceMod(start_idx, node_idx - 1, num_nodes);
1773 const unsigned this_idx = AdvanceMod(start_idx, node_idx, num_nodes);
1774
1775 const c_vector<double, SPACE_DIM>& r_last_location_added = locations_straightened.back();
1776
1777 const c_vector<double, SPACE_DIM>& r_prev_location = pLamina->GetNodeLocation(prev_idx);
1778 const c_vector<double, SPACE_DIM>& r_this_location = pLamina->GetNodeLocation(this_idx);
1779
1780 locations_straightened.emplace_back(r_last_location_added +
1781 this->GetVectorFromAtoB(r_prev_location, r_this_location));
1782 }
1783
1784 assert(locations_straightened.size() == 1 + num_nodes);
1785
1786 const bool closed_path = false;
1787 const bool permute_order = false;
1788 const std::size_t num_pts_to_place = 1 + num_nodes;
1789
1790 std::vector<c_vector<double, SPACE_DIM>> evenly_spaced_locations = EvenlySpaceAlongPath(
1791 locations_straightened,
1792 closed_path,
1793 permute_order,
1794 num_pts_to_place
1795 );
1796
1797 assert(evenly_spaced_locations.size() == 1 + num_nodes);
1798
1799 // Conform all locations to geometry
1800 for (c_vector<double, SPACE_DIM>& r_loc : evenly_spaced_locations)
1801 {
1802 ConformToGeometry(r_loc);
1803 }
1804
1805 // Update the node locations, ignoring the very last location that was added to make the path spacing even
1806 for (unsigned node_idx = 0; node_idx < num_nodes; ++node_idx)
1807 {
1808 const unsigned this_idx = AdvanceMod(start_idx, node_idx, num_nodes);
1809 pLamina->GetNode(this_idx)->rGetModifiableLocation() = evenly_spaced_locations[node_idx];
1810 }
1811 }
1812 else
1813 {
1815 }
1816}
1817
1818template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1820{
1821 for (unsigned dim = 0; dim < SPACE_DIM; ++dim)
1822 {
1823 assert(rLocation[dim] >= -2.0); // It's expected that the location is near the domain. This method is
1824 assert(rLocation[dim] < 3.0); // inefficient if the location is far away.
1825
1826 while (rLocation[dim] < 0.0)
1827 {
1828 rLocation[dim] += 1.0;
1829 }
1830
1831 while (rLocation[dim] >= 1.0)
1832 {
1833 rLocation[dim] -= 1.0;
1834 }
1835 }
1836}
1837
1838template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1840 Node<SPACE_DIM>* pNodeB)
1841{
1842 // If neither are lamina nodes, we can just check equality of first containing element indices
1843 if (pNodeA->GetRegion() != LAMINA_REGION && pNodeB->GetRegion() != LAMINA_REGION)
1844 {
1845 return *(pNodeA->ContainingElementsBegin()) != *(pNodeB->ContainingElementsBegin());
1846 }
1847 else // either one or both nodes is in lamina; assume that two laminas never interact
1848 {
1849 return true;
1850 }
1851}
1852
1853template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1855{
1856 if (!mLaminas.empty())
1857 {
1858 EXCEPTION("This method does not yet work in the presence of laminas");
1859 }
1860
1861 /*
1862 * This method uses geometric information from a voronoi diagram of every node. The voronoi cell of an element
1863 * (the union of voronoi cells of nodes in the element) provides a precise boundary between elements on which the
1864 * concept of 'neighbourhood' is well-defined.
1865 *
1866 * We first update the voronoi diagram, if it is out of date.
1867 */
1868 UpdateNodeLocationsVoronoiDiagramIfOutOfDate();
1869
1870 std::set<unsigned> neighbouring_element_indices;
1871
1872 for (unsigned node_local_idx = 0; node_local_idx < this->GetElement(elemIdx)->GetNumNodes(); ++node_local_idx)
1873 {
1874 const unsigned node_global_idx = this->GetElement(elemIdx)->GetNodeGlobalIndex(node_local_idx);
1875 const unsigned voronoi_cell_id = mVoronoiCellIdsIndexedByNodeIndex[node_global_idx];
1876
1877 /*
1878 * Iterate over the edges of the voronoi cell corresponding to the current node. Each primary edge has a twin
1879 * in a voronoi cell corresponding to a different node. The element containing that node (which may be the
1880 * current element under consideration) is added to the set of neighbours precisely if the distance between
1881 * nodes is less than mNeighbourDist.
1882 */
1883 const auto voronoi_cell = mNodeLocationsVoronoiDiagram.cells()[voronoi_cell_id];
1884 auto p_edge = voronoi_cell.incident_edge();
1885
1886 do
1887 {
1888 // The global node index corresponding to a voronoi cell cell is encoded in its 'color' variable
1889 const unsigned twin_node_idx = p_edge->twin()->cell()->color();
1890 const unsigned twin_elem_idx = *this->GetNode(twin_node_idx)->ContainingElementsBegin();
1891
1892 // Only bother to check the node distances if the nodes are in different elements
1893 if (twin_elem_idx != elemIdx)
1894 {
1895 if (this->GetDistanceBetweenNodes(node_global_idx, twin_node_idx) < mNeighbourDist)
1896 {
1897 neighbouring_element_indices.insert(twin_elem_idx);
1898 }
1899 }
1900 p_edge = p_edge->next();
1901 } while (p_edge != voronoi_cell.incident_edge());
1902 }
1903
1904 return neighbouring_element_indices;
1905} // LCOV_EXCL_LINE
1906
1907template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1909{
1910 if constexpr (SPACE_DIM == 2)
1911 {
1912 std::array<unsigned, 13> polygon_dist = {{0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u, 0u}};
1913
1914 for (auto elem_it = this->GetElementIteratorBegin(); elem_it != this->GetElementIteratorEnd(); ++elem_it)
1915 {
1916 if (!elem_it->IsElementOnBoundary())
1917 {
1918 // Accumulate all 12+ sided shapes
1919 unsigned num_neighbours = std::min<unsigned>(12u, GetNeighbouringElementIndices(elem_it->GetIndex()).size());
1920 polygon_dist[num_neighbours]++;
1921 }
1922 }
1923
1924 return polygon_dist;
1925 }
1926 else
1927 {
1929 }
1930}
1931
1932template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1933double ImmersedBoundaryMesh<ELEMENT_DIM, SPACE_DIM>::CalculateLengthOfVoronoiEdge(const boost::polygon::voronoi_diagram<double>::edge_type& rEdge)
1934{
1935 assert(rEdge.is_finite());
1936
1937 const double d_x = rEdge.vertex1()->x() - rEdge.vertex0()->x();
1938 const double d_y = rEdge.vertex1()->y() - rEdge.vertex0()->y();
1939
1940 return ScaleDistanceDownFromVoronoi(std::sqrt(d_x * d_x + d_y * d_y));
1941}
1942
1943template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1945{
1946 if (this->mNodes.begin() == this->mNodes.end())
1947 {
1948 return UINT_MAX;
1949 }
1950 else
1951 {
1952 const auto max_idx_it = std::max_element(this->mNodes.begin(), this->mNodes.end(),
1953 [](Node<SPACE_DIM>* const a, Node<SPACE_DIM>* const b)
1954 {
1955 return a->GetIndex() < b->GetIndex();
1956 });
1957
1958 return (*max_idx_it)->GetIndex();
1959 }
1960}
1961
1962template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1964{
1965 if (this->mElements.begin() == this->mElements.end())
1966 {
1967 return UINT_MAX;
1968 }
1969 else
1970 {
1972
1973 const auto max_idx_it = std::max_element(this->mElements.begin(), this->mElements.end(),
1974 [](IbElem* const a, IbElem* const b)
1975 {
1976 return a->GetIndex() < b->GetIndex();
1977 });
1978
1979 return (*max_idx_it)->GetIndex();
1980 }
1981}
1982
1983template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
1985{
1986 if (this->mLaminas.begin() == this->mLaminas.end())
1987 {
1988 return UINT_MAX;
1989 }
1990 else
1991 {
1992 using IbLam = ImmersedBoundaryElement<ELEMENT_DIM - 1, SPACE_DIM>;
1993
1994 const auto max_idx_it = std::max_element(this->mLaminas.begin(), this->mLaminas.end(),
1995 [](IbLam* const a, IbLam* const b)
1996 {
1997 return a->GetIndex() < b->GetIndex();
1998 });
1999
2000 return (*max_idx_it)->GetIndex();
2001 }
2002}
2003
2004template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2006{
2007 if (!mLaminas.empty())
2008 {
2009 EXCEPTION("This method does not yet work in the presence of laminas");
2010 }
2011
2012 // Need a large double, but since we also perform calculations with it, we can't really use DBL_MAX
2013 const double LARGE_DOUBLE = 1e6;
2014
2015 /*
2016 * This method uses geometric information from a voronoi diagram of every node.
2017 *
2018 * Nodes with infinite voronoi cells must be in a boundary element.
2019 * In addition, some nodes in boundary elements will have unusually long shared-edges, and the voronoi perimeter
2020 * of boundary elements will exceed the true perimeter by an unusually large amount.
2021 * We use a threshold on the median values of these identifiers to robustly determine boundary elements.
2022 */
2023 UpdateNodeLocationsVoronoiDiagramIfOutOfDate();
2024
2025 // Get the maximum element index in the mesh
2026 const unsigned max_elem_idx = GetMaxElementIndex();
2027
2028 // Keep track of all finite shared edge lengths (between two elements), as well as (for each element) the maximum
2029 // shared edge length.
2030 std::vector<double> max_shared_lengths(1 + max_elem_idx, 0.0);
2031 std::vector<double> voronoi_perimeter(1 + max_elem_idx, 0.0);
2032
2033 // Loop over nodes and calculate values to fill the containers defined above
2034 for (const auto& p_node : this->mNodes)
2035 {
2036 const unsigned this_node_idx = p_node->GetIndex();
2037 const unsigned this_elem_idx = *p_node->ContainingElementsBegin();
2038
2039 // Get the voronoi cell that corresponds to this node
2040 const unsigned voronoi_cell_id = mVoronoiCellIdsIndexedByNodeIndex[this_node_idx];
2041 const auto& voronoi_cell = mNodeLocationsVoronoiDiagram.cells()[voronoi_cell_id];
2042
2043 // Iterate over the edges of this cell. Note that due to the halo region none of these edges should be infinite.
2044 auto p_edge = voronoi_cell.incident_edge();
2045
2046 // Loop over the voronoi edges associated with this node's voronoi cell
2047 do
2048 {
2049 // If the edge is infinite, the element containing it must be on the boundary
2050 if (p_edge->is_infinite())
2051 {
2052 max_shared_lengths[this_elem_idx] = LARGE_DOUBLE;
2053 voronoi_perimeter[this_elem_idx] += LARGE_DOUBLE;
2054 }
2055 else
2056 {
2057 // The global node index corresponding to a voronoi cell cell is encoded in its 'color' variable
2058 const unsigned twin_node_idx = p_edge->twin()->cell()->color();
2059 const unsigned twin_elem_idx = *this->GetNode(twin_node_idx)->ContainingElementsBegin();
2060
2061 // If the edge is between two elements, it counts
2062 if (this_elem_idx != twin_elem_idx)
2063 {
2064 const double edge_length = CalculateLengthOfVoronoiEdge(*p_edge);
2065 max_shared_lengths[this_elem_idx] = std::max(max_shared_lengths[this_elem_idx], edge_length);
2066 voronoi_perimeter[this_elem_idx] += edge_length;
2067 }
2068 }
2069 p_edge = p_edge->next();
2070 } while (p_edge != voronoi_cell.incident_edge());
2071 }
2072
2073 // Calculate what proportion bigger the voronoi perimeter is than the actual one
2074 std::vector<double> perimeter_multiples(voronoi_perimeter.size());
2075 for (const auto& p_elem : this->mElements)
2076 {
2077 const unsigned idx = p_elem->GetIndex();
2078 perimeter_multiples[idx] = voronoi_perimeter[idx] / this->GetSurfaceAreaOfElement(idx);
2079 }
2080
2081 // Calculate the medians of each vector. Copy the vector accounting for possible non-contiguous element indices.
2082 std::vector<double> copy_max_shared_lengths;
2083 std::vector<double> copy_perimeter_multiples;
2084
2085 for (const auto& p_elem : this->mElements)
2086 {
2087 const unsigned idx = p_elem->GetIndex();
2088 copy_max_shared_lengths.emplace_back(max_shared_lengths[idx]);
2089 copy_perimeter_multiples.emplace_back(perimeter_multiples[idx]);
2090 }
2091
2092 const std::size_t half_way = copy_max_shared_lengths.size() / 2;
2093 assert(half_way == copy_perimeter_multiples.size() / 2);
2094
2095 std::nth_element(copy_max_shared_lengths.begin(), copy_max_shared_lengths.begin() + half_way, copy_max_shared_lengths.end());
2096 std::nth_element(copy_perimeter_multiples.begin(), copy_perimeter_multiples.begin() + half_way, copy_perimeter_multiples.end());
2097
2098 double median_max_edge_length = copy_max_shared_lengths[half_way];
2099 double median_perimeter_multiple = copy_perimeter_multiples[half_way];
2100
2101 // In the event that most elements have infinite edges, for instance if there are a small number of elements,
2102 // a reduced median can be used
2103 if (median_max_edge_length == LARGE_DOUBLE || median_perimeter_multiple >= LARGE_DOUBLE)
2104 {
2105 median_max_edge_length = 0.5 * LARGE_DOUBLE;
2106 median_perimeter_multiple = 0.5 * LARGE_DOUBLE;
2107 }
2108
2109 // Finally, tag the boundary elements
2110 for (const auto& p_elem : this->mElements)
2111 {
2112 const unsigned idx = p_elem->GetIndex();
2113
2114 const bool large_shared_edge = max_shared_lengths[idx] > 1.1 * median_max_edge_length;
2115 const bool large_perimeter = perimeter_multiples[idx] > 1.1 * median_perimeter_multiple;
2116
2117 p_elem->SetIsBoundaryElement(large_shared_edge && large_perimeter);
2118 }
2119}
2120
2121template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2123{
2124 if constexpr (SPACE_DIM == 2)
2125 {
2126 using boost_point = boost::polygon::point_data<int>;
2127
2128 /*
2129 * The voronoi diagram needs updating only if the node locations have changed (i.e. not more than once per
2130 * timestep). We test a chosen summary of node locations against the cached mSummaryOfNodeLocations to determine
2131 * this.
2132 */
2133 double new_location_summary = this->mNodes.front()->rGetLocation()[0] +
2134 this->mNodes.front()->rGetLocation()[1] +
2135 this->mNodes.back()->rGetLocation()[0] +
2136 this->mNodes.back()->rGetLocation()[1];
2137
2138 bool voronoi_needs_updating = std::fabs(mSummaryOfNodeLocations - new_location_summary) > DBL_EPSILON;
2139
2140 /*
2141 * Method outline:
2142 * - Add all node locations as boost points, correctly scaled
2143 * - Add extra locations for any nodes within distance mVoronoiHalo of the domain edge for periodicity
2144 * - Calculate the voronoi diagram
2145 * - Populate mVoronoiCellIdsIndexedByNodeIndex for efficient mapping between node index and corresponding voronoi
2146 * cell
2147 * - Store the corresponding node index in each voronoi cell's "color" variable for efficient reverse-lookup
2148 * - Tag boundary elements
2149 */
2150 if (voronoi_needs_updating)
2151 {
2152 mSummaryOfNodeLocations = new_location_summary;
2153
2154 // Helper c_vectors for halo region
2155 const c_vector<double, SPACE_DIM> halo_up = Create_c_vector(0.0, 1.0);
2156 const c_vector<double, SPACE_DIM> halo_down = Create_c_vector(0.0, -1.0);
2157 const c_vector<double, SPACE_DIM> halo_left = Create_c_vector(-1.0, 0.0);
2158 const c_vector<double, SPACE_DIM> halo_right = Create_c_vector(1.0, 0.0);
2159
2160 std::vector<std::pair<unsigned, c_vector<double, SPACE_DIM>>> halo_ids_and_locations;
2161 std::vector<unsigned> node_ids_in_source_idx_order;
2162
2163 // We need to translate node locations into boost points, which are scaled to an integer grid
2164 std::vector<boost_point> points;
2165
2166 // First add a point for each node, and calculate any additional locations needed within the halo
2167 for (const auto& p_node : this->mNodes)
2168 {
2169 const double x_pos = p_node->rGetLocation()[0];
2170 const double y_pos = p_node->rGetLocation()[1];
2171
2172 // Put the integer voronoi coordinate in for this node's location
2173 points.emplace_back(boost_point(ScaleUpToVoronoiCoordinate(x_pos), ScaleUpToVoronoiCoordinate(y_pos)));
2174 node_ids_in_source_idx_order.emplace_back(p_node->GetIndex());
2175
2176 // Now, for this location, decide whether any copies of this location are needed in the halo region
2177 const bool needed_up = y_pos < mVoronoiHalo;
2178 const bool needed_down = y_pos > 1.0 - mVoronoiHalo;
2179 const bool needed_left = x_pos > 1.0 - mVoronoiHalo;
2180 const bool needed_right = x_pos < mVoronoiHalo;
2181
2182 if (needed_up)
2183 {
2184 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2185 p_node->rGetLocation() + halo_up));
2186 }
2187 if (needed_down)
2188 {
2189 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2190 p_node->rGetLocation() + halo_down));
2191 }
2192 if (needed_left)
2193 {
2194 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2195 p_node->rGetLocation() + halo_left));
2196 }
2197 if (needed_right)
2198 {
2199 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2200 p_node->rGetLocation() + halo_right));
2201 }
2202 if (needed_up && needed_left)
2203 {
2204 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2205 p_node->rGetLocation() + halo_up + halo_left));
2206 }
2207 if (needed_up && needed_right)
2208 {
2209 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2210 p_node->rGetLocation() + halo_up + halo_right));
2211 }
2212 if (needed_down && needed_left)
2213 {
2214 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2215 p_node->rGetLocation() + halo_down + halo_left));
2216 }
2217 if (needed_down && needed_right)
2218 {
2219 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2220 p_node->rGetLocation() + halo_down + halo_right));
2221 }
2222 }
2223
2224 // Next add the additional points
2225 for (const auto& pair : halo_ids_and_locations)
2226 {
2227 const unsigned node_idx = pair.first;
2228 c_vector<double, SPACE_DIM> location = pair.second;
2229
2230 const int x_coord = ScaleUpToVoronoiCoordinate(location[0]);
2231 const int y_coord = ScaleUpToVoronoiCoordinate(location[1]);
2232
2233 points.emplace_back(boost_point(x_coord, y_coord));
2234 node_ids_in_source_idx_order.emplace_back(node_idx);
2235 }
2236
2237 // Construct the voronoi diagram. This is the costly part of this method.
2238 mNodeLocationsVoronoiDiagram.clear();
2239 construct_voronoi(std::begin(points), std::end(points), &mNodeLocationsVoronoiDiagram);
2240
2241 // We need an efficient map from node global index to the voronoi cell representing it, and we can't assume
2242 // that the nodes are ordered sequentially. We first identify the largest node index.
2243 const unsigned max_node_idx = GetMaxNodeIndex();
2244
2245 mVoronoiCellIdsIndexedByNodeIndex.resize(1 + max_node_idx);
2246
2247 for (unsigned vor_cell_id = 0; vor_cell_id < mNodeLocationsVoronoiDiagram.cells().size(); ++vor_cell_id)
2248 {
2249 // Source index is incrementally given to each input point, which is in order of nodes in this->mNodes.
2250 // We need to be able to identify each vor_cell_id by the global node index
2251
2252 auto& r_this_cell = mNodeLocationsVoronoiDiagram.cells()[vor_cell_id];
2253 const auto source_idx = r_this_cell.source_index();
2254 const unsigned node_idx = node_ids_in_source_idx_order[source_idx];
2255
2256 r_this_cell.color(node_idx);
2257
2258 if (source_idx < this->mNodes.size())
2259 {
2260 mVoronoiCellIdsIndexedByNodeIndex[node_idx] = vor_cell_id;
2261 }
2262 }
2263
2264 // Finally, if the diagram was out-of-date, we will need to re-tag boundary elements.
2265 this->TagBoundaryElements();
2266 }
2267 }
2268 else
2269 {
2271 }
2272}
2273
2274template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2276{
2277 assert(location >= -mVoronoiHalo);
2278 assert(location <= 1.0 + mVoronoiHalo);
2279
2280 constexpr auto DBL_INT_MIN = static_cast<double>(INT_MIN);
2281 constexpr auto DBL_INT_RANGE = static_cast<double>(INT_MAX) - static_cast<double>(INT_MIN);
2282
2283 // To handle periodicity, we add (again) any nodes within a halo region of the unit square
2284 constexpr double scale_factor = DBL_INT_RANGE / (1.0 + 2.0 * mVoronoiHalo);
2285
2286 // By construction there is no narrowing, so this static_cast might not be necessary
2287 return static_cast<int>(std::lround(DBL_INT_MIN + (mVoronoiHalo + location) * scale_factor));
2288}
2289
2290template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
2292{
2293 constexpr auto DBL_INT_RANGE = static_cast<double>(INT_MAX) - static_cast<double>(INT_MIN);
2294 constexpr double scale_factor = (1.0 + 2.0 * mVoronoiHalo) / DBL_INT_RANGE;
2295
2296 return scale_factor * distance;
2297}
2298
2299template<unsigned int ELEMENT_DIM, unsigned int SPACE_DIM>
2300const boost::polygon::voronoi_diagram<double>& ImmersedBoundaryMesh<ELEMENT_DIM, SPACE_DIM>::rGetNodeLocationsVoronoiDiagram(bool update)
2301{
2302 if (update)
2303 {
2304 this->UpdateNodeLocationsVoronoiDiagramIfOutOfDate();
2305 }
2306
2307 return mNodeLocationsVoronoiDiagram;
2308}
2309
2310template<unsigned int ELEMENT_DIM, unsigned int SPACE_DIM>
2312{
2313 return mVoronoiCellIdsIndexedByNodeIndex;
2314}
2315
2316// Explicit instantiation
2317template class ImmersedBoundaryMesh<1, 1>;
2318template class ImmersedBoundaryMesh<1, 2>;
2319template class ImmersedBoundaryMesh<1, 3>;
2320template class ImmersedBoundaryMesh<2, 2>;
2321template class ImmersedBoundaryMesh<2, 3>;
2322template class ImmersedBoundaryMesh<3, 3>;
2323
2324// Serialization for Boost >= 1.36
const double DOUBLE_UNSET
Definition Exception.hpp:57
#define EXCEPTION(message)
#define NEVER_REACHED
#define EXPORT_TEMPLATE_CLASS_ALL_DIMS(CLASS)
void SetAttribute(double attribute)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
std::vector< double > & rGetElementAttributes()
double GetNodeLocation(unsigned localIndex, unsigned dimension) const
unsigned GetNumElementAttributes()
unsigned GetNumNodes() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
unsigned GetIndex() const
virtual bool HasNodePermutation()
bool mMeshChangesDuringSimulation
std::vector< Node< SPACE_DIM > * > mNodes
std::shared_ptr< FluidSource< SPACE_DIM > > GetFluidSource()
std::vector< Node< SPACE_DIM > * > & rGetCornerNodes()
void SetIsBoundaryElement(bool isBoundaryElement)
ImmersedBoundaryElementData GetNextImmersedBoundaryElementData()
ImmersedBoundaryElementData GetNextImmersedBoundaryLaminaData()
ChasteCuboid< SPACE_DIM > CalculateBoundingBoxOfElement(unsigned index)
std::vector< ImmersedBoundaryElement< ELEMENT_DIM, SPACE_DIM > * > mElements
void ReMesh(bool randomOrder=false)
virtual double GetVolumeOfElement(unsigned index)
double GetCharacteristicNodeSpacing() const
void SetNeighbourDist(double neighbourDist)
bool NodesInDifferentElementOrLamina(Node< SPACE_DIM > *pNodeA, Node< SPACE_DIM > *pNodeB)
unsigned DivideElement(ImmersedBoundaryElement< ELEMENT_DIM, SPACE_DIM > *pElement, unsigned nodeAIndex, unsigned nodeBIndex, c_vector< double, SPACE_DIM > centroid, c_vector< double, SPACE_DIM > axisOfDivision)
const std::vector< Node< SPACE_DIM > * > & rGetNodes() const
unsigned GetMaxElementIndex() const
multi_array< double, 3 > & rGetModifiable2dVelocityGrids()
unsigned SolveElementMapping(unsigned index) const
ImmersedBoundaryElement< ELEMENT_DIM - 1, SPACE_DIM > * GetLamina(unsigned index) const
void ReMeshLamina(ImmersedBoundaryElement< ELEMENT_DIM - 1, SPACE_DIM > *pLamina, bool randomOrder)
void SetNode(unsigned nodeIndex, ChastePoint< SPACE_DIM > point)
double GetSkewnessOfElementMassDistributionAboutAxis(unsigned elemIndex, c_vector< double, SPACE_DIM > axis)
void SetNumGridPtsXAndY(unsigned numGridPts)
double ScaleDistanceDownFromVoronoi(const double distance) const
void SetCharacteristicNodeSpacing(double nodeSpacing)
std::array< unsigned, 13 > GetPolygonDistribution()
double CalculateLengthOfVoronoiEdge(const boost::polygon::voronoi_diagram< double >::edge_type &rEdge)
void SetCellRearrangementThreshold(double cellRearrangementThreshold)
std::vector< std::shared_ptr< FluidSource< SPACE_DIM > > > mElementFluidSources
std::vector< std::shared_ptr< FluidSource< SPACE_DIM > > > & rGetElementFluidSources()
virtual unsigned GetNumNodes() const
double GetElongationShapeFactorOfElement(unsigned elementIndex)
unsigned SolveNodeMapping(unsigned index) const
std::vector< std::shared_ptr< FluidSource< SPACE_DIM > > > mBalancingFluidSources
unsigned AddNode(Node< SPACE_DIM > *pNewNode)
double GetAverageNodeSpacingOfElement(unsigned index, bool recalculate=true)
virtual c_vector< double, 3 > CalculateMomentsOfElement(unsigned index)
void ConformToGeometry(c_vector< double, SPACE_DIM > &rLocation)
const boost::polygon::voronoi_diagram< double > & rGetNodeLocationsVoronoiDiagram(bool update=true)
void ReMeshElement(ImmersedBoundaryElement< ELEMENT_DIM, SPACE_DIM > *pElement, bool randomOrder)
multi_array< double, 3 > m2dVelocityGrids
std::set< unsigned > GetNeighbouringNodeIndices(unsigned nodeIndex)
unsigned DivideElementAlongGivenAxis(ImmersedBoundaryElement< ELEMENT_DIM, SPACE_DIM > *pElement, c_vector< double, SPACE_DIM > axisOfDivision, bool placeOriginalElementBelow=false)
virtual unsigned GetNumElements() const
void SetNumGridPtsX(unsigned meshPointsX)
std::vector< std::shared_ptr< FluidSource< SPACE_DIM > > > & rGetBalancingFluidSources()
double GetAverageNodeSpacingOfLamina(unsigned index, bool recalculate=true)
int ScaleUpToVoronoiCoordinate(double location) const
unsigned SolveBoundaryElementMapping(unsigned index) const
unsigned DivideElementAlongShortAxis(ImmersedBoundaryElement< ELEMENT_DIM, SPACE_DIM > *pElement, bool placeOriginalElementBelow=false)
virtual c_vector< double, SPACE_DIM > GetCentroidOfElement(unsigned index)
void SetElementDivisionSpacing(double elementDivisionSpacing)
c_vector< double, SPACE_DIM > GetShortAxisOfElement(unsigned index)
std::set< unsigned > GetNeighbouringElementIndices(unsigned elemIdx)
const multi_array< double, 3 > & rGet2dVelocityGrids() const
std::vector< ImmersedBoundaryElement< ELEMENT_DIM - 1, SPACE_DIM > * > mLaminas
void ConstructFromMeshReader(AbstractMeshReader< ELEMENT_DIM, SPACE_DIM > &rMeshReader)
ImmersedBoundaryElement< ELEMENT_DIM, SPACE_DIM > * GetElement(unsigned index) const
double GetVoronoiSurfaceAreaOfElement(unsigned elemIdx)
virtual double GetSurfaceAreaOfElement(unsigned index)
const std::vector< unsigned int > & GetVoronoiCellIdsIndexedByNodeIndex() const
c_vector< double, SPACE_DIM > GetVectorFromAtoB(const c_vector< double, SPACE_DIM > &rLocation1, const c_vector< double, SPACE_DIM > &rLocation2)
void SetNumGridPtsY(unsigned meshPointsY)
unsigned GetNumAllElements() const
Definition Node.hpp:59
void SetIndex(unsigned index)
Definition Node.cpp:121
ContainingElementIterator ContainingElementsBegin() const
Definition Node.hpp:485
unsigned GetRegion() const
Definition Node.cpp:437
static RandomNumberGenerator * Instance()
unsigned randMod(unsigned base)