From 111c1cf1801db16a6d086da316c97aae8b549642 Mon Sep 17 00:00:00 2001 From: npaster Date: Wed, 11 Dec 2024 19:29:24 +0100 Subject: [PATCH] Prepare elastodynamics --- apps/contact/src/common.hpp | 14 +- apps/contact/src/diffusion_nitsche_solver.hpp | 4 +- apps/contact/src/signorini_newton_solver.hpp | 16 +- apps/diffusion/src/diffusion_hho_test.cpp | 2 +- apps/maxwell/src/maxwell.cpp | 8 +- apps/maxwell/src/maxwell_sip_dg.cpp | 14 +- .../src/contact_test.cpp | 26 +- .../src/elastodynamic_test.cpp | 106 +---- .../src/nonlinear_solid_mechanics.cpp | 6 +- apps/tests/curl_reconstruction.cpp | 2 +- .../include/diskpp/bases/bases_scalar.hpp | 2 +- .../include/diskpp/bases/bases_utils.hpp | 9 +- .../include/diskpp/bases/bases_vector.hpp | 49 +- .../boundary_conditions.hpp | 144 +++--- .../diskpp/common/simplicial_formula.hpp | 107 +++-- .../diskpp/mechanics/NewtonSolver/Fields.hpp | 331 ++++++++++++++ .../NewtonSolver/NewtonIteration.hpp | 422 ++++++++---------- .../mechanics/NewtonSolver/NewtonSolver.hpp | 398 +++++++---------- .../NewtonSolver/NewtonSolverComput.hpp | 55 +-- .../NewtonSolver/NewtonSolverDynamic.hpp | 234 ++++++++++ .../NewtonSolver/NewtonSolverParameters.hpp | 37 +- .../mechanics/NewtonSolver/NewtonStep.hpp | 90 +--- .../mechanics/behaviors/laws/law_qp_bones.hpp | 15 +- .../LogarithmicStrain_qp.hpp | 4 +- .../methods/implementation_hho/curl.hpp | 36 +- .../vector_hho_laplacian.hpp | 66 ++- .../quadratures/bits/quad_raw_tetra.hpp | 23 +- 27 files changed, 1305 insertions(+), 915 deletions(-) create mode 100644 libdiskpp/include/diskpp/mechanics/NewtonSolver/Fields.hpp create mode 100644 libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverDynamic.hpp diff --git a/apps/contact/src/common.hpp b/apps/contact/src/common.hpp index 7f5263d6..8c8dde1d 100644 --- a/apps/contact/src/common.hpp +++ b/apps/contact/src/common.hpp @@ -1903,7 +1903,7 @@ class diffusion_condensed_assembler_nitsche_cells { auto fb = make_scalar_monomial_basis(msh, fc, di.face_degree()); auto dirichlet_bf = m_bnd.dirichlet_boundary_func(face_id); - Matrix mass = make_mass_matrix(msh, fc, fb, di.face_degree()); + Matrix mass = make_mass_matrix(msh, fc, fb); Matrix rhs = make_rhs(msh, fc, fb, dirichlet_bf, di.face_degree()); //ret.block(face_i*fbs, 0, fbs, 1) = mass.llt().solve(rhs); } @@ -2249,7 +2249,7 @@ class diffusion_full_assembler if (m_bnd.is_dirichlet_face( face_id)) { auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree()); - Matrix mass = make_mass_matrix(msh, fc, fb, di.face_degree()); + Matrix mass = make_mass_matrix(msh, fc, fb); auto velocity = m_bnd.dirichlet_boundary_func(face_id); Matrix rhs = make_rhs(msh, fc, fb, velocity, di.face_degree()); svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs); @@ -2569,7 +2569,7 @@ class diffusion_mix_full_assembler if (m_bnd.is_dirichlet_face( face_id)) { auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree()); - Matrix mass = make_mass_matrix(msh, fc, fb, di.face_degree()); + Matrix mass = make_mass_matrix(msh, fc, fb); auto velocity = m_bnd.dirichlet_boundary_func(face_id); Matrix rhs = make_rhs(msh, fc, fb, velocity, di.face_degree()); svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs); @@ -2893,7 +2893,7 @@ class contact_full_assembler if (dirichlet) { auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree()); - Matrix mass = make_mass_matrix(msh, fc, fb, di.face_degree()); + Matrix mass = make_mass_matrix(msh, fc, fb); auto velocity = m_bnd.dirichlet_boundary_func(face_id); Matrix rhs = make_rhs(msh, fc, fb, velocity);//, di.face_degree()); svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs); @@ -3119,7 +3119,7 @@ class contact_full_assembler_new auto fb = make_scalar_monomial_basis(msh, face, di.face_degree()); auto dirichlet_fun = m_bnd.dirichlet_boundary_func(face_id); - matrix_type mass = make_mass_matrix(msh, face, fb);// di.face_degree()); + matrix_type mass = make_mass_matrix(msh, face, fb); vector_type rhs = make_rhs(msh, face, fb, dirichlet_fun);// di.face_degree()); sol.block(face_ofs, 0, fbs, 1) = mass.llt().solve(rhs); @@ -3448,7 +3448,7 @@ class contact_face_assembler_new Matrix robin = make_rhs(msh,bfc,fb,bnd.robin_boundary_func(face_id), face_degree); assert (robin.size() == num_face_dofs); - Matrix mass = make_mass_matrix(msh, bfc, fb, face_degree); + Matrix mass = make_mass_matrix(msh, bfc, fb); for (size_t i = 0; i < num_face_dofs; i++) { @@ -3518,7 +3518,7 @@ class contact_face_assembler_new auto fb = make_scalar_monomial_basis(msh, face, di.face_degree()); auto dirichlet_fun = m_bnd.dirichlet_boundary_func(face_id); - Matrix mass = make_mass_matrix(msh, face, fb, di.face_degree()); + Matrix mass = make_mass_matrix(msh, face, fb); Matrix rhs = make_rhs(msh, face, fb, dirichlet_fun, di.face_degree()); sol.block(face_ofs, 0, fb.size(), 1) = mass.llt().solve(rhs); diff --git a/apps/contact/src/diffusion_nitsche_solver.hpp b/apps/contact/src/diffusion_nitsche_solver.hpp index ab1655bd..4ed5f730 100644 --- a/apps/contact/src/diffusion_nitsche_solver.hpp +++ b/apps/contact/src/diffusion_nitsche_solver.hpp @@ -208,7 +208,7 @@ run_hho_diffusion_nitsche_faces(const Mesh& msh, auto diff = realsol - fullsol; H1_error += diff.dot(A*diff); - matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = diff.block(0, 0, cbs, 1); L2_error += u_diff.dot(mass * u_diff); @@ -374,7 +374,7 @@ run_hho_diffusion_nitsche_cells_full(const Mesh& msh, auto diff = realsol - fullsol; H1_error += diff.dot(A*diff); - matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = diff.block(0, 0, cbs, 1); L2_error += u_diff.dot(mass * u_diff); diff --git a/apps/contact/src/signorini_newton_solver.hpp b/apps/contact/src/signorini_newton_solver.hpp index 36936eac..6184c736 100644 --- a/apps/contact/src/signorini_newton_solver.hpp +++ b/apps/contact/src/signorini_newton_solver.hpp @@ -157,7 +157,7 @@ solve_faces(const Mesh& msh, const Function& rhs_fun, const Analytical& sol_fun H1_increment += du_full.dot(A * du_full); - matrix_type mass = make_mass_matrix(msh, cl, cb);//, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = du_full.block(0, 0, cbs, 1); L2_increment += u_diff.dot(mass * u_diff); @@ -211,7 +211,7 @@ solve_faces(const Mesh& msh, const Function& rhs_fun, const Analytical& sol_fun H1_error += diff.dot(Ah*diff); auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree()); - matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = diff.block(0, 0, cbs, 1); L2_error += u_diff.dot(mass * u_diff); @@ -419,7 +419,7 @@ solve_faces_hier(const Mesh& msh, const Function& rhs_fun, const Analytical& so H1_increment += du_full.dot(A * du_full); - matrix_type mass = make_mass_matrix(msh, cl, cb);//, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = du_full.block(0, 0, cbs, 1); L2_increment += u_diff.dot(mass * u_diff); @@ -473,7 +473,7 @@ solve_faces_hier(const Mesh& msh, const Function& rhs_fun, const Analytical& so H1_error += diff.dot(Ah*diff); auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree()); - matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = diff.block(0, 0, cbs, 1); L2_error += u_diff.dot(mass * u_diff); @@ -643,7 +643,7 @@ solve_cells_full(const Mesh& msh, const Function& rhs_fun, const Analytical& so diff_sol.block(cell_ofs, 0, num_total_dofs ,1) = du_full; auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree()); - matrix_type mass = make_mass_matrix(msh, cl, cb);//, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = du_full.block(0, 0, cbs, 1); H1_increment += du_full.dot(Ah * du_full); @@ -721,7 +721,7 @@ solve_cells_full(const Mesh& msh, const Function& rhs_fun, const Analytical& so #endif auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree()); - matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = diff.block(0, 0, cbs, 1); L2_error += u_diff.dot(mass * u_diff); @@ -911,7 +911,7 @@ solve_cells_full_hier(const Mesh& msh, const Function& rhs_fun, const Analytica diff_sol.block(cell_ofs, 0, num_total_dofs ,1) = du_full; auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree()); - matrix_type mass = make_mass_matrix(msh, cl, cb);//, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = du_full.block(0, 0, cbs, 1); H1_increment += du_full.dot(Ah * du_full); @@ -966,7 +966,7 @@ solve_cells_full_hier(const Mesh& msh, const Function& rhs_fun, const Analytica H1_error += diff.dot(Ah*diff); auto cb = make_scalar_monomial_basis(msh, cl, hdi.cell_degree()); - matrix_type mass = make_mass_matrix(msh, cl, cb, hdi.cell_degree()); + matrix_type mass = make_mass_matrix(msh, cl, cb); vector_type u_diff = diff.block(0, 0, cbs, 1); L2_error += u_diff.dot(mass * u_diff); diff --git a/apps/diffusion/src/diffusion_hho_test.cpp b/apps/diffusion/src/diffusion_hho_test.cpp index 25da32b8..b44aa86e 100644 --- a/apps/diffusion/src/diffusion_hho_test.cpp +++ b/apps/diffusion/src/diffusion_hho_test.cpp @@ -191,7 +191,7 @@ run_hho_diffusion_solver(const Mesh& msh, const size_t degree) matrix_type ATT = A.block(0,0, cbs, cbs); vector_type u_cell = fullsol.block(0,0, cbs, 1); vector_type res = ATT * u_cell - rhs; - matrix_type mm = make_mass_matrix(msh, cl, cb, hdi.cell_degree()); + matrix_type mm = make_mass_matrix(msh, cl, cb); error += res.dot(mm * res); //auto diff = realsol - fullsol; diff --git a/apps/maxwell/src/maxwell.cpp b/apps/maxwell/src/maxwell.cpp index 6e358264..fec94d62 100644 --- a/apps/maxwell/src/maxwell.cpp +++ b/apps/maxwell/src/maxwell.cpp @@ -663,7 +663,7 @@ vector_wave_solver(Mesh& msh, size_t order, auto qps = integrate(msh, cl, 2*chdi.reconstruction_degree()); for (auto& qp : qps) { - Matrix hphi = rb.eval_curls2(qp.point()); + Matrix hphi = rb.eval_curls(qp.point()); Matrix hphi2 = hphi.block(0,0,cb.size(),2); Matrix ephi = cb.eval_functions(qp.point()); //Matrix rphi = rb.eval_functions(qp.point()); @@ -1156,7 +1156,7 @@ vector_wave_solver_complex(Mesh& msh, parameter_loader& msh, parameter_loader f_phi = fb.eval_functions(qp.point()); - Matrix c_cphi_tmp = cb.eval_curls2(qp.point()); + Matrix c_cphi_tmp = cb.eval_curls(qp.point()); Matrix c_cphi = disk::vcross(c_cphi_tmp, n); mass += qp.weight() * f_phi * f_phi.transpose(); @@ -1217,7 +1217,7 @@ vector_wave_solver_complex(Mesh& msh, parameter_loader cphi_tmp = cb.eval_functions(fc_pts[j]); Matrix n_x_cphi_x_n = disk::vcross(n, disk::vcross(cphi_tmp, n)); - Matrix ccphi_tmp = cb.eval_curls2(fc_pts[j]); + Matrix ccphi_tmp = cb.eval_curls(fc_pts[j]); Matrix ccphi_x_n = disk::vcross(ccphi_tmp, n); Matrix imp = (1./mur)*ccphi_x_n + (jwmu0/Z)*n_x_cphi_x_n; diff --git a/apps/maxwell/src/maxwell_sip_dg.cpp b/apps/maxwell/src/maxwell_sip_dg.cpp index abd4866b..5d98c37b 100644 --- a/apps/maxwell/src/maxwell_sip_dg.cpp +++ b/apps/maxwell/src/maxwell_sip_dg.cpp @@ -222,7 +222,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ for (auto& qp : qps) { auto phi = tbasis.eval_functions(qp.point()); - auto cphi = tbasis.eval_curls2(qp.point()); + auto cphi = tbasis.eval_curls(qp.point()); M += qp.weight() * phi * phi.transpose(); K += qp.weight() * cphi * cphi.transpose(); @@ -253,7 +253,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ for (auto& fqp : f_qps) { auto tphi = tbasis.eval_functions(fqp.point()); - auto tcphi = tbasis.eval_curls2(fqp.point()); + auto tcphi = tbasis.eval_curls(fqp.point()); auto n_x_tphi = disk::vcross(n, tphi); auto n_x_tphi_x_n = disk::vcross(n_x_tphi, n); @@ -287,7 +287,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ } auto nphi = nbasis.eval_functions(fqp.point()); - auto ncphi = nbasis.eval_curls2(fqp.point()); + auto ncphi = nbasis.eval_curls(fqp.point()); auto n_x_nphi = disk::vcross(n, nphi); Atn += - fqp.weight() * eta_l * n_x_tphi * n_x_nphi.transpose(); @@ -375,7 +375,7 @@ run_maxwell_solver(Mesh& msh, size_t degree, const typename Mesh::coordinate_typ { Matrix cphi_tmp = cb.eval_functions(fc_pts[j]); Matrix n_x_cphi_x_n = disk::vcross(n, disk::vcross(cphi_tmp, n)); - Matrix ccphi_tmp = cb.eval_curls2(fc_pts[j]); + Matrix ccphi_tmp = cb.eval_curls(fc_pts[j]); Matrix ccphi_x_n = disk::vcross(ccphi_tmp, n); Matrix imp = (1./mur)*ccphi_x_n - (jwmu0/Z)*n_x_cphi_x_n; @@ -513,7 +513,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo for (auto& qp : qps) { auto phi = tbasis.eval_functions(qp.point()); - auto cphi = tbasis.eval_curls2(qp.point()); + auto cphi = tbasis.eval_curls(qp.point()); M += qp.weight() * phi * phi.transpose(); K += qp.weight() * cphi * cphi.transpose(); @@ -540,7 +540,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo for (auto& fqp : f_qps) { auto tphi = tbasis.eval_functions(fqp.point()); - auto tcphi = tbasis.eval_curls2(fqp.point()); + auto tcphi = tbasis.eval_curls(fqp.point()); auto n_x_tphi = disk::vcross(n, tphi); if (nv.second) @@ -558,7 +558,7 @@ run_maxwell_eigenvalue_solver(Mesh& msh, size_t degree, const typename Mesh::coo } auto nphi = nbasis.eval_functions(fqp.point()); - auto ncphi = nbasis.eval_curls2(fqp.point()); + auto ncphi = nbasis.eval_curls(fqp.point()); auto n_x_nphi = disk::vcross(n, nphi); Atn += - fqp.weight() * eta_l * n_x_tphi * n_x_nphi.transpose(); diff --git a/apps/nonlinear_solid_mechanics/src/contact_test.cpp b/apps/nonlinear_solid_mechanics/src/contact_test.cpp index 88792b2d..9cf9d5a5 100644 --- a/apps/nonlinear_solid_mechanics/src/contact_test.cpp +++ b/apps/nonlinear_solid_mechanics/src/contact_test.cpp @@ -204,16 +204,19 @@ run_tresca_solver(Mesh& msh, return -time * result_type{fx, fy} * exy / (3.0 * (1.0 + material_data.getLambda())); }; - auto solution = [material_data](const disk::point& p) -> result_type + auto solution = [material_data](const disk::point& p, const T time) -> result_type { T y = p.y(); T x = p.x(); T exy = std::exp(x * y); T coeff = 1.0 / (1.0 + material_data.getLambda()); - return result_type{x * exy * (1.0 + coeff), y * exy * (-1.0 + coeff)} / 6.0; + return time*result_type{x * exy * (1.0 + coeff), y * exy * (-1.0 + coeff)} / 6.0; }; + auto sol = [solution](const disk::point& p) -> auto + { return solution(p, 1.0); }; + auto s = [rp, material_data](const disk::point& p) -> T { T y = p.y(); @@ -236,7 +239,7 @@ run_tresca_solver(Mesh& msh, nl.addBehavior(disk::DeformationMeasure::SMALL_DEF, disk::LawType::ELASTIC); nl.addMaterialData(material_data); - nl.initial_guess(solution); + nl.initial_guess(sol); if (nl.verbose()) { @@ -254,8 +257,8 @@ run_tresca_solver(Mesh& msh, error.h = average_diameter(msh); error.degree = rp.m_face_degree; error.nb_dof = nl.numberOfDofs(); - error.error_L2 = nl.compute_l2_displacement_error(solution); - error.error_H1 = nl.compute_H1_error(solution); + error.error_L2 = nl.compute_l2_displacement_error(sol); + error.error_H1 = nl.compute_H1_error(sol); // nl.compute_stress_GP("stress2D_GP_test.msh"); // nl.compute_continuous_displacement("depl2D_cont_test.msh"); @@ -294,16 +297,19 @@ run_tresca_solver(const Mesh& msh, return -time * result_type{fx, 0.0, fz} * exz / (3.0 * (1.0 + material_data.getLambda())); }; - auto solution = [material_data](const disk::point& p) -> result_type + auto solution = [material_data](const disk::point& p, const T& time) -> result_type { T z = p.z(); T x = p.x(); T exz = std::exp(x * z); T coeff = 1.0 / (1.0 + material_data.getLambda()); - return result_type{x * exz * (1.0 + coeff), 0.0, z * exz * (-1.0 + coeff)} / 6.0; + return time*result_type{x * exz * (1.0 + coeff), 0.0, z * exz * (-1.0 + coeff)} / 6.0; }; + auto sol = [solution](const disk::point& p) -> auto + { return solution(p, 1.0); }; + auto s = [rp, material_data](const disk::point& p) -> T { T x = p.x(); @@ -322,7 +328,7 @@ run_tresca_solver(const Mesh& msh, nl.addBehavior(disk::DeformationMeasure::SMALL_DEF, disk::LawType::ELASTIC); nl.addMaterialData(material_data); - nl.initial_guess(solution); + nl.initial_guess(sol); SolverInfo solve_info = nl.compute(load); @@ -335,8 +341,8 @@ run_tresca_solver(const Mesh& msh, error.h = average_diameter(msh); error.degree = rp.m_face_degree; error.nb_dof = nl.numberOfDofs(); - error.error_L2 = nl.compute_l2_displacement_error(solution); - error.error_H1 = nl.compute_H1_error(solution); + error.error_L2 = nl.compute_l2_displacement_error(sol); + error.error_H1 = nl.compute_H1_error(sol); return error; } diff --git a/apps/nonlinear_solid_mechanics/src/elastodynamic_test.cpp b/apps/nonlinear_solid_mechanics/src/elastodynamic_test.cpp index 6ccb0066..2bdb04df 100644 --- a/apps/nonlinear_solid_mechanics/src/elastodynamic_test.cpp +++ b/apps/nonlinear_solid_mechanics/src/elastodynamic_test.cpp @@ -90,60 +90,23 @@ run_linear_elasticity_solver(const Mesh& msh, return -M_PI * M_PI / (lambda + 1) * result_type{fx, fy}; }; - auto displacement = [material_data](const disk::point& p) -> result_type + auto displacement = [material_data](const disk::point& p, const T& time) -> result_type { - T time = 1.0; T fx = sin(2 * M_PI * p.y()) * (cos(2 * M_PI * p.x()) - 1) + 1.0 / (1 + material_data.getLambda()) * sin(M_PI * p.x()) * sin(M_PI * p.y()); T fy = -sin(2 * M_PI * p.x()) * (cos(2 * M_PI * p.y()) - 1) + 1.0 / (1 + material_data.getLambda()) * sin(M_PI * p.x()) * sin(M_PI * p.y()); - return result_type{fx, fy}; + return time * result_type{fx, fy}; }; - auto velocity = [material_data](const disk::point& p) -> result_type + auto velocity = [material_data, displacement](const disk::point& p, const T& time) -> result_type { - T fx = sin(2 * M_PI * p.y()) * (cos(2 * M_PI * p.x()) - 1) + - 1.0 / (1 + material_data.getLambda()) * sin(M_PI * p.x()) * sin(M_PI * p.y()); - T fy = -sin(2 * M_PI * p.x()) * (cos(2 * M_PI * p.y()) - 1) + - 1.0 / (1 + material_data.getLambda()) * sin(M_PI * p.x()) * sin(M_PI * p.y()); - - return result_type{fx, fy}; + return displacement(p, 1.0); }; - auto acceleration = [material_data](const disk::point& p) -> result_type { return result_type{0., 0.}; }; - - auto sigma = [material_data](const disk::point& p) -> grad_type - { - const T lambda = material_data.getLambda(); - const T mu = material_data.getMu(); - - T g11 = - -(2 * (lambda + 1) * sin(2 * M_PI * p.x()) * sin(2 * M_PI * p.y()) - sin(M_PI * p.y()) * cos(M_PI * p.x())); - T g12 = (2 * lambda + 2) * (cos(2 * M_PI * p.x()) - 1) * cos(2 * M_PI * p.y()) + - sin(M_PI * p.x()) * cos(M_PI * p.y()); - - T g21 = (-2 * lambda + 2) * (cos(2 * M_PI * p.y()) - 1) * cos(2 * M_PI * p.x()) + - sin(M_PI * p.y()) * cos(M_PI * p.x()); - - T g22 = - 2 * (lambda + 1) * sin(2 * M_PI * p.x()) * sin(2 * M_PI * p.y()) + sin(M_PI * p.x()) * cos(M_PI * p.y()); - - grad_type g = grad_type::Zero(); - - g(0, 0) = g11; - g(0, 1) = g12; - g(1, 0) = g21; - g(1, 1) = g22; - - g *= M_PI / (lambda + 1); - - const grad_type gs = 0.5 * (g + g.transpose()); - - const T divu = gs.trace(); - - return 2 * mu * gs + lambda * divu * disk::static_matrix::Identity(); - }; + auto acceleration = [material_data, displacement](const disk::point &p, const T &time) -> result_type + { return displacement(p, 0.0); }; Bnd_type bnd(msh); bnd.addDirichletEverywhere(displacement); @@ -153,7 +116,12 @@ run_linear_elasticity_solver(const Mesh& msh, nl.addBehavior(disk::DeformationMeasure::SMALL_DEF, disk::LawType::ELASTIC); nl.addMaterialData(material_data); - nl.initial_fields(displacement, velocity, acceleration); + nl.initial_field(disk::mechanics::FieldName::DEPL, [displacement](const disk::point &p) + { return displacement(p, 0.0); }); + nl.initial_field(disk::mechanics::FieldName::VITE_CELLS, [velocity](const disk::point &p) + { return velocity(p, 0.0); }); + nl.initial_field(disk::mechanics::FieldName::ACCE_CELLS, [acceleration](const disk::point &p) + { return acceleration(p, 0.0); }); if (nl.verbose()) { @@ -171,8 +139,10 @@ run_linear_elasticity_solver(const Mesh& msh, error.h = average_diameter(msh); error.degree = rp.m_face_degree; error.nb_dof = nl.numberOfDofs(); - error.error_L2 = nl.compute_l2_displacement_error(displacement); - error.error_H1 = nl.compute_H1_error(displacement); + error.error_L2 = nl.compute_l2_displacement_error([displacement](const disk::point &p) + { return displacement(p, 1.0); }); + error.error_H1 = nl.compute_H1_error([displacement](const disk::point &p) + { return displacement(p, 1.0); }); return error; } @@ -212,46 +182,14 @@ run_linear_elasticity_solver(const Mesh& msh, return -M_PI * M_PI / (lambda + 1) * result_type{fx, fy, 0}; }; - auto displacement = [material_data](const disk::point& p) -> result_type + auto displacement = [material_data](const disk::point& p, const T& time) -> result_type { T fx = sin(2 * M_PI * p.y()) * (cos(2 * M_PI * p.x()) - 1) + 1.0 / (1 + material_data.getLambda()) * sin(M_PI * p.x()) * sin(M_PI * p.y()); T fy = -sin(2 * M_PI * p.x()) * (cos(2 * M_PI * p.y()) - 1) + 1.0 / (1 + material_data.getLambda()) * sin(M_PI * p.x()) * sin(M_PI * p.y()); - return result_type{fx, fy, 0}; - }; - - auto sigma = [material_data](const disk::point& p) -> grad_type - { - const T lambda = material_data.getLambda(); - const T mu = material_data.getMu(); - - T g11 = - -(2 * (lambda + 1) * sin(2 * M_PI * p.x()) * sin(2 * M_PI * p.y()) - sin(M_PI * p.y()) * cos(M_PI * p.x())); - T g12 = (2 * lambda + 2) * (cos(2 * M_PI * p.x()) - 1) * cos(2 * M_PI * p.y()) + - sin(M_PI * p.x()) * cos(M_PI * p.y()); - - T g21 = (-2 * lambda + 2) * (cos(2 * M_PI * p.y()) - 1) * cos(2 * M_PI * p.x()) + - sin(M_PI * p.y()) * cos(M_PI * p.x()); - - T g22 = - 2 * (lambda + 1) * sin(2 * M_PI * p.x()) * sin(2 * M_PI * p.y()) + sin(M_PI * p.x()) * cos(M_PI * p.y()); - - grad_type g = grad_type::Zero(); - - g(0, 0) = g11; - g(0, 1) = g12; - g(1, 0) = g21; - g(1, 1) = g22; - - g *= M_PI / (lambda + 1); - - const grad_type gs = 0.5 * (g + g.transpose()); - - const T divu = gs.trace(); - - return 2 * mu * gs + lambda * divu * disk::static_matrix::Identity(); + return time * result_type{fx, fy, 0}; }; Bnd_type bnd(msh); @@ -262,7 +200,8 @@ run_linear_elasticity_solver(const Mesh& msh, nl.addBehavior(disk::DeformationMeasure::SMALL_DEF, disk::LawType::ELASTIC); nl.addMaterialData(material_data); - nl.initial_guess(displacement); + nl.initial_field(disk::mechanics::FieldName::DEPL, [displacement](const disk::point &p) + { return displacement(p, 0.0); }); if (nl.verbose()) { @@ -280,8 +219,9 @@ run_linear_elasticity_solver(const Mesh& msh, error.h = average_diameter(msh); error.degree = rp.m_face_degree; error.nb_dof = nl.numberOfDofs(); - error.error_L2 = nl.compute_l2_displacement_error(displacement); - error.error_H1 = nl.compute_H1_error(displacement); + error.error_L2 = nl.compute_l2_displacement_error([displacement](const disk::point &p) + { return displacement(p, 1.0); }); + error.error_H1 = nl.compute_H1_error([displacement](const disk::point &p) { return displacement(p, 1.0); } ); return error; } diff --git a/apps/nonlinear_solid_mechanics/src/nonlinear_solid_mechanics.cpp b/apps/nonlinear_solid_mechanics/src/nonlinear_solid_mechanics.cpp index ce240262..24866b90 100644 --- a/apps/nonlinear_solid_mechanics/src/nonlinear_solid_mechanics.cpp +++ b/apps/nonlinear_solid_mechanics/src/nonlinear_solid_mechanics.cpp @@ -58,7 +58,7 @@ run_nl_solid_mechanics_solver(const Mesh& msh, // Cook with quadrilaterals auto zero = [material_data](const disk::point& p) -> result_type { return result_type{0.0, 0}; }; - auto trac = [material_data](const disk::point& p) -> result_type { return result_type{0.0, 0.3125}; }; + auto trac = [material_data](const disk::point& p, const T& time) -> result_type { return time*result_type{0.0, 0.3125}; }; bnd.addDirichletBC(disk::CLAMPED, 3, zero); bnd.addNeumannBC(disk::NEUMANN, 8, trac); @@ -132,7 +132,7 @@ run_nl_solid_mechanics_solver(const Mesh& msh, return 3 * er; }; - auto deplr = [material_data](const disk::point& p) -> result_type + auto deplr = [material_data](const disk::point& p, const T& time) -> result_type { result_type er = result_type::Zero(); @@ -142,7 +142,7 @@ run_nl_solid_mechanics_solver(const Mesh& msh, er /= er.norm(); - return 0.157 * er; + return time*0.157 * er; }; // Sphere Hpp diff --git a/apps/tests/curl_reconstruction.cpp b/apps/tests/curl_reconstruction.cpp index 9370e712..31044782 100644 --- a/apps/tests/curl_reconstruction.cpp +++ b/apps/tests/curl_reconstruction.cpp @@ -142,7 +142,7 @@ struct test_functor_curl_reconstruction, mixed> auto qps = integrate(msh, cl, 2*rd); for (auto& qp : qps) { - auto phi = rb.eval_curls2(qp.point()); + auto phi = rb.eval_curls(qp.point()); Matrix val = disk::eval(reconstr_curl, phi) - curl_f(qp.point()); error += qp.weight() * val.dot(val); } diff --git a/libdiskpp/include/diskpp/bases/bases_scalar.hpp b/libdiskpp/include/diskpp/bases/bases_scalar.hpp index be7a5936..d7605935 100644 --- a/libdiskpp/include/diskpp/bases/bases_scalar.hpp +++ b/libdiskpp/include/diskpp/bases/bases_scalar.hpp @@ -286,7 +286,7 @@ class scaled_monomial_scalar_basis, typename Mesh>& a, const Matrix template Matrix -make_mass_matrix(const Mesh& msh, const Element& elem, const Basis& basis, size_t di = 0) +make_mass_matrix(const Mesh& msh, const Element& elem, const Basis& basis) { const auto degree = basis.degree(); const auto basis_size = basis.size(); @@ -216,7 +216,7 @@ make_mass_matrix(const Mesh& msh, const Element& elem, const Basis& basis, size_ Matrix ret = Matrix::Zero(basis_size, basis_size); - const auto qps = integrate(msh, elem, 2 * (degree+di)); + const auto qps = integrate(msh, elem, 2 * degree); for (auto& qp : qps) { @@ -228,11 +228,7 @@ make_mass_matrix(const Mesh& msh, const Element& elem, const Basis& basis, size_ return ret; } -/* We have a problem here: this definition could be ambiguous with the previous - * one. */ -//#if 0 template -[[deprecated("DiSk++ issue: this declaration is ambiguous, fix is needed")]] Matrix make_mass_matrix(const Mesh& msh, const Element& elem, const Basis& basis, const MaterialField& material_tensor, size_t di = 0) { @@ -255,7 +251,6 @@ make_mass_matrix(const Mesh& msh, const Element& elem, const Basis& basis, const return ret; } -//#endif template Matrix diff --git a/libdiskpp/include/diskpp/bases/bases_vector.hpp b/libdiskpp/include/diskpp/bases/bases_vector.hpp index ac47f94e..4c539b39 100644 --- a/libdiskpp/include/diskpp/bases/bases_vector.hpp +++ b/libdiskpp/include/diskpp/bases/bases_vector.hpp @@ -196,37 +196,8 @@ class scaled_monomial_vector_basis, typename Mesh dphi_i = dphi.row(i); - // row 1 - ret(j, 0) = dphi_i(1); - ret(j, 1) = dphi_i(2); - j++; - // row 2 - ret(j, 0) = -dphi_i(0); - ret(j, 2) = dphi_i(2); - j++; - // row 2 - ret(j, 1) = -dphi_i(0); - ret(j, 2) = -dphi_i(1); - j++; - } - assert(j == basis_size); - return ret; - } - - function_type - eval_curls2(const point_type& pt) const + eval_curls(const point_type &pt) const { function_type ret = function_type::Zero(basis_size, 3); @@ -687,24 +658,6 @@ class scaled_monomial_vector_basis, typename Mesh - eval_curls(const point_type& pt) const - { - Matrix ret = Matrix::Zero(basis_size); - - const function_type dphi = scalar_basis.eval_gradients(pt); - - size_t j = 0; - for (size_t i = 0; i < scalar_basis.size(); i++) - { - Matrix dphi_i = dphi.row(i); - - ret(j++) = dphi_i(1); - ret(j++) = -dphi_i(0); - } - return ret; - } - divergence_type eval_divergences(const point_type& pt) const { diff --git a/libdiskpp/include/diskpp/boundary_conditions/boundary_conditions.hpp b/libdiskpp/include/diskpp/boundary_conditions/boundary_conditions.hpp index c2fc2ffd..bae52946 100644 --- a/libdiskpp/include/diskpp/boundary_conditions/boundary_conditions.hpp +++ b/libdiskpp/include/diskpp/boundary_conditions/boundary_conditions.hpp @@ -69,26 +69,6 @@ enum ContactType : size_t namespace priv { -template -T -bnd_product(const T& fact, const T& func) -{ - return fact * func; -} - -template -Matrix -bnd_product(const T& fact, const Matrix& func) -{ - return fact * func; -} - -template -Matrix -bnd_product(const T& fact, const Matrix& func) -{ - return fact * func; -} template struct FunctionType @@ -186,15 +166,32 @@ class BoundaryConditions typedef typename mesh_type::face_type face_type; typedef typename mesh_type::coordinate_type scalar_type; typedef point point_type; - typedef typename priv::FunctionType::function_type function_type; + typedef typename priv::FunctionType::function_type fct_result_type; + typedef typename std::function function_type; - private: +private: const mesh_type& m_msh; - std::vector> m_dirichlet_func; - std::vector> m_neumann_func; - std::vector> m_robin_func; - std::vector> m_contact_func; + std::vector m_dirichlet_func; + std::vector m_neumann_func; + std::vector m_robin_func; + std::vector> m_contact_func; + + function_type + _conv_fct(const function_type &fct) + { + return fct; + } + + function_type + _conv_fct(const std::function fct) + { + + auto rfunc = [fct](const point_type &p, const scalar_type &time) -> auto + { return fct(p); }; + + return rfunc; + } // 1) bool to know if a boundary condition is associated to the face // 2) type of boundary conditions @@ -208,7 +205,7 @@ class BoundaryConditions size_t m_dirichlet_faces, m_neumann_faces, m_robin_faces, m_contact_faces; - scalar_type m_factor; + scalar_type m_time; // search faces that have the boundary id "b_id" std::vector @@ -236,8 +233,8 @@ class BoundaryConditions public: BoundaryConditions() = delete; - BoundaryConditions(const mesh_type& msh) : - m_msh(msh), m_dirichlet_faces(0), m_neumann_faces(0), m_robin_faces(0), m_contact_faces(0), m_factor(1) + BoundaryConditions(const mesh_type &msh) : m_msh(msh), m_dirichlet_faces(0), m_neumann_faces(0), m_robin_faces(0), m_contact_faces(0), + m_time(1.) { m_faces_is_dirichlet.assign(m_msh.faces_size(), std::make_tuple(false, NOTHING, 0, 0)); m_faces_is_neumann.assign(m_msh.faces_size(), std::make_tuple(false, FREE, 0, 0)); @@ -279,7 +276,7 @@ class BoundaryConditions addDirichletEverywhere(const Function& bcf) { const size_t bcf_id = m_dirichlet_func.size(); - m_dirichlet_func.push_back(bcf); + m_dirichlet_func.push_back(_conv_fct(bcf)); for (auto itor = m_msh.boundary_faces_begin(); itor != m_msh.boundary_faces_end(); itor++) { @@ -314,7 +311,7 @@ class BoundaryConditions addDirichletBC(const size_t& btype, const size_t& b_id, const Function& bcf) { const size_t bcf_id = m_dirichlet_func.size(); - m_dirichlet_func.push_back(bcf); + m_dirichlet_func.push_back(_conv_fct(bcf)); const auto list_faces = search_faces(b_id); @@ -341,12 +338,6 @@ class BoundaryConditions } } - void - multiplyAllFunctionsByAFactor(const scalar_type& factor) - { - m_factor = factor; - } - size_t nb_faces_boundary() const { @@ -662,44 +653,82 @@ class BoundaryConditions auto dirichlet_boundary_func(const size_t face_i) const + { + return dirichlet_boundary_func(face_i, m_time); + } + + auto + dirichlet_boundary_func(const face_type &fc) const + { + return dirichlet_boundary_func(m_msh.lookup(fc)); + } + + auto + dirichlet_boundary_func(const size_t face_i, const scalar_type &time) const { if (!is_dirichlet_face(face_i)) { throw std::logic_error("You want the Dirichlet function of face which is not a Dirichlet face"); } - const auto func = m_dirichlet_func.at(std::get<3>(m_faces_is_dirichlet.at(face_i))); - const scalar_type factor = m_factor; + const auto func = m_dirichlet_func.at(std::get<3>(m_faces_is_dirichlet.at(face_i))); - auto rfunc = [ func, factor ](const point_type& p) -> auto { return priv::bnd_product(factor, func(p)); }; + auto rfunc = [func, time](const point_type &p, scalar_type t = -123456789) -> auto + { + if (std::abs(t + 123456789.) > 1e-6) + { + return func(p, t); + } + + return func(p, time); + }; return rfunc; } auto - dirichlet_boundary_func(const face_type& fc) const + dirichlet_boundary_func(const face_type &fc, const scalar_type &time) const { - return dirichlet_boundary_func(m_msh.lookup(fc)); + return dirichlet_boundary_func(m_msh.lookup(fc), time); } auto - neumann_boundary_func(const size_t face_i) const + neumann_boundary_func(const size_t face_i, const scalar_type &time) const { if (!is_neumann_face(face_i)) { throw std::logic_error("You want the Neumann function of face which is not a Neumann face"); } - const auto func = m_neumann_func.at(std::get<3>(m_faces_is_neumann.at(face_i))); - const scalar_type factor = m_factor; + const auto func = m_neumann_func.at(std::get<3>(m_faces_is_neumann.at(face_i))); - auto rfunc = [ func, factor ](const point_type& p) -> auto { return priv::bnd_product(factor, func(p)); }; + auto rfunc = [func, time](const point_type &p, scalar_type t = -123456789) -> auto + { + if (std::abs(t + 123456789.) > 1e-6) + { + return func(p, t); + } + + return func(p, time); + }; return rfunc; } auto - neumann_boundary_func(const face_type& fc) const + neumann_boundary_func(const size_t face_i) const + { + return neumann_boundary_func(face_i, m_time); + } + + auto + neumann_boundary_func(const face_type &fc, const scalar_type &time) const + { + return neumann_boundary_func(m_msh.lookup(fc), time); + } + + auto + neumann_boundary_func(const face_type &fc) const { return neumann_boundary_func(m_msh.lookup(fc)); } @@ -713,12 +742,18 @@ class BoundaryConditions } const auto func = m_robin_func.at(std::get<3>(m_faces_is_robin.at(face_i))); - const scalar_type factor = m_factor; + const scalar_type time = m_time; - auto rfunc = [ func, factor ](const point_type& p) -> auto + auto rfunc = [func, time](const point_type &p, scalar_type t = -123456789) -> auto { - return disk::priv::inner_product(factor, func(p)); + if (std::abs(t + 123456789.) > 1e-6) + { + return func(p, t); + } + + return func(p, time); }; + return rfunc; } @@ -778,6 +813,15 @@ class BoundaryConditions return 0; } + + void setTime(const scalar_type &time) + { + m_time = time; + } + + scalar_type getTime(void) const { + return m_time; + } }; template diff --git a/libdiskpp/include/diskpp/common/simplicial_formula.hpp b/libdiskpp/include/diskpp/common/simplicial_formula.hpp index c46a19ac..93a0eed0 100644 --- a/libdiskpp/include/diskpp/common/simplicial_formula.hpp +++ b/libdiskpp/include/diskpp/common/simplicial_formula.hpp @@ -34,33 +34,39 @@ namespace disk { -/** - * @brief Compute the area of a triangle by using the Kahan formula, - * which minimize the round-off error - * - * @param p0 first point of the triangle - * @param p1 second point of the triangle - * @param p2 third point of the triangle - * @return T area - */ -template -T -area_triangle_kahan(const point& p0, const point& p1, const point& p2) -{ - const T l10 = (p1 - p0).to_vector().norm(); - const T l20 = (p2 - p0).to_vector().norm(); - const T l21 = (p2 - p1).to_vector().norm(); + /** + * @brief Compute the area of a triangle by using the Kahan formula, + * which minimize the round-off error + *https://inria.hal.science/hal-00790071/document + * + * @param p0 first point of the triangle + * @param p1 second point of the triangle + * @param p2 third point of the triangle + * @return T area + */ + template + T area_triangle_kahan(const point &p0, const point &p1, const point &p2) + { + const T l10 = (p1 - p0).to_vector().norm(); + const T l20 = (p2 - p0).to_vector().norm(); + const T l21 = (p2 - p1).to_vector().norm(); - std::array length = {l10, l20, l21}; + std::array length = {l10, l20, l21}; - std::sort(length.begin(), length.end()); + std::sort(length.begin(), length.end()); - const T a = length[2]; - const T b = length[1]; - const T c = length[0]; + const T a = length[2]; + const T b = length[1]; + const T c = length[0]; - return T(0.25) * std::sqrt((a + (b + c)) * (c - (a - b)) * (c + (a - b)) * (a + (b - c))); -} + // Minimal assumption on geometry + if (a > (b + c)) + { + throw std::runtime_error("Area is negative"); + } + + return T(0.25) * std::sqrt((((a + (b + c)) * (a + (b - c))) * (c + (a - b))) * (c - (a - b))); + } /** * @brief Compute the integration basis in order to map a point from the reference to physcal frame (for a triangle) @@ -107,11 +113,11 @@ integration_basis(const point& p0, const point& p1, const point T volume_tetrahedron_kahan(const point& p0, const point& p1, const point& p2, const point& p3) { - const auto v0 = (p1 - p0).to_vector(); - const auto v1 = (p2 - p0).to_vector(); - const auto v2 = (p3 - p0).to_vector(); - - return std::abs(v0.dot(v1.cross(v2))) / T(6); + // facial difference = (u-v+w) + auto fd = [](const T &u, const T &v, const T &w) + { return (std::max(u, w) - v) + std::min(u, w); }; + + // lengths of the edges + const auto l10 = (p1 - p0).to_vector().norm(); + const auto l20 = (p2 - p0).to_vector().norm(); + const auto l30 = (p3 - p0).to_vector().norm(); + + const auto l32 = (p3 - p2).to_vector().norm(); + const auto l31 = (p3 - p1).to_vector().norm(); + const auto l21 = (p2 - p1).to_vector().norm(); + + // sort edges + std::array, 3> length = { + std::make_pair(l10, l32), + std::make_pair(l20, l31), + std::make_pair(l30, l21)}; + + std::sort(length.begin(), length.end(), [&](const std::pair &va, const std::pair &vb) + { return (va.first + va.second < vb.first + vb.second); }); + + const auto [u, U] = length[2]; + const auto [v, V] = length[1]; + const auto [w, W] = length[0]; + + // Accurate products + const auto X = fd(w, U, v) * (U + v + w); + const auto Y = fd(u, V, w) * (V + w + u); + const auto Z = fd(v, W, u) * (W + u + v); + + const auto x = fd(U, v, w) * fd(v, w, U); + const auto y = fd(V, w, u) * fd(w, u, V); + const auto z = fd(W, u, v) * fd(u, v, W); + + // elementary factors + const auto xi = std::sqrt(x * Y * Z); + const auto eta = std::sqrt(y * Z * X); + const auto zeta = std::sqrt(z * X * Y); + const auto lambda = std::sqrt(x * y * z); + + const auto det = T(192.) * u * v * w; + + return std::sqrt((xi + eta + zeta - lambda) * (lambda + xi + eta - zeta) * (eta + zeta + lambda - xi) * (zeta + lambda + xi - eta)) / det; } /** diff --git a/libdiskpp/include/diskpp/mechanics/NewtonSolver/Fields.hpp b/libdiskpp/include/diskpp/mechanics/NewtonSolver/Fields.hpp new file mode 100644 index 00000000..5c926430 --- /dev/null +++ b/libdiskpp/include/diskpp/mechanics/NewtonSolver/Fields.hpp @@ -0,0 +1,331 @@ +/* + * /\ Matteo Cicuttin (C) 2016, 2017, 2018 + * /__\ matteo.cicuttin@enpc.fr + * /_\/_\ École Nationale des Ponts et Chaussées - CERMICS + * /\ /\ + * /__\ /__\ DISK++, a template library for DIscontinuous SKeletal + * /_\/_\/_\/_\ methods. + * + * This file is copyright of the following authors: + * Nicolas Pignet (C) 2024 + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * If you use this code or parts of it for scientific publications, you + * are required to cite it as following: + * + * Hybrid High-Order methods for finite elastoplastic deformations + * within a logarithmic strain framework. + * M. Abbas, A. Ern, N. Pignet. + * International Journal of Numerical Methods in Engineering (2019) + * 120(3), 303-327 + * DOI: 10.1002/nme.6137 + */ + +#pragma once + +#include + +#include "diskpp/adaptivity/adaptivity.hpp" + +namespace disk +{ + + namespace mechanics + { + + enum FieldName + { + DEPL, + DEPL_CELLS, + DEPL_FACES, + VITE_CELLS, + ACCE_CELLS, + }; + + std::string + getFieldName(const FieldName &name) + { + switch (name) + { + case FieldName::DEPL: + return "DEPL"; + break; + case FieldName::DEPL_CELLS: + return "DEPL_CELLS"; + break; + case FieldName::DEPL_FACES: + return "DEPL_FACES"; + break; + case FieldName::VITE_CELLS: + return "VITE_CELLS"; + break; + case FieldName::ACCE_CELLS: + return "ACCE_CELLS"; + break; + default: + break; + } + + return "Unknown type"; + } + + /** + * @brief Field at one time. + * + * @tparam scalar_type + */ + template + class TimeField + { + typedef dynamic_vector vector_type; + + std::map> m_fields; + + scalar_type m_time; + + public: + TimeField() : m_time(0) {} + + /** + * @brief Get the current time + * + */ + scalar_type + getTime(void) const + { + return m_time; + } + + /** + * @brief Set the current time + * + */ + void + setTime(const scalar_type &time) + { + m_time = time; + } + + void + setField(const FieldName &name, const std::vector &field) + { + m_fields[name] = field; + } + + auto + getField(const FieldName &name) const + { + // std::cout << getFieldName(name) << std::endl; + return m_fields.at(name); + } + + void + clearField(const FieldName &name) + { + if (auto search = m_fields.find(name); search != m_fields.end()) + { + m_fields.erase(search); + } + } + + template + void createZeroField(const FieldName &name, const Mesh &mesh, const MeshDegreeInfo °ree_infos) + { + std::vector field; + + if (name == DEPL || name == DEPL_CELLS || name == VITE_CELLS || name == ACCE_CELLS) + { + field.reserve(mesh.cells_size()); + + for (auto &cl : mesh) + { + const auto di = degree_infos.cellDegreeInfo(mesh, cl); + const auto num_cell_dofs = vector_basis_size(di.cell_degree(), Mesh::dimension, Mesh::dimension); + size_t num_faces_dofs = 0; + + if (name == DEPL) + { + const auto fcs = faces(mesh, cl); + const auto fcs_di = di.facesDegreeInfo(); + + for (auto &fc_di : fcs_di) + { + if (fc_di.hasUnknowns()) + { + num_faces_dofs += vector_basis_size(fc_di.degree(), Mesh::dimension - 1, Mesh::dimension); + } + } + } + + field.push_back(vector_type::Zero(num_cell_dofs + num_faces_dofs)); + } + } + else if (name == DEPL_FACES) + { + field.reserve(mesh.faces_size()); + + for (auto itor = mesh.faces_begin(); itor != mesh.faces_end(); itor++) + { + const auto fc = *itor; + const auto di = degree_infos.degreeInfo(mesh, fc); + + size_t num_face_dofs = 0; + if (di.hasUnknowns()) + { + num_face_dofs = vector_basis_size(di.degree(), Mesh::dimension - 1, Mesh::dimension); + } + + field.push_back(vector_type::Zero(num_face_dofs)); + } + } + else + { + throw std::runtime_error("Unknown field"); + } + + this->setField(name, field); + } + + template + void createField(const FieldName &name, const Mesh &mesh, const MeshDegreeInfo °ree_infos, const vector_rhs_function func) + { + std::vector field; + + if (name == DEPL || name == DEPL_CELLS || name == VITE_CELLS || name == ACCE_CELLS) + { + field.reserve(mesh.cells_size()); + + for (auto &cl : mesh) + { + if (name == DEPL) + { + field.push_back(project_function(mesh, cl, degree_infos, func, 2)); + } + else + { + const auto di = degree_infos.cellDegreeInfo(mesh, cl); + + field.push_back(project_function(mesh, cl, di.cell_degree(), func, 2)); + } + } + } + else if (name == DEPL_FACES) + { + field.reserve(mesh.faces_size()); + + for (auto itor = mesh.faces_begin(); itor != mesh.faces_end(); itor++) + { + const auto fc = *itor; + const auto fdi = degree_infos.degreeInfo(mesh, fc); + + field.push_back(project_function(mesh, fc, fdi.degree(), func, 2)); + } + } + else + { + throw std::runtime_error("Unknown field"); + } + + this->setField(name, field); + } + }; + + template + class MultiTimeField + { + private: + typedef dynamic_vector vector_type; + typedef TimeField field_type; + + std::vector m_fields; + + public: + MultiTimeField() {}; + + MultiTimeField(const int n_fields) + { + m_fields.resize(n_fields); + }; + + void + setCurrentTime(const scalar_type& time) + { + m_fields.at(0).setTime(time); + } + + field_type getCurrentTimeField() const + { + return this->getTimeField(0); + } + + field_type getPreviousTimeField() const + { + return this->getTimeField(-1); + } + + field_type getTimeField(const int &relative_index) const + { + return m_fields.at(-relative_index); + } + + void setTimeField(const int &relative_index, const field_type &field) + { + m_fields.at(-relative_index) = field; + } + + void setCurrentTimeField(const field_type &field) + { + return this->setTimeField(0, field); + } + + auto + getField(const int &relative_index, const FieldName &name) const + { + return m_fields.at(-relative_index).getField(name); + } + + auto + getCurrentField(const FieldName &name) const + { + return m_fields.at(0).getField(name); + } + + void + setField(const int &relative_index, const FieldName &name, const std::vector &field) + { + return m_fields.at(-relative_index).setField(name, field); + } + + void setCurrentField(const FieldName &name, const std::vector &field) + { + return this->setField(0, name, field); + } + + template + void createField(const int &relative_index, const FieldName &name, const Mesh &mesh, const MeshDegreeInfo °ree_infos, const vector_rhs_function func) + { + m_fields.at(-relative_index).createField(name, mesh, degree_infos, func); + } + + auto getNumberOfTimeField() const + { + return m_fields.size(); + } + + void update() + { + + int n_fields = this->getNumberOfTimeField(); + + for (int i = n_fields - 1; i > 0; i--) + { + m_fields.at(i) = m_fields.at(i - 1); + } + } + }; + } + +} // end disk diff --git a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonIteration.hpp b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonIteration.hpp index a0d124f8..7cca0e8a 100755 --- a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonIteration.hpp +++ b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonIteration.hpp @@ -38,14 +38,19 @@ #include "diskpp/boundary_conditions/boundary_conditions.hpp" #include "diskpp/common/timecounter.hpp" #include "diskpp/mechanics/NewtonSolver/NewtonSolverComput.hpp" +#include "diskpp/mechanics/NewtonSolver/NewtonSolverDynamic.hpp" #include "diskpp/mechanics/NewtonSolver/NewtonSolverInformations.hpp" #include "diskpp/mechanics/NewtonSolver/NewtonSolverParameters.hpp" #include "diskpp/mechanics/NewtonSolver/StabilizationManager.hpp" #include "diskpp/mechanics/NewtonSolver/TimeManager.hpp" +#include "diskpp/mechanics/NewtonSolver/Fields.hpp" #include "diskpp/mechanics/behaviors/laws/behaviorlaws.hpp" #include "diskpp/methods/hho" #include "diskpp/solvers/solver.hpp" +#include "mumps.hpp" + + namespace disk { @@ -62,21 +67,23 @@ namespace mechanics * * @tparam MeshType type of the mesh */ -template +template class NewtonIteration { - typedef MeshType mesh_type; + typedef MeshType mesh_type; + typedef typename mesh_type::cell cell_type; typedef typename mesh_type::coordinate_type scalar_type; typedef dynamic_matrix matrix_type; typedef dynamic_vector vector_type; - typedef NewtonSolverParameter param_type; + typedef NewtonSolverParameter param_type; typedef vector_boundary_conditions bnd_type; - typedef Behavior behavior_type; + typedef Behavior behavior_type; typedef vector_mechanics_hho_assembler assembler_type; - typedef mechanical_computation elem_type; + typedef mechanical_computation elem_type; + typedef dynamic_computation dyna_type; vector_type m_system_displ; @@ -85,25 +92,101 @@ class NewtonIteration std::vector m_bL; std::vector m_AL; - std::vector m_postprocess_data; - std::vector m_displ, m_displ_faces; - std::vector m_velocity, m_acce; - std::vector m_velocity_p, m_acce_p, m_acce_pred; - TimeStep m_time_step; - scalar_type m_F_int; + dyna_type m_dyna; + + scalar_type m_F_int, m_sol_norm; bool m_verbose; - public: - NewtonIteration(const mesh_type& msh, - const bnd_type& bnd, - const param_type& rp, - const MeshDegreeInfo& degree_infos, - const TimeStep& current_step) : - m_verbose(rp.m_verbose), - m_time_step(current_step) + matrix_type + _gradrec(const mesh_type &msh, + const cell_type &cl, + const param_type &rp, + const MeshDegreeInfo °ree_infos, + const bool &small_def, + const std::vector &gradient_precomputed) const + { + if (rp.m_precomputation) + { + const auto cell_i = msh.lookup(cl); + return gradient_precomputed[cell_i]; + } + + if (small_def) + { + const auto gradrec_sym_full = make_matrix_hho_symmetric_gradrec(msh, cl, degree_infos); + return gradrec_sym_full.first; + } + else + { + const auto gradrec_full = make_matrix_hho_gradrec(msh, cl, degree_infos); + return gradrec_full.first; + } + + return matrix_type(); + } + + matrix_type + _stab(const mesh_type &msh, + const cell_type &cl, + const param_type &rp, + const MeshDegreeInfo °ree_infos, + const std::vector &stab_precomputed) const + { + if (rp.m_precomputation) + { + const auto cell_i = msh.lookup(cl); + return stab_precomputed.at(cell_i); + } + else + { + switch (rp.m_stab_type) + { + case HHO_SYM: + { + const auto recons = make_vector_hho_symmetric_laplacian(msh, cl, degree_infos); + return make_vector_hho_stabilization(msh, cl, recons.first, + degree_infos); + break; + } + case HHO: + { + const auto recons_scalar = make_scalar_hho_laplacian(msh, cl, degree_infos); + return make_vector_hho_stabilization_optim(msh, cl, recons_scalar.first, degree_infos); + break; + } + case HDG: + { + return make_vector_hdg_stabilization(msh, cl, degree_infos); + break; + } + case DG: + { + return make_vector_dg_stabilization(msh, cl, degree_infos); + break; + } + case NO: + { + break; + } + default: + throw std::invalid_argument("Unknown stabilization"); + } + } + + return matrix_type(); + } + +public: + NewtonIteration(const mesh_type &msh, + const bnd_type &bnd, + const param_type &rp, + const MeshDegreeInfo °ree_infos, + const TimeStep ¤t_step) : m_verbose(rp.m_verbose), + m_time_step(current_step), + m_dyna(rp) { m_AL.clear(); m_AL.resize(msh.cells_size()); @@ -127,77 +210,26 @@ class NewtonIteration } void - initialize(const mesh_type& msh, - const param_type& rp, - const std::vector& initial_displ, - const std::vector& initial_displ_faces, - const std::vector& initial_velocity, - const std::vector& initial_acce) + initialize(const mesh_type &msh, + const MultiTimeField &fields) { - m_displ_faces.clear(); - m_displ_faces = initial_displ_faces; - - m_displ.clear(); - m_displ = initial_displ; - - m_velocity.clear(); - m_velocity = initial_velocity; - - m_acce.clear(); - m_acce = initial_acce; - - m_velocity_p.clear(); - m_velocity_p = initial_velocity; - - m_acce_p.clear(); - m_acce_p = initial_acce; - - m_acce_pred.clear(); - m_acce_pred.reserve(msh.cells_size()); - - if (rp.isUnsteady()) - { - auto dyna_rp = rp.getUnsteadyParameters(); - auto beta = dyna_rp["beta"]; - scalar_type dt = m_time_step.increment_time(); - - m_acce_pred.clear(); - m_acce_pred.reserve(msh.cells_size()); - - for (auto& cl : msh) - { - const auto cl_id = msh.lookup(cl); - const auto num_cell_dofs = m_acce[cl_id].size(); - - auto uT = m_displ[cl_id].head(num_cell_dofs); - auto acce_pred = -uT / (beta * dt * dt) - m_velocity[cl_id] / (beta * dt) - - dt * dt / 2.0 * (1.0 - 2.0 * beta) * m_acce[cl_id]; - - m_acce_pred.push_back(acce_pred); - } - } - else - { - for (auto& cl : msh) - { - m_acce_pred.push_back(vector_type::Zero(1)); - } - } + m_dyna.prediction(msh, fields, m_time_step); } - template + template AssemblyInfo - assemble(const mesh_type& msh, - const bnd_type& bnd, - const param_type& rp, - const MeshDegreeInfo& degree_infos, - const LoadFunction& lf, - const std::vector& gradient_precomputed, - const std::vector& stab_precomputed, - behavior_type& behavior, - StabCoeffManager& stab_manager) + assemble(const mesh_type &msh, + const bnd_type &bnd, + const param_type &rp, + const MeshDegreeInfo °ree_infos, + const LoadFunction &lf, + const std::vector &gradient_precomputed, + const std::vector &stab_precomputed, + behavior_type &behavior, + StabCoeffManager &stab_manager, + const MultiTimeField &fields) { - elem_type elem; + elem_type elem; AssemblyInfo ai; // set RHS to zero @@ -206,42 +238,32 @@ class NewtonIteration const bool small_def = (behavior.getDeformation() == SMALL_DEF); + const auto depl = fields.getCurrentField(FieldName::DEPL); + const auto depl_faces = fields.getCurrentField(FieldName::DEPL_FACES); + + m_sol_norm = 0.0; + for (auto &uF : depl_faces) + { + m_sol_norm += uF.squaredNorm(); + } + timecounter tc, ttot; ttot.tic(); - for (auto& cl : msh) + for (auto &cl : msh) { const auto cell_i = msh.lookup(cl); - const auto di = degree_infos.cellDegreeInfo(msh, cl); + + const auto huT = depl.at(cell_i); // Gradient Reconstruction // std::cout << "Grad" << std::endl; - matrix_type GT; tc.tic(); - if (rp.m_precomputation) - { - GT = gradient_precomputed[cell_i]; - } - else - { - if (small_def) - { - const auto gradrec_sym_full = make_matrix_hho_symmetric_gradrec(msh, cl, degree_infos); - GT = gradrec_sym_full.first; - } - else - { - const auto gradrec_full = make_matrix_hho_gradrec(msh, cl, degree_infos); - GT = gradrec_full.first; - } - } + matrix_type GT = _gradrec(msh, cl, rp, degree_infos, small_def, gradient_precomputed); tc.toc(); ai.m_time_gradrec += tc.elapsed(); - // Begin Assembly - // Build rhs and lhs - // Mechanical Computation tc.tic(); @@ -253,8 +275,7 @@ class NewtonIteration degree_infos, lf, GT, - m_displ.at(cell_i), - m_acce_pred.at(cell_i), + huT, m_time_step, behavior, stab_manager, @@ -271,67 +292,29 @@ class NewtonIteration // Stabilisation Contribution // std::cout << "Stab" << std::endl; tc.tic(); - if (rp.m_stab) { - matrix_type stab; - if (rp.m_precomputation) - { - stab = stab_precomputed.at(cell_i); - } - else - { - switch (rp.m_stab_type) - { - case HHO: - { - // we do not make any difference for the displacement reconstruction - // if (small_def) - // { - // const auto recons = make_vector_hho_symmetric_laplacian(msh, cl, degree_infos); - // stab_HHO = make_vector_hho_stabilization(msh, cl, recons.first, - // degree_infos); - // } - // else - // { - const auto recons_scalar = make_scalar_hho_laplacian(msh, cl, degree_infos); - stab = make_vector_hho_stabilization_optim(msh, cl, recons_scalar.first, degree_infos); - // } - - break; - } - case HDG: - { - stab = make_vector_hdg_stabilization(msh, cl, degree_infos); - break; - } - case DG: - { - stab = make_vector_dg_stabilization(msh, cl, degree_infos); - break; - } - case NO: - { - break; - } - default: throw std::invalid_argument("Unknown stabilization"); - } - } - - assert(elem.K_int.rows() == stab.rows()); - assert(elem.K_int.cols() == stab.cols()); - assert(elem.RTF.rows() == stab.rows()); - assert(elem.RTF.cols() == m_displ.at(cell_i).cols()); - const auto beta_s = stab_manager.getValue(msh, cl); + + matrix_type stab = beta_s * _stab(msh, cl, rp, degree_infos, stab_precomputed); + // std::cout << beta_s << std::endl; - lhs += beta_s * stab; - rhs -= beta_s * stab * m_displ.at(cell_i); + lhs += stab; + rhs -= stab * huT; } tc.toc(); ai.m_time_stab += tc.elapsed(); + // Dynamic contribution + if (m_dyna.enable()) + { + m_dyna.compute(msh, cl, degree_infos, huT, m_time_step); + lhs += m_dyna.K_iner; + rhs += m_dyna.RTF; + ai.m_time_dyna += m_dyna.time_dyna; + } + // std::cout << "R: " << rhs.norm() << std::endl; // std::cout << rhs.transpose() << std::endl; @@ -346,14 +329,14 @@ class NewtonIteration tc.toc(); ai.m_time_statcond += tc.elapsed(); - const auto& lc = std::get<0>(scnp); + const auto &lc = std::get<0>(scnp); // std::cout << "rhs: " << lc.second.norm() << std::endl; // std::cout << lc.second.transpose() << std::endl; // std::cout << "lhs: " << lc.first.norm() << std::endl; // std::cout << "Assemb" << std::endl; - m_assembler.assemble_nonlinear(msh, cl, bnd, lc.first, lc.second, m_displ_faces); + m_assembler.assemble_nonlinear(msh, cl, bnd, lc.first, lc.second, depl_faces); } m_F_int = sqrt(m_F_int); @@ -363,7 +346,7 @@ class NewtonIteration m_assembler.finalize(); ttot.toc(); - ai.m_time_assembly = ttot.elapsed(); + ai.m_time_assembly = ttot.elapsed(); ai.m_linear_system_size = m_assembler.LHS.rows(); return ai; } @@ -371,7 +354,6 @@ class NewtonIteration SolveInfo solve(void) { - // std::cout << "begin solve" << std::endl; timecounter tc; tc.tic(); @@ -380,60 +362,48 @@ class NewtonIteration #ifdef HAVE_INTEL_MKL solvers::pardiso_params pparams; mkl_pardiso(pparams, m_assembler.LHS, m_assembler.RHS, m_system_displ); +#elif HAVE_MUMPS + m_system_displ = mumps_lu(m_assembler.LHS, m_assembler.RHS); #else - throw std::runtime_error("Pardiso is not installed"); + throw std::runtime_error("No linear solver avalaible."); #endif tc.toc(); // std::cout << "LHS" << m_assembler.LHS << std::endl; // std::cout << "RHS" << m_assembler.RHS << std::endl; - // std::cout << "end solve" << std::endl; - return SolveInfo(m_assembler.LHS.rows(), m_assembler.LHS.nonZeros(), tc.elapsed()); } scalar_type - postprocess(const mesh_type& msh, - const bnd_type& bnd, - const param_type& rp, - const MeshDegreeInfo& degree_infos) + postprocess(const mesh_type &msh, + const bnd_type &bnd, + const param_type &rp, + const MeshDegreeInfo °ree_infos, + MultiTimeField &fields) { - // std::cout << "begin post_process" << std::endl; timecounter tc; tc.tic(); - // std::cout << m_system_displ << std::endl; + auto depl_faces = fields.getCurrentField(FieldName::DEPL_FACES); + auto depl_cells = fields.getCurrentField(FieldName::DEPL_CELLS); + auto depl = fields.getCurrentField(FieldName::DEPL); // Update cell - for (auto& cl : msh) + for (auto &cl : msh) { const auto cell_i = msh.lookup(cl); const vector_type xdT = - m_assembler.take_local_solution_nonlinear(msh, cl, bnd, m_system_displ, m_displ_faces); - - vector_type x_cond = xdT; + m_assembler.take_local_solution_nonlinear(msh, cl, bnd, m_system_displ, depl_faces); // static decondensation - const vector_type xT = m_bL[cell_i] - m_AL[cell_i] * x_cond; - - assert(m_displ.at(cell_i).size() == xT.size() + xdT.size()); - // Update element U^{i+1} = U^i + delta U^i /// - m_displ.at(cell_i).head(xT.size()) += xT; - m_displ.at(cell_i).segment(xT.size(), xdT.size()) += xdT; + const vector_type xT = m_bL[cell_i] - m_AL[cell_i] * xdT; - if (rp.isUnsteady()) - { - auto dyna_rp = rp.getUnsteadyParameters(); - auto beta = dyna_rp["beta"]; - auto gamma = dyna_rp["gamma"]; - scalar_type dt = m_time_step.increment_time(); - - m_acce.at(cell_i) = xT / (beta * dt * dt) + m_acce_pred.at(cell_i); - m_velocity.at(cell_i) = - m_velocity_p.at(cell_i) + dt * ((1.0 - gamma) * m_acce.at(cell_i) + gamma * m_acce.at(cell_i)); - } + // Update element U^{i+1} = U^i + delta U^i + depl.at(cell_i).head(xT.size()) += xT; + depl.at(cell_i).tail(xdT.size()) += xdT; + depl_cells.at(cell_i) += xT; // std::cout << "KT_F " << m_AL[cell_i].norm() << std::endl; // std::cout << "sol_F" << std::endl; @@ -442,36 +412,35 @@ class NewtonIteration // std::cout << m_bL[cell_i].transpose() << std::endl; // std::cout << "sol_T" << std::endl; // std::cout << xT.transpose() << std::endl; - // std::cout << (m_displ.at(cell_i)).segment(0, xT.size()).transpose() << std::endl; + // std::cout << depl.at(cell_i).transpose() << std::endl; } + fields.setCurrentField(FieldName::DEPL, depl); + fields.setCurrentField(FieldName::DEPL_CELLS, depl_cells); // Update unknowns // Update face Uf^{i+1} = Uf^i + delta Uf^i - size_t face_i = 0; + + auto depl_faces_new = fields.getCurrentField(FieldName::DEPL_FACES); + int face_i = 0; for (auto itor = msh.faces_begin(); itor != msh.faces_end(); itor++) { const auto fc = *itor; - m_displ_faces.at(face_i++) += - m_assembler.take_local_solution_nonlinear(msh, fc, bnd, m_system_displ, m_displ_faces); + depl_faces_new.at(face_i++) += + m_assembler.take_local_solution_nonlinear(msh, fc, bnd, m_system_displ, depl_faces); } + fields.setCurrentField(FieldName::DEPL_FACES, depl_faces_new); + + m_dyna.postprocess(msh, fields, m_time_step); - // std::cout << "end post_process" << std::endl; tc.toc(); return tc.elapsed(); } bool - convergence(const param_type& rp, const size_t iter) + convergence(const param_type &rp, const size_t iter) { // norm of the solution - scalar_type error_un = 0; - for (size_t i = 0; i < m_displ_faces.size(); i++) - { - scalar_type norm = m_displ_faces[i].norm(); - error_un += norm * norm; - } - - error_un = std::sqrt(error_un); + auto error_un = std::sqrt(m_sol_norm); if (error_un <= scalar_type(10E-15)) { @@ -479,15 +448,15 @@ class NewtonIteration } // norm of the rhs - const scalar_type residual = m_assembler.RHS.norm(); - scalar_type max_error = 0.0; + const scalar_type residual = m_assembler.RHS.norm(); + scalar_type max_error = 0.0; for (size_t i = 0; i < m_assembler.RHS.size(); i++) max_error = std::max(max_error, std::abs(m_assembler.RHS(i))); // norm of the increment - const scalar_type error_incr = m_system_displ.norm(); - scalar_type relative_displ = error_incr / error_un; - scalar_type relative_error = residual / m_F_int; + const scalar_type error_incr = m_system_displ.norm(); + scalar_type relative_displ = error_incr / error_un; + scalar_type relative_error = residual / m_F_int; if (iter == 0) { @@ -537,29 +506,6 @@ class NewtonIteration return false; } } - - void - save_solutions(std::vector& displ, - std::vector& displ_faces, - std::vector& velocity, - std::vector& acce) - { - displ_faces.clear(); - displ_faces = m_displ_faces; - assert(m_displ_faces.size() == displ_faces.size()); - - displ.clear(); - displ = m_displ; - assert(m_displ.size() == displ.size()); - - velocity.clear(); - velocity = m_velocity; - assert(m_velocity.size() == velocity.size()); - - acce.clear(); - acce = m_acce; - assert(m_acce.size() == acce.size()); - } }; } } // end diskpp \ No newline at end of file diff --git a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolver.hpp b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolver.hpp index d95e4dd9..74c30996 100644 --- a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolver.hpp +++ b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolver.hpp @@ -37,6 +37,7 @@ #include "diskpp/mechanics/NewtonSolver/NewtonStep.hpp" #include "diskpp/mechanics/NewtonSolver/StabilizationManager.hpp" #include "diskpp/mechanics/NewtonSolver/TimeManager.hpp" +#include "diskpp/mechanics/NewtonSolver/Fields.hpp" #include "diskpp/mechanics/behaviors/laws/behaviorlaws.hpp" #include "diskpp/mechanics/behaviors/tensor_conversion.hpp" #include "diskpp/mechanics/stress_tensors.hpp" @@ -91,8 +92,7 @@ class NewtonSolver MeshDegreeInfo m_degree_infos; StabCoeffManager m_stab_manager; - std::vector m_displ, m_displ_faces; - std::vector m_velocity, m_acce; + MultiTimeField m_fields; std::vector m_gradient_precomputed, m_stab_precomputed; PostMesh m_post_mesh; @@ -134,68 +134,17 @@ class NewtonSolver void init(void) { - const auto dimension = mesh_type::dimension; - - size_t total_dof = 0; - - m_displ.clear(); - m_displ_faces.clear(); - - m_displ.reserve(m_msh.cells_size()); - m_displ_faces.reserve(m_msh.faces_size()); - + TimeField tf; + tf.createZeroField(FieldName::DEPL, m_msh, m_degree_infos); + tf.createZeroField(FieldName::DEPL_CELLS, m_msh, m_degree_infos); + tf.createZeroField(FieldName::DEPL_FACES, m_msh, m_degree_infos); if (m_rp.isUnsteady()) { - m_velocity.clear(); - m_acce.clear(); - - m_velocity.reserve(m_msh.cells_size()); - m_acce.reserve(m_msh.cells_size()); - } - - for (auto& cl : m_msh) - { - const auto di = m_degree_infos.cellDegreeInfo(m_msh, cl); - const auto num_cell_dofs = vector_basis_size(di.cell_degree(), dimension, dimension); - total_dof += num_cell_dofs; - - const auto fcs = faces(m_msh, cl); - const auto fcs_di = di.facesDegreeInfo(); - - const auto grad_degree = di.grad_degree(); - - size_t num_faces_dofs = 0; - for (auto& fc_di : fcs_di) - { - if (fc_di.hasUnknowns()) - { - num_faces_dofs += vector_basis_size(fc_di.degree(), dimension - 1, dimension); - } - } - - m_displ.push_back(vector_type::Zero(num_cell_dofs + num_faces_dofs)); - - if (m_rp.isUnsteady()) - { - m_velocity.push_back(vector_type::Zero(num_cell_dofs)); - m_acce.push_back(vector_type::Zero(num_cell_dofs)); - } - } - - for (auto itor = m_msh.faces_begin(); itor != m_msh.faces_end(); itor++) - { - const auto fc = *itor; - const auto di = m_degree_infos.degreeInfo(m_msh, fc); - - size_t num_face_dofs = 0; - if (di.hasUnknowns()) - { - num_face_dofs = vector_basis_size(di.degree(), dimension - 1, dimension); - } - total_dof += num_face_dofs; - - m_displ_faces.push_back(vector_type::Zero(num_face_dofs)); + tf.createZeroField(FieldName::VITE_CELLS, m_msh, m_degree_infos); + tf.createZeroField(FieldName::ACCE_CELLS, m_msh, m_degree_infos); } + m_fields.setCurrentTimeField(tf); + m_fields.setTimeField(-1, tf); // compute mesh for post-processing m_post_mesh = PostMesh(m_msh); @@ -205,9 +154,7 @@ class NewtonSolver std::cout << "** Numbers of cells: " << m_msh.cells_size() << std::endl; std::cout << "** Numbers of faces: " << m_msh.faces_size() << " ( boundary faces: " << m_msh.boundary_faces_size() << " )" << std::endl; - std::cout << "** Numbers of dofs: " << std::endl; - std::cout << " ** Before static condensation: " << total_dof << std::endl; - std::cout << " ** After static condensation: " << this->numberOfDofs() << std::endl; + std::cout << "** Numbers of dofs after static condensation: " << this->numberOfDofs() << std::endl; std::cout << " " << std::endl; } } @@ -246,19 +193,16 @@ class NewtonSolver { case HHO: { - // we do not make any difference for the displacement reconstruction - // if (m_behavior.getDeformation() == SMALL_DEF) - // { - // const auto recons = make_vector_hho_symmetric_laplacian(m_msh, cl, m_degree_infos); - // m_stab_precomputed.push_back( - // make_vector_hho_stabilization(m_msh, cl, recons.first, m_degree_infos)); - // } - // else - // { const auto recons_scalar = make_scalar_hho_laplacian(m_msh, cl, m_degree_infos); m_stab_precomputed.push_back( - make_vector_hho_stabilization_optim(m_msh, cl, recons_scalar.first, m_degree_infos)); - // } + make_vector_hho_stabilization_optim(m_msh, cl, recons_scalar.first, m_degree_infos)); + break; + } + case HHO_SYM: + { + const auto recons = make_vector_hho_symmetric_laplacian(m_msh, cl, m_degree_infos); + m_stab_precomputed.push_back( + make_vector_hho_stabilization(m_msh, cl, recons.first, m_degree_infos)); break; } case HDG: @@ -282,58 +226,60 @@ class NewtonSolver } public: - NewtonSolver(const mesh_type& msh, const bnd_type& bnd, const param_type& rp) : - m_msh(msh), m_verbose(rp.m_verbose), m_convergence(false), m_rp(rp), m_bnd(bnd), m_stab_manager(msh, rp.m_beta) - { - if (m_verbose) - { - std::cout << "------------------------------------------------------------------------" - "----------------------" - << std::endl; - std::cout - << "|********************** Nonlinear Newton's Solver for solid mechanics ***********************|" - << std::endl; - std::cout << "------------------------------------------------------------------------" - "----------------------" - << std::endl; - } - int face_degree = rp.m_face_degree; - if (rp.m_face_degree < 0) - { - std::cout << "'face_degree' should be > 0. Reverting to 1." << std::endl; - face_degree = 1; - } - - m_rp.m_face_degree = face_degree; - - int cell_degree = rp.m_cell_degree; - if ((face_degree - 1 > cell_degree) or (cell_degree > face_degree + 1)) - { - std::cout << "'cell_degree' should be 'face_degree + 1' =>" - << "'cell_degree' => 'face_degree -1'. Reverting to 'face_degree'." << std::endl; - cell_degree = face_degree; - } - - m_rp.m_cell_degree = cell_degree; - - int grad_degree = rp.m_grad_degree; - if (grad_degree < face_degree) - { - std::cout << "'grad_degree' should be >= 'face_degree'. Reverting to 'face_degree'." << std::endl; - grad_degree = face_degree; - } - - if (m_verbose) - { - m_rp.infos(); - } - - // Initialization - if (m_verbose) - std::cout << "Initialization ..." << std::endl; - this->init_degree(cell_degree, face_degree, grad_degree); - this->init(); - } + NewtonSolver(const mesh_type &msh, const bnd_type &bnd, const param_type &rp) : m_msh(msh), m_verbose(rp.m_verbose), + m_convergence(false), m_rp(rp), m_bnd(bnd), + m_stab_manager(msh, rp.m_beta), + m_fields(2) + { + if (m_verbose) + { + std::cout << "------------------------------------------------------------------------" + "----------------------" + << std::endl; + std::cout + << "|********************** Nonlinear Newton's Solver for solid mechanics ***********************|" + << std::endl; + std::cout << "------------------------------------------------------------------------" + "----------------------" + << std::endl; + } + int face_degree = rp.m_face_degree; + if (rp.m_face_degree < 0) + { + std::cout << "'face_degree' should be > 0. Reverting to 1." << std::endl; + face_degree = 1; + } + + m_rp.m_face_degree = face_degree; + + int cell_degree = rp.m_cell_degree; + if ((face_degree - 1 > cell_degree) or (cell_degree > face_degree + 1)) + { + std::cout << "'cell_degree' should be 'face_degree + 1' =>" + << "'cell_degree' => 'face_degree -1'. Reverting to 'face_degree'." << std::endl; + cell_degree = face_degree; + } + + m_rp.m_cell_degree = cell_degree; + + int grad_degree = rp.m_grad_degree; + if (grad_degree < face_degree) + { + std::cout << "'grad_degree' should be >= 'face_degree'. Reverting to 'face_degree'." << std::endl; + grad_degree = face_degree; + } + + if (m_verbose) + { + m_rp.infos(); + } + + // Initialization + if (m_verbose) + std::cout << "Initialization ..." << std::endl; + this->init_degree(cell_degree, face_degree, grad_degree); + this->init(); + } /** * @brief return a boolean to know if the verbosity mode is activated @@ -366,36 +312,11 @@ class NewtonSolver { size_t cell_i = 0; - for (auto& cl : m_msh) - { - m_displ.at(cell_i++) = project_function(m_msh, cl, m_degree_infos, func, 2); - } - - for (auto itor = m_msh.faces_begin(); itor != m_msh.faces_end(); itor++) - { - const auto bfc = *itor; - const auto face_id = m_msh.lookup(bfc); - const auto fdi = m_degree_infos.degreeInfo(m_msh, bfc); - - if (m_bnd.contact_boundary_type(face_id) == SIGNORINI_FACE) - { - const auto proj_bcf = project_function(m_msh, bfc, fdi.degree(), func, 2); - assert(m_displ_faces[face_id].size() == proj_bcf.size()); + m_fields.createField(0, FieldName::DEPL, m_msh, m_degree_infos, func); + m_fields.createField(0, FieldName::DEPL_CELLS, m_msh, m_degree_infos, func); + m_fields.createField(0, FieldName::DEPL_FACES, m_msh, m_degree_infos, func); - m_displ_faces[face_id] = proj_bcf; - } - else if (m_bnd.contact_boundary_type(face_id) == SIGNORINI_CELL) - { - assert(m_displ_faces[face_id].size() == 0); - } - else - { - const auto proj_bcf = project_function(m_msh, bfc, fdi.degree(), func, 2); - assert(m_displ_faces[face_id].size() == proj_bcf.size()); - - m_displ_faces[face_id] = proj_bcf; - } - } + /*TODO: look for contact.*/ } /** @@ -404,22 +325,16 @@ class NewtonSolver * @param func given function */ void - initial_fields(const vector_rhs_function depl, - const vector_rhs_function velo, - const vector_rhs_function acce) + initial_field(const FieldName &name, const vector_rhs_function func) { - if (m_rp.isUnsteady()) - { - this->initial_guess(depl); - for (auto& cl : m_msh) - { - const auto cell_id = m_msh.lookup(cl); - const auto di = m_degree_infos.cellDegreeInfo(m_msh, cl); - - m_velocity.at(cell_id) = project_function(m_msh, cl, di.cell_degree(), velo, 2); - m_acce.at(cell_id) = project_function(m_msh, cl, di.cell_degree(), acce, 2); - } + if(name == FieldName::DEPL){ + m_fields.createField(-1, FieldName::DEPL_CELLS, m_msh, m_degree_infos, func); + m_fields.createField(-1, FieldName::DEPL_FACES, m_msh, m_degree_infos, func); + m_fields.createField(0, FieldName::DEPL_CELLS, m_msh, m_degree_infos, func); + m_fields.createField(0, FieldName::DEPL_FACES, m_msh, m_degree_infos, func); } + m_fields.createField(-1, name, m_msh, m_degree_infos, func); + m_fields.createField(0, name, m_msh, m_degree_infos, func); } /** @@ -543,23 +458,22 @@ class NewtonSolver // Newton step NewtonStep newton_step(m_rp); - newton_step.initialize(m_displ, m_displ_faces, m_velocity, m_acce); // Loop on time step while (!list_time_step.empty()) { const auto current_step = list_time_step.getCurrentTimeStep(); const auto current_time = current_step.end_time(); + m_fields.setCurrentTime(current_time); if (m_verbose) { list_time_step.printCurrentTimeStep(); } - auto rlf = [&lf, ¤t_time](const point& p) -> auto + auto rlf = [&lf, ¤t_time](const point &p) -> auto { return lf(p, current_time); }; - - m_bnd.multiplyAllFunctionsByAFactor(current_time); + m_bnd.setTime(current_time); // Newton correction NewtonSolverInfo newton_info = newton_step.compute(m_msh, @@ -571,7 +485,8 @@ class NewtonSolver m_gradient_precomputed, m_stab_precomputed, m_behavior, - m_stab_manager); + m_stab_manager, + m_fields); si.updateInfo(newton_info); if (m_verbose) @@ -608,36 +523,30 @@ class NewtonSolver list_time_step.removeCurrentTimeStep(); m_behavior.update(); m_stab_manager.update(); + m_fields.update(); - if (time_saving) + if (time_saving && (m_rp.m_time_save.front() < current_time + 1E-5)) { - if (m_rp.m_time_save.front() < current_time + 1E-5) - { - newton_step.save_solutions(m_displ, m_displ_faces, m_velocity, m_acce); - - std::cout << "** Save results" << std::endl; - std::string name = - "result" + std::to_string(mesh_type::dimension) + "D_t" + std::to_string(current_time) + "_"; - - this->output_discontinuous_displacement(name + "depl_disc.msh"); - this->output_continuous_displacement(name + "depl_cont.msh"); - this->output_CauchyStress_GP(name + "CauchyStress_GP.msh"); - this->output_CauchyStress_GP(name + "CauchyStress_GP_def.msh", true); - this->output_discontinuous_deformed(name + "deformed_disc.msh"); - this->output_is_plastic_GP(name + "plastic_GP.msh"); - this->output_stabCoeff(name + "stabCoeff.msh"); - this->output_equivalentPlasticStrain_GP(name + "equivalentPlasticStrain_GP.msh"); - - m_rp.m_time_save.pop_front(); - if (m_rp.m_time_save.empty()) - time_saving = false; - } + std::cout << "** Save results" << std::endl; + std::string name = + "result" + std::to_string(mesh_type::dimension) + "D_t" + std::to_string(current_time) + "_"; + + this->output_discontinuous_displacement(name + "depl_disc.msh"); + this->output_continuous_displacement(name + "depl_cont.msh"); + this->output_CauchyStress_GP(name + "CauchyStress_GP.msh"); + this->output_CauchyStress_GP(name + "CauchyStress_GP_def.msh", true); + this->output_discontinuous_deformed(name + "deformed_disc.msh"); + this->output_is_plastic_GP(name + "plastic_GP.msh"); + this->output_stabCoeff(name + "stabCoeff.msh"); + this->output_equivalentPlasticStrain_GP(name + "equivalentPlasticStrain_GP.msh"); + + m_rp.m_time_save.pop_front(); + if (m_rp.m_time_save.empty()) + time_saving = false; } } } - // save solutions - newton_step.save_solutions(m_displ, m_displ_faces, m_velocity, m_acce); si.m_time_step = list_time_step.numberOfTimeStep(); ttot.toc(); @@ -674,13 +583,13 @@ class NewtonSolver printSolutionCell() const { size_t cell_i = 0; + const auto depl_cells = m_fields.getCurrentField(FieldName::DEPL_CELLS); + std::cout << "Solution at the cells:" << std::endl; for (auto& cl : m_msh) { - const auto di = m_degree_infos.degreeInfo(m_msh, cl); - const auto num_cell_dofs = vector_basis_size(di.cell_degree(), mesh_type::dimension, mesh_type::dimension); std::cout << "cell " << cell_i << ": " << std::endl; - std::cout << m_displ.at(cell_i++).head(num_cell_dofs).transpose() << std::endl; + std::cout << depl_cells.at(cell_i++).transpose() << std::endl; } } @@ -692,12 +601,12 @@ class NewtonSolver scalar_type err_dof = 0; size_t cell_i = 0; + const auto depl_cells = m_fields.getCurrentField(FieldName::DEPL_CELLS); for (auto& cl : m_msh) { const auto cdi = m_degree_infos.degreeInfo(m_msh, cl); - const auto num_cell_dofs = vector_basis_size(cdi.degree(), mesh_type::dimension, mesh_type::dimension); - const vector_type comp_dof = m_displ.at(cell_i++).head(num_cell_dofs); + const vector_type comp_dof = depl_cells.at(cell_i++); const vector_type true_dof = project_function(m_msh, cl, cdi.degree(), as, 2); const auto cb = make_vector_monomial_basis(m_msh, cl, cdi.degree()); @@ -719,11 +628,12 @@ class NewtonSolver scalar_type err_dof = 0; size_t cell_i = 0; + const auto depl = m_fields.getCurrentField(FieldName::DEPL); matrix_type grad; matrix_type stab; - for (auto& cl : m_msh) + for (auto &cl : m_msh) { // std::cout << m_degree_infos.cellDegreeInfo(m_msh, cl) << std::endl; /////// Gradient Reconstruction ///////// @@ -740,44 +650,42 @@ class NewtonSolver { switch (m_rp.m_stab_type) { - case HHO: - { - // we do not make any difference for thre displacement reconstruction - // if (m_behavior.getDeformation() == SMALL_DEF) - // { - // const auto recons = make_vector_hho_symmetric_laplacian(m_msh, cl, m_degree_infos); - // stab = make_vector_hho_stabilization(m_msh, cl, recons.first, - // m_degree_infos); - // } - // else - // { - const auto recons_scalar = make_scalar_hho_laplacian(m_msh, cl, m_degree_infos); - stab = make_vector_hho_stabilization_optim(m_msh, cl, recons_scalar.first, m_degree_infos); - // } - break; - } - case HDG: - { - stab = make_vector_hdg_stabilization(m_msh, cl, m_degree_infos); - break; - } - case DG: - { - stab = make_vector_dg_stabilization(m_msh, cl, m_degree_infos); - break; - } - case NO: - { - break; - stab.setZero(); - } - default: throw std::invalid_argument("Unknown stabilization"); + case HHO_SYM: + { + const auto recons = make_vector_hho_symmetric_laplacian(m_msh, cl, m_degree_infos); + stab = make_vector_hho_stabilization(m_msh, cl, recons.first, + m_degree_infos); + break; + } + case HHO: + { + const auto recons_scalar = make_scalar_hho_laplacian(m_msh, cl, m_degree_infos); + stab = make_vector_hho_stabilization_optim(m_msh, cl, recons_scalar.first, m_degree_infos); + break; + } + case HDG: + { + stab = make_vector_hdg_stabilization(m_msh, cl, m_degree_infos); + break; + } + case DG: + { + stab = make_vector_dg_stabilization(m_msh, cl, m_degree_infos); + break; + } + case NO: + { + break; + stab.setZero(); + } + default: + throw std::invalid_argument("Unknown stabilization"); } } const auto Ah = grad + stab; - const vector_type comp_dof = m_displ.at(cell_i); + const vector_type comp_dof = depl.at(cell_i); const vector_type true_dof = project_function(m_msh, cl, m_degree_infos, as, 2); const vector_type diff_dof = (true_dof - comp_dof); @@ -798,13 +706,15 @@ class NewtonSolver std::vector data; // create data (not used) const std::vector subdata; // create subdata to save soution at gauss point + const auto depl_cells = m_fields.getCurrentField(FieldName::DEPL_CELLS); + int cell_i = 0; int nb_nodes = 0; for (auto& cl : m_msh) { const auto di = m_degree_infos.cellDegreeInfo(m_msh, cl); const auto cb = make_vector_monomial_basis(m_msh, cl, di.cell_degree()); - const vector_type x = m_displ.at(cell_i++).head(cb.size()); + const vector_type x = depl_cells.at(cell_i++); auto cell_nodes = points(m_msh, cl); std::vector new_nodes; @@ -848,6 +758,7 @@ class NewtonSolver const static_vector vzero = static_vector::Zero(); + const auto depl_cells = m_fields.getCurrentField(FieldName::DEPL_CELLS); const size_t nb_nodes(gmsh.getNumberofNodes()); // first(number of data at this node), second(cumulated value) @@ -858,7 +769,7 @@ class NewtonSolver { const auto di = m_degree_infos.cellDegreeInfo(m_msh, cl); const auto cb = make_vector_monomial_basis(m_msh, cl, di.cell_degree()); - const vector_type x = m_displ.at(cell_i).head(cb.size()); + const vector_type x = depl_cells.at(cell_i); auto cell_nodes = m_post_mesh.nodes_cell(cell_i); // Loop on the nodes of the cell @@ -904,12 +815,14 @@ class NewtonSolver std::vector subdata; // create subdata to save soution at gauss point size_t nb_nodes(gmsh.getNumberofNodes()); + const auto depl = m_fields.getCurrentField(FieldName::DEPL); + int cell_i = 0; for (auto& cl : m_msh) { const auto di = m_degree_infos.cellDegreeInfo(m_msh, cl); - const auto uTF = m_displ.at(cell_i); + const auto uTF = depl.at(cell_i); matrix_type gr; if (m_rp.m_precomputation) { @@ -1074,6 +987,7 @@ class NewtonSolver gmsh::Gmesh gmsh(mesh_type::dimension); auto storage = m_msh.backend_storage(); + const auto depl_cells = m_fields.getCurrentField(FieldName::DEPL_CELLS); int cell_i = 0; size_t nb_nodes = 0; for (auto& cl : m_msh) @@ -1081,7 +995,7 @@ class NewtonSolver const auto di = m_degree_infos.cellDegreeInfo(m_msh, cl); auto cb = make_vector_monomial_basis(m_msh, cl, di.cell_degree()); - const vector_type x = m_displ.at(cell_i++).head(cb.size()); + const vector_type x = depl_cells.at(cell_i++); const auto cell_nodes = points(m_msh, cl); std::vector new_nodes; diff --git a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverComput.hpp b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverComput.hpp index b1df5a6e..f5166c6a 100644 --- a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverComput.hpp +++ b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverComput.hpp @@ -326,8 +326,7 @@ class mechanical_computation vector_type RTF; vector_type F_int; double time_law; - double time_contact; - double time_dyna; + double time_contact; mechanical_computation(void) { @@ -339,25 +338,23 @@ class mechanical_computation assert(false); } - template + template void - compute(const mesh_type& msh, - const cell_type& cl, - const bnd_type& bnd, - const param_type& rp, - const MeshDegreeInfo& degree_infos, - const Function& load, - const matrix_type& RkT, - const vector_type& uTF, - const vector_type& aT_pred, - const TimeStep& time_step, - behavior_type& behavior, - StabCoeffManager& stab_manager, - const bool small_def) + compute(const mesh_type &msh, + const cell_type &cl, + const bnd_type &bnd, + const param_type &rp, + const MeshDegreeInfo °ree_infos, + const Function &load, + const matrix_type &RkT, + const vector_type &uTF, + const TimeStep &time_step, + behavior_type &behavior, + StabCoeffManager &stab_manager, + const bool small_def) { time_law = 0.0; time_contact = 0.0; - time_dyna = 0.0; timecounter tc; const auto cell_infos = degree_infos.cellDegreeInfo(msh, cl); @@ -366,7 +363,6 @@ class mechanical_computation const auto cell_degree = cell_infos.cell_degree(); const auto grad_degree = cell_infos.grad_degree(); - const auto cb = make_vector_monomial_basis(msh, cl, cell_degree); const auto cell_basis_size = vector_basis_size(cell_degree, dimension, dimension); size_t gb_size = 0; @@ -527,29 +523,6 @@ class mechanical_computation assert(K_int.rows() == num_total_dofs); assert(K_int.cols() == num_total_dofs); assert(RTF.rows() == num_total_dofs); - - // Unsteady Computation - tc.tic(); - if (rp.isUnsteady()) - { - auto dt = time_step.increment_time(); - auto dyna_rp = rp.getUnsteadyParameters(); - auto beta = dyna_rp["beta"]; - auto rho = dyna_rp["rho"]; - - const matrix_type mass_mat = rho * make_mass_matrix(msh, cl, cb); - - K_int.topLeftCorner(cell_basis_size, cell_basis_size) += mass_mat / (beta * dt * dt); - - auto F_iner = (mass_mat / (beta * dt * dt)) * uTF.head(cell_basis_size); - F_int.head(cell_basis_size) += F_iner; - RTF.head(cell_basis_size) -= F_iner; - - // std::cout << mass_mat.cols() << ", " << aT_pred.rows() << std::endl; - RTF.head(cell_basis_size) += mass_mat * aT_pred; - } - tc.toc(); - time_dyna += tc.elapsed(); } }; } diff --git a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverDynamic.hpp b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverDynamic.hpp new file mode 100644 index 00000000..8dfb1dda --- /dev/null +++ b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverDynamic.hpp @@ -0,0 +1,234 @@ +/* + * /\ Matteo Cicuttin (C) 2016, 2017 + * /__\ matteo.cicuttin@enpc.fr + * /_\/_\ École Nationale des Ponts et Chaussées - CERMICS + * /\ /\ + * /__\ /__\ DISK++, a template library for DIscontinuous SKeletal + * /_\/_\/_\/_\ methods. + * + * This file is copyright of the following authors: + * Nicolas Pignet (C) 2019, 2025 nicolas.pignet@enpc.fr + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * If you use this code or parts of it for scientific publications, you + * are required to cite it as following: + * + * Hybrid High-Order methods for finite elastoplastic deformations + * within a logarithmic strain framework. + * M. Abbas, A. Ern, N. Pignet. + * International Journal of Numerical Methods in Engineering (2019) + * 120(3), 303-327 + * DOI: 10.1002/nme.6137 + */ + +#pragma once + +#include + +#include "diskpp/bases/bases.hpp" +#include "diskpp/common/eigen.hpp" +#include "diskpp/common/timecounter.hpp" +#include "diskpp/mechanics/NewtonSolver/Fields.hpp" +#include "diskpp/mechanics/NewtonSolver/NewtonSolverParameters.hpp" +#include "diskpp/mechanics/NewtonSolver/TimeManager.hpp" +#include "diskpp/methods/hho" +#include "diskpp/quadratures/quadratures.hpp" + +namespace disk +{ + + namespace mechanics + { + + template + class dynamic_computation + { + typedef MeshType mesh_type; + typedef typename mesh_type::coordinate_type scalar_type; + typedef typename mesh_type::cell cell_type; + + typedef dynamic_matrix matrix_type; + typedef dynamic_vector vector_type; + + bool m_enable; + + std::map m_param; + DynamicType m_scheme; + + std::vector m_acce_pred; + + public: + matrix_type K_iner; + vector_type RTF; + + double time_dyna; + + dynamic_computation(const NewtonSolverParameter &rp) + { + m_enable = rp.isUnsteady(); + m_param = rp.getUnsteadyParameters(); + m_scheme = rp.getUnsteadyScheme(); + } + + bool + enable(void) const + { + return m_enable; + } + + void + prediction(const mesh_type &mesh, + const MultiTimeField &fields, + const TimeStep &time_step) + { + if (this->enable()) + { + switch (m_scheme) + { + case DynamicType::NEWMARK: + { + m_acce_pred.clear(); + m_acce_pred.reserve(mesh.cells_size()); + + const auto depl_prev = fields.getField(-1, FieldName::DEPL_CELLS); + const auto velo_prev = fields.getField(-1, FieldName::VITE_CELLS); + const auto acce_prev = fields.getField(-1, FieldName::ACCE_CELLS); + + auto beta = m_param.at("beta"); + scalar_type dt = time_step.increment_time(); + scalar_type denom = 1.0 / (beta * dt * dt); + scalar_type nume = dt * dt / 2.0 * (1.0 - 2.0 * beta); + + for (auto &cl : mesh) + { + const auto cl_id = mesh.lookup(cl); + + auto aT_pred = denom * (depl_prev[cl_id] + dt * velo_prev[cl_id] + nume * acce_prev[cl_id]); + m_acce_pred.push_back(aT_pred); + } + break; + } + default: + break; + } + } + } + + scalar_type + postprocess(const mesh_type &msh, + MultiTimeField &fields, + const TimeStep &time_step) const + { + timecounter tc; + tc.tic(); + + if (this->enable()) + { + switch (m_scheme) + { + case DynamicType::NEWMARK: + { + std::vector vite, acce; + vite.reserve(msh.cells_size()); + acce.reserve(msh.cells_size()); + const auto vite_prev = fields.getField(-1, FieldName::VITE_CELLS); + const auto acce_prev = fields.getField(-1, FieldName::ACCE_CELLS); + + const auto depl_cells = fields.getCurrentField(FieldName::DEPL_CELLS); + + auto beta = m_param.at("beta"); + auto gamma = m_param.at("gamma"); + scalar_type dt = time_step.increment_time(); + + const scalar_type denom = 1.0 / (beta * dt * dt); + const scalar_type nume = dt * dt / 2.0 * (1.0 - 2.0 * beta); + const scalar_type g0 = (1.0 - gamma) * dt, g1 = gamma * dt; + + for (auto &cl : msh) + { + const auto cell_i = msh.lookup(cl); + + auto aT = denom * depl_cells.at(cell_i) - m_acce_pred.at(cell_i); + auto vT = vite_prev.at(cell_i) + g0 * acce_prev.at(cell_i) + g1 * aT; + + vite.push_back(vT); + acce.push_back(aT); + } + fields.setCurrentField(FieldName::VITE_CELLS, vite); + fields.setCurrentField(FieldName::ACCE_CELLS, acce); + } + default: + break; + } + } + + tc.toc(); + return tc.elapsed(); + } + + void compute(const mesh_type &msh, + const cell_type &cl, + const MeshDegreeInfo °ree_infos, + const vector_type &uTF, + const TimeStep &time_step) + { + // Unsteady Computation + + time_dyna = 0.0; + timecounter tc; + + tc.tic(); + if (this->enable()) + { + const auto cell_i = msh.lookup(cl); + + const auto cell_infos = degree_infos.cellDegreeInfo(msh, cl); + const auto cell_degree = cell_infos.cell_degree(); + + const auto cb = make_vector_monomial_basis(msh, cl, cell_degree); + + const auto faces_infos = cell_infos.facesDegreeInfo(); + + const auto num_cell_dofs = vector_basis_size(cell_degree, mesh_type::dimension, mesh_type::dimension); + const auto num_faces_dofs = vector_faces_dofs(msh, faces_infos); + const auto num_total_dofs = num_cell_dofs + num_faces_dofs; + + const vector_type uT = uTF.head(num_cell_dofs); + + RTF = vector_type::Zero(num_total_dofs); + K_iner = matrix_type::Zero(num_total_dofs, num_total_dofs); + + switch (m_scheme) + { + case DynamicType::NEWMARK: + { + + auto dt = time_step.increment_time(); + auto beta = m_param.at("beta"); + auto rho = m_param.at("rho"); + + const matrix_type mass_mat = rho * make_mass_matrix(msh, cl, cb); + + const auto coeff = 1.0 / (beta * dt * dt); + + K_iner.topLeftCorner(num_cell_dofs, num_cell_dofs) = coeff * mass_mat; + + auto F_iner = coeff * (mass_mat * uT); + RTF.head(num_cell_dofs) -= F_iner; + + RTF.head(num_cell_dofs) += mass_mat * m_acce_pred.at(cell_i); + } + default: + break; + } + } + tc.toc(); + time_dyna += tc.elapsed(); + } + }; + } + +} // end namespace diskpp \ No newline at end of file diff --git a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverParameters.hpp b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverParameters.hpp index e074c5ac..70cc80ce 100755 --- a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverParameters.hpp +++ b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonSolverParameters.hpp @@ -36,8 +36,9 @@ enum StabilizationType : int { HDG = 0, HHO = 1, - NO = 2, - DG = 3 + HHO_SYM = 4, + NO = 2, + DG = 3 }; enum FrictionType : int @@ -47,6 +48,13 @@ enum FrictionType : int COULOMB = 2, }; +enum DynamicType : int +{ + STATIC = 0, + HHT = 1, + NEWMARK = 2, +}; + template class NewtonSolverParameter { @@ -71,7 +79,8 @@ class NewtonSolverParameter bool m_stab; // stabilization yes or no bool m_adapt_stab; // adaptative stabilization - bool m_dynamic; // dynamic or static simulation + DynamicType m_dyna_type; // type of dyna + bool m_dynamic; // dynamic or static simulation std::map m_dyna_para; // list of parameters int m_n_time_save; // number of saving @@ -82,11 +91,10 @@ class NewtonSolverParameter T m_threshold; // threshol for Tesca friction int m_frot_type; // Friction type ? - NewtonSolverParameter() : - m_face_degree(1), m_cell_degree(1), m_grad_degree(1), m_sublevel(5), m_iter_max(20), m_epsilon(T(1E-6)), - m_verbose(false), m_precomputation(false), m_stab(true), m_beta(1), m_stab_type(HHO), m_n_time_save(0), - m_user_end_time(1.0), m_has_user_end_time(false), m_adapt_stab(false), m_dynamic(false), m_theta(1), m_gamma_0(1), - m_threshold(0), m_frot_type(NO_FRICTION) + NewtonSolverParameter() : m_face_degree(1), m_cell_degree(1), m_grad_degree(1), m_sublevel(5), m_iter_max(20), m_epsilon(T(1E-6)), + m_verbose(false), m_precomputation(false), m_stab(true), m_beta(1), m_stab_type(HHO), m_n_time_save(0), + m_user_end_time(1.0), m_has_user_end_time(false), m_adapt_stab(false), m_dynamic(false), m_theta(1), m_gamma_0(1), + m_threshold(0), m_frot_type(NO_FRICTION), m_dyna_type(DynamicType::STATIC) { m_time_step.push_back(std::make_pair(m_user_end_time, 1)); } @@ -228,6 +236,8 @@ class NewtonSolverParameter m_stab_type = HDG; else if (type == "HHO") m_stab_type = HHO; + else if (type == "HHO_SYM") + m_stab_type = HHO_SYM; else if (type == "DG") m_stab_type = DG; else if (type == "NO") @@ -396,6 +406,17 @@ class NewtonSolverParameter return m_dynamic; } + DynamicType getUnsteadyScheme() const + { + return m_dyna_type; + } + + void + setUnsteadyScheme(const DynamicType &scheme) + { + m_dyna_type = scheme; + } + void isUnsteady(const bool dyna) { diff --git a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonStep.hpp b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonStep.hpp index 684afb3e..a1e0bc43 100755 --- a/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonStep.hpp +++ b/libdiskpp/include/diskpp/mechanics/NewtonSolver/NewtonStep.hpp @@ -71,19 +71,12 @@ class NewtonStep typedef vector_boundary_conditions bnd_type; typedef Behavior behavior_type; - std::vector m_displ, m_displ_faces; - std::vector m_velocity, m_acce; - bool m_verbose; bool m_convergence; public: NewtonStep(const param_type& rp) : m_verbose(rp.m_verbose), m_convergence(false) { - m_displ.clear(); - m_displ_faces.clear(); - m_velocity.clear(); - m_acce.clear(); } /** @@ -107,32 +100,6 @@ class NewtonStep m_verbose = v; } - /** - * @brief Initialize the initial guess of the Newton's step with a given guess - * - * @param initial_displ solution \f$ u_T, u_{\partial T} \f$ for each cell - * @param initial_displ_faces solution \f$ u_{F} \f$ for each face - */ - - void - initialize(const std::vector& initial_displ, - const std::vector& initial_displ_faces, - const std::vector& initial_velocity, - const std::vector& initial_acce) - { - m_displ_faces.clear(); - m_displ_faces = initial_displ_faces; - - m_displ.clear(); - m_displ = initial_displ; - - m_velocity.clear(); - m_velocity = initial_velocity; - - m_acce.clear(); - m_acce = initial_acce; - } - /** * @brief Compute the Newton's step until convergence or stopped criterion * @@ -142,18 +109,19 @@ class NewtonStep * @param stab_precomputed contains the precomputed stabilization operators for HHO methods (can be empty) * @return NewtonSolverInfo Informations about the Newton's step during the computation */ - template + template NewtonSolverInfo - compute(const mesh_type& msh, - const bnd_type& bnd, - const param_type& rp, - const MeshDegreeInfo& degree_infos, - const LoadIncrement& lf, - const TimeStep& current_step, - const std::vector& gradient_precomputed, - const std::vector& stab_precomputed, - behavior_type& behavior, - StabCoeffManager& stab_manager) + compute(const mesh_type &msh, + const bnd_type &bnd, + const param_type &rp, + const MeshDegreeInfo °ree_infos, + const LoadIncrement &lf, + const TimeStep ¤t_step, + const std::vector &gradient_precomputed, + const std::vector &stab_precomputed, + behavior_type &behavior, + StabCoeffManager &stab_manager, + MultiTimeField &fields) { NewtonSolverInfo ni; timecounter tc; @@ -162,7 +130,7 @@ class NewtonStep // initialise the NewtonRaphson iteration NewtonIteration newton_iter(msh, bnd, rp, degree_infos, current_step); - newton_iter.initialize(msh, rp, m_displ, m_displ_faces, m_velocity, m_acce); + newton_iter.initialize(msh, fields); m_convergence = false; @@ -173,7 +141,7 @@ class NewtonStep try { assembly_info = newton_iter.assemble( - msh, bnd, rp, degree_infos, lf, gradient_precomputed, stab_precomputed, behavior, stab_manager); + msh, bnd, rp, degree_infos, lf, gradient_precomputed, stab_precomputed, behavior, stab_manager, fields); } catch (const std::invalid_argument& ia) { @@ -201,7 +169,6 @@ class NewtonStep if (m_convergence) { - newton_iter.save_solutions(m_displ, m_displ_faces, m_velocity, m_acce); tc.toc(); ni.m_time_newton = tc.elapsed(); return ni; @@ -211,7 +178,7 @@ class NewtonStep SolveInfo solve_info = newton_iter.solve(); ni.updateSolveInfo(solve_info); // update unknowns - ni.m_assembly_info.m_time_postpro += newton_iter.postprocess(msh, bnd, rp, degree_infos); + ni.m_assembly_info.m_time_postpro += newton_iter.postprocess(msh, bnd, rp, degree_infos, fields); ni.m_iter++; } @@ -232,33 +199,6 @@ class NewtonStep { return m_convergence; } - - /** - * @brief Save solution of the Newton's step - * - */ - void - save_solutions(std::vector& displ, - std::vector& displ_faces, - std::vector& velocity, - std::vector& acce) - { - displ_faces.clear(); - displ_faces = m_displ_faces; - assert(m_displ_faces.size() == displ_faces.size()); - - displ.clear(); - displ = m_displ; - assert(m_displ.size() == displ.size()); - - velocity.clear(); - velocity = m_velocity; - assert(m_velocity.size() == velocity.size()); - - acce.clear(); - acce = m_acce; - assert(m_acce.size() == acce.size()); - } }; } diff --git a/libdiskpp/include/diskpp/mechanics/behaviors/laws/law_qp_bones.hpp b/libdiskpp/include/diskpp/mechanics/behaviors/laws/law_qp_bones.hpp index dda24d92..0da59b7b 100644 --- a/libdiskpp/include/diskpp/mechanics/behaviors/laws/law_qp_bones.hpp +++ b/libdiskpp/include/diskpp/mechanics/behaviors/laws/law_qp_bones.hpp @@ -39,6 +39,15 @@ namespace disk // Bones for the computation of a behavior law at a quadrature point +/* pqrst is there only to squelch -fpermissive + * GCC warnings. sorry for this. */ + +template +using pqrst = point; + +template +using qprst = quadrature_point; + template class law_qp_bones { @@ -51,7 +60,7 @@ class law_qp_bones protected: // coordinat and weight of considered gauss point. - point m_point; + pqrst m_point; scalar_type m_weight; // internal variables at current step @@ -92,13 +101,13 @@ class law_qp_bones { } - quadrature_point + qprst quadrature_point() const { return make_qp(m_point, m_weight); } - point + auto point() const { return m_point; diff --git a/libdiskpp/include/diskpp/mechanics/behaviors/logarithmic_strain/LogarithmicStrain_qp.hpp b/libdiskpp/include/diskpp/mechanics/behaviors/logarithmic_strain/LogarithmicStrain_qp.hpp index c99b90de..3e9c7734 100644 --- a/libdiskpp/include/diskpp/mechanics/behaviors/logarithmic_strain/LogarithmicStrain_qp.hpp +++ b/libdiskpp/include/diskpp/mechanics/behaviors/logarithmic_strain/LogarithmicStrain_qp.hpp @@ -155,13 +155,13 @@ class LogarithmicStrain_qp Pn = static_tensor::Zero(); } - quadrature_point + auto quadrature_point() const { return make_qp(point(), weight()); } - point + auto point() const { return m_law_hpp_qp.point(); diff --git a/libdiskpp/include/diskpp/methods/implementation_hho/curl.hpp b/libdiskpp/include/diskpp/methods/implementation_hho/curl.hpp index 384c96b6..e9c29f12 100644 --- a/libdiskpp/include/diskpp/methods/implementation_hho/curl.hpp +++ b/libdiskpp/include/diskpp/methods/implementation_hho/curl.hpp @@ -1,11 +1,11 @@ /* * DISK++, a template library for DIscontinuous SKeletal methods. - * + * * Matteo Cicuttin (C) 2020 * matteo.cicuttin@uliege.be * * University of Liège - Montefiore Institute - * Applied and Computational Electromagnetics group + * Applied and Computational Electromagnetics group */ @@ -260,7 +260,7 @@ curl_reconstruction(const Mesh& msh, const auto qps = integrate(msh, cl, 2*recdeg); for (auto& qp : qps) { - const auto cphi = rb.eval_curls2(qp.point()); + const auto cphi = rb.eval_curls(qp.point()); stiff += qp.weight() * cphi * cphi.transpose(); } @@ -277,7 +277,7 @@ curl_reconstruction(const Mesh& msh, const auto qps_f = integrate(msh, fc, 2*recdeg); for (auto& qp : qps_f) { - Matrix cphi = rb.eval_curls2(qp.point()); + Matrix cphi = rb.eval_curls(qp.point()); Matrix cphi_n = vcross(cphi, n).block(1,0,rbs-1,1); Matrix f_phi = fb.eval_functions(qp.point()); Matrix c_phi = cb.eval_functions(qp.point()); @@ -330,7 +330,7 @@ curl_reconstruction(const Mesh& msh, const auto qps = integrate(msh, cl, 2*recdeg); for (auto& qp : qps) { - const auto cphi = rb.eval_curls2(qp.point()); + const auto cphi = rb.eval_curls(qp.point()); stiff += qp.weight() * cphi * cphi.transpose(); } @@ -347,7 +347,7 @@ curl_reconstruction(const Mesh& msh, const auto qps_f = integrate(msh, fc, 2*recdeg); for (auto& qp : qps_f) { - Matrix cphi = rb.eval_curls2(qp.point()); + Matrix cphi = rb.eval_curls(qp.point()); Matrix cphi_n = vcross(cphi, n).block(3,0,rbs-3,3); Matrix f_phi = fb.eval_functions(qp.point()); Matrix c_phi = cb.eval_functions(qp.point()); @@ -378,7 +378,7 @@ maxwell_prj(const Mesh& msh, const typename Mesh::cell_type& cl, const auto fbs = vector_basis_size(hdi.face_degree(), Mesh::dimension - 1, Mesh::dimension-1); const auto cbs = vector_basis_size(hdi.cell_degree(), Mesh::dimension, Mesh::dimension); - + const auto fcs = faces(msh, cl); dynamic_vector ret = dynamic_vector::Zero(cbs+fcs.size() * fbs); @@ -449,7 +449,7 @@ curl_reconstruction_pk(const Mesh& msh, const auto r_phi = rb.eval_functions(qp.point()); mass += qp.weight() * r_phi * r_phi.transpose(); - const auto r_cphi = rb.eval_curls2(qp.point()); + const auto r_cphi = rb.eval_curls(qp.point()); const auto c_phi = cb.eval_functions(qp.point()); cr_rhs.block(0, 0, rbs, cbs) += qp.weight() * r_cphi * c_phi.transpose(); } @@ -516,7 +516,7 @@ curl_reconstruction_nedelec(const Mesh& msh, const auto r_phi = rb.eval_functions(qp.point()); mass += qp.weight() * r_phi * r_phi.transpose(); - const auto r_cphi = rb.eval_curls2(qp.point()); + const auto r_cphi = rb.eval_curls(qp.point()); const auto c_phi = cb.eval_functions(qp.point()); cr_rhs.block(0, 0, rbs, cbs) += qp.weight() * r_cphi * c_phi.transpose(); } @@ -772,7 +772,7 @@ curl_hdg_stabilization(const Mesh& msh, abort(); } - rhs.block(0,0,fbs,cbs) = -llt.solve(trace); + rhs.block(0,0,fbs,cbs) = -llt.solve(trace); stab += rhs.transpose() * mass * rhs;// / ht; @@ -847,7 +847,7 @@ curl_hdg_stabilization_nedelec(const Mesh& abort(); } - rhs.block(0,0,fbs,cbs) = -llt.solve(trace); + rhs.block(0,0,fbs,cbs) = -llt.solve(trace); stab += rhs.transpose() * mass * rhs / ht; @@ -921,7 +921,7 @@ curl_Z_stabilization(const Mesh& msh, Matrix n_x_phi_x_n = disk::vcross(n, disk::vcross(phi_tmp, n)); /* Evaluate (1/mur)*(∇×Eᵗ)×n */ - //Matrix cphi_tmp = cb.eval_curls2(qp.point()); + //Matrix cphi_tmp = cb.eval_curls(qp.point()); //Matrix cphi_x_n = (1./mur)*disk::vcross(cphi_tmp, n); Matrix f_phi = fb.eval_functions(qp.point()); @@ -957,7 +957,7 @@ make_curl_curl_matrix(const Mesh& msh, const Element& elem, const Basis& basis, const auto qps = integrate(msh, elem, 2*(degree+di)); for (auto& qp : qps) { - const auto cphi = basis.eval_curls2(qp.point()); + const auto cphi = basis.eval_curls(qp.point()); ret += qp.weight() * cphi * cphi.transpose(); } @@ -1017,12 +1017,12 @@ curl_reconstruction_div(const Mesh& msh, for (auto& qp : qps) { const auto Rphi = rb.eval_functions(qp.point()); - const auto curl_Rphi = rb.eval_curls2(qp.point()); + const auto curl_Rphi = rb.eval_curls(qp.point()); const auto div_Rphi = rb.eval_divergences(qp.point()); - + const auto Dphi = db.eval_functions(qp.point()); const auto grad_Dphi = db.eval_gradients(qp.point()); - + const Matrix Cphi = Rphi.block(0, 0, cbs, 3); const Matrix curl_Cphi = curl_Rphi.block(0, 0, cbs, 3); const Matrix div_Cphi = div_Rphi.segment(0, cb.size()); @@ -1053,7 +1053,7 @@ curl_reconstruction_div(const Mesh& msh, { Matrix phi = rb.eval_functions(qp.point()); Matrix phi_n = vcross(phi, n); - Matrix curl_phi = rb.eval_curls2(qp.point()); + Matrix curl_phi = rb.eval_curls(qp.point()); Matrix curl_phi_n = vcross(curl_phi, n); Matrix f_phi = fb.eval_functions(qp.point()); Matrix c_phi = cb.eval_functions(qp.point()); @@ -1108,7 +1108,7 @@ curl_hho_stabilization(const Mesh& msh, const auto num_faces_dofs = fcs.size() * fbs; matrix_type stab = matrix_type::Zero(cbs + num_faces_dofs, cbs + num_faces_dofs); - + matrix_type Rtrace = matrix_type::Zero(cbs, rbs); auto qps = integrate(msh, cl, 2*recdeg); diff --git a/libdiskpp/include/diskpp/methods/implementation_hho/vector_hho_laplacian.hpp b/libdiskpp/include/diskpp/methods/implementation_hho/vector_hho_laplacian.hpp index efc2e1f9..457c2ad3 100644 --- a/libdiskpp/include/diskpp/methods/implementation_hho/vector_hho_laplacian.hpp +++ b/libdiskpp/include/diskpp/methods/implementation_hho/vector_hho_laplacian.hpp @@ -40,14 +40,27 @@ namespace disk namespace priv { -inline size_t -nb_lag(const size_t dim) -{ - size_t lag = 1; - if (dim == 3) - lag = 3; - return lag; -} + template + inline size_t + nb_lag() + { + static_assert(N != 1); + return 0; + } + + template <> + constexpr size_t + nb_lag<2>() + { + return 1; + } + + template <> + constexpr size_t + nb_lag<3>() + { + return 3; + } } template @@ -77,7 +90,7 @@ make_vector_hho_symmetric_laplacian(const Mesh& msh, const size_t rbs_ho = rbs - N; const size_t num_total_dofs = cbs + num_faces_dofs; - const size_t nb_lag = priv::nb_lag(N); + constexpr size_t nb_lag = priv::nb_lag(); matrix_type stiff = matrix_type::Zero(rbs, rbs); matrix_type gr_lhs = matrix_type::Zero(rbs_ho + nb_lag, rbs_ho + nb_lag); @@ -130,10 +143,43 @@ make_vector_hho_symmetric_laplacian(const Mesh& msh, const auto qps_2 = integrate(msh, cl, recdeg); + const auto lb = make_scalar_monomial_basis(msh, cl, recdeg); + + matrix_type::Zero(rbs, nb_lag); + matrix_type rot = matrix_type::Zero(rbs, nb_lag); for (auto& qp : qps_2) { - const auto rphi = rb.eval_curls(qp.point()); + const function_type dphi = lb.eval_gradients(qp.point()); + Matrix rphi = Matrix::Zero(rbs, N); + + int j = 0; + for (size_t i = 0; i < lb.size(); i++) + { + const auto dphi_i = dphi.row(i); + if (N == 2) + { + rphi(j++) = dphi_i(1); + rphi(j++) = -dphi_i(0); + } + else + { + // row 1 + rphi(j, 0) = dphi_i(1); + rphi(j, 1) = dphi_i(2); + j++; + // row 2 + rphi(j, 0) = -dphi_i(0); + rphi(j, 2) = dphi_i(2); + j++; + // row 3 + rphi(j, 1) = -dphi_i(0); + rphi(j, 2) = -dphi_i(1); + j++; + } + } + assert(rbs == j); + rot += qp.weight() * rphi; } gr_lhs.block(0, rbs_ho, rbs_ho, nb_lag) += rot.bottomLeftCorner(rbs_ho, nb_lag); diff --git a/libdiskpp/include/diskpp/quadratures/bits/quad_raw_tetra.hpp b/libdiskpp/include/diskpp/quadratures/bits/quad_raw_tetra.hpp index dc64e4cf..167a6226 100644 --- a/libdiskpp/include/diskpp/quadratures/bits/quad_raw_tetra.hpp +++ b/libdiskpp/include/diskpp/quadratures/bits/quad_raw_tetra.hpp @@ -10,9 +10,11 @@ #pragma once +#include + +#include "diskpp/common/simplicial_formula.hpp" #include "diskpp/mesh/point.hpp" #include "diskpp/quadratures/quadrature_point.hpp" -#include #include "simplex_gm_rule.hpp" #include "tetrahedron_arbq_rule.hpp" @@ -26,18 +28,9 @@ namespace quadrature namespace priv { -template -inline T -tetra_volume(const Eigen::Matrix& v0, const Eigen::Matrix& v1, const Eigen::Matrix& v2) -{ - return std::abs(v0.dot(v1.cross(v2))) / 6.; -} - -static const double arbq_to_ref_A[3][3] = { - { 0.500000000000000, -0.288675134594813, -0.204124145231932 }, - { 0.000000000000000, 0.577350269189626, -0.204124145231932 }, - { 0.000000000000000, 0.000000000000000, 0.612372435695795 } -}; +static const double arbq_to_ref_A[3][3] = {{0.500000000000000, -0.288675134594813, -0.204124145231932}, + {0.000000000000000, 0.577350269189626, -0.204124145231932}, + {0.000000000000000, 0.000000000000000, 0.612372435695795}}; static const double arbq_to_ref_b[3] = {-1.000000000000000, -0.577350269189626, -0.408248290463863}; @@ -61,7 +54,7 @@ arbq(size_t degree, const point& p0, const point& p1, const point