MFEM  v4.6.0
Finite element discretization library
dinvariants.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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  // Requires dI2b.
104  MFEM_HOST_DEVICE inline double *Get_dI1b()
105  {
106  // I1b = I1/I2b
107  // dI1b = (1/I2b)*dI1 - (I1/I2b^2)*dI2b = (2/I2b)*[J - (I1b/2)*dI2b]
108  const double c1 = 2.0/Get_I2b();
109  const double c2 = Get_I1b()/2.0;
110  Get_dI2b();
111  dI1b[0] = c1*(J[0] - c2*dI2b[0]);
112  dI1b[1] = c1*(J[1] - c2*dI2b[1]);
113  dI1b[2] = c1*(J[2] - c2*dI2b[2]);
114  dI1b[3] = c1*(J[3] - c2*dI2b[3]);
115  return dI1b;
116  }
117 
118  // Requires dI2b.
119  MFEM_HOST_DEVICE inline double *Get_dI2()
120  {
121  // I2 = I2b^2
122  // dI2 = 2*I2b*dI2b = 2*det(J)*adj(J)^T
123  const double c1 = 2*Get_I2b();
124  Get_dI2b();
125  dI2[0] = c1*dI2b[0];
126  dI2[1] = c1*dI2b[1];
127  dI2[2] = c1*dI2b[2];
128  dI2[3] = c1*dI2b[3];
129  return dI2;
130  }
131 
132  MFEM_HOST_DEVICE inline double *Get_dI2b()
133  {
134  // I2b = det(J)
135  // dI2b = adj(J)^T
136  double sign_detJ;
137  Get_I2b(sign_detJ);
138  dI2b[0] = sign_detJ*J[3];
139  dI2b[1] = -sign_detJ*J[2];
140  dI2b[2] = -sign_detJ*J[1];
141  dI2b[3] = sign_detJ*J[0];
142  return dI2b;
143  }
144 
145  // ddI1_ijkl = 2 I_ijkl = 2 δ_ik δ_jl
146  MFEM_HOST_DEVICE inline double *Get_ddI1(int i, int j)
147  {
148  // ddI1_ijkl = 2 I_ijkl = 2 δ_ik δ_jl
149  DeviceMatrix ddi1(ddI1,2,2);
150  for (int k=0; k<2; k++)
151  {
152  for (int l=0; l<2; l++)
153  {
154  ddi1(k,l) = (i==k && j==l) ? 2.0 : 0.0;
155  }
156  }
157  return ddI1;
158  }
159 
160  // Requires dI2b + ddI1.
161  // ddI1b = X1 + X2 + X3, where
162  // X1_ijkl = (I1b/I2) [ dI2b_ij dI2b_kl + dI2b_kj dI2b_il ]
163  // X2_ijkl = (1/I2b) ddI1_ijkl
164  // X3_ijkl = -(2/I2) (J_ij dI2b_kl + dI2b_ij J_kl)
165  MFEM_HOST_DEVICE inline double *Get_ddI1b(int i, int j)
166  {
167  double X1_p[4], X2_p[4], X3_p[4];
168 
169  // X1_ijkl = (I1b/I2) [ dI2b_ij dI2b_kl + dI2b_kj dI2b_il ]
170  const double I2 = Get_I2();
171  const double I1b = Get_I1b();
172  ConstDeviceMatrix di2b(Get_dI2b(),2,2);
173  const double alpha = I1b / I2;
174  DeviceMatrix X1(X1_p,2,2);
175  for (int k=0; k<2; k++)
176  {
177  for (int l=0; l<2; l++)
178  {
179  X1(k,l) = alpha * (di2b(i,j)*di2b(k,l) + di2b(k,j)*di2b(i,l));
180  }
181  }
182  // X2_ijkl = (1/I2b) ddI1_ijkl
183  DeviceMatrix X2(X2_p,2,2);
184  const double beta = 1.0 / Get_I2b();
185  ConstDeviceMatrix ddi1(Get_ddI1(i,j),2,2);
186  for (int k=0; k<2; k++)
187  {
188  for (int l=0; l<2; l++)
189  {
190  X2(k,l) = beta * ddi1(k,l);
191  }
192  }
193  // X3_ijkl = -(2/I2) (J_ij dI2b_kl + dI2b_ij J_kl)
194  DeviceMatrix X3(X3_p,2,2);
195  const double gamma = -2.0/Get_I2();
196  ConstDeviceMatrix Jpt(J,2,2);
197  for (int k=0; k<2; k++)
198  {
199  for (int l=0; l<2; l++)
200  {
201  X3(k,l) = gamma * (Jpt(i,j)*di2b(k,l) + di2b(i,j)*Jpt(k,l));
202  }
203  }
204  DeviceMatrix ddi1b(ddI1b,2,2);
205  for (int k=0; k<2; k++)
206  {
207  for (int l=0; l<2; l++)
208  {
209  ddi1b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
210  }
211  }
212  return ddI1b;
213  }
214 
215  // Requires dI2b.
216  // ddI2_ijkl = 2 dI2b_ij dI2b_kl + 2 (dI2b_ij dI2b_kl - dI2b_kj dI2b_il)
217  MFEM_HOST_DEVICE inline double *Get_ddI2(int i, int j)
218  {
219  DeviceMatrix ddi2(ddI2,2,2);
220  ConstDeviceMatrix di2b(Get_dI2b(),2,2);
221  for (int k=0; k<2; k++)
222  {
223  for (int l=0; l<2; l++)
224  {
225  ddi2(k,l) = 2*di2b(i,j)*di2b(k,l)
226  + 2*(di2b(i,j)*di2b(k,l) - di2b(k,j)*di2b(i,l));
227  }
228  }
229  return ddI2;
230  }
231 
232  // Requires dI2b.
233  // ddI2b_ijkl = (1/I2b) (δ_ks δ_it - δ_kt δ_si) dI2b_tj dI2b_sl
234  MFEM_HOST_DEVICE inline double *Get_ddI2b(int i, int j)
235  {
236  DeviceMatrix ddi2b(ddI2b,2,2);
237  const double alpha = 1.0/Get_I2b();
238  ConstDeviceMatrix di2b(Get_dI2b(),2,2);
239  for (int k=0; k<2; k++)
240  {
241  for (int l=0; l<2; l++)
242  {
243  ddi2b(k,l) = 0.0;
244  for (int s=0; s<2; s++)
245  {
246  for (int t=0; t<2; t++)
247  {
248  const double ks_it = k==s && i==t ? 1.0 : 0.0;
249  const double kt_si = k==t && s==i ? 1.0 : 0.0;
250  ddi2b(k,l) += alpha * (ks_it - kt_si) * di2b(t,j) * di2b(s,l);
251  }
252  }
253  }
254  }
255  return ddI2b;
256  }
257 };
258 
259 
261 {
262 public:
263  class Buffers
264  {
265  friend class InvariantsEvaluator3D;
266  private:
267  const double * J_ = nullptr;
268  double * B_ = nullptr;
269  double * dI1_ = nullptr;
270  double * dI1b_ = nullptr;
271  double * ddI1_ = nullptr;
272  double * ddI1b_ = nullptr;
273  double * dI2_ = nullptr;
274  double * dI2b_ = nullptr;
275  double * ddI2_ = nullptr;
276  double * ddI2b_ = nullptr;
277  double * dI3b_ = nullptr;
278  double * ddI3b_ = nullptr;
279  public:
280  MFEM_HOST_DEVICE Buffers() {}
281  MFEM_HOST_DEVICE Buffers &J(const double *b) { J_ = b; return *this; }
282  MFEM_HOST_DEVICE Buffers &B(double *b) { B_ = b; return *this; }
283  MFEM_HOST_DEVICE Buffers &dI1(double *b) { dI1_ = b; return *this; }
284  MFEM_HOST_DEVICE Buffers &dI1b(double *b) { dI1b_ = b; return *this; }
285  MFEM_HOST_DEVICE Buffers &ddI1(double *b) { ddI1_ = b; return *this; }
286  MFEM_HOST_DEVICE Buffers &ddI1b(double *b) { ddI1b_ = b; return *this; }
287  MFEM_HOST_DEVICE Buffers &dI2(double *b) { dI2_ = b; return *this; }
288  MFEM_HOST_DEVICE Buffers &dI2b(double *b) { dI2b_ = b; return *this; }
289  MFEM_HOST_DEVICE Buffers &ddI2(double *b) { ddI2_ = b; return *this; }
290  MFEM_HOST_DEVICE Buffers &ddI2b(double *b) { ddI2b_ = b; return *this; }
291  MFEM_HOST_DEVICE Buffers &dI3b(double *b) { dI3b_ = b; return *this; }
292  MFEM_HOST_DEVICE Buffers &ddI3b(double *b) { ddI3b_ = b; return *this; }
293  };
294 
295 private:
296  double const * const J;
297  double * const B;
298  double * const dI1, * const dI1b, * const ddI1, * const ddI1b;
299  double * const dI2, * const dI2b, * const ddI2, * const ddI2b;
300  double * const dI3b, * const ddI3b;
301 
302 public:
303  MFEM_HOST_DEVICE
305  J(b.J_), B(b.B_),
306  dI1(b.dI1_), dI1b(b.dI1b_), ddI1(b.ddI1_), ddI1b(b.ddI1b_),
307  dI2(b.dI2_), dI2b(b.dI2b_), ddI2(b.ddI2_), ddI2b(b.ddI2b_),
308  dI3b(b.dI3b_), ddI3b(b.ddI3b_) { }
309 
310  MFEM_HOST_DEVICE inline double Get_I3b(double &sign_detJ) // det(J) + sign
311  {
312  const double I3b = + J[0]*(J[4]*J[8] - J[7]*J[5])
313  - J[1]*(J[3]*J[8] - J[5]*J[6])
314  + J[2]*(J[3]*J[7] - J[4]*J[6]);
315  sign_detJ = I3b >= 0.0 ? 1.0 : -1.0;
316  return sign_detJ * I3b;
317  }
318 
319  MFEM_HOST_DEVICE inline double Get_I3b() // det(J)
320  {
321  const double I3b = + J[0]*(J[4]*J[8] - J[7]*J[5])
322  - J[1]*(J[3]*J[8] - J[5]*J[6])
323  + J[2]*(J[3]*J[7] - J[4]*J[6]);
324  return I3b;
325  }
326 
327  MFEM_HOST_DEVICE inline double Get_I3() // det(J)^{2}
328  {
329  const double I3b = Get_I3b();
330  return I3b * I3b;
331  }
332 
333  MFEM_HOST_DEVICE inline double Get_I3b_p() // I3b^{-2/3}
334  {
335  double sign_detJ;
336  const double i3b = Get_I3b(sign_detJ);
337  return sign_detJ * std::pow(i3b, -2./3.);
338  }
339 
340  MFEM_HOST_DEVICE inline double Get_I3b_p(double &sign_detJ) // I3b^{-2/3}
341  {
342  const double i3b = Get_I3b(sign_detJ);
343  return sign_detJ * std::pow(i3b, -2./3.);
344  }
345 
346  MFEM_HOST_DEVICE inline double Get_I1()
347  {
348  B[0] = J[0]*J[0] + J[3]*J[3] + J[6]*J[6];
349  B[1] = J[1]*J[1] + J[4]*J[4] + J[7]*J[7];
350  B[2] = J[2]*J[2] + J[5]*J[5] + J[8]*J[8];
351  const double I1 = B[0] + B[1] + B[2];
352  return I1;
353  }
354 
355  MFEM_HOST_DEVICE inline
356  double Get_I1b() // det(J)^{-2/3}*I_1 = I_1/I_3^{1/3}
357  {
358  const double I1b = Get_I1() * Get_I3b_p();
359  return I1b;
360  }
361 
362  MFEM_HOST_DEVICE inline void Get_B_offd()
363  {
364  // B = J J^t
365  // B[3]=B(0,1), B[4]=B(0,2), B[5]=B(1,2)
366  B[3] = J[0]*J[1] + J[3]*J[4] + J[6]*J[7]; // B(0,1)
367  B[4] = J[0]*J[2] + J[3]*J[5] + J[6]*J[8]; // B(0,2)
368  B[5] = J[1]*J[2] + J[4]*J[5] + J[7]*J[8]; // B(1,2)
369  }
370 
371  MFEM_HOST_DEVICE inline double Get_I2()
372  {
373  Get_B_offd();
374  const double I1 = Get_I1();
375  const double BF2 = B[0]*B[0] + B[1]*B[1] + B[2]*B[2] +
376  2*(B[3]*B[3] + B[4]*B[4] + B[5]*B[5]);
377  const double I2 = (I1*I1 - BF2)/2;
378  return I2;
379  }
380 
381  MFEM_HOST_DEVICE inline double Get_I2b() // I2b = I2*I3b^{-4/3}
382  {
383  const double I3b_p = Get_I3b_p();
384  return Get_I2() * I3b_p * I3b_p;
385  }
386 
387  MFEM_HOST_DEVICE inline double *Get_dI1()
388  {
389  for (int i = 0; i < 9; i++) { dI1[i] = 2*J[i]; }
390  return dI1;
391  }
392 
393  MFEM_HOST_DEVICE inline double *Get_dI1b()
394  {
395  // I1b = I3b^{-2/3}*I1
396  // dI1b = 2*I3b^{-2/3}*(J - (1/3)*I1/I3b*dI3b)
397  double sign_detJ;
398  const double I3b = Get_I3b(sign_detJ);
399  const double I3b_p = Get_I3b_p();
400  const double c1 = 2.0 * I3b_p;
401  const double c2 = Get_I1()/(3.0 * I3b);
402  Get_dI3b(sign_detJ);
403  for (int i = 0; i < 9; i++) { dI1b[i] = c1*(J[i] - c2*dI3b[i]); }
404  return dI1b;
405  }
406 
407  MFEM_HOST_DEVICE inline double *Get_dI2()
408  {
409  // dI2 = 2 I_1 J - 2 J J^t J = 2 (I_1 I - B) J
410  const double I1 = Get_I1();
411  Get_B_offd();
412  // B[0]=B(0,0), B[1]=B(1,1), B[2]=B(2,2)
413  // B[3]=B(0,1), B[4]=B(0,2), B[5]=B(1,2)
414  const double C[6] =
415  {
416  2*(I1 - B[0]), 2*(I1 - B[1]), 2*(I1 - B[2]),
417  -2*B[3], -2*B[4], -2*B[5]
418  };
419  // | C[0] C[3] C[4] | | J[0] J[3] J[6] |
420  // dI2 = | C[3] C[1] C[5] | | J[1] J[4] J[7] |
421  // | C[4] C[5] C[2] | | J[2] J[5] J[8] |
422  dI2[0] = C[0]*J[0] + C[3]*J[1] + C[4]*J[2];
423  dI2[1] = C[3]*J[0] + C[1]*J[1] + C[5]*J[2];
424  dI2[2] = C[4]*J[0] + C[5]*J[1] + C[2]*J[2];
425 
426  dI2[3] = C[0]*J[3] + C[3]*J[4] + C[4]*J[5];
427  dI2[4] = C[3]*J[3] + C[1]*J[4] + C[5]*J[5];
428  dI2[5] = C[4]*J[3] + C[5]*J[4] + C[2]*J[5];
429 
430  dI2[6] = C[0]*J[6] + C[3]*J[7] + C[4]*J[8];
431  dI2[7] = C[3]*J[6] + C[1]*J[7] + C[5]*J[8];
432  dI2[8] = C[4]*J[6] + C[5]*J[7] + C[2]*J[8];
433  return dI2;
434  }
435 
436  MFEM_HOST_DEVICE inline double *Get_dI2b()
437  {
438  // I2b = det(J)^{-4/3}*I2 = I3b^{-4/3}*I2
439  // dI2b = (-4/3)*I3b^{-7/3}*I2*dI3b + I3b^{-4/3}*dI2
440  // = I3b^{-4/3} * [ dI2 - (4/3)*I2/I3b*dI3b ]
441  double sign_detJ;
442  const double I2 = Get_I2();
443  const double I3b_p = Get_I3b_p();
444  const double I3b = Get_I3b(sign_detJ);
445  const double c1 = I3b_p*I3b_p;
446  const double c2 = (4*I2/I3b)/3;
447  Get_dI2();
448  Get_dI3b(sign_detJ);
449  for (int i = 0; i < 9; i++) { dI2b[i] = c1*(dI2[i] - c2*dI3b[i]); }
450  return dI2b;
451  }
452 
453  MFEM_HOST_DEVICE inline double *Get_dI3b(const double sign_detJ)
454  {
455  // I3b = det(J)
456  // dI3b = adj(J)^T
457  dI3b[0] = sign_detJ*(J[4]*J[8] - J[5]*J[7]); // 0 3 6
458  dI3b[1] = sign_detJ*(J[5]*J[6] - J[3]*J[8]); // 1 4 7
459  dI3b[2] = sign_detJ*(J[3]*J[7] - J[4]*J[6]); // 2 5 8
460  dI3b[3] = sign_detJ*(J[2]*J[7] - J[1]*J[8]);
461  dI3b[4] = sign_detJ*(J[0]*J[8] - J[2]*J[6]);
462  dI3b[5] = sign_detJ*(J[1]*J[6] - J[0]*J[7]);
463  dI3b[6] = sign_detJ*(J[1]*J[5] - J[2]*J[4]);
464  dI3b[7] = sign_detJ*(J[2]*J[3] - J[0]*J[5]);
465  dI3b[8] = sign_detJ*(J[0]*J[4] - J[1]*J[3]);
466  return dI3b;
467  }
468 
469  // ddI1_ijkl = 2 I_ijkl = 2 δ_ik δ_jl
470  MFEM_HOST_DEVICE inline double *Get_ddI1(int i, int j)
471  {
472  DeviceMatrix ddi1(ddI1,3,3);
473  for (int k=0; k<3; k++)
474  {
475  for (int l=0; l<3; l++)
476  {
477  const double I_ijkl = (i==k && j==l) ? 1.0 : 0.0;
478  ddi1(k,l) = 2.0 * I_ijkl;
479  }
480  }
481  return ddI1;
482  }
483 
484  // ddI1b = X1 + X2 + X3, where
485  // X1_ijkl = (2/3*I1b/I3) [ 2/3 dI3b_ij dI3b_kl + dI3b_kj dI3b_il ]
486  // X2_ijkl = (I3b^{-2/3}) ddI1_ijkl
487  // X3_ijkl = -(4/3*I3b^{-5/3}) (J_ij dI3b_kl + dI3b_ij J_kl)
488  MFEM_HOST_DEVICE inline double *Get_ddI1b(int i, int j)
489  {
490  // X1_ijkl = (2/3*I1b/I3) [ 2/3 dI3b_ij dI3b_kl + dI3b_kj dI3b_il ]
491  double sign_detJ;
492  Get_I3b(sign_detJ);
493  double X1_p[9], X2_p[9], X3_p[9];
494  DeviceMatrix X1(X1_p,3,3);
495  const double I3 = Get_I3();
496  const double I1b = Get_I1b();
497  const double alpha = (2./3.)*I1b/I3;
498  ConstDeviceMatrix di3b(Get_dI3b(sign_detJ),3,3);
499  for (int k=0; k<3; k++)
500  {
501  for (int l=0; l<3; l++)
502  {
503  X1(k,l) = alpha * ((2./3.)*di3b(i,j) * di3b(k,l) +
504  di3b(k,j)*di3b(i,l));
505  }
506  }
507  // ddI1_ijkl = 2 δ_ik δ_jl
508  // X2_ijkl = (I3b^{-2/3}) ddI1_ijkl
509  DeviceMatrix X2(X2_p,3,3);
510  const double beta = Get_I3b_p();
511  for (int k=0; k<3; k++)
512  {
513  for (int l=0; l<3; l++)
514  {
515  const double ddI1_ijkl = (i==k && j==l) ? 2.0 : 0.0;
516  X2(k,l) = beta * ddI1_ijkl;
517  }
518  }
519  // X3_ijkl = -(4/3*I3b^{-5/3}) (J_ij dI3b_kl + dI3b_ij J_kl)
520  DeviceMatrix X3(X3_p,3,3);
521  const double I3b = Get_I3b();
522  const double gamma = -(4./3.)*Get_I3b_p()/I3b;
523  ConstDeviceMatrix Jpt(J,3,3);
524  for (int k=0; k<3; k++)
525  {
526  for (int l=0; l<3; l++)
527  {
528  X3(k,l) = gamma * (Jpt(i,j) * di3b(k,l) + di3b(i,j) * Jpt(k,l));
529  }
530  }
531  DeviceMatrix ddi1b(ddI1b,3,3);
532  for (int k=0; k<3; k++)
533  {
534  for (int l=0; l<3; l++)
535  {
536  ddi1b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
537  }
538  }
539  return ddI1b;
540  }
541 
542  // ddI2 = x1 + x2 + x3
543  // x1_ijkl = (2 I1) δ_ik δ_jl
544  // x2_ijkl = 2 ( 2 δ_ku δ_iv - δ_ik δ_uv - δ_kv δ_iu ) J_vj J_ul
545  // x3_ijkl = -2 (J J^t)_ik δ_jl = -2 B_ik δ_jl
546  MFEM_HOST_DEVICE inline double *Get_ddI2(int i, int j)
547  {
548  double x1_p[9], x2_p[9], x3_p[9];
549  DeviceMatrix x1(x1_p,3,3), x2(x2_p,3,3), x3(x3_p,3,3);
550  // x1_ijkl = (2 I1) δ_ik δ_jl
551  const double I1 = Get_I1();
552  for (int k=0; k<3; k++)
553  {
554  for (int l=0; l<3; l++)
555  {
556  const double ik_jl = (i==k && j==l) ? 1.0 : 0.0;
557  x1(k,l) = 2.0 * I1 * ik_jl;
558  }
559  }
560  // x2_ijkl = 2 ( 2 δ_ku δ_iv - δ_ik δ_uv - δ_kv δ_iu ) J_vj J_ul
561  ConstDeviceMatrix Jpt(J,3,3);
562  for (int k=0; k<3; k++)
563  {
564  for (int l=0; l<3; l++)
565  {
566  x2(k,l) = 0.0;
567  for (int u=0; u<3; u++)
568  {
569  for (int v=0; v<3; v++)
570  {
571  const double ku_iv = k==u && i==v ? 1.0 : 0.0;
572  const double ik_uv = i==k && u==v ? 1.0 : 0.0;
573  const double kv_iu = k==v && i==u ? 1.0 : 0.0;
574  x2(k,l) += 2.0*(2.*ku_iv-ik_uv-kv_iu)*Jpt(v,j)*Jpt(u,l);
575  }
576  }
577  }
578  }
579  // x3_ijkl = -2 B_ik δ_jl
580  B[0] = J[0]*J[0] + J[3]*J[3] + J[6]*J[6];
581  B[1] = J[1]*J[1] + J[4]*J[4] + J[7]*J[7];
582  B[2] = J[2]*J[2] + J[5]*J[5] + J[8]*J[8];
583  B[3] = J[0]*J[1] + J[3]*J[4] + J[6]*J[7]; // B(0,1)
584  B[4] = J[0]*J[2] + J[3]*J[5] + J[6]*J[8]; // B(0,2)
585  B[5] = J[1]*J[2] + J[4]*J[5] + J[7]*J[8]; // B(1,2)
586  const double b_p[9] =
587  {
588  B[0], B[3], B[4],
589  B[3], B[1], B[5],
590  B[4], B[5], B[2]
591  };
592  ConstDeviceMatrix b(b_p,3,3);
593  for (int k=0; k<3; k++)
594  {
595  for (int l=0; l<3; l++)
596  {
597  const double jl = j==l ? 1.0 : 0.0;
598  x3(k,l) = -2.0 * b(i,k) * jl;
599  }
600  }
601  // ddI2 = x1 + x2 + x3
602  DeviceMatrix ddi2(ddI2,3,3);
603  for (int k=0; k<3; k++)
604  {
605  for (int l=0; l<3; l++)
606  {
607  ddi2(k,l) = x1(k,l) + x2(k,l) + x3(k,l);
608  }
609  }
610  return ddI2;
611  }
612 
613  // ddI2b = X1 + X2 + X3
614  // X1_ijkl = 16/9 det(J)^{-10/3} I2 dI3b_ij dI3b_kl +
615  // 4/3 det(J)^{-10/3} I2 dI3b_il dI3b_kj
616  // X2_ijkl = -4/3 det(J)^{-7/3} (dI2_ij dI3b_kl + dI2_kl dI3b_ij)
617  // X3_ijkl = det(J)^{-4/3} ddI2_ijkl
618  MFEM_HOST_DEVICE inline double *Get_ddI2b(int i, int j)
619  {
620  double X1_p[9], X2_p[9], X3_p[9];
621  // X1_ijkl = 16/9 det(J)^{-10/3} I2 dI3b_ij dI3b_kl +
622  // 4/3 det(J)^{-10/3} I2 dI3b_il dI3b_kj
623  double sign_detJ;
624  DeviceMatrix X1(X1_p,3,3);
625  const double I3b_p = Get_I3b_p(); // I3b^{-2/3}
626  const double I3b = Get_I3b(sign_detJ); // det(J)
627  const double I2 = Get_I2();
628  const double I3b_p43 = I3b_p*I3b_p;
629  const double I3b_p73 = I3b_p*I3b_p/I3b;
630  const double I3b_p103 = I3b_p*I3b_p/(I3b*I3b);
631  ConstDeviceMatrix di3b(Get_dI3b(sign_detJ),3,3);
632  for (int k=0; k<3; k++)
633  {
634  for (int l=0; l<3; l++)
635  {
636  const double up = (16./9.)*I3b_p103*I2*di3b(i,j)*di3b(k,l);
637  const double down = (4./3.)*I3b_p103*I2*di3b(i,l)*di3b(k,j);
638  X1(k,l) = up + down;
639  }
640  }
641  // X2_ijkl = -4/3 det(J)^{-7/3} (dI2_ij dI3b_kl + dI2_kl dI3b_ij)
642  DeviceMatrix X2(X2_p,3,3);
643  ConstDeviceMatrix di2(Get_dI2(),3,3);
644  for (int k=0; k<3; k++)
645  {
646  for (int l=0; l<3; l++)
647  {
648  X2(k,l) = -(4./3.)*I3b_p73*(di2(i,j)*di3b(k,l)+di2(k,l)*di3b(i,j));
649  }
650  }
651  // X3_ijkl = det(J)^{-4/3} ddI2_ijkl
652  DeviceMatrix X3(X3_p,3,3);
653  ConstDeviceMatrix ddi2(Get_ddI2(i,j),3,3);
654  for (int k=0; k<3; k++)
655  {
656  for (int l=0; l<3; l++)
657  {
658  X3(k,l) = I3b_p43 * ddi2(k,l);
659  }
660  }
661  // ddI2b = X1 + X2 + X3
662  DeviceMatrix ddi2b(ddI2b,3,3);
663  for (int k=0; k<3; k++)
664  {
665  for (int l=0; l<3; l++)
666  {
667  ddi2b(k,l) = X1(k,l) + X2(k,l) + X3(k,l);
668  }
669  }
670  return ddI2b;
671  }
672 
673  // dI3b = adj(J)^T
674  // ddI3b_ijkl = (1/I3b) (δ_ks δ_it - δ_kt δ_si) dI3b_tj dI3b_sl
675  MFEM_HOST_DEVICE inline double *Get_ddI3b(int i, int j)
676  {
677  const double c1 = 1./Get_I3b();
678  ConstDeviceMatrix di3b(dI3b,3,3);
679  DeviceMatrix ddi3b(ddI3b,3,3);
680  for (int k=0; k<3; k++)
681  {
682  for (int l=0; l<3; l++)
683  {
684  ddi3b(k,l) = 0.0;
685  for (int s=0; s<3; s++)
686  {
687  for (int t=0; t<3; t++)
688  {
689  const double ks_it = k==s && i==t ? 1.0 : 0.0;
690  const double kt_si = k==t && s==i ? 1.0 : 0.0;
691  ddi3b(k,l) += c1*(ks_it-kt_si)*di3b(t,j)*di3b(s,l);
692  }
693  }
694  }
695  }
696  return ddI3b;
697  }
698 };
699 
700 } // namespace kernels
701 
702 } // namespace mfem
703 
704 #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)
Vector beta
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
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:22
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