diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 10c92f181a0a53ce0e7216d595eb725dda843404..ecba28f40d2199858b79033aa5c4c9ee3fd6b394 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -491,6 +491,33 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("moleculardiffusion",
+ std::make_shared<BuiltinFunctionEmbedder<std::tuple<
+ std::shared_ptr<const IDiscreteFunction>,
+ 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>&,
+ const std::vector<std::shared_ptr<
+ 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::tuple<std::shared_ptr<const IDiscreteFunction>,
+ std::shared_ptr<const IDiscreteFunction>> {
+ return VectorDiamondSchemeHandler{alpha, lambdab, mub, lambda, mu,
+ f, bc_descriptor_list}
+ .apply();
+ }
+
+ ));
+
this->_addBuiltinFunction("heat2",
std::make_shared<BuiltinFunctionEmbedder<
void(std::shared_ptr<const IMesh>,
diff --git a/src/scheme/VectorDiamondScheme.cpp b/src/scheme/VectorDiamondScheme.cpp
index a337db06c1e5e781a1209befef7a4453c3507149..e8c392e1062c20e02b31426a12c50dee9652366e 100644
--- a/src/scheme/VectorDiamondScheme.cpp
+++ b/src/scheme/VectorDiamondScheme.cpp
@@ -158,6 +158,8 @@ class VectorDiamondSchemeHandler::IVectorDiamondScheme
public:
virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0;
virtual std::shared_ptr<const IDiscreteFunction> getDualSolution() const = 0;
+ virtual std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> apply()
+ const = 0;
IVectorDiamondScheme() = default;
virtual ~IVectorDiamondScheme() = default;
@@ -265,6 +267,12 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
return m_dual_solution;
}
+ std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+ apply() const final
+ {
+ return {m_solution, m_dual_solution};
+ }
+
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,
@@ -933,14 +941,16 @@ 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>>>& Velocity,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>& velocity,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>>& dual_velocity,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
Assert(mesh == alpha->mesh());
- Assert(mesh == Velocity->mesh());
+ Assert(mesh == velocity->mesh());
Assert(dual_lambda->mesh() == dual_mu->mesh());
Assert(dual_lambdab->mesh() == dual_mu->mesh());
Assert(dual_mub->mesh() == dual_mu->mesh());
+ Assert(dual_mub->mesh() == dual_velocity->mesh());
if (DiamondDualMeshManager::instance().getDiamondDualMesh(mesh) != dual_mu->mesh()) {
throw NormalError("diffusion coefficient is not defined on the dual mesh!");
}
@@ -1187,7 +1197,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
CellValue<const double> dual_mubj = dual_mub->cellValues();
CellValue<const double> dual_lambdabj = dual_lambdab->cellValues();
// attention, fj not in this context
- CellValue<const TinyVector<Dimension>> fj = Velocity->cellValues();
+ CellValue<const TinyVector<Dimension>> fj = velocity->cellValues();
const CellValue<const double> dual_Vj = diamond_mesh_data.Vj();
@@ -1365,7 +1375,7 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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);
- flux[face_id] = M * (Velocity[i_id] - Velocity[j_id]);
+ flux[face_id] = M * (velocity[i_id] - velocity[j_id]);
}
}
@@ -1408,7 +1418,11 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
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];
- flux[face_id] -= w_rj(node_id, j_cell) * M * (Velocity[i_id] - Velocity[j_id]);
+ flux[face_id] -= w_rj(node_id, j_cell) * M * (velocity[i_id] - velocity[j_id]);
+ if (primal_face_is_neumann[face_id]) {
+ flux[face_id] += 1.e0 * w_rj(node_id, j_cell) * Mb *
+ (velocity[i_id] - dual_velocity[face_dual_cell_id[face_id]]);
+ }
}
if (primal_node_is_on_boundary[node_id]) {
for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
@@ -1417,7 +1431,8 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
for (size_t i = 0; i < Dimension; ++i) {
for (size_t j = 0; j < Dimension; ++j) {
// Mb? + attention, bad index l_id
- flux[face_id] -= w_rl(node_id, l_face) * M(i, j) * (Velocity[i_id] - Velocity[l_id]);
+ flux[face_id] -= w_rl(node_id, l_face) * M(i, j) *
+ (velocity[i_id] - dual_velocity[face_dual_cell_id[l_id]]);
}
}
}
@@ -1490,6 +1505,12 @@ class VectorDiamondSchemeHandler::VectorDiamondScheme : public VectorDiamondSche
}
};
+std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>>
+VectorDiamondSchemeHandler::apply() const
+{
+ return m_scheme->apply();
+}
+
std::shared_ptr<const IDiscreteFunction>
VectorDiamondSchemeHandler::solution() const
{
diff --git a/src/scheme/VectorDiamondScheme.hpp b/src/scheme/VectorDiamondScheme.hpp
index 7208a091cccab1ef5f0098c8afa9a5164b8b4fd0..09df26fe4651466fe7ddd7908f989a42326ff363 100644
--- a/src/scheme/VectorDiamondScheme.hpp
+++ b/src/scheme/VectorDiamondScheme.hpp
@@ -39,6 +39,8 @@ class VectorDiamondSchemeHandler
std::shared_ptr<const IDiscreteFunction> dual_solution() const;
+ std::tuple<std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> apply() const;
+
VectorDiamondSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& alphab,
const std::shared_ptr<const IDiscreteFunction>& lambdab,