diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 45db502f0eb5b3940c4082ba686fc16cf8c5cab4..edbce1e3f20264e1a5da949314b1f2a4123c133e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -42,6 +42,7 @@
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/ScalarDiamondScheme.hpp>
+#include <scheme/ScalarNodalScheme.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/VectorDiamondScheme.hpp>
#include <utils/Socket.hpp>
@@ -601,6 +602,24 @@ SchemeModule::SchemeModule()
return ScalarDiamondSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
}
+ ));
+ this->_addBuiltinFunction(
+ "nodalheat",
+ 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 ScalarNodalSchemeHandler{alpha, mub_dual, mu_dual, f, bc_descriptor_list}.solution();
+ }
+
));
this->_addBuiltinFunction("unsteadyelasticity",
diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp
index 791e1c8ff7b94c39663bea036faff5ecc5a08b51..215cf8510ec40f3cf19268d70c561f1034fc6ac6 100644
--- a/src/scheme/ScalarNodalScheme.cpp
+++ b/src/scheme/ScalarNodalScheme.cpp
@@ -3,7 +3,6 @@
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
-
class ScalarNodalSchemeHandler::IScalarNodalScheme
{
public:
@@ -56,6 +55,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
private:
const Array<const double> m_value_list;
const Array<const FaceId> m_face_list;
+ const Array<const NodeId> m_node_list;
public:
const Array<const FaceId>&
@@ -64,14 +64,22 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
return m_face_list;
}
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_node_list;
+ }
+
const Array<const double>&
valueList() const
{
return m_value_list;
}
- NeumannBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
- : m_value_list{value_list}, m_face_list{face_list}
+ NeumannBoundaryCondition(const Array<const FaceId>& face_list,
+ const Array<const NodeId>& node_list,
+ const Array<const double>& value_list)
+ : m_value_list{value_list}, m_face_list{face_list}, m_node_list{node_list}
{
Assert(m_value_list.size() == m_face_list.size());
}
@@ -118,7 +126,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
~FourierBoundaryCondition() = default;
};
-
public:
std::shared_ptr<const IDiscreteFunction>
getSolution() const final
@@ -127,15 +134,14 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh,
- const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
- const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_bound,
- 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)
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& alpha,
+ const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_bound,
+ 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)
{
-
- using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition,
- NeumannBoundaryCondition>;
+ using BoundaryCondition =
+ std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
@@ -179,6 +185,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
if constexpr (Dimension > 1) {
MeshFaceBoundary<Dimension> mesh_face_boundary =
getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
+ MeshNodeBoundary<Dimension> mesh_node_boundary =
+ getMeshNodeBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor());
const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId();
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
@@ -189,10 +197,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
mesh_face_boundary
.faceList());
- boundary_condition_list.push_back(NeumannBoundaryCondition{mesh_face_boundary.faceList(), value_list});
+ boundary_condition_list.push_back(
+ NeumannBoundaryCondition{mesh_face_boundary.faceList(), mesh_node_boundary.nodeList(), value_list});
} else {
- throw NotImplementedError("Dirichlet BC in 1d");
+ throw NotImplementedError("Neumann BC in 1d");
}
break;
}
@@ -207,268 +216,377 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand
}
}
- if constexpr (Dimension > 1) {
- const CellValue<const size_t> cell_dof_number = [&] {
- CellValue<size_t> compute_cell_dof_number{mesh->connectivity()};
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { compute_cell_dof_number[cell_id] = cell_id; });
- return compute_cell_dof_number;
- }();
- size_t number_of_dof = mesh->numberOfCells();
-
- const FaceValue<const size_t> face_dof_number = [&] {
- FaceValue<size_t> compute_face_dof_number{mesh->connectivity()};
- compute_face_dof_number.fill(std::numeric_limits<size_t>::max());
- for (const auto& boundary_condition : boundary_condition_list) {
- std::visit(
- [&](auto&& bc) {
- using T = std::decay_t<decltype(bc)>;
- if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
- (std::is_same_v<T, FourierBoundaryCondition>) or
- (std::is_same_v<T, DirichletBoundaryCondition>)) {
- const auto& face_list = bc.faceList();
-
- for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
- const FaceId face_id = face_list[i_face];
- if (compute_face_dof_number[face_id] != std::numeric_limits<size_t>::max()) {
- std::ostringstream os;
- os << "The face " << face_id << " is used at least twice for boundary conditions";
- throw NormalError(os.str());
+ MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+ const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr();
+
+ const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+ const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
+
+ const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
+
+ const auto is_boundary_node = mesh->connectivity().isBoundaryNode();
+
+ // const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
+ const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+ const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh->connectivity().nodeLocalNumbersInTheirCells();
+ const CellValue<const double> Vj = mesh_data.Vj();
+ const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+
+ const NodeValue<const bool> node_is_neumann = [&] {
+ NodeValue<bool> compute_node_is_neumann{mesh->connectivity()};
+ compute_node_is_neumann.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ const NodeId node_id = face_nodes[i_node];
+
+ compute_node_is_neumann[node_id] = true;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_node_is_neumann;
+ }();
+
+ const NodeValue<const bool> node_is_dirichlet = [&] {
+ NodeValue<bool> compute_node_is_dirichlet{mesh->connectivity()};
+ compute_node_is_dirichlet.fill(false);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
+ const auto& face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ const NodeId node_id = face_nodes[i_node];
+
+ compute_node_is_dirichlet[node_id] = true;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_node_is_dirichlet;
+ }();
+
+ const NodeValue<const bool> node_is_corner = [&] {
+ NodeValue<bool> compute_node_is_corner{mesh->connectivity()};
+ compute_node_is_corner.fill(false);
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ if (is_boundary_node[node_id]) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ const NodeId prev_node_id = cell_to_node_matrix[cell_id][(node_id - 1) % cell_nodes.size()];
+ const NodeId next_node_id = cell_to_node_matrix[cell_id][(node_id + 1) % cell_nodes.size()];
+
+ if (is_boundary_node[prev_node_id] and is_boundary_node[next_node_id]) {
+ compute_node_is_corner[node_id] = true;
+ }
+ }
+ }
+ });
+
+ return compute_node_is_corner;
+ }();
+
+ const NodeValue<const TinyVector<Dimension>> exterior_normal = [&] {
+ NodeValue<TinyVector<Dimension>> compute_exterior_normal;
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ compute_exterior_normal[node_id] = zero;
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const unsigned int i_node_loc = node_local_number_in_its_cell[cell_id];
+ compute_exterior_normal[node_id] += Cjr(cell_id, i_node_loc);
+ }
+ compute_exterior_normal[node_id] =
+ 1. / l2Norm(compute_exterior_normal[node_id]) * compute_exterior_normal[node_id];
+ });
+ return compute_exterior_normal;
+ }();
+
+ const FaceValue<const double> mes_l = [&] {
+ if constexpr (Dimension == 1) {
+ FaceValue<double> compute_mes_l{mesh->connectivity()};
+ compute_mes_l.fill(1);
+ return compute_mes_l;
+ } else {
+ return mesh_data.ll();
+ }
+ }();
+
+ const NodeValue<const double> node_boundary_values = [&] {
+ NodeValue<double> compute_node_boundary_values;
+ compute_node_boundary_values.fill(0);
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ const NodeId node_id = face_nodes[i_node];
+ if (node_is_dirichlet[node_id] == false) {
+ compute_node_boundary_values[node_id] += 0.5 * value_list[i_face] * mes_l[face_id];
+ }
+ }
+ }
+ } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ const NodeId node_id = face_nodes[i_node];
+ if (node_is_neumann[node_id] == false) {
+ compute_node_boundary_values[node_id] += 0.5 * value_list[i_face] * mes_l[face_id];
} else {
- compute_face_dof_number[face_id] = number_of_dof++;
+ compute_node_boundary_values[node_id] = value_list[i_face] * mes_l[face_id];
}
}
}
- },
- boundary_condition);
- }
+ }
+ },
+ boundary_condition);
+ }
+ return compute_node_boundary_values;
+ }();
- return compute_face_dof_number;
+ {
+ CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
+ CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
+
+ const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
+ NodeValue<TinyMatrix<Dimension>> kappa;
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ kappa[node_id] = zero;
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ double weight = 0;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ double local_weight = 1. / l2Norm(xr[node_id] - xj[J]);
+ kappa[node_id] += local_weight * cell_kappaj[J];
+ weight += local_weight;
+ }
+ kappa[node_id] = 1. / weight * kappa[node_id];
+ });
+ return kappa;
}();
- const FaceValue<const bool> face_is_neumann = [&] {
- FaceValue<bool> face_is_neumann{mesh->connectivity()};
- face_is_neumann.fill(false);
- for (const auto& boundary_condition : boundary_condition_list) {
- std::visit(
- [&](auto&& bc) {
- using T = std::decay_t<decltype(bc)>;
- if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
- (std::is_same_v<T, FourierBoundaryCondition>)) {
- const auto& face_list = bc.faceList();
-
- for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
- const FaceId face_id = face_list[i_face];
- face_is_neumann[face_id] = true;
- }
- }
- },
- boundary_condition);
- }
+ const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] {
+ NodeValue<TinyMatrix<Dimension>> kappa;
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ kappa[node_id] = zero;
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ double weight = 0;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ double local_weight = 1. / l2Norm(xr[node_id] - xj[J]);
+ kappa[node_id] += local_weight * cell_kappajb[J];
+ weight += local_weight;
+ }
+ kappa[node_id] = 1. / weight * kappa[node_id];
+ });
+ return kappa;
+ }();
- return face_is_neumann;
+ const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
+ NodeValue<TinyMatrix<Dimension>> beta;
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ beta[node_id] = zero;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cell[j];
+ beta[node_id] += tensorProduct(Cjr(J, R), xr[node_id] - xj[J]);
+ }
+ });
+ return beta;
}();
+ const NodeValue<const TinyMatrix<Dimension>> kappar_invBetar = [&] {
+ NodeValue<TinyMatrix<Dimension>> kappa_invBeta;
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ kappa_invBeta[node_id] = node_kappar[node_id] * inverse(node_betar[node_id]);
+ });
+ return kappa_invBeta;
+ }();
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+ const Array<int> non_zeros{mesh->numberOfCells()};
+ non_zeros.fill(Dimension); // Modif pour optimiser
+ CRSMatrixDescriptor<double> S(mesh->numberOfCells(), mesh->numberOfCells(), non_zeros);
- const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
- const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj();
- const NodeValue<const TinyVector<Dimension>>& xr = mesh_data.xr();
+ for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
- const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix();
- const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+ if (not is_boundary_node[node_id]) {
+ for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
+ const CellId J = node_to_cell[cell_id_j];
+ const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
+ const size_t node_id_j = node_local_number_in_j_cell[node_id];
- const NodeValuePerCell<const TinyVector<Dimension>>& Cjr = mesh_data.Cjr();
+ for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
+ const CellId K = node_to_cell[cell_id_k];
+ const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
+ const size_t node_id_k = node_local_number_in_k_cell[node_id];
- {
+ S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), Cjr(J, node_id_j));
+ }
+ }
+ } else if (not node_is_corner[node_id]) {
+ }
+ }
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+ const auto& node_list = bc.nodeList();
- CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues();
- CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_bound->cellValues();
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
- const CellValue<const double> Vj = mesh_data.Vj();
- const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+ for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
+ const CellId J = node_to_cell[cell_id_j];
+ const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
+ const size_t node_id_j = node_local_number_in_j_cell[node_id];
- const FaceValue<const double> mes_l = [&] {
- if constexpr (Dimension == 1) {
- FaceValue<double> compute_mes_l{mesh->connectivity()};
- compute_mes_l.fill(1);
- return compute_mes_l;
- } else {
- return mesh_data.ll();
- }
- }();
+ for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
+ const CellId K = node_to_cell[cell_id_k];
+ const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
+ const size_t node_id_k = node_local_number_in_k_cell[node_id];
- const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] {
- NodeValue<const TinyMatrix<Dimension>> kappa;
- parallel_for(
- mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- kappa[node_id] = zero;
- const auto& node_to_cell = node_to_cell_matrix[node_id];
- double weight = 0;
- for (size_t j = 0; j < node_to_cell.size(); ++j) {
- const CellId J = node_to_cell[j];
- double local_weight = 1 / l2Norm(xr[node_id]-xj[J]);
- kappa[node_id] += cell_kappaj[J] * local_weight;
- weight += local_weight;
+ const TinyMatrix<Dimension> I = identity;
+ const TinyVector<Dimension> n = exterior_normal[node_id];
+ const TinyMatrix<Dimension> nxn = tensorProduct(n, n);
+ const TinyMatrix<Dimension> P = I - nxn;
+
+ S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), P * Cjr(J, node_id_j));
+ }
}
- kappa[node_id] /= weight;
- }
- );
- return kappa;
- }();
-
- const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] {
- NodeValue<const TinyMatrix<Dimension>> kappa;
- parallel_for(
- mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- kappa[node_id] = zero;
+ }
+ } else if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
+ const auto& node_list = bc.faceList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
const auto& node_to_cell = node_to_cell_matrix[node_id];
- double weight = 0;
- for (size_t j = 0; j < node_to_cell.size(); ++j) {
- const CellId J = node_to_cell[j];
- double local_weight = 1 / l2Norm(xr[node_id]-xj[J]);
- kappa[node_id] += cell_kappajb[J] * local_weight;
- weight += local_weight;
- }
- kappa[node_id] /= weight;
- });
- return kappa;
- }();
-
- const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] {
- NodeValue<const TinyMatrix<Dimension>> beta;
- parallel_for(
- mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
- const auto& node_to_cell = node_to_cell_matrix[node_id];
- beta[node_id] = zero;
- for (size_t j = 0; j < node_to_cell.size(); ++j) {
- const CellId J = node_to_cell[j];
- const unsigned int R = node_local_number_in_its_cell[j];
- beta[node_id] += tensorProduct(Cjr(J,R),xr[node_id]-xj[J]);
- }
- });
- return beta;
- }();
-
-
- FaceValue<const double> alpha_l = [&] {
- FaceValue<double> alpha_j{mesh->connectivity()};
-
- parallel_for(
- mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
- alpha_j[face_id] = 1;
- });
-
- return alpha_j;
- }();
-
- FaceValue<const double> alphab_l = [&] {
- FaceValue<double> alpha_lb{mesh->connectivity()};
-
- parallel_for(
- mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
- alpha_lb[face_id] = 1; //Refaire
- });
-
- return alpha_lb;
- }();
-
-
- const Array<int> non_zeros{number_of_dof};
- non_zeros.fill(Dimension); //Modif pour optimiser
- CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zeros);
-
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- const auto& 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 < face_to_cell.size(); ++i_cell) {
- const CellId cell_id1 = face_to_cell[i_cell];
- for (size_t j_cell = 0; j_cell < face_to_cell.size(); ++j_cell) {
- const CellId cell_id2 = face_to_cell[j_cell];
- if (i_cell == j_cell) {
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
- if (face_is_neumann[face_id]) {
- S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= betab_l;
+
+ for (size_t cell_id_j = 0; cell_id_j < node_to_cell.size(); ++cell_id_j) {
+ const CellId J = node_to_cell[cell_id_j];
+ const auto& node_local_number_in_j_cell = node_local_numbers_in_their_cells.itemArray(J);
+ const size_t node_id_j = node_local_number_in_j_cell[node_id];
+
+ for (size_t cell_id_k = 0; cell_id_k < node_to_cell.size(); ++cell_id_k) {
+ const CellId K = node_to_cell[cell_id_k];
+ const auto& node_local_number_in_k_cell = node_local_numbers_in_their_cells.itemArray(K);
+ const size_t node_id_k = node_local_number_in_k_cell[node_id];
+
+ S(J, K) += dot(kappar_invBetar[node_id] * Cjr(K, node_id_k), Cjr(J, node_id_j));
+ }
}
- } else {
- S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
}
}
- }
- }
+ },
+ boundary_condition);
+ }
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- const size_t j = cell_dof_number[cell_id];
- S(j, j) += (*alpha)[cell_id] * Vj[cell_id];
- }
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ S(cell_id, cell_id) += Vj[cell_id];
+ }
+ CellValue<const double> fj = f->cellValues();
- for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
- if (face_is_dirichlet[face_id]) {
- S(face_dof_number[face_id], face_dof_number[face_id]) += 1;
- }
- }
-
- CellValue<const double> fj = f->cellValues();
+ CRSMatrix A{S.getCRSMatrix()};
+ Vector<double> b{mesh->numberOfCells()};
+ b = zero;
+ for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+ b[cell_id] = fj[cell_id] * Vj[cell_id];
+ }
- CRSMatrix A{S.getCRSMatrix()};
- Vector<double> b{number_of_dof};
- b = zero;
- for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
- b[cell_dof_number[cell_id]] = fj[cell_id] * Vj[cell_id];
- }
- // Dirichlet on b^L_D
- {
- for (const auto& boundary_condition : boundary_condition_list) {
- std::visit(
- [&](auto&& bc) {
- using T = std::decay_t<decltype(bc)>;
- if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) {
- const auto& face_list = bc.faceList();
- const auto& value_list = bc.valueList();
- for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
- const FaceId face_id = face_list[i_face];
- b[face_dof_number[face_id]] += value_list[i_face];
- }
- }
- },
- boundary_condition);
- }
- }
- // EL b^L
- for (const auto& boundary_condition : boundary_condition_list) {
- std::visit(
- [&](auto&& bc) {
- using T = std::decay_t<decltype(bc)>;
- if constexpr ((std::is_same_v<T, NeumannBoundaryCondition>) or
- (std::is_same_v<T, FourierBoundaryCondition>)) {
- const auto& face_list = bc.faceList();
- const auto& value_list = bc.valueList();
- for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
- FaceId face_id = face_list[i_face];
- b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face];
+ for (const auto& boundary_condition : boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) { // To do
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ b[face_id] += value_list[i_face];
+ }
+ } else if constexpr (std::is_same_v<T, NeumannBoundaryCondition>) {
+ const auto& node_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ NodeId node_id = node_list[i_node];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const TinyVector<Dimension> n = exterior_normal[node_id];
+
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ CellId J = node_to_cell[i_cell];
+ const auto& node_local_number_in_cell = node_local_numbers_in_their_cells.itemArray(J);
+ const size_t node_id_j = node_local_number_in_cell[node_id];
+ b[J] -= dot(Cjr(J, node_id_j), n) * value_list[i_node];
}
}
- },
- boundary_condition);
- }
+ }
+ },
+ boundary_condition);
+ }
- Vector<double> T{number_of_dof};
- T = zero;
+ Vector<double> T{mesh->numberOfCells()};
+ T = zero;
- LinearSolver solver;
- solver.solveLocalSystem(A, T, b);
+ LinearSolver solver;
+ solver.solveLocalSystem(A, T, b);
- m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
- auto& solution = *m_solution;
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_dof_number[cell_id]]; });
- }
+ m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh);
+ auto& solution = *m_solution;
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { solution[cell_id] = T[cell_id]; });
}
}
};
@@ -481,7 +599,7 @@ ScalarNodalSchemeHandler::solution() const
ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
const std::shared_ptr<const IDiscreteFunction>& alpha,
- const std::shared_ptr<const IDiscreteFunction>& k_bound,
+ const std::shared_ptr<const IDiscreteFunction>& cell_k_bound,
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)
@@ -490,30 +608,12 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
if (not i_mesh) {
throw NormalError("primal discrete functions are not defined on the same mesh");
}
- const std::shared_ptr i_dual_mesh = getCommonMesh({cell_k_bound, cell_k});
- if (not i_dual_mesh) {
- throw NormalError("dual discrete functions are not defined on the same mesh");
- }
+
checkDiscretizationType({alpha, cell_k_bound, cell_k, f}, DiscreteFunctionType::P0);
switch (i_mesh->dimension()) {
case 1: {
- using MeshType = Mesh<Connectivity<1>>;
- using DiscreteTensorFunctionType = DiscreteFunctionP0<1, TinyMatrix<1>>;
- using DiscreteScalarFunctionType = DiscreteFunctionP0<1, double>;
-
- std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
- if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
- throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
- }
-
- m_scheme =
- std::make_unique<ScalarNodalScheme<1>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
- std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
- bc_descriptor_list);
+ throw NormalError("The scheme is not implemented in dimension 1");
break;
}
case 2: {
@@ -523,16 +623,12 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
- if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
- throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
- }
-
m_scheme =
std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
- std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
- bc_descriptor_list);
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
+ bc_descriptor_list);
break;
}
case 3: {
@@ -542,16 +638,12 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler(
std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
- if (DualMeshManager::instance().getDiamondDualMesh(*mesh) != i_dual_mesh) {
- throw NormalError("dual variables are is not defined on the diamond dual of the primal mesh");
- }
-
m_scheme =
std::make_unique<ScalarNodalScheme<3>>(mesh, std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(alpha),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
- std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
- std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
- bc_descriptor_list);
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_bound),
+ std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k),
+ std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f),
+ bc_descriptor_list);
break;
}
default: {