Merge latest changes from upstream/master.

.github/workflows/stale.yml

@@ -0,0 +1,31 @@
+# This workflow warns and then closes issues and PRs that have had no activity for a specified amount of time.
+# For more information, see: https://github.com/actions/stale
+name: Mark stale issues and pull requests
+
+on:
+  workflow_dispatch:
+  schedule:
+  - cron: '0 0 * * *'
+
+jobs:
+  stale:
+
+    runs-on: ubuntu-latest
+    permissions:
+      issues: write
+      pull-requests: write
+      actions: write
+
+    steps:
+    - uses: actions/stale@v9
+      with:
+        repo-token: ${{ secrets.GITHUB_TOKEN }}
+        stale-issue-message: ':warning: This issue has been automatically marked as stale because it has not had any activity in the last month. *If no activity occurs in the next week, it will be automatically closed.* Thank you for your contributions.'
+        stale-pr-message: ':warning: This PR has been automatically marked as stale because it has not had any activity in the last month. *If no activity occurs in the next week, it will be automatically closed.* Thank you for your contributions.'
+        days-before-stale: 30
+        days-before-close: 7
+        stale-issue-label: 'stale'
+        stale-pr-label: 'stale'
+        operations-per-run: 500
+        exempt-issue-labels: "bug,WIP,ready-for-review,in-review,in-next"
+        exempt-pr-labels: "bug,WIP,ready-for-review,in-review,in-next"

.gitignore

@@ -8,13 +8,17 @@
 
 # Object and library files
 *.o
+*.o.tmp
 /libmfem.*
 /miniapps/common/libmfem-common.*
 
 # CMake generated files
 CMakeCache.txt
 CMakeFiles/
 
+# Clangd server cache
+*.cache*
+
 # Backup files
 *~
 

CHANGELOG

@@ -10,6 +10,8 @@
 
 Version 4.7.1 (development)
 ===========================
+- Refactored ALGOIM cut integration rules. The interface is unified with
+  the interface for moment based cut integration rules.
 
 Discretization improvements
 ---------------------------
@@ -18,6 +20,8 @@ Discretization improvements
 
 - Added support for boundary constraints to the hybridization class.
 
+- Added support for external boundary submeshes with nonconformal mesh adaptation.
+
 Meshing improvements
 --------------------
 - The ExodusII reader now handles pyramid and wedge element types. Mixed meshes
@@ -38,8 +42,27 @@ GPU computing
 - Added support for GPU-accelerated batched linear algebra (using cuBLAS,
   hipBLAS, MAGMA, or native MFEM functionality) through the BatchedLinAlg class.
 
+- A new GPU kernel dispatch mechanism was introduced. Users can instantiate
+  specialized kernels for specific combinations of (for example) polynomial
+  degree and number of quadrature points using
+  `DiffusionIntegrator::AddSpecialization` and
+  `MassIntegrator::AddSpecialization` (this functionality may be added to more
+  integrators in the future).
+
+- Calls to slower fallback kernels can be reported to `mfem::err` by setting
+  the environment variable `MFEM_REPORT_KERNELS` to any value other than `NO`
+  or by explicitly calling `KernelReporter::Enable`. Users can then add
+  specializations for these kernels to achieve higher performance.
+
+- Element assembly kernels have been added for low-order refined to
+  high-order transfer operators. New kernels can be offloaded as device
+  kernels. Example usage may be found in lor-transfer.cpp under miniapps/tools.
+
 Miscellaneous
 -------------
+- Added support for SUNDIALS v7. See the section "API changes" for some small
+  changes related to this new version.
+
 - Refactored the `ARKStepSolver` class (ARKODE interface) to use
   `TimeDependentOperator::Mult` only when the associated ODE operator is
   expressed in explicit form (i.e., `TimeDependentOperator::isExplicit()`),
@@ -52,10 +75,27 @@ Miscellaneous
 
 - Added support for custom interpolation procedure in FindPointsGSLIB.
 
+- `FiniteElementSpace` has new methods to directly set prolongation and
+  restriction operators to arbitrary sparse matrices.
+
+- There are new convenience constructors for NURBS patches and knot vectors.
+
 API changes
 -----------
 - API change: in class GridFunction, 'fec' was renamed to 'fec_owned'.
 
