[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] [getfem-commits] branch master updated: Whitespace fixe
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] [getfem-commits] branch master updated: Whitespace fixes and other minor changes |
|
Date: |
Tue, 26 Dec 2023 14:40:22 -0500 |
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 <logari81@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;
}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] [getfem-commits] branch master updated: Whitespace fixes and other minor changes,
Konstantinos Poulios <=