field_evaluate.cc
Go to the documentation of this file.
//
// evaluate a field on a predefined point set: hat_x[q], q=0..nq
// See also piola_transformation.h
//
#include "rheolef/field_evaluate.h"
#include "rheolef/piola_util.h"
namespace rheolef {
// TODO: DVT_FEM_CLASS : should move in a
// fem_on_pointset & fem class
// fem: we known how to evaluate : RTk, etc
// fem_on_pointset: loop on a quadrature
//
// -----------------------------------------------------
// scalar-valued case:
// -----------------------------------------------------
template<class T, class M>
void
field_evaluate (
const field_basic<T,M>& uh,
const basis_on_pointset<T>& bops,
reference_element hat_K,
const std::vector<size_t>& dis_idof,
Eigen::Matrix<T,Eigen::Dynamic,1>& value)
{
typedef typename field_basic<T,M>::size_type size_type;
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
phij_xi = bops.template evaluate<T> (hat_K);
size_type loc_nnod = phij_xi.rows();
size_type loc_ndof = phij_xi.cols();
value.resize (loc_nnod);
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
const T& u_jdof = uh.dis_dof (dis_idof[loc_jdof]);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod] += phij_xi (loc_inod,loc_jdof)*u_jdof;
}
}
}
// -----------------------------------------------------
// vector-valued case:
// -----------------------------------------------------
template<class T, class M>
void
vector_field_evaluate (
const field_basic<T,M>& uh,
const basis_on_pointset<T>& bops,
reference_element hat_K,
const std::vector<size_t>& dis_idof_tab,
const basis_on_pointset<T>& piola_on_geo_basis,
std::vector<size_t>& dis_inod_geo,
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value)
{
typedef typename field_basic<T,M>::size_type size_type;
size_type dis_ndof = uh.dis_ndof();
size_type loc_ndof = dis_idof_tab.size();
if (uh.get_space().get_constitution().is_hierarchical()) {
// vector: as a product of scalar basis, e.g. (Pk)^d
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
phij_xi = bops.template evaluate<T> (hat_K);
size_type loc_nnod = phij_xi.rows();
size_type loc_comp_ndof = phij_xi.cols();
value.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod] = point_basic<T>(0,0,0);
}
size_type n_comp = uh.get_geo().dimension();
for (size_type loc_comp_jdof = 0; loc_comp_jdof < loc_comp_ndof; ++loc_comp_jdof) {
for (size_type k_comp = 0; k_comp < n_comp; k_comp++) {
size_type loc_jdof = loc_comp_jdof + k_comp*loc_comp_ndof;
assert_macro (loc_jdof < loc_ndof, "invalid local index "<<loc_jdof<<" out of range [0:"<<loc_ndof<<"[");
size_type dis_jdof = dis_idof_tab[loc_jdof];
assert_macro (dis_jdof < dis_ndof, "invalid distr index");
const T& u_jdof = uh.dis_dof (dis_jdof);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod][k_comp] += phij_xi(loc_inod,loc_jdof)*u_jdof;
}
}
}
} else { // non-hierarchical vector basis: e.g. Hdiv RTk requires piola vector transformation
check_macro (piola_on_geo_basis.is_set(), "un-set piola_on_geo_basis");
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>&
phij_xi = bops.template evaluate<point_basic<T>> (hat_K);
size_type loc_nnod = phij_xi.rows();
size_type loc_ndof = phij_xi.cols();
value.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod] = point_basic<T>(0,0,0);
}
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
size_type dis_jdof = dis_idof_tab[loc_jdof];
assert_macro (dis_jdof < dis_ndof, "invalid distributed index");
const T& u_jdof = uh.dis_dof (dis_jdof);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod] += phij_xi(loc_inod,loc_jdof)*u_jdof;
}
}
// piola_on_geo_basis : on evalue la base de la transf geo aux points xq
// c'est different de bops : evalue la base fem_RT aux points xq
// TODO: l'initialiser en amont de la boucle en q, dans la classe field_evaluate_base
// si on est vecteur et non-hierarchical
// basis_on_pointset<T> piola_on_geo_basis;
// piola_on_geo_basis.set (b, omega.get_piola_basis());
// std::vector<size_type> dis_inod_geo;
Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1> DF;
jacobian_piola_transformation (uh.get_geo(), piola_on_geo_basis, hat_K, dis_inod_geo, DF);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
T J = det_jacobian_piola_transformation (DF[loc_inod], uh.get_geo().dimension(), hat_K.dimension());
value[loc_inod] = (1/J)*(DF[loc_inod]*value[loc_inod]); // see RavTho-1977, eqn (3.17)
}
}
}
// -----------------------------------------------------
// tensor-valued case:
// -----------------------------------------------------
template<class T, class M>
void
tensor_field_evaluate (
const field_basic<T,M>& uh,
const basis_on_pointset<T>& bops,
reference_element hat_K,
const std::vector<size_t>& dis_idof_tab,
Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& value)
{
typedef typename field_basic<T,M>::size_type size_type;
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
phij_xi = bops.template evaluate<T> (hat_K);
size_type loc_nnod = phij_xi.rows();
size_type loc_comp_ndof = phij_xi.cols();
value.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod] = tensor_basic<T>();
}
size_type dis_ndof = uh.dis_ndof();
size_type d = uh.get_geo().dimension();
space_constant::coordinate_type sys_coord = uh.get_geo().coordinate_system();
size_type n_comp = space_constant::n_component (space_constant::tensor, d, sys_coord);
for (size_type loc_comp_jdof = 0; loc_comp_jdof < loc_comp_ndof; ++loc_comp_jdof) {
for (size_type kl_comp = 0; kl_comp < n_comp; kl_comp++) {
size_type loc_jdof = loc_comp_jdof + kl_comp*loc_comp_ndof;
assert_macro (loc_jdof < loc_ndof, "invalid local index "<<loc_jdof<<" out of range [0:"<<loc_ndof<<"[");
size_type dis_jdof = dis_idof_tab[loc_jdof];
assert_macro (dis_jdof < dis_ndof, "invalid distr index");
const T& u_jdof = uh.dis_dof (dis_jdof);
std::pair<size_type,size_type> kl
= space_constant::tensor_subscript (space_constant::tensor, sys_coord, kl_comp);
size_type k = kl.first;
size_type l = kl.second;
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
T tmp = phij_xi (loc_inod, loc_comp_jdof)*u_jdof;
value[loc_inod](k,l) += tmp;
if (k != l) value[loc_inod](l,k) += tmp;
}
}
}
}
// -----------------------------------------------------
// homogeneous multi-component case: get the i-th value
// -----------------------------------------------------
template<class T, class M>
void
field_component_evaluate (
const field_basic<T,M>& uh,
const basis_on_pointset<T>& bops,
reference_element hat_K,
const std::vector<size_t>& dis_idof_tab,
size_t k_comp,
Eigen::Matrix<T,Eigen::Dynamic,1>& value)
{
typedef typename field_basic<T,M>::size_type size_type;
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>&
phij_xi = bops.template evaluate<T> (hat_K);
size_type loc_nnod = phij_xi.rows();
size_type loc_comp_ndof = phij_xi.cols();
value.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod] = 0;
}
size_type dis_ndof = uh.dis_ndof();
size_type loc_ndof = dis_idof_tab.size();
size_type n_comp = uh.get_geo().dimension();
for (size_type loc_comp_jdof = 0; loc_comp_jdof < loc_comp_ndof; ++loc_comp_jdof) {
size_type loc_jdof = loc_comp_jdof + k_comp*loc_comp_ndof;
assert_macro (loc_jdof < loc_ndof, "invalid local index "<<loc_jdof<<" out of range [0:"<<loc_ndof<<"[");
size_type dis_jdof = dis_idof_tab[loc_jdof];
assert_macro (dis_jdof < dis_ndof, "invalid distr index");
const T& u_jdof = uh.dis_dof (dis_jdof);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
value[loc_inod] += phij_xi (loc_inod,loc_comp_jdof)*u_jdof;
}
}
}
// ----------------------------------------------------------------------------
// evaluate on a pointset: new version wih fops
// ----------------------------------------------------------------------------
template<class T, class M, class Value>
void
field_evaluate_continued (
const field_basic<T,M>& uh,
const geo_basic<T,M>& omega_K,
const geo_element& K,
const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& phij_xi,
Eigen::Matrix<Value,Eigen::Dynamic,1>& value)
{
using size_type = typename field_basic<T,M>::size_type;
// 1) extract uh dofs on K
// TODO: merge all cases in assembly_dis_idof instead of here
std::vector<size_type> dis_idof;
uh.get_space().get_constitution().assembly_dis_idof (omega_K, K, dis_idof);
size_type loc_ndof = dis_idof.size();
Eigen::Matrix<T,Eigen::Dynamic,1> udof (loc_ndof);
check_macro (size_type(phij_xi.cols()) == loc_ndof,
"invalid sizes phij_xi("<<phij_xi.rows()<<","<<phij_xi.cols()<<") and udof("<<udof.size()<<")");
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
size_type dis_jdof = dis_idof[loc_jdof];
udof [loc_jdof] = uh.dis_dof (dis_idof[loc_jdof]);
}
// 2) interpolate uh on the prescribed pointset
// TODO: DVT_EIGEN_BLAS2: value = phij_xi*udof; but Eigen compil pb when Value=point
size_type loc_nnod = phij_xi.rows();
value.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
Value sum = Value();
for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
sum += phij_xi (loc_inod,loc_jdof)*udof[loc_jdof];
}
value[loc_inod] = sum;
}
}
template<class T, class M, class Value>
void
field_evaluate (
const field_basic<T,M>& uh,
const fem_on_pointset<T>& fops,
const geo_basic<T,M>& omega_K,
const geo_element& K,
Eigen::Matrix<Value,Eigen::Dynamic,1>& value)
{
Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic> phij_xi;
details::differentiate_option none;
fops.template evaluate <M,Value,details::differentiate_option::none> (omega_K, K, none, phij_xi);
field_evaluate_continued (uh, omega_K, K, phij_xi, value);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation_value(T,M,Value) \
template void field_evaluate_continued ( \
const field_basic<T,M>& uh, \
const geo_basic<T,M>& omega_K, \
const geo_element& K, \
const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& phij_xi, \
Eigen::Matrix<Value,Eigen::Dynamic,1>& value); \
template void field_evaluate ( \
const field_basic<T,M>& uh, \
const fem_on_pointset<T>& fops, \
const geo_basic<T,M>& omega_K, \
const geo_element& K, \
Eigen::Matrix<Value,Eigen::Dynamic,1>& value); \
#define _RHEOLEF_instanciation(T,M) \
_RHEOLEF_instanciation_value(T,M,T) \
_RHEOLEF_instanciation_value(T,M,point_basic<T>) \
_RHEOLEF_instanciation_value(T,M,tensor_basic<T>) \
_RHEOLEF_instanciation_value(T,M,tensor3_basic<T>) \
_RHEOLEF_instanciation_value(T,M,tensor4_basic<T>) \
template \
void \
field_evaluate ( \
const field_basic<T,M>& uh, \
const basis_on_pointset<T>& bops, \
reference_element hat_K, \
const std::vector<size_t>& dis_idof, \
Eigen::Matrix<T,Eigen::Dynamic,1>& value); \
template \
void \
vector_field_evaluate<T,M> ( \
const field_basic<T,M>& uh, \
const basis_on_pointset<T>& bops, \
reference_element hat_K, \
const std::vector<size_t>& dis_idof_tab, \
const basis_on_pointset<T>& piola_on_geo_basis, \
std::vector<size_t>& dis_inod_geo, \
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& value); \
template \
void \
tensor_field_evaluate<T,M> ( \
const field_basic<T,M>& uh, \
const basis_on_pointset<T>& bops, \
reference_element hat_K, \
const std::vector<size_t>& dis_idof_tab, \
Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& value); \
template \
void \
field_component_evaluate<T,M> ( \
const field_basic<T,M>& uh, \
const basis_on_pointset<T>& bops, \
reference_element hat_K, \
const std::vector<size_t>& dis_idof_tab, \
size_t i_comp, \
Eigen::Matrix<T,Eigen::Dynamic,1>& value); \
_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI
} // namespace rheolef
size_type
field::size_type size_type
Definition: branch.cc:425
Float
see the Float page for the full documentation
rheolef::basis_on_pointset
Definition: basis_on_pointset.h:192
rheolef::basis_on_pointset::is_set
bool is_set() const
Definition: basis_on_pointset.h:291
rheolef::fem_on_pointset
Definition: fem_on_pointset.h:117
rheolef::field_basic< T, M >
rheolef::field_basic::get_space
const space_type & get_space() const
Definition: field.h:300
rheolef::field_basic::dis_dof
const T & dis_dof(size_type dis_idof) const
Definition: field.cc:376
rheolef::field_basic::dis_ndof
size_type dis_ndof() const
Definition: field.h:342
rheolef::field_basic::get_geo
const geo_type & get_geo() const
Definition: field.h:301
rheolef::geo_basic< T, M >
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::point_basic
Definition: point.h:87
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
rheolef::tensor_basic
Definition: tensor.h:90
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::distributed
distributed
Definition: asr.cc:228
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
rheolef::value
rheolef::std value
assert_macro
#define assert_macro(ok_condition, message)
Definition: dis_macros.h:113
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.tmp
tmp
Definition: mkgeo_ball.sh:380
mkgeo_grid.sys_coord
string sys_coord
Definition: mkgeo_grid.sh:171
rheolef::space_constant::tensor
@ tensor
Definition: space_constant.h:138
rheolef::space_constant::n_component
size_type n_component(valued_type valued_tag, size_type d, coordinate_type sys_coord)
Definition: space_constant.cc:136
rheolef::space_constant::coordinate_type
coordinate_type
Definition: space_constant.h:121
rheolef::space_constant::tensor_subscript
std::pair< size_type, size_type > tensor_subscript(valued_type valued_tag, coordinate_type sys_coord, size_type i_comp)
Definition: space_constant.cc:209
rheolef::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition: space_numbering.cc:187
rheolef::space_numbering::dis_ndof
size_type dis_ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:41
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::vector_field_evaluate
void vector_field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof_tab, const basis_on_pointset< T > &piola_on_geo_basis, std::vector< size_t > &dis_inod_geo, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &value)
Definition: field_evaluate.cc:64
rheolef::field_component_evaluate
void field_component_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof_tab, size_t k_comp, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value)
Definition: field_evaluate.cc:181
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::tensor_field_evaluate
void tensor_field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof_tab, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &value)
Definition: field_evaluate.cc:137
rheolef::field_evaluate
void field_evaluate(const field_basic< T, M > &uh, const basis_on_pointset< T > &bops, reference_element hat_K, const std::vector< size_t > &dis_idof, Eigen::Matrix< T, Eigen::Dynamic, 1 > &value)
Definition: field_evaluate.cc:39
rheolef::jacobian_piola_transformation
void jacobian_piola_transformation(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_pointset, reference_element hat_K, const std::vector< size_t > &dis_inod, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &DF)
Definition: piola_util.cc:80
rheolef::field_evaluate_continued
void field_evaluate_continued(const field_basic< T, M > &uh, const geo_basic< T, M > &omega_K, const geo_element &K, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &phij_xi, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value)
Definition: field_evaluate.cc:217
rheolef::_RHEOLEF_instanciation
_RHEOLEF_instanciation(Float, sequential, std::allocator< Float >) _RHEOLEF_instanciation(Float