GetFEM  5.5
getfem_error_estimate.cc
1 /*===========================================================================
2 
3  Copyright (C) 2013-2026 Yves Renard
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program. If not, see https://www.gnu.org/licenses/.
18 
19 ===========================================================================*/
20 
21 
22 #include <getfem/getfem_error_estimate.h>
23 #include <getfem/getfem_contact_and_friction_common.h>
24 #include <getfem/getfem_mesher.h>
25 
26 namespace getfem {
27 
28 #ifdef EXPERIMENTAL_PURPOSE_ONLY
29 
30 
31  void error_estimate_nitsche(const mesh_im & mim,
32  const mesh_fem &mf_u,
33  const base_vector &U,
34  int GAMMAC,
35  int GAMMAN,
36  scalar_type lambda,
37  scalar_type mu,
38  scalar_type gamma0,
39  scalar_type f_coeff,
40  scalar_type vertical_force,
41  base_vector &ERR) {
42 
43 
44 
45 
46  // static double lambda, mu;
47  const mesh &m = mf_u.linked_mesh();
48  size_type N = m.dim();
49  size_type qdim = mf_u.get_qdim();
50  gmm::clear(ERR);
51  std::vector<scalar_type> coeff1, coeff2;
52  base_matrix grad1(qdim, N), grad2(qdim, N), E(N, N), S1(N, N), S2(N, N);
53  base_matrix hess1(qdim, N*N);
54  base_matrix G1, G2;
55  bgeot::geotrans_inv_convex gic;
56  base_node xref2(N);
57  base_small_vector up(N), jump(N), U1(N), sig(N), sigt(N);
58  //scalar_type young_modulus = mu*(3*lambda + 2*mu)/(lambda+mu);
59  scalar_type Pr, scal, sign, Un ;
60  scalar_type eta1 = 0, eta2 = 0, eta3 = 0, eta4 = 0;
61  scalar_type force_coeff= f_coeff; // inutile pour l'instant
62  force_coeff =0;
63 
64  //vertical force
65 
66  base_small_vector F(N);
67  for (unsigned ii=0; ii < N-1; ++ii) F[ii]=0;
68  F[N-1]=-vertical_force;
69 
70  GMM_ASSERT1(!mf_u.is_reduced(), "To be adapted");
71 
72  for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
73 
74  bgeot::pgeometric_trans pgt1 = m.trans_of_convex(cv);
75  getfem::papprox_integration pai1 =
76  get_approx_im_or_fail(mim.int_method_of_element(cv));
77  getfem::pfem pf1 = mf_u.fem_of_element(cv);
78  scalar_type radius = m.convex_radius_estimate(cv);
79  m.points_of_convex(cv, G1);
80  coeff1.resize(mf_u.nb_basic_dof_of_element(cv));
81  gmm::copy(gmm::sub_vector(U, gmm::sub_index(mf_u.ind_basic_dof_of_element(cv))), coeff1);
82  getfem::fem_interpolation_context ctx1(pgt1, pf1, base_node(N), G1, cv);
83 
84  // Residual on the element
85 
86  for (unsigned ii=0; ii < pai1->nb_points_on_convex(); ++ii) {
87  base_small_vector res(N);
88  ctx1.set_xref(pai1->point(ii));
89  pf1->interpolation_hess(ctx1, coeff1, hess1, dim_type(qdim));
90 
91  for (size_type i = 0; i < N; ++i)
92  for (size_type j = 0; j < N; ++j)
93  res[i] += (lambda + mu) * hess1(j, i*N+j) + mu * hess1(i, j*N+j)+F[i];
94 
95  ERR[cv] += radius*radius*ctx1.J()*pai1->coeff(ii)*gmm::vect_norm2_sqr(res); //norme carrée
96  eta1 += (radius*radius*ctx1.J()*pai1->coeff(ii)*gmm::vect_norm2_sqr(res));
97  }
98  //if (ERR[cv] > 100)
99  //cout << "Erreur en résidu sur element " << cv << " : " << ERR[cv] << endl;
100 
101 
102  // jump of the stress between the element ant its neighbors.
103 
104  for (short_type f1=0; f1 < m.structure_of_convex(cv)->nb_faces(); ++f1) {
105 
106  size_type cvn = m.neighbor_of_convex(cv, f1);
107  if (cvn == size_type(-1)) continue;
108 
109  bgeot::pgeometric_trans pgt2 = m.trans_of_convex(cvn);
110  getfem::pfem pf2 = mf_u.fem_of_element(cvn);
111  m.points_of_convex(cvn, G2);
112  coeff2.resize(mf_u.nb_basic_dof_of_element(cvn));
113  gmm::copy(gmm::sub_vector(U, gmm::sub_index(mf_u.ind_basic_dof_of_element(cvn))), coeff2);
114  getfem::fem_interpolation_context ctx2(pgt2, pf2, base_node(N), G2, cvn);
115  gic.init(m.points_of_convex(cvn), pgt2);
116  for (unsigned ii=0; ii < pai1->nb_points_on_face(f1); ++ii) {
117 
118  ctx1.set_xref(pai1->point_on_face(f1, ii));
119  gmm::mult(ctx1.B(), pgt1->normals()[f1], up);
120  scalar_type norm = gmm::vect_norm2(up);
121  up /= norm;
122  scalar_type coefficient = pai1->coeff_on_face(f1, ii) * ctx1.J() * norm;
123  pf1->interpolation_grad(ctx1, coeff1, grad1, dim_type(qdim));
124  gmm::copy(grad1, E); gmm::add(gmm::transposed(grad1), E);
125  gmm::scale(E, 0.5);
126  scalar_type trace = gmm::mat_trace(E);
127  gmm::copy(gmm::identity_matrix(), S1);
128  gmm::scale(S1, lambda * trace);
129  gmm::add(gmm::scaled(E, 2*mu), S1);
130  bool converged;
131  gic.invert(ctx1.xreal(), xref2, converged);
132  GMM_ASSERT1(converged, "geometric transformation not well inverted ... !");
133  ctx2.set_xref(xref2);
134  pf2->interpolation_grad(ctx2, coeff2, grad2, dim_type(qdim));
135  gmm::copy(grad2, E); gmm::add(gmm::transposed(grad2), E);
136  gmm::scale(E, 0.5);
137  trace = gmm::mat_trace(E);
138  gmm::copy(gmm::identity_matrix(), S2);
139  gmm::scale(S2, lambda * trace);
140  gmm::add(gmm::scaled(E, 2*mu), S2);
141  gmm::mult(S1, up, jump);
142  gmm::mult_add(S2, gmm::scaled(up, -1.0), jump);
143 
144  ERR[cv] +=radius * coefficient * gmm::vect_norm2_sqr(jump);
145  eta2 += (radius * coefficient * gmm::vect_norm2_sqr(jump));
146 
147 
148  }
149 
150 
151  }
152 
153 
154  };
155 
156 
157  {
158 
159  int bnum = GAMMAN;
160 
161  getfem::mesh_region region = m.region(bnum);
162  for (getfem::mr_visitor v(region,m); !v.finished(); ++v) {
163 
164  bgeot::pgeometric_trans pgt1 = m.trans_of_convex(v.cv());
165  getfem::papprox_integration pai1 =
166  get_approx_im_or_fail(mim.int_method_of_element(v.cv()));
167  getfem::pfem pf1 = mf_u.fem_of_element(v.cv());
168  scalar_type radius = m.convex_radius_estimate(v.cv());
169 
170  m.points_of_convex(v.cv(), G1);
171 
172  coeff1.resize(mf_u.nb_basic_dof_of_element(v.cv()));
173  gmm::copy(gmm::sub_vector(U, gmm::sub_index(mf_u.ind_basic_dof_of_element(v.cv()))), coeff1);
174 
175  getfem::fem_interpolation_context ctx1(pgt1, pf1, base_node(N), G1, v.cv());
176 
177  short_type f = v.f();
178 
179  for (unsigned ii=0; ii < pai1->nb_points_on_face(f); ++ii) {
180 
181  ctx1.set_xref(pai1->point_on_face(f, ii));
182  gmm::mult(ctx1.B(), pgt1->normals()[f], up);
183  scalar_type norm = gmm::vect_norm2(up);
184  up /= norm;
185  scalar_type coefficient = pai1->coeff_on_face(f, ii) * ctx1.J() * norm;
186  pf1->interpolation_grad(ctx1, coeff1, grad1, dim_type(qdim));
187  gmm::copy(grad1, E); gmm::add(gmm::transposed(grad1), E);
188  gmm::scale(E, 0.5);
189  scalar_type trace = gmm::mat_trace(E);
190  gmm::copy(gmm::identity_matrix(), S1);
191  gmm::scale(S1, lambda * trace);
192  gmm::add(gmm::scaled(E, 2*mu), S1);
193  gmm::mult(S1, up, jump);
194 
195  ERR[v.cv()] +=radius * coefficient * gmm::vect_norm2_sqr(jump);
196  eta2 += (radius * coefficient * gmm::vect_norm2_sqr(jump));
197  }
198  //cout << "Erreur en neumann sur " << v.cv() << " ERR[v.cv()] = " << ERR[v.cv()]-eee << endl;
199 
200  //if (ERR[v.cv()] > 100)
201  // cout << "Erreur en neumann sur element " << v.cv() << " : " << ERR[v.cv()] << endl;
202 
203  }
204 
205  };
206 
207 
208  {
209  int bnum = GAMMAC;
210 
211  getfem::mesh_region region = m.region(bnum);
212  for (getfem::mr_visitor v(region, m); !v.finished(); ++v) {
213 
214  bgeot::pgeometric_trans pgt1 = m.trans_of_convex(v.cv());
215  getfem::papprox_integration pai1 =
216  get_approx_im_or_fail(mim.int_method_of_element(v.cv()));
217  getfem::pfem pf1 = mf_u.fem_of_element(v.cv());
218 
219  m.points_of_convex(v.cv(), G1);
220 
221  coeff1.resize(mf_u.nb_basic_dof_of_element(v.cv()));
222  gmm::copy(gmm::sub_vector(U, gmm::sub_index(mf_u.ind_basic_dof_of_element(v.cv()))), coeff1);
223 
224  getfem::fem_interpolation_context ctx1(pgt1, pf1, base_node(N), G1, v.cv());
225 
226  // computation of h for gamma = gamma0*h
227  //scalar_type emax, emin, gamma;
228  //gmm::condition_number(ctx1.K(),emax,emin);
229  //gamma = gamma0 * emax * sqrt(scalar_type(N));
230  // Test autre gamma
231  scalar_type radius = m.convex_radius_estimate(v.cv());
232  scalar_type gamma=radius*gamma0;
233  short_type f = v.f();
234  for (unsigned ii=0; ii < pai1->nb_points_on_face(f); ++ii) {
235 
236  ctx1.set_xref(pai1->point_on_face(f, ii));
237  gmm::mult(ctx1.B(), pgt1->normals()[f], up);
238  scalar_type norm = gmm::vect_norm2(up);
239  up /= norm;
240  scalar_type coefficient = pai1->coeff_on_face(f, ii) * ctx1.J() * norm;
241  pf1->interpolation_grad(ctx1, coeff1, grad1, dim_type(qdim));
242  pf1->interpolation(ctx1, coeff1, U1, dim_type(qdim));
243  gmm::copy(grad1, E); gmm::add(gmm::transposed(grad1), E);
244  gmm::scale(E, 0.5);
245  scalar_type trace = gmm::mat_trace(E);
246  gmm::copy(gmm::identity_matrix(), S1);
247  gmm::scale(S1, lambda * trace);
248  gmm::add(gmm::scaled(E, 2*mu), S1);
249  gmm::mult(S1, up, sig); // sig = sigma(u)n
250  sign = gmm::vect_sp(sig,up);// sign = sigma_n(u)
251  Un = gmm::vect_sp(U1,up);// un = u_n
252  scal = Un-gamma*sign;
253  if (scal<0)
254  Pr = sign;
255  else
256  Pr = (scal/gamma + sign);
257  ERR[v.cv()] += coefficient*radius*Pr*Pr;
258  eta4 += coefficient*radius* Pr*Pr;
259 
260  gmm::copy(up,sigt);
261  gmm::scale(sigt, - sign);
262  gmm::add(sig,sigt);
263  ERR[v.cv()] += coefficient *radius*gmm::vect_norm2_sqr(sigt);
264  eta3 += coefficient *radius*gmm::vect_norm2_sqr(sigt);
265  }
266  // cout << "Erreur en contact sur " << v.cv() << " ERR[v.cv()] = " << ERR[v.cv()]-ee << endl;
267 
268  // if (ERR[v.cv()] > 100)
269  // cout << "Erreur en contact sur element " << v.cv() << " : " << ERR[v.cv()] << endl;
270 
271  }
272 
273  }
274 
275  cout << "eta1, eta2, eta3, eta4 = " << endl;
276  cout << sqrt(eta1) << endl;
277  cout << sqrt(eta2) << endl;
278  cout << sqrt(eta3) << endl;
279  cout << sqrt(eta4) << endl;
280  cout << sqrt(eta1+eta2+eta3+eta4) << endl;
281 
282  }
283 
284 #endif
285 
286 }
287 
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:63
bgeot::geotrans_inv_convex::invert
bool invert(const base_node &n, base_node &n_ref, scalar_type IN_EPS=1e-12, bool project_into_element=false)
given the node on the real element, returns the node on the reference element (even if it is outside ...
Definition: bgeot_geotrans_inv.cc:60
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
getfem::mesh_region::visitor
"iterator" class for regions.
Definition: getfem_mesh_region.h:236
getfem::mesh_region
structure used to hold a set of convexes and/or convex faces.
Definition: getfem_mesh_region.h:54
getfem_contact_and_friction_common.h
Comomon tools for unilateral contact and Coulomb friction bricks.
getfem_error_estimate.h
Definition of a posteriori error estimates.
getfem_mesher.h
An experimental mesher.
gmm::mult_add
void mult_add(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1790
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:556
gmm::mat_trace
linalg_traits< M >::value_type mat_trace(const M &m)
Trace of a matrix.
Definition: gmm_blas.h:527
gmm::vect_norm2_sqr
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2_sqr(const V &v)
squared Euclidean norm of a vector.
Definition: gmm_blas.h:543
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm::vect_sp
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*‍/
Definition: gmm_blas.h:263
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:243
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:185
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46