diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt index 9faaa98dae7d2ebecb141d79d2b0030270ff9b62..087b6828b85e3eecde6a9b031a1dca21f0c749bb 100644 --- a/src/language/algorithms/CMakeLists.txt +++ b/src/language/algorithms/CMakeLists.txt @@ -3,7 +3,6 @@ add_library(PugsLanguageAlgorithms AcousticSolverAlgorithm.cpp ElasticityDiamondAlgorithm.cpp - ElasticityDiamondAlgorithm2.cpp Heat5PointsAlgorithm.cpp HeatDiamondAlgorithm.cpp ) diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp index e5c7d3c2786d874e00fa3f4a4247cce725fc693e..0dd1739e0cbcdcd6c85d3e1817fa5ed58de833c7 100644 --- a/src/language/algorithms/ElasticityDiamondAlgorithm.cpp +++ b/src/language/algorithms/ElasticityDiamondAlgorithm.cpp @@ -82,50 +82,17 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); const RefId& ref = ref_face_list.refId(); if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) { - MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list}; - const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId(); - const auto& node_list = mesh_node_boundary.nodeList(); - - Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>( - TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(), node_list); - - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - NodeId node_id = node_list[i_node]; - if (not is_dirichlet[node_id]) { - is_dirichlet[node_id] = true; - dirichlet_value[node_id] = value_list[i_node]; - } - } - - boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list}); - } - } - - if constexpr (Dimension > 1) { - for (size_t i_ref_node_list = 0; - i_ref_node_list < mesh->connectivity().template numberOfRefItemList<ItemType::node>(); - ++i_ref_node_list) { - const auto& ref_node_list = mesh->connectivity().template refItemList<ItemType::node>(i_ref_node_list); - const RefId& ref = ref_node_list.refId(); - if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) { - const FunctionSymbolId g_id = dirichlet_bc_descriptor.rhsSymbolId(); - - const auto& node_list = ref_node_list.list(); + if constexpr (Dimension > 1) { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + Array<const FaceId> face_list = ref_face_list.list(); Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>( - TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(), node_list); - - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - NodeId node_id = node_list[i_node]; - if (not is_dirichlet[node_id]) { - is_dirichlet[node_id] = true; - dirichlet_value[node_id] = value_list[i_node]; - } - } - - boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list}); + TinyVector<Dimension>)>::template interpolate<ItemType::face>(g_id, mesh_data.xl(), face_list); + boundary_condition_list.push_back(DirichletBoundaryCondition{face_list, value_list}); + } else { + throw NotImplementedError("Neumann conditions are not supported in 1d"); } } } @@ -182,7 +149,8 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( [&](auto&& bc) { using T = std::decay_t<decltype(bc)>; if constexpr ((std::is_same_v<T, NormalStrainBoundaryCondition>) or - (std::is_same_v<T, SymmetryBoundaryCondition>)) { + (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) { @@ -203,8 +171,30 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( return compute_face_dof_number; }(); - const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); - const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix(); + const auto& primal_face_to_node_matrix = mesh->connectivity().faceToNodeMatrix(); + const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix(); + 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, NormalStrainBoundaryCondition>)) { + 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; + }(); + NodeValue<bool> primal_node_is_on_boundary(mesh->connectivity()); if (parallel::size() > 1) { throw NotImplementedError("Calculation of node_is_on_boundary is incorrect"); @@ -232,21 +222,14 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( } } - FaceValue<bool> primal_face_is_neumann(mesh->connectivity()); + FaceValue<bool> primal_face_is_dirichlet(mesh->connectivity()); if (parallel::size() > 1) { throw NotImplementedError("Calculation of face_is_neumann is incorrect"); } - primal_face_is_neumann.fill(false); + primal_face_is_dirichlet.fill(false); for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { - if (primal_face_is_on_boundary[face_id]) { - 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]; - if (not is_dirichlet[node_id]) { - primal_face_is_neumann[face_id] = true; - } - } - } + 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); @@ -447,10 +430,10 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( return computed_dual_node_primal_node_id; }(); - TinyMatrix<Dimension, double> eye = zero; - for (size_t i = 0; i < Dimension; ++i) { - eye(i, i) = 1; - } + TinyMatrix<Dimension, double> eye = identity; + // for (size_t i = 0; i < Dimension; ++i) { + // eye(i, i) = 1; + // } const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr(); const auto& primal_face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed(); @@ -479,6 +462,9 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( for (size_t i = 0; i < Dimension; ++i) { for (size_t j = 0; j < Dimension; ++j) { S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) += M(i, j); + if (primal_face_is_neumann[face_id]) { + S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) -= M(i, j); + } } } } else { @@ -507,38 +493,70 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( 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]; - if (not primal_node_is_on_boundary[node_id]) { - const TinyVector<Dimension> nil = [&] { - if (is_face_reversed_for_cell_i) { - return -primal_nlr(face_id, i_node); - } else { - return primal_nlr(face_id, i_node); - } - }(); - - 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); - - TinyMatrix<Dimension, double> M = alpha_mu_face_id * (Clr, nil) * eye + - alpha_mu_face_id * tensorProduct(Clr, nil) + - alpha_lambda_face_id * tensorProduct(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]; - for (size_t i = 0; i < Dimension; ++i) { - for (size_t j = 0; j < Dimension; ++j) { - S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -= + // if (not primal_node_is_on_boundary[node_id]) { + const TinyVector<Dimension> nil = [&] { + if (is_face_reversed_for_cell_i) { + return -primal_nlr(face_id, i_node); + } else { + return primal_nlr(face_id, i_node); + } + }(); + + 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); + + TinyMatrix<Dimension, double> M = alpha_mu_face_id * (Clr, nil) * eye + + alpha_mu_face_id * tensorProduct(Clr, nil) + + alpha_lambda_face_id * tensorProduct(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]; + for (size_t i = 0; i < Dimension; ++i) { + for (size_t j = 0; j < Dimension; ++j) { + S((cell_dof_number[i_id] * Dimension) + i, (cell_dof_number[j_id] * Dimension) + j) -= + w_rj(node_id, j_cell) * M(i, j); + if (primal_face_is_neumann[face_id]) { + S(face_dof_number[face_id] * Dimension + i, cell_dof_number[j_id] * Dimension + j) += w_rj(node_id, j_cell) * M(i, j); } } } } + 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]) { + for (size_t i = 0; i < Dimension; ++i) { + for (size_t j = 0; j < Dimension; ++j) { + S(cell_dof_number[i_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) -= + w_rl(node_id, l_face) * M(i, j); + } + } + if (primal_face_is_neumann[face_id]) { + for (size_t i = 0; i < Dimension; ++i) { + for (size_t j = 0; j < Dimension; ++j) { + S(face_dof_number[face_id] * Dimension + i, face_dof_number[l_id] * Dimension + j) += + w_rl(node_id, l_face) * M(i, j); + } + } + } + } + } + } } } + // } + } + } + } + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + if (primal_face_is_dirichlet[face_id]) { + for (size_t i = 0; i < Dimension; ++i) { + S(face_dof_number[face_id] * Dimension + i, face_dof_number[face_id] * Dimension + i) += 1; } } } @@ -551,7 +569,6 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( } } - const auto& primal_cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix(); // Dirichlet NodeValue<bool> node_tag{mesh->connectivity()}; node_tag.fill(false); @@ -560,64 +577,32 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( [&](auto&& bc) { using T = std::decay_t<decltype(bc)>; if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) { - const auto& node_list = bc.nodeList(); + const auto& face_list = bc.faceList(); const auto& value_list = bc.valueList(); - for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - const NodeId node_id = node_list[i_node]; - if (not node_tag[node_id]) { - node_tag[node_id] = true; - - const auto& node_cells = primal_node_to_cell_matrix[node_id]; - for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) { - const CellId cell_id = node_cells[i_cell]; - - const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id]; + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; - for (size_t i_cell_face = 0; i_cell_face < primal_cell_to_face.size(); ++i_cell_face) { - FaceId face_id = primal_cell_to_face_matrix[cell_id][i_cell_face]; + for (size_t i = 0; i < Dimension; ++i) { + b[(face_dof_number[face_id] * Dimension) + i] += value_list[i_face][i]; + } + } + } + }, + boundary_condition); + } - const bool is_face_reversed_for_cell_i = [=] { - for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) { - FaceId primal_face_id = primal_cell_to_face[i_face]; - if (primal_face_id == face_id) { - return primal_face_cell_is_reversed(cell_id, i_face); - } - } - Assert(false, "cannot find cell's face"); - return false; - }(); - - const auto& primal_face_to_node = primal_face_to_node_matrix[face_id]; - - for (size_t i_face_node = 0; i_face_node < primal_face_to_node.size(); ++i_face_node) { - if (primal_face_to_node[i_face_node] == node_id) { - const TinyVector<Dimension> nil = [=] { - if (is_face_reversed_for_cell_i) { - return -primal_nlr(face_id, i_face_node); - } else { - return primal_nlr(face_id, i_face_node); - } - }(); - - 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); - TinyMatrix<Dimension, double> M = alpha_mu_l[face_id] * (Clr, nil) * eye + - alpha_mu_l[face_id] * tensorProduct(Clr, nil) + - alpha_lambda_l[face_id] * tensorProduct(nil, Clr); - for (size_t i = 0; i < Dimension; ++i) { - b[(cell_dof_number[cell_id] * Dimension) + i] += (M * value_list[i_node])[i]; - } - } - } - } - } - } - } + 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, NormalStrainBoundaryCondition>)) { + 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]; + Assert(face_to_cell_matrix[face_id].size() == 1); + for (size_t i = 0; i < Dimension; ++i) { + b[face_dof_number[face_id] * Dimension + i] += mes_l[face_id] * value_list[i_face][i]; // sign } } } @@ -626,36 +611,75 @@ ElasticityDiamondScheme<Dimension>::ElasticityDiamondScheme( } Vector<double> U{number_of_dof * Dimension}; - U = 0; + U = 1; + CellValue<TinyVector<Dimension>> Speed{mesh->connectivity()}; + FaceValue<TinyVector<Dimension>> Ul = InterpolateItemValue<TinyVector<Dimension>( + TinyVector<Dimension>)>::template interpolate<ItemType::face>(U_id, mesh_data.xl()); + FaceValue<TinyVector<Dimension>> Speed_face{mesh->connectivity()}; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + for (size_t i = 0; i < Dimension; ++i) { + U[(cell_dof_number[cell_id] * Dimension) + i] = Uj[cell_id][i]; + } + } + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + for (size_t i = 0; i < Dimension; ++i) { + if (primal_face_is_on_boundary[face_id]) { + U[(face_dof_number[face_id] * Dimension) + i] = Ul[face_id][i]; + } + } + } + Vector r = A * U - b; + + std::cout << "initial residu = " << std::sqrt((r, r)) << '\n'; BiCGStab{b, A, U, 10000, 1e-15}; - Vector r = A * U - b; + r = A * U - b; std::cout << "real residu = " << std::sqrt((r, r)) << '\n'; - CellValue<TinyVector<Dimension>> Speed{mesh->connectivity()}; - for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { for (size_t i = 0; i < Dimension; ++i) { Speed[cell_id][i] = U[(cell_dof_number[cell_id] * Dimension) + i]; } } - Vector<double> Uexacte{number_of_dof * Dimension}; - for (CellId j = 0; j < number_of_dof; ++j) { + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { + for (size_t i = 0; i < Dimension; ++i) { + if (primal_face_is_on_boundary[face_id]) { + Speed_face[face_id][i] = U[(face_dof_number[face_id] * Dimension) + i]; + } else { + Speed_face[face_id][i] = Ul[face_id][i]; + } + } + } + Vector<double> Uexacte{mesh->numberOfCells() * Dimension}; + for (CellId j = 0; j < mesh->numberOfCells(); ++j) { for (size_t l = 0; l < Dimension; ++l) { Uexacte[(cell_dof_number[j] * Dimension) + l] = Uj[j][l]; } } - Vector<double> error{number_of_dof * Dimension}; + Vector<double> error{mesh->numberOfCells() * Dimension}; for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { for (size_t i = 0; i < Dimension; ++i) { error[(cell_id * Dimension) + i] = (Speed[cell_id][i] - Uj[cell_id][i]) * sqrt(primal_Vj[cell_id]); } } + Vector<double> error_face{mesh->numberOfFaces() * Dimension}; + parallel_for( + mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { + if (primal_face_is_on_boundary[face_id]) { + for (size_t i = 0; i < Dimension; ++i) { + error_face[face_id * Dimension + i] = (Speed_face[face_id][i] - Ul[face_id][i]) * sqrt(mes_l[face_id]); + } + } else { + error_face[face_id] = 0; + } + }); - std::cout << "||Error||_2 = " << std::sqrt((error, error)) << "\n"; + std::cout << "||Error||_2 (cell)= " << std::sqrt((error, error)) << "\n"; + std::cout << "||Error||_2 (face)= " << std::sqrt((error_face, error_face)) << "\n"; + std::cout << "||Error||_2 (total)= " << std::sqrt((error, error)) + std::sqrt((error_face, error_face)) << "\n"; VTKWriter vtk_writer("Speed_" + std::to_string(Dimension), 0.01); @@ -672,13 +696,13 @@ class ElasticityDiamondScheme<Dimension>::DirichletBoundaryCondition { private: const Array<const TinyVector<Dimension>> m_value_list; - const Array<const NodeId> m_node_list; + const Array<const FaceId> m_face_list; public: - const Array<const NodeId>& - nodeList() const + const Array<const FaceId>& + faceList() const { - return m_node_list; + return m_face_list; } const Array<const TinyVector<Dimension>>& @@ -687,10 +711,10 @@ class ElasticityDiamondScheme<Dimension>::DirichletBoundaryCondition return m_value_list; } - DirichletBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list) - : m_value_list{value_list}, m_node_list{node_list} + DirichletBoundaryCondition(const Array<const FaceId>& face_list, const Array<const TinyVector<Dimension>>& value_list) + : m_value_list{value_list}, m_face_list{face_list} { - Assert(m_value_list.size() == m_node_list.size()); + Assert(m_value_list.size() == m_face_list.size()); } ~DirichletBoundaryCondition() = default; diff --git a/src/language/algorithms/ElasticityDiamondAlgorithm.hpp b/src/language/algorithms/ElasticityDiamondAlgorithm.hpp index 733aa7d580a10b983b3fe0e686c4473ca8d2a706..2d6c298eb6ea26d6ff626a1c9de3ee4162ae5721 100644 --- a/src/language/algorithms/ElasticityDiamondAlgorithm.hpp +++ b/src/language/algorithms/ElasticityDiamondAlgorithm.hpp @@ -26,4 +26,4 @@ class ElasticityDiamondScheme const FunctionSymbolId& U_id); }; -#endif // ELASTICITY_DIAMOND_ALGORITHM_HPP +#endif // ELASTICITY_DIAMOND_ALGORITHM2_HPP diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 4730fb2ae5f7250ae5f003a3d94fe805897032d6..61c90914a45536cf839e098b7b58eaddcc4cc3f1 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -2,7 +2,6 @@ #include <language/algorithms/AcousticSolverAlgorithm.hpp> #include <language/algorithms/ElasticityDiamondAlgorithm.hpp> -#include <language/algorithms/ElasticityDiamondAlgorithm2.hpp> #include <language/algorithms/Heat5PointsAlgorithm.hpp> #include <language/algorithms/HeatDiamondAlgorithm.hpp> #include <language/utils/BuiltinFunctionEmbedder.hpp> @@ -294,39 +293,6 @@ SchemeModule::SchemeModule() )); - this->_addBuiltinFunction("elasticity2", - std::make_shared<BuiltinFunctionEmbedder< - void(std::shared_ptr<const IMesh>, - const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&, - const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&, - const FunctionSymbolId&)>>( - - [](std::shared_ptr<const IMesh> p_mesh, - const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& - bc_descriptor_list, - const FunctionSymbolId& lambda_id, const FunctionSymbolId& mu_id, - const FunctionSymbolId& f_id, const FunctionSymbolId& U_id) -> void { - switch (p_mesh->dimension()) { - case 1: { - ElasticityDiamondScheme2<1>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id}; - break; - } - case 2: { - ElasticityDiamondScheme2<2>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id}; - break; - } - case 3: { - ElasticityDiamondScheme2<3>{p_mesh, bc_descriptor_list, lambda_id, mu_id, f_id, U_id}; - break; - } - default: { - throw UnexpectedError("invalid mesh dimension"); - } - } - } - - )); - this->_addBuiltinFunction("heat5points", std::make_shared<BuiltinFunctionEmbedder< void(std::shared_ptr<const IMesh>,