Merge pull request mfem#3827 from mfem/quadinterp-surface-phys-deriv

Support surface meshes in QuadratureInterpolator physical derivatives

data/diag-segment-2d.mesh

@@ -0,0 +1,37 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT = 0
+# SEGMENT = 1
+# TRIANGLE = 2
+# SQUARE = 3
+# TETRAHEDRON = 4
+# CUBE = 5
+# PRISM = 6
+#
+
+dimension
+1
+
+elements
+4
+1 1 0 1
+1 1 1 2
+1 1 2 3
+1 1 3 4
+
+boundary
+2
+1 0 0
+2 0 4
+
+vertices
+5
+2
+0 0
+0.25 0.25
+0.50 0.50
+0.75 0.75
+1 1

data/diag-segment-3d.mesh

@@ -0,0 +1,37 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT = 0
+# SEGMENT = 1
+# TRIANGLE = 2
+# SQUARE = 3
+# TETRAHEDRON = 4
+# CUBE = 5
+# PRISM = 6
+#
+
+dimension
+1
+
+elements
+4
+1 1 0 1
+1 1 1 2
+1 1 2 3
+1 1 3 4
+
+boundary
+2
+1 0 0
+2 0 4
+
+vertices
+5
+3
+0 0 0
+0.25 0.25 0.25
+0.50 0.50 0.50
+0.75 0.75 0.75
+1 1 1

fem/qinterp/grad.hpp

@@ -33,31 +33,60 @@ static void Derivatives1D(const int NE,
  const double *j_,
  const double *x_,
  double *y_,
+ const int sdim,
  const int vdim,
  const int d1d,
  const int q1d)
 {
  const auto g = Reshape(g_, q1d, d1d);
- const auto j = Reshape(j_, q1d, NE);
+ const auto j = Reshape(j_, q1d, sdim, NE);
  const auto x = Reshape(x_, d1d, vdim, NE);
  auto y = Q_LAYOUT == QVectorLayout::byNODES ?
- Reshape(y_, q1d, vdim, NE):
- Reshape(y_, vdim, q1d, NE);
+ Reshape(y_, q1d, vdim, sdim, NE):
+ Reshape(y_, vdim, sdim, q1d, NE);
 
  mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
  {
  for (int c = 0; c < vdim; c++)
  {
  for (int q = 0; q < q1d; q++)
  {
- double u = 0.0;
+ double du[3] = {0.0, 0.0, 0.0};
  for (int d = 0; d < d1d; d++)
  {
- u += g(q, d) * x(d, c, e);
+ du[0] += g(q, d) * x(d, c, e);
+ }
+ if (GRAD_PHYS)
+ {
+ if (sdim == 1) { du[0] /= j(q, 0, e); }
+ else if (sdim == 2)
+ {
+ const double Jloc[2] = {j(q,0,e), j(q,1,e)};
+ double Jinv[3];
+ kernels::CalcLeftInverse<2,1>(Jloc, Jinv);
+ const double U = Jinv[0]*du[0];
+ const double V = Jinv[1]*du[0];
+ du[0] = U;
+ du[1] = V;
+ }
+ else // sdim == 3
+ {
+ const double Jloc[3] = {j(q,0,e), j(q,1,e), j(q,2,e)};
+ double Jinv[3];
+ kernels::CalcLeftInverse<3,1>(Jloc, Jinv);
+ const double U = Jinv[0]*du[0];
+ const double V = Jinv[1]*du[0];
+ const double W = Jinv[2]*du[0];
+ du[0] = U;
+ du[1] = V;
+ du[2] = W;
+ }
+ }
+ for (int d = 0; d < sdim; ++d)
+ {
+ if (Q_LAYOUT == QVectorLayout::byVDIM) { y(c, d, q, e) = du[d]; }
+ if (Q_LAYOUT == QVectorLayout::byNODES) { y(q, c, d, e) = du[d]; }
  }
- if (GRAD_PHYS) { u /= j(q, e); }
- if (Q_LAYOUT == QVectorLayout::byVDIM) { y(c, q, e) = u; }
- if (Q_LAYOUT == QVectorLayout::byNODES) { y(q, c, e) = u; }
  }
  }
  });
