1#ifndef _RHEOLEF_WARBURTON_ICC
2#define _RHEOLEF_WARBURTON_ICC
30#include "rheolef/point.h"
31#include "rheolef/tensor.h"
32#include "rheolef/reference_element.h"
33#include "rheolef/gauss_lobatto_jacobi.h"
48 std::vector<T> zeta(degree+1), omega(degree+1);
51 for (
size_t i = 0; i < degree+1; i++) {
54 if (!map_on_reference_element)
return;
57 std::vector<size_t> perm(degree+1);
58 for (
size_t i = 0; i < degree+1; i++) {
61 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x_new (x.size());
62 for (
size_t i = 0; i < degree+1; i++) {
63 x_new[perm[i]][0] = 0.5*(1+x[i][0]);
76 std::vector<T> zeta(degree+1), omega(degree+1);
78 x.resize ((degree+1)*(degree+1));
79 for (
size_t i = 0; i < degree+1; i++) {
80 for (
size_t j = 0; j < degree+1; j++) {
94 std::vector<T> zeta(degree+1), omega(degree+1);
96 x.resize ((degree+1)*(degree+1)*(degree+1));
97 for (
size_t i = 0; i < degree+1; i++) {
98 for (
size_t j = 0; j < degree+1; j++) {
99 for (
size_t k = 0; k < degree+1; k++) {
108template<
class T,
class Matrix>
111trans_vandermonde1d (
size_t N,
const std::vector<T>& x, Matrix& trans_v)
113 for (
size_t j = 0; j <
N+1; j++) {
114 jacobi<T> p_j (j, 0, 0);
115 for (
size_t i = 0; i <
N+1; i++) {
116 trans_v(j,i) = p_j(x[i]);
124warpfactor (
size_t N,
const std::vector<T>& rout, std::vector<T>& warp)
127 std::vector<T> req(
N+1);
129 for (
size_t i = 1; i <
N; i++) {
130 req[i] = -1 + 2*
T(
int(i))/
T(
int(
N));
135 std::vector<T> zeta(
N+1), omega(
N+1);
139 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> trans_v_eq (
N+1,
N+1);
140 trans_vandermonde1d (
N, req, trans_v_eq);
143 size_t Nr = rout.size();
144 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> pmat (
N+1, Nr);
145 for (
size_t i = 0; i <
N+1; i++) {
146 jacobi<T> p_i (i, 0, 0);
147 for (
size_t j = 0; j < Nr; j++) {
148 pmat(i,j) = p_i(rout[j]);
152 Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> inv_trans_v_eq (
N+1,
N+1);
154 pmat = inv_trans_v_eq*pmat;
158 for (
size_t i = 0; i < Nr; i++) {
160 for (
size_t j = 0; j <
N+1; j++) {
161 warp[i] += pmat(j,i)*(zeta[j] - req[j]);
166 for (
size_t i = 0, n = warp.size(); i <
n; i++) {
167 T zerof = (fabs(rout[i]) < 1.0-tol);
168 T scale = 1.0 - sqr(zerof*rout[i]);
169 warp[i] = warp[i]/scale + warp[i]*(zerof-1);
179 static const T alpha_opt[] = {0.0000, 0.0000, 1.4152, 0.1001, 0.2751,
180 0.9800, 1.0999, 1.2832, 1.3648, 1.4773,
181 1.4959, 1.5743, 1.5770, 1.6223, 1.6258};
182 T alpha = (degree <=
sizeof(alpha_opt)/
sizeof(
T)) ? alpha_opt[degree-1] :
T(5)/
T(3);
185 size_t loc_ndof = (degree+1)*(degree+2)/2;
189 std::vector<T> L1 (loc_ndof), L2 (loc_ndof), L3 (loc_ndof);
190 const T sqrt_3 = sqrt(
T(3));
191 for (
size_t i = 0, loc_idof = 0; i < degree+1; i++) {
192 for (
size_t j = 0; j < degree+1-i; j++, loc_idof++) {
193 L1[loc_idof] =
T(
int(i))/
T(
int(degree));
194 L3[loc_idof] =
T(
int(j))/
T(
int(degree));
195 L2[loc_idof] = 1 - L1[loc_idof] - L3[loc_idof];
196 x[loc_idof][0] = - L2[loc_idof] + L3[loc_idof];
197 x[loc_idof][1] = (2*L1[loc_idof] - L2[loc_idof] - L3[loc_idof])/sqrt_3;
201 std::vector<T> blend1 (loc_ndof), blend2 (loc_ndof), blend3 (loc_ndof);
202 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
203 blend1[loc_idof] = 4*L2[loc_idof]*L3[loc_idof];
204 blend2[loc_idof] = 4*L1[loc_idof]*L3[loc_idof];
205 blend3[loc_idof] = 4*L1[loc_idof]*L2[loc_idof];
208 std::vector<T> warpf1 (loc_ndof, 0), warpf2 (loc_ndof, 0), warpf3 (loc_ndof, 0);
209 std::vector<T> L32 (loc_ndof), L13 (loc_ndof), L21 (loc_ndof);
210 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
211 L32[loc_idof] = L3[loc_idof] - L2[loc_idof];
212 L13[loc_idof] = L1[loc_idof] - L3[loc_idof];
213 L21[loc_idof] = L2[loc_idof] - L1[loc_idof];
215 warpfactor (degree, L32, warpf1);
216 warpfactor (degree, L13, warpf2);
217 warpfactor (degree, L21, warpf3);
220 std::vector<T> warp1 (loc_ndof), warp2 (loc_ndof), warp3 (loc_ndof);
221 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
222 warp1[loc_idof] = blend1[loc_idof]*warpf1[loc_idof]*(1 + sqr(alpha*L1[loc_idof]));
223 warp2[loc_idof] = blend2[loc_idof]*warpf2[loc_idof]*(1 + sqr(alpha*L2[loc_idof]));
224 warp3[loc_idof] = blend3[loc_idof]*warpf3[loc_idof]*(1 + sqr(alpha*L3[loc_idof]));
227 const T pi = acos(
T(-1));
228 const T cos23 = cos(2*pi/3),
232 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
233 x[loc_idof][0] += warp1[loc_idof] + cos23*warp2[loc_idof] + cos43*warp3[loc_idof];
234 x[loc_idof][1] += sin23*warp2[loc_idof] + sin43*warp3[loc_idof];
236 if (!map_on_reference_element)
return;
239 T xa = -1;
T ya = -sqrt(
T(3))/3.;
240 T xb = 1;
T yb = -sqrt(
T(3))/3.;
241 T xc = 0;
T yc = 2*sqrt(
T(3))/3.;
248 T det = a11*a22-a12*a21;
253 std::vector<size_t> perm(loc_ndof);
254 for (
size_t loc_idof = 0, j = 0; j <= degree; j++) {
255 for (
size_t i = 0; i+j <= degree; i++) {
259 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x_new (x.size());
260 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
261 T x0 = x[loc_idof][0];
262 T y0 = x[loc_idof][1];
263 size_t new_loc_idof = perm[loc_idof];
264 x_new[new_loc_idof][0] = b11*(x0-xa) + b12*(y0-ya);
265 x_new[new_loc_idof][1] = b21*(x0-xa) + b22*(y0-ya);
279 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> xt;
282 std::vector<T> zeta(degree+1), omega(degree+1);
285 x.resize (xt.size()*zeta.size());
286 for (
size_t j = 0; j <= degree; j++) {
287 for (
size_t i = 0; i+j <= degree; i++) {
289 for (
size_t k = 0; k <= degree; k++) {
291 x[loc_idof] =
point_basic<T> (xt[loc_t_idof][0], xt[loc_t_idof][1], zeta[k]);
306 size_t loc_ndof = (degree+1)*(degree+2)*(degree+3)/6;
311 for (
size_t k = 0; k < degree+1; k++)
312 for (
size_t j = 0; j+k < degree+1; j++)
313 for (
size_t i = 0; i+j+k < degree+1; i++, loc_idof++) {
315 -1 + 2*
T(
int(i))/
T(
int(degree)),
316 -1 + 2*
T(
int(j))/
T(
int(degree)),
317 -1 + 2*
T(
int(k))/
T(
int(degree)));
325 const std::vector<T>& xnodes,
326 const std::vector<T>& xout,
327 std::vector<T>& warp)
329 std::vector<T> xeq (degree+1);
330 for (
size_t i = 0; i < degree+1; i++) {
331 xeq[i] = -1 + 2*
T(
int(degree-i))/
T(
int(degree));
333 size_t loc_ndof = xout.size();
334 warp.resize (loc_ndof);
335 std::vector<T>
d (loc_ndof);
336 for (
size_t i = 0; i < degree+1; i++) {
337 std::fill (
d.begin(),
d.end(), xnodes[i] - xeq[i]);
338 for (
size_t j = 1; j+1 < degree+1; j++) {
340 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
341 d[loc_idof] *= (xout[loc_idof]-xeq[j])/(xeq[i]-xeq[j]);
346 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
347 d[loc_idof] = -
d[loc_idof]/(xeq[i]-xeq[0]);
351 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
352 d[loc_idof] =
d[loc_idof]/(xeq[i]-xeq[degree]);
355 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
356 warp[loc_idof] +=
d[loc_idof];
366 const std::vector<T>& L0,
367 const std::vector<T>& L1,
368 const std::vector<T>& L2,
373 std::vector<T> gauss_x(degree+1), omega(degree+1);
375 for (
size_t i = 0; i < degree+1; i++) {
376 gauss_x[i] = - gauss_x[i];
379 size_t loc_ndof = L0.size();
380 std::vector<T> blend1 (loc_ndof), blend2 (loc_ndof), blend3 (loc_ndof);
381 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
382 blend1[loc_idof] = L1[loc_idof]*L2[loc_idof];
383 blend2[loc_idof] = L0[loc_idof]*L2[loc_idof];
384 blend3[loc_idof] = L0[loc_idof]*L1[loc_idof];
387 std::vector<T> L21 (loc_ndof), L02 (loc_ndof), L10 (loc_ndof);
388 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
389 L21[loc_idof] = L2[loc_idof] - L1[loc_idof];
390 L02[loc_idof] = L0[loc_idof] - L2[loc_idof];
391 L10[loc_idof] = L1[loc_idof] - L0[loc_idof];
393 std::vector<T> warpfactor1 (loc_ndof), warpfactor2 (loc_ndof), warpfactor3 (loc_ndof);
394 evalwarp (degree, gauss_x, L21, warpfactor1);
395 evalwarp (degree, gauss_x, L02, warpfactor2);
396 evalwarp (degree, gauss_x, L10, warpfactor3);
397 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
398 warpfactor1 [loc_idof] *= 4;
399 warpfactor2 [loc_idof] *= 4;
400 warpfactor3 [loc_idof] *= 4;
403 std::vector<T> warp1 (loc_ndof), warp2 (loc_ndof), warp3 (loc_ndof);
404 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
405 warp1 [loc_idof] = blend1 [loc_idof]*warpfactor1 [loc_idof]*(1 + sqr(alpha*L0 [loc_idof]));
406 warp2 [loc_idof] = blend2 [loc_idof]*warpfactor2 [loc_idof]*(1 + sqr(alpha*L1 [loc_idof]));
407 warp3 [loc_idof] = blend3 [loc_idof]*warpfactor3 [loc_idof]*(1 + sqr(alpha*L2 [loc_idof]));
410 dx.resize (loc_ndof);
411 dy.resize (loc_ndof);
412 const T pi = acos(
T(-1));
413 const T cos23 = cos(2*pi/3),
417 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
418 dx[loc_idof] = warp1[loc_idof] + cos23*warp2[loc_idof] + cos43*warp3[loc_idof];
419 dy[loc_idof] = sin23*warp2[loc_idof] + sin43*warp3[loc_idof];
431 const std::vector<T>& Ldummy,
432 const std::vector<T>& L1,
433 const std::vector<T>& L2,
434 const std::vector<T>& L3,
435 std::vector<T>& warpx,
436 std::vector<T>& warpy)
438 evalshift (degree, alpha, L1, L2, L3, warpx, warpy);
446 static const T alpha_opt[] = {0, 0, 0, 0.1002, 1.13320,
447 1.5608, 1.3413, 1.2577, 1.1603, 1.10153,
448 0.6080, 0.4523, 0.8856, 0.8717, 0.96550};
449 T alpha = (degree <=
sizeof(alpha_opt)/
sizeof(
T)) ? alpha_opt[degree-1] :
T(5)/
T(3);
453 size_t loc_ndof = x.size();
456 std::array<std::vector<T>, 4> L;
457 for (
size_t i = 0; i < 4; i++) L[i].resize (loc_ndof);
459 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
460 L[0][loc_idof] = (1 + x[loc_idof][2])/2;
461 L[1][loc_idof] = (1 + x[loc_idof][1])/2;
462 L[2][loc_idof] = - (1 + x[loc_idof][0] + x[loc_idof][1] + x[loc_idof][2])/2;
463 L[3][loc_idof] = (1 + x[loc_idof][0])/2;
467 v2 ( 1, -1/sqrt(
T(3)), -1/sqrt(
T(6))),
468 v3 ( 0, 2/sqrt(
T(3)), -1/sqrt(
T(6))),
469 v4 ( 0, 0, 3/sqrt(
T(6)));
472 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
473 x[loc_idof] = L[2][loc_idof]*v1 + L[3][loc_idof]*v2 + L[1][loc_idof]*v3 + L[0][loc_idof]*v4;
476 std::array<point_basic<T>, 4> t1, t2;
482 t2[0] = v3-0.5*(v1+v2);
483 t2[1] = v4-0.5*(v1+v2);
484 t2[2] = v4-0.5*(v2+v3);
485 t2[3] = v4-0.5*(v1+v3);
488 for (
size_t i = 0; i < 4; i++) {
489 t1[i] /=
norm(t1[i]);
490 t2[i] /=
norm(t2[i]);
494 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> shift (loc_ndof);
495 std::vector<T> blend (loc_ndof),
497 for (
size_t ifac = 0; ifac < 4; ifac++) {
500 case 0: a = 0; b = 1; c = 2;
d = 3;
break;
501 case 1: a = 1; b = 0; c = 2;
d = 3;
break;
502 case 2: a = 2; b = 0; c = 3;
d = 1;
break;
504 default: a = 3; b = 0; c = 2;
d = 1;
break;
507 std::vector<T> warp1 (loc_ndof), warp2 (loc_ndof);
508 WarpShiftFace3D (degree, alpha, alpha, L[a], L[b], L[c], L[
d], warp1, warp2);
510 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
512 blend[loc_idof] = L[b][loc_idof]*L[c][loc_idof]*L[
d][loc_idof];
514 denom[loc_idof] = (L[b][loc_idof] + 0.5*L[a][loc_idof])
515 *(L[c][loc_idof] + 0.5*L[a][loc_idof])
516 *(L[
d][loc_idof] + 0.5*L[a][loc_idof]);
517 if (denom[loc_idof] > tol) {
518 blend[loc_idof] = (1 + sqr(alpha*L[a][loc_idof]))*blend[loc_idof]/denom[loc_idof];
522 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
523 shift[loc_idof] += blend[loc_idof]*warp1 [loc_idof]*t1 [ifac]
524 + blend[loc_idof]*warp2 [loc_idof]*t2 [ifac];
527 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
529 if (L[b][loc_idof] > tol) n_sup++;
530 if (L[c][loc_idof] > tol) n_sup++;
531 if (L[
d][loc_idof] > tol) n_sup++;
532 if ((L[a][loc_idof] < tol) && (n_sup < 3)) {
533 shift[loc_idof] = warp1[loc_idof]*t1[ifac] + warp2[loc_idof]*t2[ifac];
537 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
538 x[loc_idof] += shift[loc_idof];
540 if (!map_on_reference_element)
return;
543 T height = sqrt(
T(6))/3*side_length;
545 b ( 1, -1/sqrt(
T(3)), -height/4),
546 c ( 0, 2/sqrt(
T(3)), -height/4),
547 d ( 0, 0, 3*height/4);
550 A.set_column (b-a, 0);
551 A.set_column (c-a, 1);
552 A.set_column (
d-a, 2);
554 std::vector<size_t> perm(loc_ndof);
555 for (
size_t loc_idof = 0, k = 0; k <= degree; k++) {
556 for (
size_t j = 0; j+k <= degree; j++) {
557 for (
size_t i = 0; i+j+k <= degree; i++) {
562 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> x_new (x.size());
563 for (
size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
564 size_t new_loc_idof = perm[loc_idof];
565 x_new[new_loc_idof] = inv_A*(x[loc_idof]-a);
578 bool map_on_reference_element =
true)
589 warburton_e (degree, hat_xnod, map_on_reference_element);
break;
591 warburton_q (degree, hat_xnod, map_on_reference_element);
break;
593 warburton_H (degree, hat_xnod, map_on_reference_element);
break;
595 warburton_t (degree, hat_xnod, map_on_reference_element);
break;
597 warburton_P (degree, hat_xnod, map_on_reference_element);
break;
599 warburton_T (degree, hat_xnod, map_on_reference_element);
break;
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
see the reference_element page for the full documentation
static const variant_type H
static const variant_type q
static const variant_type e
static const variant_type p
variant_type variant() const
static size_type n_node(variant_type variant, size_type order)
static const variant_type T
static const variant_type P
static const variant_type t
#define error_macro(message)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
void warburton_T(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool map_on_reference_element=true)
void reference_element_barycenter(reference_element hat_K, point_basic< T > &c)
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)
void evalwarp(size_t degree, const std::vector< T > &xnodes, const std::vector< T > &xout, std::vector< T > &warp)
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)
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)
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
void warburton_H(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool dummy=true)
void warburton_t(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool map_on_reference_element=true)
void warburton_e(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool map_on_reference_element=true)
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)
void warburton_P(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool map_on_reference_element=true)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
void equispaced_T(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x)
void warburton_q(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool dummy=true)
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)