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