Remove unnecessary SetCurvature

Fixed a missing GetRajaResource() for hip

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

@@ -142,6 +142,10 @@ jobs:
 
     continue-on-error: ${{ matrix.enzyme && true || false }}
 
+    # Enable ccache for all jobs except Windows (would need sccache).
+    env:
+      USE_CCACHE: ${{ matrix.os != 'windows-latest' }}
+
     steps:
       # Fix 'No space left on device' errors for Ubuntu builds.
       - name: Run Actions Cleaner
@@ -290,6 +294,52 @@ jobs:
           echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
           echo "OMPI_CXX=$LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
 
+      # Restore the compiler cache (ccache). The key embeds the run id, so new
+      # runs save a fresh snapshot; the restore-keys prefix warm-starts from the
+      # most recent prior run (incl. the base branch for PRs).
+      - name: cache ccache
+        if: ${{ env.USE_CCACHE == 'true' }}
+        uses: actions/cache@v5
+        with:
+          path: .ccache
+          key: ccache-${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}-${{ github.run_id }}
+          restore-keys: |
+            ccache-${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}-
+
+      # Configure ccache and select how it is injected into the MFEM build:
+      # - make:  set CXX="ccache g++"; for MPI, OMPI_CXX="ccache g++" so mpicxx
+      #          runs ccache around g++ (not ccache around the mpicxx wrapper).
+      # - cmake: set CMAKE_<LANG>_COMPILER_LAUNCHER=ccache.
+      # - enzyme: wrap the brew clang++ via OMPI_CXX.
+      # The chosen options are passed through build-mfem's 'config-options'
+      # input (see the build step below).
+      - name: configure ccache
+        if: ${{ env.USE_CCACHE == 'true' }}
+        run: |
+          command -v ccache >/dev/null 2>&1 || {
+            if [[ "${{ runner.os }}" == "Linux" ]]; then
+              sudo apt-get update && sudo apt-get install -y ccache
+            else
+              brew install ccache
+            fi
+          }
+          echo "CCACHE_DIR=${{ github.workspace }}/.ccache" >> $GITHUB_ENV
+          echo "CCACHE_MAXSIZE=1G" >> $GITHUB_ENV
+          echo "CCACHE_COMPILERCHECK=content" >> $GITHUB_ENV
+          # Ignore header timestamps (restamped by each checkout) so direct mode hits.
+          echo "CCACHE_SLOPPINESS=include_file_mtime,include_file_ctime,time_macros" >> $GITHUB_ENV
+          # Hash absolute paths relative to the workspace.
+          echo "CCACHE_BASEDIR=${{ github.workspace }}" >> $GITHUB_ENV
+          if [[ "${{ matrix.enzyme }}" == "true" ]]; then
+            echo "OMPI_CXX=ccache $LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
+          elif [[ "${{ matrix.build-system }}" == "cmake" ]]; then
+            echo 'CCACHE_CONFIG_OPTS=-DCMAKE_CXX_COMPILER_LAUNCHER=ccache -DCMAKE_C_COMPILER_LAUNCHER=ccache' >> $GITHUB_ENV
+          else
+            echo "OMPI_CXX=ccache g++" >> $GITHUB_ENV
+            echo 'CCACHE_CONFIG_OPTS=CXX="ccache g++" MPICXX="mpicxx"' >> $GITHUB_ENV
+          fi
+        shell: bash
+
       # MFEM build and test
       - name: build
         uses: mfem/github-actions/build-mfem@v2.7
@@ -305,9 +355,14 @@ jobs:
           metis-dir: ${{ env.METIS_TOP_DIR }}
           mfem-dir: ${{ env.MFEM_TOP_DIR }}
           precision: ${{ matrix.precision }}
-          config-options: ${{ matrix.config-opts }}
+          config-options: ${{ matrix.config-opts }} ${{ env.CCACHE_CONFIG_OPTS }}
           library-only: ${{ matrix.target == 'dbg' && matrix.os != 'ubuntu-latest' }}
 
