Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
23#include "piola_fem.h"
24namespace rheolef {
25
26// ----------------------------------------------------------------------------
27// grad-post treatment: tensor symmetry & surface projector
28// ----------------------------------------------------------------------------
29// vector-valued grad:
30template<class T>
31static
32void
33grad_post (
34 const piola<T>& p,
35 const details::differentiate_option& gopt,
36 const T& u,
37 point_basic<T>& grad_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:
46template<class T>
47static
48void
49grad_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:
85template<class T>
86static
87void
88grad_post (
89 const piola<T>& p,
90 const details::differentiate_option& gopt,
91 const tensor_basic<T>& u,
92 tensor3_basic<T>& grad_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
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition csr.h:455
Definition sphere.icc:25
size_t d
Definition sphere.icc:32
Definition leveque.h:25