diff --git a/src/geometry/TriangleTransformation.hpp b/src/geometry/TriangleTransformation.hpp
index 9210e0bda55bb19dfdfd877e509801fa87f0a4f0..240d0c919068ceb3aebb715f89a9e360ffe5f0a8 100644
--- a/src/geometry/TriangleTransformation.hpp
+++ b/src/geometry/TriangleTransformation.hpp
@@ -19,6 +19,13 @@ class TriangleTransformation<2>
double m_jacobian_determinant;
public:
+ PUGS_INLINE
+ TinyMatrix<Dimension>
+ jacobian() const
+ {
+ return m_jacobian;
+ }
+
PUGS_INLINE
TinyVector<Dimension>
operator()(const TinyVector<2>& x) const
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 813bde93a055a514e35a1e949ca6ea2241609ad2..fb6e9600205cd4e14fc73b09566b6b2d7e4cb407 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -17,6 +17,7 @@
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshRandomizer.hpp>
#include <mesh/MeshSmoother.hpp>
+#include <mesh/MeshSmootherEscobar.hpp>
#include <scheme/AcousticSolver.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
@@ -306,6 +307,30 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("smoothMeshEscobar",
+ 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> {
+ MeshSmootherEscobarHandler handler;
+ return handler.getSmoothedMesh(p_mesh, bc_descriptor_list);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("qualityEscobar",
+ std::function(
+#warning REMOVE BEFORE MERGING
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::shared_ptr<const ItemValueVariant> {
+ MeshSmootherEscobarHandler handler;
+ return handler.getQuality(p_mesh, bc_descriptor_list);
+ }
+
+ ));
+
this->_addBuiltinFunction("fixed", std::function(
[](std::shared_ptr<const IBoundaryDescriptor> boundary)
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index ee69a9710d8706072e1fbd80eb51a6d3e0dd9b5c..54c7afbf2f19799524faf0038e944d5463e1f7c1 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -39,5 +39,6 @@ add_library(
MeshNodeInterface.cpp
MeshRandomizer.cpp
MeshSmoother.cpp
+ MeshSmootherEscobar.cpp
MeshTransformer.cpp
SynchronizerManager.cpp)
diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..793153c0abd0a169477182e34ef40072bffa7bec
--- /dev/null
+++ b/src/mesh/MeshSmootherEscobar.cpp
@@ -0,0 +1,654 @@
+#include <mesh/MeshSmootherEscobar.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <geometry/TriangleTransformation.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/ItemValueVariant.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 MeshSmootherEscobarHandler::MeshSmootherEscobar
+{
+ 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};
+
+ parallel_for(
+ m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ Rd mean_position(zero);
+ size_t number_of_neighbours = 0;
+
+ 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];
+ mean_position += given_xr[cell_node_id];
+ number_of_neighbours++;
+ }
+ }
+ }
+ mean_position = 1. / number_of_neighbours * mean_position;
+ shift_r[node_id] = mean_position - given_xr[node_id];
+ });
+
+ this->_applyBC(shift_r);
+
+ synchronize(shift_r);
+
+ return shift_r;
+ }
+
+ public:
+ std::shared_ptr<const ItemValueVariant>
+ getQuality() const
+ {
+ if constexpr (Dimension == 2) {
+ const ConnectivityType& connectivity = m_given_mesh.connectivity();
+ NodeValue<const Rd> xr = m_given_mesh.xr();
+
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+ auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+
+ auto is_boundary_node = connectivity.isBoundaryNode();
+ NodeValue<double> quality{connectivity};
+
+ constexpr double eps = 1E-15;
+ quality.fill(2);
+
+ auto f_inner = [=](const NodeId node_id, const TinyVector<Dimension>& x) -> double {
+ auto cell_list = node_to_cell_matrix[node_id];
+ auto node_number_in_cell = node_number_in_their_cells[node_id];
+
+ const double alpha = 2 * std::acos(-1) / cell_list.size();
+ const TinyMatrix<Dimension> W{1, std::cos(alpha), 0, std::sin(alpha)};
+
+ const TinyMatrix<Dimension> inv_W = inverse(W);
+
+ SmallArray<TinyMatrix<Dimension>> S(cell_list.size());
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const size_t i_cell_node = node_number_in_cell[i_cell];
+ auto cell_node_list = cell_to_node_matrix[cell_list[i_cell]];
+ const size_t cell_nb_nodes = cell_node_list.size();
+
+ TriangleTransformation<Dimension> T(x, //
+ xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]],
+ xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]]);
+ S[i_cell] = T.jacobian() * inv_W;
+ }
+
+ SmallArray<double> sigma(S.size());
+ for (size_t j = 0; j < S.size(); ++j) {
+ sigma[j] = det(S[j]);
+ }
+
+ const double sigma_min = min(sigma);
+
+ const double delta =
+ (sigma_min < eps) ? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min)) : 0;
+
+ double f = 0;
+ for (size_t j = 0; j < S.size(); ++j) {
+ const TinyMatrix<Dimension> Sigma = sigma[j] * inverse(S[j]);
+
+ const double Sj_norm = std::sqrt(trace(transpose(S[j]) * S[j]));
+ const double Sigma_norm = std::sqrt(trace(transpose(Sigma) * Sigma));
+
+ f += Sj_norm * Sigma_norm / (sigma[j] + std::sqrt(sigma[j] * sigma[j] + 4 * delta * delta));
+ }
+ return f;
+ };
+
+ parallel_for(
+ connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ auto cell_list = node_to_cell_matrix[node_id];
+ auto node_number_in_cell = node_number_in_their_cells[node_id];
+
+ if (is_boundary_node[node_id]) {
+ quality[node_id] = 1;
+ } else {
+ TinyVector x = xr[node_id];
+ quality[node_id] = f_inner(node_id, x);
+
+ TinyMatrix<Dimension> B = identity;
+ }
+ });
+
+ return std::make_shared<ItemValueVariant>(quality);
+ } else {
+ throw NotImplementedError("Dimension != 2");
+ }
+ }
+
+ 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);
+ }
+
+ MeshSmootherEscobar(const MeshSmootherEscobar&) = delete;
+ MeshSmootherEscobar(MeshSmootherEscobar&&) = delete;
+
+ MeshSmootherEscobar(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))
+ {}
+
+ ~MeshSmootherEscobar() = default;
+};
+
+template <size_t Dimension>
+class MeshSmootherEscobarHandler::MeshSmootherEscobar<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 MeshSmootherEscobarHandler::MeshSmootherEscobar<1>::AxisBoundaryCondition
+{
+ public:
+ AxisBoundaryCondition() = default;
+ ~AxisBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class MeshSmootherEscobarHandler::MeshSmootherEscobar<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 MeshSmootherEscobarHandler::MeshSmootherEscobar<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 ItemValueVariant>
+MeshSmootherEscobarHandler::getQuality(
+ 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>>;
+ MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getQuality();
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getQuality();
+ }
+ case 3: {
+ constexpr size_t Dimension = 3;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getQuality();
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherEscobarHandler::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>>;
+ MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh();
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh();
+ }
+ case 3: {
+ constexpr size_t Dimension = 3;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), bc_descriptor_list);
+ return smoother.getSmoothedMesh();
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
+
+std::shared_ptr<const IMesh>
+MeshSmootherEscobarHandler::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>>;
+ MeshSmootherEscobar 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>>;
+ MeshSmootherEscobar 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>>;
+ MeshSmootherEscobar 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>
+MeshSmootherEscobarHandler::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>>;
+ MeshSmootherEscobar 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>>;
+ MeshSmootherEscobar 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>>;
+ MeshSmootherEscobar 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>
+MeshSmootherEscobarHandler::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>>;
+ MeshSmootherEscobar 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>>;
+ MeshSmootherEscobar 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>>;
+ MeshSmootherEscobar 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/MeshSmootherEscobar.hpp b/src/mesh/MeshSmootherEscobar.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f60ac7de08b762a1220bfc2ba2fbb4e21b0ac6e6
--- /dev/null
+++ b/src/mesh/MeshSmootherEscobar.hpp
@@ -0,0 +1,50 @@
+#ifndef MESH_SMOOTHER_ESCOBAR_HPP
+#define MESH_SMOOTHER_ESCOBAR_HPP
+
+#include <mesh/IMesh.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+class FunctionSymbolId;
+class IZoneDescriptor;
+class DiscreteFunctionVariant;
+class ItemValueVariant;
+
+class MeshSmootherEscobarHandler
+{
+ private:
+ template <size_t Dimension>
+ class MeshSmootherEscobar;
+
+ public:
+ std::shared_ptr<const ItemValueVariant> getQuality(
+ 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;
+
+ 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;
+
+ MeshSmootherEscobarHandler() = default;
+ MeshSmootherEscobarHandler(MeshSmootherEscobarHandler&&) = default;
+ ~MeshSmootherEscobarHandler() = default;
+};
+
+#endif // MESH_SMOOTHER_ESCOBAR_HPP