This is an automated email from the git hooks/post-receive script. logari81 pushed a commit to branch master in repository getfem.
The following commit(s) were added to refs/heads/master by this push: new 752f0399 Whitespace fixes and other minor changes 752f0399 is described below commit 752f03999a6c50af90bf20dc2f6224430e86fab2 Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Tue Dec 26 20:39:58 2023 +0100 Whitespace fixes and other minor changes --- contrib/crack_plate/crack_bilaplacian_problem.cc | 1248 +++++++++++----------- interface/src/gf_linsolve.cc | 4 +- src/gmm/gmm_solver_Schwarz_additive.h | 4 +- tests/schwarz_additive.cc | 10 +- 4 files changed, 635 insertions(+), 631 deletions(-) diff --git a/contrib/crack_plate/crack_bilaplacian_problem.cc b/contrib/crack_plate/crack_bilaplacian_problem.cc index d6fc76d9..f258a788 100644 --- a/contrib/crack_plate/crack_bilaplacian_problem.cc +++ b/contrib/crack_plate/crack_bilaplacian_problem.cc @@ -31,84 +31,82 @@ using std::ends; using std::cin; template <typename T> std::ostream &operator << (std::ostream &o, const std::vector<T>& m) { gmm::write(o,m); return o; } -size_type is_global_dof_type_bis(getfem::pdof_description dof){ -size_type global_dof = 0 ; - for (dim_type d = 0; d < 4 ; ++d){ - if (dof == getfem::global_dof(d)) { - global_dof = 1; - } - } -return global_dof ; +size_type is_global_dof_type_bis(getfem::pdof_description dof) { + size_type global_dof = 0; + for (dim_type d = 0; d < 4; ++d) { + if (dof == getfem::global_dof(d)) { + global_dof = 1; + } + } + return global_dof; } /******** Exact Solution *******************************/ -scalar_type D = 1. ; -scalar_type nu = 0.3 ; -scalar_type AAA = 0.1 ; // mode II -scalar_type BBB = AAA * (3. * nu + 5.)/ (3. * (nu - 1.)) ; -scalar_type DD = 0.0 ; // mode 1 -scalar_type CC = DD * (nu + 7.)/ (3. * (nu - 1.)) ; -scalar_type EE = 0.0 ; // singul 61 -scalar_type FF = 0.0 ; // singul 62 -scalar_type GG = 0.0 ; // singul 63 -scalar_type HH = 0.0 ; //3.0 // singul 6 - -scalar_type P0 = 0.0 ; -scalar_type P1 = 0.0 ; -scalar_type P2 = 0.0 ; - -scalar_type sol_u(const base_node &x){ - scalar_type r = sqrt( x[0] * x[0] + x[1] * x[1] ) ; - //scalar_type theta = 2. * atan( x[1] / ( x[0] + r ) ) ; - scalar_type theta = atan2(x[1], x[0]); - return sqrt(r*r*r)*(AAA*sin(theta/2.0)+BBB*sin(3.0/2.0*theta)+CC*cos(3.0/2.0*theta)+DD*cos(theta/2.0)) + EE * x[1] * (10. * x[1] * x[1]* x[1] + 1.) ; - +scalar_type D = 1.; +scalar_type nu = 0.3; +scalar_type AAA = 0.1; // mode II +scalar_type BBB = AAA * (3. * nu + 5.)/ (3. * (nu - 1.)); +scalar_type DD = 0.0; // mode 1 +scalar_type CC = DD * (nu + 7.)/ (3. * (nu - 1.)); +scalar_type EE = 0.0; // singul 61 +scalar_type FF = 0.0; // singul 62 +scalar_type GG = 0.0; // singul 63 +scalar_type HH = 0.0; //3.0 // singul 6 + +scalar_type P0 = 0.0; +scalar_type P1 = 0.0; +scalar_type P2 = 0.0; + +scalar_type sol_u(const base_node &x) { + scalar_type r = sqrt( x[0] * x[0] + x[1] * x[1] ); + //scalar_type theta = 2. * atan( x[1] / ( x[0] + r ) ); + scalar_type theta = atan2(x[1], x[0]); + return sqrt(r*r*r)*(AAA*sin(theta/2.0)+BBB*sin(3.0/2.0*theta)+CC*cos(3.0/2.0*theta)+DD*cos(theta/2.0)) + EE * x[1] * (10. * x[1] * x[1]* x[1] + 1.); } -scalar_type sol_F(const base_node &) -{return 1. ;//EE * D * 240. ;//256. * cos(2. * x[1]) ; +scalar_type sol_F(const base_node &) { + return 1.; //EE * D * 240.; //256. * cos(2. * x[1]); } void exact_solution_bilap::init(getfem::level_set &ls) { - std::vector<getfem::pglobal_function> cfun(11) ; + std::vector<getfem::pglobal_function> cfun(11); for (unsigned j=0; j < 4; ++j) - cfun[j] = bilaplacian_crack_singular(j, ls, nu, 0.) ; - cfun[4] = bilaplacian_crack_singular(61, ls, nu, 0.) ; - cfun[5] = bilaplacian_crack_singular(62, ls, nu, 0.) ; - cfun[6] = bilaplacian_crack_singular(63, ls, nu, 0.) ; - cfun[7] = bilaplacian_crack_singular(6, ls, nu, 0.) ; - cfun[8] = bilaplacian_crack_singular(10, ls, nu, 0.) ; - cfun[9] = bilaplacian_crack_singular(11, ls, nu, 0.) ; - cfun[10] = bilaplacian_crack_singular(12, ls, nu, 0.) ; + cfun[j] = bilaplacian_crack_singular(j, ls, nu, 0.); + cfun[4] = bilaplacian_crack_singular(61, ls, nu, 0.); + cfun[5] = bilaplacian_crack_singular(62, ls, nu, 0.); + cfun[6] = bilaplacian_crack_singular(63, ls, nu, 0.); + cfun[7] = bilaplacian_crack_singular(6, ls, nu, 0.); + cfun[8] = bilaplacian_crack_singular(10, ls, nu, 0.); + cfun[9] = bilaplacian_crack_singular(11, ls, nu, 0.); + cfun[10] = bilaplacian_crack_singular(12, ls, nu, 0.); mf.set_functions(cfun); U.resize(11); assert(mf.nb_dof() == 11); - U[0] = AAA ; - U[1] = BBB ; - U[2] = CC ; - U[3] = DD ; - U[4] = EE ; - U[5] = FF ; - U[6] = GG ; - U[7] = HH ; - U[8] = P0 ; - U[9] = P1 ; - U[10] = P2 ; - + U[0] = AAA; + U[1] = BBB; + U[2] = CC; + U[3] = DD; + U[4] = EE; + U[5] = FF; + U[6] = GG; + U[7] = HH; + U[8] = P0; + U[9] = P1; + U[10] = P2; } scalar_type eval_fem_gradient_with_finite_differences(getfem::pfem pf, - const base_vector &coeff, - size_type cv, - bgeot::pgeometric_trans pgt, - bgeot::geotrans_inv_convex &gic, - const base_matrix &G, - base_node X0, - scalar_type h, unsigned dg) { + const base_vector &coeff, + size_type cv, + bgeot::pgeometric_trans pgt, + bgeot::geotrans_inv_convex &gic, + const base_matrix &G, + base_node X0, + scalar_type h, unsigned dg) { X0[dg] -= h/2; base_node X0ref; gic.invert(X0, X0ref); getfem::fem_interpolation_context c(pgt, pf, X0ref, G, cv); @@ -127,14 +125,14 @@ scalar_type eval_fem_gradient_with_finite_differences(getfem::pfem pf, } scalar_type eval_fem_hessian_with_finite_differences(getfem::pfem pf, - const base_vector &coeff, - size_type cv, - bgeot::pgeometric_trans pgt, - bgeot::geotrans_inv_convex &gic, - const base_matrix &G, - base_node X0, - scalar_type h, - unsigned dg, unsigned dh) { + const base_vector &coeff, + size_type cv, + bgeot::pgeometric_trans pgt, + bgeot::geotrans_inv_convex &gic, + const base_matrix &G, + base_node X0, + scalar_type h, + unsigned dg, unsigned dh) { X0[dh] -= h/2; scalar_type Gr0 = eval_fem_gradient_with_finite_differences(pf, coeff, cv, pgt, gic, G, X0, h, dg); @@ -146,7 +144,7 @@ scalar_type eval_fem_hessian_with_finite_differences(getfem::pfem pf, } void validate_fem_derivatives(getfem::pfem pf, unsigned cv, - bgeot::pgeometric_trans pgt, const base_matrix &G) { + bgeot::pgeometric_trans pgt, const base_matrix &G) { unsigned N = unsigned(gmm::mat_nrows(G)); scalar_type h = 1e-5; @@ -174,8 +172,8 @@ void validate_fem_derivatives(getfem::pfem pf, unsigned cv, // avoid discontinuity lines in the HCT composite element.. if (gmm::abs(X0ref[0] + X0ref[1] - 1) > 1e-2 && - gmm::abs(X0ref[0] - X0ref[1]) > 1e-2 && - gmm::abs(X0[0]) > 1e-3 && gmm::abs(X0[1])> 1e-3) break; + gmm::abs(X0ref[0] - X0ref[1]) > 1e-2 && + gmm::abs(X0[0]) > 1e-3 && gmm::abs(X0[1])> 1e-3) break; } while (1); //cout << "testing X0 = " << X0 << " (X0ref=" << X0ref << ")\n"; @@ -190,10 +188,10 @@ void validate_fem_derivatives(getfem::pfem pf, unsigned cv, for (unsigned dg = 0; dg < N; ++dg) { grad_fd[dg] = - eval_fem_gradient_with_finite_differences(pf, coeff, cv, pgt, gic, G, X0, h, dg); + eval_fem_gradient_with_finite_differences(pf, coeff, cv, pgt, gic, G, X0, h, dg); for (unsigned dh = 0; dh < N; ++dh) { - hess_fd(0,dg*N+dh) = - eval_fem_hessian_with_finite_differences(pf, coeff, cv, pgt, gic, G, X0, h, dg, dh); + hess_fd(0,dg*N+dh) = + eval_fem_hessian_with_finite_differences(pf, coeff, cv, pgt, gic, G, X0, h, dg, dh); } } @@ -203,11 +201,11 @@ void validate_fem_derivatives(getfem::pfem pf, unsigned cv, gmm::vect_dist2((base_vector&)hess, (base_vector&)hess_fd); if (err_grad > 1e-4 || - err_hess > 1e-4) { + err_hess > 1e-4) { cout << "validate_fem_derivatives dof=" << idof << "/" << pf->nb_dof(cv) << " -- X0ref = " << X0ref << "\n"; if (gmm::vect_dist2((base_vector&)grad, (base_vector&)grad_fd) > 1e-4) - cout << "grad = " << (base_vector&)grad << "\ngrad_fd = " << (base_vector&)grad_fd << "\n"; + cout << "grad = " << (base_vector&)grad << "\ngrad_fd = " << (base_vector&)grad_fd << "\n"; cout << "hess = " << (base_vector&)hess << "\nhess_fd = " << (base_vector&)hess_fd << "\n"; if (err_grad + err_hess > 1.0) { cout << "---------> COMPLETEMENT FAUX!!!!!!!!!!!!!!!!!!!!!!!!!!!\n"; abort(); } } @@ -243,79 +241,79 @@ void bilaplacian_crack_problem::init(void) { std::string MESH_TYPE = PARAM.string_value("MESH_TYPE","Mesh type "); std::string FEM_TYPE = PARAM.string_value("FEM_TYPE","FEM name"); std::string INTEGRATION = PARAM.string_value("INTEGRATION", - "Name of integration method"); + "Name of integration method"); cout << "MESH_TYPE=" << MESH_TYPE << "\n"; cout << "FEM_TYPE=" << FEM_TYPE << "\n"; cout << "INTEGRATION=" << INTEGRATION << "\n"; std::string SIMPLEX_INTEGRATION = PARAM.string_value("SIMPLEX_INTEGRATION", - "Name of simplex integration method"); + "Name of simplex integration method"); std::string SINGULAR_INTEGRATION = PARAM.string_value("SINGULAR_INTEGRATION"); enrichment_option = unsigned(PARAM.int_value("ENRICHMENT_OPTION", - "Enrichment option")); + "Enrichment option")); enr_area_radius = PARAM.real_value("RADIUS_ENR_AREA", - "radius of the enrichment area"); + "radius of the enrichment area"); - sol_ref = PARAM.int_value("SOL_REF") ; - scalar_type angle_rot = PARAM.real_value("ANGLE_ROT") ; - size_type N = 2 ; + sol_ref = PARAM.int_value("SOL_REF"); + scalar_type angle_rot = PARAM.real_value("ANGLE_ROT"); + size_type N = 2; /* First step : build the mesh */ if (!MESH_FILE.empty()) { - mesh.read_from_file(MESH_FILE); - // printing number of elements - cout << "Number of element of the mesh: " << mesh.convex_index().card() << "\n" ; - base_small_vector tt(N); - tt[0] = PARAM.real_value("TRANSLAT_X") ; - tt[1] = PARAM.real_value("TRANSLAT_Y") ; - if (sol_ref == 1 || sol_ref == 2){ - tt[0] -= PARAM.real_value("CRACK_SEMI_LENGTH") ; - } - cout << "TRANSLAT_X = " << tt[0] << " ; TRANSLAT_Y = " << tt[1] << "\n" ; - mesh.translation(tt); - MESH_TYPE = bgeot::name_of_geometric_trans - (mesh.trans_of_convex(mesh.convex_index().first_true())); - bgeot::pgeometric_trans pgt = - bgeot::geometric_trans_descriptor(MESH_TYPE); - cout << "MESH_TYPE=" << MESH_TYPE << "\n"; - N = mesh.dim(); + mesh.read_from_file(MESH_FILE); + // printing number of elements + cout << "Number of element of the mesh: " << mesh.convex_index().card() << "\n"; + base_small_vector tt(N); + tt[0] = PARAM.real_value("TRANSLAT_X"); + tt[1] = PARAM.real_value("TRANSLAT_Y"); + if (sol_ref == 1 || sol_ref == 2) { + tt[0] -= PARAM.real_value("CRACK_SEMI_LENGTH"); + } + cout << "TRANSLAT_X = " << tt[0] << " ; TRANSLAT_Y = " << tt[1] << "\n"; + mesh.translation(tt); + MESH_TYPE = bgeot::name_of_geometric_trans + (mesh.trans_of_convex(mesh.convex_index().first_true())); + bgeot::pgeometric_trans pgt = + bgeot::geometric_trans_descriptor(MESH_TYPE); + cout << "MESH_TYPE=" << MESH_TYPE << "\n"; + N = mesh.dim(); } else { - bgeot::pgeometric_trans pgt = - bgeot::geometric_trans_descriptor(MESH_TYPE); - N = pgt->dim(); - GMM_ASSERT1(N == 2, "For a plate problem, N should be 2"); - std::vector<size_type> nsubdiv(N); - NX = PARAM.int_value("NX", "Number of space steps ") ; - std::fill(nsubdiv.begin(),nsubdiv.end(), NX); - if (sol_ref == 1) nsubdiv[0] = NX / 2 ; - if (sol_ref == 2) { - size_type NY = PARAM.int_value("NY") ; - nsubdiv[1] = NY ; - } - getfem::regular_unit_mesh(mesh, nsubdiv, pgt, PARAM.int_value("MESH_NOISED") != 0); + bgeot::pgeometric_trans pgt = + bgeot::geometric_trans_descriptor(MESH_TYPE); + N = pgt->dim(); + GMM_ASSERT1(N == 2, "For a plate problem, N should be 2"); + std::vector<size_type> nsubdiv(N); + NX = PARAM.int_value("NX", "Number of space steps "); + std::fill(nsubdiv.begin(),nsubdiv.end(), NX); + if (sol_ref == 1) nsubdiv[0] = NX / 2; + if (sol_ref == 2) { + size_type NY = PARAM.int_value("NY"); + nsubdiv[1] = NY; + } + getfem::regular_unit_mesh(mesh, nsubdiv, pgt, PARAM.int_value("MESH_NOISED") != 0); /* scale the unit mesh to [LX,LY,..] and incline it */ base_small_vector tt(N); switch (sol_ref) { case 0 : { - tt[0] = - 0.5 + PARAM.real_value("TRANSLAT_X") ; - tt[1] = - 0.5 + PARAM.real_value("TRANSLAT_Y") ; - } break ; + tt[0] = - 0.5 + PARAM.real_value("TRANSLAT_X"); + tt[1] = - 0.5 + PARAM.real_value("TRANSLAT_Y"); + } break; case 1 : { - tt[0] = - PARAM.real_value("CRACK_SEMI_LENGTH") ; - tt[1] = - 0.5 + PARAM.real_value("TRANSLAT_Y") ; - } break ; + tt[0] = - PARAM.real_value("CRACK_SEMI_LENGTH"); + tt[1] = - 0.5 + PARAM.real_value("TRANSLAT_Y"); + } break; case 2 : { - tt[0] = - PARAM.real_value("CRACK_SEMI_LENGTH"); - tt[1] = - 0.5 + PARAM.real_value("TRANSLAT_Y") ; - } break ; + tt[0] = - PARAM.real_value("CRACK_SEMI_LENGTH"); + tt[1] = - 0.5 + PARAM.real_value("TRANSLAT_Y"); + } break; case 3 : { - tt[0] = - 0.5 + PARAM.real_value("TRANSLAT_X") ; - tt[1] = - 0.5 + PARAM.real_value("TRANSLAT_Y") ; - } break ; + tt[0] = - 0.5 + PARAM.real_value("TRANSLAT_X"); + tt[1] = - 0.5 + PARAM.real_value("TRANSLAT_Y"); + } break; default : - GMM_ASSERT1(false, "SOL_REF parameter is undefined"); - break ; + GMM_ASSERT1(false, "SOL_REF parameter is undefined"); + break; } mesh.translation(tt); bgeot::base_matrix M(N,N); @@ -323,44 +321,46 @@ void bilaplacian_crack_problem::init(void) { static const char *t[] = {"LX","LY","LZ"}; M(i,i) = (i<3) ? PARAM.real_value(t[i],t[i]) : 1.0; } - if (sol_ref == 2){ - M(0, 0) = 1.5 ; - M(1, 1) = 1.0 ; + if (sol_ref == 2) { + M(0, 0) = 1.5; + M(1, 1) = 1.0; } - if (sol_ref == 0 || sol_ref == 3){ - M(0, 0) = cos(angle_rot) ; M(1, 1) = M(0, 0) ; - M(1, 0) = sin(angle_rot) ; M(0, 1) = - M(1, 0) ; + if (sol_ref == 0 || sol_ref == 3) { + M(0, 0) = cos(angle_rot); + M(1, 1) = M(0, 0); + M(1, 0) = sin(angle_rot); + M(0, 1) = - M(1, 0); } mesh.transformation(M); } - scalar_type quality = 1.0, avg_area = 0. , min_area = 1. , max_area = 0., area ; - scalar_type radius, avg_radius = 0., min_radius = 1., max_radius = 0. ; - size_type cpt = 0 ; - for (dal::bv_visitor i(mesh.convex_index()); !i.finished(); ++i){ - quality = std::min(quality, mesh.convex_quality_estimate(i)); - area = mesh.convex_area_estimate(i) ; - radius = mesh.convex_radius_estimate(i) ; - avg_area += area ; - avg_radius += radius ; - min_radius = std::min(radius, min_radius) ; - max_radius = std::max(radius, max_radius) ; - min_area = std::min(min_area, area) ; - max_area = std::max(max_area, area) ; - cpt++ ; + scalar_type quality = 1.0, avg_area = 0. , min_area = 1. , max_area = 0., area; + scalar_type radius, avg_radius = 0., min_radius = 1., max_radius = 0.; + size_type cpt = 0; + for (dal::bv_visitor i(mesh.convex_index()); !i.finished(); ++i) { + quality = std::min(quality, mesh.convex_quality_estimate(i)); + area = mesh.convex_area_estimate(i); + radius = mesh.convex_radius_estimate(i); + avg_area += area; + avg_radius += radius; + min_radius = std::min(radius, min_radius); + max_radius = std::max(radius, max_radius); + min_area = std::min(min_area, area); + max_area = std::max(max_area, area); + cpt++; } - avg_area /= scalar_type(cpt) ; - avg_radius /= scalar_type(cpt) ; + avg_area /= scalar_type(cpt); + avg_radius /= scalar_type(cpt); /* cout << "quality of mesh : " << quality << endl; cout << "average radius : " << avg_radius << endl; cout << "radius min : " << min_radius << " ; radius max : " << max_radius << endl; - cout << "average area : " << avg_area << endl ; + cout << "average area : " << avg_area << endl; cout << "area min : " << min_area << " ; area max : " << max_area << endl; */ /* read the parameters */ - epsilon = PARAM.real_value("EPSILON", "thickness") ; - nu = PARAM.real_value("NU", "nu") ; - D = PARAM.real_value("D", "Flexure modulus") ; + epsilon = PARAM.real_value("EPSILON", "thickness"); + nu = PARAM.real_value("NU", "nu"); + D = PARAM.real_value("D", "Flexure modulus"); int mv = int(PARAM.int_value("MORTAR_VERSION", "Mortar version") ); int cv = int(PARAM.int_value("CLOSING_VERSION") ); mortar_version = size_type(mv); @@ -383,8 +383,9 @@ void bilaplacian_crack_problem::init(void) { getfem::int_method_descriptor(INTEGRATION); getfem::pintegration_method sppi = getfem::int_method_descriptor(SIMPLEX_INTEGRATION); - getfem::pintegration_method sing_ppi = (SINGULAR_INTEGRATION.size() ? - getfem::int_method_descriptor(SINGULAR_INTEGRATION) : 0); + getfem::pintegration_method sing_ppi = + (SINGULAR_INTEGRATION.size() ? getfem::int_method_descriptor(SINGULAR_INTEGRATION) + : 0); mim.set_integration_method(mesh.convex_index(), ppi); mls.add_level_set(ls); @@ -392,7 +393,8 @@ void bilaplacian_crack_problem::init(void) { /* Setting the finite element on the mf_u */ mf_pre_u.set_finite_element(mesh.convex_index(), pf_u); - getfem::pfem pf_partition_of_unity = getfem::fem_descriptor(PARAM.string_value("PARTITION_OF_UNITY_FEM_TYPE")) ; + getfem::pfem pf_partition_of_unity = + getfem::fem_descriptor(PARAM.string_value("PARTITION_OF_UNITY_FEM_TYPE")); mf_partition_of_unity.set_finite_element(mesh.convex_index(), pf_partition_of_unity); // set the mesh_fem of the multipliers (for the dirichlet condition) @@ -402,7 +404,7 @@ void bilaplacian_crack_problem::init(void) { else { cout << "DIRICHLET_FEM_TYPE=" << dirichlet_fem_name << "\n"; mf_mult.set_finite_element(mesh.convex_index(), - getfem::fem_descriptor(dirichlet_fem_name)); + getfem::fem_descriptor(dirichlet_fem_name)); } std::string dirichlet_der_fem_name = PARAM.string_value("DIRICHLET_DER_FEM_TYPE", ""); // for the dirichlet condition on the derivative @@ -411,7 +413,7 @@ void bilaplacian_crack_problem::init(void) { else { cout << "DIRICHLET_DER_FEM_TYPE=" << dirichlet_der_fem_name << "\n"; mf_mult_d.set_finite_element(mesh.convex_index(), - getfem::fem_descriptor(dirichlet_der_fem_name)); + getfem::fem_descriptor(dirichlet_der_fem_name)); } /* set the finite element on mf_rhs (same as mf_u if DATA_FEM_TYPE is @@ -419,12 +421,12 @@ void bilaplacian_crack_problem::init(void) { std::string data_fem_name = PARAM.string_value("DATA_FEM_TYPE"); if (data_fem_name.size() == 0) { GMM_ASSERT1(pf_u->is_lagrange(), "You are using a non-lagrange FEM. " - << "In that case you need to set " - << "DATA_FEM_TYPE in the .param file"); + << "In that case you need to set " + << "DATA_FEM_TYPE in the .param file"); mf_rhs.set_finite_element(mesh.convex_index(), pf_u); } else { mf_rhs.set_finite_element(mesh.convex_index(), - getfem::fem_descriptor(data_fem_name)); + getfem::fem_descriptor(data_fem_name)); } // set the mesh_fem for the multipliers for the case the Integral Matching @@ -437,93 +439,93 @@ void bilaplacian_crack_problem::init(void) { cout << "Selecting Neumann and Dirichlet boundaries\n"; getfem::mesh_region border_faces; getfem::outer_faces_of_mesh(mesh, border_faces); - if (sol_ref == 0 && angle_rot == 0.){ - for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { - base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f()); - un /= gmm::vect_norm2(un); - if ( gmm::abs(un[0]) >= 0.9999 ) { // vertical edges - mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); - mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); - } - else{ // horizontal edges - unsigned id_point_1_of_face, id_point_2_of_face ; - scalar_type x1, x2 ; - id_point_1_of_face = mesh.structure_of_convex(i.cv())->ind_points_of_face(i.f())[0] ; - id_point_2_of_face = mesh.structure_of_convex(i.cv())->ind_points_of_face(i.f())[1] ; - x1 = mesh.points_of_convex(i.cv())[id_point_1_of_face][1] ; - x2 = mesh.points_of_convex(i.cv())[id_point_2_of_face][1] ; - if ( gmm::abs(x1) > 0.4999 && gmm::abs(x2) > 0.4999 ) {// on the boundary => clamped - mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); - mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); - } - else { // on the crack => free boundary condition - mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f()); - mesh.region(FORCE_BOUNDARY_NUM).add(i.cv(), i.f()); - } - } + if (sol_ref == 0 && angle_rot == 0.) { + for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { + base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f()); + un /= gmm::vect_norm2(un); + if (gmm::abs(un[0]) >= 0.9999) { // vertical edges + mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); + mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); + } + else { // horizontal edges + unsigned id_point_1_of_face, id_point_2_of_face; + scalar_type x1, x2; + id_point_1_of_face = mesh.structure_of_convex(i.cv())->ind_points_of_face(i.f())[0]; + id_point_2_of_face = mesh.structure_of_convex(i.cv())->ind_points_of_face(i.f())[1]; + x1 = mesh.points_of_convex(i.cv())[id_point_1_of_face][1]; + x2 = mesh.points_of_convex(i.cv())[id_point_2_of_face][1]; + if (gmm::abs(x1) > 0.4999 && gmm::abs(x2) > 0.4999) { // on the boundary => clamped + mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); + mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); + } + else { // on the crack => free boundary condition + mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f()); + mesh.region(FORCE_BOUNDARY_NUM).add(i.cv(), i.f()); + } } - } - - - if (sol_ref == 1 ){ - for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { - base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f()); - un /= gmm::vect_norm2(un); - if ( (un[0] <= -0.9999) || (un[0] >= 0.9999) ) { - //if ( -un[0] >= 0.999 ) { - mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); - - } - if (gmm::abs(un[1]) >= 0.999) { - mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); - mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f()); - } - } + } + } + + + if (sol_ref == 1 ) { + for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { + base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f()); + un /= gmm::vect_norm2(un); + if (un[0] <= -0.9999 || un[0] >= 0.9999) { + //if (-un[0] >= 0.999) { + mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); + } + if (gmm::abs(un[1]) >= 0.999) { + mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); + mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f()); + } + } } - if (sol_ref == 2 ){ - for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { - base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f()); - un /= gmm::vect_norm2(un); - if (un[0] < - 0.9999) { // symetry condition - mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); - } - else{ - mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); - mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f()); - } - } + if (sol_ref == 2 ) { + for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { + base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f()); + un /= gmm::vect_norm2(un); + if (un[0] < - 0.9999) { // symetry condition + mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); + } + else { + mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); + mesh.region(MOMENTUM_BOUNDARY_NUM).add(i.cv(), i.f()); + } } + } - if (sol_ref == 0 && angle_rot != 0. ){ + if (sol_ref == 0 && angle_rot != 0. ) { // does not support the case of a conformal mesh around the crack // (all the boundary, including the crack faces, will be clamped). - for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { - mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); - mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); - } - } + for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { + mesh.region(CLAMPED_BOUNDARY_NUM).add(i.cv(), i.f()); + mesh.region(SIMPLE_SUPPORT_BOUNDARY_NUM).add(i.cv(), i.f()); + } + } exact_sol.init(ls); - cout << "initialisation de la level-set : \n" ; + cout << "initialisation de la level-set : \n"; // Setting the level-set ls.reinit(); - // scalar_type a = PARAM.real_value("CRACK_SEMI_LENGTH") ; + // scalar_type a = PARAM.real_value("CRACK_SEMI_LENGTH"); for (size_type d = 0; d < ls.get_mesh_fem().nb_dof(); ++d) { scalar_type x = ls.get_mesh_fem().point_of_basic_dof(d)[0]; scalar_type y = ls.get_mesh_fem().point_of_basic_dof(d)[1]; - if (sol_ref == 0){ - ls.values(0)[d] = y ; // + 1/4.*(x + .25); - ls.values(1)[d] = x ;} - if (sol_ref == 1){ - ls.values(0)[d] = y ; - ls.values(1)[d] = x ; //x * x - a * a ; + if (sol_ref == 0) { + ls.values(0)[d] = y; // + 1/4.*(x + .25); + ls.values(1)[d] = x; + } + if (sol_ref == 1) { + ls.values(0)[d] = y; + ls.values(1)[d] = x; //x * x - a * a; } - if (sol_ref == 2){ // to modify if rotation is supported - ls.values(0)[d] = y ; - ls.values(1)[d] = x ; //x * x - a * a ; + if (sol_ref == 2) { // to modify if rotation is supported + ls.values(0)[d] = y; + ls.values(1)[d] = x; //x * x - a * a; } } //ls.simplify(0.5); @@ -551,23 +553,23 @@ void bilaplacian_crack_problem::init(void) { // << "H2 error = " << getfem::asm_H2_norm(mim, mf_rhs, V) << endl // /*<< "Linfty error = " << gmm::vect_norminf(V) << endl*/; // cout << "semi-norme H1 = " << getfem::asm_H1_semi_norm(mim, mf_rhs, V) << endl -// << "semi-norme H2 = " << getfem::asm_H2_semi_norm(mim, mf_rhs, V) << endl ; +// << "semi-norme H2 = " << getfem::asm_H2_semi_norm(mim, mf_rhs, V) << endl; // // } /* compute the error with respect to the exact solution */ void bilaplacian_crack_problem::compute_error(plain_vector &U) { - if (PARAM.real_value("RADIUS_SPLIT_DOMAIN") == 0){ - plain_vector V(gmm::vect_size(U)) ; - gmm::clear(V) ; + if (PARAM.real_value("RADIUS_SPLIT_DOMAIN") == 0) { + plain_vector V(gmm::vect_size(U)); + gmm::clear(V); cout << "L2 ERROR:" << getfem::asm_L2_dist(mim, mf_u(), U, - exact_sol.mf, exact_sol.U) << "\n"; + exact_sol.mf, exact_sol.U) << "\n"; cout << "H1 ERROR:" << getfem::asm_H1_dist(mim, mf_u(), U, - exact_sol.mf, exact_sol.U) << "\n"; + exact_sol.mf, exact_sol.U) << "\n"; cout << "H2 ERROR:" << getfem::asm_H2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U) << "\n"; @@ -577,59 +579,61 @@ void bilaplacian_crack_problem::compute_error(plain_vector &U) { // cout << "SEMI H2 ERROR:" // << getfem::asm_H2_semi_dist(mim, mf_u(), U, // exact_sol.mf, exact_sol.U) << "\n"; - if ( PARAM.int_value("NORM_EXACT") ){ + if (PARAM.int_value("NORM_EXACT")) { cout << "L2 exact:" << getfem::asm_L2_dist(mim, mf_u(), V, - exact_sol.mf, exact_sol.U) << "\n"; + exact_sol.mf, exact_sol.U) << "\n"; cout << "H1 exact:" << getfem::asm_H1_dist(mim, mf_u(), V, - exact_sol.mf, exact_sol.U) << "\n"; + exact_sol.mf, exact_sol.U) << "\n"; cout << "H2 exact:" << getfem::asm_H2_dist(mim, mf_u(), V, exact_sol.mf, exact_sol.U) << "\n"; } } else { - getfem::mesh_region r_center, r_ext ; - scalar_type radius_split_domain = PARAM.real_value("RADIUS_SPLIT_DOMAIN") ; - bool in_area ; - for (dal::bv_visitor cv(mesh.convex_index()) ; !cv.finished() ; ++cv){ - in_area = true; - /* For each element, we test all of its nodes. - If all the nodes are inside the enrichment area, - then the element is completly inside the area too */ - for (unsigned j=0; j < mesh.nb_points_of_convex(cv); ++j) { - if (gmm::sqr(mesh.points_of_convex(cv)[j][0] ) + - gmm::sqr(mesh.points_of_convex(cv)[j][1] ) > - gmm::sqr(radius_split_domain)) - in_area = false; - break; - } - if (in_area) r_center.add(cv) ; - else r_ext.add(cv) ; - } - scalar_type L2_center, H1_center, H2_center; - cout << "ERROR SPLITTED - RADIUS = " << radius_split_domain << "\n"; - cout << "Error on the crack tip zone:\n" ; - L2_center = getfem::asm_L2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_center) ; - cout << " L2 error:" << L2_center << "\n"; - H1_center = getfem::asm_H1_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_center) ; - cout << " H1 error:" << H1_center << "\n"; - H2_center = getfem::asm_H2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_center) ; - cout << " H2 error:" << H2_center << "\n"; - - cout << "Error on the remaining part of the domain:\n"; - scalar_type L2_ext, H1_ext, H2_ext; - L2_ext = getfem::asm_L2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_ext) ; - cout << " L2 error:" << L2_ext << "\n"; - H1_ext = getfem::asm_H1_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_ext) ; - cout << " H1 error:" << H1_ext << "\n"; - H2_ext = getfem::asm_H2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_ext) ; - cout << " H2 error:" << H2_ext << "\n"; - - cout << "Error on the hole domain:\n"; - cout << "L2 ERROR:" - << gmm::sqrt( gmm::sqr(L2_center) + gmm::sqr(L2_ext) ) << "\n"; + getfem::mesh_region r_center, r_ext; + scalar_type radius_split_domain = PARAM.real_value("RADIUS_SPLIT_DOMAIN"); + bool in_area; + for (dal::bv_visitor cv(mesh.convex_index()); !cv.finished(); ++cv) { + in_area = true; + /* For each element, we test all of its nodes. + If all the nodes are inside the enrichment area, + then the element is completly inside the area too */ + for (unsigned j=0; j < mesh.nb_points_of_convex(cv); ++j) { + if (gmm::sqr(mesh.points_of_convex(cv)[j][0] ) + + gmm::sqr(mesh.points_of_convex(cv)[j][1] ) > + gmm::sqr(radius_split_domain)) + in_area = false; + break; + } + if (in_area) + r_center.add(cv); + else + r_ext.add(cv); + } + scalar_type L2_center, H1_center, H2_center; + cout << "ERROR SPLITTED - RADIUS = " << radius_split_domain << "\n"; + cout << "Error on the crack tip zone:\n"; + L2_center = getfem::asm_L2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_center); + cout << " L2 error:" << L2_center << "\n"; + H1_center = getfem::asm_H1_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_center); + cout << " H1 error:" << H1_center << "\n"; + H2_center = getfem::asm_H2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_center); + cout << " H2 error:" << H2_center << "\n"; + + cout << "Error on the remaining part of the domain:\n"; + scalar_type L2_ext, H1_ext, H2_ext; + L2_ext = getfem::asm_L2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_ext); + cout << " L2 error:" << L2_ext << "\n"; + H1_ext = getfem::asm_H1_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_ext); + cout << " H1 error:" << H1_ext << "\n"; + H2_ext = getfem::asm_H2_dist(mim, mf_u(), U, exact_sol.mf, exact_sol.U, r_ext); + cout << " H2 error:" << H2_ext << "\n"; + + cout << "Error on the hole domain:\n"; + cout << "L2 ERROR:" + << gmm::sqrt( gmm::sqr(L2_center) + gmm::sqr(L2_ext) ) << "\n"; cout << "H1 ERROR:" << gmm::sqrt( gmm::sqr(H1_center) + gmm::sqr(H1_ext) ) << "\n"; @@ -647,110 +651,109 @@ void bilaplacian_crack_problem::compute_error(plain_vector &U) { bool bilaplacian_crack_problem::solve(plain_vector &U) { - if (enrichment_option == 2 || enrichment_option == 3){ - cout << "setting singularities \n" ; - if (PARAM.int_value("SING_BASE_TYPE") == 0){ - std::vector<getfem::pglobal_function> ufunc(4); - for (size_type i = 0 ; i < ufunc.size() ; ++i) { - ufunc[i] = bilaplacian_crack_singular(i, ls, nu, 0.); - } - mf_sing_u.set_functions(ufunc); - } - if (PARAM.int_value("SING_BASE_TYPE") == 1){ - std::vector<getfem::pglobal_function> ufunc(2); - for (size_type i = 0 ; i < ufunc.size() ; ++i) { - ufunc[i] = bilaplacian_crack_singular(i + 4, ls, nu, 0.); - } - mf_sing_u.set_functions(ufunc); - } + if (enrichment_option == 2 || enrichment_option == 3) { + cout << "setting singularities \n"; + if (PARAM.int_value("SING_BASE_TYPE") == 0) { + std::vector<getfem::pglobal_function> ufunc(4); + for (size_type i = 0; i < ufunc.size(); ++i) { + ufunc[i] = bilaplacian_crack_singular(i, ls, nu, 0.); + } + mf_sing_u.set_functions(ufunc); + } + if (PARAM.int_value("SING_BASE_TYPE") == 1) { + std::vector<getfem::pglobal_function> ufunc(2); + for (size_type i = 0; i < ufunc.size(); ++i) { + ufunc[i] = bilaplacian_crack_singular(i + 4, ls, nu, 0.); + } + mf_sing_u.set_functions(ufunc); + } } - else if (enrichment_option == 4){ - cout << "Setting up the singular functions for the cutoff enrichment\n"; - if (PARAM.int_value("SING_BASE_TYPE") == 0) { - std::vector<getfem::pglobal_function> vfunc(4); - for (size_type i = 0; i < vfunc.size(); ++i) { - /* use the singularity */ - getfem::pxy_function - s = std::make_shared<crack_singular_bilaplacian_xy_function> - (i, ls, nu, 0.); - /* use the product of the singularity function - with a cutoff */ - getfem::pxy_function - cc = std::make_shared<getfem::cutoff_xy_function> - (int(cutoff.fun_num), cutoff.radius, cutoff.radius1, - cutoff.radius0); - s = std::make_shared<getfem::product_of_xy_functions>(s, cc); - vfunc[i] = getfem::global_function_on_level_set(ls, s); - } - mf_sing_u.set_functions(vfunc); - } - if (PARAM.int_value("SING_BASE_TYPE") == 1){ - std::vector<getfem::pglobal_function> vfunc(2); - for (size_type i = 0; i < vfunc.size(); ++i) { - /* use the singularity */ - getfem::pxy_function - s = std::make_shared<crack_singular_bilaplacian_xy_function> - (i+4, ls, nu, 0.); - /* use the product of the singularity function - with a cutoff */ - getfem::pxy_function - cc = std::make_shared<getfem::cutoff_xy_function> - (int(cutoff.fun_num), cutoff.radius, cutoff.radius1, - cutoff.radius0); - s = std::make_shared<getfem::product_of_xy_functions>(s, cc); - vfunc[i] = getfem::global_function_on_level_set(ls, s); - } - mf_sing_u.set_functions(vfunc); - } + else if (enrichment_option == 4) { + cout << "Setting up the singular functions for the cutoff enrichment\n"; + if (PARAM.int_value("SING_BASE_TYPE") == 0) { + std::vector<getfem::pglobal_function> vfunc(4); + for (size_type i = 0; i < vfunc.size(); ++i) { + /* use the singularity */ + getfem::pxy_function + s = std::make_shared<crack_singular_bilaplacian_xy_function> + (i, ls, nu, 0.); + /* use the product of the singularity function + with a cutoff */ + getfem::pxy_function + cc = std::make_shared<getfem::cutoff_xy_function> + (int(cutoff.fun_num), cutoff.radius, cutoff.radius1, cutoff.radius0); + s = std::make_shared<getfem::product_of_xy_functions>(s, cc); + vfunc[i] = getfem::global_function_on_level_set(ls, s); + } + mf_sing_u.set_functions(vfunc); + } + if (PARAM.int_value("SING_BASE_TYPE") == 1) { + std::vector<getfem::pglobal_function> vfunc(2); + for (size_type i = 0; i < vfunc.size(); ++i) { + /* use the singularity */ + getfem::pxy_function + s = std::make_shared<crack_singular_bilaplacian_xy_function> + (i+4, ls, nu, 0.); + /* use the product of the singularity function + with a cutoff */ + getfem::pxy_function + cc = std::make_shared<getfem::cutoff_xy_function> + (int(cutoff.fun_num), cutoff.radius, cutoff.radius1, cutoff.radius0); + s = std::make_shared<getfem::product_of_xy_functions>(s, cc); + vfunc[i] = getfem::global_function_on_level_set(ls, s); + } + mf_sing_u.set_functions(vfunc); + } } // Setting the enrichment --------------------------------------------/ switch(enrichment_option) { case -1: // classical FEM { - mf_u_sum.set_mesh_fems(mf_pre_u); + mf_u_sum.set_mesh_fems(mf_pre_u); } - break ; + break; case 0 : // No enrichment { - mf_u_sum.set_mesh_fems(mfls_u); + mf_u_sum.set_mesh_fems(mfls_u); // an optional treatment : exporting a representation of the mesh - if (!PARAM.int_value("MIXED_ELEMENTS")){ - getfem::mesh_fem mf_enrich(mesh); - getfem::pfem pf_mef = getfem::classical_fem(mesh.trans_of_convex(mesh.convex_index().first_true()), 1 ); - mf_enrich.set_finite_element(mesh.convex_index(), pf_mef) ; - std::vector<scalar_type> UU(mf_enrich.nb_dof()) ; - std::fill(UU.begin(), UU.end() ,0.) ; - cout << "exporting mesh to " << "mesh_representation.vtk" << "..\n"; - getfem::vtk_export exp("mesh_representation.vtk", false); - exp.exporting(mf_enrich); - exp.write_point_data(mf_enrich, UU, "mesh"); - cout << "export done, you can view the data file with (for example)\n" - "mayavi -d mesh_representation.vtk -f " - "WarpScalar -m BandedSurfaceMap -m Outline\n";} + if (!PARAM.int_value("MIXED_ELEMENTS")) { + getfem::mesh_fem mf_enrich(mesh); + getfem::pfem pf_mef = getfem::classical_fem(mesh.trans_of_convex(mesh.convex_index().first_true()), 1 ); + mf_enrich.set_finite_element(mesh.convex_index(), pf_mef); + std::vector<scalar_type> UU(mf_enrich.nb_dof()); + std::fill(UU.begin(), UU.end() ,0.); + cout << "exporting mesh to " << "mesh_representation.vtk" << "..\n"; + getfem::vtk_export exp("mesh_representation.vtk", false); + exp.exporting(mf_enrich); + exp.write_point_data(mf_enrich, UU, "mesh"); + cout << "export done, you can view the data file with (for example)\n" + "mayavi -d mesh_representation.vtk -f " + "WarpScalar -m BandedSurfaceMap -m Outline\n"; + } } - break ; + break; case 1 : { cout << "\npointwise matching\n"; - /* first : selecting the convexes that are completly included in the enrichment area */ - for (dal::bv_visitor i(mesh.convex_index()); !i.finished(); ++i) { - pm_convexes.add(i) ; - /* For each element, we test all of its nodes. - If all the nodes are inside the enrichment area, - then the element is completly inside the area too */ - for (unsigned j=0; j < mesh.nb_points_of_convex(i); ++j) { - if (gmm::sqr(mesh.points_of_convex(i)[j][0]) + - gmm::sqr(mesh.points_of_convex(i)[j][1]) > - gmm::sqr(enr_area_radius)) - pm_convexes.sup(i); - break; - } - } + /* first : selecting the convexes that are completly included in the enrichment area */ + for (dal::bv_visitor i(mesh.convex_index()); !i.finished(); ++i) { + pm_convexes.add(i); + /* For each element, we test all of its nodes. + If all the nodes are inside the enrichment area, + then the element is completly inside the area too */ + for (unsigned j=0; j < mesh.nb_points_of_convex(i); ++j) { + if (gmm::sqr(mesh.points_of_convex(i)[j][0]) + + gmm::sqr(mesh.points_of_convex(i)[j][1]) > + gmm::sqr(enr_area_radius)) + pm_convexes.sup(i); + break; + } + } for (dal::bv_visitor cv(mf_sing_u.convex_index()); !cv.finished(); ++cv) { - if (!pm_convexes.is_in(cv)) - mf_sing_u.set_finite_element(cv, 0); + if (!pm_convexes.is_in(cv)) + mf_sing_u.set_finite_element(cv, 0); } cout << "mf_sing_u: convex_index() = " << mf_sing_u.convex_index().card() << " convexes\n"; @@ -758,100 +761,101 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { mf_u_sum.set_smart_global_dof_linking(true); mf_u_sum.set_mesh_fems(mf_pre_u, mf_sing_u); - cout << "mf_u_sum.nb_dof = " << mf_u_sum.nb_dof() << "\n"; cout << "mfls_u.convex_index = " << mfls_u.convex_index() << "\nmf_sing_u: " << mf_sing_u.convex_index() << "\n"; - } break ; + } break; case 2 : // standard XFEM on a fixed zone { dal::bit_vector enriched_dofs; plain_vector X(mf_partition_of_unity.nb_dof()); plain_vector Y(mf_partition_of_unity.nb_dof()); getfem::interpolation(ls.get_mesh_fem(), mf_partition_of_unity, - ls.values(1), X); + ls.values(1), X); getfem::interpolation(ls.get_mesh_fem(), mf_partition_of_unity, - ls.values(0), Y); + ls.values(0), Y); for (size_type j = 0; j < mf_partition_of_unity.nb_dof(); ++j) { - if (gmm::sqr(X[j]) + gmm::sqr(Y[j]) <= gmm::sqr(enr_area_radius)) - enriched_dofs.add(j); - } + if (gmm::sqr(X[j]) + gmm::sqr(Y[j]) <= gmm::sqr(enr_area_radius)) + enriched_dofs.add(j); + } //cout << "enriched_dofs: " << enriched_dofs << "\n"; if (enriched_dofs.card() < 3) - GMM_WARNING0("There is " << enriched_dofs.card() << - " enriched dofs for the crack tip"); + GMM_WARNING0("There is " << enriched_dofs.card() << + " enriched dofs for the crack tip"); mf_u_product.set_enrichment(enriched_dofs); mf_u_sum.set_mesh_fems(mf_u_product, mfls_u); - cout << "enrichment done \n" ;} - break ; - case 4 :{ - if(cutoff.fun_num == getfem::cutoff_xy_function::EXPONENTIAL_CUTOFF) - cout<<"Using exponential Cutoff..."<<endl; + cout << "enrichment done \n"; + } break; + case 4 : + { + if (cutoff.fun_num == getfem::cutoff_xy_function::EXPONENTIAL_CUTOFF) + cout<<"Using exponential Cutoff..."<<endl; else - cout<<"Using Polynomial Cutoff..."<<endl; -// dal::bit_vector enriched_dofs; -// enriched_dofs.clear() ; -// cout << "mf_cutoff.nb_dof() = " << mf_cutoff.nb_dof() << "\n" ; + cout<<"Using Polynomial Cutoff..."<<endl; +// dal::bit_vector enriched_dofs; +// enriched_dofs.clear(); +// cout << "mf_cutoff.nb_dof() = " << mf_cutoff.nb_dof() << "\n"; // -// for (size_type j = 0; j < mf_cutoff.nb_dof(); ++j) { -// enriched_dofs.add(j); -// } -// for (dal::bv_visitor j(enriched_dofs) ; !j.finished() ; ++j){ -// cout << j << " ; " ; -// } -// cout << "\n" ; -// cout << "mf_prod_cutoff.nb_dof() = " << mf_prod_cutoff.nb_dof() << "\n" ; -// cout << "mf_sing_u.nb_dof() = " << mf_sing_u.nb_dof() << "\n" ; -// mf_prod_cutoff.set_enrichment(enriched_dofs) ; +// for (size_type j = 0; j < mf_cutoff.nb_dof(); ++j) { +// enriched_dofs.add(j); +// } +// for (dal::bv_visitor j(enriched_dofs); !j.finished(); ++j) { +// cout << j << "; "; +// } +// cout << "\n"; +// cout << "mf_prod_cutoff.nb_dof() = " << mf_prod_cutoff.nb_dof() << "\n"; +// cout << "mf_sing_u.nb_dof() = " << mf_sing_u.nb_dof() << "\n"; +// mf_prod_cutoff.set_enrichment(enriched_dofs); // -// cout << "mf_prod_cutoff.nb_dof() = " << mf_prod_cutoff.nb_dof() << "\n" ; +// cout << "mf_prod_cutoff.nb_dof() = " << mf_prod_cutoff.nb_dof() << "\n"; //mf_u_sum.set_mesh_fems(mf_prod_cutoff, mfls_u); mf_u_sum.set_mesh_fems(mf_sing_u, mfls_u); } break; case 3 : // Integral matching (mortar) { - cout << "\nIntegral Matching (Mortar)\n" ; - dal::bit_vector cvlist_in_area, cvlist_out_area; - bool in_area = true; - for (dal::bv_visitor cv(mesh.convex_index()); - !cv.finished(); ++cv) { - in_area = true; - /* For each element, we test all of its nodes. - If all the nodes are inside the enrichment area, - then the element is completly inside the area too */ - for (unsigned j=0; j < mesh.nb_points_of_convex(cv); ++j) { -// if (gmm::sqr(mesh.points_of_convex(cv)[j][0] ) + -// gmm::sqr(mesh.points_of_convex(cv)[j][1] ) > -// gmm::sqr(enr_area_radius)) { + cout << "\nIntegral Matching (Mortar)\n"; + dal::bit_vector cvlist_in_area, cvlist_out_area; + bool in_area = true; + for (dal::bv_visitor cv(mesh.convex_index()); + !cv.finished(); ++cv) { + in_area = true; + /* For each element, we test all of its nodes. + If all the nodes are inside the enrichment area, + then the element is completly inside the area too */ + for (unsigned j=0; j < mesh.nb_points_of_convex(cv); ++j) { +// if (gmm::sqr(mesh.points_of_convex(cv)[j][0] ) + +// gmm::sqr(mesh.points_of_convex(cv)[j][1] ) > +// gmm::sqr(enr_area_radius)) { if ( (gmm::abs(mesh.points_of_convex(cv)[j][0] ) > enr_area_radius) || (gmm::abs(mesh.points_of_convex(cv)[j][1] ) > enr_area_radius)) { - in_area = false; - break; - } - } - - /* "remove" the global function on convexes outside the enrichment - area */ - if (!in_area) { - cvlist_out_area.add(cv); - mf_sing_u.set_finite_element(cv, 0); - mf_u().set_dof_partition(cv, 1); - } else cvlist_in_area.add(cv); + in_area = false; + break; + } + } + + /* "remove" the global function on convexes outside the enrichment + area */ + if (!in_area) { + cvlist_out_area.add(cv); + mf_sing_u.set_finite_element(cv, 0); + mf_u().set_dof_partition(cv, 1); + } else + cvlist_in_area.add(cv); } - /* extract the boundary of the enrichment area, from the - "inside" point-of-view, and from the "outside" - point-of-view */ + /* extract the boundary of the enrichment area, from the + "inside" point-of-view, and from the "outside" + point-of-view */ getfem::mesh_region r_border, r_enr_out; getfem::outer_faces_of_mesh(mesh, r_border); getfem::outer_faces_of_mesh(mesh, cvlist_in_area, - mesh.region(MORTAR_BOUNDARY_IN)); + mesh.region(MORTAR_BOUNDARY_IN)); getfem::outer_faces_of_mesh(mesh, cvlist_out_area, - mesh.region(MORTAR_BOUNDARY_OUT)); + mesh.region(MORTAR_BOUNDARY_OUT)); for (getfem::mr_visitor v(r_border); !v.finished(); ++v) { - mesh.region(MORTAR_BOUNDARY_OUT).sup(v.cv(), v.f()); + mesh.region(MORTAR_BOUNDARY_OUT).sup(v.cv(), v.f()); } if (PARAM.int_value("MORTAR_WITHOUT_SINGUL")) mf_u_sum.set_mesh_fems(mfls_u); @@ -866,18 +870,18 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { // an optional treatment : creating a representation of the enrichment area getfem::mesh_fem mf_enrich(mesh); - for (dal::bv_visitor i(mesh.convex_index()) ; !i.finished() ; ++i){ - getfem::pfem pf_mef = getfem::classical_fem(mesh.trans_of_convex(i), 1 ); - mf_enrich.set_finite_element(i, pf_mef) ; + for (dal::bv_visitor i(mesh.convex_index()); !i.finished(); ++i) { + getfem::pfem pf_mef = getfem::classical_fem(mesh.trans_of_convex(i), 1); + mf_enrich.set_finite_element(i, pf_mef); } - std::vector<scalar_type> UU(mf_enrich.nb_dof()) ; - std::fill(UU.begin(), UU.end() ,0.) ; - cout << "exporting the enrichment zone: \n" ; + std::vector<scalar_type> UU(mf_enrich.nb_dof()); + std::fill(UU.begin(), UU.end() ,0.); + cout << "exporting the enrichment zone: \n"; GMM_ASSERT1(!mf_enrich.is_reduced(), "To be adapted"); - for (dal::bv_visitor i(cvlist_in_area) ; !i.finished() ; ++i){ - for (unsigned int j = 0 ; - j < mf_enrich.ind_basic_dof_of_element(i).size() ; ++j ) - UU[mf_enrich.ind_basic_dof_of_element(i)[j]] = 1. ; + for (dal::bv_visitor i(cvlist_in_area); !i.finished(); ++i) { + for (unsigned int j = 0; + j < mf_enrich.ind_basic_dof_of_element(i).size(); ++j) + UU[mf_enrich.ind_basic_dof_of_element(i)[j]] = 1.; } cout << "exporting enrichment to " << "enrichment_zone.vtk" << "..\n"; @@ -885,37 +889,37 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { exp.exporting(mf_enrich); exp.write_point_data(mf_enrich, UU, "enrichment"); cout << "export done, you can view the data file with (for example)\n" - "mayavi -d enrichment_zone.vtk -f " - "WarpScalar -m BandedSurfaceMap -m Outline\n"; - -// // Another optional treatment : -// // Searching the elements that are both crossed by the crack -// // and with one of their faces which constitutes a part of the -// // boundary between the enriched zone and the rest of the domain. -// getfem::mesh_region &boundary = mesh.region(MORTAR_BOUNDARY_IN); -// unsigned int cpt = 0 ; -// for (dal::bv_visitor i(cvlist_in_area); !i.finished(); ++i) { -// if (mls.is_convex_cut(i)){ -// // Among the faces of the convex, we search if some are -// // part of the boundary -// cpt = 0 ; -// for (unsigned j=0; j < mesh.structure_of_convex(i) ->nb_faces(); ++j) { -// if (boundary.is_in(i,j)) -// cpt += 1; -// } -// if (cpt) { -// cout << "\n The convex number " << i << " is crossed by the crack :\n" ; -// cout << " it has : " << cpt << " face(s) among the boundary.\n \n " ; -// } -// } -// } + "mayavi -d enrichment_zone.vtk -f " + "WarpScalar -m BandedSurfaceMap -m Outline\n"; + +// // Another optional treatment : +// // Searching the elements that are both crossed by the crack +// // and with one of their faces which constitutes a part of the +// // boundary between the enriched zone and the rest of the domain. +// getfem::mesh_region &boundary = mesh.region(MORTAR_BOUNDARY_IN); +// unsigned int cpt = 0; +// for (dal::bv_visitor i(cvlist_in_area); !i.finished(); ++i) { +// if (mls.is_convex_cut(i)) { +// // Among the faces of the convex, we search if some are +// // part of the boundary +// cpt = 0; +// for (unsigned j=0; j < mesh.structure_of_convex(i) ->nb_faces(); ++j) { +// if (boundary.is_in(i,j)) +// cpt += 1; +// } +// if (cpt) { +// cout << "\n The convex number " << i << " is crossed by the crack :\n"; +// cout << " it has : " << cpt << " face(s) among the boundary.\n \n "; +// } +// } +// } } // end of "enrichment_option = 3" - break ; + break; default : - GMM_ASSERT1(false, "Enrichment_option parameter is undefined"); - break ; - } // end of switch + GMM_ASSERT1(false, "Enrichment_option parameter is undefined"); + break; + } // end of switch mesh.write_to_file("toto.mesh"); @@ -953,8 +957,8 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { //cout << "validate mf_u():\n"; validate_fem_derivatives(mf_u()); cout << "Number of dof for u --> : " << mf_u().nb_dof() << endl; - scalar_type pressure ; - pressure = PARAM.real_value("PRESSURE") ; + scalar_type pressure; + pressure = PARAM.real_value("PRESSURE"); // Bilaplacian brick. getfem::model model; @@ -973,11 +977,11 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { // Defining the volumic source term. size_type nb_dof_rhs = mf_rhs.nb_dof(); plain_vector F(nb_dof_rhs); - gmm::clear(F) ; - if (sol_ref == 2){ - getfem::interpolation_function(mf_rhs, F, sol_F); - gmm::scale(F, pressure) ; - } + gmm::clear(F); + if (sol_ref == 2) { + getfem::interpolation_function(mf_rhs, F, sol_F); + gmm::scale(F, pressure); + } //Volumic source term brick. model.add_initialized_fem_data("VolumicData", mf_rhs, F); getfem::add_source_term_brick(model, mim, "u", "VolumicData"); @@ -999,18 +1003,15 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { (model, mim, "u", mf_mult, SIMPLE_SUPPORT_BOUNDARY_NUM, "DData"); - - - sparse_matrix H(1, mf_u().nb_dof()); if (enrichment_option == 3 ) { - /* add a constraint brick for the mortar junction between + /* add a constraint brick for the mortar junction between the enriched area and the rest of the mesh */ // calcul des matrices de contraintes - plain_vector R(1) ; -// sparse_matrix H(1, mf_u().nb_dof()); - this->set_matrix_mortar(H) ; + plain_vector R(1); +// sparse_matrix H(1, mf_u().nb_dof()); + this->set_matrix_mortar(H); /* because of the discontinuous partition of mf_u(), some levelset enriched functions do not contribute any more to the @@ -1020,22 +1021,21 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { sparse_matrix M2(mf_u().nb_dof(), mf_u().nb_dof()); getfem::asm_mass_matrix(M2, mim, mf_u(), mf_u()); //gmm::HarwellBoeing_IO::write("M2.hb", M2); - cout << "PARAM.real_value(\"SEUIL\") : " << PARAM.real_value("SEUIL") << "\n" ; + cout << "PARAM.real_value(\"SEUIL\") : " << PARAM.real_value("SEUIL") << "\n"; for (size_type d = 0; d < mf_u().nb_dof(); ++d) { - // if (M2(d,d) < 1e-7) cout << " weak mf_u() dof " << d << " @ " << - // mf_u().point_of_dof(d) << " M2(d,d) = " << M2(d,d) << "\n"; + // if (M2(d,d) < 1e-7) cout << " weak mf_u() dof " << d << " @ " << + // mf_u().point_of_dof(d) << " M2(d,d) = " << M2(d,d) << "\n"; if (M2(d,d) < PARAM.real_value("SEUIL")) { - cout << "removed : " << mf_u().point_of_basic_dof(d) << "\n"; - size_type n = gmm::mat_nrows(H); - gmm::resize(H, n+1, gmm::mat_ncols(H)); - H(n, d) = 1; + cout << "removed : " << mf_u().point_of_basic_dof(d) << "\n"; + size_type n = gmm::mat_nrows(H); + gmm::resize(H, n+1, gmm::mat_ncols(H)); + H(n, d) = 1; } } gmm::resize(R, gmm::mat_nrows(H)); model.add_fixed_size_variable("mult_mo", gmm::mat_nrows(H)); getfem::add_constraint_with_multipliers(model, "u", "mult_mo", H, R); gmm::Harwell_Boeing_save("H.hb", H); - } if ( PARAM.real_value("SEUIL_FINAL")!=0 ) { @@ -1048,19 +1048,19 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { getfem::asm_stiffness_matrix_for_bilaplacian(M2, mim, mf_u(), mf_rhs, RR); - //cout << "stiffness_matrix_for_bilaplacian : " << M2 << "\n" ; - cout << "termes diagonaux, de la matrice de rigidite, inferieurs a 1e-10 : " ; + //cout << "stiffness_matrix_for_bilaplacian : " << M2 << "\n"; + cout << "termes diagonaux, de la matrice de rigidite, inferieurs a 1e-10 : "; for (size_type d = 0; d < mf_u().nb_dof(); ++d) { - if (M2(d,d) < 1e-10) cout << M2(d,d) << " ; " ; + if (M2(d,d) < 1e-10) cout << M2(d,d) << " ; "; } - cout << "\n" ; - cout << "SEUIL_FINAL = " << PARAM.real_value("SEUIL_FINAL") << "\n" ; + cout << "\n"; + cout << "SEUIL_FINAL = " << PARAM.real_value("SEUIL_FINAL") << "\n"; for (size_type d = 0; d < mf_u().nb_dof(); ++d) { if (M2(d,d) < PARAM.real_value("SEUIL_FINAL")) { - cout << "OULALA " << d << " @ " << mf_u().point_of_basic_dof(d) << " : " << M2(d,d) << "\n"; + cout << "OULALA " << d << " @ " << mf_u().point_of_basic_dof(d) << " : " << M2(d,d) << "\n"; size_type n = gmm::mat_nrows(H); - gmm::resize(H1, n+1, gmm::mat_ncols(H)); - H1(n, d) = 1; + gmm::resize(H1, n+1, gmm::mat_ncols(H)); + H1(n, d) = 1; } } base_vector R(gmm::mat_nrows(H1)); @@ -1069,149 +1069,145 @@ bool bilaplacian_crack_problem::solve(plain_vector &U) { gmm::Harwell_Boeing_save("M2.hb", M2); } - cout << "Total number of variables : " << model.nb_dof() << endl; gmm::iteration iter(residual, 1, 40000); getfem::standard_solve(model, iter); - // Solution extraction gmm::resize(U, mf_u().nb_dof()); gmm::copy(model.real_variable("u"), U); /****************************/ if (PARAM.int_value("FIC_ORTHO") ) { + sparse_matrix A = model.real_tangent_matrix(); + plain_vector b = model.real_rhs(); + gmm::scale(b, -1.); + plain_vector X(b); + scalar_type condest; + SuperLU_solve(A, X, gmm::scaled(b, scalar_type(-1)), condest, 1); + cout << "cond super LU = " << 1./condest << "\n"; + cout << "X = " << gmm::sub_vector(X, gmm::sub_interval(0, 10)) << "\n"; + cout << "U = " << gmm::sub_vector(U, gmm::sub_interval(0, 10)) << "\n"; + + unsigned q = mf_u().get_qdim(); + + base_small_vector tab_fic(4); + std::vector<size_type> ind_sing(2); + unsigned cpt = 0; + if (PARAM.int_value("ENRICHMENT_OPTION") == 3) { + // affichage des coeffs devant les singularites, avec le raccord integral + GMM_ASSERT1(!mf_u().is_reduced(), "To be adapted"); + + for (unsigned d=0; d < mf_u().nb_dof(); d += q) { + size_type cv = mf_u().first_convex_of_basic_dof(d); + getfem::pfem pf = mf_u().fem_of_element(cv); + unsigned ld = unsigned(-1); + for (unsigned dd = 0; dd < mf_u().nb_basic_dof_of_element(cv); dd += q) { + if (mf_u().ind_basic_dof_of_element(cv)[dd] == d) { + ld = dd/q; break; + } + } + if (ld == unsigned(-1)) { + cout << "DOF " << d << "NOT FOUND in " << cv << " BUG BUG\n"; + } + else { + if (is_global_dof_type_bis(pf->dof_types().at(ld))) { + cout << "coeff:" << U[d] << "\n"; + cout << "dof index:" << d << "\n"; + tab_fic[cpt] = U[d]; + ind_sing[cpt] = d; + cpt +=1; + } + } + } + } - sparse_matrix A = model.real_tangent_matrix() ; - plain_vector b = model.real_rhs() ; - gmm::scale(b, -1.) ; - plain_vector X(b) ; - scalar_type condest ; - SuperLU_solve(A, X, gmm::scaled(b, scalar_type(-1)), condest, 1) ; - cout << "cond super LU = " << 1./condest << "\n" ; - cout << "X = " << gmm::sub_vector(X, gmm::sub_interval(0, 10)) << "\n" ; - cout << "U = " << gmm::sub_vector(U, gmm::sub_interval(0, 10)) << "\n" ; - - - unsigned q = mf_u().get_qdim(); - - base_small_vector tab_fic(4); - std::vector<size_type> ind_sing(2) ; - unsigned cpt = 0; - if (PARAM.int_value("ENRICHMENT_OPTION") == 3){ - // affichage des coeffs devant les singularites, avec le raccord integral - GMM_ASSERT1(!mf_u().is_reduced(), "To be adapted"); - - for (unsigned d=0; d < mf_u().nb_dof(); d += q) { - size_type cv = mf_u().first_convex_of_basic_dof(d) ; - getfem::pfem pf = mf_u().fem_of_element(cv); - unsigned ld = unsigned(-1); - for (unsigned dd = 0; dd < mf_u().nb_basic_dof_of_element(cv); dd += q) { - if (mf_u().ind_basic_dof_of_element(cv)[dd] == d) { - ld = dd/q; break; - } - } - if (ld == unsigned(-1)) { - cout << "DOF " << d << "NOT FOUND in " << cv << " BUG BUG\n"; - } - else { - if ( is_global_dof_type_bis(pf->dof_types().at(ld)) ){ - cout << "coeff:" << U[d] << "\n" ; - cout << "dof index:" << d << "\n" ; - tab_fic[cpt] = U[d] ; - ind_sing[cpt] = d ; - cpt +=1 ; - } - } - } - } - - - plain_vector b1(gmm::mat_nrows(A)), b2(b1), X1(b1), X2(b1) ; - - scalar_type as1, as2, as1s2, bs1, bs2 ; - as1 = A(ind_sing[0], ind_sing[0]) ; - scalar_type max = 0. ; - size_type imax = 0 ; - for (unsigned i = 0 ; i < gmm::mat_nrows(A) ; i++){ - if (gmm::abs(A(ind_sing[0],i)) > max && i!= ind_sing[0] && i!= ind_sing[1]){ - max = gmm::abs(A(ind_sing[0],i)) ; - imax = i ; - } - } - cout << "imax = " << imax << endl ; - cout << "max = " << max << endl ; - if (imax < mf_u().nb_dof()) - cout << "position de imax = [" << mf_u().point_of_basic_dof(imax)[0] << " ; " << mf_u().point_of_basic_dof(imax)[1] << "\n" ; - as2 = A(ind_sing[1], ind_sing[1]) ; - as1s2 = A(ind_sing[0], ind_sing[1]) ; - bs1 = b[ind_sing[0]] ; - bs2 = b[ind_sing[1]] ; - gmm::copy(gmm::mat_col(A, ind_sing[0]), b1) ; - gmm::copy(gmm::mat_col(A, ind_sing[1]), b2) ; - //cout << "gmm::mat_col(A, ind_sing[0]) = " << b1 << "\n" ; - //cout << "gmm::mat_col(A, ind_sing[1]) = " << b2 << "\n" ; - - for (int i=0 ; i < 2 ; i++){ - for (unsigned j=0 ; j < gmm::mat_nrows(A) ; j++){ - A(ind_sing[i],j) = 0. ; - A(j,ind_sing[i]) = 0. ; - } - A(ind_sing[i], ind_sing[i]) = 1. ; - b[ind_sing[i]] = 0. ; - b1[ind_sing[i]] = 0. ; - b2[ind_sing[i]] = 0. ; - } - - SuperLU_solve(A, X1, b1, condest, 1) ; - cout << "solving for s1 OK, cond = " << 1./condest << "\n" ; - SuperLU_solve(A, X2, b2, condest, 1) ; - cout << "solving for s2 OK, cond = " << 1./condest << "\n" ; - cout << "X1[ind_sing[0]] = " << X1[ind_sing[0]] << "\n" ; - cout << "X1 = " << gmm::sub_vector(X1, gmm::sub_interval(0, 10)) << "\n" ; - scalar_type max1 = 0., max2 = 0. ; - size_type imax1 = 0 , imax2 = 0 ; - for (unsigned i = 0 ; i < X1.size() ; i++){ - if (gmm::abs(X1[i]) > max1){ - max1 = gmm::abs(X1[i]) ; - imax1 = i ; - } - if (gmm::abs(X2[i]) > max2){ - max2 = gmm::abs(X2[i]) ; - imax2 = i ; - } - } - cout << "imax1 = " << imax1 << endl ; - cout << "max1 = " << max1 << endl ; - if (imax1 < mf_u().nb_dof()) - cout << "position de imax1 = [" << mf_u().point_of_basic_dof(imax1)[0] << " ; " << mf_u().point_of_basic_dof(imax1)[1] << "\n" ; - cout << "imax2 = " << imax2 << endl ; - cout << "max2 = " << max2 << endl ; - if (imax2 < mf_u().nb_dof()) - cout << "position de imax2 = [" << mf_u().point_of_basic_dof(imax2)[0] << " ; " << mf_u().point_of_basic_dof(imax2)[1] << "\n" ; - //cout << "X1 = " << gmm::sub_vector(X1, gmm::sub_interval(0, 100)) << "\n" ; - //cout << "X2 = " << gmm::sub_vector(X2, gmm::sub_interval(0, 100)) << "\n" ; - base_matrix M(2,2) ; - plain_vector AX1(gmm::mat_nrows(A)), AX2(AX1) ; - gmm::mult(A, X1, AX1) ; - gmm::mult(A, X2, AX2) ; - M(0,0) = as1 - 2. * gmm::vect_sp(b1, X1) + gmm::vect_sp(X1, AX1) ; - M(1,1) = as2 - 2. * gmm::vect_sp(b2, X2) + gmm::vect_sp(X2, AX2) ; - M(0,1) = as1s2 - gmm::vect_sp(b1, X2) - gmm::vect_sp(b2, X1) + gmm::vect_sp(X1, AX2) ; - M(1,0) = M(0,1) ; - plain_vector Z(2), FIC_ORTHO(2) ; - Z[0] = bs1 - gmm::vect_sp(X1, b) ; - Z[1] = bs2 - gmm::vect_sp(X2, b) ; - gmm::lu_solve(M, FIC_ORTHO, Z) ; - - cout << "[as1 as2 as1s2] = " << as1 << " ; " << as2 << " ; " << as1s2 << "\n" ; - cout << "[bs1 bs2] = " << bs1 << " ; " << bs2 << "\n" ; - cout << "M = " << M << "\n" ; - cout << "Z = " << Z << "\n" ; - - cout << "FIC1 ORTHO = " << FIC_ORTHO[0] << "\n" ; - cout << "FIC2 ORTHO = " << FIC_ORTHO[1] << "\n" ; - cout << "-----------------------\n" ; + + plain_vector b1(gmm::mat_nrows(A)), b2(b1), X1(b1), X2(b1); + + scalar_type as1, as2, as1s2, bs1, bs2; + as1 = A(ind_sing[0], ind_sing[0]); + scalar_type max = 0.; + size_type imax = 0; + for (unsigned i = 0; i < gmm::mat_nrows(A); i++) { + if (gmm::abs(A(ind_sing[0],i)) > max && i!= ind_sing[0] && i!= ind_sing[1]) { + max = gmm::abs(A(ind_sing[0],i)); + imax = i; + } + } + cout << "imax = " << imax << endl; + cout << "max = " << max << endl; + if (imax < mf_u().nb_dof()) + cout << "position de imax = [" << mf_u().point_of_basic_dof(imax)[0] << "; " << mf_u().point_of_basic_dof(imax)[1] << "\n"; + as2 = A(ind_sing[1], ind_sing[1]); + as1s2 = A(ind_sing[0], ind_sing[1]); + bs1 = b[ind_sing[0]]; + bs2 = b[ind_sing[1]]; + gmm::copy(gmm::mat_col(A, ind_sing[0]), b1); + gmm::copy(gmm::mat_col(A, ind_sing[1]), b2); + //cout << "gmm::mat_col(A, ind_sing[0]) = " << b1 << "\n"; + //cout << "gmm::mat_col(A, ind_sing[1]) = " << b2 << "\n"; + + for (int i=0; i < 2; i++) { + for (unsigned j=0; j < gmm::mat_nrows(A); j++) { + A(ind_sing[i],j) = 0.; + A(j,ind_sing[i]) = 0.; + } + A(ind_sing[i], ind_sing[i]) = 1.; + b[ind_sing[i]] = 0.; + b1[ind_sing[i]] = 0.; + b2[ind_sing[i]] = 0.; + } + + SuperLU_solve(A, X1, b1, condest, 1); + cout << "solving for s1 OK, cond = " << 1./condest << "\n"; + SuperLU_solve(A, X2, b2, condest, 1); + cout << "solving for s2 OK, cond = " << 1./condest << "\n"; + cout << "X1[ind_sing[0]] = " << X1[ind_sing[0]] << "\n"; + cout << "X1 = " << gmm::sub_vector(X1, gmm::sub_interval(0, 10)) << "\n"; + scalar_type max1 = 0., max2 = 0.; + size_type imax1 = 0 , imax2 = 0; + for (unsigned i = 0; i < X1.size(); i++) { + if (gmm::abs(X1[i]) > max1) { + max1 = gmm::abs(X1[i]); + imax1 = i; + } + if (gmm::abs(X2[i]) > max2) { + max2 = gmm::abs(X2[i]); + imax2 = i; + } + } + cout << "imax1 = " << imax1 << endl; + cout << "max1 = " << max1 << endl; + if (imax1 < mf_u().nb_dof()) + cout << "position de imax1 = [" << mf_u().point_of_basic_dof(imax1)[0] << " ; " << mf_u().point_of_basic_dof(imax1)[1] << "\n"; + cout << "imax2 = " << imax2 << endl; + cout << "max2 = " << max2 << endl; + if (imax2 < mf_u().nb_dof()) + cout << "position de imax2 = [" << mf_u().point_of_basic_dof(imax2)[0] << " ; " << mf_u().point_of_basic_dof(imax2)[1] << "\n"; + //cout << "X1 = " << gmm::sub_vector(X1, gmm::sub_interval(0, 100)) << "\n"; + //cout << "X2 = " << gmm::sub_vector(X2, gmm::sub_interval(0, 100)) << "\n"; + base_matrix M(2,2); + plain_vector AX1(gmm::mat_nrows(A)), AX2(AX1); + gmm::mult(A, X1, AX1); + gmm::mult(A, X2, AX2); + M(0,0) = as1 - 2. * gmm::vect_sp(b1, X1) + gmm::vect_sp(X1, AX1); + M(1,1) = as2 - 2. * gmm::vect_sp(b2, X2) + gmm::vect_sp(X2, AX2); + M(0,1) = as1s2 - gmm::vect_sp(b1, X2) - gmm::vect_sp(b2, X1) + gmm::vect_sp(X1, AX2); + M(1,0) = M(0,1); + plain_vector Z(2), FIC_ORTHO(2); + Z[0] = bs1 - gmm::vect_sp(X1, b); + Z[1] = bs2 - gmm::vect_sp(X2, b); + gmm::lu_solve(M, FIC_ORTHO, Z); + + cout << "[as1 as2 as1s2] = " << as1 << " ; " << as2 << " ; " << as1s2 << "\n"; + cout << "[bs1 bs2] = " << bs1 << " ; " << bs2 << "\n"; + cout << "M = " << M << "\n"; + cout << "Z = " << Z << "\n"; + + cout << "FIC1 ORTHO = " << FIC_ORTHO[0] << "\n"; + cout << "FIC2 ORTHO = " << FIC_ORTHO[1] << "\n"; + cout << "-----------------------\n"; } /************************************************/ return true; diff --git a/interface/src/gf_linsolve.cc b/interface/src/gf_linsolve.cc index 8b42282f..c4c6f063 100644 --- a/interface/src/gf_linsolve.cc +++ b/interface/src/gf_linsolve.cc @@ -97,7 +97,7 @@ superlu_solver(gsparse &gsp, out.pop().from_scalar(rcond ? 1./rcond : 0.); } -#if defined(GMM_USES_MUMPS) || defined(HAVE_DMUMPS_C_H) +#if defined(GMM_USES_MUMPS) template <typename T> static void mumps_solver(gsparse &gsp, getfemint::mexargs_in& in, getfemint::mexargs_out& out, T) { @@ -204,7 +204,7 @@ void gf_linsolve(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { else superlu_solver(gsp, in, out, scalar_type()); ); -#if defined(GMM_USES_MUMPS) || defined(HAVE_DMUMPS_C_H) +#if defined(GMM_USES_MUMPS) /*@FUNC @CELL{U, cond} = ('mumps', @tsp M, @vec b) Solve `M.U = b` using the MUMPS solver.@*/ sub_command diff --git a/src/gmm/gmm_solver_Schwarz_additive.h b/src/gmm/gmm_solver_Schwarz_additive.h index a37bbef4..2d6e9fdd 100644 --- a/src/gmm/gmm_solver_Schwarz_additive.h +++ b/src/gmm/gmm_solver_Schwarz_additive.h @@ -333,8 +333,8 @@ namespace gmm { // { // MPI_Sendrecv(&(qbis[0]), gmm::vect_size(q), MPI_DOUBLE, next, tag1, -// &(qter[0]), gmm::vect_size(q),MPI_DOUBLE,previous,tag1, -// MPI_COMM_WORLD,&status); +// &(qter[0]), gmm::vect_size(q),MPI_DOUBLE,previous,tag1, +// MPI_COMM_WORLD,&status); // gmm::copy(qter, qbis); // add(qbis,q,q); // } diff --git a/tests/schwarz_additive.cc b/tests/schwarz_additive.cc index 68040126..47f0d481 100644 --- a/tests/schwarz_additive.cc +++ b/tests/schwarz_additive.cc @@ -26,7 +26,7 @@ /* */ /**************************************************************************/ -#define GMM_USES_SUPERLU +#define GETFEM_USES_SUPERLU #include "getfem/getfem_assembling.h" #include "getfem/getfem_regular_meshes.h" @@ -77,7 +77,9 @@ struct pb_data { int solve_cg(); int solve_cg2(); +#if defined(GETFEM_USES_SUPERLU) int solve_superlu(); +#endif int solve_schwarz(int); int solve() { @@ -85,7 +87,9 @@ struct pb_data { switch (solver) { case 0 : return solve_cg(); case 1 : return solve_cg2(); +#if defined(GETFEM_USES_SUPERLU) case 2 : return solve_superlu(); +#endif default : return solve_schwarz(solver); } return 0; @@ -224,11 +228,13 @@ int pb_data::solve_cg() { return int(iter.get_iteration()); } +#if defined(GETFEM_USES_SUPERLU) int pb_data::solve_superlu() { double rcond; SuperLU_solve(RM, U, F, rcond); return 1; } +#endif int pb_data::solve_cg2() { gmm::iteration iter(residual, 1, 1000000); @@ -269,10 +275,12 @@ int pb_data::solve_schwarz(int version) { gmm::ilu_precond<general_sparse_matrix>(), vB, iter, gmm::using_gmres(), gmm::using_gmres()); break; +#if defined(GETFEM_USES_SUPERLU) case 5 : gmm::additive_schwarz(RM, U, F, gmm::ilu_precond<general_sparse_matrix>(), vB, iter, gmm::using_superlu(), gmm::using_cg()); break; +#endif } return 0; }