+      - name: ccache stats
+        if: ${{ env.USE_CCACHE == 'true' }}
+        run: ccache -s
+        shell: bash
+
       # Run checks (and only checks) on debug targets
       - name: checks
         if: matrix.build-system == 'make' && matrix.target == 'dbg'

.github/workflows/cleanup-pr-caches.yml

@@ -0,0 +1,42 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+---
+# A closed PR's caches can never be restored again, so delete them to free
+# space against the 10 GB per-repo cache limit.
+name: Cleanup PR caches
+
+on:
+  pull_request:
+    types: [closed]
+
+permissions:
+  actions: write
+
+jobs:
+  cleanup:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Delete caches for the closed PR
+        env:
+          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          GH_REPO: ${{ github.repository }}
+          PR_REF: refs/pull/${{ github.event.pull_request.number }}/merge
+        run: |
+          echo "Deleting caches for $PR_REF"
+          while :; do
+            ids=$(gh cache list --ref "$PR_REF" --limit 100 --json id --jq '.[].id')
+            [ -n "$ids" ] || break
+            echo "$ids" | while read -r id; do
+              [ -n "$id" ] || continue
+              echo "Deleting cache $id"
+              gh cache delete "$id" || echo "  (already gone)"
+            done
+          done

.github/workflows/repo-check.yml

@@ -13,6 +13,7 @@ name: "Checks"
 
 permissions:
   actions: write
+  pull-requests: read
 
 on:
   push:
@@ -29,6 +30,11 @@ concurrency:
 # by checking if the workflow trigger is 'push' ("github.event_name == 'push'")
 # and if we are in a fork ("github.event.pull_request.head.repo.full_name !=
 # github.repository").
+#
+# The logic for the branch-history check is slightly different, since that check
+# also inspects the PR's labels to allow for overriding failures. In this case,
+# we run on all 'pull_request' triggers, but only run for 'push' triggers that
+# do not correspond to any open PRs.
 
 jobs:
   file-headers-check:
@@ -128,10 +134,7 @@ jobs:
 
   branch-history:
     if: |
-      github.ref != 'refs/heads/next' &&
-      github.ref != 'refs/heads/master' &&
-      (github.event_name == 'push' ||
-       github.event.pull_request.head.repo.full_name != github.repository)
+      github.ref != 'refs/heads/next' && github.ref != 'refs/heads/master'
     runs-on: ubuntu-latest
     steps:
       - name: checkout mfem
@@ -139,7 +142,27 @@ jobs:
         with:
           fetch-depth: 0
 
+      - name: check for pull request
+        id: check_pr
+        if: github.event_name == 'push'
+        env:
+          GH_TOKEN: ${{ github.token }}
+        run: |
+          pr_exists=$(gh pr list --repo "$GITHUB_REPOSITORY" \
+            --head "$GITHUB_REF_NAME" \
+            --state open \
+            --json number \
+            --jq 'length > 0')
+          echo "pr_exists=$pr_exists" >> "$GITHUB_OUTPUT"
+
       - name: branch-history
+        id: branch_history
+        if: |
+          (github.event_name == 'pull_request' ||
+           github.event_name == 'workflow_dispatch' ||
+           steps.check_pr.outputs.pr_exists == 'false')
+        continue-on-error: ${{ contains(github.event.pull_request.labels.*.name,
+                                        'branch-history-override') }}
         run: |
           # We override origin to make sure we point to the main repo.
           # This is to have consistent test results on PRs from forks.
@@ -147,3 +170,9 @@ jobs:
           git remote add origin https://github.com/mfem/mfem.git
           git checkout -b gh-actions-branch-history
           ./config/githooks/pre-push --history
+
+      - name: report branch-history override
+        if: steps.branch_history.outcome == 'failure'
+        run: |
+          echo "::warning::branch-history check failed, but the" \
+            "'branch-history-override' label is set."

config/githooks/README.md

