12 #ifndef MFEM_DINVARIANTS_HPP
13 #define MFEM_DINVARIANTS_HPP
15 #include "../config/config.hpp"
16 #include "../general/cuda.hpp"
33 const double * J_ =
nullptr;
34 double * dI1_ =
nullptr;
35 double * dI1b_ =
nullptr;
36 double * ddI1_ =
nullptr;
37 double * ddI1b_ =
nullptr;
38 double * dI2_ =
nullptr;
39 double * dI2b_ =
nullptr;
40 double * ddI2_ =
nullptr;
41 double * ddI2b_ =
nullptr;
44 MFEM_HOST_DEVICE
Buffers &
J(
const double *
b) { J_ =
b;
return *
this; }
45 MFEM_HOST_DEVICE
Buffers &
dI1(
double *
b) { dI1_ =
b;
return *
this; }
49 MFEM_HOST_DEVICE
Buffers &
dI2(
double *
b) { dI2_ =
b;
return *
this; }
56 double const *
const J;
57 double *
const dI1, *
const dI1b, *
const ddI1, *
const ddI1b;
58 double *
const dI2, *
const dI2b, *
const ddI2, *
const ddI2b;
64 dI1(b.dI1_), dI1b(b.dI1b_), ddI1(b.ddI1_), ddI1b(b.ddI1b_),
65 dI2(b.dI2_), dI2b(b.dI2b_), ddI2(b.ddI2_), ddI2b(b.ddI2b_) { }
67 MFEM_HOST_DEVICE
inline double Get_I2b(
double &sign_detJ)
69 const double I2b = J[0]*J[3] - J[1]*J[2];
70 sign_detJ = I2b >= 0.0 ? 1.0 : -1.0;
71 return sign_detJ * I2b;
80 MFEM_HOST_DEVICE
inline double Get_I2()
86 MFEM_HOST_DEVICE
inline double Get_I1()
88 return J[0]*J[0] + J[1]*J[1] + J[2]*J[2] + J[3]*J[3];
98 dI1[0] = 2*J[0]; dI1[2] = 2*J[2];
99 dI1[1] = 2*J[1]; dI1[3] = 2*J[3];
107 const double c1 = 2.0/
Get_I2b();
108 const double c2 =
Get_I1b()/2.0;
110 dI1b[0] = c1*(J[0] - c2*dI2b[0]);
111 dI1b[1] = c1*(J[1] - c2*dI2b[1]);
112 dI1b[2] = c1*(J[2] - c2*dI2b[2]);
113 dI1b[3] = c1*(J[3] - c2*dI2b[3]);
136 dI2b[0] = sign_detJ*J[3];
137 dI2b[1] = -sign_detJ*J[2];
138 dI2b[2] = -sign_detJ*J[1];
139 dI2b[3] = sign_detJ*J[0];
144 MFEM_HOST_DEVICE
inline double *
Get_ddI1(
int i,
int j)
148 for (
int k=0; k<2; k++)
150 for (
int l=0; l<2; l++)
152 ddi1(k,l) = (i==k && j==l) ? 2.0 : 0.0;
164 double X1_p[4], X2_p[4], X3_p[4];
167 const double I2 =
Get_I2();
170 const double alpha = I1b / I2;
172 for (
int k=0; k<2; k++)
174 for (
int l=0; l<2; l++)
176 X1(k,l) = alpha * (di2b(i,j)*di2b(k,l) + di2b(k,j)*di2b(i,l));
181 const double beta = 1.0 /
Get_I2b();
183 for (
int k=0; k<2; k++)
185 for (
int l=0; l<2; l++)
187 X2(k,l) = beta * ddi1(k,l);
192 const double gamma = -2.0/
Get_I2();
194 for (
int k=0; k<2; k++)
196 for (
int l=0; l<2; l++)
198 X3(k,l) = gamma * (Jpt(i,j)*di2b(k,l) + di2b(i,j)*Jpt(k,l));
202 for (
int k=0; k<2; k++)
204 for (
int l=0; l<2; l++)
206 ddi1b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
213 MFEM_HOST_DEVICE
inline double *
Get_ddI2(
int i,
int j)
217 for (
int k=0; k<2; k++)
219 for (
int l=0; l<2; l++)
221 ddi2(k,l) = 2*di2b(i,j)*di2b(k,l)
222 + 2*(di2b(i,j)*di2b(k,l) - di2b(k,j)*di2b(i,l));
234 for (
int k=0; k<2; k++)
236 for (
int l=0; l<2; l++)
239 for (
int s=0;
s<2;
s++)
241 for (
int t=0;
t<2;
t++)
243 const double ks_it = k==
s && i==
t ? 1.0 : 0.0;
244 const double kt_si = k==
t &&
s==i ? 1.0 : 0.0;
245 ddi2b(k,l) += alpha * (ks_it - kt_si) * di2b(
t,j) * di2b(s,l);
262 const double * J_ =
nullptr;
263 double * B_ =
nullptr;
264 double * dI1_ =
nullptr;
265 double * dI1b_ =
nullptr;
266 double * ddI1_ =
nullptr;
267 double * ddI1b_ =
nullptr;
268 double * dI2_ =
nullptr;
269 double * dI2b_ =
nullptr;
270 double * ddI2_ =
nullptr;
271 double * ddI2b_ =
nullptr;
272 double * dI3b_ =
nullptr;
273 double * ddI3b_ =
nullptr;
276 MFEM_HOST_DEVICE
Buffers &
J(
const double *
b) { J_ =
b;
return *
this; }
277 MFEM_HOST_DEVICE
Buffers &
B(
double *
b) { B_ =
b;
return *
this; }
291 double const *
const J;
293 double *
const dI1, *
const dI1b, *
const ddI1, *
const ddI1b;
294 double *
const dI2, *
const dI2b, *
const ddI2, *
const ddI2b;
295 double *
const dI3b, *
const ddI3b;
301 dI1(b.dI1_), dI1b(b.dI1b_), ddI1(b.ddI1_), ddI1b(b.ddI1b_),
302 dI2(b.dI2_), dI2b(b.dI2b_), ddI2(b.ddI2_), ddI2b(b.ddI2b_),
303 dI3b(b.dI3b_), ddI3b(b.ddI3b_) { }
305 MFEM_HOST_DEVICE
inline double Get_I3b(
double &sign_detJ)
307 const double I3b = + J[0]*(J[4]*J[8] - J[7]*J[5])
308 - J[1]*(J[3]*J[8] - J[5]*J[6])
309 + J[2]*(J[3]*J[7] - J[4]*J[6]);
310 sign_detJ = I3b >= 0.0 ? 1.0 : -1.0;
311 return sign_detJ * I3b;
316 const double I3b = + J[0]*(J[4]*J[8] - J[7]*J[5])
317 - J[1]*(J[3]*J[8] - J[5]*J[6])
318 + J[2]*(J[3]*J[7] - J[4]*J[6]);
331 const double i3b =
Get_I3b(sign_detJ);
332 return sign_detJ *
std::pow(i3b, -2./3.);
335 MFEM_HOST_DEVICE
inline double Get_I3b_p(
double &sign_detJ)
337 const double i3b =
Get_I3b(sign_detJ);
338 return sign_detJ *
std::pow(i3b, -2./3.);
343 B[0] = J[0]*J[0] + J[3]*J[3] + J[6]*J[6];
344 B[1] = J[1]*J[1] + J[4]*J[4] + J[7]*J[7];
345 B[2] = J[2]*J[2] + J[5]*J[5] + J[8]*J[8];
346 const double I1 = B[0] + B[1] + B[2];
350 MFEM_HOST_DEVICE
inline
361 B[3] = J[0]*J[1] + J[3]*J[4] + J[6]*J[7];
362 B[4] = J[0]*J[2] + J[3]*J[5] + J[6]*J[8];
363 B[5] = J[1]*J[2] + J[4]*J[5] + J[7]*J[8];
369 const double I1 =
Get_I1();
370 const double BF2 = B[0]*B[0] + B[1]*B[1] + B[2]*B[2] +
371 2*(B[3]*B[3] + B[4]*B[4] + B[5]*B[5]);
372 const double I2 = (I1*I1 - BF2)/2;
379 return Get_I2() * I3b_p * I3b_p;
384 for (
int i = 0; i < 9; i++) { dI1[i] = 2*J[i]; }
393 const double I3b =
Get_I3b(sign_detJ);
395 const double c1 = 2.0 * I3b_p;
396 const double c2 =
Get_I1()/(3.0 * I3b);
398 for (
int i = 0; i < 9; i++) { dI1b[i] = c1*(J[i] - c2*dI3b[i]); }
405 const double I1 =
Get_I1();
411 2*(I1 - B[0]), 2*(I1 - B[1]), 2*(I1 - B[2]),
412 -2*B[3], -2*B[4], -2*B[5]
417 dI2[0] = C[0]*J[0] + C[3]*J[1] + C[4]*J[2];
418 dI2[1] = C[3]*J[0] + C[1]*J[1] + C[5]*J[2];
419 dI2[2] = C[4]*J[0] + C[5]*J[1] + C[2]*J[2];
421 dI2[3] = C[0]*J[3] + C[3]*J[4] + C[4]*J[5];
422 dI2[4] = C[3]*J[3] + C[1]*J[4] + C[5]*J[5];
423 dI2[5] = C[4]*J[3] + C[5]*J[4] + C[2]*J[5];
425 dI2[6] = C[0]*J[6] + C[3]*J[7] + C[4]*J[8];
426 dI2[7] = C[3]*J[6] + C[1]*J[7] + C[5]*J[8];
427 dI2[8] = C[4]*J[6] + C[5]*J[7] + C[2]*J[8];
437 const double I2 =
Get_I2();
439 const double I3b =
Get_I3b(sign_detJ);
440 const double c1 = I3b_p*I3b_p;
441 const double c2 = (4*I2/I3b)/3;
444 for (
int i = 0; i < 9; i++) { dI2b[i] = c1*(dI2[i] - c2*dI3b[i]); }
448 MFEM_HOST_DEVICE
inline double *
Get_dI3b(
const double sign_detJ)
452 dI3b[0] = sign_detJ*(J[4]*J[8] - J[5]*J[7]);
453 dI3b[1] = sign_detJ*(J[5]*J[6] - J[3]*J[8]);
454 dI3b[2] = sign_detJ*(J[3]*J[7] - J[4]*J[6]);
455 dI3b[3] = sign_detJ*(J[2]*J[7] - J[1]*J[8]);
456 dI3b[4] = sign_detJ*(J[0]*J[8] - J[2]*J[6]);
457 dI3b[5] = sign_detJ*(J[1]*J[6] - J[0]*J[7]);
458 dI3b[6] = sign_detJ*(J[1]*J[5] - J[2]*J[4]);
459 dI3b[7] = sign_detJ*(J[2]*J[3] - J[0]*J[5]);
460 dI3b[8] = sign_detJ*(J[0]*J[4] - J[1]*J[3]);
465 MFEM_HOST_DEVICE
inline double *
Get_ddI1(
int i,
int j)
468 for (
int k=0; k<3; k++)
470 for (
int l=0; l<3; l++)
472 const double I_ijkl = (i==k && j==l) ? 1.0 : 0.0;
473 ddi1(k,l) = 2.0 * I_ijkl;
488 double X1_p[9], X2_p[9], X3_p[9];
490 const double I3 =
Get_I3();
492 const double alpha = (2./3.)*I1b/I3;
494 for (
int k=0; k<3; k++)
496 for (
int l=0; l<3; l++)
498 X1(k,l) = alpha * ((2./3.)*di3b(i,j) * di3b(k,l) +
499 di3b(k,j)*di3b(i,l));
506 for (
int k=0; k<3; k++)
508 for (
int l=0; l<3; l++)
510 const double ddI1_ijkl = (i==k && j==l) ? 2.0 : 0.0;
511 X2(k,l) = beta * ddI1_ijkl;
517 const double gamma = -(4./3.)*
Get_I3b_p()/I3b;
519 for (
int k=0; k<3; k++)
521 for (
int l=0; l<3; l++)
523 X3(k,l) = gamma * (Jpt(i,j) * di3b(k,l) + di3b(i,j) * Jpt(k,l));
527 for (
int k=0; k<3; k++)
529 for (
int l=0; l<3; l++)
531 ddi1b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
541 MFEM_HOST_DEVICE
inline double *
Get_ddI2(
int i,
int j)
543 double x1_p[9], x2_p[9], x3_p[9];
546 const double I1 =
Get_I1();
547 for (
int k=0; k<3; k++)
549 for (
int l=0; l<3; l++)
551 const double ik_jl = (i==k && j==l) ? 1.0 : 0.0;
552 x1(k,l) = 2.0 * I1 * ik_jl;
557 for (
int k=0; k<3; k++)
559 for (
int l=0; l<3; l++)
562 for (
int u=0;
u<3;
u++)
564 for (
int v=0; v<3; v++)
566 const double ku_iv = k==
u && i==v ? 1.0 : 0.0;
567 const double ik_uv = i==k &&
u==v ? 1.0 : 0.0;
568 const double kv_iu = k==v && i==u ? 1.0 : 0.0;
569 x2(k,l) += 2.0*(2.*ku_iv-ik_uv-kv_iu)*Jpt(v,j)*Jpt(u,l);
575 B[0] = J[0]*J[0] + J[3]*J[3] + J[6]*J[6];
576 B[1] = J[1]*J[1] + J[4]*J[4] + J[7]*J[7];
577 B[2] = J[2]*J[2] + J[5]*J[5] + J[8]*J[8];
578 B[3] = J[0]*J[1] + J[3]*J[4] + J[6]*J[7];
579 B[4] = J[0]*J[2] + J[3]*J[5] + J[6]*J[8];
580 B[5] = J[1]*J[2] + J[4]*J[5] + J[7]*J[8];
581 const double b_p[9] =
588 for (
int k=0; k<3; k++)
590 for (
int l=0; l<3; l++)
592 const double jl = j==l ? 1.0 : 0.0;
593 x3(k,l) = -2.0 *
b(i,k) * jl;
598 for (
int k=0; k<3; k++)
600 for (
int l=0; l<3; l++)
602 ddi2(k,l) = x1(k,l) + x2(k,l) + x3(k,l);
615 double X1_p[9], X2_p[9], X3_p[9];
621 const double I3b =
Get_I3b(sign_detJ);
622 const double I2 =
Get_I2();
623 const double I3b_p43 = I3b_p*I3b_p;
624 const double I3b_p73 = I3b_p*I3b_p/I3b;
625 const double I3b_p103 = I3b_p*I3b_p/(I3b*I3b);
627 for (
int k=0; k<3; k++)
629 for (
int l=0; l<3; l++)
631 const double up = (16./9.)*I3b_p103*I2*di3b(i,j)*di3b(k,l);
632 const double down = (4./3.)*I3b_p103*I2*di3b(i,l)*di3b(k,j);
639 for (
int k=0; k<3; k++)
641 for (
int l=0; l<3; l++)
643 X2(k,l) = -(4./3.)*I3b_p73*(di2(i,j)*di3b(k,l)+di2(k,l)*di3b(i,j));
649 for (
int k=0; k<3; k++)
651 for (
int l=0; l<3; l++)
653 X3(k,l) = I3b_p43 * ddi2(k,l);
658 for (
int k=0; k<3; k++)
660 for (
int l=0; l<3; l++)
662 ddi2b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
672 const double c1 = 1./
Get_I3b();
675 for (
int k=0; k<3; k++)
677 for (
int l=0; l<3; l++)
680 for (
int s=0;
s<3;
s++)
682 for (
int t=0;
t<3;
t++)
684 const double ks_it = k==
s && i==
t ? 1.0 : 0.0;
685 const double kt_si = k==
t &&
s==i ? 1.0 : 0.0;
686 ddi3b(k,l) += c1*(ks_it-kt_si)*di3b(
t,j)*di3b(s,l);
699 #endif // MFEM_DINVARIANTS_HPP
MFEM_HOST_DEVICE Buffers & dI1(double *b)
MFEM_HOST_DEVICE double Get_I3()
MFEM_HOST_DEVICE double * Get_ddI1b(int i, int j)
MFEM_HOST_DEVICE Buffers & ddI2(double *b)
MFEM_HOST_DEVICE Buffers & ddI2b(double *b)
MFEM_HOST_DEVICE double Get_I2()
MFEM_HOST_DEVICE Buffers & B(double *b)
MFEM_HOST_DEVICE Buffers & dI2(double *b)
MFEM_HOST_DEVICE Buffers & dI2b(double *b)
MFEM_HOST_DEVICE Buffers & ddI2b(double *b)
MFEM_HOST_DEVICE double * Get_dI2()
MFEM_HOST_DEVICE double * Get_dI2()
MFEM_HOST_DEVICE double * Get_dI1()
MFEM_HOST_DEVICE Buffers & ddI1b(double *b)
MFEM_HOST_DEVICE Buffers & ddI3b(double *b)
MFEM_HOST_DEVICE Buffers()
MFEM_HOST_DEVICE double Get_I3b_p()
MFEM_HOST_DEVICE Buffers & J(const double *b)
MFEM_HOST_DEVICE double * Get_ddI1(int i, int j)
MFEM_HOST_DEVICE double * Get_ddI1b(int i, int j)
MFEM_HOST_DEVICE Buffers & dI3b(double *b)
MFEM_HOST_DEVICE double * Get_ddI1(int i, int j)
MFEM_HOST_DEVICE double Get_I3b(double &sign_detJ)
MFEM_HOST_DEVICE double Get_I3b()
MFEM_HOST_DEVICE double * Get_dI3b(const double sign_detJ)
A basic generic Tensor class, appropriate for use on the GPU.
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
MFEM_HOST_DEVICE double * Get_ddI3b(int i, int j)
MFEM_HOST_DEVICE double * Get_ddI2(int i, int j)
MFEM_HOST_DEVICE double Get_I1b()
MFEM_HOST_DEVICE void Get_B_offd()
MFEM_HOST_DEVICE double Get_I2b(double &sign_detJ)
MFEM_HOST_DEVICE double * Get_dI1b()
MFEM_HOST_DEVICE double Get_I1()
MFEM_HOST_DEVICE Buffers & ddI2(double *b)
MFEM_HOST_DEVICE double * Get_ddI2b(int i, int j)
MFEM_HOST_DEVICE double * Get_dI1b()
MFEM_HOST_DEVICE Buffers & J(const double *b)
MFEM_HOST_DEVICE Buffers & dI1(double *b)
MFEM_HOST_DEVICE double Get_I2b()
MFEM_HOST_DEVICE InvariantsEvaluator3D(Buffers &b)
MFEM_HOST_DEVICE InvariantsEvaluator2D(Buffers &b)
MFEM_HOST_DEVICE Buffers & dI2b(double *b)
MFEM_HOST_DEVICE double * Get_ddI2b(int i, int j)
MFEM_HOST_DEVICE Buffers & dI2(double *b)
MFEM_HOST_DEVICE double * Get_ddI2(int i, int j)
MFEM_HOST_DEVICE Buffers & ddI1b(double *b)
MFEM_HOST_DEVICE double Get_I2()
MFEM_HOST_DEVICE double Get_I3b_p(double &sign_detJ)
double u(const Vector &xvec)
MFEM_HOST_DEVICE Buffers & ddI1(double *b)
MFEM_HOST_DEVICE double * Get_dI2b()
MFEM_HOST_DEVICE double Get_I1b()
MFEM_HOST_DEVICE double * Get_dI2b()
MFEM_HOST_DEVICE Buffers & dI1b(double *b)
MFEM_HOST_DEVICE double Get_I1()
MFEM_HOST_DEVICE Buffers & dI1b(double *b)
MFEM_HOST_DEVICE double * Get_dI1()
MFEM_HOST_DEVICE Buffers & ddI1(double *b)
MFEM_HOST_DEVICE double Get_I2b()
MFEM_HOST_DEVICE Buffers()