Initial commit for Exaconstit using libCEED through MFEM.

system, install libceed q-function headers.

Add a runtime check for q-function header location: install/source
directories.

CHANGELOG

@@ -20,6 +20,17 @@ GPU support
   backends is: "occa-cuda", "raja-cuda", "cuda", "hip", "occa-omp", "raja-omp",
   "omp", "occa-cpu", "raja-cpu", and "cpu".
 
+libCEED support
+---------------
+- Added support for libCEED, the portable library for high-order operator
+  evaluation developed by the Center for Efficient Exascale Discretizations in
+  the Exascale Computing Project, https://github.com/CEED/libCEED. This initial
+  integration includes Mass and Diffusion integrators. libCEED GPU backends can
+  be used without specific MFEM configuration, however it is highly recommended
+  to use the "cuda" build option to minimize memory transfers. Both CPU and GPU
+  modes are available as MFEM device backends (ceed-cpu and ceed-cuda), using
+  some of the best performing CPU and GPU backends from libCEED.
+
 Miscellaneous
 -------------
 - Upgraded the SUNDIALS interface to utilize SUNDIALS version 5.0.  This 

CMakeLists.txt

@@ -248,6 +248,10 @@ if (MFEM_USE_MPFR)
   find_package(MPFR REQUIRED)
 endif()
 
+if (MFEM_USE_CEED)
+  find_package(libCEED REQUIRED)
+endif()
+
 if (MFEM_USE_CONDUIT)
     find_package(Conduit REQUIRED conduit relay blueprint )
 endif()
@@ -323,7 +327,7 @@ endif()
 #    be before SuiteSparse.
 set(MFEM_TPLS MPI_CXX OPENMP BLAS LAPACK METIS HYPRE SuiteSparse SUNDIALS PETSC
     MESQUITE SuperLUDist STRUMPACK AXOM CONDUIT GECKO GNUTLS NETCDF MPFR PUMI
-    POSIXCLOCKS MFEMBacktrace ZLIB OCCA RAJA)
+    POSIXCLOCKS MFEMBacktrace ZLIB OCCA RAJA CEED)
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")
@@ -535,6 +539,13 @@ if (MFEM_USE_OCCA)
     FILES_MATCHING PATTERN "*.okl")
 endif()
 
+# Install the libCEED files
+if (MFEM_USE_CEED)
+  install(DIRECTORY ${MFEM_SOURCE_DIRS}
+    DESTINATION ${INSTALL_INCLUDE_DIR}/mfem
+    FILES_MATCHING PATTERN "fem/libceed/*.h")
+endif()
+
 # Install ${HEADERS}
 # ---
 # foreach (HDR ${HEADERS})

INSTALL

@@ -442,6 +442,11 @@ MFEM_USE_OCCA = YES/NO
    backends. In order to use the OCCA CUDA backend, CUDA support must be enabled
    in MFEM as well, i.e. MFEM_USE_CUDA=YES must be set.
 
+MFEM_USE_CEED = YES/NO
+   Enables support for the libCEED library in MFEM. libCEED is a portable
+   library for performant high-order operator evaluation developed by the Center
+   for Efficient Exascale Discretizations in the Exascale Computing Project.
+
 MFEM_BUILD_TAG = (any value)
    An optional tag to characterize the build. Exported to config/config.mk.
    Can be used to identify the MFEM build from other makefiles.
@@ -570,6 +575,11 @@ The specific libraries and their options are:
   URL: https://libocca.org
   Options: OCCA_DIR, OCCA_OPT, OCCA_LIB.
 
+- CEED, used when MFEM_USE_CEED = YES.
+  URL: https://github.com/CEED/libCEED
+       https://ceed.exascaleproject.org/libceed
+  Options: CEED_DIR, CEED_OPT, CEED_LIB
+
 - RAJA, used when MFEM_USE_RAJA = YES.
   URL: https://github.com/LLNL/RAJA
   Options: RAJA_DIR, RAJA_OPT, RAJA_LIB.
@@ -712,6 +722,7 @@ MFEM_USE_GZSTREAM
 MFEM_USE_PUMI
 MFEM_USE_CUDA
 MFEM_USE_OCCA
+MFEM_USE_CEED
 MFEM_USE_RAJA
 
 The following options are CMake specific:

config/cmake/MFEMConfig.cmake.in

@@ -44,6 +44,7 @@ set(MFEM_USE_PUMI @MFEM_USE_PUMI@)
 set(MFEM_USE_CUDA @MFEM_USE_CUDA@)
 set(MFEM_USE_OCCA @MFEM_USE_OCCA@)
 set(MFEM_USE_RAJA @MFEM_USE_RAJA@)
+set(MFEM_USE_CEED @MFEM_USE_CEED@)
 
 set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
 set(MFEM_CXX_FLAGS "@CMAKE_CXX_FLAGS@")

config/cmake/config.hpp.in

@@ -120,6 +120,9 @@
 // Enable MFEM functionality based on the OCCA library
 #cmakedefine MFEM_USE_OCCA
 
+// Enable MFEM functionality based on the libCEED library
+#cmakedefine MFEM_USE_CEED
+
 // Which library functions to use in class StopWatch for measuring time.
 // For a list of the available options, see INSTALL.
 // If not defined, an option is selected automatically.

config/cmake/modules/FindlibCEED.cmake

@@ -0,0 +1,19 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Defines the following variables:
+#   - CEED_FOUND
+#   - CEED_LIBRARIES
+#   - CEED_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(libCEED CEED CEED_DIR "include" ceed.h "lib" ceed
+  "Paths to headers required by libCEED." "Libraries required by libCEED.")

config/config.hpp.in

@@ -135,6 +135,9 @@
 // Enable functionality based on the OCCA library.
 // #define MFEM_USE_OCCA
 
+// Enable functionality based on the CEED library.
+// #define MFEM_USE_CEED
+
 // Version of HYPRE used for building MFEM.
 // #define MFEM_HYPRE_VERSION @MFEM_HYPRE_VERSION@
 

config/config.mk.in

@@ -45,6 +45,7 @@ MFEM_USE_CUDA          = @MFEM_USE_CUDA@
 MFEM_USE_HIP           = @MFEM_USE_HIP@
 MFEM_USE_RAJA          = @MFEM_USE_RAJA@
 MFEM_USE_OCCA          = @MFEM_USE_OCCA@
+MFEM_USE_CEED          = @MFEM_USE_CEED@
 
 # Compiler, compile options, and link options
 MFEM_CXX       = @MFEM_CXX@

config/defaults.cmake

@@ -45,6 +45,7 @@ option(MFEM_USE_PUMI "Enable PUMI" OFF)
 option(MFEM_USE_CUDA "Enable CUDA" OFF)
 option(MFEM_USE_OCCA "Enable OCCA" OFF)
 option(MFEM_USE_RAJA "Enable RAJA" OFF)
+option(MFEM_USE_CEED "Enable CEED" OFF)
 
 set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 
@@ -162,6 +163,7 @@ set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
 
 set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
 set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
+set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")
 
 set(BLAS_INCLUDE_DIRS "" CACHE STRING "Path to BLAS headers.")
 set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")

config/defaults.mk

@@ -133,6 +133,7 @@ MFEM_USE_CUDA          = NO
 MFEM_USE_HIP           = NO
 MFEM_USE_RAJA          = NO
 MFEM_USE_OCCA          = NO
+MFEM_USE_CEED          = NO
 
 # Compile and link options for zlib.
 ZLIB_DIR =
@@ -321,6 +322,11 @@ OCCA_DIR = @MFEM_DIR@/../occa
 OCCA_OPT = -I$(OCCA_DIR)/include
 OCCA_LIB = $(XLINKER)-rpath,$(OCCA_DIR)/lib -L$(OCCA_DIR)/lib -locca
 
+# CEED library configuration
+CEED_DIR ?= @MFEM_DIR@/../libCEED
+CEED_OPT = -I$(CEED_DIR)/include
+CEED_LIB = $(XLINKER)-rpath,$(CEED_DIR)/lib -L$(CEED_DIR)/lib -lceed
+
 # RAJA library configuration
 RAJA_DIR = @MFEM_DIR@/../raja
 RAJA_OPT = -I$(RAJA_DIR)/include

fem/CMakeLists.txt

@@ -93,6 +93,17 @@ if (MFEM_USE_MPI)
     pnonlinearform.hpp)
 endif()
 
+if (MFEM_USE_CEED)
+   list(APPEND SRCS
+     libceed/ceed.cpp
+     libceed/diffusion.cpp
+     libceed/mass.cpp)
+   list(APPEND HDRS
+     libceed/ceed.hpp
+     libceed/diffusion.hpp
+     libceed/mass.hpp)
+endif()
+
 convert_filenames_to_full_paths(SRCS)
 convert_filenames_to_full_paths(HDRS)
 

fem/bilinearform_ext.cpp

@@ -102,25 +102,35 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
    Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
 
    const int iSz = integrators.Size();
-   if (elem_restrict_lex)
+   if (Device::Allows(Backend::CEED_MASK))
    {
-      elem_restrict_lex->Mult(x, localX);
-      localY = 0.0;
-      for (int i = 0; i < iSz; ++i)
-      {
-         integrators[i]->AddMultPA(localX, localY);
-      }
-      elem_restrict_lex->MultTranspose(localY, y);
-   }
-   else
-   {
-      y.UseDevice(true); // typically this is a large vector, so store on device
-      y = 0.0;
       for (int i = 0; i < iSz; ++i)
       {
          integrators[i]->AddMultPA(x, y);
       }
    }
+   else
+   {
+      if (elem_restrict_lex)
+      {
+         elem_restrict_lex->Mult(x, localX);
+         localY = 0.0;
+         for (int i = 0; i < iSz; ++i)
+         {
+            integrators[i]->AddMultPA(localX, localY);
+         }
+         elem_restrict_lex->MultTranspose(localY, y);
+      }
+      else
+      {
+         y.UseDevice(true); // typically this is a large vector, so store on device
+         y = 0.0;
+         for (int i = 0; i < iSz; ++i)
+         {
+            integrators[i]->AddMultPA(x, y);
+         }
+      }
+   }
 }
 
 void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const

fem/bilininteg.hpp

@@ -19,6 +19,10 @@
 namespace mfem
 {
 
+#ifdef MFEM_USE_CEED
+struct CeedData;
+#endif
+
 /// Abstract base class BilinearFormIntegrator
 class BilinearFormIntegrator : public NonlinearFormIntegrator
 {
@@ -1679,6 +1683,11 @@ private:
    int dim, ne, dofs1D, quad1D;
    Vector pa_data;
 
+#ifdef MFEM_USE_CEED
+   // CEED extension
+   CeedData* ceedDataPtr;
+#endif
+
 public:
    /// Construct a diffusion integrator with coefficient Q = 1
    DiffusionIntegrator() { Q = NULL; MQ = NULL; maps = NULL; geom = NULL; }
@@ -1725,6 +1734,79 @@ public:
                                          const FiniteElement &test_fe);
 };
 