+- API change: support for SUNDIALS v7:
+  * the SUNDIALS types `realtype` and `booleantype` are no longer defined by v7
+    and therefore MFEM now uses the new type names `sunrealtype` and
+    `sunbooleantype`, respectively, which MFEM defines when using SUNDIALS < v6
+    where these types were not defined.
+  * The SUNDIALS macro `SUNLS_SUCCESS` and some other `*_SUCCESS` macros were
+    removed and replaced by `SUN_SUCCESS` in v7, so to avoid tedious checks for
+    SUNDIALS versions, MFEM now defines and uses the constant `SUN_SUCCESS` when
+    using SUNDIALS < v7.
+  * The constants `SUN_PREC_*`, introduced by SUNDIALS v6 are now introduced by
+    MFEM when using SUNDIALS < v6 to avoid tedious version checks.
+
 
 Version 4.7, released on May 7, 2024
 ====================================
@@ -142,6 +182,15 @@ New and updated examples and miniapps
 - Added two new example codes: 38 and 39/39p described above. Substantially
   updated Example 18/18p.
 
+- Added ODE solvers selection routines. This creates a uniformity across examples,
+  miniapps and other executables in regard to ODE(time-integrator) selection.
+
+- Added new mechanism for retrieving and setting state vectors in ODE solvers.
+  This is relevant for AB/AM and gen-alpha solvers.
+
+- Added ODEsolver/ODEsolver2 unit tests to verify order of convergence and
+  read/write functionality.
+
 Miscellaneous
 -------------
 - Updated the Doxygen documentation style, which now requires Doxygen version

CMakeLists.txt

@@ -340,7 +340,10 @@ if (MFEM_USE_SUNDIALS)
   if (MFEM_USE_HIP)
     list(APPEND SUNDIALS_COMPONENTS NVector_Hip)
   endif()
-  find_package(SUNDIALS REQUIRED ${SUNDIALS_COMPONENTS})
+  # The Core component was added in SUNDIALS v7, so we treat it as optional in
+  # order to support older versions.
+  find_package(SUNDIALS REQUIRED ${SUNDIALS_COMPONENTS}
+    OPTIONAL_COMPONENTS Core)
 endif()
 
 # SuperLU_DIST can only be enabled in parallel

INSTALL

@@ -502,10 +502,14 @@ MFEM_USE_CODIPACK = YES/NO
 MFEM_USE_ALGOIM = YES/NO
   Enable the usage of Algoim - a collection of high-order accurate numerical
   methods and C++ algorithms for working with implicitly-defined geometry and
-  level set methods. The Algoim library requires the Blitz++ library. The MFEM
-  provides interface to Algoim v1. Thus, to check out the specific state use:
+  level set methods, see https://algoim.github.io. MFEM provides interface to
+  Algoim v1. To check out the specific Algoim state use:
+  https://github.com/algoim/algoim
   git checkout 9c9ca0ef094d8ab0390ed36367a1151b459bbe0a
-  https://algoim.github.io
+  The Algoim library requires the Blitz++ library. To use the latest state of
+  Blitz++ that has been tested with MFEM, use:
+  https://github.com/blitzpp/blitz
+  git checkout f24a250a43dff88c31ad92916da828b7ea9a98b7
 
 MFEM_USE_ADFORWARD = YES/NO
   Enable forward mode for AD packages. This option is valid

config/cmake/modules/FindSUNDIALS.cmake

@@ -31,4 +31,5 @@ mfem_find_package(SUNDIALS SUNDIALS SUNDIALS_DIR
   ADD_COMPONENT CVODE "include" cvode/cvode.h "lib" sundials_cvode
   ADD_COMPONENT CVODES "include" cvodes/cvodes.h "lib" sundials_cvodes
   ADD_COMPONENT ARKODE "include" arkode/arkode.h "lib" sundials_arkode
-  ADD_COMPONENT KINSOL "include" kinsol/kinsol.h "lib" sundials_kinsol)
+  ADD_COMPONENT KINSOL "include" kinsol/kinsol.h "lib" sundials_kinsol
+  ADD_COMPONENT Core "include" sundials/sundials_core.h "lib" sundials_core)

