piola_on_pointset.cc
Go to the documentation of this file.
//
#include "rheolef/piola_on_pointset.h"
#include "rheolef/piola_util.h"
namespace rheolef {
// ----------------------------------------------------------------------------
// initializers
// ----------------------------------------------------------------------------
template<class T>
void
piola_on_pointset_rep<T>::initialize (
const basis_basic<T>& piola_basis,
const quadrature<T>& quad,
const integrate_option& iopt)
{
_bops.set (quad, piola_basis);
_ignore_sys_coord = iopt.ignore_sys_coord;
_last_visited_geo.fill ("");
_last_visited_dis_ie.fill (std::numeric_limits<size_type>::max());
}
template<class T>
void
piola_on_pointset_rep<T>::initialize (
const basis_basic<T>& piola_basis,
const basis_basic<T>& nodal_basis,
const integrate_option& iopt)
{
_bops.set (nodal_basis, piola_basis);
_ignore_sys_coord = iopt.ignore_sys_coord;
_last_visited_geo.fill ("");
_last_visited_dis_ie.fill (std::numeric_limits<size_type>::max());
}
// ----------------------------------------------------------------------------
// piola transformation on a pointset
// ----------------------------------------------------------------------------
template<class T>
template<class M>
void
piola_on_pointset_rep<T>::_update (const geo_basic<T,M>& omega_K, const geo_element& K) const
{
// memorize the last K call: avoid multiple _piola evaluation
reference_element hat_K = K;
if (_last_visited_geo [hat_K.variant()] == omega_K.name() &&
_last_visited_dis_ie [hat_K.variant()] == K.dis_ie()) {
return;
}
_last_visited_geo [hat_K.variant()] = omega_K.name();
_last_visited_dis_ie [hat_K.variant()] = K.dis_ie();
// here, compute piola and weight on K:
Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = _piola [hat_K.variant()];
Eigen::Matrix<T,Eigen::Dynamic,1>& weight = _weight [hat_K.variant()];
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
s_phij_xi = _bops.template evaluate<T> (hat_K);
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>&
grad_phij_xi = _bops.template grad_evaluate<point_basic<T>> (hat_K);
size_type loc_nnod = s_phij_xi.rows();
size_type loc_ndof = s_phij_xi.cols();
size_type d = omega_K.dimension();
size_type map_d = hat_K.dimension();
space_constant::coordinate_type sys_coord = omega_K.coordinate_system();
// pass 0: reset
piola .resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
piola[loc_inod].clear();
piola[loc_inod].d = d;
piola[loc_inod].map_d = map_d;
piola[loc_inod].sys_coord = sys_coord;
piola[loc_inod].ignore_sys_coord = _ignore_sys_coord;
}
// pass 1: F and DF
omega_K.dis_inod (K, _dis_inod_K);
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
// dis_node: in outer loop: could require more time with external node
const point_basic<T>& xjnod = omega_K.dis_node (_dis_inod_K[loc_jdof]);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
// step 1.1: F(hat_x)
for (size_type alpha = 0; alpha < d; alpha++) {
piola[loc_inod].F[alpha] += s_phij_xi (loc_inod,loc_jdof)*xjnod[alpha];
}
// step 1.2: DF(hat_x)
cumul_otimes (piola[loc_inod].DF, xjnod, grad_phij_xi(loc_inod,loc_jdof), d, map_d);
// step 1.3: surface projector
if (map_d > 0 && map_d + 1 == d) {
map_projector (piola[loc_inod].DF, d, map_d, piola[loc_inod].P);
}
}
}
// pass 2: invDF and detDF
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
piola[loc_inod].invDF = pseudo_inverse_jacobian_piola_transformation (piola[loc_inod].DF, d, map_d);
piola[loc_inod].detDF = det_jacobian_piola_transformation (piola[loc_inod].DF, d, map_d);
}
// pass 3: quadrature weight, when pointset is a quadrature
if (!_bops.has_quadrature()) return;
weight.resize (loc_nnod);
const quadrature<T>& quad = _bops.get_quadrature();
size_type i_comp_axi = (omega_K.coordinate_system() == space_constant::axisymmetric_rz) ? 0 : 1;
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
T weight_q = quad (hat_K, loc_inod).w;
weight[loc_inod] = piola[loc_inod].detDF*weight_q;
if (omega_K.coordinate_system() != space_constant::cartesian && !_ignore_sys_coord) {
weight[loc_inod] *= piola[loc_inod].F [i_comp_axi];
}
}
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T) \
template class piola_on_pointset_rep<T>; \
#define _RHEOLEF_instanciation_update(T,M) \
template void piola_on_pointset_rep<T>::_update ( \
const geo_basic<T,M>& omega_K, \
const geo_element& K) const; \
_RHEOLEF_instanciation(Float)
_RHEOLEF_instanciation_update(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation_update(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
rheolef::geo_basic< T, M >
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
@ 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
@ 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