4.3/bilininteg__convection__pa_8cpp_source.html

 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced

 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files

 // LICENSE and NOTICE for details. LLNL-CODE-806117.

 //

 // This file is part of the MFEM library. For more information and source code

 // availability visit https://mfem.org.

 //

 // MFEM is free software; you can redistribute it and/or modify it under the

 // terms of the BSD-3 license. We welcome feedback and contributions, see file

 // CONTRIBUTING.md for details.


 #include "../general/forall.hpp"

 #include "bilininteg.hpp"

 #include "gridfunc.hpp"

 #include "ceed/convection.hpp"

 #include "quadinterpolator.hpp"


 namespace mfem

 {


 // PA Convection Integrator


 // PA Convection Assemble 2D kernel

 static void PAConvectionSetup2D(const int NQ,

                                 const int NE,

                                 const Array<double> &w,

                                 const Vector &j,

                                 const Vector &vel,

                                 const double alpha,

                                 Vector &op)

 {

    constexpr int DIM = 2;


    const bool const_v = vel.Size() == DIM;


    const auto W = w.Read();

    const auto J = Reshape(j.Read(), NQ,DIM,DIM,NE);

    const auto V = const_v ?

                   Reshape(vel.Read(), DIM,1,1) :

                   Reshape(vel.Read(), DIM,NQ,NE);

    auto y = Reshape(op.Write(), NQ,DIM,NE);


    MFEM_FORALL(q_global, NE*NQ,

    {

       const int e = q_global / NQ;

       const int q = q_global % NQ;

       const double J11 = J(q,0,0,e);

       const double J21 = J(q,1,0,e);

       const double J12 = J(q,0,1,e);

       const double J22 = J(q,1,1,e);

       const double w = alpha * W[q];

       const double v0 = const_v ? V(0,0,0) : V(0,q,e);

       const double v1 = const_v ? V(1,0,0) : V(1,q,e);

       const double wx = w * v0;

       const double wy = w * v1;

       // y = alpha * W * det(J) * J^{-1} . v = adj(J) . { wx, wy }

       y(q,0,e) =  wx * J22 - wy * J12; // 1

       y(q,1,e) = -wx * J21 + wy * J11; // 2

    });

 }


 // PA Convection Assemble 3D kernel

 static void PAConvectionSetup3D(const int NQ,

                                 const int NE,

                                 const Array<double> &w,

                                 const Vector &j,

                                 const Vector &vel,

                                 const double alpha,

                                 Vector &op)

 {

    constexpr int DIM = 3;

    constexpr int SDIM = DIM;

    const auto W = Reshape(w.Read(), NQ);

    const auto J = Reshape(j.Read(), NQ,SDIM,DIM,NE);

    const bool const_v = vel.Size() == DIM;

    const auto V = const_v ?

                   Reshape(vel.Read(), 3,1,1) :

                   Reshape(vel.Read(), 3,NQ,NE);

    auto y = Reshape(op.Write(), NQ,3,NE);

    MFEM_FORALL(q_global, NE*NQ,

    {

       const int e = q_global / NQ;

       const int q = q_global % NQ;

       const double J11 = J(q,0,0,e);

       const double J12 = J(q,0,1,e);

       const double J13 = J(q,0,2,e);

       const double J21 = J(q,1,0,e);

       const double J22 = J(q,1,1,e);

       const double J23 = J(q,1,2,e);

       const double J31 = J(q,2,0,e);

       const double J32 = J(q,2,1,e);

       const double J33 = J(q,2,2,e);

       const double w = alpha * W(q);

       const double v0 = const_v ? V(0,0,0) : V(0,q,e);

       const double v1 = const_v ? V(1,0,0) : V(1,q,e);

       const double v2 = const_v ? V(2,0,0) : V(2,q,e);

       const double wx = w * v0;

       const double wy = w * v1;

       const double wz = w * v2;

       // A = adj(J)

       const double A11 = (J22 * J33) - (J23 * J32);

       const double A12 = (J32 * J13) - (J12 * J33);

       const double A13 = (J12 * J23) - (J22 * J13);

       const double A21 = (J31 * J23) - (J21 * J33);

       const double A22 = (J11 * J33) - (J13 * J31);

       const double A23 = (J21 * J13) - (J11 * J23);

       const double A31 = (J21 * J32) - (J31 * J22);

       const double A32 = (J31 * J12) - (J11 * J32);

       const double A33 = (J11 * J22) - (J12 * J21);

       // y = alpha * W * det(J) * J^{-1} . v = adj(J) . { wx, wy, wz }

       y(q,0,e) = wx * A11 + wy * A12 + wz * A13;

       y(q,1,e) = wx * A21 + wy * A22 + wz * A23;

       y(q,2,e) = wx * A31 + wy * A32 + wz * A33;

    });

 }


 static void PAConvectionSetup(const int dim,

                               const int NQ,

                               const int NE,

                               const Array<double> &W,

                               const Vector &J,

                               const Vector &coeff,

                               const double alpha,

                               Vector &op)

 {

    if (dim == 1) { MFEM_ABORT("dim==1 not supported in PAConvectionSetup"); }

    if (dim == 2)

    {

       PAConvectionSetup2D(NQ, NE, W, J, coeff, alpha, op);

    }

    if (dim == 3)

    {

       PAConvectionSetup3D(NQ, NE, W, J, coeff, alpha, op);

    }

 }


 // PA Convection Apply 2D kernel

 template<int T_D1D = 0, int T_Q1D = 0> static

 void PAConvectionApply2D(const int ne,

                          const Array<double> &b,

                          const Array<double> &g,

                          const Array<double> &bt,

                          const Array<double> &gt,

                          const Vector &op_,

                          const Vector &x_,

                          Vector &y_,

                          const int d1d = 0,

                          const int q1d = 0)

 {

    const int NE = ne;

    const int D1D = T_D1D ? T_D1D : d1d;

    const int Q1D = T_Q1D ? T_Q1D : q1d;

    MFEM_VERIFY(D1D <= MAX_D1D, "");

    MFEM_VERIFY(Q1D <= MAX_Q1D, "");

    auto B = Reshape(b.Read(), Q1D, D1D);

    auto G = Reshape(g.Read(), Q1D, D1D);

    auto Bt = Reshape(bt.Read(), D1D, Q1D);

    auto op = Reshape(op_.Read(), Q1D, Q1D, 2, NE);

    auto x = Reshape(x_.Read(), D1D, D1D, NE);

    auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);

    MFEM_FORALL(e, NE,

    {

       const int D1D = T_D1D ? T_D1D : d1d;

       const int Q1D = T_Q1D ? T_Q1D : q1d;

       // the following variables are evaluated at compile time

       constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;

       constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;


       double u[max_D1D][max_D1D];

       for (int dy = 0; dy < D1D; ++dy)

       {

          for (int dx = 0; dx < D1D; ++dx)

          {

             u[dy][dx] = x(dx,dy,e);

          }

       }

       double Bu[max_D1D][max_Q1D];

       double Gu[max_D1D][max_Q1D];

       for (int dy = 0; dy < D1D; ++dy)

       {

          for (int qx = 0; qx < Q1D; ++qx)

          {

             Bu[dy][qx] = 0.0;

             Gu[dy][qx] = 0.0;

             for (int dx = 0; dx < D1D; ++dx)

             {

                const double bx = B(qx,dx);

                const double gx = G(qx,dx);

                const double x = u[dy][dx];

                Bu[dy][qx] += bx * x;

                Gu[dy][qx] += gx * x;

             }

          }

       }

       double GBu[max_Q1D][max_Q1D];

       double BGu[max_Q1D][max_Q1D];

       for (int qx = 0; qx < Q1D; ++qx)

       {

          for (int qy = 0; qy < Q1D; ++qy)

          {

             GBu[qy][qx] = 0.0;

             BGu[qy][qx] = 0.0;

             for (int dy = 0; dy < D1D; ++dy)

             {

                const double bx = B(qy,dy);

                const double gx = G(qy,dy);

                GBu[qy][qx] += gx * Bu[dy][qx];

                BGu[qy][qx] += bx * Gu[dy][qx];

             }

          }

       }

       // Calculate Dxy, xDy in plane

       double DGu[max_Q1D][max_Q1D];

       for (int qy = 0; qy < Q1D; ++qy)

       {

          for (int qx = 0; qx < Q1D; ++qx)

          {

             const double O1 = op(qx,qy,0,e);

             const double O2 = op(qx,qy,1,e);


             const double gradX = BGu[qy][qx];

             const double gradY = GBu[qy][qx];


             DGu[qy][qx] = (O1 * gradX) + (O2 * gradY);

          }

       }

       double BDGu[max_D1D][max_Q1D];

       for (int qx = 0; qx < Q1D; ++qx)

       {

          for (int dy = 0; dy < D1D; ++dy)

          {

             BDGu[dy][qx] = 0.0;

             for (int qy = 0; qy < Q1D; ++qy)

             {

                const double w = Bt(dy,qy);

                BDGu[dy][qx] += w * DGu[qy][qx];

             }

          }

       }

       for (int dx = 0; dx < D1D; ++dx)

       {

          for (int dy = 0; dy < D1D; ++dy)

          {

             double BBDGu = 0.0;

             for (int qx = 0; qx < Q1D; ++qx)

             {

                const double w = Bt(dx,qx);

                BBDGu += w * BDGu[dy][qx];

             }

             y(dx,dy,e) += BBDGu;

          }

       }

    });

 }


 // Optimized PA Convection Apply 2D kernel

 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static

 void SmemPAConvectionApply2D(const int ne,

                              const Array<double> &b,

                              const Array<double> &g,

                              const Array<double> &bt,

                              const Array<double> &gt,

                              const Vector &op_,

                              const Vector &x_,

                              Vector &y_,

                              const int d1d = 0,

                              const int q1d = 0)

 {

    const int NE = ne;

    const int D1D = T_D1D ? T_D1D : d1d;

    const int Q1D = T_Q1D ? T_Q1D : q1d;

    constexpr int NBZ = T_NBZ ? T_NBZ : 1;

    MFEM_VERIFY(D1D <= MAX_D1D, "");

    MFEM_VERIFY(Q1D <= MAX_Q1D, "");

    auto B = Reshape(b.Read(), Q1D, D1D);

    auto G = Reshape(g.Read(), Q1D, D1D);

    auto Bt = Reshape(bt.Read(), D1D, Q1D);

    auto op = Reshape(op_.Read(), Q1D, Q1D, 2, NE);

    auto x = Reshape(x_.Read(), D1D, D1D, NE);

    auto y = Reshape(y_.ReadWrite(), D1D, D1D, NE);

    MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,

    {

       const int tidz = MFEM_THREAD_ID(z);

       const int D1D = T_D1D ? T_D1D : d1d;

       const int Q1D = T_Q1D ? T_Q1D : q1d;

       // the following variables are evaluated at compile time

       constexpr int NBZ = T_NBZ ? T_NBZ : 1;

       constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;

       constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;

       // constexpr int MDQ = (max_Q1D > max_D1D) ? max_Q1D : max_D1D;

       MFEM_SHARED double u[NBZ][max_D1D][max_D1D];

       MFEM_FOREACH_THREAD(dy,y,D1D)

       {

          MFEM_FOREACH_THREAD(dx,x,D1D)

          {

             // e is really equal to e+tidz

             u[tidz][dy][dx] = x(dx,dy,e);

          }

       }

       MFEM_SYNC_THREAD;

       MFEM_SHARED double Bu[NBZ][max_D1D][max_Q1D];

       MFEM_SHARED double Gu[NBZ][max_D1D][max_Q1D];

       MFEM_FOREACH_THREAD(dy,y,D1D)

       {

          MFEM_FOREACH_THREAD(qx,x,Q1D)

          {

             Bu[tidz][dy][qx] = 0.0;

             Gu[tidz][dy][qx] = 0.0;

             for (int dx = 0; dx < D1D; ++dx)

             {

                const double bx = B(qx,dx);

                const double gx = G(qx,dx);

                const double x = u[tidz][dy][dx];

                Bu[tidz][dy][qx] += bx * x;

                Gu[tidz][dy][qx] += gx * x;

             }

          }

       }

       MFEM_SYNC_THREAD;

       MFEM_SHARED double GBu[NBZ][max_Q1D][max_Q1D];

       MFEM_SHARED double BGu[NBZ][max_Q1D][max_Q1D];

       MFEM_FOREACH_THREAD(qx,x,Q1D)

       {

          MFEM_FOREACH_THREAD(qy,y,Q1D)

          {

             GBu[tidz][qy][qx] = 0.0;

             BGu[tidz][qy][qx] = 0.0;

             for (int dy = 0; dy < D1D; ++dy)

             {

                const double bx = B(qy,dy);

                const double gx = G(qy,dy);

                GBu[tidz][qy][qx] += gx * Bu[tidz][dy][qx];

                BGu[tidz][qy][qx] += bx * Gu[tidz][dy][qx];

             }

          }

       }

       MFEM_SYNC_THREAD;

       // Calculate Dxy, xDy in plane

       MFEM_SHARED double DGu[NBZ][max_Q1D][max_Q1D];

       MFEM_FOREACH_THREAD(qy,y,Q1D)

       {

          MFEM_FOREACH_THREAD(qx,x,Q1D)

          {

             const double O1 = op(qx,qy,0,e);

             const double O2 = op(qx,qy,1,e);


             const double gradX = BGu[tidz][qy][qx];

             const double gradY = GBu[tidz][qy][qx];


             DGu[tidz][qy][qx] = (O1 * gradX) + (O2 * gradY);

          }

       }

       MFEM_SYNC_THREAD;

       MFEM_SHARED double BDGu[NBZ][max_D1D][max_Q1D];

       MFEM_FOREACH_THREAD(qx,x,Q1D)

       {

          MFEM_FOREACH_THREAD(dy,y,D1D)

          {

             BDGu[tidz][dy][qx] = 0.0;

             for (int qy = 0; qy < Q1D; ++qy)

             {

                const double w = Bt(dy,qy);

                BDGu[tidz][dy][qx] += w * DGu[tidz][qy][qx];

             }

          }

       }

       MFEM_SYNC_THREAD;

       MFEM_FOREACH_THREAD(dx,x,D1D)

       {

          MFEM_FOREACH_THREAD(dy,y,D1D)

          {

             double BBDGu = 0.0;

             for (int qx = 0; qx < Q1D; ++qx)

             {

                const double w = Bt(dx,qx);

                BBDGu += w * BDGu[tidz][dy][qx];

             }

             y(dx,dy,e) += BBDGu;

          }

       }

    });

 }


 // PA Convection Apply 3D kernel

 template<int T_D1D = 0, int T_Q1D = 0> static

 void PAConvectionApply3D(const int ne,

                          const Array<double> &b,

                          const Array<double> &g,

                          const Array<double> &bt,

                          const Array<double> &gt,

                          const Vector &op_,

                          const Vector &x_,

                          Vector &y_,

                          const int d1d = 0,

                          const int q1d = 0)

 {

    const int NE = ne;

    const int D1D = T_D1D ? T_D1D : d1d;

    const int Q1D = T_Q1D ? T_Q1D : q1d;

    MFEM_VERIFY(D1D <= MAX_D1D, "");

    MFEM_VERIFY(Q1D <= MAX_Q1D, "");

    auto B = Reshape(b.Read(), Q1D, D1D);

    auto G = Reshape(g.Read(), Q1D, D1D);

    auto Bt = Reshape(bt.Read(), D1D, Q1D);

    auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, 3, NE);

    auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);

    auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);

    MFEM_FORALL(e, NE,

    {

       const int D1D = T_D1D ? T_D1D : d1d;

       const int Q1D = T_Q1D ? T_Q1D : q1d;

       // the following variables are evaluated at compile time

       constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;

       constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;


       double u[max_D1D][max_D1D][max_D1D];

       for (int dz = 0; dz < D1D; ++dz)

       {

          for (int dy = 0; dy < D1D; ++dy)

          {

             for (int dx = 0; dx < D1D; ++dx)

             {

                u[dz][dy][dx] = x(dx,dy,dz,e);

             }

          }

       }

       double Bu[max_D1D][max_D1D][max_Q1D];

       double Gu[max_D1D][max_D1D][max_Q1D];

       for (int dz = 0; dz < D1D; ++dz)

       {

          for (int dy = 0; dy < D1D; ++dy)

          {

             for (int qx = 0; qx < Q1D; ++qx)

             {

                Bu[dz][dy][qx] = 0.0;

                Gu[dz][dy][qx] = 0.0;

                for (int dx = 0; dx < D1D; ++dx)

                {

                   const double bx = B(qx,dx);

                   const double gx = G(qx,dx);

                   const double x = u[dz][dy][dx];

                   Bu[dz][dy][qx] += bx * x;

                   Gu[dz][dy][qx] += gx * x;

                }

             }

          }

       }

       double BBu[max_D1D][max_Q1D][max_Q1D];

       double GBu[max_D1D][max_Q1D][max_Q1D];

       double BGu[max_D1D][max_Q1D][max_Q1D];

       for (int dz = 0; dz < D1D; ++dz)

       {

          for (int qx = 0; qx < Q1D; ++qx)

          {

             for (int qy = 0; qy < Q1D; ++qy)

             {

                BBu[dz][qy][qx] = 0.0;

                GBu[dz][qy][qx] = 0.0;

                BGu[dz][qy][qx] = 0.0;

                for (int dy = 0; dy < D1D; ++dy)

                {

                   const double bx = B(qy,dy);

                   const double gx = G(qy,dy);

                   BBu[dz][qy][qx] += bx * Bu[dz][dy][qx];

                   GBu[dz][qy][qx] += gx * Bu[dz][dy][qx];

                   BGu[dz][qy][qx] += bx * Gu[dz][dy][qx];

                }

             }

          }

       }

       double GBBu[max_Q1D][max_Q1D][max_Q1D];

       double BGBu[max_Q1D][max_Q1D][max_Q1D];

       double BBGu[max_Q1D][max_Q1D][max_Q1D];

       for (int qx = 0; qx < Q1D; ++qx)

       {

          for (int qy = 0; qy < Q1D; ++qy)

          {

             for (int qz = 0; qz < Q1D; ++qz)

             {

                GBBu[qz][qy][qx] = 0.0;

                BGBu[qz][qy][qx] = 0.0;

                BBGu[qz][qy][qx] = 0.0;

                for (int dz = 0; dz < D1D; ++dz)

                {

                   const double bx = B(qz,dz);

                   const double gx = G(qz,dz);

                   GBBu[qz][qy][qx] += gx * BBu[dz][qy][qx];

                   BGBu[qz][qy][qx] += bx * GBu[dz][qy][qx];

                   BBGu[qz][qy][qx] += bx * BGu[dz][qy][qx];

                }

             }

          }

       }

       // Calculate Dxy, xDy in plane

       double DGu[max_Q1D][max_Q1D][max_Q1D];

       for (int qz = 0; qz < Q1D; ++qz)

       {

          for (int qy = 0; qy < Q1D; ++qy)

          {

             for (int qx = 0; qx < Q1D; ++qx)

             {

                const double O1 = op(qx,qy,qz,0,e);

                const double O2 = op(qx,qy,qz,1,e);

                const double O3 = op(qx,qy,qz,2,e);


                const double gradX = BBGu[qz][qy][qx];

                const double gradY = BGBu[qz][qy][qx];

                const double gradZ = GBBu[qz][qy][qx];


                DGu[qz][qy][qx] = (O1 * gradX) + (O2 * gradY) + (O3 * gradZ);

             }

          }

       }

       double BDGu[max_D1D][max_Q1D][max_Q1D];

       for (int qx = 0; qx < Q1D; ++qx)

       {

          for (int qy = 0; qy < Q1D; ++qy)

          {

             for (int dz = 0; dz < D1D; ++dz)

             {

                BDGu[dz][qy][qx] = 0.0;

                for (int qz = 0; qz < Q1D; ++qz)

                {

                   const double w = Bt(dz,qz);

                   BDGu[dz][qy][qx] += w * DGu[qz][qy][qx];

                }

             }

          }

       }

       double BBDGu[max_D1D][max_D1D][max_Q1D];

       for (int dz = 0; dz < D1D; ++dz)

       {

          for (int qx = 0; qx < Q1D; ++qx)

          {

             for (int dy = 0; dy < D1D; ++dy)

             {

                BBDGu[dz][dy][qx] = 0.0;

                for (int qy = 0; qy < Q1D; ++qy)

                {

                   const double w = Bt(dy,qy);

                   BBDGu[dz][dy][qx] += w * BDGu[dz][qy][qx];

                }

             }

          }

       }

       for (int dz = 0; dz < D1D; ++dz)

       {

          for (int dy = 0; dy < D1D; ++dy)

          {

             for (int dx = 0; dx < D1D; ++dx)

             {

                double BBBDGu = 0.0;

                for (int qx = 0; qx < Q1D; ++qx)

                {

                   const double w = Bt(dx,qx);

                   BBBDGu += w * BBDGu[dz][dy][qx];

                }

                y(dx,dy,dz,e) += BBBDGu;

             }

          }

       }

    });

 }


 // Optimized PA Convection Apply 3D kernel

 template<int T_D1D = 0, int T_Q1D = 0> static

 void SmemPAConvectionApply3D(const int ne,

                              const Array<double> &b,

                              const Array<double> &g,

                              const Array<double> &bt,

                              const Array<double> &gt,

                              const Vector &op_,

                              const Vector &x_,

                              Vector &y_,

                              const int d1d = 0,

                              const int q1d = 0)

 {

    const int NE = ne;

    const int D1D = T_D1D ? T_D1D : d1d;

    const int Q1D = T_Q1D ? T_Q1D : q1d;

    MFEM_VERIFY(D1D <= MAX_D1D, "");

    MFEM_VERIFY(Q1D <= MAX_Q1D, "");

    auto B = Reshape(b.Read(), Q1D, D1D);

    auto G = Reshape(g.Read(), Q1D, D1D);

    auto Bt = Reshape(bt.Read(), D1D, Q1D);

    auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, 3, NE);

    auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);

    auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);

    MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,

    {

       const int D1D = T_D1D ? T_D1D : d1d;

       const int Q1D = T_Q1D ? T_Q1D : q1d;

       // the following variables are evaluated at compile time

       constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;

       constexpr int max_Q1D = T_Q1D ? T_Q1D : MAX_Q1D;

       constexpr int max_DQ = (max_Q1D > max_D1D) ? max_Q1D : max_D1D;

       MFEM_SHARED double sm0[max_DQ*max_DQ*max_DQ];

       MFEM_SHARED double sm1[max_DQ*max_DQ*max_DQ];

       MFEM_SHARED double sm2[max_DQ*max_DQ*max_DQ];

       MFEM_SHARED double sm3[max_DQ*max_DQ*max_DQ];

       MFEM_SHARED double sm4[max_DQ*max_DQ*max_DQ];

       MFEM_SHARED double sm5[max_DQ*max_DQ*max_DQ];


       double (*u)[max_D1D][max_D1D] = (double (*)[max_D1D][max_D1D]) sm0;

       MFEM_FOREACH_THREAD(dz,z,D1D)

       {

          MFEM_FOREACH_THREAD(dy,y,D1D)

          {

             MFEM_FOREACH_THREAD(dx,x,D1D)

             {

                u[dz][dy][dx] = x(dx,dy,dz,e);

             }

          }

       }

       MFEM_SYNC_THREAD;

       double (*Bu)[max_D1D][max_Q1D] = (double (*)[max_D1D][max_Q1D])sm1;

       double (*Gu)[max_D1D][max_Q1D] = (double (*)[max_D1D][max_Q1D])sm2;

       MFEM_FOREACH_THREAD(dz,z,D1D)

       {

          MFEM_FOREACH_THREAD(dy,y,D1D)

          {

             MFEM_FOREACH_THREAD(qx,x,Q1D)

             {

                double Bu_ = 0.0;

                double Gu_ = 0.0;

                for (int dx = 0; dx < D1D; ++dx)

                {

                   const double bx = B(qx,dx);

                   const double gx = G(qx,dx);

                   const double x = u[dz][dy][dx];

                   Bu_ += bx * x;

                   Gu_ += gx * x;

                }

                Bu[dz][dy][qx] = Bu_;

                Gu[dz][dy][qx] = Gu_;

             }

          }

       }

       MFEM_SYNC_THREAD;

       double (*BBu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm3;

       double (*GBu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm4;

       double (*BGu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm5;

       MFEM_FOREACH_THREAD(dz,z,D1D)

       {

          MFEM_FOREACH_THREAD(qx,x,Q1D)

          {

             MFEM_FOREACH_THREAD(qy,y,Q1D)

             {

                double BBu_ = 0.0;

                double GBu_ = 0.0;

                double BGu_ = 0.0;

                for (int dy = 0; dy < D1D; ++dy)

                {

                   const double bx = B(qy,dy);

                   const double gx = G(qy,dy);

                   BBu_ += bx * Bu[dz][dy][qx];

                   GBu_ += gx * Bu[dz][dy][qx];

                   BGu_ += bx * Gu[dz][dy][qx];

                }

                BBu[dz][qy][qx] = BBu_;

                GBu[dz][qy][qx] = GBu_;

                BGu[dz][qy][qx] = BGu_;

             }

          }

       }

       MFEM_SYNC_THREAD;

       double (*GBBu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm0;

       double (*BGBu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm1;

       double (*BBGu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm2;

       MFEM_FOREACH_THREAD(qx,x,Q1D)

       {

          MFEM_FOREACH_THREAD(qy,y,Q1D)

          {

             MFEM_FOREACH_THREAD(qz,z,Q1D)

             {

                double GBBu_ = 0.0;

                double BGBu_ = 0.0;

                double BBGu_ = 0.0;

                for (int dz = 0; dz < D1D; ++dz)

                {

                   const double bx = B(qz,dz);

                   const double gx = G(qz,dz);

                   GBBu_ += gx * BBu[dz][qy][qx];

                   BGBu_ += bx * GBu[dz][qy][qx];

                   BBGu_ += bx * BGu[dz][qy][qx];

                }

                GBBu[qz][qy][qx] = GBBu_;

                BGBu[qz][qy][qx] = BGBu_;

                BBGu[qz][qy][qx] = BBGu_;

             }

          }

       }

       MFEM_SYNC_THREAD;

       double (*DGu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm3;

       MFEM_FOREACH_THREAD(qz,z,Q1D)

       {

          MFEM_FOREACH_THREAD(qy,y,Q1D)

          {

             MFEM_FOREACH_THREAD(qx,x,Q1D)

             {

                const double O1 = op(qx,qy,qz,0,e);

                const double O2 = op(qx,qy,qz,1,e);

                const double O3 = op(qx,qy,qz,2,e);


                const double gradX = BBGu[qz][qy][qx];

                const double gradY = BGBu[qz][qy][qx];

                const double gradZ = GBBu[qz][qy][qx];


                DGu[qz][qy][qx] = (O1 * gradX) + (O2 * gradY) + (O3 * gradZ);

             }

          }

       }

       MFEM_SYNC_THREAD;

       double (*BDGu)[max_Q1D][max_Q1D] = (double (*)[max_Q1D][max_Q1D])sm4;

       MFEM_FOREACH_THREAD(qx,x,Q1D)

       {

          MFEM_FOREACH_THREAD(qy,y,Q1D)

          {

             MFEM_FOREACH_THREAD(dz,z,D1D)

             {

                double BDGu_ = 0.0;

                for (int qz = 0; qz < Q1D; ++qz)

                {

                   const double w = Bt(dz,qz);

                   BDGu_ += w * DGu[qz][qy][qx];

                }

                BDGu[dz][qy][qx] = BDGu_;

             }

          }

       }

       MFEM_SYNC_THREAD;

       double (*BBDGu)[max_D1D][max_Q1D] = (double (*)[max_D1D][max_Q1D])sm5;

       MFEM_FOREACH_THREAD(dz,z,D1D)

       {

          MFEM_FOREACH_THREAD(qx,x,Q1D)

          {

             MFEM_FOREACH_THREAD(dy,y,D1D)

             {

                double BBDGu_ = 0.0;

                for (int qy = 0; qy < Q1D; ++qy)

                {

                   const double w = Bt(dy,qy);

                   BBDGu_ += w * BDGu[dz][qy][qx];

                }

                BBDGu[dz][dy][qx] = BBDGu_;

             }

          }

       }

       MFEM_SYNC_THREAD;

       MFEM_FOREACH_THREAD(dz,z,D1D)

       {

          MFEM_FOREACH_THREAD(dy,y,D1D)

          {

             MFEM_FOREACH_THREAD(dx,x,D1D)

             {

                double BBBDGu = 0.0;

                for (int qx = 0; qx < Q1D; ++qx)

                {

                   const double w = Bt(dx,qx);

                   BBBDGu += w * BBDGu[dz][dy][qx];

                }

                y(dx,dy,dz,e) = BBBDGu;

             }

          }

       }

    });

 }


 void ConvectionIntegrator::AssemblePA(const FiniteElementSpace &fes)

 {

    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?

                          Device::GetDeviceMemoryType() : pa_mt;

    // Assumes tensor-product elements

    Mesh *mesh = fes.GetMesh();

    const FiniteElement &el = *fes.GetFE(0);

    ElementTransformation &Trans = *fes.GetElementTransformation(0);

    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, Trans);

    if (DeviceCanUseCeed())

    {

       delete ceedOp;

       ceedOp = new ceed::PAConvectionIntegrator(fes, *ir, Q, alpha);

       return;

    }

    const int dims = el.GetDim();

    const int symmDims = dims;

    const int nq = ir->GetNPoints();

    dim = mesh->Dimension();

    ne = fes.GetNE();

    geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS, mt);

    maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);

    dofs1D = maps->ndof;

    quad1D = maps->nqpt;

    pa_data.SetSize(symmDims * nq * ne, mt);

    Vector vel;

    if (VectorConstantCoefficient *cQ =

           dynamic_cast<VectorConstantCoefficient*>(Q))

    {

       vel = cQ->GetVec();

    }

    else if (VectorGridFunctionCoefficient *vgfQ =

                dynamic_cast<VectorGridFunctionCoefficient*>(Q))

    {

       vel.SetSize(dim * nq * ne, mt);


       const GridFunction *gf = vgfQ->GetGridFunction();

       const FiniteElementSpace &gf_fes = *gf->FESpace();

       const QuadratureInterpolator *qi(gf_fes.GetQuadratureInterpolator(*ir));

       const bool use_tensor_products = UsesTensorBasis(gf_fes);

       const ElementDofOrdering ordering = use_tensor_products ?

                                           ElementDofOrdering::LEXICOGRAPHIC :

                                           ElementDofOrdering::NATIVE;

       const Operator *R = gf_fes.GetElementRestriction(ordering);


       Vector xe(R->Height(), mt);

       xe.UseDevice(true);


       R->Mult(*gf, xe);

       qi->SetOutputLayout(QVectorLayout::byVDIM);

       qi->DisableTensorProducts(!use_tensor_products);

       qi->Values(xe,vel);

    }

    else if (VectorQuadratureFunctionCoefficient* cQ =

                dynamic_cast<VectorQuadratureFunctionCoefficient*>(Q))

    {

       const QuadratureFunction &qFun = cQ->GetQuadFunction();

       MFEM_VERIFY(qFun.Size() == dim * nq * ne,

                   "Incompatible QuadratureFunction dimension \n");


       MFEM_VERIFY(ir == &qFun.GetSpace()->GetElementIntRule(0),

                   "IntegrationRule used within integrator and in"

                   " QuadratureFunction appear to be different");


       qFun.Read();

       vel.MakeRef(const_cast<QuadratureFunction &>(qFun),0);

    }

    else

    {

       vel.SetSize(dim * nq * ne);

       auto C = Reshape(vel.HostWrite(), dim, nq, ne);

       DenseMatrix Q_ir;

       for (int e = 0; e < ne; ++e)

       {

          ElementTransformation& T = *fes.GetElementTransformation(e);

          Q->Eval(Q_ir, T, *ir);

          for (int q = 0; q < nq; ++q)

          {

             for (int i = 0; i < dim; ++i)

             {

                C(i,q,e) = Q_ir(i,q);

             }

          }

       }

    }

    PAConvectionSetup(dim, nq, ne, ir->GetWeights(), geom->J,

                      vel, alpha, pa_data);

 }


 static void PAConvectionApply(const int dim,

                               const int D1D,

                               const int Q1D,

                               const int NE,

                               const Array<double> &B,

                               const Array<double> &G,

                               const Array<double> &Bt,

                               const Array<double> &Gt,

                               const Vector &op,

                               const Vector &x,

                               Vector &y)

 {

    if (dim == 2)

    {

       switch ((D1D << 4 ) | Q1D)

       {

          case 0x22: return SmemPAConvectionApply2D<2,2,8>(NE,B,G,Bt,Gt,op,x,y);

          case 0x33: return SmemPAConvectionApply2D<3,3,4>(NE,B,G,Bt,Gt,op,x,y);

          case 0x34: return SmemPAConvectionApply2D<3,4,4>(NE,B,G,Bt,Gt,op,x,y);

          case 0x44: return SmemPAConvectionApply2D<4,4,4>(NE,B,G,Bt,Gt,op,x,y);

          case 0x46: return SmemPAConvectionApply2D<4,6,4>(NE,B,G,Bt,Gt,op,x,y);

          case 0x55: return SmemPAConvectionApply2D<5,5,2>(NE,B,G,Bt,Gt,op,x,y);

          case 0x58: return SmemPAConvectionApply2D<5,8,2>(NE,B,G,Bt,Gt,op,x,y);

          case 0x66: return SmemPAConvectionApply2D<6,6,1>(NE,B,G,Bt,Gt,op,x,y);

          case 0x77: return SmemPAConvectionApply2D<7,7,1>(NE,B,G,Bt,Gt,op,x,y);

          case 0x88: return SmemPAConvectionApply2D<8,8,1>(NE,B,G,Bt,Gt,op,x,y);

          case 0x99: return SmemPAConvectionApply2D<9,9,1>(NE,B,G,Bt,Gt,op,x,y);

          default:   return PAConvectionApply2D(NE,B,G,Bt,Gt,op,x,y,D1D,Q1D);

       }

    }

    else if (dim == 3)

    {

       switch ((D1D << 4 ) | Q1D)

       {

          case 0x23: return SmemPAConvectionApply3D<2,3>(NE,B,G,Bt,Gt,op,x,y);

          case 0x24: return SmemPAConvectionApply3D<2,4>(NE,B,G,Bt,Gt,op,x,y);

          case 0x26: return SmemPAConvectionApply3D<2,6>(NE,B,G,Bt,Gt,op,x,y);

          case 0x34: return SmemPAConvectionApply3D<3,4>(NE,B,G,Bt,Gt,op,x,y);

          case 0x35: return SmemPAConvectionApply3D<3,5>(NE,B,G,Bt,Gt,op,x,y);

          case 0x45: return SmemPAConvectionApply3D<4,5>(NE,B,G,Bt,Gt,op,x,y);

          case 0x48: return SmemPAConvectionApply3D<4,8>(NE,B,G,Bt,Gt,op,x,y);

          case 0x56: return SmemPAConvectionApply3D<5,6>(NE,B,G,Bt,Gt,op,x,y);

          case 0x67: return SmemPAConvectionApply3D<6,7>(NE,B,G,Bt,Gt,op,x,y);

          case 0x78: return SmemPAConvectionApply3D<7,8>(NE,B,G,Bt,Gt,op,x,y);

          case 0x89: return SmemPAConvectionApply3D<8,9>(NE,B,G,Bt,Gt,op,x,y);

          default:   return PAConvectionApply3D(NE,B,G,Bt,Gt,op,x,y,D1D,Q1D);

       }

    }

    MFEM_ABORT("Unknown kernel.");

 }


 // PA Convection Apply kernel

 void ConvectionIntegrator::AddMultPA(const Vector &x, Vector &y) const

 {

    if (DeviceCanUseCeed())

    {

       ceedOp->AddMult(x, y);

    }

    else

    {

       PAConvectionApply(dim, dofs1D, quad1D, ne,

                         maps->B, maps->G, maps->Bt, maps->Gt,

                         pa_data, x, y);

    }

 }


 void ConvectionIntegrator::AssembleDiagonalPA(Vector &diag)

 {

    if (DeviceCanUseCeed())

    {

       ceedOp->GetDiagonal(diag);

    }

    else

    {

       MFEM_ABORT("AssembleDiagonalPA not yet implemented for"

                  " ConvectionIntegrator.");

    }

 }


 } // namespace mfem

mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247

mfem::FiniteElement
Abstract class for all finite elements.
Definition: fe.hpp:243

mfem::ConvectionIntegrator::alpha
double alpha
Definition: bilininteg.hpp:2133

mfem::FiniteElement::GetDim
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:317

mfem::Mesh
Definition: mesh.hpp:52

mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90

mfem::GridFunction
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30

mfem::DofToQuad::TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe.hpp:165

mfem::Mesh::GetGeometricFactors
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition: mesh.cpp:784

mfem::QuadratureSpace::GetElementIntRule
const IntegrationRule & GetElementIntRule(int idx) const
Get the IntegrationRule associated with mesh element idx.
Definition: fespace.hpp:950

mfem::VectorCoefficient::Eval
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...

mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513

mfem::ConvectionIntegrator::geom
const GeometricFactors * geom
Not owned.
Definition: bilininteg.hpp:2137

mfem::VectorQuadratureFunctionCoefficient
Vector quadrature function coefficient which requires that the quadrature rules used for this vector ...
Definition: coefficient.hpp:1666

mfem::VectorConstantCoefficient
Vector coefficient that is constant in space and time.
Definition: coefficient.hpp:389

mfem::ConvectionIntegrator::ne
int ne
Definition: bilininteg.hpp:2138

mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23

mfem::DofToQuad::nqpt
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe.hpp:177

mfem::Vector::HostWrite
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:438

mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190

mfem::DofToQuad::Gt
Array< double > Gt
Transpose of G.
Definition: fe.hpp:216

mfem::DofToQuad::ndof
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe.hpp:173

mfem::Operator::Mult
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).