@@ -39,3 +39,8 @@ when a picture was added for documentation.
 If that is the case, make sure the failure is indeed justified, and rerun the
 push command with the `--no-verify` option. This will skip the hooks, allowing
 you to push those changes.
+
+The `branch-history` check is run automatically through GitHub Actions. If a
+branch is known to have a large number of changes that are legitimate, the
+check can be overridden by setting  the label 'branch-history-override' on the
+pull request.

doc/CodeDocumentation.dox

@@ -119,6 +119,8 @@ namespace mfem {
  * - <a class="el" href="ex40p_8cpp_source.html">Example 40p</a>: parallel eikonal equation
  * - <a class="el" href="ex41_8cpp_source.html">Example 41</a>: DG/CG IMEX time-dependent advection-diffusion
  * - <a class="el" href="ex41p_8cpp_source.html">Example 41p</a>: parallel DG/CG IMEX time-dependent advection-diffusion
+ * - <a class="el" href="ex42_8cpp_source.html">Example 42</a>: clamped biharmonic equation
+ * - <a class="el" href="ex42p_8cpp_source.html">Example 42p</a>: parallel clamped biharmonic equation
  *
  * <H4>AmgX Examples</H4>
  * - Variants of Examples

examples/CMakeLists.txt

@@ -47,6 +47,7 @@ list(APPEND ALL_EXE_SRCS
   ex39.cpp
   ex40.cpp
   ex41.cpp
+  ex42.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -91,6 +92,7 @@ if (MFEM_USE_MPI)
     ex39p.cpp
     ex40p.cpp
     ex41p.cpp
+    ex42p.cpp
   )
 endif()
 

examples/ex42.cpp

@@ -0,0 +1,311 @@
+//                                MFEM Example 42
+//
+// Compile with: make ex42p
+//
+// Sample runs: ex42 -r 3
+//              ex42 -m ../data/hexagon.mesh -r 3 -o 3
+//              ex42 -m ../data/square-mixed.mesh -r 2 -eta 5
+//              ex42 -m ../data/l-shape.mesh -r 3
+//
+// Description: This example solves the clamped biharmonic equation,
+//
+//                      ∇⁴u = f in Ω,  u = 0 and ∇u⋅n = 0 on ∂Ω,
+//
+//              in 2D using just H¹-conforming finite elements by employing the interior penalty
+//              method outlined in [1]. This example demonstrates an approach to solving higher-order
+//              PDEs in MFEM and implementation of custom domain and face integrators to solve the
+//              weak form
+//
+//                      (H(u), H(v))_D - <{{n^T⋅H(u)⋅n}}, [[∇v⋅n]]>_F
+//                                     - <{{n^T⋅H(v)⋅n}}, [[∇u⋅n]]>_F
+//                                     + (η/h_e)<[[∇u⋅n]], [[∇v⋅n]]>_F = (f,v)_D ,
+//
+//               where (⋅,⋅)_D is domain integration, <⋅,⋅>_F is face
+//               integration, and H(⋅) is the Hessian.
+//
+//              [1] Brenner, Susanne & Sung, Li-yeng. (2005). C0 Interior Penalty Methods
+//                  for Fourth Order Elliptic Boundary Value Problems on Polygonal Domains.
+//                  Journal of Scientific Computing. 22-23. 83-118. 10.1007/s10915-004-4135-7.
+
+#include <mfem.hpp>
+
+using namespace mfem;
+using namespace std;
+
+class BiharmonicIntegrator : public BilinearFormIntegrator
+{
+private:
+   Coefficient &D;
+
+   inline static const Vector factors_2D{1.0, 2.0, 1.0};
+   mutable DenseMatrix hessian;
+   mutable Vector factors;
+public:
+   BiharmonicIntegrator(Coefficient &D_) : D(D_) {}
+
+   void AssembleElementMatrix(const FiniteElement &el,
+                              ElementTransformation &Trans, DenseMatrix &elmat) override;
+};
+
+class C0InteriorPenaltyIntegrator : public BilinearFormIntegrator
+{
+private:
+   const double eta;
+
+   mutable Vector normal[2], dnshape[2], nv[2], nd2nshape[2];
+   mutable DenseMatrix dshape[2], hessian[2], blockJ[2][2], blockC[2][2], elmatJ_p,
+           elmatC_p;
+public:
+   C0InteriorPenaltyIntegrator(double eta_) : eta(eta_) {};
+
+   void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2,
+                           FaceElementTransformations &Trans, DenseMatrix &elmat) override;
+};
+
+int main(int argc, char *argv[])
+{
+   // Parse command line args
+   const char *mesh_file = "../data/star.mesh";
+   int order = 2;
+   int ref_levels = 0;
+   real_t eta = 10;
+   int max_it = 10000;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&eta, "-eta", "--penalty-coeff",
+                  "Penalty coefficient.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // Read the mesh file
+   Mesh mesh(mesh_file);
+   int dim = mesh.Dimension();
+   MFEM_ASSERT(dim == 2, "This example only supports 2D meshes.");
+
+   // Refine the mesh
+   for (int i = 0; i < ref_levels; i++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // Initialize the FE collection and FiniteElementSpace
+   H1_FECollection fe_coll(order, dim);
+   FiniteElementSpace fespace(&mesh, &fe_coll, 1);
+
+   // Get the degrees-of-freedom (DOFs) associated with the sides of the panel
+   Array<int> all_bdr_marker(mesh.bdr_attributes.Size());
+   all_bdr_marker = 1; // Mark all sides
+   Array<int> ess_tdof_list;
+   fespace.GetEssentialTrueDofs(all_bdr_marker, ess_tdof_list);
+
+   ConstantCoefficient one(1.0);
+
+   // Initialize the bilinear form
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new BiharmonicIntegrator(one));
+   a.AddInteriorFaceIntegrator(new C0InteriorPenaltyIntegrator(eta));
+   a.AddBdrFaceIntegrator(new C0InteriorPenaltyIntegrator(eta));
+   a.Assemble();
+
+   // Initialize the linear form f=1.0
+   LinearForm b(&fespace);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // Form the linear system
+   GridFunction x(&fespace);
+   x = 0.0; // initial guess
+   SparseMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   // Solve the system using CG with symmetric Gauss-Seidel preconditioner
+   GSSmoother M(A);
+   PCG(A, M, B, X, 1, max_it, 1e-12, 0.0);
+
+   // Recover solution and visualize
+   a.RecoverFEMSolution(X, B, x);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport    = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << mesh << x << flush;
+   }
+
+   return 0;
+}
+
+
+void BiharmonicIntegrator::AssembleElementMatrix(const FiniteElement &el,
+                                                 ElementTransformation &Trans, DenseMatrix &elmat)
+{
+   int ndof = el.GetDof();
+   int dim = el.GetDim();
+
+   MFEM_ASSERT(dim == 2, "Dimension must be 2.");
+
+   hessian.SetSize(ndof, dim * (dim + 1) / 2);
+   elmat.SetSize(ndof);
+   factors.SetSize(dim * (dim + 1) / 2);
+
+   elmat = 0.0;
+
+   const IntegrationRule *ir = GetIntegrationRule(el, Trans);
+   if (ir == NULL)
+   {
+      int order = 2*el.GetOrder();
+      ir = &IntRules.Get(el.GetGeomType(), order);
+   }
+
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint(&ip);
+
+      el.CalcPhysHessian(Trans, hessian);
+
+      factors = factors_2D;
+      factors *= D.Eval(Trans, ip) * ip.weight * Trans.Weight();
+
+      AddMultADAt(hessian, factors, elmat);
+   }
+}
+
+void C0InteriorPenaltyIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
+                                                     const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
+{
+   int dim = el1.GetDim();
+   MFEM_ASSERT(dim == 2, "Dimension must be 2.");
+
+   int ndof[2] = {el1.GetDof(), 0};
+   int num_elems = 1;
+   if (Trans.Elem2No >= 0)
+   {
+      ndof[1] = el2.GetDof();
+      num_elems++;
+   }
+
+   for (int i = 0; i < num_elems; i++)
+   {
+      normal[i].SetSize(dim);
+      dshape[i].SetSize(ndof[i], dim);
+      hessian[i].SetSize(ndof[i], dim * (dim + 1) / 2);
+      nv[i].SetSize(dim * (dim + 1) / 2);
+      dnshape[i].SetSize(ndof[i]);
+      nd2nshape[i].SetSize(ndof[i]);
+   }
+
+   for (int i = 0; i < num_elems; i++)
+   {
+      for (int j = 0; j < num_elems; j++)
+      {
+         blockJ[i][j].SetSize(ndof[i], ndof[j]);
+         blockC[i][j].SetSize(ndof[i], ndof[j]);
+      }
+   }
+
+   elmatJ_p.SetSize(ndof[0] + ndof[1]);
+   elmatC_p.SetSize(ndof[0] + ndof[1]);
+   elmat.SetSize(ndof[0] + ndof[1]);
+   elmat = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order = 2 * max(el1.GetOrder(), ndof[1] ? el2.GetOrder() : 0);
+      ir = &IntRules.Get(Trans.GetGeometryType(), order);
+   }
+
+   // Compute edge length
+   double h_e  = 0.0;
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+      Trans.SetAllIntPoints(&ip);
+      h_e += ip.weight * Trans.Weight();
+   }
+
+   const FiniteElement *els[2] = {&el1, &el2};
+   ElementTransformation *el_trans[2] = {Trans.Elem1, Trans.Elem2};
+
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      elmatJ_p = 0.0;
+      elmatC_p = 0.0;
+
+      const IntegrationPoint &ip = ir->IntPoint(p);
+
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Compute normal gradients + Hessians
+      for (int i = 0; i < num_elems; i++)
+      {
+         if (i == 0)
+         {
+            CalcOrtho(Trans.Jacobian(), normal[i]);
+            normal[i] /= normal[i].Norml2();
+         }
+         else
+         {
+            normal[i] = normal[0];
+            normal[i] *= -1;
+         }
+         els[i]->CalcPhysDShape(*el_trans[i], dshape[i]);
+         els[i]->CalcPhysHessian(*el_trans[i], hessian[i]);
+         dshape[i].Mult(normal[i], dnshape[i]);
+         nv[i][0] = normal[i][0]*normal[i][0];
+         nv[i][1] = 2*normal[i][0]*normal[i][1];
+         nv[i][2] = normal[i][1]*normal[i][1];
+         hessian[i].Mult(nv[i], nd2nshape[i]);
+      }
+
+      // Compute blocks
+      for (int i = 0; i < num_elems; i++)
+      {
+         for (int j = 0; j < num_elems; j++)
+         {
+            blockJ[i][j] = 0.0;
+            blockC[i][j] = 0.0;
+            AddMult_a_VWt(-1.0, dnshape[i], nd2nshape[j], blockJ[i][j]);
+            elmatJ_p.SetSubMatrix(i*ndof[0], j*ndof[0], blockJ[i][j]);
+
+            AddMult_a_VWt(eta/h_e, dnshape[i], dnshape[j], blockC[i][j]);
+            elmatC_p.SetSubMatrix(i*ndof[0], j*ndof[0], blockC[i][j]);
+         }
+      }
+
+      // Symmetrize the jump term
+      elmatJ_p.Symmetrize();
+      if (!ndof[1])
+      {
+         elmatJ_p *= 2;
+      }
+
+      // Add penalty term
+      elmatJ_p += elmatC_p;
+      elmatJ_p *= ip.weight * Trans.Weight();
+      elmat += elmatJ_p;
+   }
+}

