From 45736307c568e1c94ab2990e1b9407178628bae0 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Thu, 1 Jun 2023 12:34:10 +0200
Subject: [PATCH] Add contraction operator to TinyMatrix and an implicit mesh
smoother
---
src/mesh/ImplicitMeshSmoother.cpp | 463 ++++++++++++++++++++++++++++++
src/mesh/ImplicitMeshSmoother.hpp | 65 +++++
2 files changed, 528 insertions(+)
create mode 100644 src/mesh/ImplicitMeshSmoother.cpp
create mode 100644 src/mesh/ImplicitMeshSmoother.hpp
diff --git a/src/mesh/ImplicitMeshSmoother.cpp b/src/mesh/ImplicitMeshSmoother.cpp
new file mode 100644
index 000000000..9a13e8098
--- /dev/null
+++ b/src/mesh/ImplicitMeshSmoother.cpp
@@ -0,0 +1,463 @@
+#include <mesh/ImplicitMeshSmoother.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 ImplicitMeshSmootherHandler::ImplicitMeshSmoother
+{
+ 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 BoundaryCondition = std::variant<FixedBoundaryCondition>;
+
+ 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
+ _browseBC(NodeValue<bool>& is_fixed) 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;
+ is_fixed[node_id] = true;
+ });
+
+ } else {
+ throw UnexpectedError("invalid boundary condition type");
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _computeMatrixSize(const NodeValue<bool>& is_fixed,
+ NodeValue<int>& id_fixed,
+ NodeValue<int>& id_free,
+ size_t& nb_free,
+ size_t& nb_fixed) const
+ {
+ nb_free = 0;
+ nb_fixed = 0;
+ for (NodeId n_id = 0; n_id < m_given_mesh.numberOfNodes(); n_id++) {
+ if (is_fixed[n_id]) {
+ id_fixed[n_id] = nb_fixed;
+ nb_fixed += 1;
+ } else {
+ id_free[n_id] = nb_free;
+ nb_free += 1;
+ }
+ }
+ }
+ void
+ _findLocalIds(const NodeValue<bool>& is_fixed,
+ Array<NodeId>& gid_node_free,
+ Array<NodeId>& gid_node_fixed,
+ Array<int>& non_zeros_free,
+ Array<int>& non_zeros_fixed) const
+ {
+ const auto& node_to_face_matrix = m_given_mesh.connectivity().nodeToFaceMatrix();
+ const auto& face_to_node_matrix = m_given_mesh.connectivity().faceToNodeMatrix();
+
+ int local_free_id = 0;
+ int local_fixed_id = 0;
+ non_zeros_free.fill(1);
+ non_zeros_fixed.fill(0);
+ for (NodeId n_id = 0; n_id < m_given_mesh.numberOfNodes(); n_id++) {
+ std::cout << " n_id " << n_id << "\n";
+ if (is_fixed[n_id]) {
+ gid_node_fixed[local_fixed_id] = n_id;
+ for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[n_id][i_face];
+ for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = face_to_node_matrix[face_id][i_node];
+ if ((node_id == n_id) or (is_fixed[node_id])) {
+ continue;
+ } else {
+ non_zeros_fixed[local_fixed_id] += 1;
+ }
+ }
+ }
+ local_fixed_id++;
+ } else {
+ gid_node_free[local_free_id] = n_id;
+ for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[n_id][i_face];
+ for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = face_to_node_matrix[face_id][i_node];
+ if ((node_id == n_id) or (is_fixed[node_id])) {
+ continue;
+ } else {
+ non_zeros_free[local_free_id] += 1;
+ }
+ }
+ }
+ local_free_id++;
+ }
+ }
+ }
+ void
+ _fillMatrix(CRSMatrixDescriptor<double>& Afree,
+ CRSMatrixDescriptor<double>& Afixed,
+ const NodeValue<bool>& is_fixed,
+ const NodeValue<int>& id_free,
+ const NodeValue<int>& id_fixed) const
+ // ,
+ // const Array<NodeId>& gid_node_free,
+ // const Array<NodeId>& gid_node_fixed) const
+ {
+ const auto& node_to_face_matrix = m_given_mesh.connectivity().nodeToFaceMatrix();
+ const auto& face_to_node_matrix = m_given_mesh.connectivity().faceToNodeMatrix();
+ int local_free_id = 0;
+ for (NodeId n_id = 0; n_id < m_given_mesh.numberOfNodes(); n_id++) {
+ if (is_fixed[n_id]) {
+ continue;
+ } else {
+ Assert(id_free[n_id] == local_free_id, "bad matrix definition");
+ for (size_t i_face = 0; i_face < node_to_face_matrix[n_id].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[n_id][i_face];
+ for (size_t i_node = 0; i_node < face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = face_to_node_matrix[face_id][i_node];
+ if (node_id == n_id) {
+ continue;
+ } else if ((is_fixed[node_id])) {
+ Afree(local_free_id, local_free_id) += 1;
+ Afixed(local_free_id, id_fixed[node_id]) -= 1;
+ } else {
+ Afree(local_free_id, local_free_id) += 1;
+ Afree(local_free_id, id_free[node_id]) -= 1;
+ }
+ }
+ }
+ local_free_id++;
+ }
+ }
+ }
+
+ NodeValue<Rd>
+ _getPosition() const
+ {
+ const ConnectivityType& connectivity = m_given_mesh.connectivity();
+ NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+ NodeValue<Rd> pos_r{connectivity};
+ 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]; });
+ _browseBC(is_fixed);
+
+ NodeValue<int> node_dof_id{connectivity};
+ size_t nb_free, nb_fixed;
+ 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_fixed};
+ 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";
+
+ 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()};
+ for (size_t dir = 0; dir < Dimension; ++dir) {
+ for (size_t lid_node = 0; lid_node < nb_fixed; lid_node++) {
+ F[lid_node] = -given_xr[gid_node_fixed[lid_node]][dir];
+ }
+ 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) {
+ if (!is_fixed[node_id]) {
+ pos_r[node_id][dir] = X[id_free[node_id]];
+ }
+ });
+ }
+ return pos_r;
+ }
+
+ public:
+ std::shared_ptr<const IMesh>
+ getSmoothedMesh() const
+ {
+ NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+ NodeValue<Rd> xr = this->_getPosition();
+
+ return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+ }
+
+ ImplicitMeshSmoother(const ImplicitMeshSmoother&) = delete;
+ ImplicitMeshSmoother(ImplicitMeshSmoother&&) = delete;
+
+ ImplicitMeshSmoother(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))
+ {}
+
+ ~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
+{
+ 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 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(
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) 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();
+ }
+ 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();
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
diff --git a/src/mesh/ImplicitMeshSmoother.hpp b/src/mesh/ImplicitMeshSmoother.hpp
new file mode 100644
index 000000000..b5ca016bc
--- /dev/null
+++ b/src/mesh/ImplicitMeshSmoother.hpp
@@ -0,0 +1,65 @@
+#ifndef IMPLICIT_MESH_SMOOTHER_HPP
+#define IMPLICIT_MESH_SMOOTHER_HPP
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
+#include <algebra/LinearSolver.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/DualConnectivityManager.hpp>
+#include <mesh/DualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/FourierBoundaryConditionDescriptor.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
+#include <mesh/IMesh.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+class FunctionSymbolId;
+class IZoneDescriptor;
+class DiscreteFunctionVariant;
+
+class ImplicitMeshSmootherHandler
+{
+ private:
+ template <size_t Dimension>
+ class ImplicitMeshSmoother;
+
+ 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;
+
+ ImplicitMeshSmootherHandler() = default;
+ ImplicitMeshSmootherHandler(ImplicitMeshSmootherHandler&&) = default;
+ ~ImplicitMeshSmootherHandler() = default;
+};
+
+#endif // IMPLICIT_MESH_SMOOTHER_HPP
--
GitLab