+/** Class for integrating the bilinear form a(u,v) := (Q grad u, grad v) where Q
+    can be a scalar or a matrix coefficient. */
+class MechanicsIntegrator: public BilinearFormIntegrator
+{
+protected:
+   Coefficient *Q;
+   MatrixCoefficient *MQ;
+
+private:
+   Vector vec, pointflux, shape;
+#ifndef MFEM_THREAD_SAFE
+   DenseMatrix dshape, dshapedxt, invdfdx, mq;
+   DenseMatrix te_dshape, te_dshapedxt;
+#endif
+
+   // PA extension
+   const DofToQuad *maps;         ///< Not owned
+   const GeometricFactors *geom;  ///< Not owned
+   int dim, ne, dofs1D, quad1D;
+   Vector pa_data;
+   
+#ifdef MFEM_USE_CEED
+   // CEED extension
+   CeedData* ceedDataPtr;
+   double* ktan_ptr;
+#endif
+
+public:
+   /// Construct a diffusion integrator with coefficient Q = 1
+   MechanicsIntegrator() { Q = NULL; MQ = NULL; maps = NULL; geom = NULL; }
+
+   /// Construct a diffusion integrator with a scalar coefficient q
+   MechanicsIntegrator(Coefficient &q)
+      : Q(&q) { MQ = NULL; maps = NULL; geom = NULL; }
+
+   /// Construct a diffusion integrator with a matrix coefficient q
+   MechanicsIntegrator(MatrixCoefficient &q)
+      : MQ(&q) { Q = NULL; maps = NULL; geom = NULL; }
+
+   /** Given a particular Finite Element
+       computes the element stiffness matrix elmat. */
+   virtual void AssembleElementMatrix(const FiniteElement &el,
+                                      ElementTransformation &Trans,
+                                      DenseMatrix &elmat);
+   /** Given a trial and test Finite Element computes the element stiffness
+       matrix elmat. */
+   virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                                       const FiniteElement &test_fe,
+                                       ElementTransformation &Trans,
+                                       DenseMatrix &elmat);
+
+   /// Perform the local action of the BilinearFormIntegrator
+   virtual void AssembleElementVector(const FiniteElement &el,
+                                      ElementTransformation &Tr,
+                                      const Vector &elfun, Vector &elvect);
+
+   virtual void ComputeElementFlux(const FiniteElement &el,
+                                   ElementTransformation &Trans,
+                                   Vector &u, const FiniteElement &fluxelem,
+                                   Vector &flux, int with_coef = 1);
+
+   virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
+                                    ElementTransformation &Trans,
+                                    Vector &flux, Vector *d_energy = NULL);
+
+   virtual void AssemblePA(const FiniteElementSpace&);
+
+   virtual void AddMultPA(const Vector&, Vector&) const;
+
+   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
+                                         const FiniteElement &test_fe);
+};
+
 /** Class for local mass matrix assembling a(u,v) := (Q u, v) */
 class MassIntegrator: public BilinearFormIntegrator
 {
@@ -1738,6 +1820,11 @@ protected:
    const DofToQuad *maps;         ///< Not owned
    const GeometricFactors *geom;  ///< Not owned
    int dim, ne, nq, dofs1D, quad1D;
+   
+#ifdef MFEM_USE_CEED
+   // CEED extension
+   CeedData* ceedDataPtr;
+#endif
 
 public:
    MassIntegrator(const IntegrationRule *ir = NULL)

fem/bilininteg_diffusion.cpp

@@ -12,6 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
+#include "libceed/diffusion.hpp"
 
 using namespace std;
 
@@ -191,25 +192,37 @@ static void PADiffusionSetup(const int dim,
 
 void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
-   // Assumes tensor-product elements
-   Mesh *mesh = fes.GetMesh();
    const FiniteElement &el = *fes.GetFE(0);
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el);
-   const int dims = el.GetDim();
-   const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
-   const int nq = ir->GetNPoints();
-   dim = mesh->Dimension();
-   ne = fes.GetNE();
-   geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
-   maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
-   dofs1D = maps->ndof;
-   quad1D = maps->nqpt;
-   pa_data.SetSize(symmDims * nq * ne, Device::GetMemoryType());
-   ConstantCoefficient *cQ = dynamic_cast<ConstantCoefficient*>(Q);
-   MFEM_VERIFY(cQ != NULL, "only ConstantCoefficient is supported!");
-   const double coeff = cQ->constant;
-   PADiffusionSetup(dim, dofs1D, quad1D, ne, ir->GetWeights(), geom->J,
-                    coeff, pa_data);
+#ifdef MFEM_USE_CEED
+   if (Device::Allows(Backend::CEED_MASK))
+   {
+      CeedData* ptr = new CeedData();
+      ceedDataPtr = ptr;
+      InitCeedCoeff(Q, ptr);
+      CeedPADiffusionAssemble(fes, *ir, *ptr);
+   }
+   else
+#endif
+   {
+      // Assumes tensor-product elements
+      Mesh *mesh = fes.GetMesh();
+      const int dims = el.GetDim();
+      const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
+      const int nq = ir->GetNPoints();
+      dim = mesh->Dimension();
+      ne = fes.GetNE();
+      geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
+      maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
+      dofs1D = maps->ndof;
+      quad1D = maps->nqpt;
+      pa_data.SetSize(symmDims * nq * ne, Device::GetMemoryType());
+      ConstantCoefficient *cQ = dynamic_cast<ConstantCoefficient*>(Q);
+      MFEM_VERIFY(cQ != NULL, "only ConstantCoefficient is supported!");
+      const double coeff = cQ->constant;
+      PADiffusionSetup(dim, dofs1D, quad1D, ne, ir->GetWeights(), geom->J,
+                       coeff, pa_data);
+   }
 }
 
 #ifdef MFEM_USE_OCCA
@@ -1088,9 +1101,40 @@ static void PADiffusionApply(const int dim,
 // PA Diffusion Apply kernel
 void DiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
-   PADiffusionApply(dim, dofs1D, quad1D, ne,
-                    maps->B, maps->G, maps->Bt, maps->Gt,
-                    pa_data, x, y);
+#ifdef MFEM_USE_CEED
+   if (Device::Allows(Backend::CEED_MASK))
+   {
+      const CeedScalar *x_ptr;
+      CeedScalar *y_ptr;
+      CeedMemType mem;
+      CeedGetPreferredMemType(internal::ceed, &mem);
+      if ( Device::Allows(Backend::CUDA) && mem==CEED_MEM_DEVICE )
+      {
+         x_ptr = x.Read();
+         y_ptr = y.ReadWrite();
+      }
+      else
+      {
+         x_ptr = x.HostRead();
+         y_ptr = y.HostReadWrite();
+         mem = CEED_MEM_HOST;
+      }
+      CeedVectorSetArray(ceedDataPtr->u, mem, CEED_USE_POINTER,
+                         const_cast<CeedScalar*>(x_ptr));
+      CeedVectorSetArray(ceedDataPtr->v, mem, CEED_USE_POINTER, y_ptr);
+
+      CeedOperatorApply(ceedDataPtr->oper, ceedDataPtr->u, ceedDataPtr->v,
+                        CEED_REQUEST_IMMEDIATE);
+
+      CeedVectorSyncArray(ceedDataPtr->v, mem);
+   }
+   else
+#endif
+   {
+      PADiffusionApply(dim, dofs1D, quad1D, ne,
+                       maps->B, maps->G, maps->Bt, maps->Gt,
+                       pa_data, x, y);
+   }
 }
 
 } // namespace mfem

fem/bilininteg_mass.cpp

@@ -12,6 +12,7 @@
 #include "../general/forall.hpp"
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
+#include "libceed/mass.hpp"
 
 using namespace std;
 
@@ -29,76 +30,88 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    const FiniteElement &el = *fes.GetFE(0);
    ElementTransformation *T = mesh->GetElementTransformation(0);
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, *T);
-   dim = mesh->Dimension();
-   ne = fes.GetMesh()->GetNE();
-   nq = ir->GetNPoints();
-   geom = mesh->GetGeometricFactors(*ir, GeometricFactors::COORDINATES |
-                                    GeometricFactors::JACOBIANS);
-   maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
-   dofs1D = maps->ndof;
-   quad1D = maps->nqpt;
-   pa_data.SetSize(ne*nq, Device::GetMemoryType());
-   ConstantCoefficient *const_coeff = dynamic_cast<ConstantCoefficient*>(Q);
-   // TODO: other types of coefficients ...
-   if (dim==1) { MFEM_ABORT("Not supported yet... stay tuned!"); }
-   if (dim==2)
+#ifdef MFEM_USE_CEED
+   if (Device::Allows(Backend::CEED_MASK))
    {
-      double constant = 0.0;
-      if (const_coeff)
-      {
-         constant = const_coeff->constant;
-      }
-      else
-      {
-         MFEM_ABORT("Coefficient type not supported");
-      }
-      const int NE = ne;
-      const int NQ = nq;
-      auto w = ir->GetWeights().Read();
-      auto J = Reshape(geom->J.Read(), NQ,2,2,NE);
-      auto v = Reshape(pa_data.Write(), NQ, NE);
-      MFEM_FORALL(e, NE,
-      {
-         for (int q = 0; q < NQ; ++q)
-         {
-            const double J11 = J(q,0,0,e);
-            const double J12 = J(q,1,0,e);
-            const double J21 = J(q,0,1,e);
-            const double J22 = J(q,1,1,e);
-            const double detJ = (J11*J22)-(J21*J12);
-            v(q,e) =  w[q] * constant * detJ;
-         }
-      });
+      CeedData* ptr = new CeedData();
+      ceedDataPtr = ptr;
+      InitCeedCoeff(Q, ptr);
+      CeedPAMassAssemble(fes, *ir, *ptr);
    }
