diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 8800b37b4ba8dd1d69387f3046e8e5cced712e77..abb60970a238402f360412b380ea05092f360740 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -603,23 +603,25 @@ SchemeModule::SchemeModule() } )); - this->_addBuiltinFunction("nodalheat", - std::make_shared<BuiltinFunctionEmbedder< - std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&, - const std::shared_ptr<const IDiscreteFunction>&, - const std::shared_ptr<const IDiscreteFunction>&, - const std::vector<std::shared_ptr< - const IBoundaryConditionDescriptor>>&)>>( - - [](const std::shared_ptr<const IDiscreteFunction>& mub_dual, - const std::shared_ptr<const IDiscreteFunction>& mu_dual, - const std::shared_ptr<const IDiscreteFunction>& f, - const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& - bc_descriptor_list) -> const std::shared_ptr<const IDiscreteFunction> { - return ScalarNodalSchemeHandler{mub_dual, mu_dual, f, bc_descriptor_list}.solution(); - } + this + ->_addBuiltinFunction("nodalheat", + std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr< + const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&, + const std::shared_ptr<const IDiscreteFunction>&, + const std::vector< + std::shared_ptr<const IBoundaryConditionDescriptor>>&, + const std::shared_ptr<const IDiscreteFunction>&, const double&)>>( + + [](const std::shared_ptr<const IDiscreteFunction>& kappa, + const std::shared_ptr<const IDiscreteFunction>& f, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list, + const std::shared_ptr<const IDiscreteFunction>& P, + const double& dt) -> const std::shared_ptr<const IDiscreteFunction> { + return ScalarNodalSchemeHandler{kappa, f, bc_descriptor_list, P, dt}.solution(); + } - )); + )); this->_addBuiltinFunction("unsteadyelasticity", std::make_shared<BuiltinFunctionEmbedder< diff --git a/src/scheme/ScalarNodalScheme.cpp b/src/scheme/ScalarNodalScheme.cpp index 2837089cda13b740a418518febcca49b6a4962c4..72936fe9222f2a9ac273f2190f37d79cf1513b39 100644 --- a/src/scheme/ScalarNodalScheme.cpp +++ b/src/scheme/ScalarNodalScheme.cpp @@ -143,10 +143,11 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand } ScalarNodalScheme(const std::shared_ptr<const MeshType>& mesh, - const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k_b, const std::shared_ptr<const DiscreteFunctionP0<Dimension, TinyMatrix<Dimension>>>& cell_k, const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& f, - const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>& P, + const double& dt) { using BoundaryCondition = std::variant<DirichletBoundaryCondition, FourierBoundaryCondition, NeumannBoundaryCondition>; @@ -408,8 +409,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand }(); { - CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues(); - CellValue<const TinyMatrix<Dimension>> cell_kappajb = cell_k_b->cellValues(); + CellValue<const TinyMatrix<Dimension>> cell_kappaj = cell_k->cellValues(); const NodeValue<const TinyMatrix<Dimension>> node_kappar = [&] { NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()}; @@ -429,24 +429,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand return kappa; }(); - const NodeValue<const TinyMatrix<Dimension>> node_kapparb = [&] { - NodeValue<TinyMatrix<Dimension>> kappa{mesh->connectivity()}; - parallel_for( - mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - kappa[node_id] = zero; - const auto& node_to_cell = node_to_cell_matrix[node_id]; - double weight = 0; - for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { - const CellId cell_id = node_to_cell[i_cell]; - double local_weight = 1. / l2Norm(xr[node_id] - xj[cell_id]); - kappa[node_id] += local_weight * cell_kappajb[cell_id]; - weight += local_weight; - } - kappa[node_id] = 1. / weight * kappa[node_id]; - }); - return kappa; - }(); - const NodeValue<const TinyMatrix<Dimension>> node_betar = [&] { NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()}; parallel_for( @@ -474,17 +456,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand return kappa_invBeta; }(); - const NodeValue<const TinyMatrix<Dimension>> kapparb_invBetar = [&] { - NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()}; - parallel_for( - mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - if (not node_is_corner[node_id]) { - kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(node_betar[node_id]); - } - }); - return kappa_invBeta; - }(); - const NodeValue<const TinyMatrix<Dimension>> corner_betar = [&] { NodeValue<TinyMatrix<Dimension>> beta{mesh->connectivity()}; parallel_for( @@ -538,17 +509,6 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand return kappa_invBeta; }(); - const NodeValue<const TinyMatrix<Dimension>> corner_kapparb_invBetar = [&] { - NodeValue<TinyMatrix<Dimension>> kappa_invBeta{mesh->connectivity()}; - parallel_for( - mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - if (node_is_corner[node_id]) { - kappa_invBeta[node_id] = node_kapparb[node_id] * inverse(corner_betar[node_id]); - } - }); - return kappa_invBeta; - }(); - const NodeValue<const TinyVector<Dimension>> sum_Cjr = [&] { NodeValue<TinyVector<Dimension>> compute_sum_Cjr{mesh->connectivity()}; parallel_for( @@ -570,8 +530,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand parallel_for( mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { if (is_boundary_node[node_id]) { - compute_v[node_id] = 1. / l2Norm(node_kapparb[node_id] * exterior_normal[node_id]) * - node_kapparb[node_id] * exterior_normal[node_id]; + compute_v[node_id] = 1. / l2Norm(node_kappar[node_id] * exterior_normal[node_id]) * node_kappar[node_id] * + exterior_normal[node_id]; } }); return compute_v; @@ -634,7 +594,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const TinyMatrix<Dimension> I = identity; const TinyVector<Dimension> n = exterior_normal[node_id]; const TinyMatrix<Dimension> nxn = tensorProduct(n, n); - const TinyMatrix<Dimension> P = I - nxn; + const TinyMatrix<Dimension> Q = I - nxn; for (size_t i_cell_j = 0; i_cell_j < node_to_cell.size(); ++i_cell_j) { const CellId cell_id_j = node_to_cell[i_cell_j]; @@ -647,7 +607,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand S(cell_id_j, cell_id_k) += dot(kappar_invBetar[node_id] * (Cjr(cell_id_k, i_node_k) - theta(cell_id_k, i_node_k) / sum_theta[node_id] * sum_Cjr[node_id]), - P * Cjr(cell_id_k, i_node_j)); + Q * Cjr(cell_id_k, i_node_j)); } } } else if ((node_is_dirichlet[node_id]) && (not node_is_corner[node_id])) { @@ -679,18 +639,25 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand } } + for (CellId cell_id_j = 0; cell_id_j < mesh->numberOfCells(); ++cell_id_j) { + for (CellId cell_id_k = 0; cell_id_k < mesh->numberOfCells(); ++cell_id_k) { + S(cell_id_j, cell_id_k) *= dt / Vj[cell_id_j]; + } + }; + for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - S(cell_id, cell_id) += Vj[cell_id]; + S(cell_id, cell_id) += 1; } CellValue<const double> fj = f->cellValues(); + CellValue<const double> Ph = P->cellValues(); CRSMatrix A{S.getCRSMatrix()}; Vector<double> b{mesh->numberOfCells()}; b = zero; for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) { - b[cell_id] = fj[cell_id] * Vj[cell_id]; + b[cell_id] = fj[cell_id]; }; for (NodeId node_id = 0; node_id < mesh->numberOfNodes(); ++node_id) { @@ -699,16 +666,16 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const TinyMatrix<Dimension> I = identity; const TinyVector<Dimension> n = exterior_normal[node_id]; const TinyMatrix<Dimension> nxn = tensorProduct(n, n); - const TinyMatrix<Dimension> P = I - nxn; + const TinyMatrix<Dimension> Q = I - nxn; if ((node_is_neumann[node_id]) && (not node_is_dirichlet[node_id])) { if ((is_boundary_node[node_id]) and (not node_is_corner[node_id])) { for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) { const CellId cell_id = node_to_cell[i_cell]; const size_t i_node = node_local_number_in_its_cells[i_cell]; - b[cell_id] += node_boundary_values[node_id] * - (1. / (sum_theta[node_id] * l2Norm(node_kapparb[node_id] * n)) * - dot(kapparb_invBetar[node_id] * sum_Cjr[node_id], P * Cjr(cell_id, i_node)) + + b[cell_id] += 1. / Vj[cell_id] * node_boundary_values[node_id] * + (1. / (sum_theta[node_id] * l2Norm(node_kappar[node_id] * n)) * + dot(kappar_invBetar[node_id] * sum_Cjr[node_id], Q * Cjr(cell_id, i_node)) + dot(Cjr(cell_id, i_node), n)); } } else if (node_is_corner[node_id]) { @@ -719,7 +686,7 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand FaceId face_id = node_to_face[i_face]; sum_mes_l += mes_l[face_id]; } - b[cell_id] += 0.5 * node_boundary_values[node_id] * sum_mes_l; + b[cell_id] += 1. / Vj[cell_id] * 0.5 * node_boundary_values[node_id] * sum_mes_l; } } else if (node_is_dirichlet[node_id]) { if (not node_is_corner[node_id]) { @@ -731,9 +698,9 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const CellId cell_id_k = node_to_cell[i_cell_k]; const size_t i_node_k = node_local_number_in_its_cells[i_cell_k]; - b[cell_id_j] += - dot(node_boundary_values[node_id] * kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k), - Cjr(cell_id_j, i_node_j)); + b[cell_id_j] += 1. / Vj[cell_id_j] * + dot(node_boundary_values[node_id] * kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), + Cjr(cell_id_j, i_node_j)); } } @@ -747,7 +714,8 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand const size_t i_node_k = node_local_number_in_its_cells[i_cell_k]; b[cell_id_j] += - dot(node_boundary_values[node_id] * corner_kapparb_invBetar[node_id] * Cjr(cell_id_k, i_node_k), + 1. / Vj[cell_id_j] * + dot(node_boundary_values[node_id] * corner_kappar_invBetar[node_id] * Cjr(cell_id_k, i_node_k), Cjr(cell_id_j, i_node_j)); } } @@ -757,8 +725,12 @@ class ScalarNodalSchemeHandler::ScalarNodalScheme : public ScalarNodalSchemeHand Vector<double> T{mesh->numberOfCells()}; T = zero; + Vector<double> B{mesh->numberOfCells()}; + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { B[cell_id] = dt * b[cell_id] + Ph[cell_id]; }); + LinearSolver solver; - solver.solveLocalSystem(A, T, b); + solver.solveLocalSystem(A, T, B); m_solution = std::make_shared<DiscreteFunctionP0<Dimension, double>>(mesh); auto& solution = *m_solution; @@ -776,17 +748,18 @@ ScalarNodalSchemeHandler::solution() const } ScalarNodalSchemeHandler::ScalarNodalSchemeHandler( - const std::shared_ptr<const IDiscreteFunction>& cell_k_b, const std::shared_ptr<const IDiscreteFunction>& cell_k, const std::shared_ptr<const IDiscreteFunction>& f, - const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const std::shared_ptr<const IDiscreteFunction>& P, + const double& dt) { const std::shared_ptr i_mesh = getCommonMesh({cell_k, f}); if (not i_mesh) { throw NormalError("primal discrete functions are not defined on the same mesh"); } - checkDiscretizationType({cell_k_b, cell_k, f}, DiscreteFunctionType::P0); + checkDiscretizationType({cell_k, f}, DiscreteFunctionType::P0); switch (i_mesh->dimension()) { case 1: { @@ -801,26 +774,14 @@ ScalarNodalSchemeHandler::ScalarNodalSchemeHandler( std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); m_scheme = - std::make_unique<ScalarNodalScheme<2>>(mesh, - std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b), - std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k), + std::make_unique<ScalarNodalScheme<2>>(mesh, std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k), std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f), - bc_descriptor_list); + bc_descriptor_list, + std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(P), dt); 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); - - m_scheme = - std::make_unique<ScalarNodalScheme<3>>(mesh, - std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k_b), - std::dynamic_pointer_cast<const DiscreteTensorFunctionType>(cell_k), - std::dynamic_pointer_cast<const DiscreteScalarFunctionType>(f), - bc_descriptor_list); + throw NormalError("The scheme is not implemented in dimension 3"); break; } default: { diff --git a/src/scheme/ScalarNodalScheme.hpp b/src/scheme/ScalarNodalScheme.hpp index 27ec1818efd0c17dc861aaf7b1c00413b5ec542d..a43ff3b4ca7b46b45948faeae252bdff62caa450 100644 --- a/src/scheme/ScalarNodalScheme.hpp +++ b/src/scheme/ScalarNodalScheme.hpp @@ -40,10 +40,11 @@ class ScalarNodalSchemeHandler std::shared_ptr<const IDiscreteFunction> solution() const; - ScalarNodalSchemeHandler(const std::shared_ptr<const IDiscreteFunction>& cell_k_bound, - const std::shared_ptr<const IDiscreteFunction>& cell_k, + ScalarNodalSchemeHandler(const std::shared_ptr<const IDiscreteFunction>& cell_k, const std::shared_ptr<const IDiscreteFunction>& f, - const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list); + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const std::shared_ptr<const IDiscreteFunction>& P, + const double& dt); ~ScalarNodalSchemeHandler(); };