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