piola_util.cc
Go to the documentation of this file.
#include "rheolef/piola_util.h"
#include "rheolef/geo_domain.h"
#include "rheolef/inv_piola.h"
#include "rheolef/damped_newton.h"
namespace rheolef {
// =========================================================================
// part 1. piola transformation
// =========================================================================
// F_K : hat_K --> K
// hat_x --> x
//
// pre-evaluation of the piola basis on a predefined point set
// e.g. quadrature nodes hat_x[q], q=0..nq on hat_K
// then, fast transformation of all hat_x[q] into xq on any K
//
// x = F(hat_x) = sum_j phi_j(hat_x)*xjnod
// for all hat_x in the pointset : hat_xi
// where xjnod: j-th node of element K in mesh omega
//
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)
{
typedef typename geo_basic<T,M>::size_type size_type;
const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& phij_xi = piola_on_pointset.template evaluate<T> (hat_K);
size_type loc_nnod = phij_xi.rows();
size_type loc_ndof = phij_xi.cols();
x.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
x[loc_inod] = point_basic<T>(0,0,0);
}
size_type d = omega.dimension();
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.dis_node (dis_inod[loc_jdof]);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
for (size_type alpha = 0; alpha < d; alpha++) {
x[loc_inod][alpha] += phij_xi (loc_inod,loc_jdof)*xjnod[alpha];
}
}
}
}
// ------------------------------------------
// jacobian of the piola transformation
// at a all quadrature points
// ------------------------------------------
//
// DF(hat_x) = sum_j grad_phi_j(hat_x)*xjnod
// for all hat_x in the pointset : hat_xi
// where xjnod: j-th node of element K in mesh omega
//
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)
{
typedef typename geo_basic<T,M>::size_type size_type;
const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>&
grad_phij_xi = piola_on_pointset.template grad_evaluate<point_basic<T>>(hat_K);
size_type loc_nnod = grad_phij_xi.rows();
size_type loc_ndof = grad_phij_xi.cols();
DF.resize (loc_nnod);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
DF[loc_inod].fill (0);
}
size_type d = omega.dimension();
size_type map_d = hat_K.dimension();
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.dis_node (dis_inod[loc_jdof]);
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
cumul_otimes (DF[loc_inod], xjnod, grad_phij_xi(loc_inod,loc_jdof), d, map_d);
}
}
}
// ------------------------------------------
// jacobian of the piola transformation
// at an aritrarily point hat_x
// ------------------------------------------
//
// DF(hat_x) = sum_j grad_phi_j(hat_x)*xjnod
// where xjnod: j-th node of element K in mesh omega
//
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)
{
typedef typename geo_basic<T,M>::size_type size_type;
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> grad_phi_x;
piola_basis.grad_evaluate (hat_K, hat_x, grad_phi_x);
DF.fill (0);
size_type d = omega.dimension();
size_type map_d = hat_K.dimension();
for (size_type loc_jdof = 0, loc_ndof = dis_inod.size(); loc_jdof < loc_ndof; loc_jdof++) {
const point_basic<T>& xjnod = omega.dis_node (dis_inod[loc_jdof]);
cumul_otimes (DF, xjnod, grad_phi_x [loc_jdof], d, map_d);
}
}
template <class T>
T
det_jacobian_piola_transformation (const tensor_basic<T>& DF, size_t d , size_t map_d)
{
if (d == map_d) {
return DF.determinant (map_d);
}
/* surface jacobian: references:
* Spectral/hp element methods for CFD
* G. E. M. Karniadakis and S. J. Sherwin
* Oxford university press
* 1999
* page 165
*/
switch (map_d) {
case 0: return 1;
case 1: return norm(DF.col(0));
case 2: return norm(vect(DF.col(0), DF.col(1)));
default:
error_macro ("det_jacobian_piola_transformation: unsupported element dimension "
<< map_d << " in " << d << "D mesh.");
return 0;
}
}
//
// application:
//
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)
{
switch (d) {
case 1: { // point in 1D: DF[1][0] is empty, so scan S[0] node connectivity
typedef typename geo_basic<T,M>::size_type size_type;
if (S.dimension() + 1 == omega.map_dimension()) {
size_type dis_ie = S.master(0);
check_macro (dis_ie != std::numeric_limits<size_type>::max(), "normal: requires neighbours initialization");
const geo_element& K = omega.dis_get_geo_element (S.dimension()+1, dis_ie);
Float sign = (S[0] == K[1]) ? 1 : -1;
return point_basic<T>(sign);
}
// omega is a domain of sides, as "boundary" or "internal_sides":
// for the side orient, we need to go back to its backgound volumic mesh
if (omega.variant() == geo_abstract_base_rep<T>::geo_domain_indirect) {
size_type dis_isid = S.dis_ie();
size_type first_dis_isid = omega.sizes().ownership_by_dimension[S.dimension()].first_index();
size_type isid = dis_isid - first_dis_isid;
check_macro (dis_isid >= first_dis_isid, "unexpected dis_index "<<dis_isid<<": out of local range");
const geo_basic<T,M>* ptr_bgd_omega = 0;
if (omega.get_background_geo().map_dimension() == 1) {
const geo_basic<T,M>& bgd_omega = omega.get_background_geo();
const geo_element& bgd_S = bgd_omega.get_geo_element(0, isid);
size_type bgd_dis_ie = bgd_S.master(0);
check_macro (bgd_dis_ie != std::numeric_limits<size_type>::max(),
"normal: bgd_S.dis_ie={"<<bgd_S.dis_ie()<<"} without neighbours; requires neighbours initialization for mesh " << bgd_omega.name());
const geo_element& bgd_K = bgd_omega.dis_get_geo_element (bgd_S.dimension()+1, bgd_dis_ie);
Float sign = (bgd_S[0] == bgd_K[1]) ? 1 : -1;
return point_basic<T>(sign);
} else {
// get a 1D geometry at a higher depth, e.g. for HDG 0D "boundary" domain in "square[sides]" 0D geo_domain
const geo_basic<T,M>& bgd_omega = omega.get_background_geo();
const geo_basic<T,M>& bgd2_omega = bgd_omega.get_background_geo();
check_macro (bgd2_omega.dimension() == 1, "unsupported depth for "<<omega.name()<<" in background domain "<<bgd_omega.name());
const geo_element& bgd_S = bgd_omega.get_geo_element(0, isid);
const geo_element& bgd2_S = bgd_omega.dom2bgd_geo_element (bgd_S);
size_type bgd2_dis_ie = bgd2_S.master(0);
check_macro (bgd2_dis_ie != std::numeric_limits<size_type>::max(),
"normal: bgd2_S.dis_ie={"<<bgd2_S.dis_ie()<<"} without neighbours; requires neighbours initialization for mesh " << bgd2_omega.name());
const geo_element& bgd2_K = bgd2_omega.dis_get_geo_element (bgd2_S.dimension()+1, bgd2_dis_ie);
Float sign = (bgd2_S[0] == bgd2_K[1]) ? 1 : -1;
return point_basic<T>(sign);
}
} else {
size_type dis_isid = S.dis_ie();
size_type first_dis_isid = omega.sizes().ownership_by_dimension[S.dimension()].first_index();
size_type isid = dis_isid - first_dis_isid;
check_macro (dis_isid >= first_dis_isid, "unexpected dis_index "<<dis_isid<<": out of local range");
const geo_basic<T,M>& bgd_omega = omega.get_background_geo();
const geo_basic<T,M>& bgd_dom = omega.get_background_domain();
const geo_element& bgd_S = bgd_dom[isid]; // TODO: pas clair, differe du cas precedent ?
size_type bgd_dis_ie = bgd_S.master(0);
check_macro (bgd_dis_ie != std::numeric_limits<size_type>::max(),
"normal: bgd_S.dis_ie={"<<bgd_S.dis_ie()<<"} without neighbours; requires neighbours initialization for mesh " << bgd_omega.name());
const geo_element& bgd_K = bgd_omega.dis_get_geo_element (bgd_S.dimension()+1, bgd_dis_ie);
Float sign = (bgd_S[0] == bgd_K[1]) ? 1 : -1;
return point_basic<T>(sign);
}
}
case 2: { // edge in 2D
// 2D: S=edge(a,b) then t=b-a, DF=[t] and n = (t1,-t0) => det(n,1)=1 : (n,t) is direct
// and the normal goes outside on a boundary edge S, when the associated element K is well oriented
point_basic<T> t = DF.col(0);
t /= norm(t);
return point_basic<T>(t[1], -t[0]);
}
case 3: { // 3D: S=triangle(a,b,c) then t0=b-a, t1=c-a, DF=[t0,t1] and n = t0^t1/|t0^t1|.
point_basic<T> t0 = DF.col(0);
point_basic<T> t1 = DF.col(1);
point_basic<T> n = vect(t0,t1);
n /= norm(n);
return n;
}
default:
error_macro ("normal: unsupported " << d << "D mesh.");
return point_basic<T>();
}
}
// 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)
{
if (d == map_d) {
return inv(DF, map_d);
}
tensor_basic<T> invDF;
switch (map_d) {
case 0: { // point in 1D
invDF(0,0) = 1;
return invDF;
}
case 1: { // segment in 2D
point_basic<T> t = DF.col(0);
invDF.set_row (t/norm2(t), 0, d);
return invDF;
}
case 2: {
point_basic<T> t0 = DF.col(0);
point_basic<T> t1 = DF.col(1);
point_basic<T> n = vect(t0,t1);
T det2 = norm2(n);
point_basic<T> v0 = vect(t1,n)/det2;
point_basic<T> v1 = - vect(t0,n)/det2;
invDF.set_row (v0, 0, d);
invDF.set_row (v1, 1, d);
return invDF;
}
default:
error_macro ("pseudo_inverse_jacobian_piola_transformation: unsupported element dimension "
<< map_d << " in " << d << "D mesh.");
return invDF;
}
}
// 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)
{
switch (sys_coord) {
case space_constant::axisymmetric_rz: return xq[0];
case space_constant::axisymmetric_zr: return xq[1];
case space_constant::cartesian: return 1;
default: {
fatal_macro ("unsupported coordinate system `"
<< space_constant::coordinate_system_name (sys_coord) << "'");
return 0;
}
}
}
// -------------------------------------------
// 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_quad,
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)
{
typedef typename geo_basic<T,M>::size_type size_type;
jacobian_piola_transformation (omega, piola_on_quad, hat_K, dis_inod, DF);
size_type loc_nnod = piola_on_quad.nnod (hat_K);
w.resize (loc_nnod);
if (omega.coordinate_system() == space_constant::cartesian || ignore_sys_coord) {
w.fill (T(1));
} else {
piola_transformation (omega, piola_on_quad, hat_K, dis_inod, x);
size_t k = (omega.coordinate_system() == space_constant::axisymmetric_rz) ? 0 : 1;
for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
w[loc_inod] = x[loc_inod][k];
}
}
size_type d = omega.dimension();
size_type map_d = hat_K.dimension();
const quadrature<T>& quad = piola_on_quad.get_quadrature();
typename quadrature<T>::const_iterator
first_quad = quad.begin(hat_K),
last_quad = quad.end (hat_K);
for (size_type q = 0; first_quad != last_quad; ++first_quad, ++q) {
T det_DF = det_jacobian_piola_transformation (DF[q], d, map_d);
T wq = det_DF*(*first_quad).w;
if (! ignore_sys_coord) {
w[q] = wq*w[q];
} else {
w[q] = wq;
}
}
}
// calcul P = I - nxn
template<class T>
void
map_projector (
const tensor_basic<T>& DF,
size_t d,
size_t map_d,
tensor_basic<T>& P)
{
point_basic<Float> n;
switch (map_d) {
case 1: {
point_basic<Float> t = DF.col(0);
check_macro (d == map_d+1, "unexpected dimension map_d="<<map_d<<" and d="<<d);
n = point_basic<T>(t[1],-t[0]);
break;
}
case 2: {
point_basic<Float> t0 = DF.col(0);
point_basic<Float> t1 = DF.col(1);
n = vect(t0,t1);
break;
}
default:
error_macro ("unexpected dimension "<<map_d);
}
n = n/norm(n);
for (size_t l = 0; l < d; l++) {
for (size_t m = 0; m < d; m++) {
P(l,m) = - n[l]*n[m];
}
P(l,l) += 1;
}
}
// =========================================================================
// part 2. inverse piola transformation
// =========================================================================
// F_K^{-1} : K --> hat(K)
// x --> hat(x)
// TODO: non-linear case
template<class T>
static
inline
point_basic<T>
inv_piola_e (
const point_basic<T>& x,
const point_basic<T>& a,
const point_basic<T>& b)
{
return point_basic<T>((x[0]-a[0])/(b[0]-a[0]));
}
template<class T>
static
inline
point_basic<T>
inv_piola_t (
const point_basic<T>& x,
const point_basic<T>& a,
const point_basic<T>& b,
const point_basic<T>& c)
{
T t9 = 1/(-b[0]*c[1]+b[0]*a[1]+a[0]*c[1]+c[0]*b[1]-c[0]*a[1]-a[0]*b[1]);
T t11 = -a[0]+x[0];
T t15 = -a[1]+x[1];
return point_basic<T>((-c[1]+a[1])*t9*t11-(-c[0]+a[0])*t9*t15,
( b[1]-a[1])*t9*t11-( b[0]-a[0])*t9*t15);
}
template<class T>
static
inline
point_basic<T>
inv_piola_T (
const point_basic<T>& x,
const point_basic<T>& a,
const point_basic<T>& b,
const point_basic<T>& c,
const point_basic<T>& d)
{
tensor_basic<T> A;
point_basic<T> ax;
for (size_t i = 0; i < 3; i++) {
ax[i] = x[i]-a[i];
A(i,0) = b[i]-a[i];
A(i,1) = c[i]-a[i];
A(i,2) = d[i]-a[i];
}
tensor_basic<T> inv_A;
bool is_singular = ! invert_3x3 (A, inv_A);
check_macro(!is_singular, "inv_piola: singular transformation in a tetrahedron");
point_basic<T> hat_x = inv_A*ax;
return hat_x;
}
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)
{
check_macro (omega.order() == 1, "inverse piola: mesh order > 1: not yet");
if (omega.order() == 1) {
switch (hat_K.variant()) {
case reference_element::e: return inv_piola_e (x, omega.dis_node(dis_inod [0]),
omega.dis_node(dis_inod [1]));
case reference_element::t: return inv_piola_t (x, omega.dis_node(dis_inod [0]),
omega.dis_node(dis_inod [1]),
omega.dis_node(dis_inod [2]));
case reference_element::T: return inv_piola_T (x, omega.dis_node(dis_inod [0]),
omega.dis_node(dis_inod [1]),
omega.dis_node(dis_inod [2]),
omega.dis_node(dis_inod [3]));
}
}
// non-linear transformation: q,P,H or high order > 1 => use Newton
inv_piola<T> F;
F.reset (omega, hat_K, dis_inod);
F.set_x (x);
point_basic<T> hat_x = F.initial();
size_t max_iter = 500, n_iter = max_iter;
T tol = std::numeric_limits<Float>::epsilon(), r = tol;
int status = damped_newton (F, hat_x, r, n_iter);
check_macro (status == 0, "inv_piola: newton failed (residue="<<r<<", n_iter="<<n_iter<<")");
return hat_x;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation1(T) \
template \
T \
det_jacobian_piola_transformation ( \
const tensor_basic<T>& DF, \
size_t d, \
size_t map_d); \
template \
tensor_basic<T> \
pseudo_inverse_jacobian_piola_transformation ( \
const tensor_basic<T>& DF, \
size_t d, \
size_t map_d); \
template \
T \
weight_coordinate_system ( \
space_constant::coordinate_type sys_coord, \
const point_basic<T>& xq); \
template \
void \
map_projector ( \
const tensor_basic<T>& DF, \
size_t d, \
size_t map_d, \
tensor_basic<T>& P); \
#define _RHEOLEF_instanciation2(T,M) \
template \
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); \
template \
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); \
template \
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); \
template \
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 \
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); \
template \
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); \
_RHEOLEF_instanciation1(Float)
_RHEOLEF_instanciation2(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation2(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_basic
Definition: basis.h:532
rheolef::basis_basic::grad_evaluate
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
Definition: basis.h:875
rheolef::basis_on_pointset
Definition: basis_on_pointset.h:192
rheolef::basis_on_pointset::nnod
size_type nnod(reference_element hat_K) const
Definition: basis_on_pointset.h:263
rheolef::basis_on_pointset::get_quadrature
const quadrature< T > & get_quadrature() const
Definition: basis_on_pointset.h:277
rheolef::geo_abstract_base_rep
abstract base interface class
Definition: geo.h:248
rheolef::geo_basic< T, M >
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::dimension
size_type dimension() const
Definition: geo_element.h:167
rheolef::geo_element::master
size_type master(bool i) const
Definition: geo_element.h:165
rheolef::geo_element::dis_ie
size_type dis_ie() const
Definition: geo_element.h:163
rheolef::inv_piola
Definition: inv_piola.h:33
rheolef::inv_piola::set_x
void set_x(const value_type &x1)
Definition: inv_piola.h:41
rheolef::inv_piola::reset
void reset(const geo_basic< T, M > &omega, const reference_element &hat_K, const std::vector< size_t > &dis_inod)
Definition: inv_piola.h:77
rheolef::inv_piola::initial
value_type initial() const
Definition: inv_piola.h:89
rheolef::piola_on_pointset
Definition: piola_on_pointset.h:142
rheolef::point_basic
Definition: point.h:87
rheolef::quadrature
Definition: quadrature.h:186
rheolef::quadrature::const_iterator
rep::const_iterator const_iterator
Definition: quadrature.h:195
rheolef::quadrature::end
const_iterator end(reference_element hat_K) const
Definition: quadrature.h:219
rheolef::quadrature::begin
const_iterator begin(reference_element hat_K) const
Definition: quadrature.h:218
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
rheolef::reference_element::e
static const variant_type e
Definition: reference_element.h:76
rheolef::reference_element::dimension
size_type dimension() const
Definition: reference_element.h:101
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::tensor_basic
Definition: tensor.h:90
rheolef::tensor_basic::col
point_basic< T > col(size_type i) const
Definition: tensor.cc:323
rheolef::tensor_basic::determinant
T determinant(size_type d=3) const
Definition: tensor.cc:288
rheolef::tensor_basic::fill
void fill(const T &init_val)
Definition: tensor.h:252
rheolef::tensor_basic::set_row
void set_row(const point_basic< T > &r, size_t i, size_t d=3)
Definition: tensor.h:435
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
fatal_macro
#define fatal_macro(message)
Definition: dis_macros.h:33
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)")
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
mkgeo_ball.b
int b
Definition: mkgeo_ball.sh:152
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
mkgeo_contraction.status
status
Definition: mkgeo_contraction.sh:290
mkgeo_grid.sys_coord
string sys_coord
Definition: mkgeo_grid.sh:171
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef::space_constant::coordinate_system_name
std::string coordinate_system_name(coordinate_type i)
Definition: space_constant.cc:69
rheolef::space_constant::coordinate_type
coordinate_type
Definition: space_constant.h:121
rheolef::space_constant::cartesian
@ cartesian
Definition: space_constant.h:122
rheolef::space_constant::axisymmetric_zr
@ axisymmetric_zr
Definition: space_constant.h:124
rheolef::space_constant::axisymmetric_rz
@ axisymmetric_rz
Definition: space_constant.h:123
rheolef::space_numbering::dis_inod
void dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
Definition: space_numbering.cc:177
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::inv
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition: tensor.cc:219
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::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
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::vect
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition: point.h:265
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::norm2
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition: vec.h:379
rheolef::inverse_piola_transformation
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)
Definition: piola_util.cc:459
rheolef::invert_3x3
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition: tensor.cc:333
rheolef::_RHEOLEF_instanciation1
_RHEOLEF_instanciation1(Float) _RHEOLEF_instanciation2(Float
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::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::weight_coordinate_system
T weight_coordinate_system(space_constant::coordinate_type sys_coord, const point_basic< T > &xq)
Definition: piola_util.cc:307
rheolef::piola_transformation_and_weight_integration
void piola_transformation_and_weight_integration(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_quad, 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)
Definition: piola_util.cc:326
rheolef::normal_from_piola_transformation
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)
Definition: piola_util.cc:164
rheolef::piola_transformation
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)
Definition: piola_util.cc:43
rheolef::damped_newton
int damped_newton(const Problem &P, const Preconditioner &T, Field &u, Real &tol, Size &max_iter, odiststream *p_derr=0)
see the damped_newton page for the full documentation
Definition: damped-newton-generic.h:29
_RHEOLEF_instanciation2
#define _RHEOLEF_instanciation2(T, M)
Definition: piola_util.cc:520
a
Definition: diffusion_isotropic.h:25
rheolef::distributed
Definition: distributed.h:88
rheolef::sequential
Definition: distributed.h:28
epsilon
Float epsilon
Definition: transmission_error.cc:25