Rheolef  7.1
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
3 //
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
rheolef::inv_piola::map_dim
size_t map_dim
Definition: inv_piola.h:51
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::inv_piola::size_type
value_type::size_type size_type
Definition: inv_piola.h:37
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
rheolef::reference_element::H
static const variant_type H
Definition: reference_element.h:81
rheolef::point_basic
Definition: point.h:87
rheolef::inv_piola::derivative_solve
value_type derivative_solve(const value_type &r) const
Definition: inv_piola.h:122
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
rheolef::dot
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
rheolef::inv_piola::reset
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
rheolef::point_basic::size_type
size_t size_type
Definition: point.h:92
rheolef::inv_piola::grad_value
Eigen::Matrix< value_type, Eigen::Dynamic, 1 > grad_value
Definition: inv_piola.h:57
rheolef::inv_piola::hat_K
reference_element hat_K
Definition: inv_piola.h:53
rheolef::value
rheolef::std value
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::tensor_basic
Definition: tensor.h:90
rheolef::inv_piola::inv_DF
tensor_basic< T > inv_DF
Definition: inv_piola.h:58
rheolef::float_type
typename float_traits< value_type >::type float_type
Definition: field_expr_recursive.h:501
rheolef::inv_piola::dim
size_t dim
Definition: inv_piola.h:51
rheolef::inv_piola::dual_space_norm
float_type dual_space_norm(const value_type &r) const
Definition: inv_piola.h:127
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
rheolef::inv_piola::initial
value_type initial() const
Definition: inv_piola.h:89
rheolef::inv_piola::update_derivative
void update_derivative(const value_type &hat_xh) const
Definition: inv_piola.h:112
rheolef::basis_basic
Definition: basis.h:532
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::inv_piola::space_norm
float_type space_norm(const value_type &hat_xh) const
Definition: inv_piola.h:132
rheolef::inv_piola::value
Eigen::Matrix< float_type, Eigen::Dynamic, 1 > value
Definition: inv_piola.h:56
rheolef::inv_piola::residue
value_type residue(const value_type &hat_xh) const
Definition: inv_piola.h:102
rheolef::inv_piola::b
basis_basic< T > b
Definition: inv_piola.h:52
rheolef::inv_piola::DF
tensor_basic< T > DF
Definition: inv_piola.h:58
rheolef::inv_piola::inv_piola
inv_piola()
Definition: inv_piola.h:61
rheolef::value_type
result_type value_type
Definition: field_expr_recursive.h:499
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::space_numbering::dis_inod
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)
Definition: space_numbering.cc:177
rheolef::inv_piola::derivative_trans_mult
value_type derivative_trans_mult(const value_type &r) const
Definition: inv_piola.h:142
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
rheolef::inv_piola::node
std::vector< value_type > node
Definition: inv_piola.h:54
rheolef::inv_piola::value_type
point_basic< T > value_type
Definition: inv_piola.h:36
rheolef::inv_piola::duality_product
float_type duality_product(const value_type &r, const value_type &s) const
Definition: inv_piola.h:137
rheolef::cumul_otimes
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
rheolef::pseudo_inverse_jacobian_piola_transformation
tensor_basic< T > pseudo_inverse_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
Definition: piola_util.cc:265
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::inv_piola::float_type
T float_type
Definition: inv_piola.h:35
rheolef::inv_piola
Definition: inv_piola.h:33
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
mkgeo_ball.dim
int dim
Definition: mkgeo_ball.sh:307
T
Expr1::float_type T
Definition: field_expr.h:261
rheolef::inv_piola::x0
value_type x0
Definition: inv_piola.h:55
rheolef::inv_piola::set_x
void set_x(const value_type &x1)
Definition: inv_piola.h:41