config/defaults.mk

@@ -289,6 +289,13 @@ endif
 ifeq ($(MFEM_USE_HIP),YES)
    SUNDIALS_LIB += -lsundials_nvechip
 endif
+SUNDIALS_CORE_PAT = $(subst\
+ @MFEM_DIR@,$(MFEM_DIR),$(SUNDIALS_DIR))/lib*/libsundials_core.*
+ifeq ($(MFEM_USE_SUNDIALS),YES)
+   ifneq ($(wildcard $(SUNDIALS_CORE_PAT)),)
+      SUNDIALS_LIB += -lsundials_core
+   endif
+endif
 # If SUNDIALS was built with KLU:
 # MFEM_USE_SUITESPARSE = YES
 
@@ -533,8 +540,10 @@ ifdef GOTCHA_DIR
 endif
 
 # BLITZ library configuration
-BLITZ_DIR = @MFEM_DIR@/../blitz
+# BLITZ_DIR must be the custom installation folder (-DCMAKE_INSTALL_PREFIX).
+BLITZ_DIR = @MFEM_DIR@/../blitz/install
 BLITZ_OPT = -I$(BLITZ_DIR)/include
+# On intel machines, use /lib64 instead of /lib.
 BLITZ_LIB = $(XLINKER)-rpath,$(BLITZ_DIR)/lib -L$(BLITZ_DIR)/lib -lblitz
 
 # ALGOIM library configuration

examples/ex10.cpp

@@ -3,14 +3,14 @@
 // Compile with: make ex10
 //
 // Sample runs:
-//    ex10 -m ../data/beam-quad.mesh -s 3 -r 2 -o 2 -dt 3
-//    ex10 -m ../data/beam-tri.mesh -s 3 -r 2 -o 2 -dt 3
-//    ex10 -m ../data/beam-hex.mesh -s 2 -r 1 -o 2 -dt 3
-//    ex10 -m ../data/beam-tet.mesh -s 2 -r 1 -o 2 -dt 3
-//    ex10 -m ../data/beam-wedge.mesh -s 2 -r 1 -o 2 -dt 3
-//    ex10 -m ../data/beam-quad.mesh -s 14 -r 2 -o 2 -dt 0.03 -vs 20
-//    ex10 -m ../data/beam-hex.mesh -s 14 -r 1 -o 2 -dt 0.05 -vs 20
-//    ex10 -m ../data/beam-quad-amr.mesh -s 3 -r 2 -o 2 -dt 3
+//    ex10 -m ../data/beam-quad.mesh -s 23 -r 2 -o 2 -dt 3
+//    ex10 -m ../data/beam-tri.mesh -s 23 -r 2 -o 2 -dt 3
+//    ex10 -m ../data/beam-hex.mesh -s 22 -r 1 -o 2 -dt 3
+//    ex10 -m ../data/beam-tet.mesh -s 22 -r 1 -o 2 -dt 3
+//    ex10 -m ../data/beam-wedge.mesh -s 22 -r 1 -o 2 -dt 3
+//    ex10 -m ../data/beam-quad.mesh -s 4 -r 2 -o 2 -dt 0.03 -vs 20
+//    ex10 -m ../data/beam-hex.mesh -s 4 -r 1 -o 2 -dt 0.05 -vs 20
+//    ex10 -m ../data/beam-quad-amr.mesh -s 23 -r 2 -o 2 -dt 3
 //
 // Description:  This examples solves a time dependent nonlinear elasticity
 //               problem of the form dv/dt = H(x) + S v, dx/dt = v, where H is a
@@ -87,16 +87,16 @@ class HyperelasticOperator : public TimeDependentOperator
                         real_t visc, real_t mu, real_t K);
 
    /// Compute the right-hand side of the ODE system.
-   virtual void Mult(const Vector &vx, Vector &dvx_dt) const;
+   void Mult(const Vector &vx, Vector &dvx_dt) const override;
    /** Solve the Backward-Euler equation: k = f(x + dt*k, t), for the unknown k.
        This is the only requirement for high-order SDIRK implicit integration.*/
