diff --git a/src/language/algorithms/HeatDiamondAlgorithm.cpp b/src/language/algorithms/HeatDiamondAlgorithm.cpp index 5cf0b0c5d6cdc73719a38b12a1e235933afbbf0f..aef3cd1aeae0ae664e6cc1e18347d416f4c49539 100644 --- a/src/language/algorithms/HeatDiamondAlgorithm.cpp +++ b/src/language/algorithms/HeatDiamondAlgorithm.cpp @@ -28,8 +28,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id, - const FunctionSymbolId& f_id, - const FunctionSymbolId& g_id) + const FunctionSymbolId& f_id) { using ConnectivityType = Connectivity<Dimension>; using MeshType = Mesh<ConnectivityType>; @@ -51,7 +50,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( BoundaryConditionList boundary_condition_list; std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n'; - std::shared_ptr m_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); + 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; @@ -67,20 +66,19 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor = dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor); - for (size_t i_ref_face_list = 0; - i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { - const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); + for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); + ++i_ref_face_list) { + 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{m_mesh, ref_face_list}; + 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 double> value_list = - InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, - m_mesh->xr(), + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, mesh->xr(), node_list); boundary_condition_list.push_back(DirichletBoundaryCondition{node_list, value_list}); @@ -112,18 +110,27 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( } if constexpr (Dimension > 1) { - MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh); + 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; + }(); + + auto dof_number = [&](CellId cell_id) { return cell_dof_number[cell_id]; }; + + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); CellValue<double> Tj = InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj()); - NodeValue<double> Tr(m_mesh->connectivity()); - const NodeValue<const TinyVector<Dimension>>& xr = m_mesh->xr(); + NodeValue<double> Tr(mesh->connectivity()); + const NodeValue<const TinyVector<Dimension>>& xr = mesh->xr(); const CellValue<const TinyVector<Dimension>>& xj = mesh_data.xj(); - const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix(); - CellValuePerNode<double> w_rj{m_mesh->connectivity()}; + const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix(); + CellValuePerNode<double> w_rj{mesh->connectivity()}; - for (NodeId i_node = 0; i_node < m_mesh->numberOfNodes(); ++i_node) { + for (NodeId i_node = 0; i_node < mesh->numberOfNodes(); ++i_node) { Vector<double> b{Dimension + 1}; b[0] = 1; for (size_t i = 1; i < Dimension + 1; i++) { @@ -146,7 +153,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( LocalRectangularMatrix B = transpose(A) * A; Vector y = transpose(A) * b; - PCG{y, B, B, x, 10, 1e-12}; + PCG<false>{y, B, B, x, 10, 1e-12}; Tr[i_node] = 0; for (size_t j = 0; j < node_to_cell.size(); j++) { @@ -158,22 +165,22 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( { VTKWriter vtk_writer("T_" + std::to_string(Dimension), 0.01); - vtk_writer.write(m_mesh, + vtk_writer.write(mesh, {NamedItemValue{"Tr", Tr}, NamedItemValue{"temperature", Tj}, - NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()}, - NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()}, - NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}}, + NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()}, + NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()}, + NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}}, 0, true); // forces last output } { - std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(m_mesh); + std::shared_ptr diamond_mesh = DiamondDualMeshManager::instance().getDiamondDualMesh(mesh); MeshDataType& diamond_mesh_data = MeshDataManager::instance().getMeshData(*diamond_mesh); std::shared_ptr mapper = DiamondDualConnectivityManager::instance().getConnectivityToDiamondDualConnectivityDataMapper( - m_mesh->connectivity()); + mesh->connectivity()); NodeValue<double> Trd{diamond_mesh->connectivity()}; @@ -199,7 +206,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( const FaceValue<const double> mes_l = [&] { if constexpr (Dimension == 1) { - FaceValue<double> compute_mes_l{m_mesh->connectivity()}; + FaceValue<double> compute_mes_l{mesh->connectivity()}; compute_mes_l.fill(1); return compute_mes_l; } else { @@ -223,7 +230,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( dual_mes_l_j[diamond_cell_id] * dual_kappaj[diamond_cell_id] / dual_Vj[diamond_cell_id]; }); - FaceValue<double> computed_alpha_l{m_mesh->connectivity()}; + FaceValue<double> computed_alpha_l{mesh->connectivity()}; mapper->fromDualCell(alpha_j, computed_alpha_l); VTKWriter vtk_writer("DD_" + std::to_string(Dimension), 0.01); @@ -237,16 +244,16 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( return computed_alpha_l; }(); - const auto& primal_face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix(); - const auto& face_to_cell_matrix = m_mesh->connectivity().faceToCellMatrix(); + 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(m_mesh->connectivity()); + 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 < m_mesh->numberOfFaces(); ++face_id) { + 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]; @@ -255,8 +262,8 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( } } - SparseMatrixDescriptor S{m_mesh->numberOfCells()}; - for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) { + SparseMatrixDescriptor S{mesh->numberOfCells()}; + 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 = 0.5 * alpha_l[face_id] * mes_l[face_id]; @@ -265,16 +272,16 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( 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_id1, cell_id2) += beta_l; + S(dof_number(cell_id1), dof_number(cell_id2)) += beta_l; } else { - S(cell_id1, cell_id2) -= beta_l; + S(dof_number(cell_id1), dof_number(cell_id2)) -= beta_l; } } } } FaceValue<const CellId> face_dual_cell_id = [=]() { - FaceValue<CellId> computed_face_dual_cell_id{m_mesh->connectivity()}; + 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; }); @@ -285,13 +292,13 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( }(); NodeValue<const NodeId> dual_node_primal_node_id = [=]() { - CellValue<NodeId> cell_ignored_id{m_mesh->connectivity()}; + CellValue<NodeId> cell_ignored_id{mesh->connectivity()}; cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()}); - NodeValue<NodeId> node_primal_id{m_mesh->connectivity()}; + NodeValue<NodeId> node_primal_id{mesh->connectivity()}; parallel_for( - m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; }); + 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()}; @@ -305,11 +312,11 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix(); const NodeValuePerFace<const TinyVector<Dimension>> primal_nlr = mesh_data.nlr(); - const auto& primal_face_cell_is_reversed = m_mesh->connectivity().cellFaceIsReversed(); - const auto& face_local_numbers_in_their_cells = m_mesh->connectivity().faceLocalNumbersInTheirCells(); - const auto& primal_node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix(); + const auto& primal_face_cell_is_reversed = mesh->connectivity().cellFaceIsReversed(); + const auto& face_local_numbers_in_their_cells = mesh->connectivity().faceLocalNumbersInTheirCells(); + const auto& primal_node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix(); - for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) { + for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) { if (face_to_cell_matrix[face_id].size() > 1) { const double alpha_face_id = alpha_l[face_id]; @@ -342,7 +349,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( 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(i_id, j_id) -= w_rj(node_id, j_cell) * a_ir; + S(dof_number(i_id), dof_number(j_id)) -= w_rj(node_id, j_cell) * a_ir; } } } @@ -355,75 +362,99 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( CellValue<double> fj = InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj()); - NodeValue<double> gr = - InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::node>(g_id, m_mesh->xr()); - const CellValue<const double> primal_Vj = mesh_data.Vj(); CRSMatrix A{S}; - Vector<double> b{m_mesh->numberOfCells()}; - - const auto& primal_cell_to_face_matrix = m_mesh->connectivity().cellToFaceMatrix(); - for (CellId cell_id = 0; cell_id < m_mesh->numberOfCells(); ++cell_id) { - b[cell_id] = fj[cell_id] * primal_Vj[cell_id]; - - const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id]; - - 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]; - - 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_node = 0; i_node < primal_face_to_node.size(); ++i_node) { - 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]; + Vector<double> b{mesh->numberOfCells()}; - NodeId face_node_id = primal_face_to_node[i_node]; + const auto& primal_cell_to_face_matrix = mesh->connectivity().cellToFaceMatrix(); + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { + b[dof_number(cell_id)] = fj[cell_id] * primal_Vj[cell_id]; + } - if (primal_node_is_on_boundary[face_node_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] == face_node_id) { - const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node); - b[cell_id] += alpha_l[face_id] * gr[face_node_id] * (nil, Clr); + { // Dirichlet + NodeValue<bool> node_tag{mesh->connectivity()}; + node_tag.fill(false); + for (const auto& boundary_condition : boundary_condition_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<T, DirichletBoundaryCondition>) { + const auto& node_list = bc.nodeList(); + 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_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]; + + 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); + b[dof_number(cell_id)] += alpha_l[face_id] * value_list[i_node] * (nil, Clr); + } + } + } + } + } + } + } } } - } - } + }, + boundary_condition); } + // std::exit(0); } - Vector<double> T{m_mesh->numberOfCells()}; + Vector<double> T{mesh->numberOfCells()}; T = 0; BiCGStab{b, A, T, 1000, 1e-15}; - CellValue<double> Temperature{m_mesh->connectivity()}; + CellValue<double> Temperature{mesh->connectivity()}; parallel_for( - m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[cell_id]; }); + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { Temperature[cell_id] = T[dof_number(cell_id)]; }); - Vector<double> error{m_mesh->numberOfCells()}; + Vector<double> error{mesh->numberOfCells()}; parallel_for( - m_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { error[cell_id] = (Temperature[cell_id] - Tj[cell_id]) * sqrt(primal_Vj[cell_id]); }); @@ -432,7 +463,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme( { VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01); - vtk_writer.write(m_mesh, + vtk_writer.write(mesh, {NamedItemValue{"T", Temperature}, NamedItemValue{"Vj", primal_Vj}, NamedItemValue{"Texact", Tj}}, 0, @@ -497,7 +528,6 @@ template HeatDiamondScheme<1>::HeatDiamondScheme( const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&, const FunctionSymbolId&, const FunctionSymbolId&, - const FunctionSymbolId&, const FunctionSymbolId&); template HeatDiamondScheme<2>::HeatDiamondScheme( @@ -505,7 +535,6 @@ template HeatDiamondScheme<2>::HeatDiamondScheme( const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&, const FunctionSymbolId&, const FunctionSymbolId&, - const FunctionSymbolId&, const FunctionSymbolId&); template HeatDiamondScheme<3>::HeatDiamondScheme( @@ -513,5 +542,4 @@ template HeatDiamondScheme<3>::HeatDiamondScheme( const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&, const FunctionSymbolId&, const FunctionSymbolId&, - const FunctionSymbolId&, const FunctionSymbolId&); diff --git a/src/language/algorithms/HeatDiamondAlgorithm.hpp b/src/language/algorithms/HeatDiamondAlgorithm.hpp index c3b3d2e84ff6e13c9fefb32303e7eb2d4b564c01..9e29d91bc1bc9f37842b8e993a698a4d26d9f68f 100644 --- a/src/language/algorithms/HeatDiamondAlgorithm.hpp +++ b/src/language/algorithms/HeatDiamondAlgorithm.hpp @@ -23,8 +23,7 @@ class HeatDiamondScheme const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, const FunctionSymbolId& T_id, const FunctionSymbolId& kappa_id, - const FunctionSymbolId& f_id, - const FunctionSymbolId& g_id); + const FunctionSymbolId& f_id); }; #endif // HEAT_DIAMOND_ALGORITHM_HPP diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index d2d25382cbb376c6fab9ae67c8044bf719e47b42..9ac6a6d09b48510c6ade818a13bb6e81e5dc034e 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -218,25 +218,24 @@ SchemeModule::SchemeModule() this->_addBuiltinFunction("heat", 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&)>>( + 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& T_id, const FunctionSymbolId& kappa_id, - const FunctionSymbolId& f_id, const FunctionSymbolId& g_id) -> void { + const FunctionSymbolId& f_id) -> void { switch (p_mesh->dimension()) { case 1: { - HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id, g_id}; + HeatDiamondScheme<1>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id}; break; } case 2: { - HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id, g_id}; + HeatDiamondScheme<2>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id}; break; } case 3: { - HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id, g_id}; + HeatDiamondScheme<3>{p_mesh, bc_descriptor_list, T_id, kappa_id, f_id}; break; } default: {