1 #include "rheolef/piola_util.h"
22 #include "rheolef/geo_domain.h"
23 #include "rheolef/inv_piola.h"
24 #include "rheolef/damped_newton.h"
41 template<
class T,
class M>
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) {
64 x[loc_inod][
alpha] += phij_xi (loc_inod,loc_jdof)*xjnod[
alpha];
78 template<
class T,
class M>
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) {
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);
114 template<
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++) {
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.");
162 template<
class T,
class M>
173 if (S.
dimension() + 1 == omega.map_dimension()) {
175 check_macro (dis_ie != std::numeric_limits<size_type>::max(),
"normal: requires neighbours initialization");
177 Float sign = (S[0] == K[1]) ? 1 : -1;
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");
188 if (omega.get_background_geo().map_dimension() == 1) {
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());
195 Float sign = (bgd_S[0] == bgd_K[1]) ? 1 : -1;
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);
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());
208 Float sign = (bgd2_S[0] == bgd2_K[1]) ? 1 : -1;
213 size_type first_dis_isid = omega.sizes().ownership_by_dimension[S.
dimension()].first_index();
214 size_type isid = dis_isid - first_dis_isid;
215 check_macro (dis_isid >= first_dis_isid,
"unexpected dis_index "<<dis_isid<<
": out of local range");
220 check_macro (bgd_dis_ie != std::numeric_limits<size_type>::max(),
221 "normal: bgd_S.dis_ie={"<<bgd_S.
dis_ie()<<
"} without neighbours; requires neighbours initialization for mesh " << bgd_omega.name());
223 Float sign = (bgd_S[0] == bgd_K[1]) ? 1 : -1;
271 return inv(DF, map_d);
296 error_macro (
"pseudo_inverse_jacobian_piola_transformation: unsupported element dimension "
297 << map_d <<
" in " <<
d <<
"D mesh.");
324 template<
class T,
class M>
330 const std::vector<size_t>&
dis_inod,
331 bool ignore_sys_coord,
334 Eigen::Matrix<T,Eigen::Dynamic,1>& w)
345 for (
size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
346 w[loc_inod] = x[loc_inod][k];
353 first_quad = quad.
begin(hat_K),
354 last_quad = quad.
end (hat_K);
355 for (
size_type q = 0; first_quad != last_quad; ++first_quad, ++q) {
357 T wq = det_DF*(*first_quad).w;
358 if (! ignore_sys_coord) {
378 check_macro (
d == map_d+1,
"unexpected dimension map_d="<<map_d<<
" and d="<<
d);
392 for (
size_t l = 0; l <
d; l++) {
393 for (
size_t m = 0;
m <
d;
m++) {
394 P(l,
m) = -
n[l]*
n[
m];
426 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]);
430 (
b[1]-
a[1])*t9*t11-(
b[0]-
a[0])*t9*t15);
445 for (
size_t i = 0; i < 3; i++) {
451 tensor_basic<T> inv_A;
453 check_macro(!is_singular,
"inv_piola: singular transformation in a tetrahedron");
457 template<
class T,
class M>
462 const std::vector<size_t>&
dis_inod,
465 check_macro (omega.order() == 1,
"inverse piola: mesh order > 1: not yet");
466 if (omega.order() == 1) {
484 size_t max_iter = 500, n_iter = max_iter;
487 check_macro (
status == 0,
"inv_piola: newton failed (residue="<<r<<
", n_iter="<<n_iter<<
")");
493 #define _RHEOLEF_instanciation1(T) \
496 det_jacobian_piola_transformation ( \
497 const tensor_basic<T>& DF, \
502 pseudo_inverse_jacobian_piola_transformation ( \
503 const tensor_basic<T>& DF, \
508 weight_coordinate_system ( \
509 space_constant::coordinate_type sys_coord, \
510 const point_basic<T>& xq); \
514 const tensor_basic<T>& DF, \
517 tensor_basic<T>& P); \
520 #define _RHEOLEF_instanciation2(T,M) \
523 piola_transformation ( \
524 const geo_basic<T,M>& omega, \
525 const basis_on_pointset<T>& piola_on_pointset, \
526 reference_element hat_K, \
527 const std::vector<size_t>& dis_inod, \
528 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x); \
531 inverse_piola_transformation ( \
532 const geo_basic<T,M>& omega, \
533 const reference_element& hat_K, \
534 const std::vector<size_t>& dis_inod, \
535 const point_basic<T>& x); \
538 jacobian_piola_transformation ( \
539 const geo_basic<T,M>& omega, \
540 const basis_basic<T>& piola_basis, \
541 reference_element hat_K, \
542 const std::vector<size_t>& dis_inod, \
543 const point_basic<T>& hat_x, \
544 tensor_basic<T>& DF); \
547 jacobian_piola_transformation ( \
548 const geo_basic<T,M>& omega, \
549 const basis_on_pointset<T>& piola_on_pointset, \
550 reference_element hat_K, \
551 const std::vector<size_t>& dis_inod, \
552 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& DF); \
555 normal_from_piola_transformation ( \
556 const geo_basic<T,M>& omega, \
557 const geo_element& S, \
558 const tensor_basic<T>& DF, \
562 piola_transformation_and_weight_integration ( \
563 const geo_basic<T,M>& omega, \
564 const basis_on_pointset<T>& piola_on_pointset, \
565 reference_element hat_K, \
566 const std::vector<size_t>& dis_inod, \
567 bool ignore_sys_coord, \
568 Eigen::Matrix<tensor_basic<T>,Eigen::Dynamic,1>& DF, \
569 Eigen::Matrix<point_basic<T>,Eigen::Dynamic,1>& x, \
570 Eigen::Matrix<T,Eigen::Dynamic,1>& w); \
575 #ifdef _RHEOLEF_HAVE_MPI
577 #endif // _RHEOLEF_HAVE_MPI