Rheolef  7.1
an efficient C++ finite element environment
piola_on_pointset.cc
Go to the documentation of this file.
1 //
21 #include "rheolef/piola_on_pointset.h"
22 #include "rheolef/piola_util.h"
23 
24 namespace rheolef {
25 
26 // ----------------------------------------------------------------------------
27 // initializers
28 // ----------------------------------------------------------------------------
29 template<class T>
30 void
32  const basis_basic<T>& piola_basis,
33  const quadrature<T>& quad,
34  const integrate_option& iopt)
35 {
36  _bops.set (quad, piola_basis);
37  _ignore_sys_coord = iopt.ignore_sys_coord;
38  _last_visited_geo.fill ("");
39  _last_visited_dis_ie.fill (std::numeric_limits<size_type>::max());
40 }
41 template<class T>
42 void
44  const basis_basic<T>& piola_basis,
45  const basis_basic<T>& nodal_basis,
46  const integrate_option& iopt)
47 {
48  _bops.set (nodal_basis, piola_basis);
49  _ignore_sys_coord = iopt.ignore_sys_coord;
50  _last_visited_geo.fill ("");
51  _last_visited_dis_ie.fill (std::numeric_limits<size_type>::max());
52 }
53 // ----------------------------------------------------------------------------
54 // piola transformation on a pointset
55 // ----------------------------------------------------------------------------
56 template<class T>
57 template<class M>
58 void
60 {
61  // memorize the last K call: avoid multiple _piola evaluation
62  reference_element hat_K = K;
63  if (_last_visited_geo [hat_K.variant()] == omega_K.name() &&
64  _last_visited_dis_ie [hat_K.variant()] == K.dis_ie()) {
65  return;
66  }
67  _last_visited_geo [hat_K.variant()] = omega_K.name();
68  _last_visited_dis_ie [hat_K.variant()] = K.dis_ie();
69 
70  // here, compute piola and weight on K:
71  Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = _piola [hat_K.variant()];
72  Eigen::Matrix<T,Eigen::Dynamic,1>& weight = _weight [hat_K.variant()];
73 
74  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
75  s_phij_xi = _bops.template evaluate<T> (hat_K);
76  const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>&
77  grad_phij_xi = _bops.template grad_evaluate<point_basic<T>> (hat_K);
78  size_type loc_nnod = s_phij_xi.rows();
79  size_type loc_ndof = s_phij_xi.cols();
80  size_type d = omega_K.dimension();
81  size_type map_d = hat_K.dimension();
82  space_constant::coordinate_type sys_coord = omega_K.coordinate_system();
83 
84  // pass 0: reset
85  piola .resize (loc_nnod);
86  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
87  piola[loc_inod].clear();
88  piola[loc_inod].d = d;
89  piola[loc_inod].map_d = map_d;
90  piola[loc_inod].sys_coord = sys_coord;
91  piola[loc_inod].ignore_sys_coord = _ignore_sys_coord;
92  }
93  // pass 1: F and DF
94  omega_K.dis_inod (K, _dis_inod_K);
95  for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
96  // dis_node: in outer loop: could require more time with external node
97  const point_basic<T>& xjnod = omega_K.dis_node (_dis_inod_K[loc_jdof]);
98  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
99  // step 1.1: F(hat_x)
100  for (size_type alpha = 0; alpha < d; alpha++) {
101  piola[loc_inod].F[alpha] += s_phij_xi (loc_inod,loc_jdof)*xjnod[alpha];
102  }
103  // step 1.2: DF(hat_x)
104  cumul_otimes (piola[loc_inod].DF, xjnod, grad_phij_xi(loc_inod,loc_jdof), d, map_d);
105 
106  // step 1.3: surface projector
107  if (map_d > 0 && map_d + 1 == d) {
108  map_projector (piola[loc_inod].DF, d, map_d, piola[loc_inod].P);
109  }
110  }
111  }
112  // pass 2: invDF and detDF
113  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
114  piola[loc_inod].invDF = pseudo_inverse_jacobian_piola_transformation (piola[loc_inod].DF, d, map_d);
115  piola[loc_inod].detDF = det_jacobian_piola_transformation (piola[loc_inod].DF, d, map_d);
116  }
117  // pass 3: quadrature weight, when pointset is a quadrature
118  if (!_bops.has_quadrature()) return;
119  weight.resize (loc_nnod);
120  const quadrature<T>& quad = _bops.get_quadrature();
121  size_type i_comp_axi = (omega_K.coordinate_system() == space_constant::axisymmetric_rz) ? 0 : 1;
122  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
123  T weight_q = quad (hat_K, loc_inod).w;
124  weight[loc_inod] = piola[loc_inod].detDF*weight_q;
125  if (omega_K.coordinate_system() != space_constant::cartesian && !_ignore_sys_coord) {
126  weight[loc_inod] *= piola[loc_inod].F [i_comp_axi];
127  }
128  }
129 }
130 // ----------------------------------------------------------------------------
131 // instanciation in library
132 // ----------------------------------------------------------------------------
133 
134 #define _RHEOLEF_instanciation(T) \
135 template class piola_on_pointset_rep<T>; \
136 
137 #define _RHEOLEF_instanciation_update(T,M) \
138 template void piola_on_pointset_rep<T>::_update ( \
139  const geo_basic<T,M>& omega_K, \
140  const geo_element& K) const; \
141 
144 #ifdef _RHEOLEF_HAVE_MPI
146 #endif // _RHEOLEF_HAVE_MPI
147 
148 } // namespace rheolef
rheolef::geo_basic
generic mesh with rerefence counting
Definition: domain_indirect.h:64
rheolef::piola::ignore_sys_coord
bool ignore_sys_coord
Definition: piola.h:78
rheolef::piola::map_d
size_type map_d
Definition: piola.h:76
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
rheolef::point_basic
Definition: point.h:87
rheolef::piola_on_pointset_rep::initialize
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
Definition: piola_on_pointset.cc:31
rheolef::piola::invDF
tensor_basic< T > invDF
Definition: piola.h:80
rheolef::piola::F
point_basic< T > F
Definition: piola.h:79
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float
rheolef::space_constant::cartesian
Definition: space_constant.h:122
rheolef::piola::sys_coord
coordinate_type sys_coord
Definition: piola.h:77
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::space_constant::axisymmetric_rz
Definition: space_constant.h:123
rheolef::det_jacobian_piola_transformation
T det_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
Definition: piola_util.cc:137
rheolef::piola_on_pointset_rep::_update
void _update(const geo_basic< T, M > &omega, const geo_element &K) const
Definition: piola_on_pointset.cc:59
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::basis_basic
Definition: basis.h:206
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::map_projector
void map_projector(const tensor_basic< T > &DF, size_t d, size_t map_d, tensor_basic< T > &P)
Definition: piola_util.cc:368
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::integrate_option::ignore_sys_coord
bool ignore_sys_coord
Definition: integrate_option.h:168
rheolef::piola::detDF
T detDF
Definition: piola.h:81
Float
see the Float page for the full documentation
mkgeo_grid.sys_coord
sys_coord
Definition: mkgeo_grid.sh:171
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
_RHEOLEF_instanciation_update
#define _RHEOLEF_instanciation_update(T, M)
Definition: piola_on_pointset.cc: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::piola::clear
void clear()
Definition: piola.h:104
rheolef::piola::d
size_type d
Definition: piola.h:76
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::space_constant::coordinate_type
coordinate_type
Definition: space_constant.h:121
rheolef::quadrature
Definition: quadrature.h:185
rheolef::piola_on_pointset_rep::size_type
reference_element::size_type size_type
Definition: piola_on_pointset.h:42
T
Expr1::float_type T
Definition: field_expr.h:218
rheolef::piola
Definition: piola.h:67