From 03c2721760b4bea294531c229223f3dbb50119b6 Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Wed, 23 Dec 2020 19:43:15 +0100
Subject: [PATCH] Put the calculation of the weights out of the schemes and add
an unsteady elasticity scheme
---
src/language/algorithms/CMakeLists.txt | 1 +
src/language/algorithms/ParabolicHeat.cpp | 148 +++-
src/language/algorithms/ParabolicHeat.hpp | 4 +-
.../algorithms/UnsteadyElasticity.cpp | 663 +++++++++++++++
.../algorithms/UnsteadyElasticity.hpp | 37 +
src/language/modules/SchemeModule.cpp | 55 +-
src/scheme/ScalarDiamondScheme.hpp | 106 +--
src/scheme/VectorDiamondScheme.hpp | 785 ++++++++++++++++++
8 files changed, 1697 insertions(+), 102 deletions(-)
create mode 100644 src/language/algorithms/UnsteadyElasticity.cpp
create mode 100644 src/language/algorithms/UnsteadyElasticity.hpp
create mode 100644 src/scheme/VectorDiamondScheme.hpp
diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt
index c1089d014..6292a5830 100644
--- a/src/language/algorithms/CMakeLists.txt
+++ b/src/language/algorithms/CMakeLists.txt
@@ -3,6 +3,7 @@
add_library(PugsLanguageAlgorithms
AcousticSolverAlgorithm.cpp
ElasticityDiamondAlgorithm.cpp
+ UnsteadyElasticity.cpp
Heat5PointsAlgorithm.cpp
HeatDiamondAlgorithm.cpp
HeatDiamondAlgorithm2.cpp
diff --git a/src/language/algorithms/ParabolicHeat.cpp b/src/language/algorithms/ParabolicHeat.cpp
index e4527a682..6ae956ec0 100644
--- a/src/language/algorithms/ParabolicHeat.cpp
+++ b/src/language/algorithms/ParabolicHeat.cpp
@@ -29,6 +29,7 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
+ const FunctionSymbolId& T_init_id,
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id,
const double& Tf,
@@ -237,18 +238,19 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
CellValue<double> Tj =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
CellValue<double> Temperature =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_init_id,
+ mesh_data.xj());
//{mesh->connectivity()};
FaceValue<double> Tl =
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
FaceValue<double> Temperature_face =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_id, mesh_data.xl());
- NodeValue<double> Tr(mesh->connectivity());
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(T_init_id,
+ mesh_data.xl());
{
VTKWriter vtk_writer("Tanalytic_" + std::to_string(Dimension), 0.01);
vtk_writer.write(mesh,
- {NamedItemValue{"Tr", Tr}, NamedItemValue{"temperature", Tj},
+ {NamedItemValue{"exact_temperature", Tj}, NamedItemValue{"init_temperature", Temperature},
NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
@@ -258,6 +260,11 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
Assert(dt > 0, "The time-step have to be positive!");
const CellValue<const double> primal_Vj = mesh_data.Vj();
double output_time = 0;
+
+ InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary);
+ iwm.compute();
+ CellValuePerNode<double> w_rj = iwm.wrj();
+ FaceValuePerNode<double> w_rl = iwm.wrl();
do {
{
if (time >= output_time) {
@@ -271,7 +278,7 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
double deltat = std::min(dt, Tf - time);
std::cout << "Current time = " << time << " time-step = " << deltat << " final time = " << Tf << "\n";
ScalarDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, kappa_id, f_id, Temperature, Temperature_face, Tf,
- deltat);
+ deltat, w_rj, w_rl);
time += deltat;
} while (time < Tf && std::abs(time - Tf) > 1e-15);
{
@@ -447,12 +454,141 @@ class ParabolicHeatScheme<Dimension>::SymmetryBoundaryCondition
~SymmetryBoundaryCondition() = default;
};
+template <size_t Dimension>
+class ParabolicHeatScheme<Dimension>::InterpolationWeightsManager
+{
+ private:
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh;
+ FaceValue<bool> m_primal_face_is_on_boundary;
+ NodeValue<bool> m_primal_node_is_on_boundary;
+ CellValuePerNode<double> m_w_rj;
+ FaceValuePerNode<double> m_w_rl;
+
+ public:
+ InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh,
+ FaceValue<bool> primal_face_is_on_boundary,
+ NodeValue<bool> primal_node_is_on_boundary)
+ : m_mesh(mesh),
+ m_primal_face_is_on_boundary(primal_face_is_on_boundary),
+ m_primal_node_is_on_boundary(primal_node_is_on_boundary)
+ {}
+ ~InterpolationWeightsManager() = default;
+ CellValuePerNode<double>&
+ wrj()
+ {
+ return m_w_rj;
+ }
+ FaceValuePerNode<double>&
+ wrl()
+ {
+ return m_w_rl;
+ }
+ void
+ compute()
+ {
+ using MeshDataType = MeshData<Dimension>;
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
+
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+ const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+ const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_face_matrix = m_mesh->connectivity().nodeToFaceMatrix();
+ CellValuePerNode<double> w_rj{m_mesh->connectivity()};
+ FaceValuePerNode<double> w_rl{m_mesh->connectivity()};
+
+ for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
+ w_rl[i] = std::numeric_limits<double>::signaling_NaN();
+ }
+
+ for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
+ Vector<double> b{Dimension + 1};
+ b[0] = 1;
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ b[i] = xr[i_node][i - 1];
+ }
+ const auto& node_to_cell = node_to_cell_matrix[i_node];
+
+ if (not m_primal_node_is_on_boundary[i_node]) {
+ LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size()};
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ A(0, j) = 1;
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+ Vector<double> x{node_to_cell.size()};
+ x = 0;
+
+ LeastSquareSolver ls_solver;
+ ls_solver.solveLocalSystem(A, x, b);
+
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ w_rj(i_node, j) = x[j];
+ }
+ } else {
+ int nb_face_used = 0;
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ nb_face_used++;
+ }
+ }
+ LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
+ for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
+ A(0, j) = 1;
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ int cpt_face = 0;
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
+ cpt_face++;
+ }
+ }
+ }
+
+ Vector<double> x{node_to_cell.size() + nb_face_used};
+ x = 0;
+
+ LeastSquareSolver ls_solver;
+ ls_solver.solveLocalSystem(A, x, b);
+
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ w_rj(i_node, j) = x[j];
+ }
+ int cpt_face = node_to_cell.size();
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ w_rl(i_node, i_face) = x[cpt_face++];
+ }
+ }
+ }
+ }
+ m_w_rj = w_rj;
+ m_w_rl = w_rl;
+ }
+};
+
template ParabolicHeatScheme<1>::ParabolicHeatScheme(
std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const double&,
const double&,
const double&);
@@ -463,6 +599,7 @@ template ParabolicHeatScheme<2>::ParabolicHeatScheme(
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const double&,
const double&,
const double&);
@@ -473,6 +610,7 @@ template ParabolicHeatScheme<3>::ParabolicHeatScheme(
const FunctionSymbolId&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const double&,
const double&,
const double&);
diff --git a/src/language/algorithms/ParabolicHeat.hpp b/src/language/algorithms/ParabolicHeat.hpp
index b5a2db881..1f8b46e88 100644
--- a/src/language/algorithms/ParabolicHeat.hpp
+++ b/src/language/algorithms/ParabolicHeat.hpp
@@ -17,16 +17,18 @@ class ParabolicHeatScheme
class FourierBoundaryCondition;
class NeumannBoundaryCondition;
class SymmetryBoundaryCondition;
+ class InterpolationWeightsManager;
public:
ParabolicHeatScheme(std::shared_ptr<const IMesh> i_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const FunctionSymbolId& T_id,
+ const FunctionSymbolId& T_init_id,
const FunctionSymbolId& kappa_id,
const FunctionSymbolId& f_id,
const double& final_time,
const double& dt,
- const double&);
+ const double& output_frequency);
};
#endif // HEAT_DIAMOND_ALGORITHM_HPP
diff --git a/src/language/algorithms/UnsteadyElasticity.cpp b/src/language/algorithms/UnsteadyElasticity.cpp
new file mode 100644
index 000000000..d95134359
--- /dev/null
+++ b/src/language/algorithms/UnsteadyElasticity.cpp
@@ -0,0 +1,663 @@
+#include <language/algorithms/UnsteadyElasticity.hpp>
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/LeastSquareSolver.hpp>
+#include <algebra/LinearSolver.hpp>
+#include <algebra/LocalRectangularMatrix.hpp>
+#include <algebra/SparseMatrixDescriptor.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <output/VTKWriter.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/VectorDiamondScheme.hpp>
+
+template <size_t Dimension>
+UnsteadyElasticity<Dimension>::UnsteadyElasticity(
+ std::shared_ptr<const IMesh> i_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& lambda_id,
+ const FunctionSymbolId& mu_id,
+ const FunctionSymbolId& f_id,
+ const FunctionSymbolId& U_id,
+ const FunctionSymbolId& U_init_id,
+ const double& Tf,
+ const double& dt,
+ const double& output_frequency)
+{
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<Dimension>;
+
+ constexpr ItemType FaceType = [] {
+ if constexpr (Dimension > 1) {
+ return ItemType::face;
+ } else {
+ return ItemType::node;
+ }
+ }();
+
+ using BoundaryCondition =
+ std::variant<DirichletBoundaryCondition, NormalStrainBoundaryCondition, SymmetryBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ BoundaryConditionList boundary_condition_list;
+
+ std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ NodeValue<bool> is_dirichlet{mesh->connectivity()};
+ is_dirichlet.fill(false);
+ NodeValue<TinyVector<Dimension>> dirichlet_value{mesh->connectivity()};
+ {
+ TinyVector<Dimension> nan_tiny_vector;
+ for (size_t i = 0; i < Dimension; ++i) {
+ nan_tiny_vector[i] = std::numeric_limits<double>::signaling_NaN();
+ }
+ dirichlet_value.fill(nan_tiny_vector);
+ }
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+ dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>();
+ ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+ if constexpr (Dimension > 1) {
+ Array<const FaceId> face_list = ref_face_list.list();
+ boundary_condition_list.push_back(SymmetryBoundaryCondition{face_list});
+ } else {
+ throw NotImplementedError("Symmetry conditions are not supported in 1d");
+ }
+ }
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (dirichlet_bc_descriptor.name() == "dirichlet") {
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension > 1) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const FaceId> face_list = ref_face_list.list();
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(), face_list);
+ boundary_condition_list.push_back(DirichletBoundaryCondition{face_list, value_list});
+ } else {
+ throw NotImplementedError("Neumann conditions are not supported in 1d");
+ }
+ }
+ }
+ } else if (dirichlet_bc_descriptor.name() == "normal_strain") {
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension > 1) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const FaceId> face_list = ref_face_list.list();
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(), face_list);
+ boundary_condition_list.push_back(NormalStrainBoundaryCondition{face_list, value_list});
+ } else {
+ throw NotImplementedError("Neumann conditions are not supported in 1d");
+ }
+ }
+ }
+ } else {
+ is_valid_boundary_condition = false;
+ }
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for elasticity equation";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ if constexpr (Dimension > 1) {
+ const CellValue<const size_t> cell_dof_number = [&] {
+ CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
+ return compute_cell_dof_number;
+ }();
+ size_t number_of_dof = mesh->numberOfCells();
+
+ const FaceValue<const size_t> face_dof_number = [&] {
+ FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
+ compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>) or
+ (std::is_same_v<T, SymmetryBoundaryCondition>) or
+ (std::is_same_v<T, DirichletBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
+ std::ostringstream os;
+ os << "The face " << face_id << " is used at least twice for boundary conditions";
+ throw NormalError(os.str());
+ } else {
+ compute_face_dof_number[face_id] = number_of_dof++;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return compute_face_dof_number;
+ }();
+
+ const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
+ const FaceValue<const bool> primal_face_is_neumann = [&] {
+ FaceValue<bool> face_is_neumann{mesh->connectivity()};
+ face_is_neumann.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ face_is_neumann[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return face_is_neumann;
+ }();
+
+ FaceValue<bool> primal_face_is_symmetry = [&] {
+ FaceValue<bool> face_is_symmetry{mesh->connectivity()};
+ face_is_symmetry.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, SymmetryBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ face_is_symmetry[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return face_is_symmetry;
+ }();
+
+ NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
+ }
+
+ primal_node_is_on_boundary.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (face_to_cell_matrix[face_id].size() == 1) {
+ for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
+ primal_node_is_on_boundary[node_id] = true;
+ }
+ }
+ }
+
+ FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
+ }
+
+ primal_face_is_on_boundary.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (face_to_cell_matrix[face_id].size() == 1) {
+ primal_face_is_on_boundary[face_id] = true;
+ }
+ }
+
+ FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of face_is_neumann is incorrect");
+ }
+
+ primal_face_is_dirichlet.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id]) &&
+ (!primal_face_is_symmetry[face_id]));
+ }
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+ const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
+
+ {
+ CellValue<TinyVector<Dimension>> velocity = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::cell>(U_id, mesh_data.xj());
+ CellValue<TinyVector<Dimension>> Uj = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::cell>(U_init_id, mesh_data.xj());
+
+ CellValue<TinyVector<Dimension>> fj = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
+ FaceValue<TinyVector<Dimension>> velocity_face = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::face>(U_id, mesh_data.xl());
+ FaceValue<TinyVector<Dimension>> Ul = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::face>(U_init_id, mesh_data.xl());
+
+ {
+ VTKWriter vtk_writer("Uanalytic_" + std::to_string(Dimension), 0.01);
+
+ vtk_writer.write(mesh,
+ {NamedItemValue{"exact_velocity", velocity}, NamedItemValue{"init_velocity", Uj},
+ NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+ NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+ 0, true); // forces last output
+ }
+ double time = 0;
+ Assert(dt > 0, "The time-step have to be positive!");
+ const CellValue<const double> primal_Vj = mesh_data.Vj();
+ double output_time = 0;
+ InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary,
+ primal_face_is_symmetry);
+ iwm.compute();
+ CellValuePerNode<double> w_rj = iwm.wrj();
+ FaceValuePerNode<double> w_rl = iwm.wrl();
+ do {
+ {
+ if (time >= output_time) {
+ VTKWriter vtk_writer("Velocity_" + std::to_string(Dimension) + "_time_" + std::to_string(time), 0.01);
+
+ vtk_writer.write(mesh, {NamedItemValue{"U", Uj}, NamedItemValue{"Vj", primal_Vj}}, 0,
+ true); // forces last output
+ output_time += output_frequency;
+ }
+ }
+ double deltat = std::min(dt, Tf - time);
+ std::cout << "Current time = " << time << " time-step = " << deltat << " final time = " << Tf << "\n";
+ VectorDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, Uj, Ul, Tf, deltat, w_rj,
+ w_rl);
+ time += deltat;
+ } while (time < Tf && std::abs(time - Tf) > 1e-15);
+
+ const FaceValue<const double> mes_l = [&] {
+ if constexpr (Dimension == 1) {
+ FaceValue<double> compute_mes_l{mesh->connectivity()};
+ compute_mes_l.fill(1);
+ return compute_mes_l;
+ } else {
+ return mesh_data.ll();
+ }
+ }();
+
+ Vector<double> Uexacte{mesh->numberOfCells() * Dimension};
+ for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+ for (size_t l = 0; l < Dimension; ++l) {
+ Uexacte[(cell_dof_number[j] * Dimension) + l] = velocity[j][l];
+ }
+ }
+
+ Vector<double> error{mesh->numberOfCells() * Dimension};
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ error[(cell_id * Dimension) + i] = (Uj[cell_id][i] - velocity[cell_id][i]) * sqrt(primal_Vj[cell_id]);
+ }
+ }
+ Vector<double> error_face{mesh->numberOfFaces() * Dimension};
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ if (primal_face_is_on_boundary[face_id]) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ error_face[face_id * Dimension + i] = (Ul[face_id][i] - velocity_face[face_id][i]) * sqrt(mes_l[face_id]);
+ }
+ } else {
+ error_face[face_id] = 0;
+ }
+ });
+
+ std::cout << "||Error||_2 (cell)= " << std::sqrt((error, error)) << "\n";
+ std::cout << "||Error||_2 (face)= " << std::sqrt((error_face, error_face)) << "\n";
+ std::cout << "||Error||_2 (total)= " << std::sqrt((error, error)) + std::sqrt((error_face, error_face)) << "\n";
+
+ VTKWriter vtk_writer("Speed_" + std::to_string(Dimension), 0.01);
+
+ NodeValue<TinyVector<3>> ur3d{mesh->connectivity()};
+ ur3d.fill(zero);
+
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ TinyVector<Dimension> x = zero;
+ const auto node_cells = node_to_cell_matrix[node_id];
+ for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
+ CellId cell_id = node_cells[i_cell];
+ x += w_rj(node_id, i_cell) * Uj[cell_id];
+ }
+ const auto node_faces = node_to_face_matrix[node_id];
+ for (size_t i_face = 0; i_face < node_faces.size(); ++i_face) {
+ FaceId face_id = node_faces[i_face];
+ if (primal_face_is_on_boundary[face_id]) {
+ x += w_rl(node_id, i_face) * Ul[face_id];
+ }
+ }
+ for (size_t i = 0; i < Dimension; ++i) {
+ ur3d[node_id][i] = x[i];
+ }
+ });
+
+ vtk_writer.write(mesh,
+ {NamedItemValue{"U", Uj}, NamedItemValue{"Ur3d", ur3d}, NamedItemValue{"Uexacte", velocity}}, 0,
+ true); // forces last output
+ }
+ } else {
+ throw NotImplementedError("not done in 1d");
+ }
+}
+
+template <size_t Dimension>
+class UnsteadyElasticity<Dimension>::DirichletBoundaryCondition
+{
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const TinyVector<Dimension>>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~DirichletBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class UnsteadyElasticity<Dimension>::NormalStrainBoundaryCondition
+{
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ NormalStrainBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const TinyVector<Dimension>>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~NormalStrainBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class UnsteadyElasticity<Dimension>::SymmetryBoundaryCondition
+{
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ public:
+ SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class UnsteadyElasticity<Dimension>::InterpolationWeightsManager
+{
+ private:
+ std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh;
+ FaceValue<bool> m_primal_face_is_on_boundary;
+ NodeValue<bool> m_primal_node_is_on_boundary;
+ FaceValue<bool> m_primal_face_is_symmetry;
+ CellValuePerNode<double> m_w_rj;
+ FaceValuePerNode<double> m_w_rl;
+
+ public:
+ InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh,
+ FaceValue<bool> primal_face_is_on_boundary,
+ NodeValue<bool> primal_node_is_on_boundary,
+ FaceValue<bool> primal_face_is_symmetry)
+ : m_mesh(mesh),
+ m_primal_face_is_on_boundary(primal_face_is_on_boundary),
+ m_primal_node_is_on_boundary(primal_node_is_on_boundary),
+ m_primal_face_is_symmetry(primal_face_is_symmetry)
+ {}
+ ~InterpolationWeightsManager() = default;
+ CellValuePerNode<double>&
+ wrj()
+ {
+ return m_w_rj;
+ }
+ FaceValuePerNode<double>&
+ wrl()
+ {
+ return m_w_rl;
+ }
+ void
+ compute()
+ {
+ using MeshDataType = MeshData<Dimension>;
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+ const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr();
+
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+ const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+ const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_face_matrix = m_mesh->connectivity().nodeToFaceMatrix();
+ const auto& face_to_cell_matrix = m_mesh->connectivity().faceToCellMatrix();
+
+ CellValuePerNode<double> w_rj{m_mesh->connectivity()};
+ FaceValuePerNode<double> w_rl{m_mesh->connectivity()};
+
+ const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
+ auto project_to_face = [&](const TinyVector<Dimension>& x, const FaceId face_id) -> const TinyVector<Dimension> {
+ TinyVector<Dimension> proj;
+ const TinyVector<Dimension> nil = primal_nlr(face_id, 0);
+ proj = x - ((x - xl[face_id]), nil) * nil;
+ return proj;
+ };
+
+ for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
+ w_rl[i] = std::numeric_limits<double>::signaling_NaN();
+ }
+
+ for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) {
+ Vector<double> b{Dimension + 1};
+ b[0] = 1;
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ b[i] = xr[i_node][i - 1];
+ }
+ const auto& node_to_cell = node_to_cell_matrix[i_node];
+
+ if (not m_primal_node_is_on_boundary[i_node]) {
+ LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size()};
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ A(0, j) = 1;
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+
+ Vector<double> x{node_to_cell.size()};
+ x = 0;
+
+ LeastSquareSolver ls_solver;
+ ls_solver.solveLocalSystem(A, x, b);
+
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ w_rj(i_node, j) = x[j];
+ }
+ } else {
+ int nb_face_used = 0;
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ nb_face_used++;
+ }
+ }
+ LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
+ for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
+ A(0, j) = 1;
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ const CellId J = node_to_cell[j];
+ A(i, j) = xj[J][i - 1];
+ }
+ }
+ for (size_t i = 1; i < Dimension + 1; i++) {
+ int cpt_face = 0;
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ if (m_primal_face_is_symmetry[face_id]) {
+ for (size_t j = 0; j < face_to_cell_matrix[face_id].size(); ++j) {
+ const CellId cell_id = face_to_cell_matrix[face_id][j];
+ TinyVector<Dimension> xproj = project_to_face(xj[cell_id], face_id);
+ A(i, node_to_cell.size() + cpt_face) = xproj[i - 1];
+ }
+ } else {
+ A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
+ }
+ cpt_face++;
+ }
+ }
+ }
+
+ Vector<double> x{node_to_cell.size() + nb_face_used};
+ x = 0;
+
+ LeastSquareSolver ls_solver;
+ ls_solver.solveLocalSystem(A, x, b);
+
+ for (size_t j = 0; j < node_to_cell.size(); j++) {
+ w_rj(i_node, j) = x[j];
+ }
+ int cpt_face = node_to_cell.size();
+ for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
+ FaceId face_id = node_to_face_matrix[i_node][i_face];
+ if (m_primal_face_is_on_boundary[face_id]) {
+ w_rl(i_node, i_face) = x[cpt_face++];
+ }
+ }
+ }
+ }
+ m_w_rj = w_rj;
+ m_w_rl = w_rl;
+ }
+};
+
+template UnsteadyElasticity<1>::UnsteadyElasticity(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const double&,
+ const double&,
+ const double&);
+
+template UnsteadyElasticity<2>::UnsteadyElasticity(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const double&,
+ const double&,
+ const double&);
+
+template UnsteadyElasticity<3>::UnsteadyElasticity(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const double&,
+ const double&,
+ const double&);
diff --git a/src/language/algorithms/UnsteadyElasticity.hpp b/src/language/algorithms/UnsteadyElasticity.hpp
new file mode 100644
index 000000000..adbb1ad95
--- /dev/null
+++ b/src/language/algorithms/UnsteadyElasticity.hpp
@@ -0,0 +1,37 @@
+#ifndef UNSTEADY_ELASTICITY_HPP
+#define UNSTEADY_ELASTICITY_HPP
+
+#include <algebra/TinyVector.hpp>
+#include <language/utils/FunctionSymbolId.hpp>
+#include <memory>
+#include <mesh/IMesh.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <variant>
+#include <vector>
+
+template <size_t Dimension>
+class UnsteadyElasticity
+{
+ private:
+ class DirichletBoundaryCondition;
+ class NormalStrainBoundaryCondition;
+ class SymmetryBoundaryCondition;
+ class InterpolationWeightsManager;
+
+ public:
+ UnsteadyElasticity(std::shared_ptr<const IMesh> i_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& lambda_id,
+ const FunctionSymbolId& mu_id,
+ const FunctionSymbolId& f_id,
+ const FunctionSymbolId& U_id,
+ const FunctionSymbolId& U_init_id,
+ const double& Tf,
+ const double& dt,
+ const double& output_frequency);
+};
+
+#endif // ELASTICITY_DIAMOND_ALGORITHM2_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index e32f667e7..714d9ae1e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -6,6 +6,7 @@
#include <language/algorithms/HeatDiamondAlgorithm.hpp>
#include <language/algorithms/HeatDiamondAlgorithm2.hpp>
#include <language/algorithms/ParabolicHeat.hpp>
+#include <language/algorithms/UnsteadyElasticity.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/Mesh.hpp>
@@ -267,27 +268,27 @@ SchemeModule::SchemeModule()
void(std::shared_ptr<const IMesh>,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
- const double&, const double&, const double&)>>(
+ const FunctionSymbolId&, const double&, const double&, const double&)>>(
[](std::shared_ptr<const IMesh> p_mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
- const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id,
- const FunctionSymbolId& f_id, const double& final_time, const double& dt,
- const double& output_frequency) -> void {
+ const FunctionSymbolId& T_id, const FunctionSymbolId& T_init_id,
+ const FunctionSymbolId& kappa_id, const FunctionSymbolId& f_id,
+ const double& final_time, const double& dt, const double& output_frequency) -> void {
switch (p_mesh->dimension()) {
case 1: {
- ParabolicHeatScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id,
+ ParabolicHeatScheme<1>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
f_id, final_time, dt, output_frequency};
break;
}
case 2: {
- ParabolicHeatScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id,
+ ParabolicHeatScheme<2>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
f_id, final_time, dt, output_frequency};
break;
}
case 3: {
- ParabolicHeatScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id,
+ ParabolicHeatScheme<3>{p_mesh, bc_descriptor_list, T_id, T_init_id, kappa_id,
f_id, final_time, dt, output_frequency};
break;
}
@@ -299,6 +300,46 @@ SchemeModule::SchemeModule()
));
+ this
+ ->_addBuiltinFunction("unsteadyelasticity",
+ std::make_shared<BuiltinFunctionEmbedder<
+ void(std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&,
+ const FunctionSymbolId&, const FunctionSymbolId&, const double&, const double&,
+ const double&)>>(
+
+ [](std::shared_ptr<const IMesh> p_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id,
+ const FunctionSymbolId& f_id, const FunctionSymbolId& U_id,
+ const FunctionSymbolId& U_init_id, const double& final_time, const double& dt,
+ const double& output_frequency) -> void {
+ switch (p_mesh->dimension()) {
+ case 1: {
+ UnsteadyElasticity<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id,
+ U_id, U_init_id, final_time, dt, output_frequency};
+ break;
+ }
+ case 2: {
+ UnsteadyElasticity<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id,
+ U_id, U_init_id, final_time, dt, output_frequency};
+ break;
+ }
+ case 3: {
+ UnsteadyElasticity<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id,
+ U_id, U_init_id, final_time, dt, output_frequency};
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+ }
+
+ ));
+
this->_addBuiltinFunction("heat2",
std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
diff --git a/src/scheme/ScalarDiamondScheme.hpp b/src/scheme/ScalarDiamondScheme.hpp
index 7259c8442..70c375b3d 100644
--- a/src/scheme/ScalarDiamondScheme.hpp
+++ b/src/scheme/ScalarDiamondScheme.hpp
@@ -1,3 +1,6 @@
+#ifndef SCALAR_DIAMOND_SCHEME_HPP
+#define SCALAR_DIAMOND_SCHEME_HPP
+
#include <algebra/CRSMatrix.hpp>
#include <algebra/LeastSquareSolver.hpp>
#include <algebra/LinearSolver.hpp>
@@ -146,7 +149,9 @@ class ScalarDiamondScheme
CellValue<double>& Temperature,
FaceValue<double>& Temperature_face,
const double& Tf,
- const double& dt)
+ const double& dt,
+ const CellValuePerNode<double>& w_rj,
+ const FaceValuePerNode<double>& w_rl)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
@@ -343,94 +348,9 @@ class ScalarDiamondScheme
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
- const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
-
const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
- const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
- CellValuePerNode<double> w_rj{mesh->connectivity()};
- FaceValuePerNode<double> w_rl{mesh->connectivity()};
-
- for (size_t i = 0; i < w_rl.numberOfValues(); ++i) {
- w_rl[i] = std::numeric_limits<double>::signaling_NaN();
- }
-
- for (NodeId i_node = 0; i_node < mesh->numberOfNodes(); ++i_node) {
- Vector<double> b{Dimension + 1};
- b[0] = 1;
- for (size_t i = 1; i < Dimension + 1; i++) {
- b[i] = xr[i_node][i - 1];
- }
- const auto& node_to_cell = node_to_cell_matrix[i_node];
-
- if (not primal_node_is_on_boundary[i_node]) {
- LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size()};
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- A(0, j) = 1;
- }
- for (size_t i = 1; i < Dimension + 1; i++) {
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- const CellId J = node_to_cell[j];
- A(i, j) = xj[J][i - 1];
- }
- }
- Vector<double> x{node_to_cell.size()};
- x = 0;
-
- LeastSquareSolver ls_solver;
- ls_solver.solveLocalSystem(A, x, b);
-
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- w_rj(i_node, j) = x[j];
- }
- } else {
- int nb_face_used = 0;
- for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
- FaceId face_id = node_to_face_matrix[i_node][i_face];
- if (primal_face_is_on_boundary[face_id]) {
- nb_face_used++;
- }
- }
- LocalRectangularMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used};
- for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) {
- A(0, j) = 1;
- }
- for (size_t i = 1; i < Dimension + 1; i++) {
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- const CellId J = node_to_cell[j];
- A(i, j) = xj[J][i - 1];
- }
- }
- for (size_t i = 1; i < Dimension + 1; i++) {
- int cpt_face = 0;
- for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
- FaceId face_id = node_to_face_matrix[i_node][i_face];
- if (primal_face_is_on_boundary[face_id]) {
- A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1];
- cpt_face++;
- }
- }
- }
-
- Vector<double> x{node_to_cell.size() + nb_face_used};
- x = 0;
-
- LeastSquareSolver ls_solver;
- ls_solver.solveLocalSystem(A, x, b);
-
- for (size_t j = 0; j < node_to_cell.size(); j++) {
- w_rj(i_node, j) = x[j];
- }
- int cpt_face = node_to_cell.size();
- for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) {
- FaceId face_id = node_to_face_matrix[i_node][i_face];
- if (primal_face_is_on_boundary[face_id]) {
- w_rl(i_node, i_face) = x[cpt_face++];
- }
- }
- }
- }
{
std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(mesh);
@@ -722,7 +642,9 @@ template ScalarDiamondScheme<1>::ScalarDiamondScheme(
CellValue<double>&,
FaceValue<double>&,
const double&,
- const double&);
+ const double&,
+ const CellValuePerNode<double>& w_rj,
+ const FaceValuePerNode<double>& w_rl);
template ScalarDiamondScheme<2>::ScalarDiamondScheme(
std::shared_ptr<const IMesh>,
@@ -732,7 +654,9 @@ template ScalarDiamondScheme<2>::ScalarDiamondScheme(
CellValue<double>&,
FaceValue<double>&,
const double&,
- const double&);
+ const double&,
+ const CellValuePerNode<double>& w_rj,
+ const FaceValuePerNode<double>& w_rl);
template ScalarDiamondScheme<3>::ScalarDiamondScheme(
std::shared_ptr<const IMesh>,
@@ -742,4 +666,8 @@ template ScalarDiamondScheme<3>::ScalarDiamondScheme(
CellValue<double>&,
FaceValue<double>&,
const double&,
- const double&);
+ const double&,
+ const CellValuePerNode<double>& w_rj,
+ const FaceValuePerNode<double>& w_rl);
+
+#endif // SCALAR_DIAMOND_SCHEME_HPP
diff --git a/src/scheme/VectorDiamondScheme.hpp b/src/scheme/VectorDiamondScheme.hpp
new file mode 100644
index 000000000..bb67cd15d
--- /dev/null
+++ b/src/scheme/VectorDiamondScheme.hpp
@@ -0,0 +1,785 @@
+#ifndef VECTOR_DIAMOND_SCHEME_HPP
+#define VECTOR_DIAMOND_SCHEME_HPP
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/LeastSquareSolver.hpp>
+#include <algebra/LinearSolver.hpp>
+#include <algebra/LocalRectangularMatrix.hpp>
+#include <algebra/SparseMatrixDescriptor.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityToDiamondDualConnectivityDataMapper.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <output/VTKWriter.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
+template <size_t Dimension>
+class VectorDiamondScheme
+{
+ private:
+ class DirichletBoundaryCondition
+ {
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ DirichletBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const TinyVector<Dimension>>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~DirichletBoundaryCondition() = default;
+ };
+
+ class NormalStrainBoundaryCondition
+ {
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ NormalStrainBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const TinyVector<Dimension>>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}
+ {
+ Assert(m_value_list.size() == m_face_list.size());
+ }
+
+ ~NormalStrainBoundaryCondition() = default;
+ };
+
+ class SymmetryBoundaryCondition
+ {
+ private:
+ const Array<const TinyVector<Dimension>> m_value_list;
+ const Array<const FaceId> m_face_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_face_list;
+ }
+
+ public:
+ SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {}
+
+ ~SymmetryBoundaryCondition() = default;
+ };
+
+ public:
+ VectorDiamondScheme(std::shared_ptr<const IMesh> i_mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const FunctionSymbolId& lambda_id,
+ const FunctionSymbolId& mu_id,
+ const FunctionSymbolId& f_id,
+ CellValue<TinyVector<Dimension>>& Uj,
+ FaceValue<TinyVector<Dimension>>& Ul,
+ const double& Tf,
+ const double& dt,
+ CellValuePerNode<double>& w_rj,
+ FaceValuePerNode<double>& w_rl)
+ {
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<Dimension>;
+
+ constexpr ItemType FaceType = [] {
+ if constexpr (Dimension > 1) {
+ return ItemType::face;
+ } else {
+ return ItemType::node;
+ }
+ }();
+
+ using BoundaryCondition =
+ std::variant<DirichletBoundaryCondition, NormalStrainBoundaryCondition, SymmetryBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ BoundaryConditionList boundary_condition_list;
+
+ std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
+ std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+ NodeValue<bool> is_dirichlet{mesh->connectivity()};
+ is_dirichlet.fill(false);
+ NodeValue<TinyVector<Dimension>> dirichlet_value{mesh->connectivity()};
+ {
+ TinyVector<Dimension> nan_tiny_vector;
+ for (size_t i = 0; i < Dimension; ++i) {
+ nan_tiny_vector[i] = std::numeric_limits<double>::signaling_NaN();
+ }
+ dirichlet_value.fill(nan_tiny_vector);
+ }
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+ dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+ if constexpr (Dimension > 1) {
+ Array<const FaceId> face_list = ref_face_list.list();
+ boundary_condition_list.push_back(SymmetryBoundaryCondition{face_list});
+ } else {
+ throw NotImplementedError("Symmetry conditions are not supported in 1d");
+ }
+ }
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (dirichlet_bc_descriptor.name() == "dirichlet") {
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension > 1) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const FaceId> face_list = ref_face_list.list();
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(), face_list);
+ boundary_condition_list.push_back(DirichletBoundaryCondition{face_list, value_list});
+ } else {
+ throw NotImplementedError("Neumann conditions are not supported in 1d");
+ }
+ }
+ }
+ } else if (dirichlet_bc_descriptor.name() == "normal_strain") {
+ for (size_t i_ref_face_list = 0;
+ i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+ const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+ const RefId& ref = ref_face_list.refId();
+ if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+ const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+ if constexpr (Dimension > 1) {
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ Array<const FaceId> face_list = ref_face_list.list();
+ Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(), face_list);
+ boundary_condition_list.push_back(NormalStrainBoundaryCondition{face_list, value_list});
+ } else {
+ throw NotImplementedError("Neumann conditions are not supported in 1d");
+ }
+ }
+ }
+ } else {
+ is_valid_boundary_condition = false;
+ }
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for elasticity equation";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ if constexpr (Dimension > 1) {
+ const CellValue<const size_t> cell_dof_number = [&] {
+ CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
+ return compute_cell_dof_number;
+ }();
+ size_t number_of_dof = mesh->numberOfCells();
+
+ const FaceValue<const size_t> face_dof_number = [&] {
+ FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
+ compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>) or
+ (std::is_same_v<T, SymmetryBoundaryCondition>) or
+ (std::is_same_v<T, DirichletBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
+ std::ostringstream os;
+ os << "The face " << face_id << " is used at least twice for boundary conditions";
+ throw NormalError(os.str());
+ } else {
+ compute_face_dof_number[face_id] = number_of_dof++;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return compute_face_dof_number;
+ }();
+
+ const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
+ const FaceValue<const bool> primal_face_is_neumann = [&] {
+ FaceValue<bool> face_is_neumann{mesh->connectivity()};
+ face_is_neumann.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ face_is_neumann[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return face_is_neumann;
+ }();
+
+ const FaceValue<const bool> primal_face_is_symmetry = [&] {
+ FaceValue<bool> face_is_symmetry{mesh->connectivity()};
+ face_is_symmetry.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, SymmetryBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ face_is_symmetry[face_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ return face_is_symmetry;
+ }();
+
+ NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of node_is_on_boundary is incorrect");
+ }
+
+ primal_node_is_on_boundary.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (face_to_cell_matrix[face_id].size() == 1) {
+ for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
+ primal_node_is_on_boundary[node_id] = true;
+ }
+ }
+ }
+
+ FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of face_is_on_boundary is incorrect");
+ }
+
+ primal_face_is_on_boundary.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (face_to_cell_matrix[face_id].size() == 1) {
+ primal_face_is_on_boundary[face_id] = true;
+ }
+ }
+
+ FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity());
+ if (parallel::size() > 1) {
+ throw NotImplementedError("Calculation of face_is_neumann is incorrect");
+ }
+
+ primal_face_is_dirichlet.fill(false);
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] &&
+ (!primal_face_is_neumann[face_id]) && (!primal_face_is_symmetry[face_id]));
+ }
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+ const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+ const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
+ const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
+
+ {
+ std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(mesh);
+
+ MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh);
+
+ std::shared_ptr mapper =
+ DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
+ mesh->connectivity());
+
+ CellValue<double> dual_muj =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(mu_id,
+ diamond_mesh_data
+ .xj());
+
+ CellValue<double> dual_lambdaj =
+ InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(lambda_id,
+ diamond_mesh_data
+ .xj());
+
+ CellValue<TinyVector<Dimension>> fj = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
+
+ const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
+
+ const FaceValue<const double> mes_l = [&] {
+ if constexpr (Dimension == 1) {
+ FaceValue<double> compute_mes_l{mesh->connectivity()};
+ compute_mes_l.fill(1);
+ return compute_mes_l;
+ } else {
+ return mesh_data.ll();
+ }
+ }();
+
+ const CellValue<const double> dual_mes_l_j = [=] {
+ CellValue<double> compute_mes_j{diamond_mesh->connectivity()};
+ mapper->toDualCell(mes_l, compute_mes_j);
+
+ return compute_mes_j;
+ }();
+
+ const CellValue<const double> primal_Vj = mesh_data.Vj();
+ FaceValue<const CellId> face_dual_cell_id = [=]() {
+ FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()};
+ CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()};
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
+
+ mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id);
+
+ return computed_face_dual_cell_id;
+ }();
+
+ NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
+ CellValue<NodeId> cell_ignored_id{mesh->connectivity()};
+ cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
+
+ NodeValue<NodeId> node_primal_id{mesh->connectivity()};
+
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
+
+ NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
+
+ mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
+
+ return computed_dual_node_primal_node_id;
+ }();
+
+ CellValue<NodeId> primal_cell_dual_node_id = [=]() {
+ CellValue<NodeId> cell_id{mesh->connectivity()};
+ NodeValue<NodeId> node_ignored_id{mesh->connectivity()};
+ node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
+
+ NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; });
+
+ CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()};
+
+ mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id);
+
+ return cell_id;
+ }();
+ const auto& dual_Cjr = diamond_mesh_data.Cjr();
+ FaceValue<TinyVector<Dimension>> dualClj = [&] {
+ FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()};
+ const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix();
+ const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
+ parallel_for(
+ mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+ for (size_t i = 0; i < primal_face_to_cell.size(); i++) {
+ CellId cell_id = primal_face_to_cell[i];
+ const NodeId dual_node_id = primal_cell_dual_node_id[cell_id];
+ for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size();
+ i_dual_cell++) {
+ const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell];
+ if (face_dual_cell_id[face_id] == dual_cell_id) {
+ for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
+ i_dual_node++) {
+ const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+ if (final_dual_node_id == dual_node_id) {
+ computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node);
+ }
+ }
+ }
+ }
+ }
+ });
+ return computedClj;
+ }();
+
+ FaceValue<TinyVector<Dimension>> nlj = [&] {
+ FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()};
+ parallel_for(
+ mesh->numberOfFaces(),
+ PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; });
+ return computedNlj;
+ }();
+
+ FaceValue<const double> alpha_lambda_l = [&] {
+ CellValue<double> alpha_j{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+ alpha_j[diamond_cell_id] = dual_lambdaj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ });
+
+ FaceValue<double> computed_alpha_l{mesh->connectivity()};
+ mapper->fromDualCell(alpha_j, computed_alpha_l);
+ return computed_alpha_l;
+ }();
+
+ FaceValue<const double> alpha_mu_l = [&] {
+ CellValue<double> alpha_j{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+ alpha_j[diamond_cell_id] = dual_muj[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ });
+
+ FaceValue<double> computed_alpha_l{mesh->connectivity()};
+ mapper->fromDualCell(alpha_j, computed_alpha_l);
+ return computed_alpha_l;
+ }();
+
+ double unsteady = (Tf == 0) ? 0 : 1;
+ double time_factor = (unsteady == 0) ? 1 : dt;
+
+ TinyMatrix<Dimension, double> eye = identity;
+
+ SparseMatrixDescriptor S{number_of_dof * Dimension};
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ const double beta_mu_l = l2Norm(dualClj[face_id]) * alpha_mu_l[face_id] * mes_l[face_id];
+ const double beta_lambda_l = l2Norm(dualClj[face_id]) * alpha_lambda_l[face_id] * mes_l[face_id];
+ const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
+ for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
+ const CellId i_id = primal_face_to_cell[i_cell];
+ const bool is_face_reversed_for_cell_i = ((dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
+
+ const TinyVector<Dimension> nil = [&] {
+ if (is_face_reversed_for_cell_i) {
+ return -nlj[face_id];
+ } else {
+ return nlj[face_id];
+ }
+ }();
+ for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
+ const CellId j_id = primal_face_to_cell[j_cell];
+ TinyMatrix<Dimension, double> M =
+ beta_mu_l * eye + beta_mu_l * tensorProduct(nil, nil) + beta_lambda_l * tensorProduct(nil, nil);
+ TinyMatrix<Dimension, double> N = tensorProduct(nil, nil);
+
+ if (i_cell == j_cell) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) +=
+ (time_factor * M(i, j) + unsteady * primal_Vj[i_id]);
+ if (primal_face_is_neumann[face_id]) {
+ S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) -= M(i, j);
+ }
+ if (primal_face_is_symmetry[face_id]) {
+ S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
+ -((i == j) ? 1 : 0) + N(i, j);
+ S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + j) +=
+ (i == j) ? 1 : 0;
+ }
+ }
+ }
+ } else {
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -=
+ time_factor * M(i, j);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
+ const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ const double alpha_mu_face_id = mes_l[face_id] * alpha_mu_l[face_id];
+ const double alpha_lambda_face_id = mes_l[face_id] * alpha_lambda_l[face_id];
+
+ for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
+ CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
+ const bool is_face_reversed_for_cell_i = ((dualClj[face_id], xl[face_id] - xj[i_id]) < 0);
+
+ for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
+ NodeId node_id = primal_face_to_node_matrix[face_id][i_node];
+
+ const TinyVector<Dimension> nil = [&] {
+ if (is_face_reversed_for_cell_i) {
+ return -nlj[face_id];
+ } else {
+ return nlj[face_id];
+ }
+ }();
+
+ CellId dual_cell_id = face_dual_cell_id[face_id];
+
+ for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
+ const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
+ if (dual_node_primal_node_id[dual_node_id] == node_id) {
+ const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
+
+ TinyMatrix<Dimension, double> M = alpha_mu_face_id * (Clr, nil) * eye +
+ alpha_mu_face_id * tensorProduct(Clr, nil) +
+ alpha_lambda_face_id * tensorProduct(nil, Clr);
+
+ for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
+ CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -=
+ time_factor * w_rj(node_id, j_cell) * M(i, j);
+ if (primal_face_is_neumann[face_id]) {
+ S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
+ w_rj(node_id, j_cell) * M(i, j);
+ }
+ }
+ }
+ }
+ if (primal_node_is_on_boundary[node_id]) {
+ for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
+ FaceId l_id = node_to_face_matrix[node_id][l_face];
+ if (primal_face_is_on_boundary[l_id]) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S(cell_dof_number[i_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) -=
+ time_factor * w_rl(node_id, l_face) * M(i, j);
+ }
+ }
+ if (primal_face_is_neumann[face_id]) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ S(face_dof_number[face_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) +=
+ w_rl(node_id, l_face) * M(i, j);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ // }
+ }
+ }
+ }
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ if (primal_face_is_dirichlet[face_id]) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + i) += 1;
+ }
+ }
+ }
+
+ Vector<double> b{number_of_dof * Dimension};
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ b[(cell_dof_number[cell_id] * Dimension) + i] =
+ primal_Vj[cell_id] * (time_factor * fj[cell_id][i] + unsteady * Uj[cell_id][i]);
+ }
+ }
+
+ // Dirichlet
+ NodeValue<bool> node_tag{mesh->connectivity()};
+ node_tag.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ for (size_t i = 0; i < Dimension; ++i) {
+ b[(face_dof_number[face_id] * Dimension) + i] += value_list[i_face][i];
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>)) {
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ FaceId face_id = face_list[i_face];
+ for (size_t i = 0; i < Dimension; ++i) {
+ b[face_dof_number[face_id] * Dimension + i] += mes_l[face_id] * value_list[i_face][i]; // sign
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ CRSMatrix A{S};
+ Vector<double> U{number_of_dof * Dimension};
+ U = 1;
+ Vector r = A * U - b;
+ std::cout << "initial (real) residu = " << std::sqrt((r, r)) << '\n';
+
+ LinearSolver solver;
+ solver.solveLocalSystem(A, U, b);
+
+ r = A * U - b;
+
+ std::cout << "final (real) residu = " << std::sqrt((r, r)) << '\n';
+
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ Uj[cell_id][i] = U[(cell_dof_number[cell_id] * Dimension) + i];
+ }
+ }
+ for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ if (primal_face_is_on_boundary[face_id]) {
+ Ul[face_id][i] = U[(face_dof_number[face_id] * Dimension) + i];
+ }
+ }
+ }
+ NodeValue<TinyVector<3>> ur3d{mesh->connectivity()};
+ ur3d.fill(zero);
+
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ TinyVector<Dimension> x = zero;
+ const auto node_cells = node_to_cell_matrix[node_id];
+ for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
+ CellId cell_id = node_cells[i_cell];
+ x += w_rj(node_id, i_cell) * Uj[cell_id];
+ }
+ const auto node_faces = node_to_face_matrix[node_id];
+ for (size_t i_face = 0; i_face < node_faces.size(); ++i_face) {
+ FaceId face_id = node_faces[i_face];
+ if (primal_face_is_on_boundary[face_id]) {
+ x += w_rl(node_id, i_face) * Ul[face_id];
+ }
+ }
+ for (size_t i = 0; i < Dimension; ++i) {
+ ur3d[node_id][i] = x[i];
+ }
+ });
+ }
+ } else {
+ throw NotImplementedError("not done in 1d");
+ }
+ }
+};
+template VectorDiamondScheme<1>::VectorDiamondScheme(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ CellValue<TinyVector<1>>&,
+ FaceValue<TinyVector<1>>&,
+ const double&,
+ const double&,
+ CellValuePerNode<double>&,
+ FaceValuePerNode<double>&);
+
+template VectorDiamondScheme<2>::VectorDiamondScheme(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ CellValue<TinyVector<2>>&,
+ FaceValue<TinyVector<2>>&,
+ const double&,
+ const double&,
+ CellValuePerNode<double>&,
+ FaceValuePerNode<double>&);
+
+template VectorDiamondScheme<3>::VectorDiamondScheme(
+ std::shared_ptr<const IMesh>,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ const FunctionSymbolId&,
+ CellValue<TinyVector<3>>&,
+ FaceValue<TinyVector<3>>&,
+ const double&,
+ const double&,
+ CellValuePerNode<double>&,
+ FaceValuePerNode<double>&);
+
+#endif // VECTOR_DIAMOND_SCHEME_HPP
--
GitLab