Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
piola_util.cc
Go to the documentation of this file.
1
21#include "rheolef/piola_util.h"
22#include "rheolef/geo_domain.h"
23#include "rheolef/inv_piola.h"
24#include "rheolef/damped_newton.h"
25namespace rheolef {
26
27// =========================================================================
28// part 1. piola transformation
29// =========================================================================
30// F_K : hat_K --> K
31// hat_x --> x
32//
33// pre-evaluation of the piola basis on a predefined point set
34// e.g. quadrature nodes hat_x[q], q=0..nq on hat_K
35// then, fast transformation of all hat_x[q] into xq on any K
36//
37// x = F(hat_x) = sum_j phi_j(hat_x)*xjnod
38// for all hat_x in the pointset : hat_xi
39// where xjnod: j-th node of element K in mesh omega
40//
41template<class T, class M>
42void
44 const geo_basic<T,M>& omega,
47 const std::vector<size_t>& dis_inod,
48 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x)
49{
50 typedef typename geo_basic<T,M>::size_type size_type;
51 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& phij_xi = piola_on_pointset.template evaluate<T> (hat_K);
52 size_type loc_nnod = phij_xi.rows();
53 size_type loc_ndof = phij_xi.cols();
54 x.resize (loc_nnod);
55 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
56 x[loc_inod] = point_basic<T>(0,0,0);
57 }
58 size_type d = omega.dimension();
59 for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
60 // dis_node: in outer loop: could require more time with external node
61 const point_basic<T>& xjnod = omega.dis_node (dis_inod[loc_jdof]);
62 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
63 for (size_type alpha = 0; alpha < d; alpha++) {
64 x[loc_inod][alpha] += phij_xi (loc_inod,loc_jdof)*xjnod[alpha];
65 }
66 }
67 }
68}
69// ------------------------------------------
70// jacobian of the piola transformation
71// at a all quadrature points
72// ------------------------------------------
73//
74// DF(hat_x) = sum_j grad_phi_j(hat_x)*xjnod
75// for all hat_x in the pointset : hat_xi
76// where xjnod: j-th node of element K in mesh omega
77//
78template<class T, class M>
79void
81 const geo_basic<T,M>& omega,
84 const std::vector<size_t>& dis_inod,
85 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& DF)
86{
87 typedef typename geo_basic<T,M>::size_type size_type;
88 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>&
89 grad_phij_xi = piola_on_pointset.template grad_evaluate<point_basic<T>>(hat_K);
90 size_type loc_nnod = grad_phij_xi.rows();
91 size_type loc_ndof = grad_phij_xi.cols();
92 DF.resize (loc_nnod);
93 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
94 DF[loc_inod].fill (0);
95 }
96 size_type d = omega.dimension();
97 size_type map_d = hat_K.dimension();
98 for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
99 // dis_node: in outer loop: could require more time with external node
100 const point_basic<T>& xjnod = omega.dis_node (dis_inod[loc_jdof]);
101 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
102 cumul_otimes (DF[loc_inod], xjnod, grad_phij_xi(loc_inod,loc_jdof), d, map_d);
103 }
104 }
105}
106// ------------------------------------------
107// jacobian of the piola transformation
108// at an aritrarily point hat_x
109// ------------------------------------------
110//
111// DF(hat_x) = sum_j grad_phi_j(hat_x)*xjnod
112// where xjnod: j-th node of element K in mesh omega
113//
114template<class T, class M>
115void
117 const geo_basic<T,M>& omega,
118 const basis_basic<T>& piola_basis,
119 reference_element hat_K,
120 const std::vector<size_t>& dis_inod,
121 const point_basic<T>& hat_x,
122 tensor_basic<T>& DF)
123{
124 typedef typename geo_basic<T,M>::size_type size_type;
125 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> grad_phi_x;
126 piola_basis.grad_evaluate (hat_K, hat_x, grad_phi_x);
127 DF.fill (0);
128 size_type d = omega.dimension();
129 size_type map_d = hat_K.dimension();
130 for (size_type loc_jdof = 0, loc_ndof = dis_inod.size(); loc_jdof < loc_ndof; loc_jdof++) {
131 const point_basic<T>& xjnod = omega.dis_node (dis_inod[loc_jdof]);
132 cumul_otimes (DF, xjnod, grad_phi_x [loc_jdof], d, map_d);
133 }
134}
135template <class T>
136T
137det_jacobian_piola_transformation (const tensor_basic<T>& DF, size_t d , size_t map_d)
138{
139 if (d == map_d) {
140 return DF.determinant (map_d);
141 }
142 /* surface jacobian: references:
143 * Spectral/hp element methods for CFD
144 * G. E. M. Karniadakis and S. J. Sherwin
145 * Oxford university press
146 * 1999
147 * page 165
148 */
149 switch (map_d) {
150 case 0: return 1;
151 case 1: return norm(DF.col(0));
152 case 2: return norm(vect(DF.col(0), DF.col(1)));
153 default:
154 error_macro ("det_jacobian_piola_transformation: unsupported element dimension "
155 << map_d << " in " << d << "D mesh.");
156 return 0;
157 }
158}
159// ------------
160// normal
161// ------------
162template<class T, class M>
163static
164point_basic<T>
165normal_from_piola_transformation_1d (
166 const geo_basic<T,M>& omega,
167 const geo_element& S,
168 const tensor_basic<T>& DF,
169 size_t d)
170{
171 // point in 1D: DF[1][0] is empty, so scan S[0] node connectivity
172 typedef typename geo_basic<T,M>::size_type size_type;
173 if (S.dimension() + 1 == omega.map_dimension()) {
174 size_type dis_ie = S.master(0);
175 check_macro (dis_ie != std::numeric_limits<size_type>::max(), "normal: requires neighbours initialization");
176 const geo_element& K = omega.dis_get_geo_element (S.dimension()+1, dis_ie);
177 Float sign = (S[0] == K[1]) ? 1 : -1;
178 return point_basic<T>(sign);
179 }
180 // omega is a domain of sides, as "boundary" or "internal_sides":
181 // for the side orient, we need to go back to its backgound volumic mesh
182 if (omega.variant() == geo_abstract_base_rep<T>::geo_domain_indirect) {
183 size_type dis_isid = S.dis_ie();
184 size_type first_dis_isid = omega.sizes().ownership_by_dimension[S.dimension()].first_index();
185 size_type isid = dis_isid - first_dis_isid;
186 check_macro (dis_isid >= first_dis_isid, "unexpected dis_index "<<dis_isid<<": out of local range");
187 const geo_basic<T,M>* ptr_bgd_omega = 0;
188 if (omega.get_background_geo().map_dimension() == 1) {
189 const geo_basic<T,M>& bgd_omega = omega.get_background_geo();
190 const geo_element& bgd_S = bgd_omega.get_geo_element(0, isid);
191 size_type bgd_dis_ie = bgd_S.master(0);
192 check_macro (bgd_dis_ie != std::numeric_limits<size_type>::max(),
193 "normal: bgd_S.dis_ie={"<<bgd_S.dis_ie()<<"} without neighbours; requires neighbours initialization for mesh " << bgd_omega.name());
194 const geo_element& bgd_K = bgd_omega.dis_get_geo_element (bgd_S.dimension()+1, bgd_dis_ie);
195 Float sign = (bgd_S[0] == bgd_K[1]) ? 1 : -1;
196 return point_basic<T>(sign);
197 } else {
198 // get a 1D geometry at a higher depth, e.g. for HDG 0D "boundary" domain in "square[sides]" 0D geo_domain
199 const geo_basic<T,M>& bgd_omega = omega.get_background_geo();
200 const geo_basic<T,M>& bgd2_omega = bgd_omega.get_background_geo();
201 check_macro (bgd2_omega.dimension() == 1, "unsupported depth for "<<omega.name()<<" in background domain "<<bgd_omega.name());
202 const geo_element& bgd_S = bgd_omega.get_geo_element(0, isid);
203 const geo_element& bgd2_S = bgd_omega.dom2bgd_geo_element (bgd_S);
204 size_type bgd2_dis_ie = bgd2_S.master(0);
205 check_macro (bgd2_dis_ie != std::numeric_limits<size_type>::max(),
206 "normal: bgd2_S.dis_ie={"<<bgd2_S.dis_ie()<<"} without neighbours; requires neighbours initialization for mesh " << bgd2_omega.name());
207 const geo_element& bgd2_K = bgd2_omega.dis_get_geo_element (bgd2_S.dimension()+1, bgd2_dis_ie);
208 Float sign = (bgd2_S[0] == bgd2_K[1]) ? 1 : -1;
209 return point_basic<T>(sign);
210 }
211 }
212 // omega.variant() != geo_abstract_base_rep<T>::geo_domain_indirect
213 size_type dis_isid = S.dis_ie();
214 size_type first_dis_isid = omega.sizes().ownership_by_dimension[S.dimension()].first_index();
215 size_type isid = dis_isid - first_dis_isid;
216 check_macro (dis_isid >= first_dis_isid, "unexpected dis_index "<<dis_isid<<": out of local range");
217 const geo_basic<T,M>& bgd_omega = omega.get_background_geo();
218 const geo_basic<T,M>& bgd_dom = omega.get_background_domain();
219 const geo_element& bgd_S = bgd_dom[isid]; // TODO: pas clair, differe du cas precedent ?
220 size_type bgd_dis_ie = bgd_S.master(0);
221 check_macro (bgd_dis_ie != std::numeric_limits<size_type>::max(),
222 "normal: bgd_S.dis_ie={"<<bgd_S.dis_ie()<<"} without neighbours; requires neighbours initialization for mesh " << bgd_omega.name());
223 const geo_element& bgd_K = bgd_omega.dis_get_geo_element (bgd_S.dimension()+1, bgd_dis_ie);
224 Float sign = (bgd_S[0] == bgd_K[1]) ? 1 : -1;
225 return point_basic<T>(sign);
226}
227template<class T, class M>
228point_basic<T>
230 const geo_basic<T,M>& omega,
231 const geo_element& S,
232 const tensor_basic<T>& DF,
233 size_t d)
234{
235 switch (d) {
236 case 1: {
237 // special case:
238 return normal_from_piola_transformation_1d (omega, S, DF, d);
239 }
240 case 2: { // edge in 2D
241 // 2D: S=edge(a,b) then t=b-a, DF=[t] and n = (t1,-t0) => det(n,1)=1 : (n,t) is direct
242 // and the normal goes outside on a boundary edge S, when the associated element K is well oriented
243 point_basic<T> t = DF.col(0);
244 t /= norm(t);
245 return point_basic<T>(t[1], -t[0]);
246 }
247 case 3: { // 3D: S=triangle(a,b,c) then t0=b-a, t1=c-a, DF=[t0,t1] and n = t0^t1/|t0^t1|.
248 point_basic<T> t0 = DF.col(0);
249 point_basic<T> t1 = DF.col(1);
250 point_basic<T> n = vect (t0,t1);
251 n /= norm(n);
252 return n;
253 }
254 default: {
255 error_macro ("normal: unsupported " << d << "D mesh.");
256 return point_basic<T>();
257 }
258 }
259}
260// The pseudo inverse extend inv(DF) for face in 3d or edge in 2d
261// i.e. useful for Laplacian-Beltrami and others surfacic forms.
262//
263// pinvDF (hat_xq) = inv DF, if tetra in 3d, tri in 2d, etc
264// = pseudo-invese, when tri in 3d, edge in 2 or 3d
265// e.g. on 3d face : pinvDF*DF = [1, 0, 0; 0, 1, 0; 0, 0, 0]
266//
267// let DF = [u, v, w], where u, v, w are the column vectors of DF
268// then det(DF) = mixt(u,v,w)
269// and det(DF)*inv(DF)^T = [v^w, w^u, u^v] where u^v = vect(u,v)
270//
271// application:
272// if K=triangle(a,b,c) then u=ab=b-a, v=ac=c-a and w = n = u^v/|u^v|.
273// Thus DF = [ab,ac,n] and det(DF)=|ab^ac|
274// and inv(DF)^T = [ac^n/|ab^ac|, -ab^n/|ab^ac|, n]
275// The pseudo-inverse is obtained by remplacing the last column n by zero.
276//
277template<class T>
278tensor_basic<T>
280 const tensor_basic<T>& DF,
281 size_t d,
282 size_t map_d)
283{
284 if (d == map_d) {
285 return inv(DF, map_d);
286 }
287 tensor_basic<T> invDF;
288 switch (map_d) {
289 case 0: { // point in 1D
290 invDF(0,0) = 1;
291 return invDF;
292 }
293 case 1: { // segment in 2D
294 point_basic<T> t = DF.col(0);
295 invDF.set_row (t/norm2(t), 0, d);
296 return invDF;
297 }
298 case 2: {
299 point_basic<T> t0 = DF.col(0);
300 point_basic<T> t1 = DF.col(1);
301 point_basic<T> n = vect(t0,t1);
302 T det2 = norm2(n);
303 point_basic<T> v0 = vect(t1,n)/det2;
304 point_basic<T> v1 = - vect(t0,n)/det2;
305 invDF.set_row (v0, 0, d);
306 invDF.set_row (v1, 1, d);
307 return invDF;
308 }
309 default:
310 error_macro ("pseudo_inverse_jacobian_piola_transformation: unsupported element dimension "
311 << map_d << " in " << d << "D mesh.");
312 return invDF;
313 }
314}
315
316
317// axisymetric weight ?
318// point_basic<T> xq = rheolef::piola_transformation (_omega, _piola_table, K, dis_inod, q);
319template<class T>
320T
322{
323 switch (sys_coord) {
324 case space_constant::axisymmetric_rz: return xq[0];
325 case space_constant::axisymmetric_zr: return xq[1];
326 case space_constant::cartesian: return 1;
327 default: {
328 fatal_macro ("unsupported coordinate system `"
329 << space_constant::coordinate_system_name (sys_coord) << "'");
330 return 0;
331 }
332 }
333}
334// -------------------------------------------
335// weight integration: w = det_DF*wq
336// with optional axisymmetric r*dr factor
337// -------------------------------------------
338template<class T, class M>
339void
341 const geo_basic<T,M>& omega,
342 const basis_on_pointset<T>& piola_on_quad,
343 reference_element hat_K,
344 const std::vector<size_t>& dis_inod,
345 bool ignore_sys_coord,
346 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& DF,
347 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x,
348 Eigen::Matrix<T,Eigen::Dynamic,1>& w)
349{
350 typedef typename geo_basic<T,M>::size_type size_type;
351 jacobian_piola_transformation (omega, piola_on_quad, hat_K, dis_inod, DF);
352 size_type loc_nnod = piola_on_quad.nnod (hat_K);
353 w.resize (loc_nnod);
354 if (omega.coordinate_system() == space_constant::cartesian || ignore_sys_coord) {
355 w.fill (T(1));
356 } else {
357 piola_transformation (omega, piola_on_quad, hat_K, dis_inod, x);
358 size_t k = (omega.coordinate_system() == space_constant::axisymmetric_rz) ? 0 : 1;
359 for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
360 w[loc_inod] = x[loc_inod][k];
361 }
362 }
363 size_type d = omega.dimension();
364 size_type map_d = hat_K.dimension();
365 const quadrature<T>& quad = piola_on_quad.get_quadrature();
367 first_quad = quad.begin(hat_K),
368 last_quad = quad.end (hat_K);
369 for (size_type q = 0; first_quad != last_quad; ++first_quad, ++q) {
370 T det_DF = det_jacobian_piola_transformation (DF[q], d, map_d);
371 T wq = det_DF*(*first_quad).w;
372 if (! ignore_sys_coord) {
373 w[q] = wq*w[q];
374 } else {
375 w[q] = wq;
376 }
377 }
378}
379// calcul P = I - nxn
380template<class T>
381void
383 const tensor_basic<T>& DF,
384 size_t d,
385 size_t map_d,
387{
389 switch (map_d) {
390 case 1: {
391 point_basic<Float> t = DF.col(0);
392 check_macro (d == map_d+1, "unexpected dimension map_d="<<map_d<<" and d="<<d);
393 n = point_basic<T>(t[1],-t[0]);
394 break;
395 }
396 case 2: {
397 point_basic<Float> t0 = DF.col(0);
398 point_basic<Float> t1 = DF.col(1);
399 n = vect(t0,t1);
400 break;
401 }
402 default:
403 error_macro ("unexpected dimension "<<map_d);
404 }
405 n = n/norm(n);
406 for (size_t l = 0; l < d; l++) {
407 for (size_t m = 0; m < d; m++) {
408 P(l,m) = - n[l]*n[m];
409 }
410 P(l,l) += 1;
411 }
412}
413// =========================================================================
414// part 2. inverse piola transformation
415// =========================================================================
416// F_K^{-1} : K --> hat(K)
417// x --> hat(x)
418// TODO: non-linear case
419template<class T>
420static
421inline
422point_basic<T>
423inv_piola_e (
424 const point_basic<T>& x,
425 const point_basic<T>& a,
426 const point_basic<T>& b)
427{
428 return point_basic<T>((x[0]-a[0])/(b[0]-a[0]));
429}
430template<class T>
431static
432inline
433point_basic<T>
434inv_piola_t (
435 const point_basic<T>& x,
436 const point_basic<T>& a,
437 const point_basic<T>& b,
438 const point_basic<T>& c)
439{
440 T t9 = 1/(-b[0]*c[1]+b[0]*a[1]+a[0]*c[1]+c[0]*b[1]-c[0]*a[1]-a[0]*b[1]);
441 T t11 = -a[0]+x[0];
442 T t15 = -a[1]+x[1];
443 return point_basic<T>((-c[1]+a[1])*t9*t11-(-c[0]+a[0])*t9*t15,
444 ( b[1]-a[1])*t9*t11-( b[0]-a[0])*t9*t15);
445}
446template<class T>
447static
448inline
449point_basic<T>
450inv_piola_T (
451 const point_basic<T>& x,
452 const point_basic<T>& a,
453 const point_basic<T>& b,
454 const point_basic<T>& c,
455 const point_basic<T>& d)
456{
457 tensor_basic<T> A;
458 point_basic<T> ax;
459 for (size_t i = 0; i < 3; i++) {
460 ax[i] = x[i]-a[i];
461 A(i,0) = b[i]-a[i];
462 A(i,1) = c[i]-a[i];
463 A(i,2) = d[i]-a[i];
464 }
465 tensor_basic<T> inv_A;
466 bool is_singular = ! invert_3x3 (A, inv_A);
467 check_macro(!is_singular, "inv_piola: singular transformation in a tetrahedron");
468 point_basic<T> hat_x = inv_A*ax;
469 return hat_x;
470}
471template<class T, class M>
472point_basic<T>
474 const geo_basic<T,M>& omega,
475 const reference_element& hat_K,
476 const std::vector<size_t>& dis_inod,
477 const point_basic<T>& x)
478{
479 check_macro (omega.order() == 1, "inverse piola: mesh order > 1: not yet");
480 if (omega.order() == 1) {
481 switch (hat_K.variant()) {
482 case reference_element::e: return inv_piola_e (x, omega.dis_node(dis_inod [0]),
483 omega.dis_node(dis_inod [1]));
484 case reference_element::t: return inv_piola_t (x, omega.dis_node(dis_inod [0]),
485 omega.dis_node(dis_inod [1]),
486 omega.dis_node(dis_inod [2]));
487 case reference_element::T: return inv_piola_T (x, omega.dis_node(dis_inod [0]),
488 omega.dis_node(dis_inod [1]),
489 omega.dis_node(dis_inod [2]),
490 omega.dis_node(dis_inod [3]));
491 }
492 }
493 // non-linear transformation: q,P,H or high order > 1 => use Newton
494 inv_piola<T> F;
495 F.reset (omega, hat_K, dis_inod);
496 F.set_x (x);
497 point_basic<T> hat_x = F.initial();
498 size_t max_iter = 500, n_iter = max_iter;
499 T tol = std::numeric_limits<Float>::epsilon(), r = tol;
500 int status = damped_newton (F, hat_x, r, n_iter);
501 check_macro (status == 0, "inv_piola: newton failed (residue="<<r<<", n_iter="<<n_iter<<")");
502 return hat_x;
503}
504// ----------------------------------------------------------------------------
505// instanciation in library
506// ----------------------------------------------------------------------------
507#define _RHEOLEF_instanciation1(T) \
508template \
509T \
510det_jacobian_piola_transformation ( \
511 const tensor_basic<T>& DF, \
512 size_t d, \
513 size_t map_d); \
514template \
515tensor_basic<T> \
516pseudo_inverse_jacobian_piola_transformation ( \
517 const tensor_basic<T>& DF, \
518 size_t d, \
519 size_t map_d); \
520template \
521T \
522weight_coordinate_system ( \
523 space_constant::coordinate_type sys_coord, \
524 const point_basic<T>& xq); \
525template \
526void \
527map_projector ( \
528 const tensor_basic<T>& DF, \
529 size_t d, \
530 size_t map_d, \
531 tensor_basic<T>& P); \
532
533
534#define _RHEOLEF_instanciation2(T,M) \
535template \
536void \
537piola_transformation ( \
538 const geo_basic<T,M>& omega, \
539 const basis_on_pointset<T>& piola_on_pointset, \
540 reference_element hat_K, \
541 const std::vector<size_t>& dis_inod, \
542 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x); \
543template \
544point_basic<T> \
545inverse_piola_transformation ( \
546 const geo_basic<T,M>& omega, \
547 const reference_element& hat_K, \
548 const std::vector<size_t>& dis_inod, \
549 const point_basic<T>& x); \
550template \
551void \
552jacobian_piola_transformation ( \
553 const geo_basic<T,M>& omega, \
554 const basis_basic<T>& piola_basis, \
555 reference_element hat_K, \
556 const std::vector<size_t>& dis_inod, \
557 const point_basic<T>& hat_x, \
558 tensor_basic<T>& DF); \
559template \
560void \
561jacobian_piola_transformation ( \
562 const geo_basic<T,M>& omega, \
563 const basis_on_pointset<T>& piola_on_pointset, \
564 reference_element hat_K, \
565 const std::vector<size_t>& dis_inod, \
566 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& DF); \
567template \
568point_basic<T> \
569normal_from_piola_transformation ( \
570 const geo_basic<T,M>& omega, \
571 const geo_element& S, \
572 const tensor_basic<T>& DF, \
573 size_t d); \
574template \
575void \
576piola_transformation_and_weight_integration ( \
577 const geo_basic<T,M>& omega, \
578 const basis_on_pointset<T>& piola_on_pointset, \
579 reference_element hat_K, \
580 const std::vector<size_t>& dis_inod, \
581 bool ignore_sys_coord, \
582 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& DF, \
583 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, \
584 Eigen::Matrix<T,Eigen::Dynamic,1>& w); \
585
586
589#ifdef _RHEOLEF_HAVE_MPI
591#endif // _RHEOLEF_HAVE_MPI
592
593} // namespace rheolef
field::size_type size_type
Definition branch.cc:430
see the Float page for the full documentation
void grad_evaluate(reference_element hat_K, const point_basic< T > &hat_x, Eigen::Matrix< Value, Eigen::Dynamic, 1 > &value) const
Definition basis.h:942
size_type nnod(reference_element hat_K) const
const quadrature< T > & get_quadrature() const
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
void set_x(const value_type &x1)
Definition inv_piola.h:41
void reset(const geo_basic< T, M > &omega, const reference_element &hat_K, const std::vector< size_t > &dis_inod)
Definition inv_piola.h:77
value_type initial() const
Definition inv_piola.h:89
rep::const_iterator const_iterator
Definition quadrature.h:195
const_iterator end(reference_element hat_K) const
Definition quadrature.h:219
const_iterator begin(reference_element hat_K) const
Definition quadrature.h:218
see the reference_element page for the full documentation
static const variant_type e
variant_type variant() const
static const variant_type T
static const variant_type t
point_basic< T > col(size_type i) const
Definition tensor.cc:323
T determinant(size_type d=3) const
Definition tensor.cc:288
void fill(const T &init_val)
Definition tensor.h:252
void set_row(const point_basic< T > &r, size_t i, size_t d=3)
Definition tensor.h:435
#define error_macro(message)
Definition dis_macros.h:49
#define fatal_macro(message)
Definition dis_macros.h:33
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)")
std::string coordinate_system_name(coordinate_type i)
This file is part of Rheolef.
point_basic< T > normal_from_piola_transformation(const geo_basic< T, M > &omega, const geo_element &S, const tensor_basic< T > &DF, size_t d)
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition tensor.cc:219
void piola_transformation(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_pointset, reference_element hat_K, const std::vector< size_t > &dis_inod, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x)
Definition piola_util.cc:43
tensor_basic< T > pseudo_inverse_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
void map_projector(const tensor_basic< T > &DF, size_t d, size_t map_d, tensor_basic< T > &P)
T weight_coordinate_system(space_constant::coordinate_type sys_coord, const point_basic< T > &xq)
point_basic< T > vect(const point_basic< T > &v, const point_basic< T > &w)
Definition point.h:264
void piola_transformation_and_weight_integration(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_quad, reference_element hat_K, const std::vector< size_t > &dis_inod, bool ignore_sys_coord, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &DF, Eigen::Matrix< point_basic< T >, Eigen::Dynamic, 1 > &x, Eigen::Matrix< T, Eigen::Dynamic, 1 > &w)
int damped_newton(const Problem &P, const Preconditioner &T, Field &u, Real &tol, Size &max_iter, odiststream *p_derr=0)
see the damped_newton page for the full documentation
T norm2(const vec< T, M > &x)
norm2(x): see the expression page for the full documentation
Definition vec.h:379
void jacobian_piola_transformation(const geo_basic< T, M > &omega, const basis_on_pointset< T > &piola_on_pointset, reference_element hat_K, const std::vector< size_t > &dis_inod, Eigen::Matrix< tensor_basic< T >, Eigen::Dynamic, 1 > &DF)
Definition piola_util.cc:80
T det_jacobian_piola_transformation(const tensor_basic< T > &DF, size_t d, size_t map_d)
point_basic< T > inverse_piola_transformation(const geo_basic< T, M > &omega, const reference_element &hat_K, const std::vector< size_t > &dis_inod, const point_basic< T > &x)
void cumul_otimes(tensor_basic< T > &t, const point_basic< T > &a, const point_basic< T > &b, size_t na, size_t nb)
Definition tensor.cc:305
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
Definition vec.h:387
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
Definition tensor.cc:333
#define _RHEOLEF_instanciation1(T)
#define _RHEOLEF_instanciation2(T, M)