-   virtual void ImplicitSolve(const real_t dt, const Vector &x, Vector &k);
+   void ImplicitSolve(const real_t dt, const Vector &x, Vector &k) override;
 
    real_t ElasticEnergy(const Vector &x) const;
    real_t KineticEnergy(const Vector &v) const;
    void GetElasticEnergyDensity(const GridFunction &x, GridFunction &w) const;
 
-   virtual ~HyperelasticOperator();
+   ~HyperelasticOperator() override;
 };
 
 /** Nonlinear operator of the form:
@@ -120,12 +120,12 @@ class ReducedSystemOperator : public Operator
    void SetParameters(real_t dt_, const Vector *v_, const Vector *x_);
 
    /// Compute y = H(x + dt (v + dt k)) + M k + S (v + dt k).
-   virtual void Mult(const Vector &k, Vector &y) const;
+   void Mult(const Vector &k, Vector &y) const override;
 
    /// Compute J = M + dt S + dt^2 grad_H(x + dt (v + dt k)).
-   virtual Operator &GetGradient(const Vector &k) const;
+   Operator &GetGradient(const Vector &k) const override;
 
-   virtual ~ReducedSystemOperator();
+   ~ReducedSystemOperator() override;
 };
 
 
@@ -141,8 +141,8 @@ class ElasticEnergyCoefficient : public Coefficient
 public:
    ElasticEnergyCoefficient(HyperelasticModel &m, const GridFunction &x_)
       : model(m), x(x_) { }
-   virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
-   virtual ~ElasticEnergyCoefficient() { }
+   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
+   ~ElasticEnergyCoefficient() override { }
 };
 
 void InitialDeformation(const Vector &x, Vector &y);
@@ -160,7 +160,7 @@ int main(int argc, char *argv[])
    const char *mesh_file = "../data/beam-quad.mesh";
    int ref_levels = 2;
    int order = 2;
-   int ode_solver_type = 3;
+   int ode_solver_type = 23;
    real_t t_final = 300.0;
    real_t dt = 3.0;
    real_t visc = 1e-2;
@@ -177,11 +177,7 @@ int main(int argc, char *argv[])
    args.AddOption(&order, "-o", "--order",
                   "Order (degree) of the finite elements.");
    args.AddOption(&ode_solver_type, "-s", "--ode-solver",
-                  "ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
-                  "            11 - Forward Euler, 12 - RK2,\n\t"
-                  "            13 - RK3 SSP, 14 - RK4."
-                  "            22 - Implicit Midpoint Method,\n\t"
-                  "            23 - SDIRK23 (A-stable), 24 - SDIRK34");
+                  ODESolver::Types.c_str());
    args.AddOption(&t_final, "-tf", "--t-final",
                   "Final time; start time is 0.");
    args.AddOption(&dt, "-dt", "--time-step",
@@ -213,28 +209,7 @@ int main(int argc, char *argv[])
    // 3. Define the ODE solver used for time integration. Several implicit
    //    singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
    //    explicit Runge-Kutta methods are available.
-   ODESolver *ode_solver;
-   switch (ode_solver_type)
-   {
-      // Implicit L-stable methods
-      case 1: ode_solver = new BackwardEulerSolver; break;
-      case 2: ode_solver = new SDIRK23Solver(2); break;
-      case 3: ode_solver = new SDIRK33Solver; break;
-      // Explicit methods
-      case 11: ode_solver = new ForwardEulerSolver; break;
-      case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
-      case 13: ode_solver = new RK3SSPSolver; break;
-      case 14: ode_solver = new RK4Solver; break;
-      case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
-      // Implicit A-stable methods (not L-stable)
-      case 22: ode_solver = new ImplicitMidpointSolver; break;
-      case 23: ode_solver = new SDIRK23Solver; break;
-      case 24: ode_solver = new SDIRK34Solver; break;
-      default:
-         cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
-         delete mesh;
-         return 3;
-   }
+   unique_ptr<ODESolver> ode_solver = ODESolver::Select(ode_solver_type);
 
    // 4. Refine the mesh to increase the resolution. In this example we do
    //    'ref_levels' of uniform refinement, where 'ref_levels' is a
@@ -371,7 +346,6 @@ int main(int argc, char *argv[])
    }
 
    // 10. Free the used memory.
-   delete ode_solver;
    delete mesh;
 
    return 0;