MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
dinvariants.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_DINVARIANTS_HPP
13 #define MFEM_DINVARIANTS_HPP
14 
15 #include "../config/config.hpp"
16 #include "../general/cuda.hpp"
17 #include "dtensor.hpp"
18 #include <cmath>
19 
20 namespace mfem
21 {
22 
23 namespace kernels
24 {
25 
27 {
28 public:
29  class Buffers
30  {
31  friend class InvariantsEvaluator2D;
32  private:
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;
42  public:
43  MFEM_HOST_DEVICE Buffers() {}
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; }
46  MFEM_HOST_DEVICE Buffers &dI1b(double *b) { dI1b_ = b; return *this; }
47  MFEM_HOST_DEVICE Buffers &ddI1(double *b) { ddI1_ = b; return *this; }
48  MFEM_HOST_DEVICE Buffers &ddI1b(double *b) { ddI1b_ = b; return *this; }
49  MFEM_HOST_DEVICE Buffers &dI2(double *b) { dI2_ = b; return *this; }
50  MFEM_HOST_DEVICE Buffers &dI2b(double *b) { dI2b_ = b; return *this; }
51  MFEM_HOST_DEVICE Buffers &ddI2(double *b) { ddI2_ = b; return *this; }
52  MFEM_HOST_DEVICE Buffers &ddI2b(double *b) { ddI2b_ = b; return *this; }
53  };
54 
55 private:
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;
59 
60 public:
61  MFEM_HOST_DEVICE
63  J(b.J_),
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_) { }
66 
67  MFEM_HOST_DEVICE inline double Get_I2b(double &sign_detJ) // det(J) + sign
68  {
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;
72  }
73 
74  MFEM_HOST_DEVICE inline double Get_I2b() // det(J)
75  {
76  double sign_detJ;
77  return Get_I2b(sign_detJ);
78  }
79 
80  MFEM_HOST_DEVICE inline double Get_I2() // det(J)^{2}
81  {
82  const double I2b = Get_I2b();
83  return I2b * I2b;
84  }
85 
86  MFEM_HOST_DEVICE inline double Get_I1() // I1 = ||J||_F^2
87  {
88  return J[0]*J[0] + J[1]*J[1] + J[2]*J[2] + J[3]*J[3];
89  }
90 
91  MFEM_HOST_DEVICE inline double Get_I1b() // I1b = I1/det(J)
92  {
93  return Get_I1() / Get_I2b();
94  }
95 
96  MFEM_HOST_DEVICE inline double *Get_dI1()
97  {
98  dI1[0] = 2*J[0]; dI1[2] = 2*J[2];
99  dI1[1] = 2*J[1]; dI1[3] = 2*J[3];
100  return dI1;
101  }
102 
103  MFEM_HOST_DEVICE inline double *Get_dI1b()
104  {
105  // I1b = I1/I2b
106  // dI1b = (1/I2b)*dI1 - (I1/I2b^2)*dI2b = (2/I2b)*[J - (I1b/2)*dI2b]
107  const double c1 = 2.0/Get_I2b();
108  const double c2 = Get_I1b()/2.0;
109  Get_dI2b();
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]);
114  return dI1b;
115  }
116 
117  MFEM_HOST_DEVICE inline double *Get_dI2()
118  {
119  // I2 = I2b^2
120  // dI2 = 2*I2b*dI2b = 2*det(J)*adj(J)^T
121  const double c1 = 2*Get_I2b();
122  Get_dI2b();
123  dI2[0] = c1*dI2b[0];
124  dI2[1] = c1*dI2b[1];
125  dI2[2] = c1*dI2b[2];
126  dI2[3] = c1*dI2b[3];
127  return dI2;
128  }
129 
130  MFEM_HOST_DEVICE inline double *Get_dI2b()
131  {
132  // I2b = det(J)
133  // dI2b = adj(J)^T
134  double sign_detJ;
135  Get_I2b(sign_detJ);
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];
140  return dI2b;
141  }
142 
143  // ddI1_ijkl = 2 I_ijkl = 2 δ_ik δ_jl
144  MFEM_HOST_DEVICE inline double *Get_ddI1(int i, int j)
145  {
146  // ddI1_ijkl = 2 I_ijkl = 2 δ_ik δ_jl
147  DeviceMatrix ddi1(ddI1,2,2);
148  for (int k=0; k<2; k++)
149  {
150  for (int l=0; l<2; l++)
151  {
152  ddi1(k,l) = (i==k && j==l) ? 2.0 : 0.0;
153  }
154  }
155  return ddI1;
156  }
157 
158  // ddI1b = X1 + X2 + X3, where
159  // X1_ijkl = (I1b/I2) [ dI2b_ij dI2b_kl + dI2b_kj dI2b_il ]
160  // X2_ijkl = (1/I2b) ddI1_ijkl
161  // X3_ijkl = -(2/I2) (J_ij dI2b_kl + dI2b_ij J_kl)
162  MFEM_HOST_DEVICE inline double *Get_ddI1b(int i, int j)
163  {
164  double X1_p[4], X2_p[4], X3_p[4];
165 
166  // X1_ijkl = (I1b/I2) [ dI2b_ij dI2b_kl + dI2b_kj dI2b_il ]
167  const double I2 = Get_I2();
168  const double I1b = Get_I1b();
169  ConstDeviceMatrix di2b(Get_dI2b(),2,2);
170  const double alpha = I1b / I2;
171  DeviceMatrix X1(X1_p,2,2);
172  for (int k=0; k<2; k++)
173  {
174  for (int l=0; l<2; l++)
175  {
176  X1(k,l) = alpha * (di2b(i,j)*di2b(k,l) + di2b(k,j)*di2b(i,l));
177  }
178  }
179  // X2_ijkl = (1/I2b) ddI1_ijkl
180  DeviceMatrix X2(X2_p,2,2);
181  const double beta = 1.0 / Get_I2b();
182  ConstDeviceMatrix ddi1(Get_ddI1(i,j),2,2);
183  for (int k=0; k<2; k++)
184  {
185  for (int l=0; l<2; l++)
186  {
187  X2(k,l) = beta * ddi1(k,l);
188  }
189  }
190  // X3_ijkl = -(2/I2) (J_ij dI2b_kl + dI2b_ij J_kl)
191  DeviceMatrix X3(X3_p,2,2);
192  const double gamma = -2.0/Get_I2();
193  ConstDeviceMatrix Jpt(J,2,2);
194  for (int k=0; k<2; k++)
195  {
196  for (int l=0; l<2; l++)
197  {
198  X3(k,l) = gamma * (Jpt(i,j)*di2b(k,l) + di2b(i,j)*Jpt(k,l));
199  }
200  }
201  DeviceMatrix ddi1b(ddI1b,2,2);
202  for (int k=0; k<2; k++)
203  {
204  for (int l=0; l<2; l++)
205  {
206  ddi1b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
207  }
208  }
209  return ddI1b;
210  }
211 
212  // ddI2_ijkl = 2 dI2b_ij dI2b_kl + 2 (dI2b_ij dI2b_kl - dI2b_kj dI2b_il)
213  MFEM_HOST_DEVICE inline double *Get_ddI2(int i, int j)
214  {
215  DeviceMatrix ddi2(ddI2,2,2);
216  ConstDeviceMatrix di2b(Get_dI2b(),2,2);
217  for (int k=0; k<2; k++)
218  {
219  for (int l=0; l<2; l++)
220  {
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));
223  }
224  }
225  return ddI2;
226  }
227 
228  // ddI2b_ijkl = (1/I2b) (δ_ks δ_it - δ_kt δ_si) dI2b_tj dI2b_sl
229  MFEM_HOST_DEVICE inline double *Get_ddI2b(int i, int j)
230  {
231  DeviceMatrix ddi2b(ddI2b,2,2);
232  const double alpha = 1.0/Get_I2b();
233  ConstDeviceMatrix di2b(Get_dI2b(),2,2);
234  for (int k=0; k<2; k++)
235  {
236  for (int l=0; l<2; l++)
237  {
238  ddi2b(k,l) = 0.0;
239  for (int s=0; s<2; s++)
240  {
241  for (int t=0; t<2; t++)
242  {
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);
246  }
247  }
248  }
249  }
250  return ddI2b;
251  }
252 };
253 
254 
256 {
257 public:
258  class Buffers
259  {
260  friend class InvariantsEvaluator3D;
261  private:
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;
274  public:
275  MFEM_HOST_DEVICE Buffers() {}
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; }
278  MFEM_HOST_DEVICE Buffers &dI1(double *b) { dI1_ = b; return *this; }
279  MFEM_HOST_DEVICE Buffers &dI1b(double *b) { dI1b_ = b; return *this; }
280  MFEM_HOST_DEVICE Buffers &ddI1(double *b) { ddI1_ = b; return *this; }
281  MFEM_HOST_DEVICE Buffers &ddI1b(double *b) { ddI1b_ = b; return *this; }
282  MFEM_HOST_DEVICE Buffers &dI2(double *b) { dI2_ = b; return *this; }
283  MFEM_HOST_DEVICE Buffers &dI2b(double *b) { dI2b_ = b; return *this; }
284  MFEM_HOST_DEVICE Buffers &ddI2(double *b) { ddI2_ = b; return *this; }
285  MFEM_HOST_DEVICE Buffers &ddI2b(double *b) { ddI2b_ = b; return *this; }
286  MFEM_HOST_DEVICE Buffers &dI3b(double *b) { dI3b_ = b; return *this; }
287  MFEM_HOST_DEVICE Buffers &ddI3b(double *b) { ddI3b_ = b; return *this; }
288  };
289 
290 private:
291  double const * const J;
292  double * const B;
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;
296 
297 public:
298  MFEM_HOST_DEVICE
300  J(b.J_), B(b.B_),
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_) { }
304 
305  MFEM_HOST_DEVICE inline double Get_I3b(double &sign_detJ) // det(J) + sign
306  {
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;
312  }
313 
314  MFEM_HOST_DEVICE inline double Get_I3b() // det(J)
315  {
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]);
319  return I3b;
320  }
321 
322  MFEM_HOST_DEVICE inline double Get_I3() // det(J)^{2}
323  {
324  const double I3b = Get_I3b();
325  return I3b * I3b;
326  }
327 
328  MFEM_HOST_DEVICE inline double Get_I3b_p() // I3b^{-2/3}
329  {
330  double sign_detJ;
331  const double i3b = Get_I3b(sign_detJ);
332  return sign_detJ * std::pow(i3b, -2./3.);
333  }
334 
335  MFEM_HOST_DEVICE inline double Get_I3b_p(double &sign_detJ) // I3b^{-2/3}
336  {
337  const double i3b = Get_I3b(sign_detJ);
338  return sign_detJ * std::pow(i3b, -2./3.);
339  }
340 
341  MFEM_HOST_DEVICE inline double Get_I1()
342  {
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];
347  return I1;
348  }
349 
350  MFEM_HOST_DEVICE inline
351  double Get_I1b() // det(J)^{-2/3}*I_1 = I_1/I_3^{1/3}
352  {
353  const double I1b = Get_I1() * Get_I3b_p();
354  return I1b;
355  }
356 
357  MFEM_HOST_DEVICE inline void Get_B_offd()
358  {
359  // B = J J^t
360  // B[3]=B(0,1), B[4]=B(0,2), B[5]=B(1,2)
361  B[3] = J[0]*J[1] + J[3]*J[4] + J[6]*J[7]; // B(0,1)
362  B[4] = J[0]*J[2] + J[3]*J[5] + J[6]*J[8]; // B(0,2)
363  B[5] = J[1]*J[2] + J[4]*J[5] + J[7]*J[8]; // B(1,2)
364  }
365 
366  MFEM_HOST_DEVICE inline double Get_I2()
367  {
368  Get_B_offd();
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;
373  return I2;
374  }
375 
376  MFEM_HOST_DEVICE inline double Get_I2b() // I2b = I2*I3b^{-4/3}
377  {
378  const double I3b_p = Get_I3b_p();
379  return Get_I2() * I3b_p * I3b_p;
380  }
381 
382  MFEM_HOST_DEVICE inline double *Get_dI1()
383  {
384  for (int i = 0; i < 9; i++) { dI1[i] = 2*J[i]; }
385  return dI1;
386  }
387 
388  MFEM_HOST_DEVICE inline double *Get_dI1b()
389  {
390  // I1b = I3b^{-2/3}*I1
391  // dI1b = 2*I3b^{-2/3}*(J - (1/3)*I1/I3b*dI3b)
392  double sign_detJ;
393  const double I3b = Get_I3b(sign_detJ);
394  const double I3b_p = Get_I3b_p();
395  const double c1 = 2.0 * I3b_p;
396  const double c2 = Get_I1()/(3.0 * I3b);
397  Get_dI3b(sign_detJ);
398  for (int i = 0; i < 9; i++) { dI1b[i] = c1*(J[i] - c2*dI3b[i]); }
399  return dI1b;
400  }
401 
402  MFEM_HOST_DEVICE inline double *Get_dI2()
403  {
404  // dI2 = 2 I_1 J - 2 J J^t J = 2 (I_1 I - B) J
405  const double I1 = Get_I1();
406  Get_B_offd();
407  // B[0]=B(0,0), B[1]=B(1,1), B[2]=B(2,2)
408  // B[3]=B(0,1), B[4]=B(0,2), B[5]=B(1,2)
409  const double C[6] =
410  {
411  2*(I1 - B[0]), 2*(I1 - B[1]), 2*(I1 - B[2]),
412  -2*B[3], -2*B[4], -2*B[5]
413  };
414  // | C[0] C[3] C[4] | | J[0] J[3] J[6] |
415  // dI2 = | C[3] C[1] C[5] | | J[1] J[4] J[7] |
416  // | C[4] C[5] C[2] | | J[2] J[5] J[8] |
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];
420 
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];
424 
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];
428  return dI2;
429  }
430 
431  MFEM_HOST_DEVICE inline double *Get_dI2b()
432  {
433  // I2b = det(J)^{-4/3}*I2 = I3b^{-4/3}*I2
434  // dI2b = (-4/3)*I3b^{-7/3}*I2*dI3b + I3b^{-4/3}*dI2
435  // = I3b^{-4/3} * [ dI2 - (4/3)*I2/I3b*dI3b ]
436  double sign_detJ;
437  const double I2 = Get_I2();
438  const double I3b_p = Get_I3b_p();
439  const double I3b = Get_I3b(sign_detJ);
440  const double c1 = I3b_p*I3b_p;
441  const double c2 = (4*I2/I3b)/3;
442  Get_dI2();
443  Get_dI3b(sign_detJ);
444  for (int i = 0; i < 9; i++) { dI2b[i] = c1*(dI2[i] - c2*dI3b[i]); }
445  return dI2b;
446  }
447 
448  MFEM_HOST_DEVICE inline double *Get_dI3b(const double sign_detJ)
449  {
450  // I3b = det(J)
451  // dI3b = adj(J)^T
452  dI3b[0] = sign_detJ*(J[4]*J[8] - J[5]*J[7]); // 0 3 6
453  dI3b[1] = sign_detJ*(J[5]*J[6] - J[3]*J[8]); // 1 4 7
454  dI3b[2] = sign_detJ*(J[3]*J[7] - J[4]*J[6]); // 2 5 8
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]);
461  return dI3b;
462  }
463 
464  // ddI1_ijkl = 2 I_ijkl = 2 δ_ik δ_jl
465  MFEM_HOST_DEVICE inline double *Get_ddI1(int i, int j)
466  {
467  DeviceMatrix ddi1(ddI1,3,3);
468  for (int k=0; k<3; k++)
469  {
470  for (int l=0; l<3; l++)
471  {
472  const double I_ijkl = (i==k && j==l) ? 1.0 : 0.0;
473  ddi1(k,l) = 2.0 * I_ijkl;
474  }
475  }
476  return ddI1;
477  }
478 
479  // ddI1b = X1 + X2 + X3, where
480  // X1_ijkl = (2/3*I1b/I3) [ 2/3 dI3b_ij dI3b_kl + dI3b_kj dI3b_il ]
481  // X2_ijkl = (I3b^{-2/3}) ddI1_ijkl
482  // X3_ijkl = -(4/3*I3b^{-5/3}) (J_ij dI3b_kl + dI3b_ij J_kl)
483  MFEM_HOST_DEVICE inline double *Get_ddI1b(int i, int j)
484  {
485  // X1_ijkl = (2/3*I1b/I3) [ 2/3 dI3b_ij dI3b_kl + dI3b_kj dI3b_il ]
486  double sign_detJ;
487  Get_I3b(sign_detJ);
488  double X1_p[9], X2_p[9], X3_p[9];
489  DeviceMatrix X1(X1_p,3,3);
490  const double I3 = Get_I3();
491  const double I1b = Get_I1b();
492  const double alpha = (2./3.)*I1b/I3;
493  ConstDeviceMatrix di3b(Get_dI3b(sign_detJ),3,3);
494  for (int k=0; k<3; k++)
495  {
496  for (int l=0; l<3; l++)
497  {
498  X1(k,l) = alpha * ((2./3.)*di3b(i,j) * di3b(k,l) +
499  di3b(k,j)*di3b(i,l));
500  }
501  }
502  // ddI1_ijkl = 2 δ_ik δ_jl
503  // X2_ijkl = (I3b^{-2/3}) ddI1_ijkl
504  DeviceMatrix X2(X2_p,3,3);
505  const double beta = Get_I3b_p();
506  for (int k=0; k<3; k++)
507  {
508  for (int l=0; l<3; l++)
509  {
510  const double ddI1_ijkl = (i==k && j==l) ? 2.0 : 0.0;
511  X2(k,l) = beta * ddI1_ijkl;
512  }
513  }
514  // X3_ijkl = -(4/3*I3b^{-5/3}) (J_ij dI3b_kl + dI3b_ij J_kl)
515  DeviceMatrix X3(X3_p,3,3);
516  const double I3b = Get_I3b();
517  const double gamma = -(4./3.)*Get_I3b_p()/I3b;
518  ConstDeviceMatrix Jpt(J,3,3);
519  for (int k=0; k<3; k++)
520  {
521  for (int l=0; l<3; l++)
522  {
523  X3(k,l) = gamma * (Jpt(i,j) * di3b(k,l) + di3b(i,j) * Jpt(k,l));
524  }
525  }
526  DeviceMatrix ddi1b(ddI1b,3,3);
527  for (int k=0; k<3; k++)
528  {
529  for (int l=0; l<3; l++)
530  {
531  ddi1b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
532  }
533  }
534  return ddI1b;
535  }
536 
537  // ddI2 = x1 + x2 + x3
538  // x1_ijkl = (2 I1) δ_ik δ_jl
539  // x2_ijkl = 2 ( 2 δ_ku δ_iv - δ_ik δ_uv - δ_kv δ_iu ) J_vj J_ul
540  // x3_ijkl = -2 (J J^t)_ik δ_jl = -2 B_ik δ_jl
541  MFEM_HOST_DEVICE inline double *Get_ddI2(int i, int j)
542  {
543  double x1_p[9], x2_p[9], x3_p[9];
544  DeviceMatrix x1(x1_p,3,3), x2(x2_p,3,3), x3(x3_p,3,3);
545  // x1_ijkl = (2 I1) δ_ik δ_jl
546  const double I1 = Get_I1();
547  for (int k=0; k<3; k++)
548  {
549  for (int l=0; l<3; l++)
550  {
551  const double ik_jl = (i==k && j==l) ? 1.0 : 0.0;
552  x1(k,l) = 2.0 * I1 * ik_jl;
553  }
554  }
555  // x2_ijkl = 2 ( 2 δ_ku δ_iv - δ_ik δ_uv - δ_kv δ_iu ) J_vj J_ul
556  ConstDeviceMatrix Jpt(J,3,3);
557  for (int k=0; k<3; k++)
558  {
559  for (int l=0; l<3; l++)
560  {
561  x2(k,l) = 0.0;
562  for (int u=0; u<3; u++)
563  {
564  for (int v=0; v<3; v++)
565  {
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);
570  }
571  }
572  }
573  }
574  // x3_ijkl = -2 B_ik δ_jl
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]; // B(0,1)
579  B[4] = J[0]*J[2] + J[3]*J[5] + J[6]*J[8]; // B(0,2)
580  B[5] = J[1]*J[2] + J[4]*J[5] + J[7]*J[8]; // B(1,2)
581  const double b_p[9] =
582  {
583  B[0], B[3], B[4],
584  B[3], B[1], B[5],
585  B[4], B[5], B[2]
586  };
587  ConstDeviceMatrix b(b_p,3,3);
588  for (int k=0; k<3; k++)
589  {
590  for (int l=0; l<3; l++)
591  {
592  const double jl = j==l ? 1.0 : 0.0;
593  x3(k,l) = -2.0 * b(i,k) * jl;
594  }
595  }
596  // ddI2 = x1 + x2 + x3
597  DeviceMatrix ddi2(ddI2,3,3);
598  for (int k=0; k<3; k++)
599  {
600  for (int l=0; l<3; l++)
601  {
602  ddi2(k,l) = x1(k,l) + x2(k,l) + x3(k,l);
603  }
604  }
605  return ddI2;
606  }
607 
608  // ddI2b = X1 + X2 + X3
609  // X1_ijkl = 16/9 det(J)^{-10/3} I2 dI3b_ij dI3b_kl +
610  // 4/3 det(J)^{-10/3} I2 dI3b_il dI3b_kj
611  // X2_ijkl = -4/3 det(J)^{-7/3} (dI2_ij dI3b_kl + dI2_kl dI3b_ij)
612  // X3_ijkl = det(J)^{-4/3} ddI2_ijkl
613  MFEM_HOST_DEVICE inline double *Get_ddI2b(int i, int j)
614  {
615  double X1_p[9], X2_p[9], X3_p[9];
616  // X1_ijkl = 16/9 det(J)^{-10/3} I2 dI3b_ij dI3b_kl +
617  // 4/3 det(J)^{-10/3} I2 dI3b_il dI3b_kj
618  double sign_detJ;
619  DeviceMatrix X1(X1_p,3,3);
620  const double I3b_p = Get_I3b_p(); // I3b^{-2/3}
621  const double I3b = Get_I3b(sign_detJ); // det(J)
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);
626  ConstDeviceMatrix di3b(Get_dI3b(sign_detJ),3,3);
627  for (int k=0; k<3; k++)
628  {
629  for (int l=0; l<3; l++)
630  {
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);
633  X1(k,l) = up + down;
634  }
635  }
636  // X2_ijkl = -4/3 det(J)^{-7/3} (dI2_ij dI3b_kl + dI2_kl dI3b_ij)
637  DeviceMatrix X2(X2_p,3,3);
638  ConstDeviceMatrix di2(Get_dI2(),3,3);
639  for (int k=0; k<3; k++)
640  {
641  for (int l=0; l<3; l++)
642  {
643  X2(k,l) = -(4./3.)*I3b_p73*(di2(i,j)*di3b(k,l)+di2(k,l)*di3b(i,j));
644  }
645  }
646  // X3_ijkl = det(J)^{-4/3} ddI2_ijkl
647  DeviceMatrix X3(X3_p,3,3);
648  ConstDeviceMatrix ddi2(Get_ddI2(i,j),3,3);
649  for (int k=0; k<3; k++)
650  {
651  for (int l=0; l<3; l++)
652  {
653  X3(k,l) = I3b_p43 * ddi2(k,l);
654  }
655  }
656  // ddI2b = X1 + X2 + X3
657  DeviceMatrix ddi2b(ddI2b,3,3);
658  for (int k=0; k<3; k++)
659  {
660  for (int l=0; l<3; l++)
661  {
662  ddi2b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
663  }
664  }
665  return ddI2b;
666  }
667 
668  // dI3b = adj(J)^T
669  // ddI3b_ijkl = (1/I3b) (δ_ks δ_it - δ_kt δ_si) dI3b_tj dI3b_sl
670  MFEM_HOST_DEVICE inline double *Get_ddI3b(int i, int j)
671  {
672  const double c1 = 1./Get_I3b();
673  ConstDeviceMatrix di3b(dI3b,3,3);
674  DeviceMatrix ddi3b(ddI3b,3,3);
675  for (int k=0; k<3; k++)
676  {
677  for (int l=0; l<3; l++)
678  {
679  ddi3b(k,l) = 0.0;
680  for (int s=0; s<3; s++)
681  {
682  for (int t=0; t<3; t++)
683  {
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);
687  }
688  }
689  }
690  }
691  return ddI3b;
692  }
693 };
694 
695 } // namespace kernels
696 
697 } // namespace mfem
698 
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)
Definition: dinvariants.hpp:51
MFEM_HOST_DEVICE Buffers & ddI2b(double *b)
Definition: dinvariants.hpp:52
MFEM_HOST_DEVICE double Get_I2()
MFEM_HOST_DEVICE Buffers & B(double *b)
MFEM_HOST_DEVICE Buffers & dI2(double *b)
Definition: dinvariants.hpp:49
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)
Definition: dinvariants.hpp:48
MFEM_HOST_DEVICE Buffers & ddI3b(double *b)
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)
double b
Definition: lissajous.cpp:42
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.
Definition: dtensor.hpp:81
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
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()
Definition: dinvariants.hpp:91
MFEM_HOST_DEVICE void Get_B_offd()
MFEM_HOST_DEVICE double Get_I2b(double &sign_detJ)
Definition: dinvariants.hpp:67
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)
Definition: dinvariants.hpp:44
MFEM_HOST_DEVICE Buffers & dI1(double *b)
Definition: dinvariants.hpp:45
MFEM_HOST_DEVICE double Get_I2b()
MFEM_HOST_DEVICE InvariantsEvaluator3D(Buffers &b)
MFEM_HOST_DEVICE InvariantsEvaluator2D(Buffers &b)
Definition: dinvariants.hpp:62
MFEM_HOST_DEVICE Buffers & dI2b(double *b)
Definition: dinvariants.hpp:50
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()
Definition: dinvariants.hpp:80
RefCoord t[3]
const double alpha
Definition: ex15.cpp:369
MFEM_HOST_DEVICE double Get_I3b_p(double &sign_detJ)
RefCoord s[3]
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
MFEM_HOST_DEVICE Buffers & ddI1(double *b)
Definition: dinvariants.hpp:47
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()
Definition: dinvariants.hpp:86
MFEM_HOST_DEVICE Buffers & dI1b(double *b)
Definition: dinvariants.hpp:46
MFEM_HOST_DEVICE double * Get_dI1()
Definition: dinvariants.hpp:96
MFEM_HOST_DEVICE Buffers & ddI1(double *b)
MFEM_HOST_DEVICE double Get_I2b()
Definition: dinvariants.hpp:74