MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
invariants.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_INVARIANTS_HPP
13 #define MFEM_INVARIANTS_HPP
14 
15 #include "../config/config.hpp"
16 #include "../general/error.hpp"
17 #include <cmath>
18 
19 namespace mfem
20 {
21 
22 // Matrix invariants and their derivatives for 2x2 and 3x3 matrices.
23 
24 /** @brief Auxiliary class used as the default for the second template parameter
25  in the classes InvariantsEvaluator2D and InvariantsEvaluator3D. */
26 template <typename scalar_t>
27 struct ScalarOps
28 {
29  static scalar_t sign(const scalar_t &a)
30  { return (a >= scalar_t(0)) ? scalar_t(1) : scalar_t(-1); }
31 
32  static scalar_t pow(const scalar_t &x, int m, int n)
33  { return std::pow(x, scalar_t(m)/n); }
34 };
35 
36 
37 /** @brief Auxiliary class for evaluating the 2x2 matrix invariants and their
38  first and second derivatives. */
39 /**
40  The type `scalar_t` must support the standard operations:
41 
42  =, +=, -=, +, -, *, /, unary -, int*scalar_t, int/scalar_t, scalar_t/int
43 
44  The type `scalar_ops` must define the static method:
45 
46  scalar_t sign(const scalar_t &);
47 */
48 template <typename scalar_t, typename scalar_ops = ScalarOps<scalar_t> >
49 class InvariantsEvaluator2D
50 {
51 protected:
52  // Transformation Jacobian
53  const scalar_t *J;
54 
55  // Invariants: I_1 = ||J||_F^2, \bar{I}_1 = I_1/det(J), \bar{I}_2 = det(J).
56  scalar_t I1, I1b, I2b;
57 
58  // Derivatives of I1, I1b, I2, and I2b using column-major storage.
59  scalar_t dI1[4], dI1b[4], dI2[4], dI2b[4];
60 
61  int D_height, alloc_height;
62  const scalar_t *D; // Always points to external data or is empty
63  scalar_t *DaJ, *DJt, *DXt, *DYt;
64 
65  enum EvalMasks
66  {
67  HAVE_I1 = 1,
68  HAVE_I1b = 2,
69  HAVE_I2b = 4,
70  HAVE_dI1 = 8,
71  HAVE_dI1b = 16,
72  HAVE_dI2 = 32,
73  HAVE_dI2b = 64,
74  HAVE_DaJ = 128, // D adj(J) = D dI2b^t
75  HAVE_DJt = 256 // D J^t
76  };
77 
78  // Bitwise OR of EvalMasks
79  int eval_state;
80 
81  bool dont(int have_mask) const { return !(eval_state & have_mask); }
82 
83  void Eval_I1()
84  {
85  eval_state |= HAVE_I1;
86  I1 = J[0]*J[0] + J[1]*J[1] + J[2]*J[2] + J[3]*J[3];
87  }
88  void Eval_I1b()
89  {
90  eval_state |= HAVE_I1b;
91  I1b = Get_I1()/Get_I2b();
92  }
93  void Eval_I2b()
94  {
95  eval_state |= HAVE_I2b;
96  const scalar_t det = J[0]*J[3] - J[1]*J[2];
97  I2b = det;
98  }
99  void Eval_dI1()
100  {
101  eval_state |= HAVE_dI1;
102  dI1[0] = 2*J[0]; dI1[2] = 2*J[2];
103  dI1[1] = 2*J[1]; dI1[3] = 2*J[3];
104  }
105  void Eval_dI1b()
106  {
107  eval_state |= HAVE_dI1b;
108  // I1b = I1/I2b
109  // dI1b = (1/I2b)*dI1 - (I1/I2b^2)*dI2b = (2/I2b)*[J - (I1b/2)*dI2b]
110  const scalar_t c1 = 2/Get_I2b();
111  const scalar_t c2 = Get_I1b()/2;
112  Get_dI2b();
113  dI1b[0] = c1*(J[0] - c2*dI2b[0]);
114  dI1b[1] = c1*(J[1] - c2*dI2b[1]);
115  dI1b[2] = c1*(J[2] - c2*dI2b[2]);
116  dI1b[3] = c1*(J[3] - c2*dI2b[3]);
117  }
118  void Eval_dI2()
119  {
120  eval_state |= HAVE_dI2;
121  // I2 = I2b^2
122  // dI2 = 2*I2b*dI2b = 2*det(J)*adj(J)^T
123  const scalar_t 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  }
130  void Eval_dI2b()
131  {
132  eval_state |= HAVE_dI2b;
133  // I2b = det(J)
134  // dI2b = adj(J)^T
135  Get_I2b();
136  dI2b[0] = J[3];
137  dI2b[1] = -J[2];
138  dI2b[2] = -J[1];
139  dI2b[3] = J[0];
140  }
141  void Eval_DaJ() // D adj(J) = D dI2b^t
142  {
143  eval_state |= HAVE_DaJ;
144  Get_dI2b();
145  Eval_DZt(dI2b, &DaJ);
146  }
147  void Eval_DJt() // D J^t
148  {
149  eval_state |= HAVE_DJt;
150  Eval_DZt(J, &DJt);
151  }
152  void Eval_DZt(const scalar_t *Z, scalar_t **DZt_ptr)
153  {
154  MFEM_ASSERT(D != NULL, "");
155  const int nd = D_height;
156  scalar_t *DZt = *DZt_ptr;
157  if (DZt == NULL) { *DZt_ptr = DZt = new scalar_t[2*alloc_height]; }
158  for (int i = 0; i < nd; i++)
159  {
160  const int i0 = i+nd*0, i1 = i+nd*1;
161  DZt[i0] = D[i0]*Z[0] + D[i1]*Z[2];
162  DZt[i1] = D[i0]*Z[1] + D[i1]*Z[3];
163  }
164  }
165 
166 public:
167  /// The Jacobian should use column-major storage.
168  InvariantsEvaluator2D(const scalar_t *Jac = NULL)
169  : J(Jac), D_height(), alloc_height(), D(), DaJ(), DJt(), DXt(), DYt(),
170  eval_state(0) { }
171 
172  ~InvariantsEvaluator2D()
173  {
174  delete [] DYt;
175  delete [] DXt;
176  delete [] DJt;
177  delete [] DaJ;
178  }
179 
180  /// The Jacobian should use column-major storage.
181  void SetJacobian(const scalar_t *Jac) { J = Jac; eval_state = 0; }
182 
183  /// The @a Deriv matrix is `dof x 2`, using column-major storage.
184  void SetDerivativeMatrix(int height, const scalar_t *Deriv)
185  {
186  eval_state &= ~(HAVE_DaJ | HAVE_DJt);
187  if (alloc_height < height)
188  {
189  delete [] DYt; DYt = NULL;
190  delete [] DXt; DXt = NULL;
191  delete [] DJt; DJt = NULL;
192  delete [] DaJ; DaJ = NULL;
193  alloc_height = height;
194  }
195  D_height = height;
196  D = Deriv;
197  }
198 
199  scalar_t Get_I1() { if (dont(HAVE_I1 )) { Eval_I1(); } return I1; }
200  scalar_t Get_I1b() { if (dont(HAVE_I1b)) { Eval_I1b(); } return I1b; }
201  scalar_t Get_I2() { if (dont(HAVE_I2b)) { Eval_I2b(); } return I2b*I2b; }
202  scalar_t Get_I2b() { if (dont(HAVE_I2b)) { Eval_I2b(); } return I2b; }
203 
204  const scalar_t *Get_dI1()
205  {
206  if (dont(HAVE_dI1 )) { Eval_dI1(); } return dI1;
207  }
208  const scalar_t *Get_dI1b()
209  {
210  if (dont(HAVE_dI1b)) { Eval_dI1b(); } return dI1b;
211  }
212  const scalar_t *Get_dI2()
213  {
214  if (dont(HAVE_dI2)) { Eval_dI2(); } return dI2;
215  }
216  const scalar_t *Get_dI2b()
217  {
218  if (dont(HAVE_dI2b)) { Eval_dI2b(); } return dI2b;
219  }
220 
221  // Assemble operation for tensor X with components X_jslt:
222  // A(i+nd*j,k+nd*l) += (\sum_st w D_is X_jslt D_kt)
223  // 0 <= i,k < nd, 0 <= j,l,s,t < 2
224  // where nd is the height of D, i.e. the number of DOFs in one component.
225 
226  void Assemble_ddI1(scalar_t w, scalar_t *A)
227  {
228  // ddI1_jslt = 2 I_jslt = 2 δ_jl δ_st
229  // A(i+nd*j,k+nd*l) += (\sum_st 2 w D_is δ_jl δ_st D_kt)
230  // or
231  // A(i+nd*j,k+nd*l) += (2 w) (\sum_s D_is D_ks) δ_jl
232  // A(i+nd*j,k+nd*l) += (2 w) (D D^t)_ik δ_jl
233 
234  const int nd = D_height;
235  const int ah = 2*nd;
236  const scalar_t a = 2*w;
237  for (int i = 0; i < nd; i++)
238  {
239  const int i0 = i+nd*0, i1 = i+nd*1;
240  const scalar_t aDi[2] = { a*D[i0], a*D[i1] };
241  // k == i
242  const scalar_t aDDt_ii = aDi[0]*D[i0] + aDi[1]*D[i1];
243  A[i0+ah*i0] += aDDt_ii;
244  A[i1+ah*i1] += aDDt_ii;
245  // 0 <= k < i
246  for (int k = 0; k < i; k++)
247  {
248  const int k0 = k+nd*0, k1 = k+nd*1;
249  const scalar_t aDDt_ik = aDi[0]*D[k0] + aDi[1]*D[k1];
250  A[i0+ah*k0] += aDDt_ik;
251  A[k0+ah*i0] += aDDt_ik;
252  A[i1+ah*k1] += aDDt_ik;
253  A[k1+ah*i1] += aDDt_ik;
254  }
255  }
256  }
257  void Assemble_ddI1b(scalar_t w, scalar_t *A)
258  {
259  // ddI1b = X1 + X2 + X3, where
260  // X1_ijkl = (I1b/I2) [ (δ_ks δ_it + δ_kt δ_si) dI2b_tj dI2b_sl ]
261  // = (I1b/I2) [ dI2b_ij dI2b_kl + dI2b_kj dI2b_il ]
262  // X2_ijkl = (2/I2b) δ_ik δ_jl = (1/I2b) ddI1_ijkl
263  // X3_ijkl = -(2/I2) (δ_ks δ_it) (J_tj dI2b_sl + dI2b_tj J_sl)
264  // = -(2/I2) (J_ij dI2b_kl + dI2b_ij J_kl)
265  //
266  // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI1b_jslt D_kt)
267  // or
268  // A(i+nd*j,k+nd*l) +=
269  // w (I1b/I2) [(D dI2b^t)_ij (D dI2b^t)_kl +
270  // (D dI2b^t)_il (D dI2b^t)_kj]
271  // + w (2/I2b) δ_jl (D D^t)_ik
272  // - w (2/I2) [(D J^t)_ij (D dI2b^t)_kl + (D dI2b^t)_ij (D J^t)_kl]
273 
274  if (dont(HAVE_DaJ)) { Eval_DaJ(); }
275  if (dont(HAVE_DJt)) { Eval_DJt(); }
276  const int nd = D_height;
277  const int ah = 2*nd;
278  const scalar_t a = w*Get_I1b()/Get_I2();
279  const scalar_t b = 2*w/Get_I2b();
280  const scalar_t c = -2*w/Get_I2();
281  for (int i = 0; i < nd; i++)
282  {
283  const int i0 = i+nd*0, i1 = i+nd*1;
284  const scalar_t aDaJ_i[2] = { a*DaJ[i0], a*DaJ[i1] };
285  const scalar_t bD_i[2] = { b*D[i0], b*D[i1] };
286  const scalar_t cDJt_i[2] = { c*DJt[i0], c*DJt[i1] };
287  const scalar_t cDaJ_i[2] = { c*DaJ[i0], c*DaJ[i1] };
288  // k == i
289  {
290  // Symmetries: A2_ii_00 = A2_ii_11
291  const scalar_t A2_ii = bD_i[0]*D[i0] + bD_i[1]*D[i1];
292 
293  A[i0+ah*i0] += 2*(aDaJ_i[0] + cDJt_i[0])*DaJ[i0] + A2_ii;
294 
295  // Symmetries: A_ii_01 = A_ii_10
296  const scalar_t A_ii_01 =
297  (2*aDaJ_i[0] + cDJt_i[0])*DaJ[i1] + cDaJ_i[0]*DJt[i1];
298  A[i0+ah*i1] += A_ii_01;
299  A[i1+ah*i0] += A_ii_01;
300 
301  A[i1+ah*i1] += 2*(aDaJ_i[1] + cDJt_i[1])*DaJ[i1] + A2_ii;
302  }
303  // 0 <= k < i
304  for (int k = 0; k < i; k++)
305  {
306  const int k0 = k+nd*0, k1 = k+nd*1;
307  // Symmetries: A1_ik_01 = A1_ik_10 = A1_ki_01 = A1_ki_10
308  const scalar_t A1_ik_01 = aDaJ_i[0]*DaJ[k1] + aDaJ_i[1]*DaJ[k0];
309 
310  // Symmetries: A2_ik_00 = A2_ik_11 = A2_ki_00 = A2_ki_11
311  const scalar_t A2_ik = bD_i[0]*D[k0] + bD_i[1]*D[k1];
312 
313  const scalar_t A_ik_00 =
314  (2*aDaJ_i[0] + cDJt_i[0])*DaJ[k0] + A2_ik + cDaJ_i[0]*DJt[k0];
315  A[i0+ah*k0] += A_ik_00;
316  A[k0+ah*i0] += A_ik_00;
317 
318  const scalar_t A_ik_01 =
319  A1_ik_01 + cDJt_i[0]*DaJ[k1] + cDaJ_i[0]*DJt[k1];
320  A[i0+ah*k1] += A_ik_01;
321  A[k1+ah*i0] += A_ik_01;
322 
323  const scalar_t A_ik_10 =
324  A1_ik_01 + cDJt_i[1]*DaJ[k0] + cDaJ_i[1]*DJt[k0];
325  A[i1+ah*k0] += A_ik_10;
326  A[k0+ah*i1] += A_ik_10;
327 
328  const scalar_t A_ik_11 =
329  (2*aDaJ_i[1] + cDJt_i[1])*DaJ[k1] + A2_ik + cDaJ_i[1]*DJt[k1];
330  A[i1+ah*k1] += A_ik_11;
331  A[k1+ah*i1] += A_ik_11;
332  }
333  }
334  }
335  void Assemble_ddI2(scalar_t w, scalar_t *A)
336  {
337  // ddI2_ijkl = 2 (2 δ_ks δ_it - δ_kt δ_si) dI2b_tj dI2b_sl
338  // = 4 dI2b_ij dI2b_kl - 2 dI2b_kj dI2b_il
339  // = 2 dI2b_ij dI2b_kl + 2 (dI2b_ij dI2b_kl - dI2b_kj dI2b_il)
340  //
341  // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI2_jslt D_kt)
342  // or
343  // A(i+nd*j,k+nd*l) +=
344  // (\sum_st w D_is (4 dI2b_js dI2b_lt - 2 dI2b_ls dI2b_jt) D_kt)
345  // A(i+nd*j,k+nd*l) +=
346  // 2 w [2 (D dI2b^t)_ij (D dI2b^t)_kl - (D dI2b^t)_il (D dI2b^t)_kj]
347  //
348  // Note: the expression
349  // (D dI2b^t)_ij (D dI2b^t)_kl - (D dI2b^t)_il (D dI2b^t)_kj
350  // is the determinant of the 2x2 matrix formed by rows {i,k} and columns
351  // {j,l} from the matrix (D dI2b^t).
352 
353  if (dont(HAVE_DaJ)) { Eval_DaJ(); }
354  const int nd = D_height;
355  const int ah = 2*nd;
356  const scalar_t a = 2*w;
357  for (int i = 0; i < ah; i++)
358  {
359  const scalar_t avi = a*DaJ[i];
360  A[i+ah*i] += avi*DaJ[i];
361  for (int j = 0; j < i; j++)
362  {
363  const scalar_t aVVt_ij = avi*DaJ[j];
364  A[i+ah*j] += aVVt_ij;
365  A[j+ah*i] += aVVt_ij;
366  }
367  }
368  const int j = 1, l = 0;
369  for (int i = 0; i < nd; i++)
370  {
371  const int ij = i+nd*j, il = i+nd*l;
372  const scalar_t aDaJ_ij = a*DaJ[ij], aDaJ_il = a*DaJ[il];
373  for (int k = 0; k < i; k++)
374  {
375  const int kj = k+nd*j, kl = k+nd*l;
376  const scalar_t A_ijkl = aDaJ_ij*DaJ[kl] - aDaJ_il*DaJ[kj];
377  A[ij+ah*kl] += A_ijkl;
378  A[kl+ah*ij] += A_ijkl;
379  A[kj+ah*il] -= A_ijkl;
380  A[il+ah*kj] -= A_ijkl;
381  }
382  }
383  }
384  void Assemble_ddI2b(scalar_t w, scalar_t *A)
385  {
386  // ddI2b_ijkl = (1/I2b) (δ_ks δ_it - δ_kt δ_si) dI2b_tj dI2b_sl
387  // [j -> u], [l -> v], [i -> j], [k -> l]
388  // ddI2b_julv = (1/I2b) (δ_ls δ_jt - δ_lt δ_sj) dI2b_tu dI2b_sv
389  //
390  // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI2b_jslt D_kt)
391  // or
392  // A(i+nd*j,k+nd*l) += (\sum_uv w D_iu ddI2b_julv D_kv)
393  // A(i+nd*j,k+nd*l) +=
394  // (\sum_uvst (w/I2b)
395  // D_iu (δ_ls δ_jt - δ_lt δ_sj) dI2b_tu dI2b_sv D_kv)
396  // A(i+nd*j,k+nd*l) +=
397  // (\sum_st (w/I2b)
398  // (D dI2b^t)_it (δ_ls δ_jt - δ_lt δ_sj) (D dI2b^t)_ks)
399  // A(i+nd*j,k+nd*l) += (w/I2b)
400  // [ (D dI2b^t)_ij (D dI2b^t)_kl - (D dI2b^t)_il (D dI2b^t)_kj ]
401 
402  if (dont(HAVE_DaJ)) { Eval_DaJ(); }
403  const int nd = D_height;
404  const int ah = 2*nd;
405  const int j = 1, l = 0;
406  const scalar_t a = w/Get_I2b();
407  for (int i = 0; i < nd; i++)
408  {
409  const int ij = i+nd*j, il = i+nd*l;
410  const scalar_t aDaJ_ij = a*DaJ[ij], aDaJ_il = a*DaJ[il];
411  for (int k = 0; k < i; k++)
412  {
413  const int kj = k+nd*j, kl = k+nd*l;
414  const scalar_t A_ijkl = aDaJ_ij*DaJ[kl] - aDaJ_il*DaJ[kj];
415  A[ij+ah*kl] += A_ijkl;
416  A[kl+ah*ij] += A_ijkl;
417  A[kj+ah*il] -= A_ijkl;
418  A[il+ah*kj] -= A_ijkl;
419  }
420  }
421  }
422  // Assemble the contribution from the term: T_ijkl = X_ij Y_kl + Y_ij X_kl,
423  // where X and Y are pointers to 2x2 matrices stored in column-major layout.
424  //
425  // The contribution to the matrix A is given by:
426  // A(i+nd*j,k+nd*l) += \sum_st w D_is T_jslt D_kt
427  // or
428  // A(i+nd*j,k+nd*l) += \sum_st w D_is (X_js Y_lt + Y_js X_lt) D_kt
429  // or
430  // A(i+nd*j,k+nd*l) +=
431  // \sum_st w [ (D X^t)_ij (D Y^t)_kl + (D Y^t)_ij (D X^t)_kl ]
432  void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y,
433  scalar_t *A)
434  {
435  Eval_DZt(X, &DXt);
436  Eval_DZt(Y, &DYt);
437  const int nd = D_height;
438  const int ah = 2*nd;
439 
440  for (int i = 0; i < ah; i++)
441  {
442  const scalar_t axi = w*DXt[i], ayi = w*DYt[i];
443  A[i+ah*i] += 2*axi*DYt[i];
444  for (int j = 0; j < i; j++)
445  {
446  const scalar_t A_ij = axi*DYt[j] + ayi*DXt[j];
447  A[i+ah*j] += A_ij;
448  A[j+ah*i] += A_ij;
449  }
450  }
451  }
452 
453  // Assemble the contribution from the term: T_ijkl = X_ij X_kl, where X is a
454  // pointer to a 2x2 matrix stored in column-major layout.
455  //
456  // The contribution to the matrix A is given by:
457  // A(i+nd*j,k+nd*l) += \sum_st w D_is X_js X_lt D_kt
458  // or
459  // A(i+nd*j,k+nd*l) += \sum_st w [ (D X^t)_ij (D X^t)_kl ]
460  void Assemble_TProd(scalar_t w, const scalar_t *X, scalar_t *A)
461  {
462  Eval_DZt(X, &DXt);
463  const int nd = D_height;
464  const int ah = 2*nd;
465 
466  for (int i = 0; i < ah; i++)
467  {
468  const scalar_t axi = w*DXt[i];
469  A[i+ah*i] += axi*DXt[i];
470  for (int j = 0; j < i; j++)
471  {
472  const scalar_t A_ij = axi*DXt[j];
473  A[i+ah*j] += A_ij;
474  A[j+ah*i] += A_ij;
475  }
476  }
477  }
478 };
479 
480 
481 /** @brief Auxiliary class for evaluating the 3x3 matrix invariants and their
482  first and second derivatives. */
483 /**
484  The type `scalar_t` must support the standard operations:
485 
486  =, +=, -=, +, -, *, /, unary -, int*scalar_t, int/scalar_t, scalar_t/int
487 
488  The type `scalar_ops` must define the static methods:
489 
490  scalar_t sign(const scalar_t &);
491  scalar_t pow(const scalar_t &x, int a, int b); // x^(a/b)
492 */
493 template <typename scalar_t, typename scalar_ops = ScalarOps<scalar_t> >
494 class InvariantsEvaluator3D
495 {
496 protected:
497  // Transformation Jacobian
498  const scalar_t *J;
499 
500  // Invatiants: I1b = det(J)^{-2/3} * ||J||_F^2
501  // I2b = det(J)^{ 2/3} * ||J^{-1}||_F^2
502  // I3b = det(J)
503  // Computed as:
504  // I_1 = ||J||_F^2, \bar{I}_1 = det(J)^{-2/3}*I_1,
505  // I_2 = (1/2)*(||J||_F^4-||J J^t||_F^2) = (1/2)*(I_1^2-||J J^t||_F^2),
506  // \bar{I}_2 = det(J)^{-4/3}*I_2,
507  // I_3 = det(J)^2, \bar{I}_3 = det(J).
508  scalar_t I1, I1b, I2, I2b, I3b;
509  scalar_t I3b_p; // I3b^{-2/3}
510 
511  // Derivatives of I1, I1b, I2, I2b, I3, and I3b using column-major storage.
512  scalar_t dI1[9], dI1b[9], dI2[9], dI2b[9], dI3[9], dI3b[9];
513  scalar_t B[6]; // B = J J^t (diagonal entries first, then off-diagonal)
514 
515  int D_height, alloc_height;
516  const scalar_t *D; // Always points to external data or is empty
517  scalar_t *DaJ, *DJt, *DdI2t, *DXt, *DYt;
518 
519  enum EvalMasks
520  {
521  HAVE_I1 = 1,
522  HAVE_I1b = 2,
523  HAVE_B_offd = 4,
524  HAVE_I2 = 8,
525  HAVE_I2b = 16,
526  HAVE_I3b = 1<<5,
527  HAVE_I3b_p = 1<<6,
528  HAVE_dI1 = 1<<7,
529  HAVE_dI1b = 1<<8,
530  HAVE_dI2 = 1<<9,
531  HAVE_dI2b = 1<<10,
532  HAVE_dI3 = 1<<11,
533  HAVE_dI3b = 1<<12,
534  HAVE_DaJ = 1<<13, // D adj(J) = D dI3b^t
535  HAVE_DJt = 1<<14, // D J^t
536  HAVE_DdI2t = 1<<15 // D dI2^t
537  };
538 
539  // Bitwise OR of EvalMasks
540  int eval_state;
541 
542  bool dont(int have_mask) const { return !(eval_state & have_mask); }
543 
544  void Eval_I1()
545  {
546  eval_state |= HAVE_I1;
547  B[0] = J[0]*J[0] + J[3]*J[3] + J[6]*J[6];
548  B[1] = J[1]*J[1] + J[4]*J[4] + J[7]*J[7];
549  B[2] = J[2]*J[2] + J[5]*J[5] + J[8]*J[8];
550  I1 = B[0] + B[1] + B[2];
551  }
552  void Eval_I1b() // det(J)^{-2/3}*I_1 = I_1/I_3^{1/3}
553  {
554  eval_state |= HAVE_I1b;
555  I1b = Get_I1()*Get_I3b_p();
556  }
557  void Eval_B_offd()
558  {
559  eval_state |= HAVE_B_offd;
560  // B = J J^t
561  // B[3]=B(0,1), B[4]=B(0,2), B[5]=B(1,2)
562  B[3] = J[0]*J[1] + J[3]*J[4] + J[6]*J[7]; // B(0,1)
563  B[4] = J[0]*J[2] + J[3]*J[5] + J[6]*J[8]; // B(0,2)
564  B[5] = J[1]*J[2] + J[4]*J[5] + J[7]*J[8]; // B(1,2)
565  }
566  void Eval_I2()
567  {
568  eval_state |= HAVE_I2;
569  Get_I1();
570  if (dont(HAVE_B_offd)) { Eval_B_offd(); }
571  const scalar_t BF2 = B[0]*B[0] + B[1]*B[1] + B[2]*B[2] +
572  2*(B[3]*B[3] + B[4]*B[4] + B[5]*B[5]);
573  I2 = (I1*I1 - BF2)/2;
574  }
575  void Eval_I2b() // I2b = I2*I3b^{-4/3}
576  {
577  eval_state |= HAVE_I2b;
578  Get_I3b_p();
579  I2b = Get_I2()*I3b_p*I3b_p;
580  }
581  void Eval_I3b() // det(J)
582  {
583  eval_state |= HAVE_I3b;
584  I3b = J[0]*(J[4]*J[8] - J[7]*J[5]) - J[1]*(J[3]*J[8] - J[5]*J[6]) +
585  J[2]*(J[3]*J[7] - J[4]*J[6]);
586  }
587  scalar_t Get_I3b_p() // I3b^{-2/3}
588  {
589  if (dont(HAVE_I3b_p))
590  {
591  eval_state |= HAVE_I3b_p;
592  const scalar_t i3b = Get_I3b();
593  I3b_p = scalar_ops::pow(i3b, -2, 3);
594  }
595  return I3b_p;
596  }
597  void Eval_dI1()
598  {
599  eval_state |= HAVE_dI1;
600  for (int i = 0; i < 9; i++)
601  {
602  dI1[i] = 2*J[i];
603  }
604  }
605  void Eval_dI1b()
606  {
607  eval_state |= HAVE_dI1b;
608  // I1b = I3b^{-2/3}*I1
609  // dI1b = 2*I3b^{-2/3}*(J - (1/3)*I1/I3b*dI3b)
610  const scalar_t c1 = 2*Get_I3b_p();
611  const scalar_t c2 = Get_I1()/(3*I3b);
612  Get_dI3b();
613  for (int i = 0; i < 9; i++)
614  {
615  dI1b[i] = c1*(J[i] - c2*dI3b[i]);
616  }
617  }
618  void Eval_dI2()
619  {
620  eval_state |= HAVE_dI2;
621  // dI2 = 2 I_1 J - 2 J J^t J = 2 (I_1 I - B) J
622  Get_I1();
623  if (dont(HAVE_B_offd)) { Eval_B_offd(); }
624  // B[0]=B(0,0), B[1]=B(1,1), B[2]=B(2,2)
625  // B[3]=B(0,1), B[4]=B(0,2), B[5]=B(1,2)
626  const scalar_t C[6] =
627  {
628  2*(I1 - B[0]), 2*(I1 - B[1]), 2*(I1 - B[2]),
629  -2*B[3], -2*B[4], -2*B[5]
630  };
631  // | C[0] C[3] C[4] | | J[0] J[3] J[6] |
632  // dI2 = | C[3] C[1] C[5] | | J[1] J[4] J[7] |
633  // | C[4] C[5] C[2] | | J[2] J[5] J[8] |
634  dI2[0] = C[0]*J[0] + C[3]*J[1] + C[4]*J[2];
635  dI2[1] = C[3]*J[0] + C[1]*J[1] + C[5]*J[2];
636  dI2[2] = C[4]*J[0] + C[5]*J[1] + C[2]*J[2];
637 
638  dI2[3] = C[0]*J[3] + C[3]*J[4] + C[4]*J[5];
639  dI2[4] = C[3]*J[3] + C[1]*J[4] + C[5]*J[5];
640  dI2[5] = C[4]*J[3] + C[5]*J[4] + C[2]*J[5];
641 
642  dI2[6] = C[0]*J[6] + C[3]*J[7] + C[4]*J[8];
643  dI2[7] = C[3]*J[6] + C[1]*J[7] + C[5]*J[8];
644  dI2[8] = C[4]*J[6] + C[5]*J[7] + C[2]*J[8];
645  }
646  void Eval_dI2b()
647  {
648  eval_state |= HAVE_dI2b;
649  // I2b = det(J)^{-4/3}*I2 = I3b^{-4/3}*I2
650  // dI2b = (-4/3)*I3b^{-7/3}*I2*dI3b + I3b^{-4/3}*dI2
651  // = I3b^{-4/3} * [ dI2 - (4/3)*I2/I3b*dI3b ]
652  Get_I3b_p();
653  const scalar_t c1 = I3b_p*I3b_p;
654  const scalar_t c2 = (4*Get_I2()/I3b)/3;
655  Get_dI2();
656  Get_dI3b();
657  for (int i = 0; i < 9; i++)
658  {
659  dI2b[i] = c1*(dI2[i] - c2*dI3b[i]);
660  }
661  }
662  void Eval_dI3()
663  {
664  eval_state |= HAVE_dI3;
665  // I3 = I3b^2
666  // dI3 = 2*I3b*dI3b = 2*det(J)*adj(J)^T
667  const scalar_t c1 = 2*Get_I3b();
668  Get_dI3b();
669  for (int i = 0; i < 9; i++)
670  {
671  dI3[i] = c1*dI3b[i];
672  }
673  }
674  void Eval_dI3b()
675  {
676  eval_state |= HAVE_dI3b;
677  // I3b = det(J)
678  // dI3b = adj(J)^T
679  dI3b[0] = J[4]*J[8] - J[5]*J[7]; // 0 3 6
680  dI3b[1] = J[5]*J[6] - J[3]*J[8]; // 1 4 7
681  dI3b[2] = J[3]*J[7] - J[4]*J[6]; // 2 5 8
682  dI3b[3] = J[2]*J[7] - J[1]*J[8];
683  dI3b[4] = J[0]*J[8] - J[2]*J[6];
684  dI3b[5] = J[1]*J[6] - J[0]*J[7];
685  dI3b[6] = J[1]*J[5] - J[2]*J[4];
686  dI3b[7] = J[2]*J[3] - J[0]*J[5];
687  dI3b[8] = J[0]*J[4] - J[1]*J[3];
688  }
689  void Eval_DZt(const scalar_t *Z, scalar_t **DZt_ptr)
690  {
691  MFEM_ASSERT(D != NULL, "");
692  const int nd = D_height;
693  scalar_t *DZt = *DZt_ptr;
694  if (DZt == NULL) { *DZt_ptr = DZt = new scalar_t[3*alloc_height]; }
695  for (int i = 0; i < nd; i++)
696  {
697  const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
698  DZt[i0] = D[i0]*Z[0] + D[i1]*Z[3] + D[i2]*Z[6];
699  DZt[i1] = D[i0]*Z[1] + D[i1]*Z[4] + D[i2]*Z[7];
700  DZt[i2] = D[i0]*Z[2] + D[i1]*Z[5] + D[i2]*Z[8];
701  }
702  }
703  void Eval_DaJ() // DaJ = D adj(J) = D dI3b^t
704  {
705  eval_state |= HAVE_DaJ;
706  Get_dI3b();
707  Eval_DZt(dI3b, &DaJ);
708  }
709  void Eval_DJt() // DJt = D J^t
710  {
711  eval_state |= HAVE_DJt;
712  Eval_DZt(J, &DJt);
713  }
714  void Eval_DdI2t() // DdI2t = D dI2^t
715  {
716  eval_state |= HAVE_DdI2t;
717  Get_dI2();
718  Eval_DZt(dI2, &DdI2t);
719  }
720 
721 public:
722  /// The Jacobian should use column-major storage.
723  InvariantsEvaluator3D(const scalar_t *Jac = NULL)
724  : J(Jac), D_height(), alloc_height(),
725  D(), DaJ(), DJt(), DdI2t(), DXt(), DYt(), eval_state(0) { }
726 
727  ~InvariantsEvaluator3D()
728  {
729  delete [] DYt;
730  delete [] DXt;
731  delete [] DdI2t;
732  delete [] DJt;
733  delete [] DaJ;
734  }
735 
736  /// The Jacobian should use column-major storage.
737  void SetJacobian(const scalar_t *Jac) { J = Jac; eval_state = 0; }
738 
739  /// The @a Deriv matrix is `dof x 3`, using column-major storage.
740  void SetDerivativeMatrix(int height, const scalar_t *Deriv)
741  {
742  eval_state &= ~(HAVE_DaJ | HAVE_DJt | HAVE_DdI2t);
743  if (alloc_height < height)
744  {
745  delete [] DYt; DYt = NULL;
746  delete [] DXt; DXt = NULL;
747  delete [] DdI2t; DdI2t = NULL;
748  delete [] DJt; DJt = NULL;
749  delete [] DaJ; DaJ = NULL;
750  alloc_height = height;
751  }
752  D_height = height;
753  D = Deriv;
754  }
755 
756  scalar_t Get_I1() { if (dont(HAVE_I1 )) { Eval_I1(); } return I1; }
757  scalar_t Get_I1b() { if (dont(HAVE_I1b)) { Eval_I1b(); } return I1b; }
758  scalar_t Get_I2() { if (dont(HAVE_I2 )) { Eval_I2(); } return I2; }
759  scalar_t Get_I2b() { if (dont(HAVE_I2b)) { Eval_I2b(); } return I2b; }
760  scalar_t Get_I3() { if (dont(HAVE_I3b)) { Eval_I3b(); } return I3b*I3b; }
761  scalar_t Get_I3b() { if (dont(HAVE_I3b)) { Eval_I3b(); } return I3b; }
762 
763  const scalar_t *Get_dI1()
764  {
765  if (dont(HAVE_dI1 )) { Eval_dI1(); } return dI1;
766  }
767  const scalar_t *Get_dI1b()
768  {
769  if (dont(HAVE_dI1b)) { Eval_dI1b(); } return dI1b;
770  }
771  const scalar_t *Get_dI2()
772  {
773  if (dont(HAVE_dI2)) { Eval_dI2(); } return dI2;
774  }
775  const scalar_t *Get_dI2b()
776  {
777  if (dont(HAVE_dI2b)) { Eval_dI2b(); } return dI2b;
778  }
779  const scalar_t *Get_dI3()
780  {
781  if (dont(HAVE_dI3)) { Eval_dI3(); } return dI3;
782  }
783  const scalar_t *Get_dI3b()
784  {
785  if (dont(HAVE_dI3b)) { Eval_dI3b(); } return dI3b;
786  }
787 
788  // Assemble operation for tensor X with components X_jslt:
789  // A(i+nd*j,k+nd*l) += (\sum_st w D_is X_jslt D_kt)
790  // 0 <= i,k < nd, 0 <= j,l,s,t < 3
791  // where nd is the height of D, i.e. the number of DOFs in one component.
792 
793  void Assemble_ddI1(scalar_t w, scalar_t *A)
794  {
795  // ddI1_jslt = 2 I_jslt = 2 δ_jl δ_st
796  // A(i+nd*j,k+nd*l) += (\sum_st 2 w D_is δ_jl δ_st D_kt)
797  // or
798  // A(i+nd*j,k+nd*l) += (2 w) (\sum_s D_is D_ks) δ_jl
799  // A(i+nd*j,k+nd*l) += (2 w) (D D^t)_ik δ_jl
800 
801  const int nd = D_height;
802  const int ah = 3*nd;
803  const scalar_t a = 2*w;
804  for (int i = 0; i < nd; i++)
805  {
806  const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
807  const scalar_t aDi[3] = { a*D[i0], a*D[i1], a*D[i2] };
808  // k == i
809  const scalar_t aDDt_ii = aDi[0]*D[i0] + aDi[1]*D[i1] + aDi[2]*D[i2];
810  A[i0+ah*i0] += aDDt_ii;
811  A[i1+ah*i1] += aDDt_ii;
812  A[i2+ah*i2] += aDDt_ii;
813  // 0 <= k < i
814  for (int k = 0; k < i; k++)
815  {
816  const int k0 = k+nd*0, k1 = k+nd*1, k2 = k+nd*2;
817  const scalar_t aDDt_ik = aDi[0]*D[k0] + aDi[1]*D[k1] + aDi[2]*D[k2];
818  A[i0+ah*k0] += aDDt_ik;
819  A[k0+ah*i0] += aDDt_ik;
820  A[i1+ah*k1] += aDDt_ik;
821  A[k1+ah*i1] += aDDt_ik;
822  A[i2+ah*k2] += aDDt_ik;
823  A[k2+ah*i2] += aDDt_ik;
824  }
825  }
826  }
827  void Assemble_ddI1b(scalar_t w, scalar_t *A)
828  {
829  // Similar to InvariantsEvaluator2D::Assemble_ddI1b():
830  //
831  // ddI1b = X1 + X2 + X3, where
832  // X1_ijkl = (2/3*I1b/I3) [ (2/3 δ_ks δ_it + δ_kt δ_si) dI3b_tj dI3b_sl ]
833  // = (2/3*I1b/I3) [ 2/3 dI3b_ij dI3b_kl + dI3b_kj dI3b_il ]
834  // X2_ijkl = (2*I3b^{-2/3}) δ_ik δ_jl = (I3b^{-2/3}) ddI1_ijkl
835  // X3_ijkl = -(4/3*I3b^{-5/3}) (δ_ks δ_it) (J_tj dI3b_sl + dI3b_tj J_sl)
836  // = -(4/3*I3b^{-5/3}) (J_ij dI3b_kl + dI3b_ij J_kl)
837  //
838  // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI1b_jslt D_kt)
839  // or
840  // A(i+nd*j,k+nd*l) +=
841  // w (2/3*I1b/I3) [ 2/3 DaJ_ij DaJ_kl + DaJ_il DaJ_kj ]
842  // + w (2*I3b^{-2/3}) (D D^t)_ik δ_jl
843  // - w (4/3*I3b^{-5/3}) [ DJt_ij DaJ_kl + DaJ_ij DJt_kl ]
844 
845  if (dont(HAVE_DaJ)) { Eval_DaJ(); }
846  if (dont(HAVE_DJt)) { Eval_DJt(); }
847  const int nd = D_height;
848  const int ah = 3*nd;
849  const scalar_t r23 = scalar_t(2)/3;
850  const scalar_t r53 = scalar_t(5)/3;
851  const scalar_t a = r23*w*Get_I1b()/Get_I3();
852  const scalar_t b = 2*w*Get_I3b_p();
853  const scalar_t c = -r23*b/I3b;
854  for (int i = 0; i < nd; i++)
855  {
856  // A1a_ik_jl = 2/3 a DaJ_ij DaJ_kl, A1b_ik_jl = a DaJ_il DaJ_kj
857  // Symmetries: A1a_ik_jl = A1a_ki_lj = 2/3 A1b_ik_lj = 2/3 A1b_ki_jl
858  // A1_ik_jl = A1_ki_lj = A1b_ik_jl + 2/3 A1b_ik_lj
859  // A1_ik_lj = A1_ki_jl = 2/3 A1b_ik_jl + A1b_ik_lj
860  // k == i:
861  // A1_ii_jl = A1_ii_lj = (5/3) a DaJ_ij DaJ_il
862  // l == j:
863  // A1_ik_jj = A1_ki_jj = (5/3) a DaJ_ij DaJ_kj
864  // k == i && l == j:
865  // A1_ii_jj = (5/3) a DaJ_ij^2
866 
867  // A2_ik_jl = b (D D^t)_ik δ_jl
868  // Symmetries:
869 
870  // A3_ik_jl = c [ DJt_ij DaJ_kl + DaJ_ij DJt_kl ]
871  // Symmetries:
872  // A3_ik_jl = A3_ki_lj = c [ DJt_ij DaJ_kl + DaJ_ij DJt_kl ]
873  // A3_ik_lj = A3_ki_jl = c [ DJt_il DaJ_kj + DaJ_il DJt_kj ]
874  // k == i:
875  // A3_ii_jl = A3_ii_lj = c [ DJt_ij DaJ_il + DaJ_ij DJt_il ]
876  // l == j:
877  // A3_ik_jj = A3_ki_jj = c [ DJt_ij DaJ_kj + DaJ_ij DJt_kj ]
878  // k == i && l == j:
879  // A3_ii_jj = 2 c DJt_ij DaJ_ij
880 
881  const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
882  const scalar_t aDaJ_i[3] = { a*DaJ[i0], a*DaJ[i1], a*DaJ[i2] };
883  const scalar_t bD_i[3] = { b*D[i0], b*D[i1], b*D[i2] };
884  const scalar_t cDJt_i[3] = { c*DJt[i0], c*DJt[i1], c*DJt[i2] };
885  const scalar_t cDaJ_i[3] = { c*DaJ[i0], c*DaJ[i1], c*DaJ[i2] };
886  // k == i
887  {
888  // Symmetries: A2_ii_00 = A2_ii_11 = A2_ii_22
889  const scalar_t A2_ii = bD_i[0]*D[i0]+bD_i[1]*D[i1]+bD_i[2]*D[i2];
890  A[i0+ah*i0] += (r53*aDaJ_i[0] + 2*cDJt_i[0])*DaJ[i0] + A2_ii;
891  A[i1+ah*i1] += (r53*aDaJ_i[1] + 2*cDJt_i[1])*DaJ[i1] + A2_ii;
892  A[i2+ah*i2] += (r53*aDaJ_i[2] + 2*cDJt_i[2])*DaJ[i2] + A2_ii;
893 
894  // Symmetries: A_ii_jl = A_ii_lj
895  for (int j = 1; j < 3; j++)
896  {
897  const int ij = i+nd*j;
898  for (int l = 0; l < j; l++)
899  {
900  const int il = i+nd*l;
901  const scalar_t A_ii_jl =
902  (r53*aDaJ_i[j] + cDJt_i[j])*DaJ[il] + cDaJ_i[j]*DJt[il];
903  A[ij+ah*il] += A_ii_jl;
904  A[il+ah*ij] += A_ii_jl;
905  }
906  }
907  }
908  // 0 <= k < i
909  for (int k = 0; k < i; k++)
910  {
911  const int k0 = k+nd*0, k1 = k+nd*1, k2 = k+nd*2;
912  // Symmetries: A2_ik_jj = A2_ki_ll
913  const scalar_t A2_ik = bD_i[0]*D[k0]+bD_i[1]*D[k1]+bD_i[2]*D[k2];
914 
915  // l == j
916  for (int j = 0; j < 3; j++)
917  {
918  const int ij = i+nd*j, kj = k+nd*j;
919  const scalar_t A_ik_jj = (r53*aDaJ_i[j] + cDJt_i[j])*DaJ[kj] +
920  cDaJ_i[j]*DJt[kj] + A2_ik;
921  A[ij+ah*kj] += A_ik_jj;
922  A[kj+ah*ij] += A_ik_jj;
923  }
924 
925  // 0 <= l < j
926  for (int j = 1; j < 3; j++)
927  {
928  const int ij = i+nd*j, kj = k+nd*j;
929  for (int l = 0; l < j; l++)
930  {
931  const int il = i+nd*l, kl = k+nd*l;
932  // A1b_ik_jl = a DaJ_il DaJ_kj
933  const scalar_t A1b_ik_jl = aDaJ_i[l]*DaJ[kj];
934  const scalar_t A1b_ik_lj = aDaJ_i[j]*DaJ[kl];
935  // A1_ik_jl = A1_ki_lj = A1b_ik_jl + 2/3 A1b_ik_lj
936  // A1_ik_lj = A1_ki_jl = 2/3 A1b_ik_jl + A1b_ik_lj
937  // A3_ik_jl = c [ DJt_ij DaJ_kl + DaJ_ij DJt_kl ]
938  const scalar_t A_ik_jl = A1b_ik_jl + r23*A1b_ik_lj +
939  cDJt_i[j]*DaJ[kl]+cDaJ_i[j]*DJt[kl];
940  A[ij+ah*kl] += A_ik_jl;
941  A[kl+ah*ij] += A_ik_jl;
942  const scalar_t A_ik_lj = r23*A1b_ik_jl + A1b_ik_lj +
943  cDJt_i[l]*DaJ[kj]+cDaJ_i[l]*DJt[kj];
944  A[il+ah*kj] += A_ik_lj;
945  A[kj+ah*il] += A_ik_lj;
946  }
947  }
948  }
949  }
950  }
951  void Assemble_ddI2(scalar_t w, scalar_t *A)
952  {
953  // dI2 = 2 I_1 J - 2 J J^t J = 2 (I_1 I - B) J
954  //
955  // ddI2 = X1 + X2 + X3
956  // X1_ijkl = (2 I_1) δ_ik δ_jl
957  // X2_ijkl = 2 ( 2 δ_ku δ_iv - δ_ik δ_uv - δ_kv δ_iu ) J_vj J_ul
958  // X3_ijkl = -2 (J J^t)_ik δ_jl = -2 B_ik δ_jl
959  //
960  // Apply: j->s, i->j, l->t, k->l
961  // X2_jslt = 2 ( δ_lu δ_jv - δ_jl δ_uv +
962  // δ_lu δ_jv - δ_lv δ_ju ) J_vs J_ut
963  //
964  // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI2_jslt D_kt)
965  //
966  // \sum_st w D_is X1_jslt D_kt =
967  // \sum_st w D_is [ (2 I_1) δ_jl δ_st ] D_kt =
968  // (2 w I_1) D_is δ_jl D_ks = (2 w I_1) (D D^t)_ik δ_jl
969  //
970  // \sum_st w D_is X2_jslt D_kt =
971  // \sum_stuv w D_is [ 2 ( δ_lu δ_jv - δ_jl δ_uv +
972  // δ_lu δ_jv - δ_lv δ_ju ) J_vs J_ut ] D_kt =
973  // \sum_uv 2 w [ δ_lu δ_jv - δ_jl δ_uv +
974  // δ_lu δ_jv - δ_lv δ_ju ] (D J^t)_iv (D J^t)_ku =
975  // 2 w ( DJt_ij DJt_kl - δ_jl (DJt DJt^t)_ik ) +
976  // 2 w ( DJt_ij DJt_kl - DJt_il DJt_kj )
977  //
978  // \sum_st w D_is X3_jslt D_kt = \sum_st w D_is [ -2 B_jl δ_st ] D_kt =
979  // -2 w (D D^t)_ik B_jl
980  //
981  // A(i+nd*j,k+nd*l) +=
982  // (2 w I_1) (D D^t)_ik δ_jl - 2 w (D D^t)_ik B_jl +
983  // 2 w DJt_ij DJt_kl - 2 w (DJt DJt^t)_ik δ_jl +
984  // 2 w ( DJt_ij DJt_kl - DJt_il DJt_kj )
985  //
986  // The last term is a determinant: rows {i,k} and columns {j,l} of DJt:
987  // | DJt_ij DJt_il |
988  // | DJt_kj DJt_kl | = DJt_ij DJt_kl - DJt_il DJt_kj
989 
990  if (dont(HAVE_DJt)) { Eval_DJt(); }
991  Get_I1(); // evaluates I1 and the diagonal of B
992  if (dont(HAVE_B_offd)) { Eval_B_offd(); }
993  const int nd = D_height;
994  const int ah = 3*nd;
995  const scalar_t a = 2*w;
996  for (int i = 0; i < ah; i++)
997  {
998  const scalar_t avi = a*DJt[i];
999  A[i+ah*i] += avi*DJt[i];
1000  for (int j = 0; j < i; j++)
1001  {
1002  const scalar_t aVVt_ij = avi*DJt[j];
1003  A[i+ah*j] += aVVt_ij;
1004  A[j+ah*i] += aVVt_ij;
1005  }
1006  }
1007 
1008  for (int i = 0; i < nd; i++)
1009  {
1010  const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
1011  const scalar_t aD_i[3] = { a*D[i0], a*D[i1], a*D[i2] };
1012  const scalar_t aDJt_i[3] = { a*DJt[i0], a*DJt[i1], a*DJt[i2] };
1013  // k == i
1014  {
1015  const scalar_t aDDt_ii =
1016  aD_i[0]*D[i0] + aD_i[1]*D[i1] + aD_i[2]*D[i2];
1017  const scalar_t Z1_ii =
1018  I1*aDDt_ii - (aDJt_i[0]*DJt[i0] + aDJt_i[1]*DJt[i1] +
1019  aDJt_i[2]*DJt[i2]);
1020  // l == j
1021  for (int j = 0; j < 3; j++)
1022  {
1023  const int ij = i+nd*j;
1024  A[ij+ah*ij] += Z1_ii - aDDt_ii*B[j];
1025  }
1026  // l != j
1027  const scalar_t Z2_ii_01 = aDDt_ii*B[3];
1028  const scalar_t Z2_ii_02 = aDDt_ii*B[4];
1029  const scalar_t Z2_ii_12 = aDDt_ii*B[5];
1030  A[i0+ah*i1] -= Z2_ii_01;
1031  A[i1+ah*i0] -= Z2_ii_01;
1032  A[i0+ah*i2] -= Z2_ii_02;
1033  A[i2+ah*i0] -= Z2_ii_02;
1034  A[i1+ah*i2] -= Z2_ii_12;
1035  A[i2+ah*i1] -= Z2_ii_12;
1036  }
1037  // 0 <= k < i
1038  for (int k = 0; k < i; k++)
1039  {
1040  const int k0 = k+nd*0, k1 = k+nd*1, k2 = k+nd*2;
1041  const scalar_t aDDt_ik =
1042  aD_i[0]*D[k0] + aD_i[1]*D[k1] + aD_i[2]*D[k2];
1043  const scalar_t Z1_ik =
1044  I1*aDDt_ik - (aDJt_i[0]*DJt[k0] + aDJt_i[1]*DJt[k1] +
1045  aDJt_i[2]*DJt[k2]);
1046  // l == j
1047  for (int j = 0; j < 3; j++)
1048  {
1049  const int ij = i+nd*j, kj = k+nd*j;
1050  const scalar_t Z2_ik_jj = Z1_ik - aDDt_ik*B[j];
1051  A[ij+ah*kj] += Z2_ik_jj;
1052  A[kj+ah*ij] += Z2_ik_jj;
1053  }
1054  // l != j
1055  {
1056  const scalar_t Z2_ik_01 = aDDt_ik*B[3];
1057  A[i0+ah*k1] -= Z2_ik_01;
1058  A[i1+ah*k0] -= Z2_ik_01;
1059  A[k0+ah*i1] -= Z2_ik_01;
1060  A[k1+ah*i0] -= Z2_ik_01;
1061  const scalar_t Z2_ik_02 = aDDt_ik*B[4];
1062  A[i0+ah*k2] -= Z2_ik_02;
1063  A[i2+ah*k0] -= Z2_ik_02;
1064  A[k0+ah*i2] -= Z2_ik_02;
1065  A[k2+ah*i0] -= Z2_ik_02;
1066  const scalar_t Z2_ik_12 = aDDt_ik*B[5];
1067  A[i1+ah*k2] -= Z2_ik_12;
1068  A[i2+ah*k1] -= Z2_ik_12;
1069  A[k1+ah*i2] -= Z2_ik_12;
1070  A[k2+ah*i1] -= Z2_ik_12;
1071  }
1072  // 0 <= l < j
1073  for (int j = 1; j < 3; j++)
1074  {
1075  const int ij = i+nd*j, kj = k+nd*j;
1076  for (int l = 0; l < j; l++)
1077  {
1078  const int il = i+nd*l, kl = k+nd*l;
1079  const scalar_t Z3_ik_jl =
1080  aDJt_i[j]*DJt[kl] - aDJt_i[l]*DJt[kj];
1081  A[ij+ah*kl] += Z3_ik_jl;
1082  A[kl+ah*ij] += Z3_ik_jl;
1083  A[il+ah*kj] -= Z3_ik_jl;
1084  A[kj+ah*il] -= Z3_ik_jl;
1085  }
1086  }
1087  }
1088  }
1089  }
1090  void Assemble_ddI2b(scalar_t w, scalar_t *A)
1091  {
1092  // dI2b = (-4/3)*I3b^{-7/3}*I2*dI3b + I3b^{-4/3}*dI2
1093  // = I3b^{-4/3} * [ dI2 - (4/3)*I2/I3b*dI3b ]
1094  //
1095  // ddI2b = X1 + X2 + X3
1096  // X1_ijkl = 16/9 det(J)^{-10/3} I2 dI3b_ij dI3b_kl +
1097  // 4/3 det(J)^{-10/3} I2 dI3b_il dI3b_kj
1098  // X2_ijkl = -4/3 det(J)^{-7/3} (dI2_ij dI3b_kl + dI2_kl dI3b_ij)
1099  // X3_ijkl = det(J)^{-4/3} ddI2_ijkl
1100  //
1101  // Apply: j->s, i->j, l->t, k->l
1102  // X1_jslt = 16/9 det(J)^{-10/3} I2 dI3b_js dI3b_lt +
1103  // 4/3 det(J)^{-10/3} I2 dI3b_jt dI3b_ls
1104  // X2_jslt = -4/3 det(J)^{-7/3} (dI2_js dI3b_lt + dI2_lt dI3b_js)
1105  //
1106  // A(i+nd*j,k+nd*l) += (\sum_st w D_is ddI2b_jslt D_kt)
1107  //
1108  // (\sum_st w D_is X1_jslt D_kt) =
1109  // 16/9 w det(J)^{-10/3} I2 DaJ_ij DaJ_kl +
1110  // 4/3 w det(J)^{-10/3} I2 DaJ_il DaJ_kj
1111  //
1112  // (\sum_st w D_is X1_jslt D_kt) =
1113  // -4/3 w det(J)^{-7/3} D_is (dI2_js dI3b_lt + dI2_lt dI3b_js) D_kt =
1114  // -4/3 w det(J)^{-7/3} [ (D dI2^t)_ij DaJ_kl + DaJ_ij (D dI2^t)_kl ]
1115  //
1116  // A(i+nd*j,k+nd*l) +=
1117  // 16/9 w det(J)^{-10/3} I2 DaJ_ij DaJ_kl +
1118  // 4/3 w det(J)^{-10/3} I2 DaJ_il DaJ_kj -
1119  // 4/3 w det(J)^{-7/3} [ DdI2t_ij DaJ_kl + DaJ_ij DdI2t_kl ] +
1120  // w det(J)^{-4/3} D_is D_kt ddI2_jslt
1121 
1122  Get_I3b_p(); // = det(J)^{-2/3}, evaluates I3b
1123  if (dont(HAVE_DaJ)) { Eval_DaJ(); }
1124  if (dont(HAVE_DdI2t)) { Eval_DdI2t(); }
1125  const int nd = D_height;
1126  const int ah = 3*nd;
1127  const scalar_t a = w*I3b_p*I3b_p;
1128  const scalar_t b = (-4*a)/(3*I3b);
1129  const scalar_t c = -b*Get_I2()/I3b;
1130  const scalar_t d = (4*c)/3;
1131 
1132  for (int i = 0; i < ah; i++)
1133  {
1134  const scalar_t dvi = d*DaJ[i];
1135  A[i+ah*i] += dvi*DaJ[i];
1136  for (int j = 0; j < i; j++)
1137  {
1138  const scalar_t dVVt_ij = dvi*DaJ[j];
1139  A[i+ah*j] += dVVt_ij;
1140  A[j+ah*i] += dVVt_ij;
1141  }
1142  }
1143  Assemble_ddI2(a, A);
1144  for (int i = 0; i < nd; i++)
1145  {
1146  const int i0 = i+nd*0, i1 = i+nd*1, i2 = i+nd*2;
1147  const scalar_t cDaJ_i[3] = { c*DaJ[i0], c*DaJ[i1], c*DaJ[i2] };
1148  const scalar_t bDaJ_i[3] = { b*DaJ[i0], b*DaJ[i1], b*DaJ[i2] };
1149  const scalar_t bDdI2t_i[3] = { b*DdI2t[i0], b*DdI2t[i1], b*DdI2t[i2] };
1150  // k == i
1151  {
1152  // l == j
1153  for (int j = 0; j < 3; j++)
1154  {
1155  const int ij = i+nd*j;
1156  A[ij+ah*ij] += (cDaJ_i[j] + 2*bDdI2t_i[j])*DaJ[ij];
1157  }
1158  // 0 <= l < j
1159  for (int j = 1; j < 3; j++)
1160  {
1161  const int ij = i+nd*j;
1162  for (int l = 0; l < j; l++)
1163  {
1164  const int il = i+nd*l;
1165  const scalar_t Z_ii_jl =
1166  (cDaJ_i[l] + bDdI2t_i[l])*DaJ[ij] + bDdI2t_i[j]*DaJ[il];
1167  A[ij+ah*il] += Z_ii_jl;
1168  A[il+ah*ij] += Z_ii_jl;
1169  }
1170  }
1171  }
1172  // 0 <= k < i
1173  for (int k = 0; k < i; k++)
1174  {
1175  // l == j
1176  for (int j = 0; j < 3; j++)
1177  {
1178  const int ij = i+nd*j, kj = k+nd*j;
1179  const scalar_t Z_ik_jj =
1180  (cDaJ_i[j] + bDdI2t_i[j])*DaJ[kj] + bDaJ_i[j]*DdI2t[kj];
1181  A[ij+ah*kj] += Z_ik_jj;
1182  A[kj+ah*ij] += Z_ik_jj;
1183  }
1184  // 0 <= l < j
1185  for (int j = 1; j < 3; j++)
1186  {
1187  const int ij = i+nd*j, kj = k+nd*j;
1188  for (int l = 0; l < j; l++)
1189  {
1190  const int il = i+nd*l, kl = k+nd*l;
1191  const scalar_t Z_ik_jl = cDaJ_i[l]*DaJ[kj] +
1192  bDdI2t_i[j]*DaJ[kl] +
1193  bDaJ_i[j]*DdI2t[kl];
1194  A[ij+ah*kl] += Z_ik_jl;
1195  A[kl+ah*ij] += Z_ik_jl;
1196  const scalar_t Z_ik_lj = cDaJ_i[j]*DaJ[kl] +
1197  bDdI2t_i[l]*DaJ[kj] +
1198  bDaJ_i[l]*DdI2t[kj];
1199  A[il+ah*kj] += Z_ik_lj;
1200  A[kj+ah*il] += Z_ik_lj;
1201  }
1202  }
1203  }
1204  }
1205  }
1206  void Assemble_ddI3(scalar_t w, scalar_t *A)
1207  {
1208  // Similar to InvariantsEvaluator2D::Assemble_ddI2():
1209  //
1210  // A(i+nd*j,k+nd*l) += 2 w [ 2 DaJ_ij DaJ_kl - DaJ_il DaJ_kj ]
1211  //
1212  // Note: the expression ( DaJ_ij DaJ_kl - DaJ_il DaJ_kj ) is the
1213  // determinant of the 2x2 matrix formed by rows {i,k} and columns {j,l}
1214  // from the matrix DaJ = D dI3b^t.
1215 
1216  if (dont(HAVE_DaJ)) { Eval_DaJ(); }
1217  const int nd = D_height;
1218  const int ah = 3*nd;
1219  const scalar_t a = 2*w;
1220 
1221  for (int i = 0; i < ah; i++)
1222  {
1223  const scalar_t avi = a*DaJ[i];
1224  A[i+ah*i] += avi*DaJ[i];
1225  for (int j = 0; j < i; j++)
1226  {
1227  const scalar_t aVVt_ij = avi*DaJ[j];
1228  A[i+ah*j] += aVVt_ij;
1229  A[j+ah*i] += aVVt_ij;
1230  }
1231  }
1232  for (int j = 1; j < 3; j++)
1233  {
1234  for (int l = 0; l < j; l++)
1235  {
1236  for (int i = 0; i < nd; i++)
1237  {
1238  const int ij = i+nd*j, il = i+nd*l;
1239  const scalar_t aDaJ_ij = a*DaJ[ij], aDaJ_il = a*DaJ[il];
1240  for (int k = 0; k < i; k++)
1241  {
1242  const int kj = k+nd*j, kl = k+nd*l;
1243  const scalar_t A_ijkl = aDaJ_ij*DaJ[kl] - aDaJ_il*DaJ[kj];
1244  A[ij+ah*kl] += A_ijkl;
1245  A[kl+ah*ij] += A_ijkl;
1246  A[kj+ah*il] -= A_ijkl;
1247  A[il+ah*kj] -= A_ijkl;
1248  }
1249  }
1250  }
1251  }
1252  }
1253  void Assemble_ddI3b(scalar_t w, scalar_t *A)
1254  {
1255  // Similar to InvariantsEvaluator2D::Assemble_ddI2b():
1256  //
1257  // A(i+nd*j,k+nd*l) += (w/I3b) [ DaJ_ij DaJ_kl - DaJ_il DaJ_kj ]
1258  //
1259  // | DaJ_ij DaJ_il | = determinant of rows {i,k}, columns {j,l} from DaJ
1260  // | DaJ_kj DaJ_kl |
1261 
1262  if (dont(HAVE_DaJ)) { Eval_DaJ(); }
1263  const int nd = D_height;
1264  const int ah = 3*nd;
1265  const scalar_t a = w/Get_I3b();
1266  for (int j = 1; j < 3; j++)
1267  {
1268  for (int l = 0; l < j; l++)
1269  {
1270  for (int i = 0; i < nd; i++)
1271  {
1272  const int ij = i+nd*j, il = i+nd*l;
1273  const scalar_t aDaJ_ij = a*DaJ[ij], aDaJ_il = a*DaJ[il];
1274  for (int k = 0; k < i; k++)
1275  {
1276  const int kj = k+nd*j, kl = k+nd*l;
1277  const scalar_t A_ijkl = aDaJ_ij*DaJ[kl] - aDaJ_il*DaJ[kj];
1278  A[ij+ah*kl] += A_ijkl;
1279  A[kl+ah*ij] += A_ijkl;
1280  A[kj+ah*il] -= A_ijkl;
1281  A[il+ah*kj] -= A_ijkl;
1282  }
1283  }
1284  }
1285  }
1286  }
1287  // Assemble the contribution from the term: T_ijkl = X_ij Y_kl + Y_ij X_kl,
1288  // where X and Y are pointers to 3x3 matrices stored in column-major layout.
1289  //
1290  // The contribution to the matrix A is given by:
1291  // A(i+nd*j,k+nd*l) += \sum_st w D_is T_jslt D_kt
1292  // or
1293  // A(i+nd*j,k+nd*l) += \sum_st w D_is (X_js Y_lt + Y_js X_lt) D_kt
1294  // or
1295  // A(i+nd*j,k+nd*l) +=
1296  // \sum_st w [ (D X^t)_ij (D Y^t)_kl + (D Y^t)_ij (D X^t)_kl ]
1297  void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y,
1298  scalar_t *A)
1299  {
1300  Eval_DZt(X, &DXt);
1301  Eval_DZt(Y, &DYt);
1302  const int nd = D_height;
1303  const int ah = 3*nd;
1304 
1305  for (int i = 0; i < ah; i++)
1306  {
1307  const scalar_t axi = w*DXt[i], ayi = w*DYt[i];
1308  A[i+ah*i] += 2*axi*DYt[i];
1309  for (int j = 0; j < i; j++)
1310  {
1311  const scalar_t A_ij = axi*DYt[j] + ayi*DXt[j];
1312  A[i+ah*j] += A_ij;
1313  A[j+ah*i] += A_ij;
1314  }
1315  }
1316  }
1317  // Assemble the contribution from the term: T_ijkl = X_ij X_kl, where X is a
1318  // pointer to a 3x3 matrix stored in column-major layout.
1319  //
1320  // The contribution to the matrix A is given by:
1321  // A(i+nd*j,k+nd*l) += \sum_st w D_is X_js X_lt D_kt
1322  // or
1323  // A(i+nd*j,k+nd*l) += \sum_st w [ (D X^t)_ij (D X^t)_kl ]
1324  void Assemble_TProd(scalar_t w, const scalar_t *X, scalar_t *A)
1325  {
1326  Eval_DZt(X, &DXt);
1327  const int nd = D_height;
1328  const int ah = 3*nd;
1329 
1330  for (int i = 0; i < ah; i++)
1331  {
1332  const scalar_t axi = w*DXt[i];
1333  A[i+ah*i] += axi*DXt[i];
1334  for (int j = 0; j < i; j++)
1335  {
1336  const scalar_t A_ij = axi*DXt[j];
1337  A[i+ah*j] += A_ij;
1338  A[j+ah*i] += A_ij;
1339  }
1340  }
1341  }
1342 };
1343 
1344 }
1345 
1346 #endif
mfem::InvariantsEvaluator3D::dont
bool dont(int have_mask) const
Definition: invariants.hpp:542
mfem::InvariantsEvaluator3D::Assemble_ddI3
void Assemble_ddI3(scalar_t w, scalar_t *A)
Definition: invariants.hpp:1206
mfem::InvariantsEvaluator2D::Assemble_ddI2b
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:384
mfem::InvariantsEvaluator2D::InvariantsEvaluator2D
InvariantsEvaluator2D(const scalar_t *Jac=NULL)
The Jacobian should use column-major storage.
Definition: invariants.hpp:168
mfem::ScalarOps::pow
static scalar_t pow(const scalar_t &x, int m, int n)
Definition: invariants.hpp:32
mfem::InvariantsEvaluator3D::Assemble_ddI2b
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:1090
mfem::InvariantsEvaluator2D::SetDerivativeMatrix
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 2, using column-major storage.
Definition: invariants.hpp:184
mfem::InvariantsEvaluator3D::Get_dI3b
const scalar_t * Get_dI3b()
Definition: invariants.hpp:783
mfem::ScalarOps::sign
static scalar_t sign(const scalar_t &a)
Definition: invariants.hpp:29
mfem::InvariantsEvaluator3D::Eval_DZt
void Eval_DZt(const scalar_t *Z, scalar_t **DZt_ptr)
Definition: invariants.hpp:689
mfem::InvariantsEvaluator2D::dont
bool dont(int have_mask) const
Definition: invariants.hpp:81
mfem::InvariantsEvaluator2D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition: invariants.hpp:432
mfem::InvariantsEvaluator2D::Assemble_ddI1
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:226
mfem::InvariantsEvaluator2D
Auxiliary class for evaluating the 2x2 matrix invariants and their first and second derivatives...
Definition: invariants.hpp:49
mfem::InvariantsEvaluator2D::Assemble_ddI2
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:335
mfem::InvariantsEvaluator2D::Eval_DZt
void Eval_DZt(const scalar_t *Z, scalar_t **DZt_ptr)
Definition: invariants.hpp:152
mfem::InvariantsEvaluator2D::SetJacobian
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:181
mfem::InvariantsEvaluator3D::Assemble_ddI3b
void Assemble_ddI3b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:1253
mfem::ScalarOps
Auxiliary class used as the default for the second template parameter in the classes InvariantsEvalua...
Definition: invariants.hpp:27
b
double b
Definition: lissajous.cpp:42
mfem::InvariantsEvaluator2D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, scalar_t *A)
Definition: invariants.hpp:460
mfem::InvariantsEvaluator2D::Get_dI2b
const scalar_t * Get_dI2b()
Definition: invariants.hpp:216
mfem::InvariantsEvaluator3D::SetDerivativeMatrix
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 3, using column-major storage.
Definition: invariants.hpp:740
mfem::InvariantsEvaluator3D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, scalar_t *A)
Definition: invariants.hpp:1324
mfem::InvariantsEvaluator2D::Assemble_ddI1b
void Assemble_ddI1b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:257
mfem::InvariantsEvaluator3D::Get_dI3
const scalar_t * Get_dI3()
Definition: invariants.hpp:779
mfem::InvariantsEvaluator3D::Assemble_ddI2
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:951
mfem::InvariantsEvaluator3D::Get_dI1
const scalar_t * Get_dI1()
Definition: invariants.hpp:763
mfem::InvariantsEvaluator2D::Get_dI1
const scalar_t * Get_dI1()
Definition: invariants.hpp:204
mfem::InvariantsEvaluator2D::Get_dI1b
const scalar_t * Get_dI1b()
Definition: invariants.hpp:208
mfem::InvariantsEvaluator3D::Get_dI1b
const scalar_t * Get_dI1b()
Definition: invariants.hpp:767
a
double a
Definition: lissajous.cpp:41
mfem::InvariantsEvaluator3D::Get_dI2b
const scalar_t * Get_dI2b()
Definition: invariants.hpp:775
mfem::InvariantsEvaluator3D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition: invariants.hpp:1297
mfem::InvariantsEvaluator3D::Get_dI2
const scalar_t * Get_dI2()
Definition: invariants.hpp:771
mfem::InvariantsEvaluator3D::Assemble_ddI1
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:793
mfem::InvariantsEvaluator3D
Auxiliary class for evaluating the 3x3 matrix invariants and their first and second derivatives...
Definition: invariants.hpp:494
mfem::InvariantsEvaluator2D::Get_dI2
const scalar_t * Get_dI2()
Definition: invariants.hpp:212
mfem::InvariantsEvaluator3D::Assemble_ddI1b
void Assemble_ddI1b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:827
mfem::InvariantsEvaluator3D::InvariantsEvaluator3D
InvariantsEvaluator3D(const scalar_t *Jac=NULL)
The Jacobian should use column-major storage.
Definition: invariants.hpp:723
mfem::InvariantsEvaluator3D::SetJacobian
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:737
Generated on Wed Dec 7 2022 10:41:47 for MFEM by   doxygen 1.8.5