diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 232f5b18f45f8d9301dc8eb92422167bd6f4d158..fb41b96a941f948aaee30ae72006ac4518d57f09 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -404,23 +404,24 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("parabolicheat",
- std::make_shared<BuiltinFunctionEmbedder<
- std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
- const std::shared_ptr<const IDiscreteFunction>&,
- const std::shared_ptr<const IDiscreteFunction>&,
- const std::vector<std::shared_ptr<
- const IBoundaryConditionDescriptor>>&)>>(
-
- [](const std::shared_ptr<const IDiscreteFunction>& alpha,
- const std::shared_ptr<const IDiscreteFunction>& mu_dual,
- const std::shared_ptr<const IDiscreteFunction>& f,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
- bc_descriptor_list) -> const std::shared_ptr<const IDiscreteFunction> {
- return ScalarDiamondSchemeHandler{alpha, mu_dual, f, bc_descriptor_list}.solution();
- }
+ this->_addBuiltinFunction(
+ "parabolicheat",
+ std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+ const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&)>>(
+
+ [](const std::shared_ptr<const IDiscreteFunction>& alpha,
+ const std::shared_ptr<const IDiscreteFunction>& mub_dual,
+ const std::shared_ptr<const IDiscreteFunction>& mu_dual, const std::shared_ptr<const IDiscreteFunction>& f,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ -> const std::shared_ptr<const IDiscreteFunction> {
+ return ScalarDiamondSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
+ }
- ));
+ ));
this->_addBuiltinFunction("unsteadyelasticity",
std::make_shared<BuiltinFunctionEmbedder<
@@ -463,6 +464,8 @@ SchemeModule::SchemeModule()
this->_addBuiltinFunction("unsteadyelasticity",
std::make_shared<BuiltinFunctionEmbedder<
std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
+ const std::shared_ptr<const IDiscreteFunction>&,
const std::shared_ptr<const IDiscreteFunction>&,
const std::shared_ptr<const IDiscreteFunction>&,
const std::shared_ptr<const IDiscreteFunction>&,
@@ -470,12 +473,16 @@ SchemeModule::SchemeModule()
const IBoundaryConditionDescriptor>>&)>>(
[](const std::shared_ptr<const IDiscreteFunction> alpha,
+ const std::shared_ptr<const IDiscreteFunction> lambdab,
+ const std::shared_ptr<const IDiscreteFunction> mub,
const std::shared_ptr<const IDiscreteFunction> lambda,
const std::shared_ptr<const IDiscreteFunction> mu,
const std::shared_ptr<const IDiscreteFunction> f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list) -> const std::shared_ptr<const IDiscreteFunction> {
- return VectorDiamondSchemeHandler{alpha, lambda, mu, f, bc_descriptor_list}.solution();
+ return VectorDiamondSchemeHandler{alpha, lambdab, mub, lambda, mu,
+ f, bc_descriptor_list}
+ .solution();
}
));
diff --git a/src/scheme/ScalarDiamondScheme.cpp b/src/scheme/ScalarDiamondScheme.cpp
index 93538514fc93ac7a566bddc270228beee34479c0..bd9feb7d919a8605794e3c4146a2715f1063f6ee 100644
--- a/src/scheme/ScalarDiamondScheme.cpp
+++ b/src/scheme/ScalarDiamondScheme.cpp
@@ -276,6 +276,7 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
ScalarDiamondScheme(const std::shared_ptr<const MeshType>& mesh,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mub,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mu,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
@@ -284,6 +285,10 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
throw NormalError("diffusion coefficient is not defined on the dual mesh!");
}
+ if (DiamondDualMeshManager::instance().getDiamondDualMesh(mesh) != dual_mub->mesh()) {
+ throw NormalError("boundary diffusion coefficient is not defined on the dual mesh!");
+ }
+
constexpr ItemType FaceType = [] {
if constexpr (Dimension > 1) {
return ItemType::face;
@@ -491,7 +496,8 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
mesh->connectivity());
- CellValue<const double> dual_kappaj = dual_mu->cellValues();
+ CellValue<const double> dual_kappaj = dual_mu->cellValues();
+ CellValue<const double> dual_kappajb = dual_mub->cellValues();
const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
@@ -605,12 +611,26 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
return computed_alpha_l;
}();
+ FaceValue<const double> alphab_l = [&] {
+ CellValue<double> alpha_jb{diamond_mesh->connectivity()};
+
+ parallel_for(
+ diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId diamond_cell_id) {
+ alpha_jb[diamond_cell_id] = dual_kappajb[diamond_cell_id] / dual_Vj[diamond_cell_id];
+ });
+
+ FaceValue<double> computed_alpha_lb{mesh->connectivity()};
+ mapper->fromDualCell(alpha_jb, computed_alpha_lb);
+ return computed_alpha_lb;
+ }();
+
const CellValue<const double> primal_Vj = mesh_data.Vj();
SparseMatrixDescriptor S{number_of_dof};
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
const double beta_l = l2Norm(dualClj[face_id]) * alpha_l[face_id] * mes_l[face_id];
+ const double betab_l = l2Norm(dualClj[face_id]) * alphab_l[face_id] * mes_l[face_id];
for (size_t i_cell = 0; i_cell < primal_face_to_cell.size(); ++i_cell) {
const CellId cell_id1 = primal_face_to_cell[i_cell];
for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
@@ -618,7 +638,7 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
if (i_cell == j_cell) {
S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
if (primal_face_is_neumann[face_id]) {
- S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= beta_l;
+ S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= betab_l;
}
} else {
S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
@@ -636,7 +656,8 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- const double alpha_face_id = mes_l[face_id] * alpha_l[face_id];
+ const double alpha_face_id = mes_l[face_id] * alpha_l[face_id];
+ const double alphab_face_id = mes_l[face_id] * alphab_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];
@@ -660,13 +681,14 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
if (dual_node_primal_node_id[dual_node_id] == node_id) {
const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
- const double a_ir = alpha_face_id * (nil, Clr);
+ const double a_ir = alpha_face_id * (nil, Clr);
+ const double ab_ir = alphab_face_id * (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];
S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
if (primal_face_is_neumann[face_id]) {
- S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * a_ir;
+ S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * ab_ir;
}
}
if (primal_node_is_on_boundary[node_id]) {
@@ -675,7 +697,7 @@ class ScalarDiamondSchemeHandler::ScalarDiamondScheme : public ScalarDiamondSche
if (primal_face_is_on_boundary[l_id]) {
S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
if (primal_face_is_neumann[face_id]) {
- S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * a_ir;
+ S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * ab_ir;
}
}
}
@@ -758,12 +780,13 @@ ScalarDiamondSchemeHandler::solution() const
ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& alpha,
+ const std::shared_ptr<const IDiscreteFunction>& dual_mub,
const std::shared_ptr<const IDiscreteFunction>& dual_mu,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
- checkDiscretizationType({alpha, dual_mu, f}, DiscreteFunctionType::P0);
+ checkDiscretizationType({alpha, dual_mub, dual_mu, f}, DiscreteFunctionType::P0);
switch (i_mesh->dimension()) {
case 1: {
@@ -773,6 +796,7 @@ ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
m_scheme =
std::make_unique<ScalarDiamondScheme<1>>(std::dynamic_pointer_cast<const MeshType>(i_mesh),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
@@ -785,6 +809,7 @@ ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
m_scheme =
std::make_unique<ScalarDiamondScheme<2>>(std::dynamic_pointer_cast<const MeshType>(i_mesh),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
@@ -797,6 +822,7 @@ ScalarDiamondSchemeHandler::ScalarDiamondSchemeHandler(
m_scheme =
std::make_unique<ScalarDiamondScheme<3>>(std::dynamic_pointer_cast<const MeshType>(i_mesh),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
diff --git a/src/scheme/ScalarDiamondScheme.hpp b/src/scheme/ScalarDiamondScheme.hpp
index efdd2937b0c4127b170a96fd780a493db41b8178..ab5303915060ea25acf6891dc6f5bcb5bd853535 100644
--- a/src/scheme/ScalarDiamondScheme.hpp
+++ b/src/scheme/ScalarDiamondScheme.hpp
@@ -42,6 +42,7 @@ class ScalarDiamondSchemeHandler
ScalarDiamondSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& alpha,
+ const std::shared_ptr<const IDiscreteFunction>& mu_dualb,
const std::shared_ptr<const IDiscreteFunction>& mu_dual,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
diff --git a/src/scheme/VectorDiamondScheme.cpp b/src/scheme/VectorDiamondScheme.cpp
index 92a262cb7928b2c3b4fbfc6518aa8ff4ee83b326..0bb3acf84527cf05e8b45d595a5bce576ca14ebd 100644
--- a/src/scheme/VectorDiamondScheme.cpp
+++ b/src/scheme/VectorDiamondScheme.cpp
@@ -258,6 +258,8 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
VectorDiamondScheme(const std::shared_ptr<const MeshType>& mesh,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
+ 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>>>& f,
@@ -266,6 +268,8 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
Assert(mesh == alpha->mesh());
Assert(mesh == f->mesh());
Assert(dual_lambda->mesh() == dual_mu->mesh());
+ Assert(dual_lambdab->mesh() == dual_mu->mesh());
+ Assert(dual_mub->mesh() == dual_mu->mesh());
if (DiamondDualMeshManager::instance().getDiamondDualMesh(mesh) != dual_mu->mesh()) {
throw NormalError("diffusion coefficient is not defined on the dual mesh!");
}
@@ -527,8 +531,10 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper(
mesh->connectivity());
- CellValue<const double> dual_muj = dual_mu->cellValues();
- CellValue<const double> dual_lambdaj = dual_lambda->cellValues();
+ CellValue<const double> dual_muj = dual_mu->cellValues();
+ CellValue<const double> dual_lambdaj = dual_lambda->cellValues();
+ CellValue<const double> dual_mubj = dual_mub->cellValues();
+ CellValue<const double> dual_lambdabj = dual_lambdab->cellValues();
CellValue<const TinyVector<Dimension>> fj = f->cellValues();
@@ -658,12 +664,40 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
return computed_alpha_l;
}();
+ FaceValue<const double> alpha_lambdab_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_lambdabj[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_mub_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_mubj[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;
+ }();
+
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 double beta_mub_l = l2Norm(dualClj[face_id]) * alpha_mub_l[face_id] * mes_l[face_id];
+ const double beta_lambdab_l = l2Norm(dualClj[face_id]) * alpha_lambdab_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];
@@ -680,6 +714,8 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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> Mb =
+ beta_mub_l * eye + beta_mub_l * tensorProduct(nil, nil) + beta_lambdab_l * tensorProduct(nil, nil);
TinyMatrix<Dimension, double> N = tensorProduct(nil, nil);
if (i_cell == j_cell) {
@@ -687,7 +723,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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);
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);
+ S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) -= Mb(i, j);
}
if (primal_face_is_symmetry[face_id]) {
S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) +=
@@ -718,8 +754,10 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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];
+ 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];
+ const double alpha_mub_face_id = mes_l[face_id] * alpha_mub_l[face_id];
+ const double alpha_lambdab_face_id = mes_l[face_id] * alpha_lambdab_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];
@@ -746,6 +784,9 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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);
+ TinyMatrix<Dimension, double> Mb = alpha_mub_face_id * (Clr, nil) * eye +
+ alpha_mub_face_id * tensorProduct(Clr, nil) +
+ alpha_lambdab_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];
@@ -755,7 +796,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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);
+ w_rj(node_id, j_cell) * Mb(i, j);
}
}
}
@@ -774,7 +815,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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);
+ w_rl(node_id, l_face) * Mb(i, j);
}
}
}
@@ -901,6 +942,8 @@ VectorDiamondSchemeHandler::solution() const
VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& alpha,
+ const std::shared_ptr<const IDiscreteFunction>& dual_lambdab,
+ const std::shared_ptr<const IDiscreteFunction>& dual_mub,
const std::shared_ptr<const IDiscreteFunction>& dual_lambda,
const std::shared_ptr<const IDiscreteFunction>& dual_mu,
const std::shared_ptr<const IDiscreteFunction>& f,
@@ -908,7 +951,7 @@ VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
{
const std::shared_ptr i_mesh = getCommonMesh({alpha, f});
const std::shared_ptr i_dual_mesh = getCommonMesh({dual_lambda, dual_mu});
- checkDiscretizationType({alpha, dual_lambda, dual_mu, f}, DiscreteFunctionType::P0);
+ checkDiscretizationType({alpha, dual_lambdab, dual_mub, dual_lambda, dual_mu, f}, DiscreteFunctionType::P0);
switch (i_mesh->dimension()) {
case 1: {
@@ -924,6 +967,9 @@ VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
m_scheme =
std::make_unique<VectorDiamondScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(
+ dual_lambdab),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambda),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(f),
@@ -943,6 +989,9 @@ VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
m_scheme =
std::make_unique<VectorDiamondScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(
+ dual_lambdab),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambda),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(f),
@@ -962,6 +1011,9 @@ VectorDiamondSchemeHandler::VectorDiamondSchemeHandler(
m_scheme =
std::make_unique<VectorDiamondScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(
+ dual_lambdab),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mub),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_lambda),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu),
std::dynamic_pointer_cast<const DiscreteVectorFunctionType>(f),
diff --git a/src/scheme/VectorDiamondScheme.hpp b/src/scheme/VectorDiamondScheme.hpp
index cc2bca603b2f1be2d6af19b842f28b0dc896c418..17870b2fd3c81598e487c94137b7c9ce70de07e7 100644
--- a/src/scheme/VectorDiamondScheme.hpp
+++ b/src/scheme/VectorDiamondScheme.hpp
@@ -38,6 +38,8 @@ class VectorDiamondSchemeHandler
std::shared_ptr<const IDiscreteFunction> solution() const;
VectorDiamondSchemeHandler(
+ const std::shared_ptr<const IDiscreteFunction>& alphab,
+ const std::shared_ptr<const IDiscreteFunction>& lambdab,
const std::shared_ptr<const IDiscreteFunction>& alpha,
const std::shared_ptr<const IDiscreteFunction>& lambda,
const std::shared_ptr<const IDiscreteFunction>& mu,