@@ -73,22 +102,24 @@ static void Derivatives2D(const int NE,
  const double *j_,
  const double *x_,
  double *y_,
+ const int sdim = 2,
  const int vdim = 0,
  const int d1d = 0,
  const int q1d = 0)
 {
  const int D1D = T_D1D ? T_D1D : d1d;
  const int Q1D = T_Q1D ? T_Q1D : q1d;
  const int VDIM = T_VDIM ? T_VDIM : vdim;
+ const int SDIM = GRAD_PHYS ? sdim : 2;
  static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
 
  const auto b = Reshape(b_, Q1D, D1D);
  const auto g = Reshape(g_, Q1D, D1D);
- const auto j = Reshape(j_, Q1D, Q1D, 2, 2, NE);
+ const auto j = Reshape(j_, Q1D, Q1D, SDIM, 2, NE);
  const auto x = Reshape(x_, D1D, D1D, VDIM, NE);
  auto y = Q_LAYOUT == QVectorLayout:: byNODES ?
- Reshape(y_, Q1D, Q1D, VDIM, 2, NE):
- Reshape(y_, VDIM, 2, Q1D, Q1D, NE);
+ Reshape(y_, Q1D, Q1D, VDIM, SDIM, NE):
+ Reshape(y_, VDIM, SDIM, Q1D, Q1D, NE);
 
  mfem::forall_2D_batch(NE, Q1D, Q1D, NBZ, [=] MFEM_HOST_DEVICE (int e)
  {
@@ -135,34 +166,55 @@ static void Derivatives2D(const int NE,
  {
  MFEM_FOREACH_THREAD(qx,x,Q1D)
  {
- double u = 0.0;
- double v = 0.0;
+ double du[3] = {0.0, 0.0, 0.0};
  for (int dy = 0; dy < D1D; ++dy)
  {
- u += DQ1(dy,qx) * B(dy,qy);
- v += DQ0(dy,qx) * G(dy,qy);
+ du[0] += DQ1(dy,qx) * B(dy,qy);
+ du[1] += DQ0(dy,qx) * G(dy,qy);
  }
  if (GRAD_PHYS)
  {
- double Jloc[4], Jinv[4];
- Jloc[0] = j(qx,qy,0,0,e);
- Jloc[1] = j(qx,qy,1,0,e);
- Jloc[2] = j(qx,qy,0,1,e);
- Jloc[3] = j(qx,qy,1,1,e);
- kernels::CalcInverse<2>(Jloc, Jinv);
- const double U = Jinv[0]*u + Jinv[1]*v;
- const double V = Jinv[2]*u + Jinv[3]*v;
- u = U; v = V;
- }
- if (Q_LAYOUT == QVectorLayout::byVDIM)
- {
- y(c,0,qx,qy,e) = u;
- y(c,1,qx,qy,e) = v;
+ if (SDIM == 2)
+ {
+ double Jloc[4], Jinv[4];
+ Jloc[0] = j(qx,qy,0,0,e);
+ Jloc[1] = j(qx,qy,1,0,e);
+ Jloc[2] = j(qx,qy,0,1,e);
+ Jloc[3] = j(qx,qy,1,1,e);
+ kernels::CalcInverse<2>(Jloc, Jinv);
+ const double U = Jinv[0]*du[0] + Jinv[1]*du[1];
+ const double V = Jinv[2]*du[0] + Jinv[3]*du[1];
+ du[0] = U;
+ du[1] = V;
+ }
+ else
+ {
+ double Jloc[6], Jinv[6];
+ Jloc[0] = j(qx,qy,0,0,e);
+ Jloc[1] = j(qx,qy,1,0,e);
+ Jloc[2] = j(qx,qy,2,0,e);
+ Jloc[3] = j(qx,qy,0,1,e);
+ Jloc[4] = j(qx,qy,1,1,e);
+ Jloc[5] = j(qx,qy,2,1,e);
+ kernels::CalcLeftInverse<3,2>(Jloc, Jinv);
+ const double U = Jinv[0]*du[0] + Jinv[1]*du[1];
+ const double V = Jinv[2]*du[0] + Jinv[3]*du[1];
+ const double W = Jinv[4]*du[0] + Jinv[5]*du[1];
+ du[0] = U;
+ du[1] = V;
+ du[2] = W;
+ }
  }
- if (Q_LAYOUT == QVectorLayout::byNODES)
+ for (int d = 0; d < SDIM; ++d)
  {
- y(qx,qy,c,0,e) = u;
- y(qx,qy,c,1,e) = v;
+ if (Q_LAYOUT == QVectorLayout::byVDIM)
+ {
+ y(c,d,qx,qy,e) = du[d];
+ }
+ else // Q_LAYOUT == QVectorLayout::byNODES
+ {
+ y(qx,qy,c,d,e) = du[d];
+ }
  }
  }
  }

fem/qinterp/grad_by_nodes.cpp

@@ -47,7 +47,7 @@ void TensorDerivatives<QVectorLayout::byNODES>(const int NE,
 
  if (dim == 1)
  {
- return Derivatives1D<L,P>(NE,G,J,X,Y,vdim,D1D,Q1D);
+ return Derivatives1D<L,P>(NE,G,J,X,Y,dim,vdim,D1D,Q1D);
  }
  if (dim == 2)
  {
@@ -83,7 +83,7 @@ void TensorDerivatives<QVectorLayout::byNODES>(const int NE,
  {
  MFEM_ABORT("");
  }
- Derivatives2D<L,P,0,0,0,0,MD,MQ>(NE,B,G,J,X,Y,vdim,D1D,Q1D);
+ Derivatives2D<L,P,0,0,0,0,MD,MQ>(NE,B,G,J,X,Y,dim,vdim,D1D,Q1D);
  return;
  }
  }

fem/qinterp/grad_by_vdim.cpp

@@ -47,7 +47,7 @@ void TensorDerivatives<QVectorLayout::byVDIM>(const int NE,
 
  if (dim == 1)
  {
- return Derivatives1D<L,P>(NE,G,J,X,Y,vdim,D1D,Q1D);
+ return Derivatives1D<L,P>(NE,G,J,X,Y,dim,vdim,D1D,Q1D);
  }
  if (dim == 2)
  {
@@ -68,7 +68,7 @@ void TensorDerivatives<QVectorLayout::byVDIM>(const int NE,
  << " are not supported!");
  MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
  << MQ << " 1D points are not supported!");
- Derivatives2D<L,P,0,0,0,0,MD,MQ>(NE,B,G,J,X,Y,vdim,D1D,Q1D);
+ Derivatives2D<L,P,0,0,0,0,MD,MQ>(NE,B,G,J,X,Y,dim,vdim,D1D,Q1D);
  return;
  }
  }

fem/qinterp/grad_phys_by_nodes.cpp

@@ -37,7 +37,7 @@ void TensorPhysDerivatives<QVectorLayout::byNODES>(const int NE,
  const int D1D = maps.ndof;
  const int Q1D = maps.nqpt;
 
- MFEM_ASSERT(geom.mesh->SpaceDimension() == dim, "");
+ const int sdim = geom.mesh->SpaceDimension();
 
  const double *B = maps.B.Read();
  const double *G = maps.G.Read();
@@ -52,25 +52,25 @@ void TensorPhysDerivatives<QVectorLayout::byNODES>(const int NE,
 
  if (dim == 1)
  {
- return Derivatives1D<L,P>(NE,G,J,X,Y,vdim,D1D,Q1D);
+ return Derivatives1D<L,P>(NE,G,J,X,Y,sdim,vdim,D1D,Q1D);
  }
  if (dim == 2)
  {
  switch (id)
  {
- case 0x133: return Derivatives2D<L,P,1,3,3,8>(NE,B,G,J,X,Y);
- case 0x134: return Derivatives2D<L,P,1,3,4,8>(NE,B,G,J,X,Y);
- case 0x143: return Derivatives2D<L,P,1,4,3,4>(NE,B,G,J,X,Y);
- case 0x144: return Derivatives2D<L,P,1,4,4,4>(NE,B,G,J,X,Y);
- case 0x146: return Derivatives2D<L,P,1,4,6,4>(NE,B,G,J,X,Y);
- case 0x158: return Derivatives2D<L,P,1,5,8,2>(NE,B,G,J,X,Y);
+ case 0x133: return Derivatives2D<L,P,1,3,3,8>(NE,B,G,J,X,Y,sdim);
+ case 0x134: return Derivatives2D<L,P,1,3,4,8>(NE,B,G,J,X,Y,sdim);
+ case 0x143: return Derivatives2D<L,P,1,4,3,4>(NE,B,G,J,X,Y,sdim);
+ case 0x144: return Derivatives2D<L,P,1,4,4,4>(NE,B,G,J,X,Y,sdim);
+ case 0x146: return Derivatives2D<L,P,1,4,6,4>(NE,B,G,J,X,Y,sdim);
+ case 0x158: return Derivatives2D<L,P,1,5,8,2>(NE,B,G,J,X,Y,sdim);
 
- case 0x233: return Derivatives2D<L,P,2,3,3,8>(NE,B,G,J,X,Y);
- case 0x234: return Derivatives2D<L,P,2,3,4,8>(NE,B,G,J,X,Y);
- case 0x243: return Derivatives2D<L,P,2,4,3,4>(NE,B,G,J,X,Y);
- case 0x244: return Derivatives2D<L,P,2,4,4,4>(NE,B,G,J,X,Y);
- case 0x246: return Derivatives2D<L,P,2,4,6,4>(NE,B,G,J,X,Y);
- case 0x258: return Derivatives2D<L,P,2,5,8,2>(NE,B,G,J,X,Y);
+ case 0x233: return Derivatives2D<L,P,2,3,3,8>(NE,B,G,J,X,Y,sdim);
+ case 0x234: return Derivatives2D<L,P,2,3,4,8>(NE,B,G,J,X,Y,sdim);
+ case 0x243: return Derivatives2D<L,P,2,4,3,4>(NE,B,G,J,X,Y,sdim);
+ case 0x244: return Derivatives2D<L,P,2,4,4,4>(NE,B,G,J,X,Y,sdim);
+ case 0x246: return Derivatives2D<L,P,2,4,6,4>(NE,B,G,J,X,Y,sdim);
+ case 0x258: return Derivatives2D<L,P,2,5,8,2>(NE,B,G,J,X,Y,sdim);
  default:
  {
  constexpr int MD = MAX_D1D;
@@ -79,7 +79,7 @@ void TensorPhysDerivatives<QVectorLayout::byNODES>(const int NE,
  << " are not supported!");
  MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
  << MQ << " 1D points are not supported!");
- Derivatives2D<L,P,0,0,0,0,MD,MQ>(NE,B,G,J,X,Y,vdim,D1D,Q1D);
+ Derivatives2D<L,P,0,0,0,0,MD,MQ>(NE,B,G,J,X,Y,sdim,vdim,D1D,Q1D);
  return;
  }
  }

fem/qinterp/grad_phys_by_vdim.cpp

@@ -37,7 +37,7 @@ void TensorPhysDerivatives<QVectorLayout::byVDIM>(const int NE,
  const int D1D = maps.ndof;
  const int Q1D = maps.nqpt;
 
- MFEM_ASSERT(geom.mesh->SpaceDimension() == dim, "");
+ const int sdim = geom.mesh->SpaceDimension();
 
  const double *B = maps.B.Read();
  const double *G = maps.G.Read();
@@ -52,20 +52,20 @@ void TensorPhysDerivatives<QVectorLayout::byVDIM>(const int NE,
 
  if (dim == 1)
  {
- return Derivatives1D<L,P>(NE,G,J,X,Y,vdim,D1D,Q1D);
+ return Derivatives1D<L,P>(NE,G,J,X,Y,sdim,vdim,D1D,Q1D);
  }
  if (dim == 2)
  {
  switch (id)
  {
- case 0x134: return Derivatives2D<L,P,1,3,4,8>(NE,B,G,J,X,Y);
- case 0x146: return Derivatives2D<L,P,1,4,6,4>(NE,B,G,J,X,Y);
- case 0x158: return Derivatives2D<L,P,1,5,8,2>(NE,B,G,J,X,Y);
+ case 0x134: return Derivatives2D<L,P,1,3,4,8>(NE,B,G,J,X,Y,sdim);
+ case 0x146: return Derivatives2D<L,P,1,4,6,4>(NE,B,G,J,X,Y,sdim);
+ case 0x158: return Derivatives2D<L,P,1,5,8,2>(NE,B,G,J,X,Y,sdim);
 
- case 0x233: return Derivatives2D<L,P,2,3,3,8>(NE,B,G,J,X,Y);
- case 0x234: return Derivatives2D<L,P,2,3,4,8>(NE,B,G,J,X,Y);
- case 0x246: return Derivatives2D<L,P,2,4,6,4>(NE,B,G,J,X,Y);
- case 0x258: return Derivatives2D<L,P,2,5,8,2>(NE,B,G,J,X,Y);
+ case 0x233: return Derivatives2D<L,P,2,3,3,8>(NE,B,G,J,X,Y,sdim);
+ case 0x234: return Derivatives2D<L,P,2,3,4,8>(NE,B,G,J,X,Y,sdim);
+ case 0x246: return Derivatives2D<L,P,2,4,6,4>(NE,B,G,J,X,Y,sdim);
+ case 0x258: return Derivatives2D<L,P,2,5,8,2>(NE,B,G,J,X,Y,sdim);
  default:
  {
  constexpr int MD = MAX_D1D;
@@ -74,7 +74,7 @@ void TensorPhysDerivatives<QVectorLayout::byVDIM>(const int NE,
  << " are not supported!");
  MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
  << MQ << " 1D points are not supported!");
- Derivatives2D<L,P,0,0,0,0,MD,MQ>(NE,B,G,J,X,Y,vdim,D1D,Q1D);
+ Derivatives2D<L,P,0,0,0,0,MD,MQ>(NE,B,G,J,X,Y,sdim,vdim,D1D,Q1D);
  return;
  }
  }