Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
24namespace rheolef {
25
26// ----------------------------------------------------------------------------
27// initializers
28// ----------------------------------------------------------------------------
29template<class T>
30void
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}
41template<class T>
42void
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// ----------------------------------------------------------------------------
56template<class T>
57template<class M>
58void
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) \
135template class piola_on_pointset_rep<T>; \
136
137#define _RHEOLEF_instanciation_update(T,M) \
138template 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
#define _RHEOLEF_instanciation(T, M, A)
Definition asr.cc:223
see the Float page for the full documentation
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
size_type dis_ie() const
see the integrate_option page for the full documentation
reference_element::size_type size_type
void initialize(const basis_basic< T > &piola_basis, const quadrature< T > &quad, const integrate_option &iopt)
void _update(const geo_basic< T, M > &omega, const geo_element &K) const
see the reference_element page for the full documentation
variant_type variant() const
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)
void map_projector(const tensor_basic< T > &DF, size_t d, size_t map_d, tensor_basic< T > &P)
T det_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
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
#define _RHEOLEF_instanciation_update(T, M)
size_type d
Definition piola.h:76
coordinate_type sys_coord
Definition piola.h:77
size_type map_d
Definition piola.h:76
tensor_basic< T > invDF
Definition piola.h:80
bool ignore_sys_coord
Definition piola.h:78
void clear()
Definition piola.h:104
point_basic< T > F
Definition piola.h:79