793 [[maybe_unused]]
unsigned index)
796 if constexpr (SPACE_DIM == 2)
801 c_vector<double, SPACE_DIM> centroid = zero_vector<double>(SPACE_DIM);
803 double centroid_x = 0;
804 double centroid_y = 0;
807 double element_signed_area = 0.0;
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);
814 for (
unsigned local_index = 0; local_index < num_nodes; local_index++)
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);
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];
824 double signed_area_term = this_x * next_y - this_y * next_x;
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;
833 assert(element_signed_area != 0.0);
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);
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);
1370 unsigned nodeAIndex,
1371 unsigned nodeBIndex,
1372 c_vector<double, SPACE_DIM> centroid,
1373 c_vector<double, SPACE_DIM> axisOfDivision)
1375 if constexpr (SPACE_DIM == 2 && ELEMENT_DIM == 2)
1379 EXCEPTION(
"The value of mElementDivisionSpacing has not been set.");
1395 double half_spacing = 0.5 * mElementDivisionSpacing;
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];
1415 unsigned start_a = (nodeAIndex + 1) % num_nodes;
1416 unsigned end_a = nodeBIndex;
1418 unsigned start_b = (nodeBIndex + 1) % num_nodes;
1419 unsigned end_b = nodeAIndex;
1422 bool no_node_satisfied_condition_1 =
true;
1423 for (
unsigned i = start_a; i != end_a;)
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);
1428 if (fabs(perpendicular_dist) >= half_spacing)
1430 no_node_satisfied_condition_1 =
false;
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);
1442 i = (i + 1) % num_nodes;
1446 bool no_node_satisfied_condition_2 =
true;
1447 for (
unsigned i = end_a; i != start_a;)
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);
1452 if (fabs(perpendicular_dist) >= half_spacing)
1454 no_node_satisfied_condition_2 =
false;
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);
1466 i = (i + num_nodes - 1) % num_nodes;
1470 bool no_node_satisfied_condition_3 =
true;
1471 for (
unsigned i = start_b; i != end_b;)
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);
1476 if (fabs(perpendicular_dist) >= half_spacing)
1478 no_node_satisfied_condition_3 =
false;
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);
1490 i = (i + 1) % num_nodes;
1494 bool no_node_satisfied_condition_4 =
true;
1495 for (
unsigned i = end_b; i != start_b;)
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);
1500 if (fabs(perpendicular_dist) >= half_spacing)
1502 no_node_satisfied_condition_4 =
false;
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);
1514 i = (i + num_nodes - 1) % num_nodes;
1517 if (no_node_satisfied_condition_1 || no_node_satisfied_condition_2 || no_node_satisfied_condition_3 || no_node_satisfied_condition_4)
1519 EXCEPTION(
"Could not space elements far enough apart during cell division. Cannot currently handle this case");
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;)
1528 daughter_a_location_stencil.push_back(c_vector<double, SPACE_DIM>(pElement->
GetNode(node_idx)->rGetLocation()));
1531 node_idx = (node_idx + 1) % num_nodes;
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;)
1537 daughter_b_location_stencil.push_back(c_vector<double, SPACE_DIM>(pElement->
GetNode(node_idx)->rGetLocation()));
1540 node_idx = (node_idx + 1) % num_nodes;
1543 assert(daughter_a_location_stencil.size() > 1);
1544 assert(daughter_b_location_stencil.size() > 1);
1547 daughter_a_location_stencil.push_back(daughter_a_location_stencil[0]);
1548 daughter_b_location_stencil.push_back(daughter_b_location_stencil[0]);
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++)
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])));
1559 for (
unsigned loc_idx = 1; loc_idx < daughter_b_location_stencil.size(); loc_idx++)
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])));
1565 double target_spacing_a = cumulative_distances_a.back() / (
double)num_nodes;
1566 double target_spacing_b = cumulative_distances_b.back() / (
double)num_nodes;
1569 unsigned last_idx_used = 0;
1570 for (
unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
1572 double location_along_arc = (
double)node_idx * target_spacing_a;
1574 while (location_along_arc > cumulative_distances_a[last_idx_used + 1])
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]);
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]);
1585 c_vector<double, SPACE_DIM> new_location_a = daughter_a_location_stencil[last_idx_used] + interpolant * this_to_next;
1592 std::vector<Node<SPACE_DIM>*> new_nodes_vec;
1593 for (
unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
1595 double location_along_arc = (
double)node_idx * target_spacing_b;
1597 while (location_along_arc > cumulative_distances_b[last_idx_used + 1])
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]);
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]);
1608 c_vector<double, SPACE_DIM> new_location_b = daughter_b_location_stencil[last_idx_used] + interpolant * this_to_next;
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());
1616 for (
unsigned node_idx = 0; node_idx < num_nodes; node_idx++)
1618 new_nodes_vec[node_idx]->SetRegion(pElement->
GetNode(node_idx)->GetRegion());
1620 for (
unsigned node_attribute = 0; node_attribute < pElement->
GetNode(node_idx)->GetNumNodeAttributes(); node_attribute++)
1622 new_nodes_vec[node_idx]->AddNodeAttribute(pElement->
GetNode(node_idx)->rGetNodeAttributes()[node_attribute]);
1627 unsigned new_elem_idx = this->mElements.size();
1629 this->mElements.back()->RegisterWithNodes();
1638 for (
unsigned corner = 0; corner < pElement->
rGetCornerNodes().size(); corner++)
1640 this->mElements.back()->rGetCornerNodes().push_back(pElement->
rGetCornerNodes()[corner]);
1647 c_vector<double, SPACE_DIM> new_centroid = this->GetCentroidOfElement(new_elem_idx);
1651 mElementFluidSources.back()->SetAssociatedElementIndex(new_elem_idx);
1652 mElementFluidSources.back()->SetStrength(0.0);
1655 mElements[new_elem_idx]->SetFluidSource(mElementFluidSources.back());
1657 return new_elem_idx;
1691 if constexpr (SPACE_DIM == 2)
1693 const unsigned num_nodes = pElement->
GetNumNodes();
1699 std::vector<c_vector<double, SPACE_DIM>> locations_straightened;
1700 locations_straightened.reserve(num_nodes);
1702 locations_straightened.emplace_back(pElement->
GetNodeLocation(start_idx));
1704 for (
unsigned node_idx = 1; node_idx < num_nodes; ++node_idx)
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);
1709 const c_vector<double, SPACE_DIM>& r_last_location_added = locations_straightened.back();
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);
1714 locations_straightened.emplace_back(r_last_location_added +
1715 this->GetVectorFromAtoB(r_prev_location, r_this_location));
1718 assert(locations_straightened.size() == num_nodes);
1720 const bool closed_path =
true;
1721 const bool permute_order =
false;
1722 const std::size_t num_pts_to_place = num_nodes;
1724 std::vector<c_vector<double, SPACE_DIM>> evenly_spaced_locations = EvenlySpaceAlongPath(
1725 locations_straightened,
1731 assert(evenly_spaced_locations.size() == num_nodes);
1734 for (c_vector<double, SPACE_DIM>& r_loc : evenly_spaced_locations)
1736 ConformToGeometry(r_loc);
1740 for (
unsigned node_idx = 0; node_idx < num_nodes; ++node_idx)
1742 const unsigned this_idx = AdvanceMod(node_idx, start_idx, num_nodes);
1743 pElement->
GetNode(this_idx)->rGetModifiableLocation() = evenly_spaced_locations[node_idx];
1756 if constexpr (SPACE_DIM == 2)
1758 const unsigned num_nodes = pLamina->
GetNumNodes();
1764 std::vector<c_vector<double, SPACE_DIM>> locations_straightened;
1765 locations_straightened.reserve(1 + num_nodes);
1767 locations_straightened.emplace_back(pLamina->
GetNodeLocation(start_idx));
1770 for (
unsigned node_idx = 1; node_idx < 1 + num_nodes; ++node_idx)
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);
1775 const c_vector<double, SPACE_DIM>& r_last_location_added = locations_straightened.back();
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);
1780 locations_straightened.emplace_back(r_last_location_added +
1781 this->GetVectorFromAtoB(r_prev_location, r_this_location));
1784 assert(locations_straightened.size() == 1 + num_nodes);
1786 const bool closed_path =
false;
1787 const bool permute_order =
false;
1788 const std::size_t num_pts_to_place = 1 + num_nodes;
1790 std::vector<c_vector<double, SPACE_DIM>> evenly_spaced_locations = EvenlySpaceAlongPath(
1791 locations_straightened,
1797 assert(evenly_spaced_locations.size() == 1 + num_nodes);
1800 for (c_vector<double, SPACE_DIM>& r_loc : evenly_spaced_locations)
1802 ConformToGeometry(r_loc);
1806 for (
unsigned node_idx = 0; node_idx < num_nodes; ++node_idx)
1808 const unsigned this_idx = AdvanceMod(start_idx, node_idx, num_nodes);
1809 pLamina->
GetNode(this_idx)->rGetModifiableLocation() = evenly_spaced_locations[node_idx];
2007 if (!mLaminas.empty())
2009 EXCEPTION(
"This method does not yet work in the presence of laminas");
2013 const double LARGE_DOUBLE = 1e6;
2023 UpdateNodeLocationsVoronoiDiagramIfOutOfDate();
2026 const unsigned max_elem_idx = GetMaxElementIndex();
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);
2034 for (
const auto& p_node : this->mNodes)
2036 const unsigned this_node_idx = p_node->GetIndex();
2037 const unsigned this_elem_idx = *p_node->ContainingElementsBegin();
2040 const unsigned voronoi_cell_id = mVoronoiCellIdsIndexedByNodeIndex[this_node_idx];
2041 const auto& voronoi_cell = mNodeLocationsVoronoiDiagram.cells()[voronoi_cell_id];
2044 auto p_edge = voronoi_cell.incident_edge();
2050 if (p_edge->is_infinite())
2052 max_shared_lengths[this_elem_idx] = LARGE_DOUBLE;
2053 voronoi_perimeter[this_elem_idx] += LARGE_DOUBLE;
2058 const unsigned twin_node_idx = p_edge->twin()->cell()->color();
2059 const unsigned twin_elem_idx = *this->GetNode(twin_node_idx)->ContainingElementsBegin();
2062 if (this_elem_idx != twin_elem_idx)
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;
2069 p_edge = p_edge->next();
2070 }
while (p_edge != voronoi_cell.incident_edge());
2074 std::vector<double> perimeter_multiples(voronoi_perimeter.size());
2075 for (
const auto& p_elem : this->mElements)
2077 const unsigned idx = p_elem->
GetIndex();
2078 perimeter_multiples[idx] = voronoi_perimeter[idx] / this->GetSurfaceAreaOfElement(idx);
2082 std::vector<double> copy_max_shared_lengths;
2083 std::vector<double> copy_perimeter_multiples;
2085 for (
const auto& p_elem : this->mElements)
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]);
2092 const std::size_t half_way = copy_max_shared_lengths.size() / 2;
2093 assert(half_way == copy_perimeter_multiples.size() / 2);
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());
2098 double median_max_edge_length = copy_max_shared_lengths[half_way];
2099 double median_perimeter_multiple = copy_perimeter_multiples[half_way];
2103 if (median_max_edge_length == LARGE_DOUBLE || median_perimeter_multiple >= LARGE_DOUBLE)
2105 median_max_edge_length = 0.5 * LARGE_DOUBLE;
2106 median_perimeter_multiple = 0.5 * LARGE_DOUBLE;
2110 for (
const auto& p_elem : this->mElements)
2112 const unsigned idx = p_elem->
GetIndex();
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;
2124 if constexpr (SPACE_DIM == 2)
2126 using boost_point = boost::polygon::point_data<int>;
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];
2138 bool voronoi_needs_updating = std::fabs(mSummaryOfNodeLocations - new_location_summary) > DBL_EPSILON;
2150 if (voronoi_needs_updating)
2152 mSummaryOfNodeLocations = new_location_summary;
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);
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;
2164 std::vector<boost_point> points;
2167 for (
const auto& p_node : this->mNodes)
2169 const double x_pos = p_node->rGetLocation()[0];
2170 const double y_pos = p_node->rGetLocation()[1];
2173 points.emplace_back(boost_point(ScaleUpToVoronoiCoordinate(x_pos), ScaleUpToVoronoiCoordinate(y_pos)));
2174 node_ids_in_source_idx_order.emplace_back(p_node->GetIndex());
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;
2184 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2185 p_node->rGetLocation() + halo_up));
2189 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2190 p_node->rGetLocation() + halo_down));
2194 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2195 p_node->rGetLocation() + halo_left));
2199 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2200 p_node->rGetLocation() + halo_right));
2202 if (needed_up && needed_left)
2204 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2205 p_node->rGetLocation() + halo_up + halo_left));
2207 if (needed_up && needed_right)
2209 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2210 p_node->rGetLocation() + halo_up + halo_right));
2212 if (needed_down && needed_left)
2214 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2215 p_node->rGetLocation() + halo_down + halo_left));
2217 if (needed_down && needed_right)
2219 halo_ids_and_locations.emplace_back(std::make_pair(p_node->GetIndex(),
2220 p_node->rGetLocation() + halo_down + halo_right));
2225 for (
const auto& pair : halo_ids_and_locations)
2227 const unsigned node_idx = pair.first;
2228 c_vector<double, SPACE_DIM> location = pair.second;
2230 const int x_coord = ScaleUpToVoronoiCoordinate(location[0]);
2231 const int y_coord = ScaleUpToVoronoiCoordinate(location[1]);
2233 points.emplace_back(boost_point(x_coord, y_coord));
2234 node_ids_in_source_idx_order.emplace_back(node_idx);
2238 mNodeLocationsVoronoiDiagram.clear();
2239 construct_voronoi(std::begin(points), std::end(points), &mNodeLocationsVoronoiDiagram);
2243 const unsigned max_node_idx = GetMaxNodeIndex();
2245 mVoronoiCellIdsIndexedByNodeIndex.resize(1 + max_node_idx);
2247 for (
unsigned vor_cell_id = 0; vor_cell_id < mNodeLocationsVoronoiDiagram.cells().size(); ++vor_cell_id)
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];
2256 r_this_cell.color(node_idx);
2258 if (source_idx < this->mNodes.size())
2260 mVoronoiCellIdsIndexedByNodeIndex[node_idx] = vor_cell_id;
2265 this->TagBoundaryElements();