diff --git a/src/scheme/Order2LocalDtHyperelasticSolver.cpp b/src/scheme/Order2LocalDtHyperelasticSolver.cpp index 824790cecc0a15899b494f92a8d9272a1d71001d..ff0c44b0fffcf57b97c71c87970f342984d94942 100644 --- a/src/scheme/Order2LocalDtHyperelasticSolver.cpp +++ b/src/scheme/Order2LocalDtHyperelasticSolver.cpp @@ -450,6 +450,228 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi return std::make_tuple(map1, map2); } + void + _vector_limiter(const MeshType& mesh, + const DiscreteFunctionP0<const Rd>& f, + DiscreteFunctionDPk<Dimension, Rd>& DPk_fh) const + { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list; + StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes}; + auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list); + + const auto xj = mesh_data.xj(); + + parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){ + const Rd fj = f[cell_id]; + + Rd f_min = fj; + Rd f_max = fj; + + const auto cell_stencil = stencil[cell_id]; + for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){ + for(size_t dim = 0; dim < Dimension; ++dim){ + f_min[dim] = std::min(f_min[dim], f[cell_stencil[i_cell]][dim]); + f_max[dim] = std::max(f_max[dim], f[cell_stencil[i_cell]][dim]); + } + } + + Rd f_bar_min = fj; + Rd f_bar_max = fj; + + for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){ + const CellId cell_k_id = cell_stencil[i_cell]; + const Rd f_xk = DPk_fh[cell_id](xj[cell_k_id]); + + for(size_t dim = 0; dim < Dimension; ++dim){ + f_bar_min[dim] = std::min(f_bar_min[dim], f_xk[dim]); + f_bar_max[dim] = std::max(f_bar_max[dim], f_xk[dim]); + } + } + + const double eps = 1E-14; + Rd coef1; + for(size_t dim = 0; dim < Dimension; ++dim){ + coef1[dim] = 1; + if (std::abs(f_bar_max[dim] - fj[dim]) > eps) { + coef1[dim] = (f_max[dim] - fj[dim]) / ((f_bar_max[dim] - fj[dim])); + } + } + + Rd coef2; + for(size_t dim = 0; dim < Dimension; ++dim){ + coef2[dim] = 1; + if (std::abs(f_bar_min[dim] - fj[dim]) > eps) { + coef2[dim] = (f_min[dim] - fj[dim]) / ((f_bar_min[dim] - fj[dim])); + } + } + + double min_coef1 = coef1[0]; + double min_coef2 = coef2[0]; + for(size_t dim = 1; dim < Dimension; ++dim){ + min_coef1 = std::min(min_coef1,coef1[dim]); + min_coef2 = std::min(min_coef2,coef2[dim]); + } + + const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2))); + + auto coefficients = DPk_fh.coefficients(cell_id); + + coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0]; + for (size_t i = 1; i < coefficients.size(); ++i) { + coefficients[i] *= lambda; + } + }); + } + + CellValue<double> + _scalar_limiter_coefficient(const MeshType& mesh, + const DiscreteFunctionP0<const double>& f, + const DiscreteFunctionDPk<Dimension, const double>& DPk_fh) const + { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list; + StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes}; + auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list); + + const auto xj = mesh_data.xj(); + + CellValue<double> lambda{mesh.connectivity()}; + + parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){ + const double fj = f[cell_id]; + + double f_min = fj; + double f_max = fj; + + const auto cell_stencil = stencil[cell_id]; + for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){ + f_min = std::min(f_min, f[cell_stencil[i_cell]]); + f_max = std::max(f_max, f[cell_stencil[i_cell]]); + } + + double f_bar_min = fj; + double f_bar_max = fj; + + for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){ + const CellId cell_k_id = cell_stencil[i_cell]; + const double f_xk = DPk_fh[cell_id](xj[cell_k_id]); + + f_bar_min = std::min(f_bar_min, f_xk); + f_bar_max = std::max(f_bar_max, f_xk); + } + + const double eps = 1E-14; + double coef1 = 1; + if (std::abs(f_bar_max - fj) > eps) { + coef1 = (f_max - fj) / ((f_bar_max - fj)); + } + + double coef2 = 1.; + if (std::abs(f_bar_min - fj) > eps) { + coef2 = (f_min - fj) / ((f_bar_min - fj)); + } + + lambda[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2))); + }); + + return lambda; + } + + NodeValuePerCell<Rdxd> + _compute_tensor_reconstruction(const std::shared_ptr<const MeshVariant>& v_mesh, + DiscreteTensorFunction sigma, + PolynomialReconstructionDescriptor reconstruction_descriptor) const + { + const MeshType& mesh = *v_mesh->get<MeshType>(); + + CellValue<double> limiter_j{mesh.connectivity()}; + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j){ + limiter_j[j] = 1; + } + ); + + const CellValue<const Rdxd>& sigma_j = copy(sigma.cellValues()); + std::vector<DiscreteFunctionDPk<Dimension, double>> sigma_coef; + std::vector<size_t> row; + std::vector<size_t> col; + for(size_t i = 0; i < Dimension; ++i){ + for(size_t k = i; k < Dimension; ++k){ + + CellValue<double> coef{mesh.connectivity()}; + for(CellId j = 0; j <mesh.numberOfCells(); ++j){ + coef[j] = sigma_j[j](i,k); + } + + const auto& coef_scalar_function = DiscreteScalarFunction(v_mesh, coef); + auto coef_v = std::make_shared<DiscreteFunctionVariant>(coef_scalar_function); + auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(coef_v); + auto DPk_coef = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const double>>(); + + const CellValue<const double>& limiter_ik = _scalar_limiter_coefficient(mesh,coef_scalar_function,DPk_coef); + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + limiter_j[j] = std::min(limiter_j[j],limiter_ik[j]); + }); + + sigma_coef.push_back(copy(DPk_coef)); + row.push_back(i); + col.push_back(k); + } + } + + NodeValuePerCell<Rdxd> sigma_lim{mesh.connectivity()} ; + const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j){ + const auto& cell_nodes = cell_to_node_matrix[j]; + for(size_t i = 0; i < Dimension; ++i){ + for(size_t k = 0; k < Dimension; ++k){ + for(size_t R = 0; R < cell_nodes.size(); ++R){ + sigma_lim(j,R)(i,k) = 0; + } + } + } + } + ); + + const NodeValue<const Rd>& xr = copy(mesh.xr()); + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + for(size_t R = 0; R < cell_nodes.size(); ++R){ + const NodeId r = cell_nodes[R]; + + for(size_t l = 0; l < sigma_coef.size(); ++l){ + const size_t& i = row[l]; + const size_t& k = col[l]; + + auto coefficients = sigma_coef[l].coefficients(j); + + coefficients[0] = (1-limiter_j[j])*sigma_j[j](i,k) + limiter_j[j] * coefficients[0]; + for (size_t indice = 1; indice < coefficients.size(); ++indice) { + coefficients[indice] *= limiter_j[j]; + } + + sigma_lim(j,R)(i,k) = sigma_coef[l][j](xr[r]); + if(i != k){ + sigma_lim(j,R)(i,k) *= 2; + } + } + + sigma_lim(j,R) += transpose(sigma_lim(j,R)); + sigma_lim(j,R) *= 0.5; + } + + } + ); + + return sigma_lim; + } + public: std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>> compute_fluxes(const SolverType& solver_type, @@ -476,10 +698,28 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>(); const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>(); + std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list; + + for(auto&& bc_descriptor : bc_descriptor_list){ + if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){ + symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared()); + } + } + + PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1, + symmetry_boundary_descriptor_list); + + NodeValuePerCell<Rdxd> sigma_lim = _compute_tensor_reconstruction(i_mesh,sigma,reconstruction_descriptor); + auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v); + auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>(); + + DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh); + this->_vector_limiter(mesh,u,u_lim); + NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT); NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr); - NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, sigma); + NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, sigma_lim); const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list); this->_applyBoundaryConditions(bc_list, mesh, Ar, br); @@ -488,7 +728,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi synchronize(br); NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br); - NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma); + NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim); return std::make_tuple(std::make_shared<const ItemValueVariant>(ur), std::make_shared<const SubItemValuePerItemVariant>(Fjr)); @@ -547,13 +787,45 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi const DiscreteScalarFunction& aT2 = aT2_v->get<DiscreteScalarFunction>(); const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>(); + std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1; + std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2; + + for(auto&& bc_descriptor : bc_descriptor_list1){ + if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){ + symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared()); + } + } + for(auto&& bc_descriptor : bc_descriptor_list2){ + if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){ + symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared()); + } + } + + PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1, + symmetry_boundary_descriptor_list1); + PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1, + symmetry_boundary_descriptor_list2); + + NodeValuePerCell<Rdxd> sigma1_lim = _compute_tensor_reconstruction(i_mesh1,sigma1,reconstruction_descriptor1); + NodeValuePerCell<Rdxd> sigma2_lim = _compute_tensor_reconstruction(i_mesh2,sigma2,reconstruction_descriptor2); + + auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v); + auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>(); + auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v); + auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>(); + + DiscreteFunctionDPk<Dimension, Rd> u1_lim = copy(DPk_uh1); + this->_vector_limiter(mesh1,u1,u1_lim); + DiscreteFunctionDPk<Dimension, Rd> u2_lim = copy(DPk_uh2); + this->_vector_limiter(mesh2,u2,u2_lim); + NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1); NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2); NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1); - NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1, sigma1); + NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1_lim, sigma1_lim); NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2); - NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2, sigma2); + NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2_lim, sigma2_lim); this->_applyCouplingBC(Ar1, Ar2, br1, br2, map1, map2); const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1); @@ -567,11 +839,11 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi synchronize(br2); NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1); - NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1); + NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, sigma1_lim); NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2); - NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2); + NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim); NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2); - NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2); + NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim); return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1), std::make_shared<const SubItemValuePerItemVariant>(Fjr1), @@ -607,10 +879,28 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi const DiscreteScalarFunction& aT = aT_v->get<DiscreteScalarFunction>(); const DiscreteTensorFunction& sigma = sigma_v->get<DiscreteTensorFunction>(); + std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list; + + for(auto&& bc_descriptor : bc_descriptor_list){ + if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){ + symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared()); + } + } + + PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1, + symmetry_boundary_descriptor_list); + + NodeValuePerCell<Rdxd> sigma_lim = _compute_tensor_reconstruction(i_mesh,sigma,reconstruction_descriptor); + auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(u_v); + auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>(); + + DiscreteFunctionDPk<Dimension, Rd> u_lim = copy(DPk_uh); + this->_vector_limiter(mesh,u,u_lim); + NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT); NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr); - NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, sigma); + NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u_lim, sigma_lim); const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list); this->_applyBoundaryConditions(bc_list, mesh, Ar, br); @@ -619,7 +909,7 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi synchronize(br); NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br); - NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma); + NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, sigma_lim); this->_applyCouplingBC(mesh, ur, CR_ur, Fjr, CR_Fjr, map2); @@ -680,13 +970,45 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi const DiscreteScalarFunction& aT2 = aT2_v->get<DiscreteScalarFunction>(); const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>(); + std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list1; + std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list2; + + for(auto&& bc_descriptor : bc_descriptor_list1){ + if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){ + symmetry_boundary_descriptor_list1.push_back(bc_descriptor->boundaryDescriptor_shared()); + } + } + for(auto&& bc_descriptor : bc_descriptor_list2){ + if(bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry){ + symmetry_boundary_descriptor_list2.push_back(bc_descriptor->boundaryDescriptor_shared()); + } + } + + PolynomialReconstructionDescriptor reconstruction_descriptor1(IntegrationMethodType::cell_center, 1, + symmetry_boundary_descriptor_list1); + PolynomialReconstructionDescriptor reconstruction_descriptor2(IntegrationMethodType::cell_center, 1, + symmetry_boundary_descriptor_list2); + + NodeValuePerCell<Rdxd> sigma1_lim = _compute_tensor_reconstruction(i_mesh1,sigma1,reconstruction_descriptor1); + NodeValuePerCell<Rdxd> sigma2_lim = _compute_tensor_reconstruction(i_mesh2,sigma2,reconstruction_descriptor2); + + auto reconstruction1 = PolynomialReconstruction{reconstruction_descriptor1}.build(u1_v); + auto DPk_uh1 = reconstruction1[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>(); + auto reconstruction2 = PolynomialReconstruction{reconstruction_descriptor2}.build(u2_v); + auto DPk_uh2 = reconstruction2[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>(); + + DiscreteFunctionDPk<Dimension, Rd> u1_lim = copy(DPk_uh1); + this->_vector_limiter(mesh1,u1,u1_lim); + DiscreteFunctionDPk<Dimension, Rd> u2_lim = copy(DPk_uh2); + this->_vector_limiter(mesh2,u2,u2_lim); + NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1); NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2); NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1); - NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1, sigma1); + NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1_lim, sigma1_lim); NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2); - NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2, sigma2); + NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2_lim, sigma2_lim); const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1); const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2); @@ -701,8 +1023,8 @@ class Order2LocalDtHyperelasticSolverHandler::Order2LocalDtHyperelasticSolver fi NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1); NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2); this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2); - NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1); - NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2); + NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1_lim, sigma1_lim); + NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2_lim, sigma2_lim); return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1), std::make_shared<const SubItemValuePerItemVariant>(Fjr1),