Add DG vector diffusion example

bugfix for 4935

.github/workflows/builds-and-tests.yml

@@ -132,12 +132,14 @@ jobs:
             hypre-target: int32
             precision: fp64
             enzyme: true
-            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme)
+            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme) LDFLAGS=-L$LLVM_PREFIX/lib/c++
 
     name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}
 
     runs-on: ${{ matrix.os }}
 
+    continue-on-error: ${{ matrix.enzyme && true || false }}
+
     steps:
       # Fix 'No space left on device' errors for Ubuntu builds.
       - name: Run Actions Cleaner
@@ -292,10 +294,12 @@ jobs:
         run: |
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
-          brew install llvm@20 enzyme
-          echo "LLVM_PREFIX=$(brew --prefix llvm@20)" >> $GITHUB_ENV
-          echo "OMPI_CC=$(brew --prefix llvm@20)/bin/clang" >> $GITHUB_ENV
-          echo "OMPI_CXX=$(brew --prefix llvm@20)/bin/clang++" >> $GITHUB_ENV
+          brew install enzyme
+          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
+          LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
+          echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
+          echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
+          echo "OMPI_CXX=$LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
 
       # MFEM build and test
       - name: build

.gitignore

@@ -63,6 +63,7 @@ examples/ex3[0-9]
 examples/ex3[0-9]p
 examples/ex4[0-9]
 examples/ex4[0-9]p
+examples/vector-dg-diffusion
 
 examples/refined.mesh
 examples/displaced.mesh

examples/vector-dg-diffusion.cpp

@@ -0,0 +1,164 @@
+#include "mfem.hpp"
+#include "vector-dg-diffusion.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class RepeatedCoefficient : public VectorCoefficient
+{
+   Coefficient &coeff;
+public:
+   RepeatedCoefficient(int dim, Coefficient &coeff_)
+      : VectorCoefficient(dim), coeff(coeff_)
+   { }
+   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      V.SetSize(vdim);
+      V = coeff.Eval(T, ip);
+   }
+};
+
+real_t u_fn(const Vector &xvec);
+real_t f_fn(const Vector &xvec);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int ref_levels = 0;
+   int order = 1;
+   real_t sigma = -1.0;
+   real_t kappa = -1.0;
+   const char *device_config = "cpu";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the three DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the three DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.ParseCheck();
+
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+
+   Device device(device_config);
+   device.Print();
+
+   Mesh mesh(mesh_file);
+   const int dim = mesh.Dimension();
+
+   {
+      if (ref_levels < 0)
+      {
+         ref_levels = (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
+      }
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec, dim);
+   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
+
+
+   FunctionCoefficient scalar_f_coeff(f_fn);
+   FunctionCoefficient scalar_u_coeff(u_fn);
+   RepeatedCoefficient f_coeff(dim, scalar_f_coeff);
+   RepeatedCoefficient u_coeff(dim, scalar_u_coeff);
+
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(5.0);
+   RepeatedCoefficient zero_vec(dim, zero);
+
+   LinearForm b(&fespace);
+   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f_coeff));
+   b.AddBdrFaceIntegrator(
+      new VectorDGDirichletLFIntegrator(u_coeff, one, sigma, kappa));
+   b.Assemble();
+
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new VectorDiffusionIntegrator(one));
+   a.AddInteriorFaceIntegrator(new VectorDGDiffusionIntegrator(
+                                  one, sigma, kappa, dim));
+   a.AddBdrFaceIntegrator(new VectorDGDiffusionIntegrator(
+                             one, sigma, kappa, dim));
+   a.Assemble();
+   a.Finalize();
+
+   const SparseMatrix &A = a.SpMat();
+#ifndef MFEM_USE_SUITESPARSE
+   GSSmoother M(A);
+   if (sigma == -1.0)
+   {
+      PCG(A, M, b, x, 1, 500, 1e-12, 0.0);
+   }
+   else
+   {
+      GMRES(A, M, b, x, 1, 500, 10, 1e-12, 0.0);
+   }
+#else
+   UMFPackSolver umf_solver;
+   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+   umf_solver.SetOperator(A);
+   umf_solver.Mult(b, x);
+#endif
+
+   ParaViewDataCollection pv("DGDiffusion", &mesh);
+   pv.SetPrefixPath("ParaView");
+   pv.SetHighOrderOutput(true);
+   pv.SetLevelsOfDetail(order);
+   pv.RegisterField("u", &x);
+   pv.SetCycle(0);
+   pv.SetTime(0.0);
+   pv.Save();
+
+   cout << "L2 error: " << x.ComputeL2Error(u_coeff) << '\n';
+
+   return 0;
+}
+
+constexpr real_t pi = M_PI;
+constexpr real_t pi2 = pi*pi;
+
+real_t u_fn(const Vector &xvec)
+{
+   int dim = xvec.Size();
+   real_t x = pi*xvec[0], y = pi*xvec[1];
+   if (dim == 2) { return sin(x)*sin(y); }
+   else { real_t z = pi*xvec[2]; return sin(x)*sin(y)*sin(z); }
+}
+
+real_t f_fn(const Vector &xvec)
+{
+   int dim = xvec.Size();
+   real_t x = pi*xvec[0], y = pi*xvec[1];
+
+   if (dim == 2)
+   {
+      return 2*pi2*sin(x)*sin(y);
+   }
+   else // dim == 3
+   {
+      real_t z = pi*xvec[2];
+      return 3*pi2*sin(x)*sin(y)*sin(z);
+   }
+}