-   if (dim==3)
+   else
+#endif
    {
-      double constant = 0.0;
-      if (const_coeff)
+      dim = mesh->Dimension();
+      ne = fes.GetMesh()->GetNE();
+      nq = ir->GetNPoints();
+      geom = mesh->GetGeometricFactors(*ir, GeometricFactors::COORDINATES |
+                                       GeometricFactors::JACOBIANS);
+      maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
+      dofs1D = maps->ndof;
+      quad1D = maps->nqpt;
+      pa_data.SetSize(ne*nq, Device::GetMemoryType());
+      ConstantCoefficient *const_coeff = dynamic_cast<ConstantCoefficient*>(Q);
+      // TODO: other types of coefficients ...
+      if (dim==1) { MFEM_ABORT("Not supported yet... stay tuned!"); }
+      if (dim==2)
       {
-         constant = const_coeff->constant;
-      }
-      else
-      {
-         MFEM_ABORT("Coefficient type not supported");
-      }
-      const int NE = ne;
-      const int NQ = nq;
-      auto W = ir->GetWeights().Read();
-      auto J = Reshape(geom->J.Read(), NQ,3,3,NE);
-      auto v = Reshape(pa_data.Write(), NQ,NE);
-      MFEM_FORALL(e, NE,
-      {
-         for (int q = 0; q < NQ; ++q)
+         double constant = 0.0;
+         if (const_coeff)
          {
-            const double J11 = J(q,0,0,e), J12 = J(q,0,1,e), J13 = J(q,0,2,e);
-            const double J21 = J(q,1,0,e), J22 = J(q,1,1,e), J23 = J(q,1,2,e);
-            const double J31 = J(q,2,0,e), J32 = J(q,2,1,e), J33 = J(q,2,2,e);
-            const double detJ = J11 * (J22 * J33 - J32 * J23) -
-            /* */               J21 * (J12 * J33 - J32 * J13) +
-            /* */               J31 * (J12 * J23 - J22 * J13);
-            v(q,e) = W[q] * constant * detJ;
+            constant = const_coeff->constant;
          }
-      });
+         else
+         {
+            MFEM_ABORT("Coefficient type not supported");
+         }
+         const int NE = ne;
+         const int NQ = nq;
+         auto w = ir->GetWeights().Read();
+         auto J = Reshape(geom->J.Read(), NQ,2,2,NE);
+         auto v = Reshape(pa_data.Write(), NQ, NE);
+         MFEM_FORALL(e, NE,
+         {
+            for (int q = 0; q < NQ; ++q)
+            {
+               const double J11 = J(q,0,0,e);
+               const double J12 = J(q,1,0,e);
+               const double J21 = J(q,0,1,e);
+               const double J22 = J(q,1,1,e);
+               const double detJ = (J11*J22)-(J21*J12);
+               v(q,e) =  w[q] * constant * detJ;
+            }
+         });
+      }
+      if (dim==3)
+      {
+         double constant = 0.0;
+         if (const_coeff)
+         {
+            constant = const_coeff->constant;
+         }
+         else
+         {
+            MFEM_ABORT("Coefficient type not supported");
+         }
+         const int NE = ne;
+         const int NQ = nq;
+         auto W = ir->GetWeights().Read();
+         auto J = Reshape(geom->J.Read(), NQ,3,3,NE);
+         auto v = Reshape(pa_data.Write(), NQ,NE);
+         MFEM_FORALL(e, NE,
+         {
+            for (int q = 0; q < NQ; ++q)
+            {
+               const double J11 = J(q,0,0,e), J12 = J(q,0,1,e), J13 = J(q,0,2,e);
+               const double J21 = J(q,1,0,e), J22 = J(q,1,1,e), J23 = J(q,1,2,e);
+               const double J31 = J(q,2,0,e), J32 = J(q,2,1,e), J33 = J(q,2,2,e);
+               const double detJ = J11 * (J22 * J33 - J32 * J23) -
+               /* */               J21 * (J12 * J33 - J32 * J13) +
+               /* */               J31 * (J12 * J23 - J22 * J13);
+               v(q,e) = W[q] * constant * detJ;
+            }
+         });
+      }
    }
 }
 
@@ -783,7 +796,38 @@ static void PAMassApply(const int dim,
 
 void MassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
-   PAMassApply(dim, dofs1D, quad1D, ne, maps->B, maps->Bt, pa_data, x, y);
+#ifdef MFEM_USE_CEED
+   if (Device::Allows(Backend::CEED_MASK))
+   {
+      const CeedScalar *x_ptr;
+      CeedScalar *y_ptr;
+      CeedMemType mem;
+      CeedGetPreferredMemType(internal::ceed, &mem);
+      if ( Device::Allows(Backend::CUDA) && mem==CEED_MEM_DEVICE )
+      {
+         x_ptr = x.Read();
+         y_ptr = y.ReadWrite();
+      }
+      else
+      {
+         x_ptr = x.HostRead();
+         y_ptr = y.HostReadWrite();
+         mem = CEED_MEM_HOST;
+      }
+      CeedVectorSetArray(ceedDataPtr->u, mem, CEED_USE_POINTER,
+                         const_cast<CeedScalar*>(x_ptr));
+      CeedVectorSetArray(ceedDataPtr->v, mem, CEED_USE_POINTER, y_ptr);
+
+      CeedOperatorApply(ceedDataPtr->oper, ceedDataPtr->u, ceedDataPtr->v,
+                        CEED_REQUEST_IMMEDIATE);
+
+      CeedVectorSyncArray(ceedDataPtr->v, mem);
+   }
+   else
+#endif
+   {
+      PAMassApply(dim, dofs1D, quad1D, ne, maps->B, maps->Bt, pa_data, x, y);
+   }
 }
 
 } // namespace mfem

fem/libceed/ceed.cpp

@@ -0,0 +1,170 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "ceed.hpp"
+#include "../../general/device.hpp"
+
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#ifndef _WIN32
+typedef struct stat struct_stat;
+#else
+#define stat(dir, buf) _stat(dir, buf)
+#define S_ISDIR(mode) _S_IFDIR(mode)
+typedef struct _stat struct_stat;
+#endif
+
+namespace mfem
+{
+
+namespace internal
+{
+
+extern Ceed ceed;
+
+std::string ceed_path;
+
+}
+
+#ifdef MFEM_USE_CEED
+void InitCeedCoeff(Coefficient* Q, CeedData* ptr)
+{
+   if (ConstantCoefficient* coeff = dynamic_cast<ConstantCoefficient*>(Q))
+   {
+      CeedConstCoeff* ceedCoeff = new CeedConstCoeff{coeff->constant};
+      ptr->coeff_type = CeedCoeff::Const;
+      ptr->coeff = (void*)ceedCoeff;
+   }
+   else if (GridFunctionCoefficient* coeff =
+               dynamic_cast<GridFunctionCoefficient*>(Q))
+   {
+      CeedGridCoeff* ceedCoeff = new CeedGridCoeff;
+      ceedCoeff->coeff = coeff->GetGridFunction();
+      ptr->coeff_type = CeedCoeff::Grid;
+      ptr->coeff = (void*)ceedCoeff;
+   }
+   else
+   {
+      MFEM_ABORT("This type of Coefficient is not supported.");
+   }
+}
+
+void InitCeedTensorBasisAndRestriction(const mfem::FiniteElementSpace &fes,
+                                       const mfem::IntegrationRule &ir,
+                                       Ceed ceed, CeedBasis *basis,
+                                       CeedElemRestriction *restr)
+{
+   mfem::Mesh *mesh = fes.GetMesh();
+   const mfem::FiniteElement *fe = fes.GetFE(0);
+   const int order = fes.GetOrder(0);
+   mfem::Array<int> dof_map;
+   switch (mesh->Dimension())
+   {
+      case 1:
+      {
+         const mfem::H1_SegmentElement *h1_fe =
+            dynamic_cast<const mfem::H1_SegmentElement *>(fe);
+         MFEM_VERIFY(h1_fe, "invalid FE");
+         h1_fe->GetDofMap().Copy(dof_map);
+         break;
+      }
+      case 2:
+      {
+         const mfem::H1_QuadrilateralElement *h1_fe =
+            dynamic_cast<const mfem::H1_QuadrilateralElement *>(fe);
+         MFEM_VERIFY(h1_fe, "invalid FE");
+         h1_fe->GetDofMap().Copy(dof_map);
+         break;
+      }
+      case 3:
+      {
+         const mfem::H1_HexahedronElement *h1_fe =
+            dynamic_cast<const mfem::H1_HexahedronElement *>(fe);
+         MFEM_VERIFY(h1_fe, "invalid FE");
+         h1_fe->GetDofMap().Copy(dof_map);
+         break;
+      }
+   }
+   const mfem::FiniteElement *fe1d =
+      fes.FEColl()->FiniteElementForGeometry(mfem::Geometry::SEGMENT);
+   mfem::DenseMatrix shape1d(fe1d->GetDof(), ir.GetNPoints());
+   mfem::DenseMatrix grad1d(fe1d->GetDof(), ir.GetNPoints());
+   mfem::Vector qref1d(ir.GetNPoints()), qweight1d(ir.GetNPoints());
+   mfem::Vector shape_i(shape1d.Height());
+   mfem::DenseMatrix grad_i(grad1d.Height(), 1);
+   const mfem::H1_SegmentElement *h1_fe1d =
+      dynamic_cast<const mfem::H1_SegmentElement *>(fe1d);
+   MFEM_VERIFY(h1_fe1d, "invalid FE");
+   const mfem::Array<int> &dof_map_1d = h1_fe1d->GetDofMap();
+   for (int i = 0; i < ir.GetNPoints(); i++)
+   {
+      const mfem::IntegrationPoint &ip = ir.IntPoint(i);
+      qref1d(i) = ip.x;
+      qweight1d(i) = ip.weight;
+      fe1d->CalcShape(ip, shape_i);
+      fe1d->CalcDShape(ip, grad_i);
+      for (int j = 0; j < shape1d.Height(); j++)
+      {
+         shape1d(j, i) = shape_i(dof_map_1d[j]);
+         grad1d(j, i) = grad_i(dof_map_1d[j], 0);
+      }
+   }
+   CeedBasisCreateTensorH1(ceed, mesh->Dimension(), fes.GetVDim(), order + 1,
+                           ir.GetNPoints(), shape1d.GetData(),
+                           grad1d.GetData(), qref1d.GetData(),
+                           qweight1d.GetData(), basis);
+
+   const mfem::Table &el_dof = fes.GetElementToDofTable();
+   mfem::Array<int> tp_el_dof(el_dof.Size_of_connections());
+   for (int i = 0; i < mesh->GetNE(); i++)
+   {
+      const int el_offset = fe->GetDof() * i;
+      for (int j = 0; j < fe->GetDof(); j++)
+      {
+         tp_el_dof[j + el_offset] = el_dof.GetJ()[dof_map[j] + el_offset];
+      }
+   }
+   CeedElemRestrictionCreate(ceed, mesh->GetNE(), fe->GetDof(),
+                             fes.GetNDofs(), fes.GetVDim(), CEED_MEM_HOST, CEED_COPY_VALUES,
+                             tp_el_dof.GetData(), restr);
+}
+
+const std::string &GetCeedPath()
+{
+   if (internal::ceed_path.empty())
+   {
+      const char *install_dir = MFEM_INSTALL_DIR "/include/mfem/fem/libceed";
+      const char *source_dir = MFEM_SOURCE_DIR "/fem/libceed";
+      struct_stat m_stat;
+      if (stat(install_dir, &m_stat) == 0 && S_ISDIR(m_stat.st_mode))
+      {
+         internal::ceed_path = install_dir;
+      }
+      else if (stat(source_dir, &m_stat) == 0 && S_ISDIR(m_stat.st_mode))
+      {
+         internal::ceed_path = source_dir;
+      }
+      else
+      {
+         MFEM_ABORT("Cannot find libCEED kernels in MFEM_INSTALL_DIR or "
+                    "MFEM_SOURCE_DIR");
+      }
+   }
+#ifdef MFEM_DEBUG
+   mfem::out << "Using libCEED dir: " << internal::ceed_path << std::endl;
+#endif
+   return internal::ceed_path;
+}
+
+#endif
+
+} // namespace mfem

fem/libceed/ceed.hpp

