Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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"
30namespace rheolef {
31
32template<class T>
33class inv_piola {
34public:
35 typedef T float_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;
50protected:
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};
60template<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}
74template<class T>
75template<class M>
76void
77inv_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}
87template<class T>
90 switch (hat_K.variant()) {
91 case reference_element::e : return value_type(0.5);
93 case reference_element::q : return value_type(0,0);
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}
100template<class T>
102inv_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}
110template<class T>
111void
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 }
118 inv_DF = pseudo_inverse_jacobian_piola_transformation (DF, dim, map_dim);
119}
120template<class T>
123 return inv_DF*r;
124}
125template<class T>
128 return norm(r);
129}
130template<class T>
133 return norm(hat_xh);
134}
135template<class T>
138 return dot (r, s);
139}
140template<class T>
143 return DF.trans_mult(r);
144}
145
146}// namespace rheolef
147#endif // _RHEOLEF_PIOLA_H
generic mesh with rerefence counting
Definition geo.h:1089
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
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
Expr1::float_type T
Definition field_expr.h:230
This file is part of Rheolef.
tensor_basic< T > pseudo_inverse_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
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
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387