MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
tmop_pa_h3s.cpp
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 #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_321 = ddI1 + (-2/I3b^3)*(dI2 x dI3b + dI3b x dI2)
129 // + (1/I3)*ddI2
130 // + (6*I2/I3b^4)*(dI3b x dI3b)
131 // + (-2*I2/I3b^3)*ddI3b
132 static MFEM_HOST_DEVICE inline
133 void EvalH_321(const int e, const int qx, const int qy, const int qz,
134  const double weight, const double *J, DeviceTensor<8,double> dP,
135  double *B, double *dI1b, double *ddI1, double *ddI1b,
136  double *dI2, double *dI2b, double *ddI2, double *ddI2b,
137  double *dI3b, double *ddI3b)
138 {
139  constexpr int DIM = 3;
140  kernels::InvariantsEvaluator3D ie(Args()
141  .J(J).B(B)
142  .dI1b(dI1b).ddI1(ddI1).ddI1b(ddI1b)
143  .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
144  .dI3b(dI3b).ddI3b(ddI3b));
145  double sign_detJ;
146  const double I2 = ie.Get_I2();
147  const double I3b = ie.Get_I3b(sign_detJ);
148 
149  ConstDeviceMatrix di2(ie.Get_dI2(),DIM,DIM);
150  ConstDeviceMatrix di3b(ie.Get_dI3b(sign_detJ),DIM,DIM);
151 
152  const double c0 = 1.0/I3b;
153  const double c1 = weight*c0*c0;
154  const double c2 = -2*c0*c1;
155  const double c3 = c2*I2;
156 
157  for (int i = 0; i < DIM; i++)
158  {
159  for (int j = 0; j < DIM; j++)
160  {
161  ConstDeviceMatrix ddi1(ie.Get_ddI1(i,j),DIM,DIM);
162  ConstDeviceMatrix ddi2(ie.Get_ddI2(i,j),DIM,DIM);
163  ConstDeviceMatrix ddi3b(ie.Get_ddI3b(i,j),DIM,DIM);
164  for (int r = 0; r < DIM; r++)
165  {
166  for (int c = 0; c < DIM; c++)
167  {
168  const double dp =
169  weight * ddi1(r,c)
170  + c1 * ddi2(r,c)
171  + c3 * ddi3b(r,c)
172  + c2 * ((di2(r,c)*di3b(i,j) + di3b(r,c)*di2(i,j)))
173  -3*c0*c3 * di3b(r,c)*di3b(i,j);
174  dP(r,c,i,j,qx,qy,qz,e) = dp;
175  }
176  }
177  }
178  }
179 }
180 
181 // H_332 = (1-gamma) H_302 + gamma H_315
182 static MFEM_HOST_DEVICE inline
183 void EvalH_332(const int e, const int qx, const int qy, const int qz,
184  const double weight, const double gamma,
185  const double *J, DeviceTensor<8,double> dP,
186  double *B, double *dI1b, double *ddI1b,
187  double *dI2, double *dI2b, double *ddI2, double *ddI2b,
188  double *dI3b, double *ddI3b)
189 {
190  constexpr int DIM = 3;
191  kernels::InvariantsEvaluator3D ie(Args()
192  .J(J).B(B)
193  .dI1b(dI1b).ddI1b(ddI1b)
194  .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
195  .dI3b(dI3b).ddI3b(ddI3b));
196  double sign_detJ;
197  const double c1 = weight/9.;
198  const double I1b = ie.Get_I1b();
199  const double I2b = ie.Get_I2b();
200  const double I3b = ie.Get_I3b(sign_detJ);
201  ConstDeviceMatrix di1b(ie.Get_dI1b(),DIM,DIM);
202  ConstDeviceMatrix di2b(ie.Get_dI2b(),DIM,DIM);
203  ConstDeviceMatrix di3b(ie.Get_dI3b(sign_detJ),DIM,DIM);
204  for (int i = 0; i < DIM; i++)
205  {
206  for (int j = 0; j < DIM; j++)
207  {
208  ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
209  ConstDeviceMatrix ddi2b(ie.Get_ddI2b(i,j),DIM,DIM);
210  ConstDeviceMatrix ddi3b(ie.Get_ddI3b(i,j),DIM,DIM);
211  for (int r = 0; r < DIM; r++)
212  {
213  for (int c = 0; c < DIM; c++)
214  {
215  const double dp_302 =
216  (di2b(r,c)*di1b(i,j) + di1b(r,c)*di2b(i,j))
217  + ddi2b(r,c)*I1b
218  + ddi1b(r,c)*I2b;
219  const double dp_315 = 2.0 * weight * (I3b - 1.0) * ddi3b(r,c) +
220  2.0 * weight * di3b(r,c) * di3b(i,j);
221  dP(r,c,i,j,qx,qy,qz,e) = (1.0 - gamma) * c1 * dp_302 +
222  gamma * dp_315;
223  }
224  }
225  }
226  }
227 }
228 
229 MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_3D,
230  const double metric_normal,
231  const double metric_param,
232  const int mid,
233  const Vector &x_,
234  const int NE,
235  const Array<double> &w_,
236  const Array<double> &b_,
237  const Array<double> &g_,
238  const DenseTensor &j_,
239  Vector &h_,
240  const int d1d,
241  const int q1d)
242 {
243  MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 ||
244  mid == 321 || mid == 332, "3D metric not yet implemented!");
245 
246  constexpr int DIM = 3;
247  const int D1D = T_D1D ? T_D1D : d1d;
248  const int Q1D = T_Q1D ? T_Q1D : q1d;
249 
250  const auto b = Reshape(b_.Read(), Q1D, D1D);
251  const auto g = Reshape(g_.Read(), Q1D, D1D);
252  const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
253  const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
254  const auto X = Reshape(x_.Read(), D1D, D1D, D1D, DIM, NE);
255  auto H = Reshape(h_.Write(), DIM, DIM, DIM, DIM, Q1D, Q1D, Q1D, NE);
256 
257  MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
258  {
259  const int D1D = T_D1D ? T_D1D : d1d;
260  const int Q1D = T_Q1D ? T_Q1D : q1d;
261  constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
262  constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
263 
264  MFEM_SHARED double s_BG[2][MQ1*MD1];
265  MFEM_SHARED double s_DDD[3][MD1*MD1*MD1];
266  MFEM_SHARED double s_DDQ[9][MD1*MD1*MQ1];
267  MFEM_SHARED double s_DQQ[9][MD1*MQ1*MQ1];
268  MFEM_SHARED double s_QQQ[9][MQ1*MQ1*MQ1];
269 
270  kernels::internal::LoadX<MD1>(e,D1D,X,s_DDD);
271  kernels::internal::LoadBG<MD1,MQ1>(D1D,Q1D,b,g,s_BG);
272 
273  kernels::internal::GradX<MD1,MQ1>(D1D,Q1D,s_BG,s_DDD,s_DDQ);
274  kernels::internal::GradY<MD1,MQ1>(D1D,Q1D,s_BG,s_DDQ,s_DQQ);
275  kernels::internal::GradZ<MD1,MQ1>(D1D,Q1D,s_BG,s_DQQ,s_QQQ);
276 
277  MFEM_FOREACH_THREAD(qz,z,Q1D)
278  {
279  MFEM_FOREACH_THREAD(qy,y,Q1D)
280  {
281  MFEM_FOREACH_THREAD(qx,x,Q1D)
282  {
283  const double *Jtr = &J(0,0,qx,qy,qz,e);
284  const double detJtr = kernels::Det<3>(Jtr);
285  const double weight = metric_normal * W(qx,qy,qz) * detJtr;
286 
287  // Jrt = Jtr^{-1}
288  double Jrt[9];
289  kernels::CalcInverse<3>(Jtr, Jrt);
290 
291  // Jpr = X^T.DSh
292  double Jpr[9];
293  kernels::internal::PullGrad<MQ1>(Q1D,qx,qy,qz, s_QQQ, Jpr);
294 
295  // Jpt = X^T . DS = (X^T.DSh) . Jrt = Jpr . Jrt
296  double Jpt[9];
297  kernels::Mult(3,3,3, Jpr, Jrt, Jpt);
298 
299  // InvariantsEvaluator3D buffers used for the metrics
300  double B[9];
301  double dI1b[9], ddI1[9], ddI1b[9];
302  double dI2[9], dI2b[9], ddI2[9], ddI2b[9];
303  // reuse local arrays, to help register allocation
304  double *dI3b=Jrt, *ddI3b=Jpr;
305 
306  // metric->AssembleH
307  if (mid == 302)
308  {
309  EvalH_302(e,qx,qy,qz,weight,Jpt,H,
310  B,dI1b,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b);
311  }
312  if (mid == 303)
313  {
314  EvalH_303(e,qx,qy,qz,weight,Jpt,H,
315  B,dI1b,ddI1,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
316  }
317  if (mid == 315)
318  {
319  EvalH_315(e,qx,qy,qz,weight,Jpt,H, dI3b,ddI3b);
320  }
321  if (mid == 321)
322  {
323  EvalH_321(e,qx,qy,qz,weight,Jpt,H,
324  B,dI1b,ddI1,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
325  }
326  if (mid == 332)
327  {
328  EvalH_332(e,qx,qy,qz,weight,metric_param,Jpt,H,
329  B,dI1b,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
330  }
331  } // qx
332  } // qy
333  } // qz
334  });
335 }
336 
337 void TMOP_Integrator::AssembleGradPA_3D(const Vector &X) const
338 {
339  const int N = PA.ne;
340  const int D1D = PA.maps->ndof;
341  const int Q1D = PA.maps->nqpt;
342  const int M = metric->Id();
343  const int id = (D1D << 4 ) | Q1D;
344  const double mn = metric_normal;
345  const DenseTensor &J = PA.Jtr;
346  const Array<double> &W = PA.ir->GetWeights();
347  const Array<double> &B = PA.maps->B;
348  const Array<double> &G = PA.maps->G;
349  Vector &H = PA.H;
350 
351  double mp = 0.0;
352  if (auto m = dynamic_cast<TMOP_Metric_332 *>(metric)) { mp = m->GetGamma(); }
353 
354  MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_3D,id,mn,mp,M,X,N,W,B,G,J,H);
355 }
356 
357 } // namespace mfem
mfem::TMOP_QualityMetric::Id
virtual int Id() const
Return the metric ID.
Definition: tmop.hpp:78
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:1671
mfem::TMOP_Integrator::metric
TMOP_QualityMetric * metric
Definition: tmop.hpp:1572
mfem::DenseTensor::Read
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:1087
mfem::TMOP_Integrator::metric_normal
double metric_normal
Definition: tmop.hpp:1582
DIM
constexpr int DIM
Definition: minimal-surface.cpp:71
tmop_pa.hpp
mfem::Vector::Write
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:457
b
double b
Definition: lissajous.cpp:42
mfem::MFEM_REGISTER_TMOP_KERNELS
MFEM_REGISTER_TMOP_KERNELS(void, DatcSize, const int NE, const int ncomp, const int sizeidx, const DenseMatrix &w_, const Array< double > &b_, const Vector &x_, DenseTensor &j_, const int d1d, const int q1d)
Definition: tmop_pa_da3.cpp:20
mfem::Array::Read
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:304
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::TMOP_Integrator::AssembleGradPA_3D
void AssembleGradPA_3D(const Vector &) const
Definition: tmop_pa_h3s.cpp:337
mfem::Args
kernels::InvariantsEvaluator2D::Buffers Args
Definition: tmop_pa_h2s.cpp:21
mfem::Vector
Vector data type.
Definition: vector.hpp:60
mfem::DenseTensor
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:953
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:449
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
Generated on Wed Dec 7 2022 10:41:46 for MFEM by   doxygen 1.8.5