Loading...
Searching...
No Matches
warburton.icc
Go to the documentation of this file.
#ifndef _RHEOLEF_WARBURTON_ICC
#define _RHEOLEF_WARBURTON_ICC
//
// compute the Warburton warp and blend point set (see War-2006 HesWar-2008)
//
// author: Pierre.Saramito@imag.fr
//
// date: 2 september 2017
//
#include "rheolef/point.h"
#include "rheolef/tensor.h"
#include "rheolef/reference_element.h"
#include "rheolef/gauss_lobatto_jacobi.h"
#include "reference_element_aux.icc"
#include "eigen_util.h"
namespace rheolef {
// ---------------------------------------------------------------------------
// 1d : edge
// ---------------------------------------------------------------------------
// Compute nodes in [-1:1] edge for polynomial of order "degree"
template<class T>
void
warburton_e (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool map_on_reference_element=true)
{
// Compute Gauss-Lobatto Legendre node distribution
std::vector<T> zeta(degree+1), omega(degree+1);
gauss_lobatto_jacobi (degree+1, T(0), T(0), zeta.begin(), omega.begin());
x.resize (degree+1);
for (size_t i = 0; i < degree+1; i++) {
x[i] = point_basic<T>(zeta[i],0,0);
}
if (!map_on_reference_element) return;
// permutation: by decreasing subgeo dimension, for rheolef
// transform: to reference_element
std::vector<size_t> perm(degree+1);
for (size_t i = 0; i < degree+1; i++) {
perm[i] = reference_element_e::ilat2loc_inod (degree, point_basic<size_t>(i));
}
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x_new (x.size());
for (size_t i = 0; i < degree+1; i++) {
x_new[perm[i]][0] = 0.5*(1+x[i][0]);
}
x = x_new;
}
// ---------------------------------------------------------------------------
// 2d : quadrangle
// ---------------------------------------------------------------------------
// Compute nodes in [-1:1]^2 quadrangle for polynomial of order "degree"
template<class T>
void
warburton_q (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool dummy=true)
{
// Compute Gauss-Lobatto Legendre node distribution
std::vector<T> zeta(degree+1), omega(degree+1);
gauss_lobatto_jacobi (degree+1, T(0), T(0), zeta.begin(), omega.begin());
x.resize ((degree+1)*(degree+1));
for (size_t i = 0; i < degree+1; i++) {
for (size_t j = 0; j < degree+1; j++) {
size_t loc_idof = reference_element_q::ilat2loc_inod (degree, point_basic<size_t>(i,j));
x[loc_idof] = point_basic<T>(zeta[i], zeta[j],0);
}}
}
// ---------------------------------------------------------------------------
// 3d : hexahedron
// ---------------------------------------------------------------------------
// Compute nodes in [-1:1]^3 quadrangle for polynomial of order "degree"
template<class T>
void
warburton_H (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool dummy=true)
{
// Compute Gauss-Lobatto Legendre node distribution
std::vector<T> zeta(degree+1), omega(degree+1);
gauss_lobatto_jacobi (degree+1, T(0), T(0), zeta.begin(), omega.begin());
x.resize ((degree+1)*(degree+1)*(degree+1));
for (size_t i = 0; i < degree+1; i++) {
for (size_t j = 0; j < degree+1; j++) {
for (size_t k = 0; k < degree+1; k++) {
size_t loc_idof = reference_element_H::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
x[loc_idof] = point_basic<T>(zeta[i], zeta[j], zeta[k]);
}}}
}
// ---------------------------------------------------------------------------
// 2d : triangle
// ---------------------------------------------------------------------------
// initialize the transpose of the 1d Vandermonde matrix: v_ij = p_j(x_i);
template<class T, class Matrix>
static
void
trans_vandermonde1d (size_t N, const std::vector<T>& x, Matrix& trans_v)
{
for (size_t j = 0; j < N+1; j++) {
jacobi<T> p_j (j, 0, 0);
for (size_t i = 0; i < N+1; i++) {
trans_v(j,i) = p_j(x[i]);
}
}
}
// Compute scaled warp function at order N based on rout interpolation nodes
template<class T>
static
void
warpfactor (size_t N, const std::vector<T>& rout, std::vector<T>& warp)
{
// Compute equidistant node distribution
std::vector<T> req(N+1);
req[0] = -1;
for (size_t i = 1; i < N; i++) {
req[i] = -1 + 2*T(int(i))/T(int(N));
}
req[N] = 1;
// Compute Gauss-Lobatto Legendre node distribution
std::vector<T> zeta(N+1), omega(N+1);
gauss_lobatto_jacobi (N+1, T(0), T(0), zeta.begin(), omega.begin());
// Compute V based on req
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> trans_v_eq (N+1, N+1);
trans_vandermonde1d (N, req, trans_v_eq);
// Evaluate Lagrange polynomial at rout: pmat(i,j) = p_i(rout[j]);
size_t Nr = rout.size();
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> pmat (N+1, Nr);
for (size_t i = 0; i < N+1; i++) {
jacobi<T> p_i (i, 0, 0);
for (size_t j = 0; j < Nr; j++) {
pmat(i,j) = p_i(rout[j]);
}
}
// pmat := inv(v_eq^T)*pmat
Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> inv_trans_v_eq (N+1, N+1);
check_macro (invert(trans_v_eq, inv_trans_v_eq), "matrix inversion failed");
pmat = inv_trans_v_eq*pmat;
// Compute warp factor: warp = pmat^T*(LGLr - req);
warp.resize (Nr);
for (size_t i = 0; i < Nr; i++) {
warp[i] = T(0);
for (size_t j = 0; j < N+1; j++) {
warp[i] += pmat(j,i)*(zeta[j] - req[j]);
}
}
// Scale factor
const T tol = 1e-10;
for (size_t i = 0, n = warp.size(); i < n; i++) {
T zerof = (fabs(rout[i]) < 1.0-tol);
T scale = 1.0 - sqr(zerof*rout[i]);
warp[i] = warp[i]/scale + warp[i]*(zerof-1);
}
}
// Compute nodes in equilateral triangle for polynomial of order "degree"
// note: equilateral triangle bbox = [-1:1]*[-1/sqrt(3):2/sqrt(3)]
template<class T>
void
warburton_t (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool map_on_reference_element=true)
{
// Set optimized parameter, alpha, depending on order "degree"
static const T alpha_opt[] = {0.0000, 0.0000, 1.4152, 0.1001, 0.2751,
0.9800, 1.0999, 1.2832, 1.3648, 1.4773,
1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
T alpha = (degree <= sizeof(alpha_opt)/sizeof(T)) ? alpha_opt[degree-1] : T(5)/T(3);
// total number of nodes
size_t loc_ndof = (degree+1)*(degree+2)/2;
x.resize(loc_ndof);
// Create equidistributed nodes on equilateral triangle
std::vector<T> L1 (loc_ndof), L2 (loc_ndof), L3 (loc_ndof);
const T sqrt_3 = sqrt(T(3));
for (size_t i = 0, loc_idof = 0; i < degree+1; i++) {
for (size_t j = 0; j < degree+1-i; j++, loc_idof++) {
L1[loc_idof] = T(int(i))/T(int(degree));
L3[loc_idof] = T(int(j))/T(int(degree));
L2[loc_idof] = 1 - L1[loc_idof] - L3[loc_idof];
x[loc_idof][0] = - L2[loc_idof] + L3[loc_idof];
x[loc_idof][1] = (2*L1[loc_idof] - L2[loc_idof] - L3[loc_idof])/sqrt_3;
}
}
// Compute blending function at each node for each edge
std::vector<T> blend1 (loc_ndof), blend2 (loc_ndof), blend3 (loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
blend1[loc_idof] = 4*L2[loc_idof]*L3[loc_idof];
blend2[loc_idof] = 4*L1[loc_idof]*L3[loc_idof];
blend3[loc_idof] = 4*L1[loc_idof]*L2[loc_idof];
}
// Amount of warp for each node, for each edge
std::vector<T> warpf1 (loc_ndof, 0), warpf2 (loc_ndof, 0), warpf3 (loc_ndof, 0);
std::vector<T> L32 (loc_ndof), L13 (loc_ndof), L21 (loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
L32[loc_idof] = L3[loc_idof] - L2[loc_idof];
L13[loc_idof] = L1[loc_idof] - L3[loc_idof];
L21[loc_idof] = L2[loc_idof] - L1[loc_idof];
}
warpfactor (degree, L32, warpf1);
warpfactor (degree, L13, warpf2);
warpfactor (degree, L21, warpf3);
// Combine blend & warp
std::vector<T> warp1 (loc_ndof), warp2 (loc_ndof), warp3 (loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
warp1[loc_idof] = blend1[loc_idof]*warpf1[loc_idof]*(1 + sqr(alpha*L1[loc_idof]));
warp2[loc_idof] = blend2[loc_idof]*warpf2[loc_idof]*(1 + sqr(alpha*L2[loc_idof]));
warp3[loc_idof] = blend3[loc_idof]*warpf3[loc_idof]*(1 + sqr(alpha*L3[loc_idof]));
}
// Accumulate deformations associated with each edge
const T pi = acos(T(-1));
const T cos23 = cos(2*pi/3),
sin23 = sin(2*pi/3),
cos43 = cos(4*pi/3),
sin43 = sin(4*pi/3);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
x[loc_idof][0] += warp1[loc_idof] + cos23*warp2[loc_idof] + cos43*warp3[loc_idof];
x[loc_idof][1] += sin23*warp2[loc_idof] + sin43*warp3[loc_idof];
}
if (!map_on_reference_element) return;
// transform: to reference_element, from the isocele triangle used by Warburton
// permutation: by decreasing subgeo dimension, for rheolef
T xa = -1; T ya = -sqrt(T(3))/3.;
T xb = 1; T yb = -sqrt(T(3))/3.;
T xc = 0; T yc = 2*sqrt(T(3))/3.;
// introduce the Piola transformation from ref elem to isocele
T a11 = xb-xa;
T a12 = xc-xa;
T a21 = yb-ya;
T a22 = yc-ya;
// invert the Piola transformation
T det = a11*a22-a12*a21;
T b11 = a22/det;
T b12 = -a12/det;
T b21 = -a21/det;
T b22 = a11/det;
std::vector<size_t> perm(loc_ndof);
for (size_t loc_idof = 0, j = 0; j <= degree; j++) {
for (size_t i = 0; i+j <= degree; i++) {
perm[loc_idof++] = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
}
}
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x_new (x.size());
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
T x0 = x[loc_idof][0];
T y0 = x[loc_idof][1];
size_t new_loc_idof = perm[loc_idof];
x_new[new_loc_idof][0] = b11*(x0-xa) + b12*(y0-ya);
x_new[new_loc_idof][1] = b21*(x0-xa) + b22*(y0-ya);
}
x = x_new;
}
// ---------------------------------------------------------------------------
// 3d : prism
// ---------------------------------------------------------------------------
// compute the equidistributed nodes on the reference prism
// note: reference prism bbox = ref_triangle*[-1:1]
template<class T>
void
warburton_P (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool map_on_reference_element=true)
{
// for x,y: reference triangle, but unpermuted nodes
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> xt;
warburton_t (degree, xt);
// for z: compute Gauss-Lobatto Legendre node distribution
std::vector<T> zeta(degree+1), omega(degree+1);
gauss_lobatto_jacobi (degree+1, T(0), T(0), zeta.begin(), omega.begin());
// tensorial product: ((x,y),z)
x.resize (xt.size()*zeta.size());
for (size_t j = 0; j <= degree; j++) {
for (size_t i = 0; i+j <= degree; i++) {
size_t loc_t_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
for (size_t k = 0; k <= degree; k++) {
size_t loc_idof = reference_element_P::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
x[loc_idof] = point_basic<T> (xt[loc_t_idof][0], xt[loc_t_idof][1], zeta[k]);
}
}
}
}
// ---------------------------------------------------------------------------
// 3d : tetrahedron
// ---------------------------------------------------------------------------
// compute the equidistributed nodes on the reference tetrahedron
// note: reference tetrahedron bbox = [-1:1]^3
template<class T>
void
equispaced_T (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x)
{
// total number of nodes
size_t loc_ndof = (degree+1)*(degree+2)*(degree+3)/6;
x.resize (loc_ndof);
// create equidistributed nodes on equilateral triangle
size_t loc_idof = 0;
for (size_t k = 0; k < degree+1; k++)
for (size_t j = 0; j+k < degree+1; j++)
for (size_t i = 0; i+j+k < degree+1; i++, loc_idof++) {
x[loc_idof] = point_basic<T>(
-1 + 2*T(int(i))/T(int(degree)),
-1 + 2*T(int(j))/T(int(degree)),
-1 + 2*T(int(k))/T(int(degree)));
}
}
// compute one-dimensional edge warping function
template<class T>
void
evalwarp (
size_t degree,
const std::vector<T>& xnodes,
const std::vector<T>& xout,
std::vector<T>& warp)
{
std::vector<T> xeq (degree+1);
for (size_t i = 0; i < degree+1; i++) {
xeq[i] = -1 + 2*T(int(degree-i))/T(int(degree));
}
size_t loc_ndof = xout.size();
warp.resize (loc_ndof);
std::vector<T> d (loc_ndof);
for (size_t i = 0; i < degree+1; i++) {
std::fill (d.begin(), d.end(), xnodes[i] - xeq[i]);
for (size_t j = 1; j+1 < degree+1; j++) {
if (i != j) {
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
d[loc_idof] *= (xout[loc_idof]-xeq[j])/(xeq[i]-xeq[j]);
}
}
}
if (i != 0) {
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
d[loc_idof] = - d[loc_idof]/(xeq[i]-xeq[0]);
}
}
if (i != degree) {
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
d[loc_idof] = d[loc_idof]/(xeq[i]-xeq[degree]);
}
}
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
warp[loc_idof] += d[loc_idof];
}
}
}
// compute two-dimensional Warp & Blend transform
template<class T>
void
evalshift (
size_t degree,
const T& alpha,
const std::vector<T>& L0,
const std::vector<T>& L1,
const std::vector<T>& L2,
std::vector<T>& dx,
std::vector<T>& dy)
{
// 1) compute Gauss-Lobatto-Legendre node distribution
std::vector<T> gauss_x(degree+1), omega(degree+1);
gauss_lobatto_jacobi (degree+1, T(0), T(0), gauss_x.begin(), omega.begin());
for (size_t i = 0; i < degree+1; i++) {
gauss_x[i] = - gauss_x[i];
}
// 2) compute blending function at each node for each edge
size_t loc_ndof = L0.size();
std::vector<T> blend1 (loc_ndof), blend2 (loc_ndof), blend3 (loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
blend1[loc_idof] = L1[loc_idof]*L2[loc_idof];
blend2[loc_idof] = L0[loc_idof]*L2[loc_idof];
blend3[loc_idof] = L0[loc_idof]*L1[loc_idof];
}
// 3) amount of warp for each node, for each edge
std::vector<T> L21 (loc_ndof), L02 (loc_ndof), L10 (loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
L21[loc_idof] = L2[loc_idof] - L1[loc_idof];
L02[loc_idof] = L0[loc_idof] - L2[loc_idof];
L10[loc_idof] = L1[loc_idof] - L0[loc_idof];
}
std::vector<T> warpfactor1 (loc_ndof), warpfactor2 (loc_ndof), warpfactor3 (loc_ndof);
evalwarp (degree, gauss_x, L21, warpfactor1);
evalwarp (degree, gauss_x, L02, warpfactor2);
evalwarp (degree, gauss_x, L10, warpfactor3);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
warpfactor1 [loc_idof] *= 4;
warpfactor2 [loc_idof] *= 4;
warpfactor3 [loc_idof] *= 4;
}
// 4) combine blend & warp
std::vector<T> warp1 (loc_ndof), warp2 (loc_ndof), warp3 (loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
warp1 [loc_idof] = blend1 [loc_idof]*warpfactor1 [loc_idof]*(1 + sqr(alpha*L0 [loc_idof]));
warp2 [loc_idof] = blend2 [loc_idof]*warpfactor2 [loc_idof]*(1 + sqr(alpha*L1 [loc_idof]));
warp3 [loc_idof] = blend3 [loc_idof]*warpfactor3 [loc_idof]*(1 + sqr(alpha*L2 [loc_idof]));
}
// 5) evaluate shift in equilateral triangle
dx.resize (loc_ndof);
dy.resize (loc_ndof);
const T pi = acos(T(-1));
const T cos23 = cos(2*pi/3),
sin23 = sin(2*pi/3),
cos43 = cos(4*pi/3),
sin43 = sin(4*pi/3);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
dx[loc_idof] = warp1[loc_idof] + cos23*warp2[loc_idof] + cos43*warp3[loc_idof];
dy[loc_idof] = sin23*warp2[loc_idof] + sin43*warp3[loc_idof];
// les indices 6 & 7 differents de matlab:
}
}
// compute two-dimensional Warp & Blend transform
template<class T>
inline
void
WarpShiftFace3D (
size_t degree,
const T& alpha,
const T& dummy,
const std::vector<T>& Ldummy,
const std::vector<T>& L1,
const std::vector<T>& L2,
const std::vector<T>& L3,
std::vector<T>& warpx,
std::vector<T>& warpy)
{
evalshift (degree, alpha, L1, L2, L3, warpx, warpy);
}
// compute warp & blend nodes
template<class T>
void
warburton_T (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool map_on_reference_element=true)
{
// choose optimized blending parameter
static const T alpha_opt[] = {0, 0, 0, 0.1002, 1.13320,
1.5608, 1.3413, 1.2577, 1.1603, 1.10153,
0.6080, 0.4523, 0.8856, 0.8717, 0.96550};
T alpha = (degree <= sizeof(alpha_opt)/sizeof(T)) ? alpha_opt[degree-1] : T(5)/T(3);
// create equidistributed nodes in the reference tetrahedron
equispaced_T (degree, x);
size_t loc_ndof = x.size();
// compute barycentric coordinates
std::array<std::vector<T>, 4> L;
for (size_t i = 0; i < 4; i++) L[i].resize (loc_ndof);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
L[0][loc_idof] = (1 + x[loc_idof][2])/2;
L[1][loc_idof] = (1 + x[loc_idof][1])/2;
L[2][loc_idof] = - (1 + x[loc_idof][0] + x[loc_idof][1] + x[loc_idof][2])/2;
L[3][loc_idof] = (1 + x[loc_idof][0])/2;
}
// set vertices of the equilateral tetrahedron
point_basic<T> v1 (-1, -1/sqrt(T(3)), -1/sqrt(T(6))),
v2 ( 1, -1/sqrt(T(3)), -1/sqrt(T(6))),
v3 ( 0, 2/sqrt(T(3)), -1/sqrt(T(6))),
v4 ( 0, 0, 3/sqrt(T(6)));
// create equidistributed nodes in the equilateral tetrahedron
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
x[loc_idof] = L[2][loc_idof]*v1 + L[3][loc_idof]*v2 + L[1][loc_idof]*v3 + L[0][loc_idof]*v4;
}
// orthogonal axis tangents on faces 1-4
std::array<point_basic<T>, 4> t1, t2;
t1[0] = v2-v1;
t1[1] = v2-v1;
t1[2] = v3-v2;
t1[3] = v3-v1;
t2[0] = v3-0.5*(v1+v2);
t2[1] = v4-0.5*(v1+v2);
t2[2] = v4-0.5*(v2+v3);
t2[3] = v4-0.5*(v1+v3);
// normalize tangents
for (size_t i = 0; i < 4; i++) {
t1[i] /= norm(t1[i]);
t2[i] /= norm(t2[i]);
}
// Warp and blend for each face (accumulated in shift)
const T tol = 1e-10;
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> shift (loc_ndof);
std::vector<T> blend (loc_ndof),
denom (loc_ndof);
for (size_t ifac = 0; ifac < 4; ifac++) {
size_t a, b, c, d;
switch (ifac) {
case 0: a = 0; b = 1; c = 2; d = 3; break;
case 1: a = 1; b = 0; c = 2; d = 3; break;
case 2: a = 2; b = 0; c = 3; d = 1; break;
case 3:
default: a = 3; b = 0; c = 2; d = 1; break;
}
// compute warp tangential to face
std::vector<T> warp1 (loc_ndof), warp2 (loc_ndof);
WarpShiftFace3D (degree, alpha, alpha, L[a], L[b], L[c], L[d], warp1, warp2);
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
// compute volume blending
blend[loc_idof] = L[b][loc_idof]*L[c][loc_idof]*L[d][loc_idof];
// modify linear blend
denom[loc_idof] = (L[b][loc_idof] + 0.5*L[a][loc_idof])
*(L[c][loc_idof] + 0.5*L[a][loc_idof])
*(L[d][loc_idof] + 0.5*L[a][loc_idof]);
if (denom[loc_idof] > tol) {
blend[loc_idof] = (1 + sqr(alpha*L[a][loc_idof]))*blend[loc_idof]/denom[loc_idof];
}
}
// compute warp & blend
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
shift[loc_idof] += blend[loc_idof]*warp1 [loc_idof]*t1 [ifac]
+ blend[loc_idof]*warp2 [loc_idof]*t2 [ifac];
}
// fix face warp
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
size_t n_sup = 0;
if (L[b][loc_idof] > tol) n_sup++;
if (L[c][loc_idof] > tol) n_sup++;
if (L[d][loc_idof] > tol) n_sup++;
if ((L[a][loc_idof] < tol) && (n_sup < 3)) {
shift[loc_idof] = warp1[loc_idof]*t1[ifac] + warp2[loc_idof]*t2[ifac];
}
}
}
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
x[loc_idof] += shift[loc_idof];
}
if (!map_on_reference_element) return;
// define the isocele tetra used by Warburton
T side_length = 2;
T height = sqrt(T(6))/3*side_length;
point_basic<T> a (-1, -1/sqrt(T(3)), -height/4),
b ( 1, -1/sqrt(T(3)), -height/4),
c ( 0, 2/sqrt(T(3)), -height/4),
d ( 0, 0, 3*height/4);
// define and invert the Piola transform
tensor_basic<T> A, inv_A;
A.set_column (b-a, 0);
A.set_column (c-a, 1);
A.set_column (d-a, 2);
invert_3x3 (A,inv_A);
std::vector<size_t> perm(loc_ndof);
for (size_t loc_idof = 0, k = 0; k <= degree; k++) {
for (size_t j = 0; j+k <= degree; j++) {
for (size_t i = 0; i+j+k <= degree; i++) {
perm[loc_idof++] = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
}
}
}
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x_new (x.size());
for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
size_t new_loc_idof = perm[loc_idof];
x_new[new_loc_idof] = inv_A*(x[loc_idof]-a);
}
x = x_new;
}
// ---------------------------------------------------------------------------
// all cases
// ---------------------------------------------------------------------------
template<class T>
void
pointset_lagrange_warburton (
reference_element hat_K,
size_t degree,
Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_xnod,
bool map_on_reference_element = true)
{
hat_xnod.resize (reference_element::n_node(hat_K.variant(), degree));
if (degree == 0) {
reference_element_barycenter (hat_K, hat_xnod[0]);
return;
}
switch (hat_K.variant()) {
case reference_element::p:
hat_xnod [0] = point_basic<T> (T(0)); break;
case reference_element::e:
warburton_e (degree, hat_xnod, map_on_reference_element); break;
case reference_element::q:
warburton_q (degree, hat_xnod, map_on_reference_element); break;
case reference_element::H:
warburton_H (degree, hat_xnod, map_on_reference_element); break;
case reference_element::t:
warburton_t (degree, hat_xnod, map_on_reference_element); break;
case reference_element::P:
warburton_P (degree, hat_xnod, map_on_reference_element); break;
case reference_element::T:
warburton_T (degree, hat_xnod, map_on_reference_element); break;
default: error_macro ("unexpected element type `"<<hat_K.name()<<"'");
}
}
} // namespace rheolef
#endif // _RHEOLEF_WARBURTON_ICC
rheolef::point_basic
Definition point.h:86
rheolef::reference_element_H::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:790
rheolef::reference_element_P::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:695
rheolef::reference_element_T::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:582
rheolef::reference_element_e::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:350
rheolef::reference_element_q::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:481
rheolef::reference_element_t::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition reference_element.cc:413
rheolef::reference_element
see the reference_element page for the full documentation
Definition reference_element.h:66
rheolef::reference_element::H
static const variant_type H
Definition reference_element.h:81
rheolef::reference_element::q
static const variant_type q
Definition reference_element.h:78
rheolef::reference_element::e
static const variant_type e
Definition reference_element.h:76
rheolef::reference_element::p
static const variant_type p
Definition reference_element.h:75
rheolef::reference_element::variant
variant_type variant() const
Definition reference_element.h:99
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition reference_element.cc:201
rheolef::reference_element::T
static const variant_type T
Definition reference_element.h:79
rheolef::reference_element::P
static const variant_type P
Definition reference_element.h:80
rheolef::reference_element::name
char name() const
Definition reference_element.h:100
rheolef::reference_element::t
static const variant_type t
Definition reference_element.h:77
rheolef::tensor_basic
Definition tensor.h:90
error_macro
#define error_macro(message)
Definition dis_macros.h:49
eigen_util.h
T
Expr1::float_type T
Definition field_expr.h:230
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
reference_element_aux.icc
mkgeo_ball.n
int n
Definition mkgeo_ball.sh:150
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::warburton_T
void warburton_T(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool map_on_reference_element=true)
Definition warburton.icc:443
rheolef::reference_element_barycenter
void reference_element_barycenter(reference_element hat_K, point_basic< T > &c)
Definition reference_element_aux.icc:38
rheolef::pointset_lagrange_warburton
void pointset_lagrange_warburton(reference_element hat_K, size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &hat_xnod, bool map_on_reference_element=true)
Definition warburton.icc:574
rheolef::evalwarp
void evalwarp(size_t degree, const std::vector< T > &xnodes, const std::vector< T > &xout, std::vector< T > &warp)
Definition warburton.icc:323
rheolef::gauss_lobatto_jacobi
void gauss_lobatto_jacobi(Size R, typename std::iterator_traits< OutputIterator1 >::value_type alpha, typename std::iterator_traits< OutputIterator1 >::value_type beta, OutputIterator1 zeta, OutputIterator2 omega)
Definition gauss_lobatto_jacobi.icc:28
rheolef::WarpShiftFace3D
void WarpShiftFace3D(size_t degree, const T &alpha, const T &dummy, const std::vector< T > &Ldummy, const std::vector< T > &L1, const std::vector< T > &L2, const std::vector< T > &L3, std::vector< T > &warpx, std::vector< T > &warpy)
Definition warburton.icc:427
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition tiny_lu.h:127
rheolef::warburton_H
void warburton_H(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool dummy=true)
Definition warburton.icc:91
rheolef::warburton_t
void warburton_t(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool map_on_reference_element=true)
Definition warburton.icc:176
rheolef::warburton_e
void warburton_e(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool map_on_reference_element=true)
Definition warburton.icc:45
rheolef::evalshift
void evalshift(size_t degree, const T &alpha, const std::vector< T > &L0, const std::vector< T > &L1, const std::vector< T > &L2, std::vector< T > &dx, std::vector< T > &dy)
Definition warburton.icc:363
rheolef::warburton_P
void warburton_P(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool map_on_reference_element=true)
Definition warburton.icc:276
rheolef::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
rheolef::equispaced_T
void equispaced_T(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x)
Definition warburton.icc:303
rheolef::warburton_q
void warburton_q(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool dummy=true)
Definition warburton.icc:73
rheolef::invert_3x3
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition tensor.cc:333
N
size_t N
Definition neumann-laplace-check.cc:24
d
size_t d
Definition neumann-laplace-lambda.cc:28
A
Definition diffusion_isotropic.h:25