Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
warburton.icc
Go to the documentation of this file.
1#ifndef _RHEOLEF_WARBURTON_ICC
2#define _RHEOLEF_WARBURTON_ICC
23//
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
37namespace rheolef {
38
39// ---------------------------------------------------------------------------
40// 1d : edge
41// ---------------------------------------------------------------------------
42// Compute nodes in [-1:1] edge for polynomial of order "degree"
43template<class T>
44void
45warburton_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"
71template<class T>
72void
73warburton_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"
89template<class T>
90void
91warburton_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);
108template<class T, class Matrix>
109static
110void
111trans_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
121template<class T>
122static
123void
124warpfactor (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)]
174template<class T>
175void
176warburton_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]
274template<class T>
275void
276warburton_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
301template<class T>
302void
303equispaced_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
321template<class T>
322void
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
361template<class T>
362void
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
424template<class T>
425inline
426void
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
441template<class T>
442void
443warburton_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// ---------------------------------------------------------------------------
572template<class T>
573void
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
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)
Definition dis_macros.h:49
Expr1::float_type T
Definition field_expr.h:230
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)
Definition tiny_lu.h:127
void warburton_H(size_t degree, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, bool dummy=true)
Definition warburton.icc:91
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)
Definition warburton.icc:45
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
Definition vec.h:387
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)
Definition warburton.icc:73
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition tensor.cc:333
size_t N