@@ -0,0 +1,81 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_LIBCEED_HPP
+#define MFEM_LIBCEED_HPP
+
+#include "../gridfunc.hpp"
+#include "../fespace.hpp"
+
+#ifdef MFEM_USE_CEED
+#include <ceed.h>
+#else
+typedef void* Ceed;
+typedef int CeedInt;
+typedef double CeedScalar;
+#define CEED_QFUNCTION(name) int name
+#endif
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+namespace internal { extern Ceed ceed; }
+
+/// A structure used to pass additional data to f_build_diff and f_apply_diff
+struct BuildContext { CeedInt dim, space_dim; CeedScalar coeff; };
+
+enum class CeedCoeff { Const, Grid };
+
+struct CeedConstCoeff
+{
+   double val;
+};
+
+struct CeedGridCoeff
+{
+   GridFunction* coeff;
+   CeedBasis basis;
+   CeedElemRestriction restr;
+   CeedVector coeffVector;
+};
+
+struct CeedData
+{
+   CeedOperator build_oper, oper;
+   CeedBasis basis, mesh_basis;
+   CeedElemRestriction restr, mesh_restr, restr_i, mesh_restr_i;
+   CeedQFunction apply_qfunc, build_qfunc;
+   CeedVector node_coords, rho;
+   CeedCoeff coeff_type;
+   void* coeff;
+   BuildContext build_ctx;
+
+   CeedVector u, v;
+};
+
+/// Identifies the type of coefficient of the Integrator to initialize accordingly the CeedData
+void InitCeedCoeff(Coefficient* Q, CeedData* ptr);
+
+/// Initialize a tensor CeedBasis and a CeedElemRestriction
+void InitCeedTensorBasisAndRestriction(const mfem::FiniteElementSpace &fes,
+                                       const mfem::IntegrationRule &ir,
+                                       Ceed ceed, CeedBasis *basis,
+                                       CeedElemRestriction *restr);
+
+const std::string &GetCeedPath();
+
+#endif
+
+}
+
+#endif // MFEM_LIBCEED_HPP

fem/libceed/diffusion.cpp

@@ -0,0 +1,145 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "ceed.hpp"
+#include "../../general/device.hpp"
+
+#include "diffusion.h"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+void CeedPADiffusionAssemble(const FiniteElementSpace &fes,
+                             const mfem::IntegrationRule &irm, CeedData& ceedData)
+{
+   Ceed ceed(internal::ceed);
+   mfem::Mesh *mesh = fes.GetMesh();
+   const int ir_order = irm.GetOrder();
+   const mfem::IntegrationRule &ir =
+      mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order);
+   CeedInt nqpts, nelem = mesh->GetNE(), dim = mesh->SpaceDimension();
+   mesh->EnsureNodes();
+   InitCeedTensorBasisAndRestriction(fes, ir, ceed, &ceedData.basis, &ceedData.restr);
+
+   const mfem::FiniteElementSpace *mesh_fes = mesh->GetNodalFESpace();
+   MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
+   InitCeedTensorBasisAndRestriction(*mesh_fes, ir, ceed, &ceedData.mesh_basis, &ceedData.mesh_restr);
+   CeedBasisGetNumQuadraturePoints(ceedData.basis, &nqpts);
+
+   CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts,
+                                     nqpts * nelem, dim * (dim + 1) / 2, &ceedData.restr_i);
+   CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts,
+                                     nqpts * nelem, 1, &ceedData.mesh_restr_i);
+
+   CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &ceedData.node_coords);
+   CeedVectorSetArray(ceedData.node_coords, CEED_MEM_HOST, CEED_USE_POINTER,
+                      mesh->GetNodes()->GetData());
+
+   CeedVectorCreate(ceed, nelem * nqpts * dim * (dim + 1) / 2, &ceedData.rho);
+
+   // Context data to be passed to the 'f_build_diff' Q-function.
+   ceedData.build_ctx.dim = mesh->Dimension();
+   ceedData.build_ctx.space_dim = mesh->SpaceDimension();
+
+   std::string diff_qf_file = GetCeedPath() + "/diffusion.h";
+   std::string diff_qf;
+
+   // Create the Q-function that builds the diff operator (i.e. computes its
+   // quadrature data) and set its context data.
+   switch (ceedData.coeff_type)
+   {
+      case CeedCoeff::Const:
+         diff_qf = diff_qf_file + ":f_build_diff_const";
+         CeedQFunctionCreateInterior(ceed, 1, f_build_diff_const,
+                                     diff_qf.c_str(),
+                                     &ceedData.build_qfunc);
+         ceedData.build_ctx.coeff = ((CeedConstCoeff*)ceedData.coeff)->val;
+         break;
+      case CeedCoeff::Grid:
+         diff_qf = diff_qf_file + ":f_build_diff_grid";
+         CeedQFunctionCreateInterior(ceed, 1, f_build_diff_grid,
+                                     diff_qf.c_str(),
+                                     &ceedData.build_qfunc);
+         CeedQFunctionAddInput(ceedData.build_qfunc, "coeff", 1, CEED_EVAL_INTERP);
+         break;
+      default:
+         MFEM_ABORT("This coeff_type is not handled");
+   }
+   CeedQFunctionAddInput(ceedData.build_qfunc, "dx", dim * dim, CEED_EVAL_GRAD);
+   CeedQFunctionAddInput(ceedData.build_qfunc, "weights", 1, CEED_EVAL_WEIGHT);
+   CeedQFunctionAddOutput(ceedData.build_qfunc, "rho", dim * (dim + 1) / 2,
+                          CEED_EVAL_NONE);
+   CeedQFunctionSetContext(ceedData.build_qfunc, &ceedData.build_ctx,
+                           sizeof(ceedData.build_ctx));
+
+   // Create the operator that builds the quadrature data for the diff operator.
+   CeedOperatorCreate(ceed, ceedData.build_qfunc, NULL, NULL,
+                      &ceedData.build_oper);
+   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
+   if (mesh_fes->GetOrdering()==Ordering::byVDIM)
+   {
+      lmode = CEED_TRANSPOSE;
+   }
+   if (ceedData.coeff_type==CeedCoeff::Grid)
+   {
+      CeedGridCoeff* ceedCoeff = (CeedGridCoeff*)ceedData.coeff;
+      InitCeedTensorBasisAndRestriction(*ceedCoeff->coeff->FESpace(), ir, ceed, &ceedCoeff->basis,
+                   &ceedCoeff->restr);
+      CeedVectorCreate(ceed, ceedCoeff->coeff->FESpace()->GetNDofs(),
+                       &ceedCoeff->coeffVector);
+      CeedVectorSetArray(ceedCoeff->coeffVector, CEED_MEM_HOST, CEED_USE_POINTER,
+                         ceedCoeff->coeff->GetData());
+      CeedOperatorSetField(ceedData.build_oper, "coeff", ceedCoeff->restr,
+                           CEED_NOTRANSPOSE, ceedCoeff->basis, ceedCoeff->coeffVector);
+   }
+   CeedOperatorSetField(ceedData.build_oper, "dx", ceedData.mesh_restr, lmode,
+                        ceedData.mesh_basis, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(ceedData.build_oper, "weights", ceedData.mesh_restr_i,
+                        CEED_NOTRANSPOSE,
+                        ceedData.mesh_basis, CEED_VECTOR_NONE);
+   CeedOperatorSetField(ceedData.build_oper, "rho", ceedData.restr_i,
+                        CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
+
+   // Compute the quadrature data for the diff operator.
+   CeedOperatorApply(ceedData.build_oper, ceedData.node_coords, ceedData.rho,
+                     CEED_REQUEST_IMMEDIATE);
+
+   // Create the Q-function that defines the action of the diff operator.
+   diff_qf = diff_qf_file + ":f_apply_diff";
+   CeedQFunctionCreateInterior(ceed, 1, f_apply_diff,
+                               diff_qf.c_str(),
+                               &ceedData.apply_qfunc);
+   CeedQFunctionAddInput(ceedData.apply_qfunc, "u", dim, CEED_EVAL_GRAD);
+   CeedQFunctionAddInput(ceedData.apply_qfunc, "rho", dim * (dim + 1) / 2,
+                         CEED_EVAL_NONE);
+   CeedQFunctionAddOutput(ceedData.apply_qfunc, "v", dim, CEED_EVAL_GRAD);
+   CeedQFunctionSetContext(ceedData.apply_qfunc, &ceedData.build_ctx,
+                           sizeof(ceedData.build_ctx));
+
+   // Create the diff operator.
+   CeedOperatorCreate(ceed, ceedData.apply_qfunc, NULL, NULL, &ceedData.oper);
+   CeedOperatorSetField(ceedData.oper, "u", ceedData.restr, CEED_NOTRANSPOSE,
+                        ceedData.basis, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(ceedData.oper, "rho", ceedData.restr_i, CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, ceedData.rho);
+   CeedOperatorSetField(ceedData.oper, "v", ceedData.restr, CEED_NOTRANSPOSE,
+                        ceedData.basis, CEED_VECTOR_ACTIVE);
+
+   CeedVectorCreate(ceed, fes.GetNDofs(), &ceedData.u);
+   CeedVectorCreate(ceed, fes.GetNDofs(), &ceedData.v);
+}
+
+#endif
+
+}

fem/libceed/diffusion.h

@@ -0,0 +1,200 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+/// A structure used to pass additional data to f_build_diff and f_apply_diff
+struct BuildContext { CeedInt dim, space_dim; CeedScalar coeff; };
+
+/// libCEED Q-function for building quadrature data for a diffusion operator with a constant coefficient
+CEED_QFUNCTION(f_build_diff_const)(void *ctx, CeedInt Q,
+                                   const CeedScalar *const *in, CeedScalar *const *out)
+{
+   BuildContext *bc = (BuildContext*)ctx;
+   // in[0] is Jacobians with shape [dim, nc=dim, Q]
+   // in[1] is quadrature weights, size (Q)
+   //
+   // At every quadrature point, compute qw/det(J).adj(J).adj(J)^T and store
+   // the symmetric part of the result.
+   const CeedScalar coeff = bc->coeff;
+   const CeedScalar *J = in[0], *qw = in[1];
+   CeedScalar *qd = out[0];
+   switch (bc->dim + 10 * bc->space_dim)
+   {
+      case 11:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            qd[i] = coeff * qw[i] / J[i];
+         }
+         break;
+      case 22:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            // J: 0 2   qd: 0 1   adj(J):  J22 -J12
+            //    1 3       1 2           -J21  J11
+            const CeedScalar J11 = J[i + Q * 0];
+            const CeedScalar J21 = J[i + Q * 1];
+            const CeedScalar J12 = J[i + Q * 2];
+            const CeedScalar J22 = J[i + Q * 3];
+            const CeedScalar w = qw[i] / (J11 * J22 - J21 * J12);
+            qd[i + Q * 0] =   coeff * w * (J12 * J12 + J22 * J22);
+            qd[i + Q * 1] = - coeff * w * (J11 * J12 + J21 * J22);
+            qd[i + Q * 2] =   coeff * w * (J11 * J11 + J21 * J21);
+         }
+         break;
+      case 33:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            // J: 0 3 6   qd: 0 1 2
+            //    1 4 7       1 3 4
+            //    2 5 8       2 4 5
+            const CeedScalar J11 = J[i + Q * 0];
+            const CeedScalar J21 = J[i + Q * 1];
+            const CeedScalar J31 = J[i + Q * 2];
+            const CeedScalar J12 = J[i + Q * 3];
+            const CeedScalar J22 = J[i + Q * 4];
+            const CeedScalar J32 = J[i + Q * 5];
+            const CeedScalar J13 = J[i + Q * 6];
+            const CeedScalar J23 = J[i + Q * 7];
+            const CeedScalar J33 = J[i + Q * 8];
+            const CeedScalar A11 = J22 * J33 - J23 * J32;
+            const CeedScalar A12 = J13 * J32 - J12 * J33;
+            const CeedScalar A13 = J12 * J23 - J13 * J22;
+            const CeedScalar A21 = J23 * J31 - J21 * J33;
+            const CeedScalar A22 = J11 * J33 - J13 * J31;
+            const CeedScalar A23 = J13 * J21 - J11 * J23;
+            const CeedScalar A31 = J21 * J32 - J22 * J31;
+            const CeedScalar A32 = J12 * J31 - J11 * J32;
+            const CeedScalar A33 = J11 * J22 - J12 * J21;
+            const CeedScalar w = qw[i] / (J11 * A11 + J21 * A12 + J31 * A13);
+            qd[i + Q * 0] = coeff * w * (A11 * A11 + A12 * A12 + A13 * A13);
+            qd[i + Q * 1] = coeff * w * (A11 * A21 + A12 * A22 + A13 * A23);
+            qd[i + Q * 2] = coeff * w * (A11 * A31 + A12 * A32 + A13 * A33);
+            qd[i + Q * 3] = coeff * w * (A21 * A21 + A22 * A22 + A23 * A23);
+            qd[i + Q * 4] = coeff * w * (A21 * A31 + A22 * A32 + A23 * A33);
+            qd[i + Q * 5] = coeff * w * (A31 * A31 + A32 * A32 + A33 * A33);
+         }
+         break;
+   }
+   return 0;
+}
+
+/// libCEED Q-function for building quadrature data for a diffusion operator with a grid function coefficient
+CEED_QFUNCTION(f_build_diff_grid)(void *ctx, CeedInt Q,
+                                  const CeedScalar *const *in, CeedScalar *const *out)
+{
+   BuildContext *bc = (BuildContext *)ctx;
+   // in[1] is Jacobians with shape [dim, nc=dim, Q]
+   // in[2] is quadrature weights, size (Q)
+   //
+   // At every quadrature point, compute qw/det(J).adj(J).adj(J)^T and store
+   // the symmetric part of the result.
+   const CeedScalar *c = in[0], *J = in[1], *qw = in[2];
+   CeedScalar *qd = out[0];
+   switch (bc->dim + 10 * bc->space_dim)
+   {
+      case 11:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            qd[i] = c[i] * qw[i] / J[i];
+         }
+         break;
+      case 22:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            // J: 0 2   qd: 0 1   adj(J):  J22 -J12
+            //    1 3       1 2           -J21  J11
+            const CeedScalar coeff = c[i];
+            const CeedScalar J11 = J[i + Q * 0];
+            const CeedScalar J21 = J[i + Q * 1];
+            const CeedScalar J12 = J[i + Q * 2];
+            const CeedScalar J22 = J[i + Q * 3];
+            const CeedScalar w = qw[i] / (J11 * J22 - J21 * J12);
+            qd[i + Q * 0] =   coeff * w * (J12 * J12 + J22 * J22);
+            qd[i + Q * 1] = - coeff * w * (J11 * J12 + J21 * J22);
+            qd[i + Q * 2] =   coeff * w * (J11 * J11 + J21 * J21);
+         }
+         break;
+      case 33:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            // J: 0 3 6   qd: 0 1 2
+            //    1 4 7       1 3 4
+            //    2 5 8       2 4 5
+            const CeedScalar coeff = c[i];
+            const CeedScalar J11 = J[i + Q * 0];
+            const CeedScalar J21 = J[i + Q * 1];
+            const CeedScalar J31 = J[i + Q * 2];
+            const CeedScalar J12 = J[i + Q * 3];
+            const CeedScalar J22 = J[i + Q * 4];
+            const CeedScalar J32 = J[i + Q * 5];
+            const CeedScalar J13 = J[i + Q * 6];
+            const CeedScalar J23 = J[i + Q * 7];
+            const CeedScalar J33 = J[i + Q * 8];
+            const CeedScalar A11 = J22 * J33 - J23 * J32;
+            const CeedScalar A12 = J13 * J32 - J12 * J33;
+            const CeedScalar A13 = J12 * J23 - J13 * J22;
+            const CeedScalar A21 = J23 * J31 - J21 * J33;
+            const CeedScalar A22 = J11 * J33 - J13 * J31;
+            const CeedScalar A23 = J13 * J21 - J11 * J23;
+            const CeedScalar A31 = J21 * J32 - J22 * J31;
+            const CeedScalar A32 = J12 * J31 - J11 * J32;
+            const CeedScalar A33 = J11 * J22 - J12 * J21;
+            const CeedScalar w = qw[i] / (J11 * A11 + J21 * A12 + J31 * A13);
+            qd[i + Q * 0] = coeff * w * (A11 * A11 + A12 * A12 + A13 * A13);
+            qd[i + Q * 1] = coeff * w * (A11 * A21 + A12 * A22 + A13 * A23);
+            qd[i + Q * 2] = coeff * w * (A11 * A31 + A12 * A32 + A13 * A33);
+            qd[i + Q * 3] = coeff * w * (A21 * A21 + A22 * A22 + A23 * A23);
+            qd[i + Q * 4] = coeff * w * (A21 * A31 + A22 * A32 + A23 * A33);
+            qd[i + Q * 5] = coeff * w * (A31 * A31 + A32 * A32 + A33 * A33);
+         }
+         break;
+   }
+   return 0;
+}
+
+/// libCEED Q-function for applying a diff operator
+CEED_QFUNCTION(f_apply_diff)(void *ctx, CeedInt Q,
+                             const CeedScalar *const *in, CeedScalar *const *out)
+{
+   BuildContext *bc = (BuildContext *)ctx;
+   // in[0], out[0] have shape [dim, nc=1, Q]
+   const CeedScalar *ug = in[0], *qd = in[1];
+   CeedScalar *vg = out[0];
+   switch (bc->dim)
+   {
+      case 1:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            vg[i] = ug[i] * qd[i];
+         }
+         break;
+      case 2:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            const CeedScalar ug0 = ug[i + Q * 0];
+            const CeedScalar ug1 = ug[i + Q * 1];
+            vg[i + Q * 0] = qd[i + Q * 0] * ug0 + qd[i + Q * 1] * ug1;
+            vg[i + Q * 1] = qd[i + Q * 1] * ug0 + qd[i + Q * 2] * ug1;
+         }
+         break;
+      case 3:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            const CeedScalar ug0 = ug[i + Q * 0];
+            const CeedScalar ug1 = ug[i + Q * 1];
+            const CeedScalar ug2 = ug[i + Q * 2];
+            vg[i + Q * 0] = qd[i + Q * 0] * ug0 + qd[i + Q * 1] * ug1 + qd[i + Q * 2] * ug2;
+            vg[i + Q * 1] = qd[i + Q * 1] * ug0 + qd[i + Q * 3] * ug1 + qd[i + Q * 4] * ug2;
+            vg[i + Q * 2] = qd[i + Q * 2] * ug0 + qd[i + Q * 4] * ug1 + qd[i + Q * 5] * ug2;
+         }
+         break;
+   }
+   return 0;
+}

