piola_util.h

Go to the documentation of this file.

#ifndef _RHEOLEF_PIOLA_UTIL_H
#define _RHEOLEF_PIOLA_UTIL_H
// 
// piola transformation:
//
//      F : hat_K ---> K
//          hat_x +--> x = F(hat_x)
//
// i.e. map any geo_element K from a reference_element hat_K
//
#include "rheolef/geo.h"
#include "rheolef/basis_on_pointset.h"
namespace rheolef {
 
// ------------------------------------------
// piola transformation:
//          x[q] = F(hat_x[q])
// on all nodes of a quadrature formulae
// ------------------------------------------
template<class T, class M>
void
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<point_basic<T>,Eigen::Dynamic,1>& x);
 
// ------------------------------------------
// jacobian of the piola transformation
//         DF(hat_x[q])
// at an aritrarily point hat_x
// ------------------------------------------
template<class T, class M>
void
jacobian_piola_transformation (
  const geo_basic<T,M>&         omega,
  const basis_basic<T>&         piola_basis,
  reference_element             hat_K,
  const std::vector<size_t>&    dis_inod,
  const point_basic<T>&         hat_x,
        tensor_basic<T>&        DF);
 
// ------------------------------------------
// jacobian of the piola transformation
//         DF(hat_x[q])
// on all nodes of a quadrature formulae
// ------------------------------------------
template<class T, class M>
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);
 
 
template <class T>
T
det_jacobian_piola_transformation (const tensor_basic<T>& DF, size_t d , size_t map_d);
 
template<class T, class M>
point_basic<T>
normal_from_piola_transformation (const geo_basic<T,M>& omega, const geo_element& S, const tensor_basic<T>& DF, size_t d);
 
// The pseudo inverse extend inv(DF) for face in 3d or edge in 2d
//  i.e. useful for Laplacian-Beltrami and others surfacic forms.
//
// pinvDF (hat_xq) = inv DF, if tetra in 3d, tri in 2d, etc
//                  = pseudo-invese, when tri in 3d, edge in 2 or 3d
// e.g. on 3d face : pinvDF*DF = [1, 0, 0; 0, 1, 0; 0, 0, 0]
//
// let DF = [u, v, w], where u, v, w are the column vectors of DF
// then det(DF) = mixt(u,v,w)
// and det(DF)*inv(DF)^T = [v^w, w^u, u^v] where u^v = vect(u,v)
//
// application:
// if K=triangle(a,b,c) then u=ab=b-a, v=ac=c-a and w = n = u^v/|u^v|.
// Thus DF = [ab,ac,n] and det(DF)=|ab^ac|
// and inv(DF)^T = [ac^n/|ab^ac|, -ab^n/|ab^ac|, n]
// The pseudo-inverse is obtained by remplacing the last column n by zero.
//
template<class T>
tensor_basic<T>
pseudo_inverse_jacobian_piola_transformation (
    const tensor_basic<T>& DF,
    size_t d,
    size_t map_d);
 
// F^{-1}: x --> hat_x  on K
template<class T, class M>
point_basic<T>
inverse_piola_transformation (
  const geo_basic<T,M>&         omega,
  const reference_element&      hat_K,
  const std::vector<size_t>&    dis_inod,
  const point_basic<T>&         x);
 
// compute: P = I - nxn, the tangential projector on a map with unit normal n
template <class T>
void map_projector (const tensor_basic<T>& DF, size_t d, size_t map_d, tensor_basic<T>& P);
 
// axisymetric weight ?
// point_basic<T> xq = rheolef::piola_transformation (_omega, _piola_table, K, dis_inod, q);
template<class T>
T
weight_coordinate_system (space_constant::coordinate_type sys_coord, const point_basic<T>& xq);
 
// -------------------------------------------
// weight integration: w = det_DF*wq
// with optional axisymmetric r*dr factor
// -------------------------------------------
template<class T, class M>
void
piola_transformation_and_weight_integration (
  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,
  bool                                             ignore_sys_coord,
  Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& DF,
  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>&  x,
  Eigen::Matrix<T,Eigen::Dynamic,1>&               w);
 
}// namespace rheolef
#endif // _RHEOLEF_PIOLA_UTIL_H