MFEM  v4.6.0
Finite element discretization library
tmop_pa_h3s.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../tmop.hpp"
13 #include "tmop_pa.hpp"
14 #include "../../general/forall.hpp"
15 #include "../../linalg/kernels.hpp"
16 #include "../../linalg/dinvariants.hpp"
17 
18 namespace mfem
19 {
20 
21 using Args = kernels::InvariantsEvaluator3D::Buffers;
22 
23 // dP_302 = (dI2b*dI1b + dI1b*dI2b)/9 + (I1b/9)*ddI2b + (I2b/9)*ddI1b
24 static MFEM_HOST_DEVICE inline
25 void EvalH_302(const int e, const int qx, const int qy, const int qz,
26  const double weight, const double *J, DeviceTensor<8,double> dP,
27  double *B, double *dI1b, double *ddI1b,
28  double *dI2, double *dI2b, double *ddI2, double *ddI2b,
29  double *dI3b)
30 {
31  constexpr int DIM = 3;
32  kernels::InvariantsEvaluator3D ie(Args()
33  .J(J).B(B)
34  .dI1b(dI1b).ddI1b(ddI1b)
35  .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
36  .dI3b(dI3b));
37  const double c1 = weight/9.;
38  const double I1b = ie.Get_I1b();
39  const double I2b = ie.Get_I2b();
40  ConstDeviceMatrix di1b(ie.Get_dI1b(),DIM,DIM);
41  ConstDeviceMatrix di2b(ie.Get_dI2b(),DIM,DIM);
42  for (int i = 0; i < DIM; i++)
43  {
44  for (int j = 0; j < DIM; j++)
45  {
46  ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
47  ConstDeviceMatrix ddi2b(ie.Get_ddI2b(i,j),DIM,DIM);
48  for (int r = 0; r < DIM; r++)
49  {
50  for (int c = 0; c < DIM; c++)
51  {
52  const double dp =
53  (di2b(r,c)*di1b(i,j) + di1b(r,c)*di2b(i,j))
54  + ddi2b(r,c)*I1b
55  + ddi1b(r,c)*I2b;
56  dP(r,c,i,j,qx,qy,qz,e) = c1 * dp;
57  }
58  }
59  }
60  }
61 }
62 
63 // dP_303 = ddI1b/3
64 static MFEM_HOST_DEVICE inline
65 void EvalH_303(const int e, const int qx, const int qy, const int qz,
66  const double weight, const double *J, DeviceTensor<8,double> dP,
67  double *B, double *dI1b, double *ddI1, double *ddI1b,
68  double *dI2, double *dI2b, double *ddI2, double *ddI2b,
69  double *dI3b, double *ddI3b)
70 {
71  constexpr int DIM = 3;
72  kernels::InvariantsEvaluator3D ie(Args()
73  .J(J).B(B)
74  .dI1b(dI1b).ddI1(ddI1).ddI1b(ddI1b)
75  .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
76  .dI3b(dI3b).ddI3b(ddI3b));
77  const double c1 = weight/3.;
78  for (int i = 0; i < DIM; i++)
79  {
80  for (int j = 0; j < DIM; j++)
81  {
82  ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
83  for (int r = 0; r < DIM; r++)
84  {
85  for (int c = 0; c < DIM; c++)
86  {
87  const double dp = ddi1b(r,c);
88  dP(r,c,i,j,qx,qy,qz,e) = c1 * dp;
89  }
90  }
91  }
92  }
93 }
94 
95 // dP_315 = 2*(dI3b x dI3b) + 2*(I3b - 1)*ddI3b
96 static MFEM_HOST_DEVICE inline
97 void EvalH_315(const int e, const int qx, const int qy, const int qz,
98  const double weight, const double *J, DeviceTensor<8,double> dP,
99  double *dI3b, double *ddI3b)
100 {
101  constexpr int DIM = 3;
102  kernels::InvariantsEvaluator3D ie(Args().
103  J(J).
104  dI3b(dI3b).ddI3b(ddI3b));
105 
106  double sign_detJ;
107  const double I3b = ie.Get_I3b(sign_detJ);
108  ConstDeviceMatrix di3b(ie.Get_dI3b(sign_detJ),DIM,DIM);
109 
110  for (int i = 0; i < DIM; i++)
111  {
112  for (int j = 0; j < DIM; j++)
113  {
114  ConstDeviceMatrix ddi3b(ie.Get_ddI3b(i,j),DIM,DIM);
115  for (int r = 0; r < DIM; r++)
116  {
117  for (int c = 0; c < DIM; c++)
118  {
119  const double dp = 2.0 * weight * (I3b - 1.0) * ddi3b(r,c) +
120  2.0 * weight * di3b(r,c) * di3b(i,j);
121  dP(r,c,i,j,qx,qy,qz,e) = dp;
122  }
123  }
124  }
125  }
126 }
127 
128 // dP_318 = (I3b - 1/I3b^3)*ddI3b + (1 + 3/I3b^4)*(dI3b x dI3b)
129 // Uses the I3b form, as dI3 and ddI3 were not implemented at the time.
130 static MFEM_HOST_DEVICE inline
131 void EvalH_318(const int e, const int qx, const int qy, const int qz,
132  const double weight, const double *J, DeviceTensor<8,double> dP,
133  double *dI3b, double *ddI3b)
134 {
135  constexpr int DIM = 3;
136  kernels::InvariantsEvaluator3D ie(Args().
137  J(J).
138  dI3b(dI3b).ddI3b(ddI3b));
139 
140  double sign_detJ;
141  const double I3b = ie.Get_I3b(sign_detJ);
142  ConstDeviceMatrix di3b(ie.Get_dI3b(sign_detJ),DIM,DIM);
143 
144  for (int i = 0; i < DIM; i++)
145  {
146  for (int j = 0; j < DIM; j++)
147  {
148  ConstDeviceMatrix ddi3b(ie.Get_ddI3b(i,j),DIM,DIM);
149  for (int r = 0; r < DIM; r++)
150  {
151  for (int c = 0; c < DIM; c++)
152  {
153  const double dp =
154  weight * (I3b - 1.0/(I3b*I3b*I3b)) * ddi3b(r,c) +
155  weight * (1.0 + 3.0/(I3b*I3b*I3b*I3b)) * di3b(r,c)*di3b(i,j);
156  dP(r,c,i,j,qx,qy,qz,e) = dp;
157  }
158  }
159  }
160  }
161 }
162 
163 // dP_321 = ddI1 + (-2/I3b^3)*(dI2 x dI3b + dI3b x dI2)
164 // + (1/I3)*ddI2
165 // + (6*I2/I3b^4)*(dI3b x dI3b)
166 // + (-2*I2/I3b^3)*ddI3b
167 static MFEM_HOST_DEVICE inline
168 void EvalH_321(const int e, const int qx, const int qy, const int qz,
169  const double weight, const double *J, DeviceTensor<8,double> dP,
170  double *B, double *dI1b, double *ddI1, double *ddI1b,
171  double *dI2, double *dI2b, double *ddI2, double *ddI2b,
172  double *dI3b, double *ddI3b)
173 {
174  constexpr int DIM = 3;
175  kernels::InvariantsEvaluator3D ie(Args()
176  .J(J).B(B)
177  .dI1b(dI1b).ddI1(ddI1).ddI1b(ddI1b)
178  .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
179  .dI3b(dI3b).ddI3b(ddI3b));
180  double sign_detJ;
181  const double I2 = ie.Get_I2();
182  const double I3b = ie.Get_I3b(sign_detJ);
183 
184  ConstDeviceMatrix di2(ie.Get_dI2(),DIM,DIM);
185  ConstDeviceMatrix di3b(ie.Get_dI3b(sign_detJ),DIM,DIM);
186 
187  const double c0 = 1.0/I3b;
188  const double c1 = weight*c0*c0;
189  const double c2 = -2*c0*c1;
190  const double c3 = c2*I2;
191 
192  for (int i = 0; i < DIM; i++)
193  {
194  for (int j = 0; j < DIM; j++)
195  {
196  ConstDeviceMatrix ddi1(ie.Get_ddI1(i,j),DIM,DIM);
197  ConstDeviceMatrix ddi2(ie.Get_ddI2(i,j),DIM,DIM);
198  ConstDeviceMatrix ddi3b(ie.Get_ddI3b(i,j),DIM,DIM);
199  for (int r = 0; r < DIM; r++)
200  {
201  for (int c = 0; c < DIM; c++)
202  {
203  const double dp =
204  weight * ddi1(r,c)
205  + c1 * ddi2(r,c)
206  + c3 * ddi3b(r,c)
207  + c2 * ((di2(r,c)*di3b(i,j) + di3b(r,c)*di2(i,j)))
208  -3*c0*c3 * di3b(r,c)*di3b(i,j);
209  dP(r,c,i,j,qx,qy,qz,e) = dp;
210  }
211  }
212  }
213  }
214 }
215 
216 // H_332 = w0 H_302 + w1 H_315
217 static MFEM_HOST_DEVICE inline
218 void EvalH_332(const int e, const int qx, const int qy, const int qz,
219  const double weight, const double *w,
220  const double *J, DeviceTensor<8,double> dP,
221  double *B, double *dI1b, double *ddI1b,
222  double *dI2, double *dI2b, double *ddI2, double *ddI2b,
223  double *dI3b, double *ddI3b)
224 {
225  constexpr int DIM = 3;
226  kernels::InvariantsEvaluator3D ie(Args()
227  .J(J).B(B)
228  .dI1b(dI1b).ddI1b(ddI1b)
229  .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
230  .dI3b(dI3b).ddI3b(ddI3b));
231  double sign_detJ;
232  const double c1 = weight/9.;
233  const double I1b = ie.Get_I1b();
234  const double I2b = ie.Get_I2b();
235  const double I3b = ie.Get_I3b(sign_detJ);
236  ConstDeviceMatrix di1b(ie.Get_dI1b(),DIM,DIM);
237  ConstDeviceMatrix di2b(ie.Get_dI2b(),DIM,DIM);
238  ConstDeviceMatrix di3b(ie.Get_dI3b(sign_detJ),DIM,DIM);
239  for (int i = 0; i < DIM; i++)
240  {
241  for (int j = 0; j < DIM; j++)
242  {
243  ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
244  ConstDeviceMatrix ddi2b(ie.Get_ddI2b(i,j),DIM,DIM);
245  ConstDeviceMatrix ddi3b(ie.Get_ddI3b(i,j),DIM,DIM);
246  for (int r = 0; r < DIM; r++)
247  {
248  for (int c = 0; c < DIM; c++)
249  {
250  const double dp_302 =
251  (di2b(r,c)*di1b(i,j) + di1b(r,c)*di2b(i,j))
252  + ddi2b(r,c)*I1b
253  + ddi1b(r,c)*I2b;
254  const double dp_315 = 2.0 * weight * (I3b - 1.0) * ddi3b(r,c) +
255  2.0 * weight * di3b(r,c) * di3b(i,j);
256  dP(r,c,i,j,qx,qy,qz,e) = w[0] * c1 * dp_302 +
257  w[1] * dp_315;
258  }
259  }
260  }
261  }
262 }
263 
264 // H_338 = w0 H_302 + w1 H_318
265 static MFEM_HOST_DEVICE inline
266 void EvalH_338(const int e, const int qx, const int qy, const int qz,
267  const double weight, const double *w,
268  const double *J, DeviceTensor<8,double> dP,
269  double *B, double *dI1b, double *ddI1b,
270  double *dI2, double *dI2b, double *ddI2, double *ddI2b,
271  double *dI3b, double *ddI3b)
272 {
273  constexpr int DIM = 3;
274  kernels::InvariantsEvaluator3D ie(Args()
275  .J(J).B(B)
276  .dI1b(dI1b).ddI1b(ddI1b)
277  .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
278  .dI3b(dI3b).ddI3b(ddI3b));
279  double sign_detJ;
280  const double c1 = weight/9.;
281  const double I1b = ie.Get_I1b();
282  const double I2b = ie.Get_I2b();
283  const double I3b = ie.Get_I3b(sign_detJ);
284  ConstDeviceMatrix di1b(ie.Get_dI1b(),DIM,DIM);
285  ConstDeviceMatrix di2b(ie.Get_dI2b(),DIM,DIM);
286  ConstDeviceMatrix di3b(ie.Get_dI3b(sign_detJ),DIM,DIM);
287  for (int i = 0; i < DIM; i++)
288  {
289  for (int j = 0; j < DIM; j++)
290  {
291  ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
292  ConstDeviceMatrix ddi2b(ie.Get_ddI2b(i,j),DIM,DIM);
293  ConstDeviceMatrix ddi3b(ie.Get_ddI3b(i,j),DIM,DIM);
294  for (int r = 0; r < DIM; r++)
295  {
296  for (int c = 0; c < DIM; c++)
297  {
298  const double dp_302 =
299  (di2b(r,c)*di1b(i,j) + di1b(r,c)*di2b(i,j))
300  + ddi2b(r,c)*I1b
301  + ddi1b(r,c)*I2b;
302  const double dp_318 =
303  weight * (I3b - 1.0/(I3b*I3b*I3b)) * ddi3b(r,c) +
304  weight * (1.0 + 3.0/(I3b*I3b*I3b*I3b)) * di3b(r,c)*di3b(i,j);
305  dP(r,c,i,j,qx,qy,qz,e) = w[0] * c1 * dp_302 +
306  w[1] * dp_318;
307  }
308  }
309  }
310  }
311 }
312 
313 MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_3D,
314  const double metric_normal,
315  const Array<double> &metric_param,
316  const int mid,
317  const Vector &x_,
318  const int NE,
319  const Array<double> &w_,
320  const Array<double> &b_,
321  const Array<double> &g_,
322  const DenseTensor &j_,
323  Vector &h_,
324  const int d1d,
325  const int q1d)
326 {
327  MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 || mid == 318 ||
328  mid == 321 || mid == 332 || mid == 338,
329  "3D metric not yet implemented!");
330 
331  constexpr int DIM = 3;
332  const int D1D = T_D1D ? T_D1D : d1d;
333  const int Q1D = T_Q1D ? T_Q1D : q1d;
334 
335  const auto b = Reshape(b_.Read(), Q1D, D1D);
336  const auto g = Reshape(g_.Read(), Q1D, D1D);
337  const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
338  const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
339  const auto X = Reshape(x_.Read(), D1D, D1D, D1D, DIM, NE);
340  auto H = Reshape(h_.Write(), DIM, DIM, DIM, DIM, Q1D, Q1D, Q1D, NE);
341 
342  const double *metric_data = metric_param.Read();
343 
344  mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
345  {
346  const int D1D = T_D1D ? T_D1D : d1d;
347  const int Q1D = T_Q1D ? T_Q1D : q1d;
348  constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
349  constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
350 
351  MFEM_SHARED double s_BG[2][MQ1*MD1];
352  MFEM_SHARED double s_DDD[3][MD1*MD1*MD1];
353  MFEM_SHARED double s_DDQ[9][MD1*MD1*MQ1];
354  MFEM_SHARED double s_DQQ[9][MD1*MQ1*MQ1];
355  MFEM_SHARED double s_QQQ[9][MQ1*MQ1*MQ1];
356 
357  kernels::internal::LoadX<MD1>(e,D1D,X,s_DDD);
358  kernels::internal::LoadBG<MD1,MQ1>(D1D,Q1D,b,g,s_BG);
359 
360  kernels::internal::GradX<MD1,MQ1>(D1D,Q1D,s_BG,s_DDD,s_DDQ);
361  kernels::internal::GradY<MD1,MQ1>(D1D,Q1D,s_BG,s_DDQ,s_DQQ);
362  kernels::internal::GradZ<MD1,MQ1>(D1D,Q1D,s_BG,s_DQQ,s_QQQ);
363 
364  MFEM_FOREACH_THREAD(qz,z,Q1D)
365  {
366  MFEM_FOREACH_THREAD(qy,y,Q1D)
367  {
368  MFEM_FOREACH_THREAD(qx,x,Q1D)
369  {
370  const double *Jtr = &J(0,0,qx,qy,qz,e);
371  const double detJtr = kernels::Det<3>(Jtr);
372  const double weight = metric_normal * W(qx,qy,qz) * detJtr;
373 
374  // Jrt = Jtr^{-1}
375  double Jrt[9];
376  kernels::CalcInverse<3>(Jtr, Jrt);
377 
378  // Jpr = X^T.DSh
379  double Jpr[9];
380  kernels::internal::PullGrad<MQ1>(Q1D,qx,qy,qz, s_QQQ, Jpr);
381 
382  // Jpt = X^T . DS = (X^T.DSh) . Jrt = Jpr . Jrt
383  double Jpt[9];
384  kernels::Mult(3,3,3, Jpr, Jrt, Jpt);
385 
386  // InvariantsEvaluator3D buffers used for the metrics
387  double B[9];
388  double dI1b[9], ddI1[9], ddI1b[9];
389  double dI2[9], dI2b[9], ddI2[9], ddI2b[9];
390  // reuse local arrays, to help register allocation
391  double *dI3b=Jrt, *ddI3b=Jpr;
392 
393  // metric->AssembleH
394  if (mid == 302)
395  {
396  EvalH_302(e,qx,qy,qz,weight,Jpt,H,
397  B,dI1b,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b);
398  }
399  if (mid == 303)
400  {
401  EvalH_303(e,qx,qy,qz,weight,Jpt,H,
402  B,dI1b,ddI1,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
403  }
404  if (mid == 315)
405  {
406  EvalH_315(e,qx,qy,qz,weight,Jpt,H, dI3b,ddI3b);
407  }
408  if (mid == 318)
409  {
410  EvalH_318(e,qx,qy,qz,weight,Jpt,H, dI3b,ddI3b);
411  }
412  if (mid == 321)
413  {
414  EvalH_321(e,qx,qy,qz,weight,Jpt,H,
415  B,dI1b,ddI1,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
416  }
417  if (mid == 332)
418  {
419  EvalH_332(e,qx,qy,qz,weight,metric_data,Jpt,H,
420  B,dI1b,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
421  }
422  if (mid == 338)
423  {
424  EvalH_338(e,qx,qy,qz,weight,metric_data,Jpt,H,
425  B,dI1b,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
426  }
427  } // qx
428  } // qy
429  } // qz
430  });
431 }
432 
433 void TMOP_Integrator::AssembleGradPA_3D(const Vector &X) const
434 {
435  const int N = PA.ne;
436  const int D1D = PA.maps->ndof;
437  const int Q1D = PA.maps->nqpt;
438  const int M = metric->Id();
439  const int id = (D1D << 4 ) | Q1D;
440  const double mn = metric_normal;
441  const DenseTensor &J = PA.Jtr;
442  const Array<double> &W = PA.ir->GetWeights();
443  const Array<double> &B = PA.maps->B;
444  const Array<double> &G = PA.maps->G;
445  Vector &H = PA.H;
446 
447  Array<double> mp;
448  if (auto m = dynamic_cast<TMOP_Combo_QualityMetric *>(metric))
449  {
450  m->GetWeights(mp);
451  }
452 
453  MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_3D,id,mn,mp,M,X,N,W,B,G,J,H);
454 }
455 
456 } // namespace mfem
mfem::Array::Read
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
mfem::forall_3D
void forall_3D(int N, int X, int Y, int Z, lambda &&body)
Definition: forall.hpp:763
mfem::ConstDeviceMatrix
DeviceTensor< 2, const double > ConstDeviceMatrix
Definition: dtensor.hpp:144
mfem::TMOP_Integrator::PA
struct mfem::TMOP_Integrator::@23 PA
mfem::TMOP_Integrator::H
Vector H
Definition: tmop.hpp:1849
mfem::TMOP_Integrator::metric
TMOP_QualityMetric * metric
Definition: tmop.hpp:1740
mfem::TMOP_Integrator::metric_normal
double metric_normal
Definition: tmop.hpp:1751
mfem::Vector::Read
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
DIM
constexpr int DIM
Definition: minimal-surface.cpp:71
tmop_pa.hpp
mfem::DenseTensor::Read
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:1230
mfem::MFEM_REGISTER_TMOP_KERNELS
MFEM_REGISTER_TMOP_KERNELS(void, DatcSize, const int NE, const int ncomp, const int sizeidx, const double input_min_size, const DenseMatrix &w_, const Array< double > &b_, const Vector &x_, const Vector &nc_reduce, DenseTensor &j_, const int d1d, const int q1d)
Definition: tmop_pa_da3.cpp:20
mfem::TMOP_QualityMetric::Id
virtual int Id() const
Return the metric ID.
Definition: tmop.hpp:78
mfem::Vector::Write
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
b
double b
Definition: lissajous.cpp:42
mfem::kernels::Mult
MFEM_HOST_DEVICE void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
Matrix vector multiplication: y = A x, where the matrix A is of size height x width with given data...
Definition: kernels.hpp:163
mfem::Args
kernels::InvariantsEvaluator2D::Buffers Args
Definition: tmop_pa_h2s.cpp:21
mfem::TMOP_Integrator::AssembleGradPA_3D
void AssembleGradPA_3D(const Vector &) const
Definition: tmop_pa_h3s.cpp:433
mfem::Vector
Vector data type.
Definition: vector.hpp:58
mfem::Reshape
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
mfem::DenseTensor
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:1096
Generated on Wed Sep 27 2023 14:52:31 for MFEM by   doxygen 1.8.14