fem/libceed/diffusion.hpp

@@ -0,0 +1,31 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_LIBCEED_DIFF_HPP
+#define MFEM_LIBCEED_DIFF_HPP
+
+#include "ceed.hpp"
+#include "../fespace.hpp"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+/// Initialize a Diffusion Integrator using libCEED
+void CeedPADiffusionAssemble(const FiniteElementSpace &fes,
+                             const mfem::IntegrationRule &ir,  CeedData& ceedData);
+
+#endif
+
+}
+
+#endif // MFEM_LIBCEED_DIFF_HPP

fem/libceed/mass.cpp

@@ -0,0 +1,141 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "ceed.hpp"
+#include "../../general/device.hpp"
+
+#include "mass.h"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+void CeedPAMassAssemble(const FiniteElementSpace &fes,
+                        const mfem::IntegrationRule &irm, CeedData& ceedData)
+{
+   Ceed ceed(internal::ceed);
+   mfem::Mesh *mesh = fes.GetMesh();
+   const int ir_order = irm.GetOrder();
+   const mfem::IntegrationRule &ir =
+      mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order);
+   CeedInt nqpts, nelem = mesh->GetNE();
+
+   InitCeedTensorBasisAndRestriction(fes, ir, ceed, &ceedData.basis, &ceedData.restr);
+   mesh->EnsureNodes();
+   const mfem::FiniteElementSpace *mesh_fes = mesh->GetNodalFESpace();
+   MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
+   InitCeedTensorBasisAndRestriction(*mesh_fes, ir, ceed, &ceedData.mesh_basis, &ceedData.mesh_restr);
+   CeedBasisGetNumQuadraturePoints(ceedData.basis, &nqpts);
+
+   CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts,
+                                     nqpts*nelem, 1, &ceedData.restr_i);
+   CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts,
+                                     nqpts*nelem, 1, &ceedData.mesh_restr_i);
+
+   CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &ceedData.node_coords);
+   CeedVectorSetArray(ceedData.node_coords, CEED_MEM_HOST, CEED_USE_POINTER,
+                      mesh->GetNodes()->GetData());
+
+   CeedVectorCreate(ceed, nelem*nqpts, &ceedData.rho);
+
+   // Context data to be passed to the 'f_build_mass' Q-function.
+   ceedData.build_ctx.dim = mesh->Dimension();
+   ceedData.build_ctx.space_dim = mesh->SpaceDimension();
+
+   std::string mass_qf_file = GetCeedPath() + "/mass.h";
+   std::string mass_qf;
+
+   // Create the Q-function that builds the mass operator (i.e. computes its
+   // quadrature data) and set its context data.
+   switch (ceedData.coeff_type)
+   {
+      case CeedCoeff::Const:
+         mass_qf = mass_qf_file + ":f_build_mass_const";
+         CeedQFunctionCreateInterior(ceed, 1, f_build_mass_const,
+                                     mass_qf.c_str(),
+                                     &ceedData.build_qfunc);
+         ceedData.build_ctx.coeff = ((CeedConstCoeff*)ceedData.coeff)->val;
+         break;
+      case CeedCoeff::Grid:
+         mass_qf = mass_qf_file + ":f_build_mass_grid";
+         CeedQFunctionCreateInterior(ceed, 1, f_build_mass_grid,
+                                     mass_qf.c_str(),
+                                     &ceedData.build_qfunc);
+         CeedQFunctionAddInput(ceedData.build_qfunc, "coeff", 1, CEED_EVAL_INTERP);
+         break;
+      default:
+         MFEM_ABORT("This coeff_type is not handled");
+   }
+   CeedQFunctionAddInput(ceedData.build_qfunc, "dx", mesh->SpaceDimension()*mesh->SpaceDimension(),
+                         CEED_EVAL_GRAD);
+   CeedQFunctionAddInput(ceedData.build_qfunc, "weights", 1, CEED_EVAL_WEIGHT);
+   CeedQFunctionAddOutput(ceedData.build_qfunc, "rho", 1, CEED_EVAL_NONE);
+   CeedQFunctionSetContext(ceedData.build_qfunc, &ceedData.build_ctx,
+                           sizeof(ceedData.build_ctx));
+
+   // Create the operator that builds the quadrature data for the mass operator.
+   CeedOperatorCreate(ceed, ceedData.build_qfunc, NULL, NULL,
+                      &ceedData.build_oper);
+   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
+   if (mesh_fes->GetOrdering()==Ordering::byVDIM)
+   {
+      lmode = CEED_TRANSPOSE;
+   }
+   if (ceedData.coeff_type==CeedCoeff::Grid)
+   {
+      CeedGridCoeff* ceedCoeff = (CeedGridCoeff*)ceedData.coeff;
+      InitCeedTensorBasisAndRestriction(*ceedCoeff->coeff->FESpace(), ir, ceed, &ceedCoeff->basis,
+                   &ceedCoeff->restr);
+      CeedVectorCreate(ceed, ceedCoeff->coeff->FESpace()->GetNDofs(),
+                       &ceedCoeff->coeffVector);
+      CeedVectorSetArray(ceedCoeff->coeffVector, CEED_MEM_HOST, CEED_USE_POINTER,
+                         ceedCoeff->coeff->GetData());
+      CeedOperatorSetField(ceedData.build_oper, "coeff", ceedCoeff->restr, CEED_NOTRANSPOSE,
+                           ceedCoeff->basis, ceedCoeff->coeffVector);
+   }
+   CeedOperatorSetField(ceedData.build_oper, "dx", ceedData.mesh_restr, lmode,
+                        ceedData.mesh_basis, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(ceedData.build_oper, "weights", ceedData.mesh_restr_i,
+                        CEED_NOTRANSPOSE,
+                        ceedData.mesh_basis, CEED_VECTOR_NONE);
+   CeedOperatorSetField(ceedData.build_oper, "rho", ceedData.restr_i,
+                        CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
+
+   // Compute the quadrature data for the mass operator.
+   CeedOperatorApply(ceedData.build_oper, ceedData.node_coords, ceedData.rho,
+                     CEED_REQUEST_IMMEDIATE);
+
+   // Create the Q-function that defines the action of the mass operator.
+   mass_qf = mass_qf_file + ":f_apply_mass";
+   CeedQFunctionCreateInterior(ceed, 1, f_apply_mass,
+                               mass_qf.c_str(), &ceedData.apply_qfunc);
+   CeedQFunctionAddInput(ceedData.apply_qfunc, "u", 1, CEED_EVAL_INTERP);
+   CeedQFunctionAddInput(ceedData.apply_qfunc, "rho", 1, CEED_EVAL_NONE);
+   CeedQFunctionAddOutput(ceedData.apply_qfunc, "v", 1, CEED_EVAL_INTERP);
+
+   // Create the mass operator.
+   CeedOperatorCreate(ceed, ceedData.apply_qfunc, NULL, NULL, &ceedData.oper);
+   CeedOperatorSetField(ceedData.oper, "u", ceedData.restr, CEED_NOTRANSPOSE,
+                        ceedData.basis, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(ceedData.oper, "rho", ceedData.restr_i, CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, ceedData.rho);
+   CeedOperatorSetField(ceedData.oper, "v", ceedData.restr, CEED_NOTRANSPOSE,
+                        ceedData.basis, CEED_VECTOR_ACTIVE);
+
+   CeedVectorCreate(ceed, fes.GetNDofs(), &ceedData.u);
+   CeedVectorCreate(ceed, fes.GetNDofs(), &ceedData.v);
+}
+
+#endif
+
+}

fem/libceed/mass.h

@@ -0,0 +1,107 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+/// A structure used to pass additional data to f_build_diff and f_apply_diff
+struct BuildContext { CeedInt dim, space_dim; CeedScalar coeff; };
+
+/// libCEED Q-function for building quadrature data for a mass operator with a constant coefficient
+CEED_QFUNCTION(f_build_mass_const)(void *ctx, CeedInt Q,
+                                   const CeedScalar *const *in, CeedScalar *const *out)
+{
+   // in[0] is Jacobians with shape [dim, nc=dim, Q]
+   // in[1] is quadrature weights, size (Q)
+   BuildContext *bc = (BuildContext *)ctx;
+   const CeedScalar coeff = bc->coeff;
+   const CeedScalar *J = in[0], *qw = in[1];
+   CeedScalar *rho = out[0];
+   switch (bc->dim + 10*bc->space_dim)
+   {
+      case 11:
+         for (CeedInt i=0; i<Q; i++)
+         {
+            rho[i] = coeff * J[i] * qw[i];
+         }
+         break;
+      case 22:
+         for (CeedInt i=0; i<Q; i++)
+         {
+            // 0 2
+            // 1 3
+            rho[i] = coeff * (J[i+Q*0]*J[i+Q*3] - J[i+Q*1]*J[i+Q*2]) * qw[i];
+         }
+         break;
+      case 33:
+         for (CeedInt i=0; i<Q; i++)
+         {
+            // 0 3 6
+            // 1 4 7
+            // 2 5 8
+            rho[i] = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) -
+                      J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) +
+                      J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])) * coeff * qw[i];
+         }
+         break;
+   }
+   return 0;
+}
+
+/// libCEED Q-function for building quadrature data for a mass operator with a grid function coefficient
+CEED_QFUNCTION(f_build_mass_grid)(void *ctx, CeedInt Q,
+                                  const CeedScalar *const *in, CeedScalar *const *out)
+{
+   // in[0] is Jacobians with shape [dim, nc=dim, Q]
+   // in[1] is quadrature weights, size (Q)
+   BuildContext *bc = (BuildContext *)ctx;
+   const CeedScalar *c = in[0], *J = in[1], *qw = in[2];
+   CeedScalar *rho = out[0];
+   switch (bc->dim + 10*bc->space_dim)
+   {
+      case 11:
+         for (CeedInt i=0; i<Q; i++)
+         {
+            rho[i] = c[i] * J[i] * qw[i];
+         }
+         break;
+      case 22:
+         for (CeedInt i=0; i<Q; i++)
+         {
+            // 0 2
+            // 1 3
+            rho[i] = c[i] * (J[i+Q*0]*J[i+Q*3] - J[i+Q*1]*J[i+Q*2]) * qw[i];
+         }
+         break;
+      case 33:
+         for (CeedInt i=0; i<Q; i++)
+         {
+            // 0 3 6
+            // 1 4 7
+            // 2 5 8
+            rho[i] = (J[i+Q*0]*(J[i+Q*4]*J[i+Q*8] - J[i+Q*5]*J[i+Q*7]) -
+                      J[i+Q*1]*(J[i+Q*3]*J[i+Q*8] - J[i+Q*5]*J[i+Q*6]) +
+                      J[i+Q*2]*(J[i+Q*3]*J[i+Q*7] - J[i+Q*4]*J[i+Q*6])) * c[i] * qw[i];
+         }
+         break;
+   }
+   return 0;
+}
+
+/// libCEED Q-function for applying a mass operator
+CEED_QFUNCTION(f_apply_mass)(void *ctx, CeedInt Q,
+                             const CeedScalar *const *in, CeedScalar *const *out)
+{
+   const CeedScalar *u = in[0], *w = in[1];
+   CeedScalar *v = out[0];
+   for (CeedInt i=0; i<Q; i++)
+   {
+      v[i] = w[i] * u[i];
+   }
+   return 0;
+}

