From 4f5b355eee3c5575cb439abf01351820ca5bbef2 Mon Sep 17 00:00:00 2001
From: Julie PATELA <julie.patela.ocre@cea.fr>
Date: Thu, 20 Aug 2020 09:23:36 +0200
Subject: [PATCH] Add internal's node matrix and Dirichlet's boundary condition
---
.../algorithms/ElasticityDiamondAlgorithm.cpp | 259 +++++++++++++++---
.../algorithms/ElasticityDiamondAlgorithm.hpp | 3 +-
src/language/modules/SchemeModule.cpp | 11 +-
3 files changed, 227 insertions(+), 46 deletions(-)
diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
index f51e7678e..e7e8de72e 100644
--- a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
+++ b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
@@ -27,7 +27,8 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
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& f_id,
+ const FunctionSymbolId& U_id)
{
using ConnectivityType = Connectivity<Dimension>;
using MeshType = Mesh<ConnectivityType>;
@@ -359,6 +360,9 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
diamond_mesh_data
.xj());
+ CellValue<TinyVector<Dimension>> Uj = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::cell>(U_id, mesh_data.xj());
+
CellValue<TinyVector<Dimension>> fj = InterpolateItemValue<TinyVector<Dimension>(
TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
@@ -427,6 +431,22 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
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;
+ }();
+
TinyMatrix<Dimension, double> eye = zero;
for (size_t i = 0; i < Dimension; ++i) {
eye(i, i) = 1;
@@ -437,35 +457,84 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells();
SparseMatrixDescriptor S{number_of_dof * Dimension};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- if (not primal_face_is_on_boundary[face_id]) {
- const double beta_mu_l = (1. / Dimension) * alpha_mu_l[face_id] * mes_l[face_id];
- const double beta_lambda_l = (1. / Dimension) * 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 =
- primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_cell));
- const TinyVector<Dimension> nil = [&] {
- if (is_face_reversed_for_cell_i) {
- return -primal_nlr(face_id, 0);
- } else {
- return primal_nlr(face_id, 0);
+ const double beta_mu_l = (1. / Dimension) * alpha_mu_l[face_id] * mes_l[face_id];
+ const double beta_lambda_l = (1. / Dimension) * 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 =
+ primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_cell));
+ const TinyVector<Dimension> nil = [&] {
+ if (is_face_reversed_for_cell_i) {
+ return -primal_nlr(face_id, 0);
+ } else {
+ return primal_nlr(face_id, 0);
+ }
+ }();
+ 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);
+ 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) += M(i, j);
+ }
}
- }();
- 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);
- 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) += M(i, j);
- }
+ } 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) -= M(i, j);
}
- } 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) -= M(i, j);
+ }
+ }
+ }
+ }
+ }
+
+ const auto& dual_Cjr = diamond_mesh_data.Cjr();
+ 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 = alpha_mu_l[face_id];
+ const double alpha_lambda_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 =
+ primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
+
+ 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];
+
+ if (not primal_node_is_on_boundary[node_id]) {
+ const TinyVector<Dimension> nil = [&] {
+ if (is_face_reversed_for_cell_i) {
+ return -primal_nlr(face_id, i_node);
+ } else {
+ return primal_nlr(face_id, i_node);
+ }
+ }();
+
+ 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) -=
+ w_rj(node_id, j_cell) * M(i, j);
+ }
+ }
}
}
}
@@ -473,13 +542,7 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
}
}
}
- // for (size_t i = 0; i < S.numberOfRows(); ++i) {
- // for (size_t j = 0; j < S.numberOfRows(); ++j) {
- // std::cout << S(i, j) << " ; ";
- // }
- // std::cout << "\n";
- // }
- // std::exit(0);
+
CRSMatrix A{S};
Vector<double> b{number_of_dof * Dimension};
for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
@@ -488,6 +551,80 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
}
}
+ const auto& primal_cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix();
+ // 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& node_list = bc.nodeList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ if (not node_tag[node_id]) {
+ node_tag[node_id] = true;
+
+ const auto& node_cells = primal_node_to_cell_matrix[node_id];
+ for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
+ const CellId cell_id = node_cells[i_cell];
+
+ const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id];
+
+ for (size_t i_cell_face = 0; i_cell_face < primal_cell_to_face.size(); ++i_cell_face) {
+ FaceId face_id = primal_cell_to_face_matrix[cell_id][i_cell_face];
+
+ const bool is_face_reversed_for_cell_i = [=] {
+ for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) {
+ FaceId primal_face_id = primal_cell_to_face[i_face];
+ if (primal_face_id == face_id) {
+ return primal_face_cell_is_reversed(cell_id, i_face);
+ }
+ }
+ Assert(false, "cannot find cell's face");
+ return false;
+ }();
+
+ const auto& primal_face_to_node = primal_face_to_node_matrix[face_id];
+
+ for (size_t i_face_node = 0; i_face_node < primal_face_to_node.size(); ++i_face_node) {
+ if (primal_face_to_node[i_face_node] == node_id) {
+ const TinyVector<Dimension> nil = [=] {
+ if (is_face_reversed_for_cell_i) {
+ return -primal_nlr(face_id, i_face_node);
+ } else {
+ return primal_nlr(face_id, i_face_node);
+ }
+ }();
+
+ 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_l[face_id] * (Clr, nil) * eye +
+ alpha_mu_l[face_id] * tensorProduct(Clr, nil) +
+ alpha_lambda_l[face_id] * tensorProduct(nil, Clr);
+ for (size_t i = 0; i < Dimension; ++i) {
+ b[(cell_dof_number[cell_id] * Dimension) + i] += (M * value_list[i_node])[i];
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+
Vector<double> U{number_of_dof * Dimension};
U = 0;
@@ -495,19 +632,58 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
Vector r = A * U - b;
+ // for (CellId j = 0; j < number_of_dof; ++j) {
+ // for (size_t l = 0; l < Dimension; ++l) {
+ // std::cout << "Ucalc : " << U[(cell_dof_number[j] * Dimension) + l] << "\n";
+ // }
+ // }
+
std::cout << "real residu = " << std::sqrt((r, r)) << '\n';
- CellValue<double> Speed{mesh->connectivity()};
+ CellValue<TinyVector<Dimension>> Speed{mesh->connectivity()};
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Speed[cell_id] = U[cell_dof_number[cell_id]]; });
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ Speed[cell_id][i] = U[(cell_dof_number[cell_id] * Dimension) + i];
+ }
+ }
+ Vector<double> Uexacte{number_of_dof * Dimension};
+ for (CellId j = 0; j < number_of_dof; ++j) {
+ for (size_t l = 0; l < Dimension; ++l) {
+ Uexacte[(cell_dof_number[j] * Dimension) + l] = Uj[j][l];
+ }
+ }
+ // std::cout << "A*Uex - b =\n";
+ // for (CellId i_id = 0; i_id < number_of_dof; ++i_id) {
+ // for (size_t j = 0; j < Dimension; ++j) {
+ // std::cout << (A * Uexacte - b)[(cell_dof_number[i_id] * Dimension) + j] << "\n";
+ // }
+ // }
+ // std::cout << "Vecteur AU - b: \n";
+ // for (CellId i = 0; i < number_of_dof; ++i) {
+ // // double temp = 0;
+ // for (size_t k = 0; k < Dimension; ++k) {
+ // double temp = 0;
+ // for (CellId j = 0; j < number_of_dof; ++j) {
+ // for (size_t l = 0; l < Dimension; ++l) {
+ // temp += S((cell_dof_number[i] * Dimension) + k, (cell_dof_number[j] * Dimension) + l) * Uj[j][l];
+ // }
+ // }
+ // std::cout << temp << "; " << b[(cell_dof_number[i] * Dimension) + k] << "=> "
+ // << temp - b[(cell_dof_number[i] * Dimension) + k] << "\n";
+ // }
+ // }
+
+ // std::cout << "Vecteur Uexacte : \n";
+ // for (CellId i = 0; i < number_of_dof; ++i) {
+ // std::cout << Uj[i] << "\n";
+ // }
VTKWriter vtk_writer("Speed_" + std::to_string(Dimension), 0.01);
- vtk_writer.write(mesh, {NamedItemValue{"U", Speed}}, 0,
+ vtk_writer.write(mesh, {NamedItemValue{"U", Speed}, NamedItemValue{"Uexacte", Uj}}, 0,
true); // forces last output
}
-
} else {
throw NotImplementedError("not done in 1d");
}
@@ -596,6 +772,7 @@ template ElasticityDiamondScheme<1>::ElasticityDiamondScheme(
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const FunctionSymbolId&);
template ElasticityDiamondScheme<2>::ElasticityDiamondScheme(
@@ -603,6 +780,7 @@ template ElasticityDiamondScheme<2>::ElasticityDiamondScheme(
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const FunctionSymbolId&);
template ElasticityDiamondScheme<3>::ElasticityDiamondScheme(
@@ -610,4 +788,5 @@ template ElasticityDiamondScheme<3>::ElasticityDiamondScheme(
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
const FunctionSymbolId&,
const FunctionSymbolId&,
+ const FunctionSymbolId&,
const FunctionSymbolId&);
diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.hpp b/src/language/algorithms/ElasticityDiamondAlgorithm.hpp
index 6eea3945a..733aa7d58 100644
--- a/src/language/algorithms/ElasticityDiamondAlgorithm.hpp
+++ b/src/language/algorithms/ElasticityDiamondAlgorithm.hpp
@@ -22,7 +22,8 @@ class ElasticityDiamondScheme
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& f_id,
+ const FunctionSymbolId& U_id);
};
#endif // ELASTICITY_DIAMOND_ALGORITHM_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 3ac4f789f..61c90914a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -264,24 +264,25 @@ SchemeModule::SchemeModule()
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 FunctionSymbolId&,
+ const FunctionSymbolId&)>>(
[](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) -> void {
+ const FunctionSymbolId& f_id, const FunctionSymbolId& U_id) -> void {
switch (p_mesh->dimension()) {
case 1: {
- ElasticityDiamondScheme<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id};
+ ElasticityDiamondScheme<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
break;
}
case 2: {
- ElasticityDiamondScheme<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id};
+ ElasticityDiamondScheme<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
break;
}
case 3: {
- ElasticityDiamondScheme<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id};
+ ElasticityDiamondScheme<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id};
break;
}
default: {
--
GitLab