From 689b7e2d5fdc83c5f14b75cb243d4d1d77a3cbfc Mon Sep 17 00:00:00 2001 From: Drouard <axelle.drouard@orange.fr> Date: Wed, 20 Apr 2022 15:47:23 +0200 Subject: [PATCH] Nouveaux fichiers pour solveur nodal --- src/scheme/ScalarNodalScheme.cpp | 851 +++++++++++++++++++++++++++++++ src/scheme/ScalarNodalScheme.hpp | 688 +++++++++++++++++++++++++ 2 files changed, 1539 insertions(+) create mode 100644 src/scheme/ScalarNodalScheme.cpp create mode 100644 src/scheme/ScalarNodalScheme.hpp diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp new file mode 100644 index 000000000..fbbaf8952 --- /dev/null +++ b/src/scheme/ScalarNodalScheme.cpp @@ -0,0 +1,851 @@ +#include <scheme/ScalarNodalScheme.hpp> + +#include <scheme/DiscreteFunctionP0.hpp> +#include <scheme/DiscreteFunctionUtils.hpp> + +template <size_t Dimension> +class ScalarNodalSchemeHandler::InterpolationWeightsManager +{ + private: + std::shared_ptr<const Mesh<Connectivity<Dimension>>> m_mesh; + FaceValue<bool> m_primal_face_is_on_boundary; + NodeValue<bool> m_primal_node_is_on_boundary; + CellValuePerNode<double> m_w_rj; + FaceValuePerNode<double> m_w_rl; + + public: + CellValuePerNode<double>& + wrj() + { + return m_w_rj; + } + + FaceValuePerNode<double>& + wrl() + { + return m_w_rl; + } + + void + compute() + { + using MeshDataType = MeshData<Dimension>; + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh); + + const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr(); + + const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl(); + const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj(); + const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix(); + const auto& node_to_face_matrix = m_mesh->connectivity().nodeToFaceMatrix(); + CellValuePerNode<double> w_rj{m_mesh->connectivity()}; + FaceValuePerNode<double> w_rl{m_mesh->connectivity()}; + + for (size_t i = 0; i < w_rl.numberOfValues(); ++i) { + w_rl[i] = std::numeric_limits<double>::signaling_NaN(); + } + + for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) { + SmallVector<double> b{Dimension + 1}; + b[0] = 1; + for (size_t i = 1; i < Dimension + 1; i++) { + b[i] = xr[i_node][i - 1]; + } + const auto& node_to_cell = node_to_cell_matrix[i_node]; + + if (not m_primal_node_is_on_boundary[i_node]) { + SmallMatrix<double> A{Dimension + 1, node_to_cell.size()}; + for (size_t j = 0; j < node_to_cell.size(); j++) { + A(0, j) = 1; + } + for (size_t i = 1; i < Dimension + 1; i++) { + for (size_t j = 0; j < node_to_cell.size(); j++) { + const CellId J = node_to_cell[j]; + A(i, j) = xj[J][i - 1]; + } + } + SmallVector<double> x{node_to_cell.size()}; + x = zero; + + LeastSquareSolver ls_solver; + ls_solver.solveLocalSystem(A, x, b); + + for (size_t j = 0; j < node_to_cell.size(); j++) { + w_rj(i_node, j) = x[j]; + } + } else { + int nb_face_used = 0; + for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) { + FaceId face_id = node_to_face_matrix[i_node][i_face]; + if (m_primal_face_is_on_boundary[face_id]) { + nb_face_used++; + } + } + SmallMatrix<double> A{Dimension + 1, node_to_cell.size() + nb_face_used}; + for (size_t j = 0; j < node_to_cell.size() + nb_face_used; j++) { + A(0, j) = 1; + } + for (size_t i = 1; i < Dimension + 1; i++) { + for (size_t j = 0; j < node_to_cell.size(); j++) { + const CellId J = node_to_cell[j]; + A(i, j) = xj[J][i - 1]; + } + } + for (size_t i = 1; i < Dimension + 1; i++) { + int cpt_face = 0; + for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) { + FaceId face_id = node_to_face_matrix[i_node][i_face]; + if (m_primal_face_is_on_boundary[face_id]) { + A(i, node_to_cell.size() + cpt_face) = xl[face_id][i - 1]; + cpt_face++; + } + } + } + + SmallVector<double> x{node_to_cell.size() + nb_face_used}; + x = zero; + + LeastSquareSolver ls_solver; + ls_solver.solveLocalSystem(A, x, b); + + for (size_t j = 0; j < node_to_cell.size(); j++) { + w_rj(i_node, j) = x[j]; + } + int cpt_face = node_to_cell.size(); + for (size_t i_face = 0; i_face < node_to_face_matrix[i_node].size(); ++i_face) { + FaceId face_id = node_to_face_matrix[i_node][i_face]; + if (m_primal_face_is_on_boundary[face_id]) { + w_rl(i_node, i_face) = x[cpt_face++]; + } + } + } + } + m_w_rj = w_rj; + m_w_rl = w_rl; + } + + InterpolationWeightsManager(std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh, + FaceValue<bool> primal_face_is_on_boundary, + NodeValue<bool> primal_node_is_on_boundary) + : m_mesh(mesh), + m_primal_face_is_on_boundary(primal_face_is_on_boundary), + m_primal_node_is_on_boundary(primal_node_is_on_boundary) + {} + + ~InterpolationWeightsManager() = default; +}; + +class ScalarNodalSchemeHandler::IScalarNodalScheme +{ + public: + virtual std::shared_ptr<const IDiscreteFunction> getSolution() const = 0; + + IScalarNodalScheme() = default; + virtual ~IScalarNodalScheme() = default; +}; + +template <size_t Dimension> +class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHandler::IScalarNodalScheme +{ + private: + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + using MeshDataType = MeshData<Dimension>; + + std::shared_ptr<const DiscreteFunctionP0<Dimension, double>> m_solution; + + class DirichletBoundaryCondition + { + private: + const Array<const double> m_value_list; + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + + const Array<const double>& + valueList() const + { + return m_value_list; + } + + DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list) + : m_value_list{value_list}, m_face_list{face_list} + { + Assert(m_value_list.size() == m_face_list.size()); + } + + ~DirichletBoundaryCondition() = default; + }; + + class NeumannBoundaryCondition + { + private: + const Array<const double> m_value_list; + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_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} + { + Assert(m_value_list.size() == m_face_list.size()); + } + + ~NeumannBoundaryCondition() = default; + }; + + class FourierBoundaryCondition + { + private: + const Array<const double> m_coef_list; + const Array<const double> m_value_list; + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + + const Array<const double>& + valueList() const + { + return m_value_list; + } + + const Array<const double>& + coefList() const + { + return m_coef_list; + } + + public: + FourierBoundaryCondition(const Array<const FaceId>& face_list, + const Array<const double>& coef_list, + const Array<const double>& value_list) + : m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list} + { + Assert(m_coef_list.size() == m_face_list.size()); + Assert(m_value_list.size() == m_face_list.size()); + } + + ~FourierBoundaryCondition() = default; + }; + + class SymmetryBoundaryCondition + { + private: + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + + public: + SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {} + + ~SymmetryBoundaryCondition() = default; + }; + + public: + std::shared_ptr<const IDiscreteFunction> + getSolution() const final + { + return m_solution; + } + + ScalarNodalScheme(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) + { + Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mu->mesh(), + "diffusion coefficient is not defined on the dual mesh!"); + Assert(DualMeshManager::instance().getDiamondDualMesh(*mesh) == dual_mub->mesh(), + "boundary diffusion coefficient is not defined on the dual mesh!"); + + using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, + NeumannBoundaryCondition, SymmetryBoundaryCondition>; + + using BoundaryConditionList = std::vector<BoundaryCondition>; + + BoundaryConditionList boundary_condition_list; + + for (const auto& bc_descriptor : bc_descriptor_list) { + bool is_valid_boundary_condition = true; + + switch (bc_descriptor->type()) { + case IBoundaryConditionDescriptor::Type::symmetry: { + throw NotImplementedError("NIY"); + break; + } + case IBoundaryConditionDescriptor::Type::dirichlet: { + const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor = + dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor); + if constexpr (Dimension > 1) { + MeshFaceBoundary<Dimension> mesh_face_boundary = + getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor()); + + const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId(); + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + Array<const double> value_list = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, + mesh_data.xl(), + mesh_face_boundary + .faceList()); + + boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list}); + + } else { + throw NotImplementedError("Dirichlet BC in 1d"); + } + break; + } + case IBoundaryConditionDescriptor::Type::fourier: { + throw NotImplementedError("NIY"); + break; + } + case IBoundaryConditionDescriptor::Type::neumann: { + const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor = + dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor); + + if constexpr (Dimension > 1) { + MeshFaceBoundary<Dimension> mesh_face_boundary = + getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor()); + + const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId(); + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + Array<const double> value_list = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, + mesh_data.xl(), + mesh_face_boundary + .faceList()); + + boundary_condition_list.push_back(NeumannBoundaryCondition{mesh_face_boundary.faceList(), value_list}); + + } else { + throw NotImplementedError("Dirichlet BC in 1d"); + } + break; + } + default: { + is_valid_boundary_condition = false; + } + } + if (not is_valid_boundary_condition) { + std::ostringstream error_msg; + error_msg << *bc_descriptor << " is an invalid boundary condition for heat equation"; + throw NormalError(error_msg.str()); + } + } + + 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, SymmetryBoundaryCondition>) 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()); + } else { + compute_face_dof_number[face_id] = number_of_dof++; + } + } + } + }, + boundary_condition); + } + + return compute_face_dof_number; + }(); + + const FaceValue<const bool> primal_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>) or + (std::is_same_v<T, SymmetryBoundaryCondition>)) { + 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); + } + + return face_is_neumann; + }(); + + const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix(); + NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity()); + if (parallel::size() > 1) { + throw NotImplementedError("Calculation of node_is_on_boundary is incorrect"); + } + + primal_node_is_on_boundary.fill(false); + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + if (face_to_cell_matrix[face_id].size() == 1) { + for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) { + NodeId node_id = primal_face_to_node_matrix[face_id][i_node]; + primal_node_is_on_boundary[node_id] = true; + } + } + } + + FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity()); + if (parallel::size() > 1) { + throw NotImplementedError("Calculation of face_is_on_boundary is incorrect"); + } + + primal_face_is_on_boundary.fill(false); + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + if (face_to_cell_matrix[face_id].size() == 1) { + primal_face_is_on_boundary[face_id] = true; + } + } + + FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity()); + if (parallel::size() > 1) { + throw NotImplementedError("Calculation of face_is_neumann is incorrect"); + } + + primal_face_is_dirichlet.fill(false); + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id])); + } + + InterpolationWeightsManager iwm(mesh, primal_face_is_on_boundary, primal_node_is_on_boundary); + iwm.compute(); + CellValuePerNode<double> w_rj = iwm.wrj(); + FaceValuePerNode<double> w_rl = iwm.wrl(); + + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl(); + const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj(); + const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix(); + + { + std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh); + + MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh); + + std::shared_ptr mapper = + DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity()); + + 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(); + + 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 CellValue<const double> dual_mes_l_j = [=] { + CellValue<double> compute_mes_j{diamond_mesh->connectivity()}; + mapper->toDualCell(mes_l, compute_mes_j); + + return compute_mes_j; + }(); + + FaceValue<const CellId> face_dual_cell_id = [=]() { + FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()}; + CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()}; + parallel_for( + diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; }); + + mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id); + + return computed_face_dual_cell_id; + }(); + + NodeValue<const NodeId> dual_node_primal_node_id = [=]() { + CellValue<NodeId> cell_ignored_id{mesh->connectivity()}; + cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()}); + + NodeValue<NodeId> node_primal_id{mesh->connectivity()}; + + parallel_for( + mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; }); + + NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()}; + + mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id); + + return computed_dual_node_primal_node_id; + }(); + + CellValue<NodeId> primal_cell_dual_node_id = [=]() { + CellValue<NodeId> cell_id{mesh->connectivity()}; + NodeValue<NodeId> node_ignored_id{mesh->connectivity()}; + node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()}); + + NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()}; + + parallel_for( + diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; }); + + CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()}; + + mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id); + + return cell_id; + }(); + const auto& dual_Cjr = diamond_mesh_data.Cjr(); + FaceValue<TinyVector<Dimension>> dualClj = [&] { + FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()}; + const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix(); + const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix(); + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { + const auto& primal_face_to_cell = face_to_cell_matrix[face_id]; + for (size_t i = 0; i < primal_face_to_cell.size(); i++) { + CellId cell_id = primal_face_to_cell[i]; + const NodeId dual_node_id = primal_cell_dual_node_id[cell_id]; + for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size(); + i_dual_cell++) { + const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell]; + if (face_dual_cell_id[face_id] == dual_cell_id) { + for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); + i_dual_node++) { + const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node]; + if (final_dual_node_id == dual_node_id) { + computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node); + } + } + } + } + } + }); + return computedClj; + }(); + + FaceValue<TinyVector<Dimension>> nlj = [&] { + FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()}; + parallel_for( + mesh->numberOfFaces(), + PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; }); + return computedNlj; + }(); + + FaceValue<const double> alpha_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_kappaj[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> 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(); + + const Array<int> non_zeros{number_of_dof}; + non_zeros.fill(Dimension); + CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zeros); + + 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) { + const CellId cell_id2 = primal_face_to_cell[j_cell]; + 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]) -= betab_l; + } + } else { + S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l; + } + } + } + } + + 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] * primal_Vj[cell_id]; + } + + 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_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]; + const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0); + + for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) { + NodeId node_id = primal_face_to_node_matrix[face_id][i_node]; + + const TinyVector<Dimension> nil = [&] { + if (is_face_reversed_for_cell_i) { + return -nlj[face_id]; + } else { + return nlj[face_id]; + } + }(); + + CellId dual_cell_id = face_dual_cell_id[face_id]; + + for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) { + const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node]; + 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 * dot(nil, Clr); + const double ab_ir = alphab_face_id * dot(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) * ab_ir; + } + } + 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) { + FaceId l_id = node_to_face_matrix[node_id][l_face]; + 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) * ab_ir; + } + } + } + } + } + } + } + } + } + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + if (primal_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{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] * primal_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]; + } + } + }, + boundary_condition); + } + + Vector<double> T{number_of_dof}; + T = zero; + + 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]]; }); + } + } + } +}; + +std::shared_ptr<const IDiscreteFunction> +ScalarNodalSchemeHandler::solution() const +{ + return m_scheme->getSolution(); +} + +ScalarNodalSchemeHandler::ScalarNodalSchemeHandler( + 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}); + if (not i_mesh) { + throw NormalError("primal discrete functions are not defined on the same mesh"); + } + const std::shared_ptr i_dual_mesh = getCommonMesh({dual_mub, dual_mu}); + if (not i_dual_mesh) { + throw NormalError("dual discrete functions are not defined on the same mesh"); + } + checkDiscretizationType({alpha, dual_mub, dual_mu, f}, DiscreteFunctionType::P0); + + switch (i_mesh->dimension()) { + case 1: { + using MeshType = Mesh<Connectivity<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 DiscreteScalarFunctionType>(dual_mub), + std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(dual_mu), + std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f), + bc_descriptor_list); + break; + } + case 2: { + using MeshType = Mesh<Connectivity<2>>; + using DiscreteScalarFunctionType = DiscreteFunctionP0<2, 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<2>>(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); + break; + } + case 3: { + using MeshType = Mesh<Connectivity<3>>; + using DiscreteScalarFunctionType = DiscreteFunctionP0<3, 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<3>>(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); + break; + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +ScalarNodalSchemeHandler::~ScalarNodalSchemeHandler() = default; diff --git a/src/scheme/ScalarNodalScheme.hpp b/src/scheme/ScalarNodalScheme.hpp new file mode 100644 index 000000000..fd6aaec18 --- /dev/null +++ b/src/scheme/ScalarNodalScheme.hpp @@ -0,0 +1,688 @@ +#ifndef SCALAR_NODAL_SCHEME_HPP +#define SCALAR_NODAL_SCHEME_HPP + +#include <algebra/CRSMatrix.hpp> +#include <algebra/CRSMatrixDescriptor.hpp> +#include <algebra/LeastSquareSolver.hpp> +#include <algebra/LinearSolver.hpp> +#include <algebra/TinyVector.hpp> +#include <algebra/Vector.hpp> +#include <language/utils/InterpolateItemValue.hpp> +#include <mesh/Connectivity.hpp> +#include <mesh/DualConnectivityManager.hpp> +#include <mesh/DualMeshManager.hpp> +#include <mesh/Mesh.hpp> +#include <mesh/MeshData.hpp> +#include <mesh/MeshDataManager.hpp> +#include <mesh/MeshFaceBoundary.hpp> +#include <mesh/MeshNodeBoundary.hpp> +#include <mesh/PrimalToDiamondDualConnectivityDataMapper.hpp> +#include <mesh/SubItemValuePerItem.hpp> +#include <scheme/DirichletBoundaryConditionDescriptor.hpp> +#include <scheme/FourierBoundaryConditionDescriptor.hpp> +#include <scheme/IDiscreteFunction.hpp> +#include <scheme/NeumannBoundaryConditionDescriptor.hpp> +#include <scheme/SymmetryBoundaryConditionDescriptor.hpp> + +class ScalarNodalSchemeHandler +{ + private: + class IScalarNodalScheme; + + template <size_t Dimension> + class ScalarNodalScheme; + + template <size_t Dimension> + class InterpolationWeightsManager; + + public: + std::unique_ptr<IScalarNodalScheme> m_scheme; + + std::shared_ptr<const IDiscreteFunction> solution() const; + + ScalarNodalSchemeHandler( + 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); + + ~ScalarNodalSchemeHandler(); +}; + +template <size_t Dimension> +class LegacyScalarNodalScheme +{ + private: + class DirichletBoundaryCondition + { + private: + const Array<const double> m_value_list; + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + + const Array<const double>& + valueList() const + { + return m_value_list; + } + + DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list) + : m_value_list{value_list}, m_face_list{face_list} + { + Assert(m_value_list.size() == m_face_list.size()); + } + + ~DirichletBoundaryCondition() = default; + }; + + class NeumannBoundaryCondition + { + private: + const Array<const double> m_value_list; + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_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} + { + Assert(m_value_list.size() == m_face_list.size()); + } + + ~NeumannBoundaryCondition() = default; + }; + + class FourierBoundaryCondition + { + private: + const Array<const double> m_coef_list; + const Array<const double> m_value_list; + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + + const Array<const double>& + valueList() const + { + return m_value_list; + } + + const Array<const double>& + coefList() const + { + return m_coef_list; + } + + public: + FourierBoundaryCondition(const Array<const FaceId>& face_list, + const Array<const double>& coef_list, + const Array<const double>& value_list) + : m_coef_list{coef_list}, m_value_list{value_list}, m_face_list{face_list} + { + Assert(m_coef_list.size() == m_face_list.size()); + Assert(m_value_list.size() == m_face_list.size()); + } + + ~FourierBoundaryCondition() = default; + }; + + class SymmetryBoundaryCondition + { + private: + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + + public: + SymmetryBoundaryCondition(const Array<const FaceId>& face_list) : m_face_list{face_list} {} + + ~SymmetryBoundaryCondition() = default; + }; + + public: + LegacyScalarNodalScheme(std::shared_ptr<const IMesh> i_mesh, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const FunctionSymbolId& kappa_id, + const FunctionSymbolId& f_id, + CellValue<double>& Temperature, + FaceValue<double>& Temperature_face, + const double& Tf, + const double& dt, + const CellValuePerNode<double>& w_rj, + const FaceValuePerNode<double>& w_rl) + { + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + using MeshDataType = MeshData<Dimension>; + + using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, + NeumannBoundaryCondition, SymmetryBoundaryCondition>; + + using BoundaryConditionList = std::vector<BoundaryCondition>; + + BoundaryConditionList boundary_condition_list; + + std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); + + for (const auto& bc_descriptor : bc_descriptor_list) { + bool is_valid_boundary_condition = true; + + switch (bc_descriptor->type()) { + case IBoundaryConditionDescriptor::Type::symmetry: { + throw NotImplementedError("NIY"); + break; + } + case IBoundaryConditionDescriptor::Type::dirichlet: { + const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor = + dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor); + if constexpr (Dimension > 1) { + MeshFaceBoundary<Dimension> mesh_face_boundary = + getMeshFaceBoundary(*mesh, dirichlet_bc_descriptor.boundaryDescriptor()); + + const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId(); + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + Array<const double> value_list = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, + mesh_data.xl(), + mesh_face_boundary + .faceList()); + boundary_condition_list.push_back(DirichletBoundaryCondition{mesh_face_boundary.faceList(), value_list}); + } + break; + } + case IBoundaryConditionDescriptor::Type::fourier: { + throw NotImplementedError("NIY"); + break; + } + case IBoundaryConditionDescriptor::Type::neumann: { + const NeumannBoundaryConditionDescriptor& neumann_bc_descriptor = + dynamic_cast<const NeumannBoundaryConditionDescriptor&>(*bc_descriptor); + + if constexpr (Dimension > 1) { + MeshFaceBoundary<Dimension> mesh_face_boundary = + getMeshFaceBoundary(*mesh, neumann_bc_descriptor.boundaryDescriptor()); + + const FunctionSymbolId g_id = neumann_bc_descriptor.rhsSymbolId(); + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + Array<const double> value_list = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, + mesh_data.xl(), + mesh_face_boundary + .faceList()); + + boundary_condition_list.push_back(NeumannBoundaryCondition{mesh_face_boundary.faceList(), value_list}); + + } else { + throw NotImplementedError("Dirichlet BC in 1d"); + } + break; + } + default: { + is_valid_boundary_condition = false; + } + } + if (not is_valid_boundary_condition) { + std::ostringstream error_msg; + error_msg << *bc_descriptor << " is an invalid boundary condition for heat equation"; + throw NormalError(error_msg.str()); + } + } + + 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, SymmetryBoundaryCondition>) 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()); + } else { + compute_face_dof_number[face_id] = number_of_dof++; + } + } + } + }, + boundary_condition); + } + + return compute_face_dof_number; + }(); + + const FaceValue<const bool> primal_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>) or + (std::is_same_v<T, SymmetryBoundaryCondition>)) { + 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); + } + + return face_is_neumann; + }(); + + const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix(); + NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity()); + if (parallel::size() > 1) { + throw NotImplementedError("Calculation of node_is_on_boundary is incorrect"); + } + + primal_node_is_on_boundary.fill(false); + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + if (face_to_cell_matrix[face_id].size() == 1) { + for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) { + NodeId node_id = primal_face_to_node_matrix[face_id][i_node]; + primal_node_is_on_boundary[node_id] = true; + } + } + } + + FaceValue<bool> primal_face_is_on_boundary(mesh->connectivity()); + if (parallel::size() > 1) { + throw NotImplementedError("Calculation of face_is_on_boundary is incorrect"); + } + + primal_face_is_on_boundary.fill(false); + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + if (face_to_cell_matrix[face_id].size() == 1) { + primal_face_is_on_boundary[face_id] = true; + } + } + + FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity()); + if (parallel::size() > 1) { + throw NotImplementedError("Calculation of face_is_neumann is incorrect"); + } + + primal_face_is_dirichlet.fill(false); + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + primal_face_is_dirichlet[face_id] = (primal_face_is_on_boundary[face_id] && (!primal_face_is_neumann[face_id])); + } + + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl(); + const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj(); + const auto& node_to_face_matrix = mesh->connectivity().nodeToFaceMatrix(); + + { + std::shared_ptr diamond_mesh = DualMeshManager::instance().getDiamondDualMesh(*mesh); + + MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh); + + std::shared_ptr mapper = + DualConnectivityManager::instance().getPrimalToDiamondDualConnectivityDataMapper(mesh->connectivity()); + + CellValue<double> kappaj = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id, + mesh_data.xj()); + + CellValue<double> dual_kappaj = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(kappa_id, + diamond_mesh_data + .xj()); + + const CellValue<const double> dual_Vj = diamond_mesh_data.Vj(); + + 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 CellValue<const double> dual_mes_l_j = [=] { + CellValue<double> compute_mes_j{diamond_mesh->connectivity()}; + mapper->toDualCell(mes_l, compute_mes_j); + + return compute_mes_j; + }(); + + FaceValue<const CellId> face_dual_cell_id = [=]() { + FaceValue<CellId> computed_face_dual_cell_id{mesh->connectivity()}; + CellValue<CellId> dual_cell_id{diamond_mesh->connectivity()}; + parallel_for( + diamond_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; }); + + mapper->fromDualCell(dual_cell_id, computed_face_dual_cell_id); + + return computed_face_dual_cell_id; + }(); + + NodeValue<const NodeId> dual_node_primal_node_id = [=]() { + CellValue<NodeId> cell_ignored_id{mesh->connectivity()}; + cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()}); + + NodeValue<NodeId> node_primal_id{mesh->connectivity()}; + + parallel_for( + mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; }); + + NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()}; + + mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id); + + return computed_dual_node_primal_node_id; + }(); + + CellValue<NodeId> primal_cell_dual_node_id = [=]() { + CellValue<NodeId> cell_id{mesh->connectivity()}; + NodeValue<NodeId> node_ignored_id{mesh->connectivity()}; + node_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()}); + + NodeValue<NodeId> dual_node_id{diamond_mesh->connectivity()}; + + parallel_for( + diamond_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { dual_node_id[node_id] = node_id; }); + + CellValue<NodeId> computed_primal_cell_dual_node_id{mesh->connectivity()}; + + mapper->fromDualNode(dual_node_id, node_ignored_id, cell_id); + + return cell_id; + }(); + const auto& dual_Cjr = diamond_mesh_data.Cjr(); + FaceValue<TinyVector<Dimension>> dualClj = [&] { + FaceValue<TinyVector<Dimension>> computedClj{mesh->connectivity()}; + const auto& dual_node_to_cell_matrix = diamond_mesh->connectivity().nodeToCellMatrix(); + const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix(); + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { + const auto& primal_face_to_cell = face_to_cell_matrix[face_id]; + for (size_t i = 0; i < primal_face_to_cell.size(); i++) { + CellId cell_id = primal_face_to_cell[i]; + const NodeId dual_node_id = primal_cell_dual_node_id[cell_id]; + for (size_t i_dual_cell = 0; i_dual_cell < dual_node_to_cell_matrix[dual_node_id].size(); + i_dual_cell++) { + const CellId dual_cell_id = dual_node_to_cell_matrix[dual_node_id][i_dual_cell]; + if (face_dual_cell_id[face_id] == dual_cell_id) { + for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); + i_dual_node++) { + const NodeId final_dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node]; + if (final_dual_node_id == dual_node_id) { + computedClj[face_id] = dual_Cjr(dual_cell_id, i_dual_node); + } + } + } + } + } + }); + return computedClj; + }(); + + FaceValue<TinyVector<Dimension>> nlj = [&] { + FaceValue<TinyVector<Dimension>> computedNlj{mesh->connectivity()}; + parallel_for( + mesh->numberOfFaces(), + PUGS_LAMBDA(FaceId face_id) { computedNlj[face_id] = 1. / l2Norm(dualClj[face_id]) * dualClj[face_id]; }); + return computedNlj; + }(); + + FaceValue<const double> alpha_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_kappaj[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; + }(); + + double lambda = (Tf == 0) ? 0 : 1; + double time_factor = (lambda == 0) ? 1 : dt; + + const CellValue<const double> primal_Vj = mesh_data.Vj(); + Array<int> non_zero{number_of_dof}; + non_zero.fill(Dimension * Dimension); + CRSMatrixDescriptor<double> S(number_of_dof, number_of_dof, non_zero); + 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]; + 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) { + const CellId cell_id2 = primal_face_to_cell[j_cell]; + if (i_cell == j_cell) { + S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) += + (time_factor * beta_l + lambda * primal_Vj[cell_id1]); + if (primal_face_is_neumann[face_id]) { + S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= beta_l; + } + } else { + S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= time_factor * beta_l; + } + } + } + } + + 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_face_id = mes_l[face_id] * alpha_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]; + const bool is_face_reversed_for_cell_i = (dot(dualClj[face_id], xl[face_id] - xj[i_id]) < 0); + + for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) { + NodeId node_id = primal_face_to_node_matrix[face_id][i_node]; + + const TinyVector<Dimension> nil = [&] { + if (is_face_reversed_for_cell_i) { + return -nlj[face_id]; + } else { + return nlj[face_id]; + } + }(); + + CellId dual_cell_id = face_dual_cell_id[face_id]; + + for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) { + const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node]; + 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 * dot(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]) -= time_factor * 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; + } + } + 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) { + FaceId l_id = node_to_face_matrix[node_id][l_face]; + if (primal_face_is_on_boundary[l_id]) { + S(cell_dof_number[i_id], face_dof_number[l_id]) -= time_factor * 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; + } + } + } + } + } + } + } + } + } + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + if (primal_face_is_dirichlet[face_id]) { + S(face_dof_number[face_id], face_dof_number[face_id]) += 1; + } + } + + CellValue<double> fj = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, + mesh_data.xj()); + + 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]] = + (time_factor * fj[cell_id] + lambda * Temperature[cell_id]) * primal_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]; + } + } + }, + boundary_condition); + } + + Vector<double> T{number_of_dof}; + T = zero; + + LinearSolver solver; + solver.solveLocalSystem(A, T, b); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_dof_number[cell_id]]; }); + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { + if (primal_face_is_neumann[face_id]) { + Temperature_face[face_id] = T[face_dof_number[face_id]]; + } + }); + } + } + } +}; + +template LegacyScalarNodalScheme<1>::LegacyScalarNodalScheme( + std::shared_ptr<const IMesh>, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&, + const FunctionSymbolId&, + const FunctionSymbolId&, + CellValue<double>&, + FaceValue<double>&, + const double&, + const double&, + const CellValuePerNode<double>& w_rj, + const FaceValuePerNode<double>& w_rl); + +template LegacyScalarNodalScheme<2>::LegacyScalarNodalScheme( + std::shared_ptr<const IMesh>, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&, + const FunctionSymbolId&, + const FunctionSymbolId&, + CellValue<double>&, + FaceValue<double>&, + const double&, + const double&, + const CellValuePerNode<double>& w_rj, + const FaceValuePerNode<double>& w_rl); + +template LegacyScalarNodalScheme<3>::LegacyScalarNodalScheme( + std::shared_ptr<const IMesh>, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&, + const FunctionSymbolId&, + const FunctionSymbolId&, + CellValue<double>&, + FaceValue<double>&, + const double&, + const double&, + const CellValuePerNode<double>& w_rj, + const FaceValuePerNode<double>& w_rl); + +#endif // SCALAR_NODAL_SCHEME_HPP -- GitLab