fem/libceed/mass.hpp

@@ -0,0 +1,31 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_LIBCEED_MASS_HPP
+#define MFEM_LIBCEED_MASS_HPP
+
+#include "ceed.hpp"
+#include "../fespace.hpp"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+/// Initialize a Mass Integrator using libCEED
+void CeedPAMassAssemble(const FiniteElementSpace &fes,
+                        const mfem::IntegrationRule &ir,  CeedData& ceedData);
+
+#endif
+
+}
+
+#endif // MFEM_LIBCEED_MASS_HPP

fem/libceed/mechanics.cpp

@@ -0,0 +1,128 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "ceed.hpp"
+#include "../../general/device.hpp"
+
+#include "mechanics.h"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+void CeedPAMechanicsAssemble(const FiniteElementSpace &fes,
+                             const mfem::IntegrationRule &irm,
+                             const double* ktan_ptr,
+                             CeedData& ceedData)
+{
+   Ceed ceed(internal::ceed);
+   mfem::Mesh *mesh = fes.GetMesh();
+   const int ir_order = irm.GetOrder();
+   const mfem::IntegrationRule &ir =
+      mfem::IntRules.Get(mfem::Geometry::SEGMENT, ir_order);
+   CeedInt nqpts, nelem = mesh->GetNE(), dim = mesh->SpaceDimension();
+   mesh->EnsureNodes();
+   InitCeedTensorBasisAndRestriction(fes, ir, ceed, &ceedData.basis, &ceedData.restr);
+
+   const mfem::FiniteElementSpace *mesh_fes = mesh->GetNodalFESpace();
+   MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
+   InitCeedTensorBasisAndRestriction(*mesh_fes, ir, ceed, &ceedData.mesh_basis, &ceedData.mesh_restr);
+   CeedBasisGetNumQuadraturePoints(ceedData.basis, &nqpts);
+
+   CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts,
+                                     nqpts * nelem, dim * dim * dim * dim, &ceedData.restr_i);
+   CeedElemRestrictionCreateIdentity(ceed, nelem, nqpts,
+                                     nqpts * nelem, 1, &ceedData.mesh_restr_i);
+
+   CeedVectorCreate(ceed, mesh->GetNodes()->Size(), &ceedData.node_coords);
+   CeedVectorSetArray(ceedData.node_coords, CEED_MEM_HOST, CEED_USE_POINTER,
+                      mesh->GetNodes()->GetData());
+
+   CeedVectorCreate(ceed, nelem * nqpts * dim * dim * dim * dim, &ceedData.rho);
+
+   // Context data to be passed to the 'f_build_diff' Q-function.
+   ceedData.build_ctx.dim = mesh->Dimension();
+   ceedData.build_ctx.space_dim = mesh->SpaceDimension();
+
+   std::string mech_qf_file = GetCeedPath() + "/mechanics.h";
+   std::string mech_qf;
+
+   // Create the Q-function that builds the diff operator (i.e. computes its
+   // quadrature data) and set its context data.
+   mech_qf = mech_qf_file + ":f_build_mech";
+   CeedQFunctionCreateInterior(ceed, 1, f_build_mech,
+                               mech_qf.c_str(),
+                               &ceedData.build_qfunc);
+   CeedQFunctionAddInput(ceedData.build_qfunc, "ktan", dim * dim * dim * dim,
+                          CEED_EVAL_NONE);
+   CeedQFunctionAddInput(ceedData.build_qfunc, "dx", dim * dim, CEED_EVAL_GRAD);
+   CeedQFunctionAddInput(ceedData.build_qfunc, "weights", 1, CEED_EVAL_WEIGHT);
+   CeedQFunctionAddOutput(ceedData.build_qfunc, "qd", dim * dim * dim * dim,
+                          CEED_EVAL_NONE);
+   CeedQFunctionSetContext(ceedData.build_qfunc, &ceedData.build_ctx,
+                           sizeof(ceedData.build_ctx));
+
+   // Create the operator that builds the quadrature data for the diff operator.
+   CeedOperatorCreate(ceed, ceedData.build_qfunc, NULL, NULL,
+                      &ceedData.build_oper);
+   CeedTransposeMode lmode = CEED_NOTRANSPOSE;
+   if (mesh_fes->GetOrdering()==Ordering::byVDIM)
+   {
+      lmode = CEED_TRANSPOSE;
+   }
+   CeedVector ktan;
+   CeedVectorCreate(ceed, nelem * nqpts * dim * dim * dim * dim, &ktan);
+   CeedVectorSetArray(ktan, CEED_MEM_DEVICE, CEED_USE_POINTER, const_cast<double*>(ktan_ptr));
+   CeedOperatorSetField(ceedData.build_oper, "ktan", ceedData.restr_i,
+                        CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, ktan);
+   CeedOperatorSetField(ceedData.build_oper, "dx", ceedData.mesh_restr, lmode,
+                        ceedData.mesh_basis, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(ceedData.build_oper, "weights", ceedData.mesh_restr_i,
+                        CEED_NOTRANSPOSE,
+                        ceedData.mesh_basis, CEED_VECTOR_NONE);
+   CeedOperatorSetField(ceedData.build_oper, "qd", ceedData.restr_i,
+                        CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
+
+   // Compute the quadrature data for the diff operator.
+   CeedOperatorApply(ceedData.build_oper, ceedData.node_coords, ceedData.rho,
+                     CEED_REQUEST_IMMEDIATE);
+
+   // Create the Q-function that defines the action of the diff operator.
+   mech_qf = mech_qf_file + ":f_apply_mech";
+   CeedQFunctionCreateInterior(ceed, 1, f_apply_mech,
+                               mech_qf.c_str(),
+                               &ceedData.apply_qfunc);
+   CeedQFunctionAddInput(ceedData.apply_qfunc, "u", dim * dim, CEED_EVAL_GRAD);
+   CeedQFunctionAddInput(ceedData.apply_qfunc, "qd", dim * dim * dim * dim,
+                         CEED_EVAL_NONE);
+   CeedQFunctionAddOutput(ceedData.apply_qfunc, "v", dim * dim, CEED_EVAL_GRAD);
+   CeedQFunctionSetContext(ceedData.apply_qfunc, &ceedData.build_ctx,
+                           sizeof(ceedData.build_ctx));
+
+   // Create the diff operator.
+   CeedOperatorCreate(ceed, ceedData.apply_qfunc, NULL, NULL, &ceedData.oper);
+   CeedOperatorSetField(ceedData.oper, "u", ceedData.restr, CEED_NOTRANSPOSE,
+                        ceedData.basis, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(ceedData.oper, "qd", ceedData.restr_i, CEED_NOTRANSPOSE,
+                        CEED_BASIS_COLLOCATED, ceedData.rho);
+   CeedOperatorSetField(ceedData.oper, "v", ceedData.restr, CEED_NOTRANSPOSE,
+                        ceedData.basis, CEED_VECTOR_ACTIVE);
+
+   CeedVectorCreate(ceed, fes.GetNDofs(), &ceedData.u);
+   CeedVectorCreate(ceed, fes.GetNDofs(), &ceedData.v);
+}
+
+#endif
+
+}

fem/libceed/mechanics.h

@@ -0,0 +1,181 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+/// A structure used to pass additional data to f_build_diff and f_apply_diff
+struct BuildContext { CeedInt dim, space_dim; CeedScalar coeff; };
+
+/// libCEED Q-function for building quadrature data for a diffusion operator with a constant coefficient
+CEED_QFUNCTION(f_build_mech)(void *ctx, CeedInt Q,
+                                   const CeedScalar *const *in, CeedScalar *const *out)
+{
+   BuildContext *bc = (BuildContext*)ctx;
+   // in[0] is ktan
+   // in[1] is Jacobians with shape [dim, nc=dim, Q]
+   // in[2] is quadrature weights, size (Q)
+   //
+   // At every quadrature point, compute qw/det(J).adj(J).adj(J)^T and store
+   // the symmetric part of the result.
+   const CeedScalar *ktan = in[0], *J = in[1], *qw = in[2];
+   CeedScalar *qd = out[0];
+   switch (bc->dim + 10 * bc->space_dim)
+   {
+      case 11:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            qd[i] = ktan[i] * qw[i] / J[i];
+         }
+         break;
+      case 22:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            // J: 0 2   qd: 0 1   adj(J):  J22 -J12
+            //    1 3       1 2           -J21  J11
+            // const CeedScalar J11 = J[i + Q * 0];
+            // const CeedScalar J21 = J[i + Q * 1];
+            // const CeedScalar J12 = J[i + Q * 2];
+            // const CeedScalar J22 = J[i + Q * 3];
+            // const CeedScalar w = qw[i] / (J11 * J22 - J21 * J12);
+            // qd[i + Q * 0] =   coeff * w * (J12 * J12 + J22 * J22);
+            // qd[i + Q * 1] = - coeff * w * (J11 * J12 + J21 * J22);
+            // qd[i + Q * 2] =   coeff * w * (J11 * J11 + J21 * J21);
+            //TODO
+         }
+         break;
+      case 33:
+         for (CeedInt q = 0; q < Q; q++)
+         {
+            // J: 0 3 6   qd: 0 1 2
+            //    1 4 7       1 3 4
+            //    2 5 8       2 4 5
+            const CeedScalar J11 = J[q + Q * 0];
+            const CeedScalar J21 = J[q + Q * 1];
+            const CeedScalar J31 = J[q + Q * 2];
+            const CeedScalar J12 = J[q + Q * 3];
+            const CeedScalar J22 = J[q + Q * 4];
+            const CeedScalar J32 = J[q + Q * 5];
+            const CeedScalar J13 = J[q + Q * 6];
+            const CeedScalar J23 = J[q + Q * 7];
+            const CeedScalar J33 = J[q + Q * 8];
+            const CeedScalar A11 = J22 * J33 - J23 * J32;
+            const CeedScalar A12 = J13 * J32 - J12 * J33;
+            const CeedScalar A13 = J12 * J23 - J13 * J22;
+            const CeedScalar A21 = J23 * J31 - J21 * J33;
+            const CeedScalar A22 = J11 * J33 - J13 * J31;
+            const CeedScalar A23 = J13 * J21 - J11 * J23;
+            const CeedScalar A31 = J21 * J32 - J22 * J31;
+            const CeedScalar A32 = J12 * J31 - J11 * J32;
+            const CeedScalar A33 = J11 * J22 - J12 * J21;
+            const CeedScalar w = qw[q] / (J11 * A11 + J21 * A12 + J31 * A13);
+            // Load ktan
+            CeedScalar K[3][3][3][3];
+            for (int j = 0; j < 3; ++j) {
+               for (int k = 0; k < 3; ++k) {
+                  for (int l = 0; l < 3; ++l) {
+                     for (int m = 0; m < 3; ++m) {
+                        K[j][k][l][m] = ktan [q + (j + k*3 + l*3*3 + m*3*3*3) * Q];
+                     }
+                  }
+               }
+            }
+            // ktan*J^-1
+            CeedScalar tmp[3][3][3][3];
+            for (int j = 0; j < 3; ++j) {
+               for (int k = 0; k < 3; ++k) {
+                  for (int l = 0; l < 3; ++l) {
+                     tmp[j][k][l][0] = K[j][k][l][0] * A11 + K[j][k][l][1] * A21 + K[j][k][l][2] * A31;
+                     tmp[j][k][l][1] = K[j][k][l][0] * A12 + K[j][k][l][1] * A22 + K[j][k][l][2] * A32;
+                     tmp[j][k][l][2] = K[j][k][l][0] * A13 + K[j][k][l][1] * A23 + K[j][k][l][2] * A33;
+                  }
+               }
+            }
+            // J^-T*ktan*J^-1
+            for (int k = 0; k < 3; ++k) {
+               for (int l = 0; l < 3; ++l) {
+                  for (int n = 0; n < 3; ++n) {
+                     qd[q + (0 + k*3 + l*3*3 + n*3*3*3)*Q] = w * (A11 * tmp[0][k][l][n] + A21 * tmp[1][k][l][n] + A31 * tmp[2][k][l][n]);
+                     qd[q + (1 + k*3 + l*3*3 + n*3*3*3)*Q] = w * (A12 * tmp[0][k][l][n] + A22 * tmp[1][k][l][n] + A32 * tmp[2][k][l][n]);
+                     qd[q + (2 + k*3 + l*3*3 + n*3*3*3)*Q] = w * (A13 * tmp[0][k][l][n] + A23 * tmp[1][k][l][n] + A33 * tmp[2][k][l][n]);
+                  }
+               }
+            }
+         }
+         break;
+   }
+   return 0;
+}
+
+/// libCEED Q-function for applying a diff operator
+CEED_QFUNCTION(f_apply_mech)(void *ctx, CeedInt Q,
+                             const CeedScalar *const *in, CeedScalar *const *out)
+{
+   BuildContext *bc = (BuildContext *)ctx;
+   // in[0], out[0] have shape [dim, nc=1, Q]
+   const CeedScalar *ug = in[0], *qd = in[1];
+   CeedScalar *vg = out[0];
+   switch (bc->dim)
+   {
+      case 1:
+         for (CeedInt i = 0; i < Q; i++)
+         {
+            vg[i] = ug[i] * qd[i];
+         }
+         break;
+      case 2:
+         // for (CeedInt i = 0; i < Q; i++)
+         // {
+         //    const CeedScalar ug0 = ug[i + Q * 0];
+         //    const CeedScalar ug1 = ug[i + Q * 1];
+         //    vg[i + Q * 0] = qd[i + Q * 0] * ug0 + qd[i + Q * 1] * ug1;
+         //    vg[i + Q * 1] = qd[i + Q * 1] * ug0 + qd[i + Q * 2] * ug1;
+         // }
+         break;
+      case 3:
+         for (CeedInt q = 0; q < Q; q++)
+         {
+            // Read spatial derivatives of u components
+            const CeedScalar uJ[3][3]        = {{ug[q+(0+0*3)*Q],
+                                                 ug[q+(0+1*3)*Q],
+                                                 ug[q+(0+2*3)*Q]},
+                                                {ug[q+(1+0*3)*Q],
+                                                 ug[q+(1+1*3)*Q],
+                                                 ug[q+(1+2*3)*Q]},
+                                                {ug[q+(2+0*3)*Q],
+                                                 ug[q+(2+1*3)*Q],
+                                                 ug[q+(2+2*3)*Q]}
+                                                };
+            // Load quadrature data
+            CeedScalar K[3][3][3][3];
+            for (int j = 0; j < 3; ++j) {
+               for (int k = 0; k < 3; ++k) {
+                  for (int l = 0; l < 3; ++l) {
+                     for (int m = 0; m < 3; ++m) {
+                        K[j][k][l][m] = qd [q + (j + k*3 + l*3*3 + m*3*3*3) * Q];
+                     }
+                  }
+               }
+            }
+            // double contraction
+            for (int j = 0; j < 3; ++j) {
+               for (int k = 0; k < 3; ++k) {
+                  vg[q + (j+k*3)*Q] = 0.0;
+                  for (int l = 0; l < 3; ++l) {
+                     for (int m = 0; m < 3; ++m) {
+                        vg[q + (j+k*3)*Q] += K[j][k][l][m] * uJ[m][l];
+                     }
+                  }
+               }
+            }
+         }
+         break;
+   }
+   return 0;
+}
+

