From 24e24145ff5a7eb41aaaab8da7dd7e23eba0e51d Mon Sep 17 00:00:00 2001 From: Drouard <axelle.drouard@orange.fr> Date: Tue, 3 May 2022 10:13:11 +0200 Subject: [PATCH] Replace scalar diffusion coefficient with tensorial one --- src/scheme/CMakeLists.txt | 3 +- src/scheme/ScalarNodalScheme.cpp | 271 ++----------- src/scheme/ScalarNodalScheme.hpp | 638 +------------------------------ 3 files changed, 31 insertions(+), 881 deletions(-) diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt index 25e1405a5..b7f170fb7 100644 --- a/src/scheme/CMakeLists.txt +++ b/src/scheme/CMakeLists.txt @@ -9,4 +9,5 @@ add_library( DiscreteFunctionVectorIntegrator.cpp DiscreteFunctionVectorInterpoler.cpp ScalarDiamondScheme.cpp - VectorDiamondScheme.cpp) + VectorDiamondScheme.cpp + ScalarNodalScheme.cpp) diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index fbbaf8952..fd4e2d5bd 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -3,137 +3,6 @@ #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 { @@ -249,23 +118,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand ~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> @@ -276,18 +128,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, double>>& dual_mub, - const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& dual_mu, + const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& nod_k_bound, + const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& nod_k, 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>; + NeumannBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -297,10 +145,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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); @@ -381,7 +225,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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(); @@ -411,8 +254,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand [&](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>)) { + (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) { @@ -466,10 +308,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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); @@ -485,8 +323,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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(); + CellValue<const TinyMatrix<Dimension>> nod_kappar = nod_k->cellValues(); + CellValue<const TinyMatrix<Dimension>> nod_kapparb = nod_k_bound->cellValues(); const CellValue<const double> dual_Vj = diamond_mesh_data.Vj(); @@ -588,29 +426,25 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }(); FaceValue<const double> alpha_l = [&] { - CellValue<double> alpha_j{diamond_mesh->connectivity()}; + FaceValue<double> alpha_j{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]; + mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { + alpha_j[face_id] = 1; }); - FaceValue<double> computed_alpha_l{mesh->connectivity()}; - mapper->fromDualCell(alpha_j, computed_alpha_l); - return computed_alpha_l; + return alpha_j; }(); FaceValue<const double> alphab_l = [&] { - CellValue<double> alpha_jb{diamond_mesh->connectivity()}; + FaceValue<double> alpha_lb{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]; + mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { + alpha_lb[face_id] = 1; //Refaire }); - FaceValue<double> computed_alpha_lb{mesh->connectivity()}; - mapper->fromDualCell(alpha_jb, computed_alpha_lb); - return computed_alpha_lb; + return alpha_lb; }(); const CellValue<const double> primal_Vj = mesh_data.Vj(); @@ -644,61 +478,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand 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; @@ -772,8 +552,8 @@ ScalarNodalSchemeHandler::solution() const 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>& nod_k_bound, + const std::shared_ptr<const IDiscreteFunction>& nod_k, const std::shared_ptr<const IDiscreteFunction>& f, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) { @@ -781,15 +561,16 @@ 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({dual_mub, dual_mu}); + const std::shared_ptr i_dual_mesh = getCommonMesh({nod_k_bound, nod_k}); 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); + checkDiscretizationType({alpha, nod_k_bound, nod_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); @@ -800,14 +581,15 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler( 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 DiscreteTensorFunctionType>(nod_k_bound), + std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k), std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f), bc_descriptor_list); break; } case 2: { using MeshType = Mesh<Connectivity<2>>; + using DiscreteTensorFunctionType = DiscreteFunctionP0<2, TinyMatrix<2>>; using DiscreteScalarFunctionType = DiscreteFunctionP0<2, double>; std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); @@ -818,14 +600,15 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler( 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 DiscreteTensorFunctionType>(nod_k_bound), + std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k), std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f), bc_descriptor_list); break; } case 3: { using MeshType = Mesh<Connectivity<3>>; + using DiscreteTensorFunctionType = DiscreteFunctionP0<3, TinyMatrix<3>>; using DiscreteScalarFunctionType = DiscreteFunctionP0<3, double>; std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); @@ -836,8 +619,8 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler( 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 DiscreteTensorFunctionType>(nod_k_bound), + std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(nod_k), std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f), bc_descriptor_list); break; diff --git a/src/scheme/ScalarNodalScheme.hpp b/src/scheme/ScalarNodalScheme.hpp index fd6aaec18..253beaf71 100644 --- a/src/scheme/ScalarNodalScheme.hpp +++ b/src/scheme/ScalarNodalScheme.hpp @@ -42,647 +42,13 @@ class ScalarNodalSchemeHandler 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>& nod_k_bound, + const std::shared_ptr<const IDiscreteFunction>& nod_k, 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