21using Args = kernels::InvariantsEvaluator3D::Buffers;
24static MFEM_HOST_DEVICE
inline
25void EvalH_302(
const int e,
const int qx,
const int qy,
const int qz,
26 const real_t weight,
const real_t *J, DeviceTensor<8,real_t> dP,
31 constexpr int DIM = 3;
32 kernels::InvariantsEvaluator3D ie(
Args()
34 .dI1b(dI1b).ddI1b(ddI1b)
35 .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
37 const real_t c1 = weight/9.;
38 const real_t I1b = ie.Get_I1b();
39 const real_t I2b = ie.Get_I2b();
42 for (
int i = 0; i <
DIM; i++)
44 for (
int j = 0; j <
DIM; j++)
48 for (
int r = 0; r <
DIM; r++)
50 for (
int c = 0; c <
DIM; c++)
53 (di2b(r,c)*di1b(i,j) + di1b(r,c)*di2b(i,j))
56 dP(r,c,i,j,qx,qy,qz,e) = c1 * dp;
64static MFEM_HOST_DEVICE
inline
65void EvalH_303(
const int e,
const int qx,
const int qy,
const int qz,
66 const real_t weight,
const real_t *J, DeviceTensor<8,real_t> dP,
71 constexpr int DIM = 3;
72 kernels::InvariantsEvaluator3D ie(
Args()
74 .dI1b(dI1b).ddI1(ddI1).ddI1b(ddI1b)
75 .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
76 .dI3b(dI3b).ddI3b(ddI3b));
77 const real_t c1 = weight/3.;
78 for (
int i = 0; i <
DIM; i++)
80 for (
int j = 0; j <
DIM; j++)
83 for (
int r = 0; r <
DIM; r++)
85 for (
int c = 0; c <
DIM; c++)
87 const real_t dp = ddi1b(r,c);
88 dP(r,c,i,j,qx,qy,qz,e) = c1 * dp;
96static MFEM_HOST_DEVICE
inline
97void EvalH_315(
const int e,
const int qx,
const int qy,
const int qz,
98 const real_t weight,
const real_t *J, DeviceTensor<8,real_t> dP,
101 constexpr int DIM = 3;
102 kernels::InvariantsEvaluator3D ie(
Args().
104 dI3b(dI3b).ddI3b(ddI3b));
107 const real_t I3b = ie.Get_I3b(sign_detJ);
110 for (
int i = 0; i <
DIM; i++)
112 for (
int j = 0; j <
DIM; j++)
115 for (
int r = 0; r <
DIM; r++)
117 for (
int c = 0; c <
DIM; c++)
119 const real_t 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;
130static MFEM_HOST_DEVICE
inline
131void EvalH_318(
const int e,
const int qx,
const int qy,
const int qz,
132 const real_t weight,
const real_t *J, DeviceTensor<8,real_t> dP,
135 constexpr int DIM = 3;
136 kernels::InvariantsEvaluator3D ie(
Args().
138 dI3b(dI3b).ddI3b(ddI3b));
141 const real_t I3b = ie.Get_I3b(sign_detJ);
144 for (
int i = 0; i <
DIM; i++)
146 for (
int j = 0; j <
DIM; j++)
149 for (
int r = 0; r <
DIM; r++)
151 for (
int c = 0; c <
DIM; c++)
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;
167static MFEM_HOST_DEVICE
inline
168void EvalH_321(
const int e,
const int qx,
const int qy,
const int qz,
169 const real_t weight,
const real_t *J, DeviceTensor<8,real_t> dP,
174 constexpr int DIM = 3;
175 kernels::InvariantsEvaluator3D ie(
Args()
177 .dI1b(dI1b).ddI1(ddI1).ddI1b(ddI1b)
178 .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
179 .dI3b(dI3b).ddI3b(ddI3b));
181 const real_t I2 = ie.Get_I2();
182 const real_t I3b = ie.Get_I3b(sign_detJ);
187 const real_t c0 = 1.0/I3b;
188 const real_t c1 = weight*c0*c0;
189 const real_t c2 = -2*c0*c1;
192 for (
int i = 0; i <
DIM; i++)
194 for (
int j = 0; j <
DIM; j++)
199 for (
int r = 0; r <
DIM; r++)
201 for (
int c = 0; c <
DIM; 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;
217static MFEM_HOST_DEVICE
inline
218void EvalH_332(
const int e,
const int qx,
const int qy,
const int qz,
220 const real_t *J, DeviceTensor<8,real_t> dP,
225 constexpr int DIM = 3;
226 kernels::InvariantsEvaluator3D ie(
Args()
228 .dI1b(dI1b).ddI1b(ddI1b)
229 .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
230 .dI3b(dI3b).ddI3b(ddI3b));
232 const real_t c1 = weight/9.;
233 const real_t I1b = ie.Get_I1b();
234 const real_t I2b = ie.Get_I2b();
235 const real_t I3b = ie.Get_I3b(sign_detJ);
239 for (
int i = 0; i <
DIM; i++)
241 for (
int j = 0; j <
DIM; j++)
246 for (
int r = 0; r <
DIM; r++)
248 for (
int c = 0; c <
DIM; c++)
251 (di2b(r,c)*di1b(i,j) + di1b(r,c)*di2b(i,j))
254 const real_t 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 +
265static MFEM_HOST_DEVICE
inline
266void EvalH_338(
const int e,
const int qx,
const int qy,
const int qz,
268 const real_t *J, DeviceTensor<8,real_t> dP,
273 constexpr int DIM = 3;
274 kernels::InvariantsEvaluator3D ie(
Args()
276 .dI1b(dI1b).ddI1b(ddI1b)
277 .dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
278 .dI3b(dI3b).ddI3b(ddI3b));
280 const real_t c1 = weight/9.;
281 const real_t I1b = ie.Get_I1b();
282 const real_t I2b = ie.Get_I2b();
283 const real_t I3b = ie.Get_I3b(sign_detJ);
287 for (
int i = 0; i <
DIM; i++)
289 for (
int j = 0; j <
DIM; j++)
294 for (
int r = 0; r <
DIM; r++)
296 for (
int c = 0; c <
DIM; c++)
299 (di2b(r,c)*di1b(i,j) + di1b(r,c)*di2b(i,j))
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 +
314 const real_t metric_normal,
328 MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 || mid == 318 ||
329 mid == 321 || mid == 332 || mid == 338,
330 "3D metric not yet implemented!");
332 const bool const_m0 = mc_.
Size() == 1;
334 constexpr int DIM = 3;
335 const int D1D = T_D1D ? T_D1D : d1d;
336 const int Q1D = T_Q1D ? T_Q1D : q1d;
338 const auto MC = const_m0 ?
348 const real_t *metric_data = metric_param.
Read();
352 const int D1D = T_D1D ? T_D1D : d1d;
353 const int Q1D = T_Q1D ? T_Q1D : q1d;
354 constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
355 constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
357 MFEM_SHARED
real_t s_BG[2][MQ1*MD1];
358 MFEM_SHARED
real_t s_DDD[3][MD1*MD1*MD1];
359 MFEM_SHARED
real_t s_DDQ[9][MD1*MD1*MQ1];
360 MFEM_SHARED
real_t s_DQQ[9][MD1*MQ1*MQ1];
361 MFEM_SHARED
real_t s_QQQ[9][MQ1*MQ1*MQ1];
363 kernels::internal::LoadX<MD1>(e,D1D,X,s_DDD);
364 kernels::internal::LoadBG<MD1,MQ1>(D1D,Q1D,
b,g,s_BG);
366 kernels::internal::GradX<MD1,MQ1>(D1D,Q1D,s_BG,s_DDD,s_DDQ);
367 kernels::internal::GradY<MD1,MQ1>(D1D,Q1D,s_BG,s_DDQ,s_DQQ);
368 kernels::internal::GradZ<MD1,MQ1>(D1D,Q1D,s_BG,s_DQQ,s_QQQ);
370 MFEM_FOREACH_THREAD(qz,z,Q1D)
372 MFEM_FOREACH_THREAD(qy,y,Q1D)
374 MFEM_FOREACH_THREAD(qx,x,Q1D)
376 const real_t *Jtr = &J(0,0,qx,qy,qz,e);
378 const real_t m_coef = const_m0 ? MC(0,0,0,0) : MC(qx,qy,qz,e);
379 const real_t weight = metric_normal * m_coef *
380 W(qx,qy,qz) * detJtr;
388 kernels::internal::PullGrad<MQ1>(Q1D,qx,qy,qz, s_QQQ, Jpr);
396 real_t dI1b[9], ddI1[9], ddI1b[9];
397 real_t dI2[9], dI2b[9], ddI2[9], ddI2b[9];
399 real_t *dI3b=Jrt, *ddI3b=Jpr;
404 EvalH_302(e,qx,qy,qz,weight,Jpt,H,
405 B,dI1b,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b);
409 EvalH_303(e,qx,qy,qz,weight,Jpt,H,
410 B,dI1b,ddI1,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
414 EvalH_315(e,qx,qy,qz,weight,Jpt,H, dI3b,ddI3b);
418 EvalH_318(e,qx,qy,qz,weight,Jpt,H, dI3b,ddI3b);
422 EvalH_321(e,qx,qy,qz,weight,Jpt,H,
423 B,dI1b,ddI1,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
427 EvalH_332(e,qx,qy,qz,weight,metric_data,Jpt,H,
428 B,dI1b,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
432 EvalH_338(e,qx,qy,qz,weight,metric_data,Jpt,H,
433 B,dI1b,ddI1b,dI2,dI2b,ddI2,ddI2b,dI3b,ddI3b);
444 const int D1D =
PA.maps->ndof;
445 const int Q1D =
PA.maps->nqpt;
447 const int id = (D1D << 4 ) | Q1D;
462 MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_3D,
id,mn,
MC,mp,M,X,N,W,B,G,J,
H);
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Rank 3 tensor (array of matrices)
const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
TMOP_QualityMetric * metric
struct mfem::TMOP_Integrator::@23 PA
void AssembleGradPA_3D(const Vector &) const
virtual int Id() const
Return the metric ID.
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
int Size() const
Returns the size of the vector.
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
MFEM_HOST_DEVICE void CalcInverse(const T *data, T *inv_data)
Return the inverse of a matrix with given size and data into the matrix with data inv_data.
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,...
MFEM_HOST_DEVICE T Det(const T *data)
Compute the determinant of a square matrix of size dim with given data.
MFEM_REGISTER_TMOP_KERNELS(void, DatcSize, const int NE, const int ncomp, const int sizeidx, const real_t input_min_size, const DenseMatrix &w_, const Array< real_t > &b_, const Vector &x_, const Vector &nc_reduce, DenseTensor &j_, const int d1d, const int q1d)
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
DeviceTensor< 2, const real_t > ConstDeviceMatrix
kernels::InvariantsEvaluator2D::Buffers Args
void forall_3D(int N, int X, int Y, int Z, lambda &&body)