fem/libceed/mechanics.hpp

@@ -0,0 +1,33 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_LIBCEED_MECH_HPP
+#define MFEM_LIBCEED_MECH_HPP
+
+#include "ceed.hpp"
+#include "../fespace.hpp"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_CEED
+
+/// Initialize a Mechanics Integrator using libCEED
+void CeedPAMechanicsAssemble(const FiniteElementSpace &fes,
+                             const mfem::IntegrationRule &irm,
+                             const double* ktan_ptr,
+                             CeedData& ceedData);
+
+#endif
+
+}
+
+#endif // MFEM_LIBCEED_NECH_HPP

general/device.cpp

@@ -12,6 +12,9 @@
 #include "forall.hpp"
 #include "cuda.hpp"
 #include "occa.hpp"
+#ifdef MFEM_USE_CEED
+#include <ceed.h>
+#endif
 
 #include <string>
 #include <map>
@@ -29,20 +32,24 @@ namespace internal
 occa::device occaDevice;
 #endif
 
+#ifdef MFEM_USE_CEED
+Ceed ceed;
+#endif
+
 // Backends listed by priority, high to low:
 static const Backend::Id backend_list[Backend::NUM_BACKENDS] =
 {
-   Backend::OCCA_CUDA, Backend::RAJA_CUDA, Backend::CUDA,
+   Backend::OCCA_CUDA, Backend::RAJA_CUDA, Backend::CEED_CUDA, Backend::CUDA,
    Backend::HIP,
    Backend::OCCA_OMP, Backend::RAJA_OMP, Backend::OMP,
-   Backend::OCCA_CPU, Backend::RAJA_CPU, Backend::CPU
+   Backend::OCCA_CPU, Backend::RAJA_CPU, Backend::CEED_CPU, Backend::CPU
 };
 
 // Backend names listed by priority, high to low:
 static const char *backend_name[Backend::NUM_BACKENDS] =
 {
-   "occa-cuda", "raja-cuda", "cuda", "hip", "occa-omp", "raja-omp", "omp",
-   "occa-cpu", "raja-cpu", "cpu"
+   "occa-cuda", "raja-cuda", "ceed-cuda", "cuda", "hip", "occa-omp", "raja-omp", "omp",
+   "occa-cpu", "raja-cpu", "ceed-cpu", "cpu"
 };
 
 } // namespace mfem::internal
