Solve 2d/3d Darcy and Richards problems using H(div) FE space

backends/ref/ceed-ref-basis.c

@@ -14,7 +14,7 @@
 #include "ceed-ref.h"
 
 //------------------------------------------------------------------------------
-// Basis Apply
+// Basis Apply H1
 //------------------------------------------------------------------------------
 static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector U, CeedVector V) {
   Ceed ceed;

examples/Hdiv-mass/Makefile

@@ -0,0 +1,82 @@
+# Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
+# Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
+# All Rights reserved. See files LICENSE and NOTICE for details.
+#
+# This file is part of CEED, a collection of benchmarks, miniapps, software
+# libraries and APIs for efficient high-order finite element and spectral
+# element discretizations for exascale applications. For more information and
+# source code availability see http://github.com/ceed.
+#
+# The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+# a collaborative effort of two U.S. Department of Energy organizations (Office
+# of Science and the National Nuclear Security Administration) responsible for
+# the planning and preparation of a capable exascale ecosystem, including
+# software, applications, hardware, advanced system engineering and early
+# testbed platforms, in support of the nation's exascale computing imperative.
+
+COMMON ?= ../../common.mk
+-include $(COMMON)
+
+# Note: PETSC_ARCH can be undefined or empty for installations which do not use
+#       PETSC_ARCH - for example when using PETSc installed through Spack.
+PETSc.pc := $(PETSC_DIR)/$(PETSC_ARCH)/lib/pkgconfig/PETSc.pc
+CEED_DIR ?= ../..
+ceed.pc := $(CEED_DIR)/lib/pkgconfig/ceed.pc
+
+CC = $(call pkgconf, --variable=ccompiler $(PETSc.pc) $(ceed.pc))
+CFLAGS = -std=c99 \
+  $(call pkgconf, --variable=cflags_extra $(PETSc.pc)) \
+  $(call pkgconf, --cflags-only-other $(PETSc.pc)) \
+  $(OPT)
+CPPFLAGS = $(call pkgconf, --cflags-only-I $(PETSc.pc) $(ceed.pc)) \
+  $(call pkgconf, --variable=cflags_dep $(PETSc.pc))
+LDFLAGS = $(call pkgconf, --libs-only-L --libs-only-other $(PETSc.pc) $(ceed.pc))
+LDFLAGS += $(patsubst -L%, $(call pkgconf, --variable=ldflag_rpath $(PETSc.pc))%, $(call pkgconf, --libs-only-L $(PETSc.pc) $(ceed.pc)))
+LDLIBS = $(call pkgconf, --libs-only-l $(PETSc.pc) $(ceed.pc)) -lm
+
+OBJDIR := build
+SRCDIR := src
+PROBLEMDIR := problems
+
+src.c := main.c $(sort $(wildcard $(PROBLEMDIR)/*.c)) $(sort $(wildcard $(SRCDIR)/*.c))
+src.o = $(src.c:%.c=$(OBJDIR)/%.o)
+
+all: main
+
+main: $(src.o) | $(PETSc.pc) $(ceed.pc)
+	$(call quiet,LINK.o) $(CEED_LDFLAGS) $^ $(LOADLIBES) $(LDLIBS) -o $@
+
+.SECONDEXPANSION: # to expand $$(@D)/.DIR
+%/.DIR :
+	@mkdir -p $(@D)
+	@touch $@
+
+# Quiet, color output
+quiet ?= $($(1))
+
+$(OBJDIR)/%.o : %.c | $$(@D)/.DIR
+	$(call quiet,CC) $(CPPFLAGS) $(CFLAGS) -c -o $@ $(abspath $<)
+
+# Rules for building the examples
+#%: %.c
+
+print: $(PETSc.pc) $(ceed.pc)
+	$(info CC      : $(CC))
+	$(info CFLAGS  : $(CFLAGS))
+	$(info CPPFLAGS: $(CPPFLAGS))
+	$(info LDFLAGS : $(LDFLAGS))
+	$(info LDLIBS  : $(LDLIBS))
+	@true
+
+clean:
+	$(RM) -r $(OBJDIR) main *.vtu *.csv
+
+$(PETSc.pc):
+	$(if $(wildcard $@),,$(error \
+	  PETSc config not found at $@. Please set PETSC_DIR and PETSC_ARCH))
+
+.PHONY: all print clean
+
+pkgconf = $(shell pkg-config $1 | sed -e 's/^"//g' -e 's/"$$//g')
+
+-include $(src.o:%.o=%.d)

examples/Hdiv-mass/basis/Hdiv-hex.h

@@ -0,0 +1,156 @@
+// Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
+// Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
+// All Rights reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+// To see how the nodal basis is constructed visit:
+// https://github.com/rezgarshakeri/H-div-Tests
+int NodalHdivBasisHex(CeedScalar *x, CeedScalar *Bx, CeedScalar *By, CeedScalar *Bz) {
+  Bx[0] = 0.0625 * x[0] * x[0] - 0.0625;
+  By[0] = -0.0625 * x[0] * x[1] * x[1] + 0.0625 * x[0] + 0.0625 * x[1] * x[1] - 0.0625;
+  Bz[0] = 0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] + 0.125 * x[0] - 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] -
+          0.125;
+  Bx[1] = 0.0625 - 0.0625 * x[0] * x[0];
+  By[1] = 0.0625 * x[0] * x[1] * x[1] - 0.0625 * x[0] + 0.0625 * x[1] * x[1] - 0.0625;
+  Bz[1] = -0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] - 0.125 * x[0] - 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] -
+          0.125;
+  Bx[2] = 0.0625 * x[0] * x[0] - 0.0625;
+  By[2] = 0.0625 * x[0] * x[1] * x[1] - 0.0625 * x[0] - 0.0625 * x[1] * x[1] + 0.0625;
+  Bz[2] = -0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] + 0.125 * x[0] + 0.125 * x[1] * x[2] - 0.125 * x[1] + 0.125 * x[2] -
+          0.125;
+  Bx[3] = 0.0625 - 0.0625 * x[0] * x[0];
+  By[3] = -0.0625 * x[0] * x[1] * x[1] + 0.0625 * x[0] - 0.0625 * x[1] * x[1] + 0.0625;
+  Bz[3] = 0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] - 0.125 * x[0] + 0.125 * x[1] * x[2] - 0.125 * x[1] + 0.125 * x[2] -
+          0.125;
+  Bx[4] = 0.0625 * x[0] * x[0] - 0.0625;
+  By[4] = -0.0625 * x[0] * x[1] * x[1] + 0.0625 * x[0] + 0.0625 * x[1] * x[1] - 0.0625;
+  Bz[4] = 0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] - 0.125 * x[0] - 0.125 * x[1] * x[2] - 0.125 * x[1] + 0.125 * x[2] +
+          0.125;
+  Bx[5] = 0.0625 - 0.0625 * x[0] * x[0];
+  By[5] = 0.0625 * x[0] * x[1] * x[1] - 0.0625 * x[0] + 0.0625 * x[1] * x[1] - 0.0625;
+  Bz[5] = -0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] + 0.125 * x[0] - 0.125 * x[1] * x[2] - 0.125 * x[1] + 0.125 * x[2] +
+          0.125;
+  Bx[6] = 0.0625 * x[0] * x[0] - 0.0625;
+  By[6] = 0.0625 * x[0] * x[1] * x[1] - 0.0625 * x[0] - 0.0625 * x[1] * x[1] + 0.0625;
+  Bz[6] = -0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] - 0.125 * x[0] + 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] +
+          0.125;
+  Bx[7] = 0.0625 - 0.0625 * x[0] * x[0];
+  By[7] = -0.0625 * x[0] * x[1] * x[1] + 0.0625 * x[0] - 0.0625 * x[1] * x[1] + 0.0625;
+  Bz[7] = 0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] + 0.125 * x[0] + 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] +
+          0.125;
+  Bx[8] = 0.0625 * x[0] * x[0] * x[2] - 0.0625 * x[0] * x[0] - 0.0625 * x[2] + 0.0625;
+  By[8] = -0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] - 0.125 * x[0] - 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] -
+          0.125;
+  Bz[8] = 0.0625 * x[2] * x[2] - 0.0625;
+  Bx[9] = -0.0625 * x[0] * x[0] * x[2] + 0.0625 * x[0] * x[0] + 0.0625 * x[2] - 0.0625;
+  By[9] = 0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] + 0.125 * x[0] - 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] -
+          0.125;
+  Bz[9]  = 0.0625 * x[2] * x[2] - 0.0625;
+  Bx[10] = -0.0625 * x[0] * x[0] * x[2] - 0.0625 * x[0] * x[0] + 0.0625 * x[2] + 0.0625;
+  By[10] = 0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] - 0.125 * x[0] + 0.125 * x[1] * x[2] + 0.125 * x[1] - 0.125 * x[2] -
+           0.125;
+  Bz[10] = 0.0625 - 0.0625 * x[2] * x[2];
+  Bx[11] = 0.0625 * x[0] * x[0] * x[2] + 0.0625 * x[0] * x[0] - 0.0625 * x[2] - 0.0625;
+  By[11] = -0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] + 0.125 * x[0] + 0.125 * x[1] * x[2] + 0.125 * x[1] -
+           0.125 * x[2] - 0.125;
+  Bz[11] = 0.0625 - 0.0625 * x[2] * x[2];
+  Bx[12] = 0.0625 * x[0] * x[0] * x[2] - 0.0625 * x[0] * x[0] - 0.0625 * x[2] + 0.0625;
+  By[12] = -0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] + 0.125 * x[0] - 0.125 * x[1] * x[2] + 0.125 * x[1] -
+           0.125 * x[2] + 0.125;
+  Bz[12] = 0.0625 * x[2] * x[2] - 0.0625;
+  Bx[13] = -0.0625 * x[0] * x[0] * x[2] + 0.0625 * x[0] * x[0] + 0.0625 * x[2] - 0.0625;
+  By[13] = 0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] - 0.125 * x[0] - 0.125 * x[1] * x[2] + 0.125 * x[1] - 0.125 * x[2] +
+           0.125;
+  Bz[13] = 0.0625 * x[2] * x[2] - 0.0625;
+  Bx[14] = -0.0625 * x[0] * x[0] * x[2] - 0.0625 * x[0] * x[0] + 0.0625 * x[2] + 0.0625;
+  By[14] = 0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] + 0.125 * x[0] + 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] +
+           0.125;
+  Bz[14] = 0.0625 - 0.0625 * x[2] * x[2];
+  Bx[15] = 0.0625 * x[0] * x[0] * x[2] + 0.0625 * x[0] * x[0] - 0.0625 * x[2] - 0.0625;
+  By[15] = -0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] - 0.125 * x[0] + 0.125 * x[1] * x[2] + 0.125 * x[1] +
+           0.125 * x[2] + 0.125;
+  Bz[15] = 0.0625 - 0.0625 * x[2] * x[2];
+  Bx[16] = 0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] + 0.125 * x[0] + 0.125 * x[1] * x[2] - 0.125 * x[1] - 0.125 * x[2] +
+           0.125;
+  By[16] = 0.0625 * x[1] * x[1] - 0.0625;
+  Bz[16] = -0.0625 * x[1] * x[2] * x[2] + 0.0625 * x[1] + 0.0625 * x[2] * x[2] - 0.0625;
+  Bx[17] = -0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] + 0.125 * x[0] - 0.125 * x[1] * x[2] + 0.125 * x[1] -
+           0.125 * x[2] + 0.125;
+  By[17] = 0.0625 - 0.0625 * x[1] * x[1];
+  Bz[17] = 0.0625 * x[1] * x[2] * x[2] - 0.0625 * x[1] + 0.0625 * x[2] * x[2] - 0.0625;
+  Bx[18] = -0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] + 0.125 * x[0] - 0.125 * x[1] * x[2] - 0.125 * x[1] +
+           0.125 * x[2] + 0.125;
+  By[18] = 0.0625 * x[1] * x[1] - 0.0625;
+  Bz[18] = 0.0625 * x[1] * x[2] * x[2] - 0.0625 * x[1] - 0.0625 * x[2] * x[2] + 0.0625;
+  Bx[19] = 0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] + 0.125 * x[0] + 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] +
+           0.125;
+  By[19] = 0.0625 - 0.0625 * x[1] * x[1];
+  Bz[19] = -0.0625 * x[1] * x[2] * x[2] + 0.0625 * x[1] - 0.0625 * x[2] * x[2] + 0.0625;
+  Bx[20] = 0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] + 0.125 * x[0] - 0.125 * x[1] * x[2] + 0.125 * x[1] + 0.125 * x[2] -
+           0.125;
+  By[20] = 0.0625 * x[1] * x[1] - 0.0625;
+  Bz[20] = -0.0625 * x[1] * x[2] * x[2] + 0.0625 * x[1] + 0.0625 * x[2] * x[2] - 0.0625;
+  Bx[21] = -0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] - 0.125 * x[0] * x[2] + 0.125 * x[0] + 0.125 * x[1] * x[2] - 0.125 * x[1] +
+           0.125 * x[2] - 0.125;
+  By[21] = 0.0625 - 0.0625 * x[1] * x[1];
+  Bz[21] = 0.0625 * x[1] * x[2] * x[2] - 0.0625 * x[1] + 0.0625 * x[2] * x[2] - 0.0625;
+  Bx[22] = -0.125 * x[0] * x[1] * x[2] - 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] + 0.125 * x[0] + 0.125 * x[1] * x[2] + 0.125 * x[1] -
+           0.125 * x[2] - 0.125;
+  By[22] = 0.0625 * x[1] * x[1] - 0.0625;
+  Bz[22] = 0.0625 * x[1] * x[2] * x[2] - 0.0625 * x[1] - 0.0625 * x[2] * x[2] + 0.0625;
+  Bx[23] = 0.125 * x[0] * x[1] * x[2] + 0.125 * x[0] * x[1] + 0.125 * x[0] * x[2] + 0.125 * x[0] - 0.125 * x[1] * x[2] - 0.125 * x[1] - 0.125 * x[2] -
+           0.125;
+  By[23] = 0.0625 - 0.0625 * x[1] * x[1];
+  Bz[23] = -0.0625 * x[1] * x[2] * x[2] + 0.0625 * x[1] - 0.0625 * x[2] * x[2] + 0.0625;
+  return 0;
+}
+static void HdivBasisHex(CeedInt Q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) {
+  // Get 1D quadrature on [-1,1]
+  CeedScalar q_ref_1d[Q], q_weight_1d[Q];
+  switch (quad_mode) {
+    case CEED_GAUSS:
+      CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d);
+      break;
+    case CEED_GAUSS_LOBATTO:
+      CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d);
+      break;
+  }
+
+  // Divergence operator; Divergence of nodal basis for ref element
+  CeedScalar D = 0.125;
+  // Loop over quadrature points
+  CeedScalar Bx[24], By[24], Bz[24];
+  CeedScalar x[3];
+  for (CeedInt k = 0; k < Q; k++) {
+    for (CeedInt i = 0; i < Q; i++) {
+      for (CeedInt j = 0; j < Q; j++) {
+        CeedInt k1            = Q * Q * k + Q * i + j;
+        q_ref[k1 + 0 * Q * Q] = q_ref_1d[j];
+        q_ref[k1 + 1 * Q * Q] = q_ref_1d[i];
+        q_ref[k1 + 2 * Q * Q] = q_ref_1d[k];
+        q_weights[k1]         = q_weight_1d[j] * q_weight_1d[i] * q_weight_1d[k];
+        x[0]                  = q_ref_1d[j];
+        x[1]                  = q_ref_1d[i];
+        x[2]                  = q_ref_1d[k];
+        NodalHdivBasisHex(x, Bx, By, Bz);
+        for (CeedInt d = 0; d < 24; d++) {
+          interp[k1 * 24 + d]                  = Bx[d];
+          interp[k1 * 24 + d + 24 * Q * Q * Q] = By[d];
+          interp[k1 * 24 + d + 48 * Q * Q * Q] = Bz[d];
+          div[k1 * 24 + d]                     = D;
+        }
+      }
+    }
+  }
+}

examples/Hdiv-mass/basis/Hdiv-quad.h

@@ -0,0 +1,82 @@
+// Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
+// Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
+// All Rights reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+// Hdiv basis for quadrilateral element in 2D
+// Local numbering is as follow (each edge has 2 vector dof)
+//     b4     b5
+//    2---------3
+//  b7|         |b3
+//    |         |
+//  b6|         |b2
+//    0---------1
+//     b0     b1
+// Bx[0-->7] = b0_x-->b7_x, By[0-->7] = b0_y-->b7_y
+// To see how the nodal basis is constructed visit:
+// https://github.com/rezgarshakeri/H-div-Tests
+int NodalHdivBasisQuad(CeedScalar *x, CeedScalar *Bx, CeedScalar *By) {
+  Bx[0] = 0.125 * x[0] * x[0] - 0.125;
+  By[0] = -0.25 * x[0] * x[1] + 0.25 * x[0] + 0.25 * x[1] - 0.25;
+  Bx[1] = 0.125 - 0.125 * x[0] * x[0];
+  By[1] = 0.25 * x[0] * x[1] - 0.25 * x[0] + 0.25 * x[1] - 0.25;
+  Bx[2] = -0.25 * x[0] * x[1] + 0.25 * x[0] - 0.25 * x[1] + 0.25;
+  By[2] = 0.125 * x[1] * x[1] - 0.125;
+  Bx[3] = 0.25 * x[0] * x[1] + 0.25 * x[0] + 0.25 * x[1] + 0.25;
+  By[3] = 0.125 - 0.125 * x[1] * x[1];
+  Bx[4] = 0.125 * x[0] * x[0] - 0.125;
+  By[4] = -0.25 * x[0] * x[1] - 0.25 * x[0] + 0.25 * x[1] + 0.25;
+  Bx[5] = 0.125 - 0.125 * x[0] * x[0];
+  By[5] = 0.25 * x[0] * x[1] + 0.25 * x[0] + 0.25 * x[1] + 0.25;
+  Bx[6] = -0.25 * x[0] * x[1] + 0.25 * x[0] + 0.25 * x[1] - 0.25;
+  By[6] = 0.125 * x[1] * x[1] - 0.125;
+  Bx[7] = 0.25 * x[0] * x[1] + 0.25 * x[0] - 0.25 * x[1] - 0.25;
+  By[7] = 0.125 - 0.125 * x[1] * x[1];
+  return 0;
+}
+static void HdivBasisQuad(CeedInt Q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedScalar *div, CeedQuadMode quad_mode) {
+  // Get 1D quadrature on [-1,1]
+  CeedScalar q_ref_1d[Q], q_weight_1d[Q];
+  switch (quad_mode) {
+    case CEED_GAUSS:
+      CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d);
+      break;
+    case CEED_GAUSS_LOBATTO:
+      CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d);
+      break;
+  }
+
+  // Divergence operator; Divergence of nodal basis for ref element
+  CeedScalar D = 0.25;
+  // Loop over quadrature points
+  CeedScalar Bx[8], By[8];
+  CeedScalar x[2];
+
+  for (CeedInt i = 0; i < Q; i++) {
+    for (CeedInt j = 0; j < Q; j++) {
+      CeedInt k1        = Q * i + j;
+      q_ref[k1]         = q_ref_1d[j];
+      q_ref[k1 + Q * Q] = q_ref_1d[i];
+      q_weights[k1]     = q_weight_1d[j] * q_weight_1d[i];
+      x[0]              = q_ref_1d[j];
+      x[1]              = q_ref_1d[i];
+      NodalHdivBasisQuad(x, Bx, By);
+      for (CeedInt k = 0; k < 8; k++) {
+        interp[k1 * 8 + k]             = Bx[k];
+        interp[k1 * 8 + k + 8 * Q * Q] = By[k];
+        div[k1 * 8 + k]                = D;
+      }
+    }
+  }
+}

examples/Hdiv-mass/basis/L2-P0.h

@@ -0,0 +1,58 @@
+// Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
+// Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
+// All Rights reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+// Build L2 constant basis
+
+static void L2BasisP0(CeedInt dim, CeedInt Q, CeedScalar *q_ref, CeedScalar *q_weights, CeedScalar *interp, CeedQuadMode quad_mode) {
+  // Get 1D quadrature on [-1,1]
+  CeedScalar q_ref_1d[Q], q_weight_1d[Q];
+  switch (quad_mode) {
+    case CEED_GAUSS:
+      CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d);
+      break;
+    case CEED_GAUSS_LOBATTO:
+      CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d);
+      break;
+  }
+
+  // P0 L2 basis is just a constant
+  CeedScalar P0 = 1.0;
+  // Loop over quadrature points
+  if (dim == 2) {
+    for (CeedInt i = 0; i < Q; i++) {
+      for (CeedInt j = 0; j < Q; j++) {
+        CeedInt k1        = Q * i + j;
+        q_ref[k1]         = q_ref_1d[j];
+        q_ref[k1 + Q * Q] = q_ref_1d[i];
+        q_weights[k1]     = q_weight_1d[j] * q_weight_1d[i];
+        interp[k1]        = P0;
+      }
+    }
+  } else {
+    for (CeedInt k = 0; k < Q; k++) {
+      for (CeedInt i = 0; i < Q; i++) {
+        for (CeedInt j = 0; j < Q; j++) {
+          CeedInt k1            = Q * Q * k + Q * i + j;
+          q_ref[k1 + 0 * Q * Q] = q_ref_1d[j];
+          q_ref[k1 + 1 * Q * Q] = q_ref_1d[i];
+          q_ref[k1 + 2 * Q * Q] = q_ref_1d[k];
+          q_weights[k1]         = q_weight_1d[j] * q_weight_1d[i] * q_weight_1d[k];
+          interp[k1]            = P0;
+        }
+      }
+    }
+  }
+}