Rheolef  7.1
an efficient C++ finite element environment
piola_fem_grad_post.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_PIOLA_FEM_GRAD_POST_ICC
2 #define _RHEOLEF_PIOLA_FEM_GRAD_POST_ICC
3 #include "piola_fem.h"
24 namespace rheolef {
25 
26 // ----------------------------------------------------------------------------
27 // grad-post treatment: tensor symmetry & surface projector
28 // ----------------------------------------------------------------------------
29 // vector-valued grad:
30 template<class T>
31 static
32 void
33 grad_post (
34  const piola<T>& p,
35  const details::differentiate_option& gopt,
36  const T& u,
38 {
39  if (gopt.surfacic && p.map_d + 1 == p.d) {
40  // apply also the tangential projection P=(I-n*n)
41  // grad_s(u) = P*grad(u)
42  grad_u = p.P*grad_u; // TODO: DVT_OPTIM_2D
43  }
44 }
45 // tensor-valued grad:
46 template<class T>
47 static
48 void
49 grad_post (
50  const piola<T>& p,
51  const details::differentiate_option& gopt,
52  const point_basic<T>& u,
53  tensor_basic<T>& grad_u)
54 {
55  if (gopt.symmetrized) {
56  // compute D(u) from grad(u):
57  grad_u = (grad_u + trans(grad_u))/2; // TODO: DVT_OPTIM_2D
58  }
59  if (gopt.surfacic) {
60  // will apply also the tangential projection P=(I-n*n)
61  // surfacic grad, where P=I-nxn is the map projector
62  // grad_s(u) = grad(u)*P
63  // or
64  // Ds(u) = P*D(u)*P
65  if (! gopt.symmetrized) {
66  grad_u = grad_u*p.P; // TODO: DVT_OPTIM_2D
67  } else {
68  grad_u = p.P*grad_u*p.P; // TODO: DVT_OPTIM_2D
69  }
70  }
71  if (! p.ignore_sys_coord && p.sys_coord != space_constant::cartesian) {
72  // set the specific grad_u(theta,theta) = vr/r term in the axi case
73  size_t i_comp_r = (p.sys_coord == space_constant::axisymmetric_rz) ? 0 : 1;
74  size_t i_comp_theta = 2;
75  const T& ur = u [i_comp_r];
76  const T& r = p.F [i_comp_r];
77  // pb: r==0 when using gauss_lobatto, e.g. for the characteristic method
78  // e.g. in VF: 2D(u):D(v) dx -> (ur/r)*(vr/r)*r dr*dz = (ur*vr)/r dr*dz
79  // or for direct interpolation of grad(uh), D(uh), div(uh), etc in the axisym. case
80  check_macro (1+abs(r) != 1, "grad(): singular axisymmetric (1/r) weight (HINT: avoid interpolate() or change quadrature formulae)");
81  grad_u (i_comp_theta, i_comp_theta) = ur/r;
82  }
83 }
84 // tensor3-valued grad:
85 template<class T>
86 static
87 void
88 grad_post (
89  const piola<T>& p,
90  const details::differentiate_option& gopt,
91  const tensor_basic<T>& u,
93 {
94  // only used by u.grad(tau) = grad(tau)*u
95  // e.g. in doc/pexamples/transport_tensor_dg.cc : grad_h(sigma)*u
96  // From BirArmHas-vol1-1987 page 589 : formula for u.grad(tau)
97  // here, in zr or rz geometries, the u_theta component is zero
98  // and then there is no additional component to add here
99  // TODO: DVT_TENSOR3_GRAD_AXI otherwise, for a general grad(tau) tensor3 usage, there are additional components
100 }
101 
102 }// namespace rheolef
103 #endif // _RHEOLEF_PIOLA_FEM_GRAD_POST_ICC
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::space_constant::cartesian
@ cartesian
Definition: space_constant.h:122
p::d
size_t d
Definition: sphere.icc:32
rheolef::space_constant::axisymmetric_rz
@ axisymmetric_rz
Definition: space_constant.h:123
p
Definition: sphere.icc:25
rheolef::tensor3_basic< T >
tensor_basic< T > tensor3_basic< T >
Definition: piola_fem.h:137
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
u
Definition: leveque.h:25
grad_u
tensor grad_u
Definition: transport_tensor_exact.icc:26
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
piola_fem.h
rheolef::trans
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
grad_u
Definition: combustion_exact.icc:34
T
Expr1::float_type T
Definition: field_expr.h:218