examples/vector-dg-diffusion.hpp

@@ -0,0 +1,433 @@
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+class VectorDGDiffusionIntegrator : public BilinearFormIntegrator
+{
+protected:
+   Coefficient *Q = nullptr;
+   MatrixCoefficient *MQ = nullptr;
+   real_t sigma, kappa;
+   int vdim;
+
+   // these are not thread-safe!
+   Vector shape1, shape2, dshape1dn, dshape2dn, nor, nh, ni;
+   DenseMatrix jmat, dshape1, dshape2, mq, adjJ;
+
+public:
+   VectorDGDiffusionIntegrator(real_t s, real_t k, int vd=-1)
+      : sigma(s), kappa(k), vdim(vd) { }
+   VectorDGDiffusionIntegrator(Coefficient &q, real_t s, real_t k, int vd=-1)
+      : Q(&q), sigma(s), kappa(k), vdim(vd) { }
+   VectorDGDiffusionIntegrator(MatrixCoefficient &mq, real_t s, real_t k,
+                               int vd=-1)
+      : MQ(&mq), sigma(s), kappa(k), vdim(vd) { }
+
+   using BilinearFormIntegrator::AssembleFaceMatrix;
+   virtual void AssembleFaceMatrix(const FiniteElement &el1,
+                                   const FiniteElement &el2,
+                                   FaceElementTransformations &Trans,
+                                   DenseMatrix &full_elmat);
+};
+
+class VectorDGDirichletLFIntegrator : public LinearFormIntegrator
+{
+protected:
+   VectorCoefficient &uD;
+   Coefficient *Q = nullptr;
+   MatrixCoefficient *MQ = nullptr;
+   real_t sigma, kappa;
+   int vdim;
+
+   // these are not thread-safe!
+   Vector shape, dshape_dn, nor, nh, ni, uD_vec;
+   DenseMatrix dshape, mq, adjJ;
+
+public:
+   VectorDGDirichletLFIntegrator(VectorCoefficient &u, real_t s, real_t k,
+                                 int vd=-1)
+      : uD(u), sigma(s), kappa(k), vdim(vd) { }
+   VectorDGDirichletLFIntegrator(VectorCoefficient &u, Coefficient &q, real_t s,
+                                 real_t k, int vd=-1)
+      : uD(u), Q(&q), sigma(s), kappa(k), vdim(vd) { }
+   VectorDGDirichletLFIntegrator(VectorCoefficient &u, MatrixCoefficient &mq,
+                                 real_t s, real_t k, int vd=-1)
+      : uD(u), MQ(&mq), sigma(s), kappa(k), vdim(vd) { }
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               ElementTransformation &Tr,
+                               Vector &elvect) override
+   { MFEM_ABORT("Not implemented."); }
+
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               FaceElementTransformations &Tr,
+                               Vector &elvect) override;
+};
+
+void VectorDGDiffusionIntegrator::AssembleFaceMatrix(
+   const FiniteElement &el1, const FiniteElement &el2,
+   FaceElementTransformations &Trans, DenseMatrix &full_elmat)
+{
+   int dim, ndof1, ndof2, ndofs;
+   bool kappa_is_nonzero = (kappa != 0.);
+   real_t w, wq = 0.0;
+
+   const int sdim = Trans.GetSpaceDim();
+   if (vdim < 0) { vdim = sdim; }
+
+   dim = el1.GetDim();
+   ndof1 = el1.GetDof();
+
+   nor.SetSize(dim);
+   nh.SetSize(dim);
+   ni.SetSize(dim);
+   adjJ.SetSize(dim);
+   if (MQ)
+   {
+      mq.SetSize(dim);
+   }
+
+   shape1.SetSize(ndof1);
+   dshape1.SetSize(ndof1, dim);
+   dshape1dn.SetSize(ndof1);
+   if (Trans.Elem2No >= 0)
+   {
+      ndof2 = el2.GetDof();
+      shape2.SetSize(ndof2);
+      dshape2.SetSize(ndof2, dim);
+      dshape2dn.SetSize(ndof2);
+   }
+   else
+   {
+      ndof2 = 0;
+   }
+
+   ndofs = ndof1 + ndof2;
+   DenseMatrix elmat;
+   elmat.SetSize(ndofs);
+   elmat = 0.0;
+   if (kappa_is_nonzero)
+   {
+      jmat.SetSize(ndofs);
+      jmat = 0.;
+   }
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == nullptr)
+   {
+      // a simple choice for the integration order
+      int order;
+      if (ndof2)
+      {
+         order = 2 * std::max(el1.GetOrder(), el2.GetOrder());
+      }
+      else
+      {
+         order = 2 * el1.GetOrder();
+      }
+      ir = &IntRules.Get(Trans.GetGeometryType(), order);
+   }
+
+   // assemble: < {(Q \nabla u).n},[v] >      --> elmat
+   //           kappa < {h^{-1} Q} [u],[v] >  --> jmat
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Access the neighboring elements' integration points
+      // Note: eip2 will only contain valid data if Elem2 exists
+      const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
+      const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
+
+      if (dim == 1)
+      {
+         nor(0) = 2 * eip1.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Trans.Jacobian(), nor);
+      }
+
+      el1.CalcShape(eip1, shape1);
+      el1.CalcDShape(eip1, dshape1);
+      w = ip.weight / Trans.Elem1->Weight();
+      if (ndof2)
+      {
+         w /= 2;
+      }
+      if (!MQ)
+      {
+         if (Q)
+         {
+            w *= Q->Eval(*Trans.Elem1, eip1);
+         }
+         ni.Set(w, nor);
+      }
+      else
+      {
+         nh.Set(w, nor);
+         MQ->Eval(mq, *Trans.Elem1, eip1);
+         mq.MultTranspose(nh, ni);
+      }
+      CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
+      adjJ.Mult(ni, nh);
+      if (kappa_is_nonzero)
+      {
+         wq = ni * nor;
+      }
+
+      // Note: in the jump term, we use 1/h1 = |nor|/det(J1) which is
+      // independent of Loc1 and always gives the size of element 1 in
+      // direction perpendicular to the face. Indeed, for linear transformation
+      //
+      //      |nor|=measure(face)/measure(ref. face),
+      //
+      //      det(J1)=measure(element)/measure(ref. element),
+      //
+      //      and the ratios measure(ref. element)/measure(ref. face)
+      //      are compatible for all element/face pairs.
+      //
+      // For example: meas(ref. tetrahedron)/meas(ref. triangle) = 1/3, and
+      // for any tetrahedron vol(tet)=(1/3)*height*area(base).
+      //
+      // For interior faces: q_e/h_e=(q1/h1+q2/h2)/2.
+
+      dshape1.Mult(nh, dshape1dn);
+      for (int i = 0; i < ndof1; i++)
+         for (int j = 0; j < ndof1; j++)
+         {
+            elmat(i, j) += shape1(i) * dshape1dn(j);
+         }
+
+      if (ndof2)
+      {
+         el2.CalcShape(eip2, shape2);
+         el2.CalcDShape(eip2, dshape2);
+         w = ip.weight / 2 / Trans.Elem2->Weight();
+         if (!MQ)
+         {
+            if (Q)
+            {
+               w *= Q->Eval(*Trans.Elem2, eip2);
+            }
+            ni.Set(w, nor);
+         }
+         else
+         {
+            nh.Set(w, nor);
+            MQ->Eval(mq, *Trans.Elem2, eip2);
+            mq.MultTranspose(nh, ni);
+         }
+         CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
+         adjJ.Mult(ni, nh);
+         if (kappa_is_nonzero)
+         {
+            wq += ni * nor;
+         }
+
+         dshape2.Mult(nh, dshape2dn);
+
+         for (int i = 0; i < ndof1; i++)
+            for (int j = 0; j < ndof2; j++)
+            {
+               elmat(i, ndof1 + j) += shape1(i) * dshape2dn(j);
+            }
+
+         for (int i = 0; i < ndof2; i++)
+            for (int j = 0; j < ndof1; j++)
+            {
+               elmat(ndof1 + i, j) -= shape2(i) * dshape1dn(j);
+            }
+
+         for (int i = 0; i < ndof2; i++)
+            for (int j = 0; j < ndof2; j++)
+            {
+               elmat(ndof1 + i, ndof1 + j) -= shape2(i) * dshape2dn(j);
+            }
+      }
+
+      if (kappa_is_nonzero)
+      {
+         // only assemble the lower triangular part of jmat
+         wq *= kappa;
+         for (int i = 0; i < ndof1; i++)
+         {
+            const real_t wsi = wq * shape1(i);
+            for (int j = 0; j <= i; j++)
+            {
+               jmat(i, j) += wsi * shape1(j);
+            }
+         }
+         if (ndof2)
+         {
+            for (int i = 0; i < ndof2; i++)
+            {
+               const int i2 = ndof1 + i;
+               const real_t wsi = wq * shape2(i);
+               for (int j = 0; j < ndof1; j++)
+               {
+                  jmat(i2, j) -= wsi * shape1(j);
+               }
+               for (int j = 0; j <= i; j++)
+               {
+                  jmat(i2, ndof1 + j) += wsi * shape2(j);
+               }
+            }
+         }
+      }
+   }
+
+   // elmat := -elmat + sigma*elmat^t + jmat
+   if (kappa_is_nonzero)
+   {
+      for (int i = 0; i < ndofs; i++)
+      {
+         for (int j = 0; j < i; j++)
+         {
+            real_t aij = elmat(i, j), aji = elmat(j, i), mij = jmat(i, j);
+            elmat(i, j) = sigma * aji - aij + mij;
+            elmat(j, i) = sigma * aij - aji + mij;
+         }
+         elmat(i, i) = (sigma - 1.) * elmat(i, i) + jmat(i, i);
+      }
+   }
+   else
+   {
+      for (int i = 0; i < ndofs; i++)
+      {
+         for (int j = 0; j < i; j++)
+         {
+            real_t aij = elmat(i, j), aji = elmat(j, i);
+            elmat(i, j) = sigma * aji - aij;
+            elmat(j, i) = sigma * aij - aji;
+         }
+         elmat(i, i) *= (sigma - 1.);
+      }
+   }
+
+   // populate full matrix following github issue #2909
+   full_elmat.SetSize(vdim*(ndof1 + ndof2));
+   full_elmat = 0.0;
+   for (int d=0; d<vdim; ++d)
+   {
+      for (int j=0; j<ndofs; ++j)
+      {
+         int jj = (j < ndof1) ? j + d*ndof1 : j - ndof1 + d*ndof2 + vdim*ndof1;
+         for (int i=0; i<ndofs; ++i)
+         {
+            int ii = (i < ndof1) ? i + d*ndof1 : i - ndof1 + d*ndof2 + vdim*ndof1;
+            full_elmat(ii, jj) += elmat(i, j);
+         }
+      }
+   }
+};
+
+void VectorDGDirichletLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
+{
+   const int dim = el.GetDim();
+   const int sdim = Tr.GetSpaceDim();
+
+   if (vdim < 0) { vdim = sdim; }
+
+   const int ndof = el.GetDof();
+
+   bool kappa_is_nonzero = (kappa != 0.);
+   real_t w;
+
+   nor.SetSize(dim);
+   nh.SetSize(dim);
+   ni.SetSize(dim);
+   adjJ.SetSize(dim);
+   if (MQ)
+   {
+      mq.SetSize(dim);
+   }
+
+   shape.SetSize(ndof);
+   dshape.SetSize(ndof, dim);
+   dshape_dn.SetSize(ndof);
+
+   elvect.SetSize(vdim * ndof);
+   elvect = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      // a simple choice for the integration order; is this OK?
+      int order = 2*el.GetOrder();
+      ir = &IntRules.Get(Tr.GetGeometryType(), order);
+   }
+
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+
+      // Set the integration point in the face and the neighboring element
+      Tr.SetAllIntPoints(&ip);
+
+      // Access the neighboring element's integration point
+      const IntegrationPoint &eip = Tr.GetElement1IntPoint();
+
+      uD.Eval(uD_vec, Tr, ip);
+
+      if (dim == 1)
+      {
+         nor(0) = 2*eip.x - 1.0;
+      }
+      else
+      {
+         CalcOrtho(Tr.Jacobian(), nor);
+      }
+
+      el.CalcShape(eip, shape);
+      el.CalcDShape(eip, dshape);
+
+      // compute uD through the face transformation
+      w = ip.weight / Tr.Elem1->Weight();
+      if (!MQ)
+      {
+         if (Q)
+         {
+            w *= Q->Eval(*Tr.Elem1, eip);
+         }
+         ni.Set(w, nor);
+      }
+      else
+      {
+         nh.Set(w, nor);
+         MQ->Eval(mq, *Tr.Elem1, eip);
+         mq.MultTranspose(nh, ni);
+      }
+      CalcAdjugate(Tr.Elem1->Jacobian(), adjJ);
+      adjJ.Mult(ni, nh);
+
+      dshape.Mult(nh, dshape_dn);
+
+      for (int vd = 0; vd < vdim; ++vd)
+      {
+         for (int i = 0; i < ndof; ++i)
+         {
+            elvect[i + vd*ndof] += sigma * uD_vec[vd] * dshape_dn[i];
+         }
+      }
+      if (kappa_is_nonzero)
+      {
+         for (int vd = 0; vd < vdim; ++vd)
+         {
+            for (int i = 0; i < ndof; ++i)
+            {
+               elvect[i + vd*ndof] += kappa*(ni*nor) * uD_vec[vd] * shape[i];
+            }
+         }
+      }
+   }
+}
+
+} // namespace mfem

fem/fe_coll.cpp

@@ -401,9 +401,6 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
    {
       MFEM_ABORT("unknown FiniteElementCollection: " << name);
    }
-   MFEM_VERIFY(!strcmp(fec->Name(), name), "input name: \"" << name
-               << "\" does not match the created collection name: \""
-               << fec->Name() << '"');
 
    return fec;
 }

fem/fespace.hpp

@@ -922,6 +922,9 @@ public:
    { return mesh->GetBdrElementType(i); }
 
    /// Returns ElementTransformation for the @a i-th element.
+   /// @note The returned pointer references an object owned by the associated
+   /// @a Mesh that will be modified by other calls to `GetElementTransformation`.
+   /// As such, this pointer should @b not be deleted by the caller.
    ElementTransformation *GetElementTransformation(int i) const
    { return mesh->GetElementTransformation(i); }