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::@26 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)