From fc63eeaa576bf36c0da7ecae8b2ade778d0026eb Mon Sep 17 00:00:00 2001
From: Drouard <axelle.drouard@orange.fr>
Date: Tue, 7 Jun 2022 08:47:34 +0200
Subject: [PATCH] Fix renormalization v
---
src/scheme/ScalarNodalScheme.cpp | 29 +++++++++++++++++------------
1 file changed, 17 insertions(+), 12 deletions(-)
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index c834964f7..2837089cd 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -143,7 +143,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
- const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_bound,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_b,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k,
const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
@@ -409,7 +409,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
{
CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
- CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
+ CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_b->cellValues();
const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()};
@@ -568,8 +568,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const NodeValue<const TinyVector<Dimension>> v = [&] {
NodeValue<TinyVector<Dimension>> compute_v{mesh->connectivity()};
parallel_for(
- mesh->numberOfNodes(),
- PUGS_LAMBDA(NodeId node_id) { compute_v[node_id] = node_kapparb[node_id] * exterior_normal[node_id]; });
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (is_boundary_node[node_id]) {
+ compute_v[node_id] = 1. / l2Norm(node_kapparb[node_id] * exterior_normal[node_id]) *
+ node_kapparb[node_id] * exterior_normal[node_id];
+ }
+ });
return compute_v;
}();
@@ -704,7 +708,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
b[cell_id] += node_boundary_values[node_id] *
(1. / (sum_theta[node_id] * l2Norm(node_kapparb[node_id] * n)) *
- dot(kapparb_invBetar[node_id] * sum_Cjr[node_id], P * Cjr(cell_id, i_node)) -
+ dot(kapparb_invBetar[node_id] * sum_Cjr[node_id], P * Cjr(cell_id, i_node)) +
dot(Cjr(cell_id, i_node), n));
}
} else if (node_is_corner[node_id]) {
@@ -727,8 +731,9 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const CellId cell_id_k = node_to_cell[i_cell_k];
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
- b[cell_id_j] += dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
- Cjr(cell_id_j, i_node_j));
+ b[cell_id_j] +=
+ dot(node_boundary_values[node_id] * kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+ Cjr(cell_id_j, i_node_j));
}
}
@@ -742,7 +747,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
const size_t i_node_k = node_local_number_in_its_cells[i_cell_k];
b[cell_id_j] +=
- dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
+ dot(node_boundary_values[node_id] * corner_kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k),
Cjr(cell_id_j, i_node_j));
}
}
@@ -771,7 +776,7 @@ ScalarNodalSchemeHandler::solution() const
}
ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
- const std::shared_ptr<const IDiscreteFunction>& cell_k_bound,
+ const std::shared_ptr<const IDiscreteFunction>& cell_k_b,
const std::shared_ptr<const IDiscreteFunction>& cell_k,
const std::shared_ptr<const IDiscreteFunction>& f,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
@@ -781,7 +786,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
throw NormalError("primal discrete functions are not defined on the same mesh");
}
- checkDiscretizationType({cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
+ checkDiscretizationType({cell_k_b, cell_k, f}, DiscreteFunctionType::P0);
switch (i_mesh->dimension()) {
case 1: {
@@ -797,7 +802,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<2>>(mesh,
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
@@ -812,7 +817,7 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
m_scheme =
std::make_unique<ScalarNodalScheme<3>>(mesh,
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b),
std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
bc_descriptor_list);
--
GitLab