21#include "rheolef/piola_util.h"
22#include "rheolef/geo_domain.h"
23#include "rheolef/inv_piola.h"
24#include "rheolef/damped_newton.h"
41template<
class T,
class M>
47 const std::vector<size_t>& dis_inod,
51 const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& phij_xi =
piola_on_pointset.template evaluate<T> (hat_K);
55 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
59 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++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];
78template<
class T,
class M>
84 const std::vector<size_t>& dis_inod,
88 const Eigen::Matrix<point_basic<T>,Eigen::Dynamic,Eigen::Dynamic>&
93 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
94 DF[loc_inod].fill (0);
98 for (
size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
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);
114template<
class T,
class M>
120 const std::vector<size_t>& dis_inod,
125 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1> grad_phi_x;
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]);
151 case 1:
return norm(DF.
col(0));
154 error_macro (
"det_jacobian_piola_transformation: unsupported element dimension "
155 << map_d <<
" in " <<
d <<
"D mesh.");
162template<
class T,
class M>
165normal_from_piola_transformation_1d (
166 const geo_basic<T,M>& omega,
167 const geo_element& S,
168 const tensor_basic<T>& DF,
172 typedef typename geo_basic<T,M>::size_type
size_type;
173 if (S.dimension() + 1 == omega.map_dimension()) {
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);
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);
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);
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);
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];
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);
227template<
class T,
class M>
238 return normal_from_piola_transformation_1d (omega, S, DF,
d);
285 return inv(DF, map_d);
310 error_macro (
"pseudo_inverse_jacobian_piola_transformation: unsupported element dimension "
311 << map_d <<
" in " <<
d <<
"D mesh.");
338template<
class T,
class M>
344 const std::vector<size_t>& dis_inod,
345 bool ignore_sys_coord,
348 Eigen::Matrix<T,Eigen::Dynamic,1>& w)
359 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
360 w[loc_inod] = x[loc_inod][k];
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) {
371 T wq = det_DF*(*first_quad).w;
372 if (! ignore_sys_coord) {
392 check_macro (
d == map_d+1,
"unexpected dimension map_d="<<map_d<<
" and d="<<
d);
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];
424 const point_basic<T>& x,
425 const point_basic<T>& a,
426 const point_basic<T>& b)
428 return point_basic<T>((x[0]-a[0])/(b[0]-a[0]));
435 const point_basic<T>& x,
436 const point_basic<T>& a,
437 const point_basic<T>& b,
438 const point_basic<T>& c)
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]);
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);
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)
459 for (
size_t i = 0; i < 3; i++) {
465 tensor_basic<T> inv_A;
467 check_macro(!is_singular,
"inv_piola: singular transformation in a tetrahedron");
468 point_basic<T> hat_x = inv_A*ax;
471template<
class T,
class M>
476 const std::vector<size_t>& dis_inod,
479 check_macro (omega.order() == 1,
"inverse piola: mesh order > 1: not yet");
480 if (omega.order() == 1) {
483 omega.dis_node(dis_inod [1]));
485 omega.dis_node(dis_inod [1]),
486 omega.dis_node(dis_inod [2]));
488 omega.dis_node(dis_inod [1]),
489 omega.dis_node(dis_inod [2]),
490 omega.dis_node(dis_inod [3]));
495 F.
reset (omega, hat_K, dis_inod);
498 size_t max_iter = 500, n_iter = max_iter;
499 T tol = std::numeric_limits<Float>::epsilon(), r = tol;
501 check_macro (status == 0,
"inv_piola: newton failed (residue="<<r<<
", n_iter="<<n_iter<<
")");
507#define _RHEOLEF_instanciation1(T) \
510det_jacobian_piola_transformation ( \
511 const tensor_basic<T>& DF, \
516pseudo_inverse_jacobian_piola_transformation ( \
517 const tensor_basic<T>& DF, \
522weight_coordinate_system ( \
523 space_constant::coordinate_type sys_coord, \
524 const point_basic<T>& xq); \
528 const tensor_basic<T>& DF, \
531 tensor_basic<T>& P); \
534#define _RHEOLEF_instanciation2(T,M) \
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); \
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); \
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); \
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); \
569normal_from_piola_transformation ( \
570 const geo_basic<T,M>& omega, \
571 const geo_element& S, \
572 const tensor_basic<T>& DF, \
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); \
589#ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
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
size_type nnod(reference_element hat_K) const
const quadrature< T > & get_quadrature() const
generic mesh with rerefence counting
see the geo_element page for the full documentation
void set_x(const value_type &x1)
void reset(const geo_basic< T, M > &omega, const reference_element &hat_K, const std::vector< size_t > &dis_inod)
value_type initial() const
rep::const_iterator const_iterator
const_iterator end(reference_element hat_K) const
const_iterator begin(reference_element hat_K) const
see the reference_element page for the full documentation
static const variant_type e
size_type dimension() const
variant_type variant() const
static const variant_type T
static const variant_type t
point_basic< T > col(size_type i) const
T determinant(size_type d=3) const
void fill(const T &init_val)
void set_row(const point_basic< T > &r, size_t i, size_t d=3)
#define error_macro(message)
#define fatal_macro(message)
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)
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)
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)
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
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)
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)
T norm(const vec< T, M > &x)
norm(x): see the expression page for the full documentation
bool invert_3x3(const tensor_basic< T > &A, tensor_basic< T > &result)
#define _RHEOLEF_instanciation1(T)
#define _RHEOLEF_instanciation2(T, M)