Skip to content
Snippets Groups Projects
Commit a751f5c4 authored by PATELA Julie's avatar PATELA Julie
Browse files

Add mesh matrix for elasticity equation

parent a4daceec
No related branches found
No related tags found
No related merge requests found
...@@ -359,10 +359,13 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( ...@@ -359,10 +359,13 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
diamond_mesh_data diamond_mesh_data
.xj()); .xj());
CellValue<double> fj = CellValue<TinyVector<Dimension>> fj = InterpolateItemValue<TinyVector<Dimension>(
InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj()); 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(); const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
...@@ -384,7 +387,127 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( ...@@ -384,7 +387,127 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme(
}(); }();
const CellValue<const double> primal_Vj = mesh_data.Vj(); 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 { } else {
throw NotImplementedError("not done in 1d"); throw NotImplementedError("not done in 1d");
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment