diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 31dd68cfdd15a88cdb75b091a15a717c2d9beecd..bf6da680dee9204324a8c5d2b4b10c2dde115b2e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -268,6 +268,34 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("implicitSmoothMesh",
+ 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> {
+ ImplicitMeshSmootherHandler handler;
+ return handler.getSmoothedMesh(p_mesh, bc_descriptor_list, function_symbol_id);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("implicitSmoothMesh",
+ 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> {
+ ImplicitMeshSmootherHandler handler;
+ return handler.getSmoothedMesh(p_mesh, bc_descriptor_list,
+ discrete_function_variant_list);
+ }
+
+ ));
+
this->_addBuiltinFunction("smoothMesh", std::function(
[](std::shared_ptr<const IMesh> p_mesh,
diff --git a/src/mesh/ImplicitMeshSmoother.cpp b/src/mesh/ImplicitMeshSmoother.cpp
index 64c5b67101dbee80d89513d7f97864166948e76e..43538914a2069268254c8d76d765bdf0ce39d7bc 100644
--- a/src/mesh/ImplicitMeshSmoother.cpp
+++ b/src/mesh/ImplicitMeshSmoother.cpp
@@ -32,10 +32,10 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
// class AxisBoundaryCondition;
class FixedBoundaryCondition;
- // class SymmetryBoundaryCondition;
+ class SymmetryBoundaryCondition;
- // using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
- using BoundaryCondition = std::variant<FixedBoundaryCondition>;
+ // using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
+ using BoundaryCondition = std::variant<FixedBoundaryCondition, SymmetryBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
BoundaryConditionList m_boundary_condition_list;
@@ -48,21 +48,21 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
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::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;
@@ -240,7 +240,7 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
}
NodeValue<Rd>
- _getPosition() const
+ _getPosition(NodeValue<const bool> is_displaced) const
{
const ConnectivityType& connectivity = m_given_mesh.connectivity();
NodeValue<const Rd> given_xr = m_given_mesh.xr();
@@ -249,7 +249,12 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
NodeValue<bool> is_fixed{connectivity};
is_fixed.fill(false);
parallel_for(
- m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { pos_r[node_id] = given_xr[node_id]; });
+ m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ pos_r[node_id] = given_xr[node_id];
+ if (not is_displaced[node_id]) {
+ is_fixed[node_id] = true;
+ }
+ });
_browseBC(is_fixed);
NodeValue<int> node_dof_id{connectivity};
@@ -257,51 +262,17 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
NodeValue<int> id_free{connectivity};
NodeValue<int> id_fixed{connectivity};
_computeMatrixSize(is_fixed, id_fixed, id_free, nb_free, nb_fixed);
- // std::cout << " is_fixed: "
- // << "\n";
- // std::cout << is_fixed << "\n";
- // std::cout << " id_fixed: "
- // << "\n";
- // std::cout << id_fixed << "\n";
- // std::cout << " is_free: "
- // << "\n";
- // std::cout << id_free << "\n";
- // std::cout << "nb_free " << nb_free << " nb_fixed " << nb_fixed << "\n";
Array<int> non_zeros_free{nb_free};
Array<int> non_zeros_fixed{nb_free};
Array<NodeId> gid_node_free{nb_free};
Array<NodeId> gid_node_fixed{nb_fixed};
- // size_t local_free_id = 0;
- // size_t local_fixed_id = 0;
_findLocalIds(is_fixed, gid_node_free, gid_node_fixed, non_zeros_free, non_zeros_fixed);
- // for (NodeId n_id = 0; n_id < m_given_mesh.numberOfNodes(); n_id++) {
- // if (is_fixed[n_id]) {
- // gid_node_fixed[local_fixed_id] = n_id;
- // local_fixed_id++;
- // } else {
- // gid_node_free[local_free_id] = n_id;
- // local_free_id++;
- // }
- // }
- // std::cout << " gid_node_fixed: "
- // << "\n";
- // std::cout << gid_node_fixed << "\n";
- // std::cout << " gid_node_free: "
- // << "\n";
- // std::cout << gid_node_free << "\n";
std::cout << " nb_free " << nb_free << " nb_fixed " << nb_fixed << "\n";
CRSMatrixDescriptor<double> Afree(nb_free, nb_free, non_zeros_free);
CRSMatrixDescriptor<double> Afixed(nb_free, nb_fixed, non_zeros_fixed);
LinearSolver solver;
_fillMatrix(Afree, Afixed, is_fixed, id_free, id_fixed);
- // std::cout << " Afree :"
- // << "\n";
- // std::cout << Afree << "\n";
- // std::cout << " Afixed :"
- // << "\n";
- // std::cout << Afixed << "\n";
-
Vector<double> F{nb_fixed};
CRSMatrix Mfree{Afree.getCRSMatrix()};
CRSMatrix Mfixed{Afixed.getCRSMatrix()};
@@ -312,8 +283,6 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
Vector<double> X{nb_free};
Vector<double> b{nb_free};
b = Mfixed * F;
- // std::cout << " F " << F << "\n";
- // std::cout << " b " << b << "\n";
solver.solveLocalSystem(Mfree, X, b);
parallel_for(
m_given_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
@@ -322,6 +291,7 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
}
});
}
+ // synchronize(pos_r);
return pos_r;
}
@@ -329,9 +299,54 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
std::shared_ptr<const IMesh>
getSmoothedMesh() const
{
+ NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+ is_displaced.fill(true);
NodeValue<const Rd> given_xr = m_given_mesh.xr();
- NodeValue<Rd> xr = this->_getPosition();
+ NodeValue<Rd> xr = this->_getPosition(is_displaced);
+
+ 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->_getPosition(is_displaced);
+
+ 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;
+ }
+ });
+ }
+
+ NodeValue<Rd> xr = this->_getPosition(is_displaced);
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}
@@ -347,45 +362,6 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother
~ImplicitMeshSmoother() = default;
};
-// template <size_t Dimension>
-// class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<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 ImplicitMeshSmootherHandler::ImplicitMeshSmoother<1>::AxisBoundaryCondition
-// {
-// public:
-// AxisBoundaryCondition() = default;
-// ~AxisBoundaryCondition() = default;
-// };
-
template <size_t Dimension>
class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<Dimension>::FixedBoundaryCondition
{
@@ -404,36 +380,36 @@ class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<Dimension>::FixedBoundar
~FixedBoundaryCondition() = default;
};
-// template <size_t Dimension>
-// class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<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;
-// };
+template <size_t Dimension>
+class ImplicitMeshSmootherHandler::ImplicitMeshSmoother<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>
ImplicitMeshSmootherHandler::getSmoothedMesh(
@@ -462,3 +438,70 @@ ImplicitMeshSmootherHandler::getSmoothedMesh(
}
}
}
+std::shared_ptr<const IMesh>
+ImplicitMeshSmootherHandler::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: {
+ throw NotImplementedError("ImplicitMeshSmoother not implemented in 1D");
+ break;
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ ImplicitMeshSmoother 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>>;
+ ImplicitMeshSmoother 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>
+ImplicitMeshSmootherHandler::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: {
+ throw NotImplementedError("ImplicitMeshSmoother not implemented in 1D");
+ break;
+ }
+ case 2: {
+ constexpr size_t Dimension = 2;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+ ImplicitMeshSmoother 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>>;
+ ImplicitMeshSmoother 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/ImplicitMeshSmoother.hpp b/src/mesh/ImplicitMeshSmoother.hpp
index b5ca016bc1c5c09cc9dcafb5a237788e3587edcd..7678417a5e6274f298b54783ff9a1354bcfdabdd 100644
--- a/src/mesh/ImplicitMeshSmoother.hpp
+++ b/src/mesh/ImplicitMeshSmoother.hpp
@@ -41,6 +41,14 @@ class ImplicitMeshSmootherHandler
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 DiscreteFunctionVariant>>& discrete_function_variant_list) const;
// std::shared_ptr<const IMesh> getSmoothedMesh(
// const std::shared_ptr<const IMesh>& mesh,