Skip to content
Snippets Groups Projects
Commit 45736307 authored by Emmanuel Labourasse's avatar Emmanuel Labourasse
Browse files

Add contraction operator to TinyMatrix and an implicit mesh smoother

parent 40538c11
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/mesh/ImplicitMeshSmoother.cpp 0 → 100644
+ 463
− 0
View file @ 45736307
#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");
}
}
}
src/mesh/ImplicitMeshSmoother.hpp 0 → 100644
+ 65
− 0
View file @ 45736307
#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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment