GetFEM  5.4.2
getfem_derivatives.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2002-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_derivatives.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>,
34  @author Julien Pommier <Julien.Pommier@insa-toulouse.fr>
35  @date June 17, 2002.
36  @brief Compute the gradient of a field on a getfem::mesh_fem.
37 */
38 #ifndef GETFEM_DERIVATIVES_H__
39 #define GETFEM_DERIVATIVES_H__
40 
41 #include "getfem_mesh_fem.h"
42 #include "getfem_interpolation.h"
43 #include "gmm/gmm_dense_qr.h"
44 
45 namespace getfem
46 {
47  /** Compute the gradient of a field on a getfem::mesh_fem.
48  @param mf the source mesh_fem.
49  @param U the source field.
50  @param mf_target should be a lagrange discontinous element
51  @param V contains on output the gradient of U, evaluated on mf_target.
52 
53  mf_target may have the same Qdim than mf, or it may
54  be a scalar mesh_fem, in which case the derivatives are stored in
55  the order: DxUx,DyUx,DzUx,DxUy,DyUy,...
56 
57  in any case, the size of V should be
58  N*(mf.qdim)*(mf_target.nbdof/mf_target.qdim)
59  elements (this is not checked by the function!)
60  */
61  template<class VECT1, class VECT2>
62  void compute_gradient(const mesh_fem &mf, const mesh_fem &mf_target,
63  const VECT1 &UU, VECT2 &VV) {
64  typedef typename gmm::linalg_traits<VECT1>::value_type T;
65 
66  size_type N = mf.linked_mesh().dim();
67  size_type qdim = mf.get_qdim();
68  size_type target_qdim = mf_target.get_qdim();
69  size_type qqdimt = qdim * N / target_qdim;
70  std::vector<T> U(mf.nb_basic_dof());
71  std::vector<T> V(mf_target.nb_basic_dof() * qqdimt);
72 
73  mf.extend_vector(UU, U);
74 
75  GMM_ASSERT1(&mf.linked_mesh() == &mf_target.linked_mesh(),
76  "meshes are different.");
77  GMM_ASSERT1(target_qdim == qdim*N || target_qdim == qdim
78  || target_qdim == 1, "invalid Qdim for gradient mesh_fem");
79 
80  base_matrix G;
81  std::vector<T> coeff;
82 
83  bgeot::pgeotrans_precomp pgp = NULL;
84  pfem_precomp pfp = NULL;
85  pfem pf, pf_target, pf_old = NULL, pf_targetold = NULL;
87 
88  for (dal::bv_visitor cv(mf_target.convex_index()); !cv.finished(); ++cv) {
89  pf = mf.fem_of_element(cv);
90  pf_target = mf_target.fem_of_element(cv);
91  if (!pf) continue;
92 
93  GMM_ASSERT1(!(pf_target->need_G()) && pf_target->is_lagrange(),
94  "finite element target not convenient");
95 
96  bgeot::vectors_to_base_matrix(G, mf.linked_mesh().points_of_convex(cv));
97 
98  pgt = mf.linked_mesh().trans_of_convex(cv);
99  if (pf_targetold != pf_target) {
100  pgp = bgeot::geotrans_precomp(pgt, pf_target->node_tab(cv), pf_target);
101  }
102  pf_targetold = pf_target;
103 
104  if (pf_old != pf) {
105  pfp = fem_precomp(pf, pf_target->node_tab(cv), pf_target);
106  }
107  pf_old = pf;
108 
109  gmm::dense_matrix<T> grad(N,qdim), gradt(qdim,N);
110  fem_interpolation_context ctx(pgp, pfp, 0, G, cv, short_type(-1));
111  slice_vector_on_basic_dof_of_element(mf, U, cv, coeff);
112  // gmm::resize(coeff, mf.nb_basic_dof_of_element(cv));
113  // gmm::copy(gmm::sub_vector
114  // (U, gmm::sub_index(mf.ind_basic_dof_of_element(cv))), coeff);
115  for (size_type j = 0; j < pf_target->nb_dof(cv); ++j) {
116  size_type dof_t =
117  mf_target.ind_basic_dof_of_element(cv)[j*target_qdim] * qqdimt;
118  ctx.set_ii(j);
119  pf->interpolation_grad(ctx, coeff, gradt, dim_type(qdim));
120  gmm::copy(gmm::transposed(gradt),grad);
121  std::copy(grad.begin(), grad.end(), V.begin() + dof_t);
122  }
123  }
124 
125  mf_target.reduce_vector(V, VV);
126  }
127 
128  /** Compute the hessian of a field on a getfem::mesh_fem.
129  @param mf the source mesh_fem.
130  @param U the source field.
131  @param mf_target should be a lagrange discontinous element
132  does not work with vectorial elements. ... to be done ...
133  @param V contains on output the gradient of U, evaluated on mf_target.
134 
135  mf_target may have the same Qdim than mf, or it may
136  be a scalar mesh_fem, in which case the derivatives are stored in
137  the order: DxxUx,DxyUx, DyxUx, DyyUx, ...
138 
139  in any case, the size of V should be
140  N*N*(mf.qdim)*(mf_target.nbdof/mf_target.qdim)
141  elements (this is not checked by the function!)
142  */
143  template<class VECT1, class VECT2>
144  void compute_hessian(const mesh_fem &mf, const mesh_fem &mf_target,
145  const VECT1 &UU, VECT2 &VV) {
146  typedef typename gmm::linalg_traits<VECT1>::value_type T;
147 
148  size_type N = mf.linked_mesh().dim();
149  size_type qdim = mf.get_qdim();
150  size_type target_qdim = mf_target.get_qdim();
151  size_type qqdimt = qdim * N * N / target_qdim;
152  std::vector<T> U(mf.nb_basic_dof());
153  std::vector<T> V(mf_target.nb_basic_dof() * qqdimt);
154 
155  mf.extend_vector(UU, U);
156 
157  GMM_ASSERT1(&mf.linked_mesh() == &mf_target.linked_mesh(),
158  "meshes are different.");
159  GMM_ASSERT1(target_qdim == qdim || target_qdim == 1,
160  "invalid Qdim for gradient mesh_fem");
161  base_matrix G;
162  std::vector<T> coeff;
163 
164  bgeot::pgeotrans_precomp pgp = NULL;
165  pfem_precomp pfp = NULL;
166  pfem pf, pf_target, pf_old = NULL, pf_targetold = NULL;
168 
169  for (dal::bv_visitor cv(mf_target.convex_index()); !cv.finished(); ++cv) {
170  pf = mf.fem_of_element(cv);
171  pf_target = mf_target.fem_of_element(cv);
172  GMM_ASSERT1(!(pf_target->need_G()) && pf_target->is_lagrange(),
173  "finite element target not convenient");
174 
175  bgeot::vectors_to_base_matrix(G, mf.linked_mesh().points_of_convex(cv));
176 
177  pgt = mf.linked_mesh().trans_of_convex(cv);
178  if (pf_targetold != pf_target) {
179  pgp = bgeot::geotrans_precomp(pgt, pf_target->node_tab(cv), pf_target);
180  }
181  pf_targetold = pf_target;
182 
183  if (pf_old != pf) {
184  pfp = fem_precomp(pf, pf_target->node_tab(cv), pf_target);
185  }
186  pf_old = pf;
187 
188  gmm::dense_matrix<T> hess(N*N,qdim), hesst(qdim,N*N);
189  fem_interpolation_context ctx(pgp, pfp, 0, G, cv, short_type(-1));
190  slice_vector_on_basic_dof_of_element(mf, U, cv, coeff);
191  // gmm::resize(coeff, mf.nb_basic_dof_of_element(cv));
192  // gmm::copy(gmm::sub_vector
193  // (U, gmm::sub_index(mf.ind_basic_dof_of_element(cv))), coeff);
194  for (size_type j = 0; j < pf_target->nb_dof(cv); ++j) {
195  size_type dof_t
196  = mf_target.ind_basic_dof_of_element(cv)[j*target_qdim] * qqdimt;
197  ctx.set_ii(j);
198  pf->interpolation_hess(ctx, coeff, hesst, dim_type(qdim));
199  gmm::copy(gmm::transposed(hesst), hess);
200  std::copy(hess.begin(), hess.end(), V.begin() + dof_t);
201  }
202  }
203 
204  mf_target.reduce_vector(V, VV);
205  }
206 
207  /**Compute the Von-Mises stress of a field (only valid for
208  linearized elasticity in 3D)
209  */
210  template <typename VEC1, typename VEC2>
212  const getfem::mesh_fem &mf_vm,
213  const VEC1 &U, VEC2 &VM,
214  scalar_type mu=1) {
215  dal::bit_vector bv; bv.add(0, mf_vm.nb_dof());
216  interpolation_von_mises(mf_u, mf_vm, U, VM, bv, mu);
217  }
218 
219  template <typename VEC1, typename VEC2>
220  void interpolation_von_mises(const getfem::mesh_fem &mf_u,
221  const getfem::mesh_fem &mf_vm,
222  const VEC1 &U, VEC2 &VM,
223  const dal::bit_vector &mf_vm_dofs,
224  scalar_type mu=1) {
225 
226  assert(mf_vm.get_qdim() == 1);
227  unsigned N = mf_u.linked_mesh().dim(); assert(N == mf_u.get_qdim());
228  std::vector<scalar_type> DU(mf_vm.nb_dof() * N * N);
229 
230  getfem::compute_gradient(mf_u, mf_vm, U, DU);
231 
232  GMM_ASSERT1(!mf_vm.is_reduced(), "Sorry, to be done");
233 
234  scalar_type vm_min = 0., vm_max = 0.;
235  for (dal::bv_visitor i(mf_vm_dofs); !i.finished(); ++i) {
236  VM[i] = 0;
237  scalar_type sdiag = 0.;
238  for (unsigned j=0; j < N; ++j) {
239  sdiag += DU[i*N*N + j*N + j];
240  for (unsigned k=0; k < N; ++k) {
241  scalar_type e = .5*(DU[i*N*N + j*N + k] + DU[i*N*N + k*N + j]);
242  VM[i] += e*e;
243  }
244  }
245  VM[i] -= 1./N * sdiag * sdiag;
246  vm_min = (i == 0 ? VM[0] : std::min(vm_min, VM[i]));
247  vm_max = (i == 0 ? VM[0] : std::max(vm_max, VM[i]));
248  }
249  cout << "Von Mises : min=" << 4*mu*mu*vm_min << ", max="
250  << 4*mu*mu*vm_max << "\n";
251  gmm::scale(VM, 4*mu*mu);
252  }
253 
254 
255  /** Compute the Von-Mises stress of a field (valid for
256  linearized elasticity in 2D and 3D)
257  */
258  template <typename VEC1, typename VEC2, typename VEC3>
260  const getfem::mesh_fem &mf_vm,
261  const VEC1 &U, VEC2 &VM,
262  const getfem::mesh_fem &mf_lambda,
263  const VEC3 &lambda,
264  const getfem::mesh_fem &mf_mu,
265  const VEC3 &mu,
266  bool tresca) {
267  assert(mf_vm.get_qdim() == 1);
268  typedef typename gmm::linalg_traits<VEC1>::value_type T;
269  size_type N = mf_u.get_qdim();
270  std::vector<T> GRAD(mf_vm.nb_dof()*N*N),
271  LAMBDA(mf_vm.nb_dof()), MU(mf_vm.nb_dof());
272  base_matrix sigma(N,N);
273  base_vector eig(N);
274  if (tresca) interpolation(mf_lambda, mf_vm, lambda, LAMBDA);
275  interpolation(mf_mu, mf_vm, mu, MU);
276  compute_gradient(mf_u, mf_vm, U, GRAD);
277 
278  GMM_ASSERT1(!mf_vm.is_reduced(), "Sorry, to be done");
279  GMM_ASSERT1(N>=2 && N<=3, "Only for 2D and 3D");
280 
281  for (size_type i = 0; i < mf_vm.nb_dof(); ++i) {
282  scalar_type trE = 0, diag = 0;
283  for (unsigned j = 0; j < N; ++j)
284  trE += GRAD[i*N*N + j*N + j];
285  if (tresca)
286  diag = LAMBDA[i]*trE; // calculation of sigma
287  else
288  diag = (-2./3.)*MU[i]*trE; // for the calculation of deviator(sigma)
289  for (unsigned j = 0; j < N; ++j) {
290  for (unsigned k = 0; k < N; ++k) {
291  scalar_type eps = /* 0.5* */ (GRAD[i*N*N + j*N + k] +
292  GRAD[i*N*N + k*N + j]);
293  sigma(j,k) = /* 2* */ MU[i]*eps;
294  }
295  sigma(j,j) += diag;
296  }
297  if (!tresca) {
298  /* von mises: norm(deviator(sigma)) */
299  //gmm::add(gmm::scaled(Id, -gmm::mat_trace(sigma) / N), sigma);
300  if (N==3)
301  VM[i] = sqrt((3./2.)*gmm::mat_euclidean_norm_sqr(sigma));
302  else // for plane strains ( s_33 = -diag )
303  VM[i] = sqrt((3./2.)*(gmm::mat_euclidean_norm_sqr(sigma) + diag*diag));
304  } else {
305  /* else compute the tresca criterion */
306  gmm::symmetric_qr_algorithm(sigma, eig);
307  std::sort(eig.begin(), eig.end());
308  VM[i] = eig.back() - eig.front();
309  }
310  }
311  }
312 }
313 
314 #endif
getfem::mesh_fem::get_qdim
virtual dim_type get_qdim() const
Return the Q dimension.
Definition: getfem_mesh_fem.h:312
getfem::mesh_fem::is_reduced
bool is_reduced() const
Return true if a reduction matrix is applied to the dofs.
Definition: getfem_mesh_fem.h:195
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::compute_gradient
void compute_gradient(const mesh_fem &mf, const mesh_fem &mf_target, const VECT1 &UU, VECT2 &VV)
Compute the gradient of a field on a getfem::mesh_fem.
Definition: getfem_derivatives.h:62
getfem::fem_precomp
pfem_precomp fem_precomp(pfem pf, bgeot::pstored_point_tab pspt, dal::pstatic_stored_object dep)
Handles precomputations for FEM.
Definition: getfem_fem.cc:4333
getfem::mesh_fem::fem_of_element
virtual pfem fem_of_element(size_type cv) const
Return the basic fem associated with an element (if no fem is associated, the function will crash!...
Definition: getfem_mesh_fem.h:444
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:148
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
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
getfem::compute_hessian
void compute_hessian(const mesh_fem &mf, const mesh_fem &mf_target, const VECT1 &UU, VECT2 &VV)
Compute the hessian of a field on a getfem::mesh_fem.
Definition: getfem_derivatives.h:144
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
getfem::mesh_fem::convex_index
const dal::bit_vector & convex_index() const
Get the set of convexes where a finite element has been assigned.
Definition: getfem_mesh_fem.h:191
getfem::mesh_fem::nb_basic_dof
virtual size_type nb_basic_dof() const
Return the total number of basic degrees of freedom (before the optional reduction).
Definition: getfem_mesh_fem.h:556
getfem_mesh_fem.h
Define the getfem::mesh_fem class.
bgeot::geotrans_interpolation_context::set_ii
void set_ii(size_type ii__)
change the current point (assuming a geotrans_precomp_ is used)
Definition: bgeot_geometric_trans.h:463
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
getfem::interpolation_von_mises
void interpolation_von_mises(const getfem::mesh_fem &mf_u, const getfem::mesh_fem &mf_vm, const VEC1 &U, VEC2 &VM, scalar_type mu=1)
Compute the Von-Mises stress of a field (only valid for linearized elasticity in 3D)
Definition: getfem_derivatives.h:211
getfem::interpolation_von_mises_or_tresca
void interpolation_von_mises_or_tresca(const getfem::mesh_fem &mf_u, const getfem::mesh_fem &mf_vm, const VEC1 &U, VEC2 &VM, const getfem::mesh_fem &mf_lambda, const VEC3 &lambda, const getfem::mesh_fem &mf_mu, const VEC3 &mu, bool tresca)
Compute the Von-Mises stress of a field (valid for linearized elasticity in 2D and 3D)
Definition: getfem_derivatives.h:259
getfem::slice_vector_on_basic_dof_of_element
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
Definition: getfem_mesh_fem.h:659
getfem_interpolation.h
Interpolation of fields from a mesh_fem onto another.
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_fem::ind_basic_dof_of_element
virtual ind_dof_ct ind_basic_dof_of_element(size_type cv) const
Give an array of the dof numbers a of convex.
Definition: getfem_mesh_fem.h:450
gmm_dense_qr.h
Dense QR factorization.
getfem::mesh_fem::nb_dof
virtual size_type nb_dof() const
Return the total number of degrees of freedom.
Definition: getfem_mesh_fem.h:562