mfem::IntegrationRule::GetWeights
const Array< double > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:83

DIM
constexpr int DIM
Definition: minimal-surface.cpp:71

mfem::ConvectionIntegrator::dim
int dim
Definition: bilininteg.hpp:2138

mfem::UsesTensorBasis
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:957

mfem::NonlinearFormIntegrator::ceedOp
ceed::Operator * ceedOp
Definition: nonlininteg.hpp:33

mfem::Vector::UseDevice
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:108

mfem::ConvectionIntegrator::AssembleDiagonalPA
virtual void AssembleDiagonalPA(Vector &diag)
Assemble diagonal and add it to Vector diag.
Definition: bilininteg_convection_pa.cpp:924

mfem::MemoryType::DEFAULT

mfem::FiniteElementSpace::GetElementRestriction
const Operator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1195

mfem::Reshape
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:136

mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:567

mfem::GeometricFactors::J
Vector J
Jacobians of the element transformations at all quadrature points.
Definition: mesh.hpp:1596

mfem::ElementDofOrdering::NATIVE
Native ordering as defined by the FiniteElement.

mfem::MAX_Q1D
const int MAX_Q1D
Definition: forall.hpp:28

mfem::Operator::Height
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66

mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:410

mfem::QuadratureInterpolator
A class that performs interpolation from an E-vector to quadrature point values and/or derivatives (Q...
Definition: quadinterpolator.hpp:36

b
double b
Definition: lissajous.cpp:42

mfem::Array< double >

mfem::Device::GetDeviceMemoryType
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273

mfem::QuadratureFunction::GetSpace
QuadratureSpace * GetSpace() const
Get the associated QuadratureSpace.
Definition: gridfunc.hpp:775

mfem::GridFunction::FESpace
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:629

mfem::Mesh::Dimension
int Dimension() const
Definition: mesh.hpp:911

mfem::superlu::Trans
Trans
Definition: superlu.hpp:43

mfem::ConvectionIntegrator::AddMultPA
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
Definition: bilininteg_convection_pa.cpp:910

bilininteg.hpp

mfem::FiniteElementSpace::GetElementTransformation
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:600

mfem::DofToQuad::Bt
Array< double > Bt
Transpose of B.
Definition: fe.hpp:194

mfem::ConvectionIntegrator::AssemblePA
virtual void AssemblePA(const FiniteElementSpace &)
Method defining partial assembly.
Definition: bilininteg_convection_pa.cpp:769

mfem::ConvectionIntegrator::dofs1D
int dofs1D
Definition: bilininteg.hpp:2138

mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87

mfem::ConvectionIntegrator::pa_data
Vector pa_data
Definition: bilininteg.hpp:2135

mfem::MemoryType
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31

mfem::ConvectionIntegrator::GetRule
static const IntegrationRule & GetRule(const FiniteElement &el, ElementTransformation &Trans)
Definition: bilininteg.cpp:1380

mfem::NonlinearFormIntegrator::IntRule
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:30

mfem::DofToQuad::B
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe.hpp:188

mfem::FiniteElement::GetDofToQuad
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe.cpp:367

mfem::ElementTransformation
Definition: eltrans.hpp:23

mfem::GeometricFactors::JACOBIANS
Definition: mesh.hpp:1570

vel
void vel(const Vector &x, double t, Vector &u)
Definition: navier_3dfoc.cpp:28

mfem::FiniteElementSpace::GetQuadratureInterpolator
const QuadratureInterpolator * GetQuadratureInterpolator(const IntegrationRule &ir) const
Return a QuadratureInterpolator that interpolates E-vectors to quadrature point values and/or derivat...
Definition: fespace.cpp:1256

mfem::ElementDofOrdering
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:65

dim
int dim
Definition: ex24.cpp:53

mfem::MAX_D1D
const int MAX_D1D
Definition: forall.hpp:27

quadinterpolator.hpp

mfem::DofToQuad::G
Array< double > G
Gradients/divergences/curls of basis functions evaluated at quadrature points.
Definition: fe.hpp:209

mfem::FiniteElementSpace::GetFE
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2388

SDIM
constexpr int SDIM
Definition: minimal-surface.cpp:72

mfem::ElementDofOrdering::LEXICOGRAPHIC
Lexicographic ordering for tensor-product FiniteElements.

alpha
const double alpha
Definition: ex15.cpp:369

mfem::Vector
Vector data type.
Definition: vector.hpp:60

mfem::Vector::MakeRef
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:577

mfem::VectorGridFunctionCoefficient
Vector coefficient defined by a vector GridFunction.
Definition: coefficient.hpp:484

gridfunc.hpp

mfem::u
double u(const Vector &xvec)
Definition: lor_mms.hpp:24

mfem::Operator
Abstract operator.
Definition: operator.hpp:24

mfem::ConvectionIntegrator::nq
int nq
Definition: bilininteg.hpp:2138

mfem::ConvectionIntegrator::quad1D
int quad1D
Definition: bilininteg.hpp:2138

mfem::QuadratureFunction
Class representing a function through its values (scalar or vector) at quadrature points...
Definition: gridfunc.hpp:734

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:426

mfem::QVectorLayout::byVDIM
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)

mfem::ConvectionIntegrator::maps
const DofToQuad * maps
Not owned.
Definition: bilininteg.hpp:2136

mfem::NonlinearFormIntegrator::pa_mt
MemoryType pa_mt
Definition: nonlininteg.hpp:35

mfem::ConvectionIntegrator::Q
VectorCoefficient * Q
Definition: bilininteg.hpp:2132