GetFEM  5.4.2
getfem_crack_sif.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2007-2020 Yves Renard, Julien Pommier
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file getfem_crack_sif.h
33  @author Julien Pommier
34  @date July 2007
35  @brief crack support functions for computation of SIF
36  (stress intensity factors)
37 */
38 
39 #ifndef GETFEM_CRACK_SIF_H
40 #define GETFEM_CRACK_SIF_H
41 
42 #include "getfem/getfem_mesh.h"
46 
47 namespace getfem {
48  /* build a "ring" of convexes of given center and radius */
49  inline dal::bit_vector
50  build_sif_ring_from_mesh(const mesh &m,
51  base_node center, scalar_type r) {
52  dal::bit_vector bv;
53  scalar_type r2 = r - 1e-4;
54  unsigned in = 0;
55  for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
56  unsigned in1=0, out1=0;
57  unsigned in2=0, out2=0;
58  for (unsigned i=0; i < m.nb_points_of_convex(cv); ++i) {
59  base_node P = m.points_of_convex(cv)[i];
60  if (gmm::vect_dist2(P, center) < r) ++in1; else ++out1;
61  if (gmm::vect_dist2(P, center) < r2) ++in2; else ++out2;
62  }
63  if ((in1 && out1) || (in2 && out2)) bv.add(cv);
64  in += in1;
65  }
66  if (in < 3) GMM_WARNING1("looks like the radius is too small...");
67  return bv;
68  }
69 
70  /* return the crack tip in P,
71  and the outgoing tangent of the crack in T,
72  and the normal in N */
73  inline void get_crack_tip_and_orientation(const level_set &/* ls */,
74  base_node &P,
75  base_small_vector &T, base_small_vector &N) {
76  cerr << __PRETTY_FUNCTION__ << " IS TO BE DONE\n";
77  /* too lazy to do it now */
78  P.resize(2); P[0] = .5; P[1] = 0;
79  T.resize(2); T[0] = 1; T[1] = 0;
80  N.resize(2); N[0] = 0; N[1] = 1;
81  }
82 
83 
84  /* compute with great precision the stress intensity factors using
85  integral computed on a ring around the crack tip */
86  template <typename VECT>
87  void compute_crack_stress_intensity_factors(const level_set &ls,
88  const mesh_im &mim,
89  const mesh_fem &mf,
90  const VECT &U,
91  scalar_type ring_radius,
92  scalar_type lambda, scalar_type mu,
93  scalar_type young_modulus,
94  scalar_type &KI, scalar_type &KII) {
95  const mesh &mring = mim.linked_mesh();
96  mesh_fem_global_function mf_mode(mring, 1);
97  mesh_fem mf_q(mring,1);
98 
99  std::vector<pglobal_function> cfun(4);
100  for (unsigned j=0; j < 4; ++j) {
101  auto s = std::make_shared<crack_singular_xy_function>(j);
102  cfun[j] = global_function_on_level_set(ls, s);
103  }
104  mf_mode.set_functions(cfun);
105  mf_mode.set_qdim(2);
106 
107  mf_q.set_classical_finite_element(1);
108 
109  base_node crack_tip;
110  base_small_vector T, N;
111  get_crack_tip_and_orientation(ls, crack_tip, T, N);
112 
113  dal::bit_vector cvring = build_sif_ring_from_mesh(mring, crack_tip,
114  ring_radius);
115 
116  /* fill the "q" ring field with a approximately linear field, equal to
117  1 on the inner boundary, and equal to zero on the outer boundary */
118  std::vector<scalar_type> q(mf_q.nb_basic_dof());
119  for (unsigned d = 0; d < mf_q.nb_basic_dof(); ++d) {
120  base_node P = mf_q.point_of_basic_dof(d);
121  q[d] = (gmm::vect_dist2(P, crack_tip) > ring_radius) ? 0 : 1;
122  }
123 
124  base_vector U_mode(mf_mode.nb_dof()); assert(U_mode.size() == 8);
125 
126  /* expression for SIF computation taken from "a finite element
127  method for crack growth without remeshing", moes, dolbow &
128  belytschko */
129 
130  generic_assembly
131  assem("lambda=data$1(1); mu=data$2(1); x1=data$3(mdim(#1)); "
132  "U1=data$4(#1); U2=data$5(#2); q=data$6(#3);"
133  "t=U1(i).U2(j).q(k).comp(vGrad(#1).vGrad(#2).Grad(#3))(i,:,:,j,:,:,k,:);"
134  "e1=(t{1,2,:,:,:}+t{2,1,:,:,:})/2;"
135  "e2=(t{:,:,3,4,:}+t{:,:,4,3,:})/2;"
136  "e12=(e1{:,:,3,4,:}+e1{:,:,4,3,:})/2;"
137  "V()+=2*mu(p).e1(i,j,i,k,j).x1(k) + lambda(p).e1(i,i,j,k,j).x1(k);"
138  "V()+=2*mu(p).e2(i,k,i,j,j).x1(k) + lambda(p).e2(j,k,i,i,j).x1(k);"
139  "V()+=-2*mu(p).e12(i,j,i,j,k).x1(k) - lambda(p).e12(i,i,j,j,k).x1(k);");
140  assem.push_mf(mf);
141  assem.push_mf(mf_mode);
142  assem.push_mf(mf_q);
143  assem.push_mi(mim);
144  base_vector vlambda(1); vlambda[0] = lambda;
145  base_vector vmu(1); vmu[0] = mu;
146  assem.push_data(vlambda);
147  assem.push_data(vmu);
148  assem.push_data(T); // outgoing tangent of the crack
149  assem.push_data(U);
150  assem.push_data(U_mode);
151  assem.push_data(q);
152  base_vector V(1);
153  assem.push_vec(V);
154 
155  /* fill with the crack opening mode I or mode II */
156  for (unsigned mode = 1; mode <= 2; ++mode) {
157  base_vector::iterator it = U_mode.begin();
158  scalar_type coeff=0.;
159  switch(mode) {
160  case 1: {
161  scalar_type A=2+2*mu/(lambda+2*mu), B=-2*(lambda+mu)/(lambda+2*mu);
162  /* "colonne" 1: ux, colonne 2: uy */
163  *it++ = 0; *it++ = A-B; /* sin(theta/2) */
164  *it++ = A+B; *it++ = 0; /* cos(theta/2) */
165  *it++ = -B; *it++ = 0; /* sin(theta/2)*sin(theta) */
166  *it++ = 0; *it++ = B; /* cos(theta/2)*cos(theta) */
167  coeff = 1/sqrt(2*M_PI);
168  } break;
169  case 2: {
170  scalar_type C1 = (lambda+3*mu)/(lambda+mu);
171  *it++ = C1+2-1; *it++ = 0;
172  *it++ = 0; *it++ = -(C1-2+1);
173  *it++ = 0; *it++ = 1;
174  *it++ = 1; *it++ = 0;
175  coeff = 2*(mu+lambda)/(lambda+2*mu)/sqrt(2*M_PI);
176  } break;
177  }
178  gmm::scale(U_mode, coeff/young_modulus);
179  V[0] = 0.;
180  cout << "assemblig SIFs ..." << std::flush;
181  double time = gmm::uclock_sec();
182  assem.assembly(cvring);
183  cout << "done (" << gmm::uclock_sec()-time << " sec)\n";
184  V[0] *= young_modulus/2; /* plane stress only, should be E/(2*(1-nu)) for plane strain */
185  if (mode == 1) KI = V[0]; else KII = V[0];
186  }
187  }
188 
189  /* return a (very rough) estimate of the stress intensity factors using
190  the local displacement near the crack tip */
191  template <typename VECT>
192  void estimate_crack_stress_intensity_factors(const level_set &ls,
193  const mesh_fem &mf, const VECT &U,
194  scalar_type young_modulus,
195  scalar_type &KI,
196  scalar_type &KII, double EPS=1e-2) {
197  base_node P(2);
198  base_small_vector T(2),N(2);
199  get_crack_tip_and_orientation(ls, P, T, N);
200  base_node P1 = P - EPS*T + EPS/100.*N;
201  base_node P2 = P - EPS*T - EPS/100.*N;
202  std::vector<double> V(4);
203  getfem::mesh_trans_inv mti(mf.linked_mesh());
204  mti.add_point(P1);
205  mti.add_point(P2);
206  cout << "P1 = " << P1 << ", P2=" << P2 << "\n";
207  base_matrix M;
208  getfem::interpolation(mf, mti, U, V, false);
209  KI = (V[1] - V[3])/sqrt(EPS/(2*M_PI)) * young_modulus / 8.0;
210  KII = (V[0] - V[2])/sqrt(EPS/(2*M_PI)) * young_modulus / 8.0;
211  }
212 }
213 
214 #endif // GETFEM_CRACK_SIF_H
getfem_assembling_tensors.h
Generic assembly implementation.
gmm::vect_dist2
number_traits< typename linalg_traits< V1 >::value_type >::magnitude_type vect_dist2(const V1 &v1, const V2 &v2)
Euclidean distance between two vectors.
Definition: gmm_blas.h:597
getfem_mesh.h
Define a getfem::getfem_mesh object.
getfem::interpolation
void interpolation(const mesh_fem &mf_source, const mesh_fem &mf_target, const VECTU &U, VECTV &V, int extrapolation=0, double EPS=1E-10, mesh_region rg_source=mesh_region::all_convexes(), mesh_region rg_target=mesh_region::all_convexes())
interpolation/extrapolation of (mf_source, U) on mf_target.
Definition: getfem_interpolation.h:693
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem_level_set.h
Define level-sets.
getfem_mesh_fem_global_function.h
Define a mesh_fem with base functions which are global functions given by the user.