diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 81549c4039487f4175e920b62a7a54a77be2c6e6..1cdc1a23b016e6744c953d3e74a159f24499ea5d 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -18,6 +18,7 @@
#include <mesh/MeshRandomizer.hpp>
#include <mesh/MeshSmoother.hpp>
#include <mesh/MeshSmootherEscobar.hpp>
+#include <mesh/MeshSmootherJun.hpp>
#include <scheme/AcousticSolver.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
@@ -307,6 +308,60 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("smoothMeshJun",
+ std::function(
+
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::shared_ptr<const IMesh> {
+ MeshSmootherJunHandler handler;
+ return handler.getSmoothedMesh(p_mesh, bc_descriptor_list);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("smoothMeshJun",
+ std::function(
+
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const FunctionSymbolId& function_symbol_id) -> std::shared_ptr<const IMesh> {
+ MeshSmootherJunHandler handler;
+ return handler.getSmoothedMesh(p_mesh, bc_descriptor_list, function_symbol_id);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("smoothMeshJun",
+ std::function(
+
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list)
+ -> std::shared_ptr<const IMesh> {
+ MeshSmootherJunHandler handler;
+ return handler.getSmoothedMesh(p_mesh, bc_descriptor_list, smoothing_zone_list);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("smoothMeshJun",
+ std::function(
+
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
+ discrete_function_variant_list) -> std::shared_ptr<const IMesh> {
+ MeshSmootherJunHandler handler;
+ return handler.getSmoothedMesh(p_mesh, bc_descriptor_list,
+ discrete_function_variant_list);
+ }
+
+ ));
+
this->_addBuiltinFunction("smoothMeshEscobar",
std::function(
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 54c7afbf2f19799524faf0038e944d5463e1f7c1..63b1b000484f10031214f93d605db9d1c869a1a1 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -40,5 +40,6 @@ add_library(
MeshRandomizer.cpp
MeshSmoother.cpp
MeshSmootherEscobar.cpp
+ MeshSmootherJun.cpp
MeshTransformer.cpp
SynchronizerManager.cpp)
diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp
index eca67723dd72c733107cdd5b0caeef04ae24e117..0f32fa6417dae941b4480cd98bbf3c4ce43a1fef 100644
--- a/src/mesh/MeshSmootherEscobar.cpp
+++ b/src/mesh/MeshSmootherEscobar.cpp
@@ -369,7 +369,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const NodeId node_id = node_list[i_node];
if (m_is_smoothed_node[node_id] and node_is_owned[node_id]) {
smooth(node_id, new_xr[node_id]);
- // new_xr[node_id] = 0.5 * (new_xr[node_id] + old_xr[node_id]);
+ new_xr[node_id] = 0.5 * (new_xr[node_id] + old_xr[node_id]);
}
});
@@ -495,22 +495,18 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const double sin_alpha = std::sin(alpha);
const double cos_alpha = std::cos(alpha);
const double inv_sin_alpha = 1. / sin_alpha;
- const double inv_tan_alpha = 1. / std::tan(alpha);
+ const double inv_tan_alpha = cos_alpha / sin_alpha;
const TinyMatrix<Dimension> W{1, cos_alpha, //
0, sin_alpha};
const TinyMatrix<Dimension> inv_W = inverse(W);
- std::array<TinyMatrix<Dimension>, Dimension> //
- S_gradient = {TinyMatrix<Dimension>{-1, -1, //
+ std::array<TinyMatrix<Dimension>, Dimension> //
+ S_gradient = {TinyMatrix<Dimension>{-1, -inv_sin_alpha + inv_tan_alpha, //
+0, +0},
- TinyMatrix<Dimension>{+inv_tan_alpha, +inv_tan_alpha, //
- -inv_sin_alpha, -inv_sin_alpha}};
-
- f = 0;
- f_gradient = zero;
- f_hessian = zero;
+ TinyMatrix<Dimension>{+0, +0, //
+ -1, -inv_sin_alpha + inv_tan_alpha}};
SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
@@ -658,6 +654,10 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
auto escobar = escobarSSigma;
+ f = 0;
+ f_gradient = zero;
+ f_hessian = zero;
+
for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
const double sigma = sigma_list[i_cell];
const TinyMatrix<Dimension> S = S_list[i_cell];
@@ -674,12 +674,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]), //
trace(Sigma * S_gradient[1])};
-
- const std::array<TinyMatrix<Dimension>, Dimension> //
- Sigma_gradient{TinyMatrix<Dimension>{+0, +1, //
+ const std::array<TinyMatrix<Dimension>, Dimension> //
+ Sigma_gradient{TinyMatrix<Dimension>{+0, +inv_sin_alpha - inv_tan_alpha, //
+0, -1},
- TinyMatrix<Dimension>{-inv_sin_alpha, -inv_tan_alpha, //
- +inv_sin_alpha, +inv_tan_alpha}};
+ TinyMatrix<Dimension>{-inv_sin_alpha + inv_tan_alpha, +0, //
+ +1, +0}};
escobar(sigma, S, sigma_gradient, S_gradient, Sigma_gradient, f, f_gradient, f_hessian);
}
@@ -698,10 +697,10 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const double mr_alpha = 2 * std::acos(-1) / ((cell_nb_nodes - 2) * (cell_list.size()));
- const double mr_sin_alpha = std::sin(mr_alpha);
- const double mr_cos_alpha = std::cos(mr_alpha);
- const double mr_inv_sin_alpha = 1. / mr_sin_alpha;
- const double mr_inv_tan_alpha = 1. / std::tan(mr_alpha);
+ const double mr_sin_alpha = std::sin(mr_alpha);
+ const double mr_cos_alpha = std::cos(mr_alpha);
+ // const double mr_inv_sin_alpha = 1. / mr_sin_alpha;
+ // const double mr_inv_tan_alpha = 1. / std::tan(mr_alpha);
const TinyMatrix<Dimension> mr_W{1, mr_cos_alpha, //
0, mr_sin_alpha};
@@ -791,96 +790,164 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
auto cell_to_node_matrix = m_given_mesh.connectivity().cellToNodeMatrix();
auto node_number_in_their_cells = m_given_mesh.connectivity().nodeLocalNumbersInTheirCells();
+ constexpr bool q_sin = true;
+
for (NodeId node_id = 0; node_id < m_given_mesh.numberOfNodes(); ++node_id) {
if (not m_is_fixed_node[node_id]) {
auto node_cells = node_to_cell_matrix[node_id];
auto node_number = node_number_in_their_cells[node_id];
- double min_sin = std::numeric_limits<double>::max();
- double max_abs_sin = -std::numeric_limits<double>::max();
+ if constexpr (q_sin) { // q_area
+ double min_sin = std::numeric_limits<double>::max();
+ double max_abs_sin = -std::numeric_limits<double>::max();
- for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
- const CellId cell_id = node_cells[i_cell];
+ for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
+ const CellId cell_id = node_cells[i_cell];
- const auto cell_nodes = cell_to_node_matrix[cell_id];
+ const auto cell_nodes = cell_to_node_matrix[cell_id];
- const size_t i_A = node_number[i_cell];
- const size_t i_B = (i_A + 1) % cell_nodes.size();
- const size_t i_C = (i_A + cell_nodes.size() - 1) % cell_nodes.size();
- const size_t i_D = (i_A + 2) % cell_nodes.size();
- const size_t i_E = (i_A + cell_nodes.size() - 2) % cell_nodes.size();
+ const size_t i_A = node_number[i_cell];
+ const size_t i_B = (i_A + 1) % cell_nodes.size();
+ const size_t i_C = (i_A + cell_nodes.size() - 1) % cell_nodes.size();
+ const size_t i_D = (i_A + 2) % cell_nodes.size();
+ const size_t i_E = (i_A + cell_nodes.size() - 2) % cell_nodes.size();
- const Rd A = xr[cell_nodes[i_A]];
- const Rd B = xr[cell_nodes[i_B]];
- const Rd C = xr[cell_nodes[i_C]];
- const Rd D = xr[cell_nodes[i_D]];
- const Rd E = xr[cell_nodes[i_E]];
+ const Rd A = xr[cell_nodes[i_A]];
+ const Rd B = xr[cell_nodes[i_B]];
+ const Rd C = xr[cell_nodes[i_C]];
+ const Rd D = xr[cell_nodes[i_D]];
+ const Rd E = xr[cell_nodes[i_E]];
- const Rd AB = B - A;
- const Rd AC = C - A;
+ const Rd AB = B - A;
+ const Rd AC = C - A;
- const Rd CA = A - C;
- const Rd CE = E - C;
+ const Rd CA = A - C;
+ const Rd CE = E - C;
- const Rd BD = D - B;
- const Rd BA = A - B;
+ const Rd BD = D - B;
+ const Rd BA = A - B;
- const double l_AB = l2Norm(AB);
- const double l_AC = l2Norm(AC);
+ const double l_AB = l2Norm(AB);
+ const double l_AC = l2Norm(AC);
- const double l_CA = l2Norm(CA);
- const double l_CE = l2Norm(CE);
+ const double l_CA = l2Norm(CA);
+ const double l_CE = l2Norm(CE);
- const double l_BD = l2Norm(BD);
- const double l_BA = l2Norm(BA);
+ const double l_BD = l2Norm(BD);
+ const double l_BA = l2Norm(BA);
- if ((l_AB == 0) or (l_AC == 0) or (l_CA == 0) or (l_CE == 0) or (l_BD == 0) or (l_BA == 0)) {
- min_sin = -1;
- max_abs_sin = 1;
- break;
- }
+ if ((l_AB == 0) or (l_AC == 0) or (l_CA == 0) or (l_CE == 0) or (l_BD == 0) or (l_BA == 0)) {
+ min_sin = -1;
+ max_abs_sin = 1;
+ break;
+ }
- const Rd u_AB = 1. / l_AB * AB;
- const Rd u_AC = 1. / l_AC * AC;
- const Rd u_CA = 1. / l_CA * CA;
- const Rd u_CE = 1. / l_CE * CE;
- const Rd u_BD = 1. / l_BD * BD;
- const Rd u_BA = 1. / l_BA * BA;
+ const Rd u_AB = 1. / l_AB * AB;
+ const Rd u_AC = 1. / l_AC * AC;
+ const Rd u_CA = 1. / l_CA * CA;
+ const Rd u_CE = 1. / l_CE * CE;
+ const Rd u_BD = 1. / l_BD * BD;
+ const Rd u_BA = 1. / l_BA * BA;
- const double sin_CAB = det(TinyMatrix{u_AB, u_AC});
- const double sin_ECA = det(TinyMatrix{u_CA, u_CE});
- const double sin_ABD = det(TinyMatrix{u_BD, u_BA});
+ const double sin_CAB = det(TinyMatrix{u_AB, u_AC});
+ const double sin_ECA = det(TinyMatrix{u_CA, u_CE});
+ const double sin_ABD = det(TinyMatrix{u_BD, u_BA});
- const double abs_sin_CAB = std::abs(sin_CAB);
- const double abs_sin_ECA = std::abs(sin_ECA);
- const double abs_sin_ABD = std::abs(sin_ABD);
+ const double abs_sin_CAB = std::abs(sin_CAB);
+ const double abs_sin_ECA = std::abs(sin_ECA);
+ const double abs_sin_ABD = std::abs(sin_ABD);
- if (sin_CAB < min_sin) {
- min_sin = sin_CAB;
- }
+ if (sin_CAB < min_sin) {
+ min_sin = sin_CAB;
+ }
- if (sin_ECA < min_sin) {
- min_sin = sin_ECA;
- }
+ if (sin_ECA < min_sin) {
+ min_sin = sin_ECA;
+ }
- if (sin_ABD < min_sin) {
- min_sin = sin_ABD;
- }
+ if (sin_ABD < min_sin) {
+ min_sin = sin_ABD;
+ }
- if (abs_sin_CAB > max_abs_sin) {
- max_abs_sin = abs_sin_CAB;
- }
+ if (abs_sin_CAB > max_abs_sin) {
+ max_abs_sin = abs_sin_CAB;
+ }
- if (abs_sin_ECA > max_abs_sin) {
- max_abs_sin = abs_sin_ECA;
+ if (abs_sin_ECA > max_abs_sin) {
+ max_abs_sin = abs_sin_ECA;
+ }
+
+ if (abs_sin_ABD > max_abs_sin) {
+ max_abs_sin = abs_sin_ABD;
+ }
}
- if (abs_sin_ABD > max_abs_sin) {
- max_abs_sin = abs_sin_ABD;
+ node_quality[node_id] = min_sin / max_abs_sin;
+ } else {
+ double min_area = std::numeric_limits<double>::max();
+ double max_abs_area = -std::numeric_limits<double>::max();
+
+ for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
+ const CellId cell_id = node_cells[i_cell];
+
+ const auto cell_nodes = cell_to_node_matrix[cell_id];
+
+ const size_t i_A = node_number[i_cell];
+ const size_t i_B = (i_A + 1) % cell_nodes.size();
+ const size_t i_C = (i_A + cell_nodes.size() - 1) % cell_nodes.size();
+ const size_t i_D = (i_A + 2) % cell_nodes.size();
+ const size_t i_E = (i_A + cell_nodes.size() - 2) % cell_nodes.size();
+
+ const Rd A = xr[cell_nodes[i_A]];
+ const Rd B = xr[cell_nodes[i_B]];
+ const Rd C = xr[cell_nodes[i_C]];
+ const Rd D = xr[cell_nodes[i_D]];
+ const Rd E = xr[cell_nodes[i_E]];
+
+ const Rd AB = B - A;
+ const Rd AC = C - A;
+
+ const Rd CA = A - C;
+ const Rd CE = E - C;
+
+ const Rd BD = D - B;
+ const Rd BA = A - B;
+
+ const double area_CAB = det(TinyMatrix{AB, AC});
+ const double area_ECA = det(TinyMatrix{CA, CE});
+ const double area_ABD = det(TinyMatrix{BD, BA});
+
+ const double abs_area_CAB = std::abs(area_CAB);
+ const double abs_area_ECA = std::abs(area_ECA);
+ const double abs_area_ABD = std::abs(area_ABD);
+
+ if (area_CAB < min_area) {
+ min_area = area_CAB;
+ }
+
+ if (area_ECA < min_area) {
+ min_area = area_ECA;
+ }
+
+ if (area_ABD < min_area) {
+ min_area = area_ABD;
+ }
+
+ if (abs_area_CAB > max_abs_area) {
+ max_abs_area = abs_area_CAB;
+ }
+
+ if (abs_area_ECA > max_abs_area) {
+ max_abs_area = abs_area_ECA;
+ }
+
+ if (abs_area_ABD > max_abs_area) {
+ max_abs_area = abs_area_ABD;
+ }
}
- }
- node_quality[node_id] = min_sin / max_abs_sin;
+ node_quality[node_id] = min_area / max_abs_area;
+ }
} else {
node_quality[node_id] = 1;
}
@@ -920,11 +987,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
quality.fill(2);
- size_t maxiter = 0;
-
- auto smooth = [this, &maxiter, &old_xr](const NodeId node_id, TinyVector<Dimension>& x) -> double {
- // SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
-
+ auto smooth = [this, &old_xr](const NodeId node_id, TinyVector<Dimension>& x) -> double {
double delta = 0;
double final_f = 0;
@@ -933,7 +996,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const ConnectivityType& connectivity = m_given_mesh.connectivity();
- for (size_t i_iter = 0; i_iter < 200; ++i_iter) {
+ for (size_t i_iter = 0; i_iter < 20; ++i_iter) {
double sigma_min = std::numeric_limits<double>::max();
delta = sigma_min;
@@ -962,17 +1025,18 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
return max_lenght;
}();
- std::cout.precision(15);
-
if constexpr (use_newton) {
const TinyVector delta_x_candidate = -inverse(f_hessian) * f_gradient;
const double delta_x_norm = l2Norm(delta_x_candidate);
+ TinyVector<Dimension> old_x = x;
+
double factor = 1;
- if (delta_x_norm > 0.2 * max_edge_lenght) {
- factor *= (0.2 * max_edge_lenght / delta_x_norm);
+ if (delta_x_norm > 0.5 * max_edge_lenght) {
+ factor *= (0.5 * max_edge_lenght / delta_x_norm);
}
+
x += factor * delta_x_candidate;
if (delta_x_norm < 1E-8 * max_edge_lenght) {
@@ -1146,25 +1210,20 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
return final_f;
};
- // parallel_for(
- // connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- // if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
- // quality[node_id] = smooth(node_id, new_xr[node_id]);
- // }
- // });
-
- size_t count = 0;
- for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
- if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
- quality[node_id] = smooth(node_id, new_xr[node_id]);
- new_xr[node_id] = 0.5 * (new_xr[node_id] + old_xr[node_id]);
- // TinyVector<Dimension> x0{2, 3};
- // quality[node_id] = smooth(node_id, x0);
- count++;
+ parallel_for(
+ connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
+ quality[node_id] = smooth(node_id, new_xr[node_id]);
+ new_xr[node_id] = 0.5 * (new_xr[node_id] + old_xr[node_id]);
+ }
+ });
+ std::cout << "nb smoothed = " << [&] {
+ size_t count = 0;
+ for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+ count += (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]);
}
- }
- std::cout << "nb smoothed = " << count << '\n';
- std::cout << "maxiter = " << maxiter << '\n';
+ return count;
+ }() << '\n';
} else {
throw NotImplementedError("Dimension != 2");
}
@@ -1176,43 +1235,48 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
{
NodeValue<const Rd> xr = m_given_mesh.xr();
- size_t nb_non_convex_cells = 0;
- size_t i_convexify = 1;
+ size_t i_convexify = 1;
+ bool need_convexify = false;
+
+ std::cout << "== calling smoothMesh ==\n";
+
+ NodeValue<const double> old_node_quality;
do {
std::cout << "- convexfication step " << i_convexify++ << '\n';
NodeValue<Rd> new_xr = copy(xr);
auto node_quality = this->_getNodeQuality(xr);
- CellValue<const int> non_convex_cell_tag = [=, this] {
- CellValue<int> cell_tag(m_given_mesh.connectivity());
- cell_tag.fill(0);
+ need_convexify = min(node_quality) <= 0;
- auto cell_to_node_matrix = m_given_mesh.connectivity().cellToNodeMatrix();
+ if (need_convexify) {
+ m_is_smoothed_node = _getSmoothedNode(node_quality, m_is_fixed_node);
- for (CellId cell_id = 0; cell_id < m_given_mesh.numberOfCells(); ++cell_id) {
- for (size_t i_node = 0; i_node < cell_to_node_matrix[cell_id].size(); ++i_node) {
- const NodeId node_id = cell_to_node_matrix[cell_id][i_node];
- if (node_quality[node_id] <= m_min_quality) {
- cell_tag[cell_id] = 1;
- break;
+ if (old_node_quality.isBuilt()) {
+ NodeValue<bool> is_smoothed_first{m_given_mesh.connectivity()};
+ is_smoothed_first.fill(false);
+ bool has_new_smoothed_node = false;
+ for (NodeId node_id = 0; node_id < m_given_mesh.numberOfNodes(); ++node_id) {
+ if ((m_is_smoothed_node[node_id]) and (old_node_quality[node_id] >= 0)) {
+ is_smoothed_first[node_id] = true;
+ has_new_smoothed_node = true;
}
}
- }
-
- return cell_tag;
- }();
- // auto non_convex_cell_tag = this->_getNonConvexCellTag(xr);
- nb_non_convex_cells = sum(non_convex_cell_tag);
+ if (has_new_smoothed_node) {
+ std::cout << "has_new_smoothed_node = " << has_new_smoothed_node << '\n';
+ m_is_smoothed_node = is_smoothed_first;
+ }
+ }
- std::cout << " remaining non convex cells: " << nb_non_convex_cells << '\n';
- if (nb_non_convex_cells > 0) {
- m_is_smoothed_node = _getSmoothedNode(node_quality, m_is_fixed_node);
this->_getNewPositions(xr, new_xr);
xr = new_xr;
}
- } while ((nb_non_convex_cells > 0) and (i_convexify < 20));
+
+ old_node_quality = node_quality;
+ std::cout << "need_convexify = " << std::boolalpha << need_convexify << '\n';
+ } while ((need_convexify) and (i_convexify < 50));
+ std::cout << " --> finished\n";
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}
diff --git a/src/mesh/MeshSmootherJun.cpp b/src/mesh/MeshSmootherJun.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..20821c5713b95b22f132ea3416170d01dd74b718
--- /dev/null
+++ b/src/mesh/MeshSmootherJun.cpp
@@ -0,0 +1,627 @@
+#include <mesh/MeshSmootherJun.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshCellZone.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshLineNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <scheme/AxisBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/RandomEngine.hpp>
+
+#include <variant>
+
+template <size_t Dimension>
+class MeshSmootherJunHandler::MeshSmootherJun
+{
+ private:
+ using Rd = TinyVector<Dimension>;
+ using Rdxd = TinyMatrix<Dimension>;
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ const MeshType& m_given_mesh;
+
+ class AxisBoundaryCondition;
+ class FixedBoundaryCondition;
+ class SymmetryBoundaryCondition;
+
+ using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+ BoundaryConditionList m_boundary_condition_list;
+
+ BoundaryConditionList
+ _getBCList(const MeshType& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ {
+ BoundaryConditionList bc_list;
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::axis: {
+ if constexpr (Dimension == 1) {
+ bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+ } else {
+ bc_list.emplace_back(
+ AxisBoundaryCondition{getMeshLineNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ bc_list.emplace_back(
+ SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::fixed: {
+ bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for mesh smoother";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
+
+ return bc_list;
+ }
+
+ void
+ _applyBC(NodeValue<Rd>& shift) const
+ {
+ for (auto&& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using BCType = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+ const Rd& n = bc.outgoingNormal();
+
+ const Rdxd I = identity;
+ const Rdxd nxn = tensorProduct(n, n);
+ const Rdxd P = I - nxn;
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ shift[node_id] = P * shift[node_id];
+ });
+
+ } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
+ if constexpr (Dimension > 1) {
+ const Rd& t = bc.direction();
+
+ const Rdxd txt = tensorProduct(t, t);
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ shift[node_id] = txt * shift[node_id];
+ });
+ } else {
+ throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
+ }
+
+ } else if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ shift[node_id] = zero;
+ });
+
+ } else {
+ throw UnexpectedError("invalid boundary condition type");
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ NodeValue<Rd>
+ _getDisplacement() const
+ {
+ const ConnectivityType& connectivity = m_given_mesh.connectivity();
+ NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+
+ NodeValue<double> max_delta_xr{connectivity};
+ parallel_for(
+ connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const Rd& x0 = given_xr[node_id];
+
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ double min_distance_2 = std::numeric_limits<double>::max();
+
+ for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+ const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
+
+ const CellId cell_id = node_cell_list[i_cell];
+ const auto& cell_node_list = cell_to_node_matrix[cell_id];
+
+ for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
+ if (i_node != i_cell_node) {
+ const NodeId cell_node_id = cell_node_list[i_node];
+ const Rd delta = x0 - given_xr[cell_node_id];
+ min_distance_2 = std::min(min_distance_2, dot(delta, delta));
+ }
+ }
+ }
+ double max_delta = std::sqrt(min_distance_2);
+
+ max_delta_xr[node_id] = max_delta;
+ });
+
+ NodeValue<Rd> shift_r{connectivity};
+
+ if constexpr (Dimension == 2) {
+ parallel_for(
+ m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const auto node_cell_list = node_to_cell_matrix[node_id];
+
+ constexpr TinyMatrix<Dimension> R{0, 1, -1, 0};
+
+ double a_r = 0;
+ Rd b_r = zero;
+
+ for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+ const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
+
+ const CellId cell_id = node_cell_list[i_cell];
+
+ const auto cell_node_list = cell_to_node_matrix[cell_id];
+ const size_t cell_nb_node = cell_node_list.size();
+
+ const size_t i_next_node = (i_cell_node + 1) % cell_nb_node;
+ const size_t i_prev_node = (i_cell_node + cell_nb_node - 1) % cell_nb_node;
+
+ const NodeId next_node_id = cell_node_list[i_next_node];
+ const NodeId prev_node_id = cell_node_list[i_prev_node];
+
+ const Rd& xrm1 = given_xr[prev_node_id];
+ const Rd& xr = given_xr[node_id];
+ const Rd& xrp1 = given_xr[next_node_id];
+
+ const Rd& a = given_xr[next_node_id];
+ const Rd& b = given_xr[prev_node_id];
+
+ const Rd a_xr = a - xr;
+ const Rd b_xr = b - xr;
+
+ const double a11 = dot(b_xr, b_xr);
+ const double a22 = dot(a_xr, a_xr);
+ const double a12 = -dot(a_xr, b_xr); //(a_xr[0] + a_xr[1]) * (b_xr[0] + b_xr[1]);
+
+ TinyMatrix<Dimension + 1> M;
+
+ M(0, 0) = 0.5 * (a11 + 2 * a12 + a22);
+ M(0, 1) = 0.5 * (-a11 - a12);
+ M(0, 2) = 0.5 * (-a12 - a22);
+
+ M(1, 0) = 0.5 * (-a11 - a12);
+ M(1, 1) = 0.5 * (a11);
+ M(1, 2) = 0.5 * (a12);
+
+ M(2, 0) = 0.5 * (-a12 - a22);
+ M(2, 1) = 0.5 * (a12);
+ M(2, 2) = 0.5 * (a22);
+
+ // const Rd xrm1{1, 3};
+ // const Rd xr(1, 1);
+ // const Rd xrp1(4, 1);
+
+ // const Rd xr_xrp1_perp = R * (xrp1 - xr);
+ // const Rd xrp1_xrm1_perp = R * (xrm1 - xrp1);
+ // const Rd xrm1_xr_perp = R * (xr - xrm1);
+
+ // std::cout << "xr_xrp1_perp = " << xr_xrp1_perp << '\n';
+ // std::cout << "xrp1_xrm1_perp = " << xrp1_xrm1_perp << '\n';
+ // std::cout << "xrm1_xr_perp = " << xrm1_xr_perp << '\n';
+ // const Rd Cjr = 0.5 * (xr_xrp1_perp + xrm1_xr_perp);
+ // const Rd Cjrp1 = 0.5 * (xr_xrp1_perp + xrp1_xrm1_perp);
+ // const Rd Cjrm1 = 0.5 * (xrm1_xr_perp + xrp1_xrm1_perp);
+
+ // std::cout << "Cjr = " << Cjr << '\n';
+ // std::cout << "Cjrp1 = " << Cjrp1 << '\n';
+ // std::cout << "Cjrm1 = " << Cjrm1 << '\n';
+
+ // const double Sj = 0.5 * (dot(Cjr, xr) + dot(Cjrm1, xrm1) + dot(Cjrp1, xrp1));
+
+ // std::cout << "grad Wjr = " << 1 / Sj * Cjr << '\n';
+ // std::cout << "grad Wjrp1 = " << 1 / Sj * Cjrp1 << '\n';
+ // std::cout << "grad Wjrm1 = " << 1 / Sj * Cjrm1 << '\n';
+ // std::cout << "Sj = " << Sj << '\n';
+
+ // b_r[0] += M(i_cell_node, i_prev_node) * xrm1[0] + M(i_cell_node, i_next_node) * xrp1[0];
+ // b_r[1] += M(i_cell_node, i_prev_node) * xrm1[1] + M(i_cell_node, i_next_node) * xrp1[1];
+ // a_r += M(i_cell_node, i_cell_node);
+
+ b_r[0] += M(0, 2) * xrm1[0] + M(0, 1) * xrp1[0];
+ b_r[1] += M(0, 2) * xrm1[1] + M(0, 1) * xrp1[1];
+ a_r += M(0, 0);
+
+ // std::cout << " ===\n";
+ }
+
+ // std::cout << "ar = " << a_r << '\n';
+ // std::cout << "br = " << b_r << '\n';
+ // std::cout << "a_r = " << a_r << '\n';
+ // std::cout << "b_r = " << b_r << '\n';
+ if (a_r != 0) {
+ Rd new_x = -1. / a_r * b_r;
+ // std::cout << rang::fgB::yellow << "new_x = " << new_x << rang::fg::reset << '\n';
+ shift_r[node_id] = new_x - given_xr[node_id]; //-1. / a_r * b_r;
+ } else {
+ shift_r[node_id] = zero;
+ }
+ // std::cout << " -1. / a_r * b_r = " << -1. / a_r * b_r << '\n';
+ // std::cout << " given_xr[" << node_id << "] = " << given_xr[node_id] << '\n';
+ // std::cout << "shift[" << node_id << "] = " << shift_r[node_id] << '\n';
+ // std::cout << "-------------------\n";
+ });
+ } else {
+ throw NotImplementedError("Jun smoother in dimension != 2");
+ }
+
+ this->_applyBC(shift_r);
+
+ synchronize(shift_r);
+
+ return shift_r;
+ }
+
+ public:
+ std::shared_ptr<const IMesh>
+ getSmoothedMesh() const
+ {
+ NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+ NodeValue<Rd> xr = this->_getDisplacement();
+
+ parallel_for(
+ m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { xr[node_id] += given_xr[node_id]; });
+
+ return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+ }
+
+ std::shared_ptr<const IMesh>
+ getSmoothedMesh(const FunctionSymbolId& function_symbol_id) const
+ {
+ NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+ NodeValue<const bool> is_displaced =
+ InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr);
+
+ NodeValue<Rd> xr = this->_getDisplacement();
+
+ parallel_for(
+ m_given_mesh.numberOfNodes(),
+ PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+ return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+ }
+
+ std::shared_ptr<const IMesh>
+ getSmoothedMesh(const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_descriptor_list) const
+ {
+ NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+ auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
+
+ NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+ is_displaced.fill(false);
+
+ for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) {
+ MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
+ const auto cell_list = cell_zone.cellList();
+ CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
+ is_zone_cell.fill(false);
+ parallel_for(
+ cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; });
+ parallel_for(
+ m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ auto node_cell_list = node_to_cell_matrix[node_id];
+ bool displace = true;
+ for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+ const CellId cell_id = node_cell_list[i_node_cell];
+ displace &= is_zone_cell[cell_id];
+ }
+ if (displace) {
+ is_displaced[node_id] = true;
+ }
+ });
+ }
+ synchronize(is_displaced);
+ NodeValue<Rd> xr = this->_getDisplacement();
+
+ parallel_for(
+ m_given_mesh.numberOfNodes(),
+ PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+ return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+ }
+
+ std::shared_ptr<const IMesh>
+ getSmoothedMesh(
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+ {
+ NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+ auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
+
+ NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+ is_displaced.fill(false);
+
+ for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
+ auto is_zone_cell = discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const double>>();
+
+ parallel_for(
+ m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ auto node_cell_list = node_to_cell_matrix[node_id];
+ bool displace = true;
+ for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+ const CellId cell_id = node_cell_list[i_node_cell];
+ displace &= (is_zone_cell[cell_id] != 0);
+ }
+ if (displace) {
+ is_displaced[node_id] = true;
+ }
+ });
+ }
+ synchronize(is_displaced);
+ NodeValue<Rd> xr = this->_getDisplacement();
+
+ parallel_for(
+ m_given_mesh.numberOfNodes(),
+ PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+ return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+ }
+
+ MeshSmootherJun(const MeshSmootherJun&) = delete;
+ MeshSmootherJun(MeshSmootherJun&&) = delete;
+
+ MeshSmootherJun(const MeshType& given_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ : m_given_mesh(given_mesh), m_boundary_condition_list(this->_getBCList(given_mesh, bc_descriptor_list))
+ {}
+
+ ~MeshSmootherJun() = default;
+};
+
+template <size_t Dimension>
+class MeshSmootherJunHandler::MeshSmootherJun<Dimension>::AxisBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
+
+ public:
+ const Rd&
+ direction() const
+ {
+ return m_mesh_line_node_boundary.direction();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_line_node_boundary.nodeList();
+ }
+
+ AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
+ : m_mesh_line_node_boundary(mesh_line_node_boundary)
+ {
+ ;
+ }
+
+ ~AxisBoundaryCondition() = default;
+};
+
+template <>
+class MeshSmootherJunHandler::MeshSmootherJun<1>::AxisBoundaryCondition
+{
+ public:
+ AxisBoundaryCondition() = default;
+ ~AxisBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class MeshSmootherJunHandler::MeshSmootherJun<Dimension>::FixedBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ FixedBoundaryCondition(MeshNodeBoundary<Dimension>&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+
+ ~FixedBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class MeshSmootherJunHandler::MeshSmootherJun<Dimension>::SymmetryBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+
+ public:
+ const Rd&
+ outgoingNormal() const
+ {
+ return m_mesh_flat_node_boundary.outgoingNormal();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+std::shared_ptr<const IMesh>
+MeshSmootherJunHandler::getSmoothedMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+{
+ switch (mesh->dimension()) {
+ case 1: {
+ constexpr size_t Dimension = 1;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh();
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh();
+ }
+ case 3: {
+ constexpr size_t Dimension = 3;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh();
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherJunHandler::getSmoothedMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& function_symbol_id) const
+{
+ switch (mesh->dimension()) {
+ case 1: {
+ constexpr size_t Dimension = 1;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(function_symbol_id);
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(function_symbol_id);
+ }
+ case 3: {
+ constexpr size_t Dimension = 3;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(function_symbol_id);
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherJunHandler::getSmoothedMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const
+{
+ switch (mesh->dimension()) {
+ case 1: {
+ constexpr size_t Dimension = 1;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(smoothing_zone_list);
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(smoothing_zone_list);
+ }
+ case 3: {
+ constexpr size_t Dimension = 3;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(smoothing_zone_list);
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherJunHandler::getSmoothedMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
+{
+ if (not hasSameMesh(discrete_function_variant_list)) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ std::shared_ptr<const IMesh> common_mesh = getCommonMesh(discrete_function_variant_list);
+
+ if (common_mesh != mesh) {
+ throw NormalError("discrete functions are not defined on the smoothed mesh");
+ }
+
+ switch (mesh->dimension()) {
+ case 1: {
+ constexpr size_t Dimension = 1;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(discrete_function_variant_list);
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(discrete_function_variant_list);
+ }
+ case 3: {
+ constexpr size_t Dimension = 3;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherJun smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh(discrete_function_variant_list);
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
diff --git a/src/mesh/MeshSmootherJun.hpp b/src/mesh/MeshSmootherJun.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a9a34438fdeec2fee1108d41c319d7ca5c33acc3
--- /dev/null
+++ b/src/mesh/MeshSmootherJun.hpp
@@ -0,0 +1,45 @@
+#ifndef MESH_SMOOTHER_JUN_HPP
+#define MESH_SMOOTHER_JUN_HPP
+
+#include <mesh/IMesh.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+class FunctionSymbolId;
+class IZoneDescriptor;
+class DiscreteFunctionVariant;
+
+class MeshSmootherJunHandler
+{
+ private:
+ template <size_t Dimension>
+ class MeshSmootherJun;
+
+ public:
+ std::shared_ptr<const IMesh> getSmoothedMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const;
+
+ std::shared_ptr<const IMesh> getSmoothedMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& function_symbol_id) const;
+
+ std::shared_ptr<const IMesh> getSmoothedMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const;
+
+ std::shared_ptr<const IMesh> getSmoothedMesh(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& smoothing_zone_list) const;
+
+ MeshSmootherJunHandler() = default;
+ MeshSmootherJunHandler(MeshSmootherJunHandler&&) = default;
+ ~MeshSmootherJunHandler() = default;
+};
+
+#endif // MESH_SMOOTHER_JUN_HPP