From a751f5c42ccc21b12bbcdd7e03d79c2d12d1ad5c Mon Sep 17 00:00:00 2001
From: Julie PATELA <julie.patela.ocre@cea.fr>
Date: Tue, 18 Aug 2020 11:43:22 +0200
Subject: [PATCH] Add mesh matrix for elasticity equation
---
.../algorithms/ElasticityDiamondAlgorithm.cpp | 129 +++++++++++++++++-
1 file changed, 126 insertions(+), 3 deletions(-)
diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
index 6ec5f17b0..f51e7678e 100644
--- a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
+++ b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp
@@ -359,10 +359,13 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
diamond_mesh_data
.xj());
- CellValue<double> fj =
- InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
+ CellValue<TinyVector<Dimension>> fj = InterpolateItemValue<TinyVector<Dimension>(
+ TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
- VTKWriter vtk_writer("D_" + std::to_string(Dimension), 0.01);
+ VTKWriter vtk_writer2("f_" + std::to_string(Dimension), 0.01);
+
+ vtk_writer2.write(mesh, {NamedItemValue{"f", fj}}, 0,
+ true); // forces last output
const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
@@ -384,7 +387,127 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
}();
const CellValue<const double> primal_Vj = mesh_data.Vj();
+
+ 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_mes_l_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_mes_l_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;
+ }();
+
+ 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;
+ }();
+
+ TinyMatrix<Dimension, double> eye = zero;
+ for (size_t i = 0; i < Dimension; ++i) {
+ eye(i, i) = 1;
+ }
+
+ const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr();
+ const auto& primal_face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed();
+ 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);
+ }
+ }();
+ 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);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ // 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) {
+ for (size_t i = 0; i < Dimension; ++i) {
+ b[(cell_dof_number[cell_id] * Dimension) + i] = primal_Vj[cell_id] * fj[cell_id][i];
+ }
+ }
+
+ Vector<double> U{number_of_dof * Dimension};
+ U = 0;
+
+ BiCGStab{b, A, U, 1000, 1e-15};
+
+ Vector r = A * U - b;
+
+ std::cout << "real residu = " << std::sqrt((r, r)) << '\n';
+
+ CellValue<double> Speed{mesh->connectivity()};
+
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Speed[cell_id] = U[cell_dof_number[cell_id]]; });
+
+ VTKWriter vtk_writer("Speed_" + std::to_string(Dimension), 0.01);
+
+ vtk_writer.write(mesh, {NamedItemValue{"U", Speed}}, 0,
+ true); // forces last output
}
+
} else {
throw NotImplementedError("not done in 1d");
}
--
GitLab