From dbf0459c8346039cc9fc75602a1aee09344e0cad Mon Sep 17 00:00:00 2001
From: labourasse <labourassee@gmail.com>
Date: Mon, 14 Feb 2022 09:39:20 +0100
Subject: [PATCH] Correction
---
src/scheme/VectorDiamondScheme.cpp | 71 ++++++++++++++++--------------
1 file changed, 39 insertions(+), 32 deletions(-)
diff --git a/src/scheme/VectorDiamondScheme.cpp b/src/scheme/VectorDiamondScheme.cpp
index 87173047e..fbbe3501e 100644
--- a/src/scheme/VectorDiamondScheme.cpp
+++ b/src/scheme/VectorDiamondScheme.cpp
@@ -281,11 +281,11 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mub,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_lambda,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mu,
- const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>& f,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>& source,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
Assert(mesh == alpha->mesh());
- Assert(mesh == f->mesh());
+ Assert(mesh == source->mesh());
Assert(dual_lambda->mesh() == dual_mu->mesh());
Assert(dual_lambdab->mesh() == dual_mu->mesh());
Assert(dual_mub->mesh() == dual_mu->mesh());
@@ -533,7 +533,10 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
CellValue<const double> dual_mubj = dual_mub->cellValues();
CellValue<const double> dual_lambdabj = dual_lambdab->cellValues();
- CellValue<const TinyVector<Dimension>> fj = f->cellValues();
+ CellValue<const TinyVector<Dimension>> fj = source->cellValues();
+ // for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ // std::cout << xj[cell_id] << "-> fj[" << cell_id << "]=" << fj[cell_id] << '\n';
+ // }
const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
@@ -727,7 +730,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
TinyMatrix<Dimension> Mb =
beta_mub_l * I + beta_mub_l * tensorProduct(nil, nil) + beta_lambdab_l * tensorProduct(nil, nil);
TinyMatrix<Dimension> N = 1.e0 * tensorProduct(nil, nil);
-
+ double coef_adim = beta_mu_l + beta_lambdab_l;
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];
if (i_cell == j_cell) {
@@ -742,9 +745,9 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
}
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.e0 : 0) + N(i, j);
+ ((i == j) ? -coef_adim : 0) + coef_adim * N(i, j);
S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + j) +=
- (i == j) ? 1.e0 : 0;
+ (i == j) ? coef_adim : 0;
}
}
}
@@ -847,7 +850,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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.e60;
+ S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + i) += 1.e0;
}
}
}
@@ -874,7 +877,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
const FaceId face_id = face_list[i_face];
for (size_t i = 0; i < Dimension; ++i) {
- b[(face_dof_number[face_id] * Dimension) + i] += 1.e60 * value_list[i_face][i];
+ b[(face_dof_number[face_id] * Dimension) + i] += 1.e0 * value_list[i_face][i];
}
}
}
@@ -1231,15 +1234,15 @@ class EnergyComputerHandler::EnergyComputer : public EnergyComputerHandler::IEne
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_lambdab,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mub,
// const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_lambda,
- // const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mu,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>& U,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>& dual_U,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>& source,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
// Assert(mesh == alpha->mesh());
Assert(mesh == U->mesh());
// Assert(dual_lambda->mesh() == dual_mu->mesh());
- // Assert(dual_lambdab->mesh() == dual_mu->mesh());
+ Assert(dual_lambdab->mesh() == dual_U->mesh());
Assert(U->mesh() == source->mesh());
Assert(dual_lambdab->mesh() == dual_mub->mesh());
if (DiamondDualMeshManager::instance().getDiamondDualMesh(mesh) != dual_mub->mesh()) {
@@ -1651,28 +1654,28 @@ class EnergyComputerHandler::EnergyComputer : public EnergyComputerHandler::IEne
return computed_dual_cell_face_id;
}();
- m_dual_solution = std::make_shared<DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>(diamond_mesh);
- m_solution = std::make_shared<DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>(mesh);
- const auto& solution = *U;
- auto& dual_solution = *m_dual_solution;
- dual_solution.fill(zero);
- const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
- for (CellId cell_id = 0; cell_id < diamond_mesh->numberOfCells(); ++cell_id) {
- const FaceId face_id = dual_cell_face_id[cell_id];
- CellId cell_id1 = face_to_cell_matrix[face_id][0];
- if (primal_face_is_on_boundary[face_id]) {
- for (size_t i = 0; i < Dimension; ++i) {
- // A revoir!!
- dual_solution[cell_id][i] = solution[cell_id1][i];
- }
- } else {
- CellId cell_id1 = face_to_cell_matrix[face_id][0];
- CellId cell_id2 = face_to_cell_matrix[face_id][1];
- for (size_t i = 0; i < Dimension; ++i) {
- dual_solution[cell_id][i] = 0.5 * (solution[cell_id1][i] + solution[cell_id2][i]);
- }
- }
- }
+ m_dual_solution = std::make_shared<DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>(diamond_mesh);
+ // m_solution = std::make_shared<DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>(mesh);
+ // const auto& solution = *U;
+ auto& dual_solution = *dual_U;
+ // dual_solution.fill(zero);
+ // const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
+ // for (CellId cell_id = 0; cell_id < diamond_mesh->numberOfCells(); ++cell_id) {
+ // const FaceId face_id = dual_cell_face_id[cell_id];
+ // CellId cell_id1 = face_to_cell_matrix[face_id][0];
+ // if (primal_face_is_on_boundary[face_id]) {
+ // for (size_t i = 0; i < Dimension; ++i) {
+ // // A revoir!!
+ // dual_solution[cell_id][i] = solution[cell_id1][i];
+ // }
+ // } else {
+ // CellId cell_id1 = face_to_cell_matrix[face_id][0];
+ // CellId cell_id2 = face_to_cell_matrix[face_id][1];
+ // for (size_t i = 0; i < Dimension; ++i) {
+ // dual_solution[cell_id][i] = 0.5 * (solution[cell_id1][i] + solution[cell_id2][i]);
+ // }
+ // }
+ // }
// const Array<int> non_zeros{number_of_dof * Dimension};
// non_zeros.fill(Dimension * Dimension);
@@ -1947,6 +1950,7 @@ EnergyComputerHandler::EnergyComputerHandler(
const std::shared_ptr<const IDiscreteFunction>& dual_lambdab,
const std::shared_ptr<const IDiscreteFunction>& dual_mub,
const std::shared_ptr<const IDiscreteFunction>& U,
+ const std::shared_ptr<const IDiscreteFunction>& dual_U,
const std::shared_ptr<const IDiscreteFunction>& source,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
@@ -1971,6 +1975,7 @@ EnergyComputerHandler::EnergyComputerHandler(
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambdab),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(U),
+ std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(dual_U),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(source),
bc_descriptor_list);
break;
@@ -1991,6 +1996,7 @@ EnergyComputerHandler::EnergyComputerHandler(
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambdab),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(U),
+ std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(dual_U),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(source),
bc_descriptor_list);
break;
@@ -2011,6 +2017,7 @@ EnergyComputerHandler::EnergyComputerHandler(
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambdab),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(U),
+ std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(dual_U),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(source),
bc_descriptor_list);
break;
--
GitLab