Rheolef  7.1
an efficient C++ finite element environment
warburton.icc
Go to the documentation of this file.
1 #ifndef _RHEOLEF_WARBURTON_ICC
2 #define _RHEOLEF_WARBURTON_ICC
3 //
24 // compute the Warburton warp and blend point set (see War-2006 HesWar-2008)
25 //
26 // author: Pierre.Saramito@imag.fr
27 //
28 // date: 2 september 2017
29 //
30 #include "rheolef/point.h"
31 #include "rheolef/tensor.h"
32 #include "rheolef/reference_element.h"
33 #include "rheolef/gauss_lobatto_jacobi.h"
35 #include "eigen_util.h"
36 
37 namespace rheolef {
38 
39 // ---------------------------------------------------------------------------
40 // 1d : edge
41 // ---------------------------------------------------------------------------
42 // Compute nodes in [-1:1] edge for polynomial of order "degree"
43 template<class T>
44 void
45 warburton_e (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool map_on_reference_element=true)
46 {
47  // Compute Gauss-Lobatto Legendre node distribution
48  std::vector<T> zeta(degree+1), omega(degree+1);
49  gauss_lobatto_jacobi (degree+1, T(0), T(0), zeta.begin(), omega.begin());
50  x.resize (degree+1);
51  for (size_t i = 0; i < degree+1; i++) {
52  x[i] = point_basic<T>(zeta[i],0,0);
53  }
54  if (!map_on_reference_element) return;
55  // permutation: by decreasing subgeo dimension, for rheolef
56  // transform: to reference_element
57  std::vector<size_t> perm(degree+1);
58  for (size_t i = 0; i < degree+1; i++) {
60  }
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]);
64  }
65  x = x_new;
66 }
67 // ---------------------------------------------------------------------------
68 // 2d : quadrangle
69 // ---------------------------------------------------------------------------
70 // Compute nodes in [-1:1]^2 quadrangle for polynomial of order "degree"
71 template<class T>
72 void
73 warburton_q (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool dummy=true)
74 {
75  // Compute Gauss-Lobatto Legendre node distribution
76  std::vector<T> zeta(degree+1), omega(degree+1);
77  gauss_lobatto_jacobi (degree+1, T(0), T(0), zeta.begin(), omega.begin());
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++) {
81  size_t loc_idof = reference_element_q::ilat2loc_inod (degree, point_basic<size_t>(i,j));
82  x[loc_idof] = point_basic<T>(zeta[i], zeta[j],0);
83  }}
84 }
85 // ---------------------------------------------------------------------------
86 // 3d : hexahedron
87 // ---------------------------------------------------------------------------
88 // Compute nodes in [-1:1]^3 quadrangle for polynomial of order "degree"
89 template<class T>
90 void
91 warburton_H (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool dummy=true)
92 {
93  // Compute Gauss-Lobatto Legendre node distribution
94  std::vector<T> zeta(degree+1), omega(degree+1);
95  gauss_lobatto_jacobi (degree+1, T(0), T(0), zeta.begin(), omega.begin());
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++) {
100  size_t loc_idof = reference_element_H::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
101  x[loc_idof] = point_basic<T>(zeta[i], zeta[j], zeta[k]);
102  }}}
103 }
104 // ---------------------------------------------------------------------------
105 // 2d : triangle
106 // ---------------------------------------------------------------------------
107 // initialize the transpose of the 1d Vandermonde matrix: v_ij = p_j(x_i);
108 template<class T, class Matrix>
109 static
110 void
111 trans_vandermonde1d (size_t N, const std::vector<T>& x, Matrix& trans_v)
112 {
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]);
117  }
118  }
119 }
120 // Compute scaled warp function at order N based on rout interpolation nodes
121 template<class T>
122 static
123 void
124 warpfactor (size_t N, const std::vector<T>& rout, std::vector<T>& warp)
125 {
126  // Compute equidistant node distribution
127  std::vector<T> req(N+1);
128  req[0] = -1;
129  for (size_t i = 1; i < N; i++) {
130  req[i] = -1 + 2*T(int(i))/T(int(N));
131  }
132  req[N] = 1;
133 
134  // Compute Gauss-Lobatto Legendre node distribution
135  std::vector<T> zeta(N+1), omega(N+1);
136  gauss_lobatto_jacobi (N+1, T(0), T(0), zeta.begin(), omega.begin());
137 
138  // Compute V based on req
139  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> trans_v_eq (N+1, N+1);
140  trans_vandermonde1d (N, req, trans_v_eq);
141 
142  // Evaluate Lagrange polynomial at rout: pmat(i,j) = p_i(rout[j]);
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]);
149  }
150  }
151  // pmat := inv(v_eq^T)*pmat
152  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> inv_trans_v_eq (N+1, N+1);
153  check_macro (invert(trans_v_eq, inv_trans_v_eq), "matrix inversion failed");
154  pmat = inv_trans_v_eq*pmat;
155 
156  // Compute warp factor: warp = pmat^T*(LGLr - req);
157  warp.resize (Nr);
158  for (size_t i = 0; i < Nr; i++) {
159  warp[i] = T(0);
160  for (size_t j = 0; j < N+1; j++) {
161  warp[i] += pmat(j,i)*(zeta[j] - req[j]);
162  }
163  }
164  // Scale factor
165  const T tol = 1e-10;
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);
170  }
171 }
172 // Compute nodes in equilateral triangle for polynomial of order "degree"
173 // note: equilateral triangle bbox = [-1:1]*[-1/sqrt(3):2/sqrt(3)]
174 template<class T>
175 void
176 warburton_t (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool map_on_reference_element=true)
177 {
178  // Set optimized parameter, alpha, depending on order "degree"
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);
183 
184  // total number of nodes
185  size_t loc_ndof = (degree+1)*(degree+2)/2;
186  x.resize(loc_ndof);
187 
188  // Create equidistributed nodes on equilateral triangle
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;
198  }
199  }
200  // Compute blending function at each node for each edge
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];
206  }
207  // Amount of warp for each node, for each edge
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];
214  }
215  warpfactor (degree, L32, warpf1);
216  warpfactor (degree, L13, warpf2);
217  warpfactor (degree, L21, warpf3);
218 
219  // Combine blend & warp
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]));
225  }
226  // Accumulate deformations associated with each edge
227  const T pi = acos(T(-1));
228  const T cos23 = cos(2*pi/3),
229  sin23 = sin(2*pi/3),
230  cos43 = cos(4*pi/3),
231  sin43 = sin(4*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];
235  }
236  if (!map_on_reference_element) return;
237  // transform: to reference_element, from the isocele triangle used by Warburton
238  // permutation: by decreasing subgeo dimension, for rheolef
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.;
242  // introduce the Piola transformation from ref elem to isocele
243  T a11 = xb-xa;
244  T a12 = xc-xa;
245  T a21 = yb-ya;
246  T a22 = yc-ya;
247  // invert the Piola transformation
248  T det = a11*a22-a12*a21;
249  T b11 = a22/det;
250  T b12 = -a12/det;
251  T b21 = -a21/det;
252  T b22 = a11/det;
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++) {
256  perm[loc_idof++] = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
257  }
258  }
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);
266  }
267  x = x_new;
268 }
269 // ---------------------------------------------------------------------------
270 // 3d : prism
271 // ---------------------------------------------------------------------------
272 // compute the equidistributed nodes on the reference prism
273 // note: reference prism bbox = ref_triangle*[-1:1]
274 template<class T>
275 void
276 warburton_P (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool map_on_reference_element=true)
277 {
278  // for x,y: reference triangle, but unpermuted nodes
279  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> xt;
280  warburton_t (degree, xt);
281  // for z: compute Gauss-Lobatto Legendre node distribution
282  std::vector<T> zeta(degree+1), omega(degree+1);
283  gauss_lobatto_jacobi (degree+1, T(0), T(0), zeta.begin(), omega.begin());
284  // tensorial product: ((x,y),z)
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++) {
288  size_t loc_t_idof = reference_element_t::ilat2loc_inod (degree, point_basic<size_t>(i,j));
289  for (size_t k = 0; k <= degree; k++) {
290  size_t loc_idof = reference_element_P::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
291  x[loc_idof] = point_basic<T> (xt[loc_t_idof][0], xt[loc_t_idof][1], zeta[k]);
292  }
293  }
294  }
295 }
296 // ---------------------------------------------------------------------------
297 // 3d : tetrahedron
298 // ---------------------------------------------------------------------------
299 // compute the equidistributed nodes on the reference tetrahedron
300 // note: reference tetrahedron bbox = [-1:1]^3
301 template<class T>
302 void
303 equispaced_T (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x)
304 {
305  // total number of nodes
306  size_t loc_ndof = (degree+1)*(degree+2)*(degree+3)/6;
307  x.resize (loc_ndof);
308 
309  // create equidistributed nodes on equilateral triangle
310  size_t loc_idof = 0;
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++) {
314  x[loc_idof] = point_basic<T>(
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)));
318  }
319 }
320 // compute one-dimensional edge warping function
321 template<class T>
322 void
324  size_t degree,
325  const std::vector<T>& xnodes,
326  const std::vector<T>& xout,
327  std::vector<T>& warp)
328 {
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));
332  }
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++) {
339  if (i != 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]);
342  }
343  }
344  }
345  if (i != 0) {
346  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
347  d[loc_idof] = - d[loc_idof]/(xeq[i]-xeq[0]);
348  }
349  }
350  if (i != degree) {
351  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
352  d[loc_idof] = d[loc_idof]/(xeq[i]-xeq[degree]);
353  }
354  }
355  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
356  warp[loc_idof] += d[loc_idof];
357  }
358  }
359 }
360 // compute two-dimensional Warp & Blend transform
361 template<class T>
362 void
364  size_t degree,
365  const T& alpha,
366  const std::vector<T>& L0,
367  const std::vector<T>& L1,
368  const std::vector<T>& L2,
369  std::vector<T>& dx,
370  std::vector<T>& dy)
371 {
372  // 1) compute Gauss-Lobatto-Legendre node distribution
373  std::vector<T> gauss_x(degree+1), omega(degree+1);
374  gauss_lobatto_jacobi (degree+1, T(0), T(0), gauss_x.begin(), omega.begin());
375  for (size_t i = 0; i < degree+1; i++) {
376  gauss_x[i] = - gauss_x[i];
377  }
378  // 2) compute blending function at each node for each edge
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];
385  }
386  // 3) amount of warp for each node, for each edge
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];
392  }
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;
401  }
402  // 4) combine blend & warp
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]));
408  }
409  // 5) evaluate shift in equilateral triangle
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),
414  sin23 = sin(2*pi/3),
415  cos43 = cos(4*pi/3),
416  sin43 = sin(4*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];
420  // les indices 6 & 7 differents de matlab:
421  }
422 }
423 // compute two-dimensional Warp & Blend transform
424 template<class T>
425 inline
426 void
428  size_t degree,
429  const T& alpha,
430  const T& dummy,
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)
437 {
438  evalshift (degree, alpha, L1, L2, L3, warpx, warpy);
439 }
440 // compute warp & blend nodes
441 template<class T>
442 void
443 warburton_T (size_t degree, Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, bool map_on_reference_element=true)
444 {
445  // choose optimized blending parameter
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);
450 
451  // create equidistributed nodes in the reference tetrahedron
452  equispaced_T (degree, x);
453  size_t loc_ndof = x.size();
454 
455  // compute barycentric coordinates
456  std::array<std::vector<T>, 4> L;
457  for (size_t i = 0; i < 4; i++) L[i].resize (loc_ndof);
458 
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;
464  }
465  // set vertices of the equilateral tetrahedron
466  point_basic<T> v1 (-1, -1/sqrt(T(3)), -1/sqrt(T(6))),
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)));
470 
471  // create equidistributed nodes in the equilateral tetrahedron
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;
474  }
475  // orthogonal axis tangents on faces 1-4
476  std::array<point_basic<T>, 4> t1, t2;
477  t1[0] = v2-v1;
478  t1[1] = v2-v1;
479  t1[2] = v3-v2;
480  t1[3] = v3-v1;
481 
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);
486 
487  // normalize tangents
488  for (size_t i = 0; i < 4; i++) {
489  t1[i] /= norm(t1[i]);
490  t2[i] /= norm(t2[i]);
491  }
492  // Warp and blend for each face (accumulated in shift)
493  const T tol = 1e-10;
494  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> shift (loc_ndof);
495  std::vector<T> blend (loc_ndof),
496  denom (loc_ndof);
497  for (size_t ifac = 0; ifac < 4; ifac++) {
498  size_t a, b, c, d;
499  switch (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;
503  case 3:
504  default: a = 3; b = 0; c = 2; d = 1; break;
505  }
506  // compute warp tangential to face
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);
509 
510  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
511  // compute volume blending
512  blend[loc_idof] = L[b][loc_idof]*L[c][loc_idof]*L[d][loc_idof];
513  // modify linear blend
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];
519  }
520  }
521  // compute warp & blend
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];
525  }
526  // fix face warp
527  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
528  size_t n_sup = 0;
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];
534  }
535  }
536  }
537  for (size_t loc_idof = 0; loc_idof < loc_ndof; loc_idof++) {
538  x[loc_idof] += shift[loc_idof];
539  }
540  if (!map_on_reference_element) return;
541  // define the isocele tetra used by Warburton
542  T side_length = 2;
543  T height = sqrt(T(6))/3*side_length;
544  point_basic<T> a (-1, -1/sqrt(T(3)), -height/4),
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);
548  // define and invert the Piola transform
549  tensor_basic<T> A, inv_A;
550  A.set_column (b-a, 0);
551  A.set_column (c-a, 1);
552  A.set_column (d-a, 2);
553  invert_3x3 (A,inv_A);
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++) {
558  perm[loc_idof++] = reference_element_T::ilat2loc_inod (degree, point_basic<size_t>(i,j,k));
559  }
560  }
561  }
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);
566  }
567  x = x_new;
568 }
569 // ---------------------------------------------------------------------------
570 // all cases
571 // ---------------------------------------------------------------------------
572 template<class T>
573 void
575  reference_element hat_K,
576  size_t degree,
577  Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& hat_xnod,
578  bool map_on_reference_element = true)
579 {
580  hat_xnod.resize (reference_element::n_node(hat_K.variant(), degree));
581  if (degree == 0) {
582  reference_element_barycenter (hat_K, hat_xnod[0]);
583  return;
584  }
585  switch (hat_K.variant()) {
587  hat_xnod [0] = point_basic<T> (T(0)); break;
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;
600  default: error_macro ("unexpected element type `"<<hat_K.name()<<"'");
601  }
602 }
603 
604 } // namespace rheolef
605 #endif // _RHEOLEF_WARBURTON_ICC
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::e
static const variant_type e
Definition: reference_element.h:76
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::H
static const variant_type H
Definition: reference_element.h:81
rheolef::reference_element::n_node
static size_type n_node(variant_type variant, size_type order)
Definition: reference_element.cc:201
eigen_util.h
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
rheolef::invert_3x3
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition: tensor.cc:333
rheolef::point_basic
Definition: point.h:87
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
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
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_T::ilat2loc_inod
static size_type ilat2loc_inod(size_type order, const point_basic< size_type > &ilat)
Definition: reference_element.cc:582
N
size_t N
Definition: neumann-laplace-check.cc:24
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::dummy
static iorheo::force_initialization dummy
Definition: iorheo.cc:147
rheolef::reference_element::T
static const variant_type T
Definition: reference_element.h:79
rheolef::tensor_basic
Definition: tensor.h:90
mkgeo_ball.c
c
Definition: mkgeo_ball.sh:153
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::norm
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition: vec.h:387
reference_element_aux.icc
rheolef::reference_element
see the reference_element page for the full documentation
Definition: reference_element.h:66
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::equispaced_T
void equispaced_T(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x)
Definition: warburton.icc:303
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
a
Definition: diffusion_isotropic.h:25
rheolef::reference_element::name
char name() const
Definition: reference_element.h:100
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
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::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::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
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
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
mkgeo_obstacle.L
L
Definition: mkgeo_obstacle.sh:133
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::point_basic< T >
point_basic< T >
Definition: piola_fem.h:135
mkgeo_ball.b
b
Definition: mkgeo_ball.sh:152
mkgeo_ball.a
a
Definition: mkgeo_ball.sh:151
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::reference_element::q
static const variant_type q
Definition: reference_element.h:78
rheolef::reference_element::P
static const variant_type P
Definition: reference_element.h:80
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
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::reference_element::t
static const variant_type t
Definition: reference_element.h:77
rheolef::invert
void invert(tiny_matrix< T > &a, tiny_matrix< T > &inv_a)
Definition: tiny_lu.h:127
rheolef::reference_element_barycenter
void reference_element_barycenter(reference_element hat_K, point_basic< T > &c)
Definition: reference_element_aux.icc:38
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::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::tensor_basic::set_column
void set_column(const point_basic< T > &c, size_t j, size_t d=3)
Definition: tensor.h:427
T
Expr1::float_type T
Definition: field_expr.h:218