examples/ex42p.cpp

@@ -0,0 +1,331 @@
+//                       MFEM Example 42 - Parallel Version
+//
+// Compile with: make ex42p
+//
+// Sample runs: mpirun -np 4 ex42p -r 3
+//              mpirun -np 4 ex42p -m ../data/hexagon.mesh -r 3 -o 3
+//              mpirun -np 4 ex42p -m ../data/square-mixed.mesh -r 2 -eta 5
+//              mpirun -np 4 ex42p -m ../data/l-shape.mesh -r 3
+//
+// Description: This example solves the clamped biharmonic equation,
+//
+//                      ∇⁴u = f in Ω,  u = 0 and ∇u⋅n = 0 on ∂Ω,
+//
+//              in 2D using just H¹-conforming finite elements by employing the interior penalty
+//              method outlined in [1]. This example demonstrates an approach to solving higher-order
+//              PDEs in MFEM and implementation of custom domain and face integrators to solve the
+//              weak form
+//
+//                      (H(u), H(v))_D - <{{n^T⋅H(u)⋅n}}, [[∇v⋅n]]>_F
+//                                     - <{{n^T⋅H(v)⋅n}}, [[∇u⋅n]]>_F
+//                                     + (η/h_e)<[[∇u⋅n]], [[∇v⋅n]]>_F = (f,v)_D ,
+//
+//               where (⋅,⋅)_D is domain integration, <⋅,⋅>_F is face
+//               integration, and H(⋅) is the Hessian.
+//
+//              [1] Brenner, Susanne & Sung, Li-yeng. (2005). C0 Interior Penalty Methods
+//                  for Fourth Order Elliptic Boundary Value Problems on Polygonal Domains.
+//                  Journal of Scientific Computing. 22-23. 83-118. 10.1007/s10915-004-4135-7.
+
+#include <mfem.hpp>
+
+using namespace mfem;
+using namespace std;
+
+class BiharmonicIntegrator : public BilinearFormIntegrator
+{
+private:
+   Coefficient &D;
+
+   inline static const Vector factors_2D{1.0, 2.0, 1.0};
+   mutable DenseMatrix hessian;
+   mutable Vector factors;
+public:
+   BiharmonicIntegrator(Coefficient &D_) : D(D_) {}
+
+   void AssembleElementMatrix(const FiniteElement &el,
+                              ElementTransformation &Trans, DenseMatrix &elmat) override;
+};
+
+class C0InteriorPenaltyIntegrator : public BilinearFormIntegrator
+{
+private:
+   const double eta;
+
+   mutable Vector normal[2], dnshape[2], nv[2], nd2nshape[2];
+   mutable DenseMatrix dshape[2], hessian[2], blockJ[2][2], blockC[2][2], elmatJ_p,
+           elmatC_p;
+public:
+   C0InteriorPenaltyIntegrator(double eta_) : eta(eta_) {};
+
+   void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2,
+                           FaceElementTransformations &Trans, DenseMatrix &elmat) override;
+};
+
+int main(int argc, char *argv[])
+{
+   // Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   int size = Mpi::WorldSize();
+   int rank = Mpi::WorldRank();
+   Hypre::Init();
+
+   // Parse command line args
+   const char *mesh_file = "../data/star.mesh";
+   int order = 2;
+   int ref_levels = 0;
+   real_t eta = 10;
+   int max_it = 10000;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&eta, "-eta", "--penalty-coeff",
+                  "Penalty coefficient.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (rank == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // Read the mesh file
+   Mesh mesh(mesh_file);
+   int dim = mesh.Dimension();
+   MFEM_ASSERT(dim == 2, "This example only supports 2D meshes.");
+
+   // Refine the mesh
+   for (int i = 0; i < ref_levels; i++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // Partition the mesh
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // Initialize the FE collection and FiniteElementSpace
+   H1_FECollection fe_coll(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fe_coll, 1);
+
+   // Get the degrees-of-freedom (DOFs) associated with the sides of the panel
+   Array<int> all_bdr_marker(pmesh.bdr_attributes.Size());
+   all_bdr_marker = 1; // Mark all sides
+   Array<int> ess_tdof_list;
+   fespace.GetEssentialTrueDofs(all_bdr_marker, ess_tdof_list);
+
+   ConstantCoefficient one(1.0);
+
+   // Initialize the bilinear form
+   ParBilinearForm a(&fespace);
+   a.AddDomainIntegrator(new BiharmonicIntegrator(one));
+   a.AddInteriorFaceIntegrator(new C0InteriorPenaltyIntegrator(eta));
+   a.AddBdrFaceIntegrator(new C0InteriorPenaltyIntegrator(eta));
+   a.Assemble();
+
+   // Initialize the linear form f=1.0
+   ParLinearForm b(&fespace);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // Form the linear system
+   ParGridFunction x(&fespace);
+   x = 0.0; // initial guess
+   HypreParMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   // Solve the system using CG with hypre's BoomerAMG preconditioner
+   HypreBoomerAMG amg(A);
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(max_it);
+   cg.SetPrintLevel(1);
+   cg.SetPreconditioner(amg);
+   cg.SetOperator(A);
+   cg.Mult(B, X);
+
+   // Recover solution and visualize
+   a.RecoverFEMSolution(X, B, x);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport    = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << size << " " << rank << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+   return 0;
+}
+
+
+void BiharmonicIntegrator::AssembleElementMatrix(const FiniteElement &el,
+                                                 ElementTransformation &Trans, DenseMatrix &elmat)
+{
+   int ndof = el.GetDof();
+   int dim = el.GetDim();
+
+   MFEM_ASSERT(dim == 2, "Dimension must be 2.");
+
+   hessian.SetSize(ndof, dim * (dim + 1) / 2);
+   elmat.SetSize(ndof);
+   factors.SetSize(dim * (dim + 1) / 2);
+
+   elmat = 0.0;
+
+   const IntegrationRule *ir = GetIntegrationRule(el, Trans);
+   if (ir == NULL)
+   {
+      int order = 2*el.GetOrder();
+      ir = &IntRules.Get(el.GetGeomType(), order);
+   }
+
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const mfem::IntegrationPoint &ip = ir->IntPoint(i);
+      Trans.SetIntPoint(&ip);
+
+      el.CalcPhysHessian(Trans, hessian);
+
+      factors = factors_2D;
+      factors *= D.Eval(Trans, ip) * ip.weight * Trans.Weight();
+
+      AddMultADAt(hessian, factors, elmat);
+   }
+}
+
+void C0InteriorPenaltyIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
+                                                     const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
+{
+   int dim = el1.GetDim();
+   MFEM_ASSERT(dim == 2, "Dimension must be 2.");
+
+   int ndof[2] = {el1.GetDof(), 0};
+   int num_elems = 1;
+   if (Trans.Elem2No >= 0)
+   {
+      ndof[1] = el2.GetDof();
+      num_elems++;
+   }
+
+   for (int i = 0; i < num_elems; i++)
+   {
+      normal[i].SetSize(dim);
+      dshape[i].SetSize(ndof[i], dim);
+      hessian[i].SetSize(ndof[i], dim * (dim + 1) / 2);
+      nv[i].SetSize(dim * (dim + 1) / 2);
+      dnshape[i].SetSize(ndof[i]);
+      nd2nshape[i].SetSize(ndof[i]);
+   }
+
+   for (int i = 0; i < num_elems; i++)
+   {
+      for (int j = 0; j < num_elems; j++)
+      {
+         blockJ[i][j].SetSize(ndof[i], ndof[j]);
+         blockC[i][j].SetSize(ndof[i], ndof[j]);
+      }
+   }
+
+   elmatJ_p.SetSize(ndof[0] + ndof[1]);
+   elmatC_p.SetSize(ndof[0] + ndof[1]);
+   elmat.SetSize(ndof[0] + ndof[1]);
+   elmat = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order = 2 * max(el1.GetOrder(), ndof[1] ? el2.GetOrder() : 0);
+      ir = &IntRules.Get(Trans.GetGeometryType(), order);
+   }
+
+   // Compute edge length
+   double h_e  = 0.0;
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(p);
+      Trans.SetAllIntPoints(&ip);
+      h_e += ip.weight * Trans.Weight();
+   }
+
+   const FiniteElement *els[2] = {&el1, &el2};
+   ElementTransformation *el_trans[2] = {Trans.Elem1, Trans.Elem2};
+
+   for (int p = 0; p < ir->GetNPoints(); p++)
+   {
+      elmatJ_p = 0.0;
+      elmatC_p = 0.0;
+
+      const IntegrationPoint &ip = ir->IntPoint(p);
+
+      // Set the integration point in the face and the neighboring elements
+      Trans.SetAllIntPoints(&ip);
+
+      // Compute normal gradients + Hessians
+      for (int i = 0; i < num_elems; i++)
+      {
+         if (i == 0)
+         {
+            CalcOrtho(Trans.Jacobian(), normal[i]);
+            normal[i] /= normal[i].Norml2();
+         }
+         else
+         {
+            normal[i] = normal[0];
+            normal[i] *= -1;
+         }
+         els[i]->CalcPhysDShape(*el_trans[i], dshape[i]);
+         els[i]->CalcPhysHessian(*el_trans[i], hessian[i]);
+         dshape[i].Mult(normal[i], dnshape[i]);
+         nv[i][0] = normal[i][0]*normal[i][0];
+         nv[i][1] = 2*normal[i][0]*normal[i][1];
+         nv[i][2] = normal[i][1]*normal[i][1];
+         hessian[i].Mult(nv[i], nd2nshape[i]);
+      }
+
+      // Compute blocks
+      for (int i = 0; i < num_elems; i++)
+      {
+         for (int j = 0; j < num_elems; j++)
+         {
+            blockJ[i][j] = 0.0;
+            blockC[i][j] = 0.0;
+            AddMult_a_VWt(-1.0, dnshape[i], nd2nshape[j], blockJ[i][j]);
+            elmatJ_p.SetSubMatrix(i*ndof[0], j*ndof[0], blockJ[i][j]);
+
+            AddMult_a_VWt(eta/h_e, dnshape[i], dnshape[j], blockC[i][j]);
+            elmatC_p.SetSubMatrix(i*ndof[0], j*ndof[0], blockC[i][j]);
+         }
+      }
+
+      // Symmetrize the jump term
+      elmatJ_p.Symmetrize();
+      if (!ndof[1])
+      {
+         elmatJ_p *= 2;
+      }
+
+      // Add penalty term
+      elmatJ_p += elmatC_p;
+      elmatJ_p *= ip.weight * Trans.Weight();
+      elmat += elmatJ_p;
+   }
+}

fem/datacollection.cpp

@@ -809,7 +809,7 @@ ParaViewDataCollectionBase::ParaViewDataCollectionBase(
 
 void ParaViewDataCollectionBase::SetLevelsOfDetail(int levels_of_detail_)
 {
-   levels_of_detail = levels_of_detail_;
+   levels_of_detail = std::max(levels_of_detail_, 1);
 }
 
 void ParaViewDataCollectionBase::SetHighOrderOutput(bool high_order_output_)

general/forall.hpp

@@ -480,8 +480,8 @@ template <typename DBODY>
 void RajaHipWrap1D(const int N, DBODY &&d_body)
 {
    //true denotes asynchronous kernel
-   RAJA::forall<RAJA::hip_exec<MFEM_HIP_BLOCKS,true> >(RAJA::RangeSegment(0,N),
-                                                       d_body);
+   RAJA::forall<RAJA::hip_exec<MFEM_HIP_BLOCKS, true> >(
+      Device::GetRajaResource(), RAJA::RangeSegment(0, N), d_body);
 }
 
 template <typename DBODY>