Rheolef  7.2
an efficient C++ finite element environment
inv_piola.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_INV_PIOLA_H
2 #define _RHEOLEF_INV_PIOLA_H
23 //
24 // invert piola tranformation on nonlinear elements :
25 // quadrangles, high-order, etc
26 // by using the newton generic algorithm
27 //
28 #include "rheolef/geo.h"
29 #include "rheolef/piola_util.h"
30 namespace rheolef {
31 
32 template<class T>
33 class inv_piola {
34 public:
35  typedef T float_type;
37  typedef typename value_type::size_type size_type;
38  inv_piola();
39  template<class M>
40  void reset (const geo_basic<T,M>& omega, const reference_element& hat_K, const std::vector<size_t>& dis_inod);
41  void set_x (const value_type& x1) { x0 = x1; }
42  value_type initial() const;
43  value_type residue (const value_type& hat_xh) const;
44  void update_derivative (const value_type& hat_xh) const;
45  value_type derivative_solve (const value_type& r) const;
47  float_type space_norm (const value_type& hat_xh) const;
48  float_type dual_space_norm (const value_type& r) const;
49  float_type duality_product (const value_type& r, const value_type& s) const;
50 protected:
51  size_t dim, map_dim;
54  std::vector<value_type> node;
56  mutable Eigen::Matrix<float_type,Eigen::Dynamic,1> value;
57  mutable Eigen::Matrix<value_type,Eigen::Dynamic,1> grad_value;
59 };
60 template<class T>
62 : dim(0),
63  map_dim(0),
64  b(),
65  hat_K(),
66  node(),
67  x0(),
68  value(),
69  grad_value(),
70  DF(),
71  inv_DF()
72 {
73 }
74 template<class T>
75 template<class M>
76 void
77 inv_piola<T>::reset (const geo_basic<T,M>& omega, const reference_element& hat_K1, const std::vector<size_t>& dis_inod) {
78  dim = omega.dimension();
79  b = omega.get_piola_basis();
80  hat_K = hat_K1;
81  map_dim = hat_K1.dimension();
82  node.resize (dis_inod.size());
83  for (size_t loc_inod = 0, loc_nnod = node.size(); loc_inod < loc_nnod; ++loc_inod) {
84  node[loc_inod] = omega.dis_node (dis_inod[loc_inod]);
85  }
86 }
87 template<class T>
90  switch (hat_K.variant()) {
91  case reference_element::e : return value_type(0.5);
92  case reference_element::t : return value_type(1/float_type(3),1/float_type(3));
93  case reference_element::q : return value_type(0,0);
94  case reference_element::T : return value_type(1/float_type(3),1/float_type(3),1/float_type(3));
95  case reference_element::P : return value_type(1/float_type(3),1/float_type(3),0);
96  case reference_element::H : return value_type(0,0,0);
97  }
98  return value_type(0);
99 }
100 template<class T>
102 inv_piola<T>::residue (const value_type& hat_x) const {
103  b.evaluate (hat_K, hat_x, value);
104  value_type r;
105  for (size_t loc_inod = 0, loc_nnod = node.size(); loc_inod < loc_nnod; ++loc_inod) {
106  r = r + value[loc_inod]*node[loc_inod];
107  }
108  return r - x0;
109 }
110 template<class T>
111 void
113  b.grad_evaluate (hat_K, hat_x, grad_value);
114  DF.reset();
115  for (size_t loc_inod = 0, loc_nnod = node.size(); loc_inod < loc_nnod; ++loc_inod) {
116  cumul_otimes(DF, node[loc_inod], grad_value[loc_inod], dim, map_dim);
117  }
119 }
120 template<class T>
123  return inv_DF*r;
124 }
125 template<class T>
128  return norm(r);
129 }
130 template<class T>
132 inv_piola<T>::space_norm (const value_type& hat_xh) const {
133  return norm(hat_xh);
134 }
135 template<class T>
138  return dot (r, s);
139 }
140 template<class T>
143  return DF.trans_mult(r);
144 }
145 
146 }// namespace rheolef
147 #endif // _RHEOLEF_PIOLA_H
tensor_basic< T > inv_DF
Definition: inv_piola.h:58
value_type derivative_trans_mult(const value_type &r) const
Definition: inv_piola.h:142
tensor_basic< T > DF
Definition: inv_piola.h:58
float_type duality_product(const value_type &r, const value_type &s) const
Definition: inv_piola.h:137
basis_basic< T > b
Definition: inv_piola.h:52
value_type derivative_solve(const value_type &r) const
Definition: inv_piola.h:122
point_basic< T > value_type
Definition: inv_piola.h:36
void set_x(const value_type &x1)
Definition: inv_piola.h:41
Eigen::Matrix< float_type, Eigen::Dynamic, 1 > value
Definition: inv_piola.h:56
value_type x0
Definition: inv_piola.h:55
void update_derivative(const value_type &hat_xh) const
Definition: inv_piola.h:112
reference_element hat_K
Definition: inv_piola.h:53
value_type::size_type size_type
Definition: inv_piola.h:37
float_type dual_space_norm(const value_type &r) const
Definition: inv_piola.h:127
std::vector< value_type > node
Definition: inv_piola.h:54
float_type space_norm(const value_type &hat_xh) const
Definition: inv_piola.h:132
void reset(const geo_basic< T, M > &omega, const reference_element &hat_K, const std::vector< size_t > &dis_inod)
Definition: inv_piola.h:77
value_type initial() const
Definition: inv_piola.h:89
value_type residue(const value_type &hat_xh) const
Definition: inv_piola.h:102
Eigen::Matrix< value_type, Eigen::Dynamic, 1 > grad_value
Definition: inv_piola.h:57
size_t size_type
Definition: point.h:91
see the reference_element page for the full documentation
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type T
static const variant_type P
static const variant_type t
typename float_traits< value_type >::type float_type
rheolef::std value
result_type value_type
Expr1::float_type T
Definition: field_expr.h:230
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
This file is part of Rheolef.
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
void cumul_otimes(tensor_basic< T > &t, const point_basic< T > &a, const point_basic< T > &b, size_t na, size_t nb)
Definition: tensor.cc:305
tensor_basic< T > pseudo_inverse_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
Definition: piola_util.cc:279