@@ -54,7 +61,11 @@ Device Device::device_singleton;
 
 Device::~Device()
 {
-   if (destroy_mm) { mm.Destroy(); }
+   if (destroy_mm)
+   {
+      free(ceed_option);
+      mm.Destroy();
+   }
 }
 
 void Device::Configure(const std::string &device, const int dev)
@@ -64,15 +75,29 @@ void Device::Configure(const std::string &device, const int dev)
    {
       bmap[internal::backend_name[i]] = internal::backend_list[i];
    }
-   std::string::size_type beg = 0, end;
+   std::string::size_type beg = 0, end, option;
    while (1)
    {
       end = device.find(',', beg);
       end = (end != std::string::npos) ? end : device.size();
       const std::string bname = device.substr(beg, end - beg);
-      std::map<std::string, Backend::Id>::iterator it = bmap.find(bname);
-      MFEM_VERIFY(it != bmap.end(), "invalid backend name: '" << bname << '\'');
-      Get().MarkBackend(it->second);
+      option = bname.find(':');
+      if (option==std::string::npos)//No option
+      {
+         const std::string backend = bname;
+         std::map<std::string, Backend::Id>::iterator it = bmap.find(backend);
+         MFEM_VERIFY(it != bmap.end(), "invalid backend name: '" << backend << '\'');
+         Get().MarkBackend(it->second);
+      }
+      else
+      {
+         const std::string backend = bname.substr(0, option);
+         const std::string boption = bname.substr(option+1);
+         Get().ceed_option = strdup(boption.c_str());
+         std::map<std::string, Backend::Id>::iterator it = bmap.find(backend);
+         MFEM_VERIFY(it != bmap.end(), "invalid backend name: '" << backend << '\'');
+         Get().MarkBackend(it->second);
+      }
       if (end == device.size()) { break; }
       beg = end + 1;
    }
@@ -82,6 +107,10 @@ void Device::Configure(const std::string &device, const int dev)
    {
       Get().MarkBackend(Backend::CUDA);
    }
+   if (Allows(Backend::CEED_CUDA))
+   {
+      Get().MarkBackend(Backend::CUDA);
+   }
 
    // Perform setup.
    Get().Setup(dev);
@@ -223,6 +252,20 @@ static void OccaDeviceSetup(const int dev)
 #endif
 }
 
+static void CeedDeviceSetup(const char* ceed_spec)
+{
+#ifdef MFEM_USE_CEED
+   CeedInit(ceed_spec, &internal::ceed);
+   const char *ceed_backend;
+   CeedGetResource(internal::ceed, &ceed_backend);
+   mfem::out << "libCEED backend: " << ceed_backend << std::endl;
+   if (strcmp(ceed_spec, ceed_backend) && strcmp(ceed_spec, "/cpu/self"))
+   {
+      std::cout << std::endl << "WARNING!!!\nlibCEED is not using the requested backend!!!\nWARNING!!!\n" << std::endl;
+   }
+#endif
+}
+
 void Device::Setup(const int device)
 {
    MFEM_VERIFY(ngpu == -1, "the mfem::Device is already configured!");
@@ -251,6 +294,28 @@ void Device::Setup(const int device)
    if (Allows(Backend::RAJA_CUDA)) { RajaDeviceSetup(dev, ngpu); }
    // The check for MFEM_USE_OCCA is in the function OccaDeviceSetup().
    if (Allows(Backend::OCCA_MASK)) { OccaDeviceSetup(dev); }
+   if (Allows(Backend::CEED_CPU))
+   {
+      if (!ceed_option)
+      {
+         CeedDeviceSetup("/cpu/self");
+      }
+      else
+      {
+         CeedDeviceSetup(ceed_option);
+      }
+   }
+   if (Allows(Backend::CEED_CUDA))
+   {
+      if (!ceed_option)
+      {
+         CeedDeviceSetup("/gpu/cuda/gen");
+      }
+      else
+      {
+         CeedDeviceSetup(ceed_option);
+      }
+   }
 }
 
 } // mfem

general/device.hpp

@@ -54,7 +54,14 @@ struct Backend
       OCCA_OMP = 1 << 8,
       /** @brief [device] OCCA CUDA backend. Enabled when MFEM_USE_OCCA = YES
           and MFEM_USE_CUDA = YES. */
-      OCCA_CUDA = 1 << 9
+      OCCA_CUDA = 1 << 9,
+      /** @brief [host] CEED backend: GPU backends can still be used, but
+          with expensive memory transfers.
+          Enabled when MFEM_USE_CEED = YES. */
+      CEED_CPU  = 1 << 10,
+      /** @brief [device] Ceed backends working in colaboration with the Cuda backend.
+          Enabled when MFEM_USE_CEED = YES and MFEM_USE_CUDA = YES. */
+      CEED_CUDA = 1 << 11
    };
 
    /** @brief Additional useful constants. For example, the *_MASK constants can
@@ -62,7 +69,7 @@ struct Backend
    enum
    {
       /// Number of backends: from (1 << 0) to (1 << (NUM_BACKENDS-1)).
-      NUM_BACKENDS = 10,
+      NUM_BACKENDS = 12,
 
       /// Biwise-OR of all CPU backends
       CPU_MASK = CPU | RAJA_CPU | OCCA_CPU,
@@ -72,6 +79,8 @@ struct Backend
       HIP_MASK = HIP,
       /// Biwise-OR of all OpenMP backends
       OMP_MASK = OMP | RAJA_OMP | OCCA_OMP,
+      /// Bitwise-OR of all CEED backends
+      CEED_MASK = CEED_CPU | CEED_CUDA,
       /// Biwise-OR of all device backends
       DEVICE_MASK = CUDA_MASK | HIP_MASK,
 
@@ -116,6 +125,7 @@ private:
    MemoryType mem_type;    ///< Current Device MemoryType
    MemoryClass mem_class;  ///< Current Device MemoryClass
 
+   char *ceed_option = NULL;
    Device(Device const&);
    void operator=(Device const&);
    static Device& Get() { return device_singleton; }
@@ -187,7 +197,13 @@ public:
        * Multiple backends can be configured at the same time.
        * Only one 'occa-*' backend can be configured at a time.
        * The backend 'occa-cuda' enables the 'cuda' backend unless 'raja-cuda'
-         is already enabled. */
+         is already enabled.
+       * The backend 'ceed-cpu' delegates to a libCEED CPU backend the setup and
+         evaluation of the operator.
+       * The backend 'ceed-cuda' delegates to a libCEED CUDA backend the setup
+         and evaluation of the operator and enables the 'cuda' backend to avoid
+         transfer between host and device.
+   */
    void Configure(const std::string &device, const int dev = 0);
 
    /// Print the configuration of the MFEM virtual device object.
@@ -220,7 +236,7 @@ public:
    /** @brief Get the current Device MemoryClass. This is the MemoryClass used
        by most MFEM device kernels to access Memory objects. */
    static inline MemoryClass GetMemoryClass() { return Get().mem_class; }
-
+   
    static void SetGPUAwareMPI(const bool force = true)
    { Get().mpi_gpu_aware = force; }
 

makefile

@@ -258,7 +258,7 @@ endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
- SUITESPARSE STRUMPACK GECKO GNUTLS NETCDF PETSC MPFR PUMI OCCA RAJA
+ SUITESPARSE STRUMPACK GECKO GNUTLS NETCDF PETSC MPFR PUMI OCCA CEED RAJA
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 
 define mfem_check_dependency
@@ -319,7 +319,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE MFEM_USE_GECKO\
  MFEM_USE_SUPERLU MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_NETCDF\
  MFEM_USE_PETSC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT MFEM_USE_PUMI\
- MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_OCCA MFEM_USE_RAJA\
+ MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_RAJA\
  MFEM_SOURCE_DIR MFEM_INSTALL_DIR
 
 # List of makefile variables that will be written to config.mk:
@@ -387,7 +387,7 @@ ifneq (,$(filter install,$(MAKECMDGOALS)))
 endif
 
 # Source dirs in logical order
-DIRS = general linalg mesh fem
+DIRS = general linalg mesh fem fem/libceed
 SOURCE_FILES = $(foreach dir,$(DIRS),$(wildcard $(SRC)$(dir)/*.cpp))
 RELSRC_FILES = $(patsubst $(SRC)%,%,$(SOURCE_FILES))
 OBJECT_FILES = $(patsubst $(SRC)%,$(BLD)%,$(SOURCE_FILES:.cpp=.o))
@@ -550,6 +550,9 @@ install: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT))
 	   mkdir -p $(PREFIX_INC)/mfem/$$dir && \
 	   $(INSTALL) -m 640 $(SRC)$$dir/*.okl $(PREFIX_INC)/mfem/$$dir; \
 	done
+# install libCEED q-function headers
+	mkdir -p $(PREFIX_INC)/mfem/fem/libceed
+	$(INSTALL) -m 640 $(SRC)fem/libceed/*.h $(PREFIX_INC)/mfem/fem/libceed
 # install config.mk in $(PREFIX_SHARE)
 	mkdir -p $(PREFIX_SHARE)
 	$(MAKE) -C $(BLD)config config-mk CONFIG_MK=config-install.mk
@@ -629,6 +632,7 @@ status info:
 	$(info MFEM_USE_HIP           = $(MFEM_USE_HIP))
 	$(info MFEM_USE_RAJA          = $(MFEM_USE_RAJA))
 	$(info MFEM_USE_OCCA          = $(MFEM_USE_OCCA))
+	$(info MFEM_USE_CEED          = $(MFEM_USE_CEED))
 	$(info MFEM_CXX               = $(value MFEM_CXX))
 	$(info MFEM_CPPFLAGS          = $(value MFEM_CPPFLAGS))
 	$(info MFEM_CXXFLAGS          = $(value MFEM_CXXFLAGS))