Silenced tee.

Range and curl dimension in physical space

.github/actions/sanitize/mfem/action.yml

@@ -25,7 +25,7 @@ runs:
   steps:
     - uses: ./.github/actions/sanitize/config
 
-    - uses: actions/cache@v5
+    - uses: actions/cache@v4
       if: ${{env.DEBUG == 'true'}}
       id: debug
       with:

.github/actions/sanitize/restore/action.yml

@@ -36,7 +36,7 @@ runs:
   steps:
     - uses: ./.github/actions/sanitize/config
 
-    - uses: actions/cache@v5
+    - uses: actions/cache@v4
       if: ${{env.DEBUG == 'true' && inputs.cache-skip != 'true'}}
       id: debug
       with:

.github/actions/sanitize/setup/action.yml

@@ -23,7 +23,7 @@ inputs:
 runs:
   using: 'composite'
   steps:
-    - uses: actions/cache/restore@v5  # Cache for LLVM libcxx
+    - uses: actions/cache/restore@v4  # Cache for LLVM libcxx
       with:
         path: ${{env.LLVM_DIR}}
         fail-on-cache-miss: true
@@ -32,14 +32,14 @@ runs:
     - uses: ./.github/actions/sanitize/mpi
       if: ${{inputs.par == 'true'}}
 
-    - uses: actions/cache/restore@v5  # Cache for Hypre
+    - uses: actions/cache/restore@v4  # Cache for Hypre
       if: ${{inputs.par == 'true'}}
       with:
         path: ${{env.HYPRE_DIR}}
         fail-on-cache-miss: true
         key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5
 
-    - uses: actions/cache/restore@v5  # Cache for Metis
+    - uses: actions/cache/restore@v4  # Cache for Metis
       if: ${{inputs.par == 'true'}}
       with:
         path: ${{env.METIS_DIR}}
@@ -51,13 +51,13 @@ runs:
       run: ln -s -f ${{env.HYPRE_DIR}} hypre && ln -s -f ${{env.METIS_DIR}} metis-4.0
       shell: bash
 
-    - uses: actions/cache/restore@v5  # Cache for LSAN suppression file
+    - uses: actions/cache/restore@v4  # Cache for LSAN suppression file
       with:
         path: ${{env.LSAN_DIR}}
         fail-on-cache-miss: true
         key: build-lsan-suppression-file
 
-    - uses: actions/checkout@v6  # Checkout the repository
+    - uses: actions/checkout@v4  # Checkout the repository
       with:
         path: mfem
         # ref: ${{env.BRANCH}}

.github/workflows/build-container.yml

@@ -43,7 +43,7 @@ jobs:
           remove-docker-images: 'true'
 
       - name: Checkout
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
 
       # It's easier to reference named variables than indexes of the matrix
       - name: Set Environment

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

@@ -153,7 +153,7 @@ jobs:
       # /home/runner/work/mfem/mfem/mfem
       # Note: Done now to access "install-hypre" and "install-metis" actions.
       - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
         with:
           path: ${{ env.MFEM_TOP_DIR }}
           # Fetch the complete history for codecov to access commits ID
@@ -225,7 +225,7 @@ jobs:
       - name: cache hypre
         id: hypre-cache
         if: matrix.mpi == 'par'
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: ${{ env.HYPRE_TOP_DIR }}
           key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-${{ matrix.precision }}-v2.5
@@ -255,7 +255,7 @@ jobs:
       - name: cache metis
         id: metis-cache
         if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: ${{ env.METIS_TOP_DIR }}
           key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5
@@ -270,7 +270,7 @@ jobs:
       - name: cache vcpkg (Windows)
         id: vcpkg-cache
         if: matrix.os == 'windows-latest'
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: vcpkg_cache
           key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1

.github/workflows/codeql-analysis.yml

@@ -40,11 +40,11 @@ jobs:
 
     steps:
       - name: Checkout repository
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
 
       # Initializes the CodeQL tools for scanning.
       - name: Initialize CodeQL
-        uses: github/codeql-action/init@v4
+        uses: github/codeql-action/init@v2
         with:
           languages: ${{ matrix.language }}
           # If you wish to specify custom queries, you can do so here or in a config file.
@@ -57,7 +57,7 @@ jobs:
       # Autobuild attempts to build any compiled languages  (C/C++, C#, or Java).
       # If this step fails, then you should remove it and run the build manually (see below)
       - name: Autobuild
-        uses: github/codeql-action/autobuild@v4
+        uses: github/codeql-action/autobuild@v2
 
       # ℹ️ Command-line programs to run using the OS shell.
       # 📚 See https://docs.github.com/en/actions/using-workflows/workflow-syntax-for-github-actions#jobsjob_idstepsrun
@@ -70,4 +70,4 @@ jobs:
       #   ./location_of_script_within_repo/buildscript.sh
 
       - name: Perform CodeQL Analysis
-        uses: github/codeql-action/analyze@v4
+        uses: github/codeql-action/analyze@v2

.github/workflows/mfem-analysis.yml

@@ -39,7 +39,7 @@ jobs:
 
     steps:
       - name: checkout MFEM
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
         with:
           path: mfem
 
@@ -50,7 +50,7 @@ jobs:
 
       - name: Cache Hypre Install
         id: hypre-cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: ${{ env.HYPRE_TOP_DIR }}
           key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-v2.5
@@ -65,7 +65,7 @@ jobs:
 
       - name: Cache Metis Install
         id: metis-cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: ${{ env.METIS_TOP_DIR }}
           key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5

.github/workflows/repo-check.yml

@@ -38,7 +38,7 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
       - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
 
       - name: copyright check
         id: copyright
@@ -93,7 +93,7 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
       - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
 
       - name: get astyle
         run: |
@@ -110,7 +110,7 @@ jobs:
        github.event.pull_request.head.repo.full_name != github.repository)
     steps:
       - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
 
       - name: get doxygen and graphviz
         run: |
@@ -135,7 +135,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
         with:
           fetch-depth: 0
 

.github/workflows/sanitize-build-hypre.yml

@@ -17,11 +17,11 @@ jobs:
     runs-on: ubuntu-latest
     name: 2.19.0
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/config
       - name: Cache
         id: cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: ${{env.HYPRE_DIR}}
           key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5

.github/workflows/sanitize-build-libcxx.yml

@@ -27,13 +27,13 @@ jobs:
             llvm_use_sanitizer: "Undefined"
     name: ${{matrix.sanitizer}}
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/config
         with:
           NO_FLAGS: true
       - name: Cache
         id: cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: ${{env.LLVM_DIR}}
           key: build-libcxx-${{env.LLVM_VER}}-${{matrix.sanitizer}}

.github/workflows/sanitize-build-lsan.yml

@@ -17,11 +17,11 @@ jobs:
     runs-on: ubuntu-latest
     name: lsan.supp
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/config
       - name: Cache
         id: cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: ${{env.LSAN_DIR}}
           key: build-lsan-suppression-file

.github/workflows/sanitize-build-metis.yml

@@ -17,11 +17,11 @@ jobs:
     runs-on: ubuntu-latest
     name: 4.0.3
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/config
       - name: Cache
         id: cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
         with:
           path: ${{env.METIS_DIR}}
           key: ${{runner.os}}-build-${{env.METIS_DIR}}-v2.5

.github/workflows/sanitize-tests.yml

@@ -28,7 +28,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/mfem
         with:
           par: ${{inputs.par}}
@@ -40,7 +40,7 @@ jobs:
     env:
       ex: ${{inputs.par && 'ex1p' || 'ex1'}}
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -58,7 +58,7 @@ jobs:
     env:
       exclude: ${{inputs.par && '-E "_ser"' || ''}}
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -82,7 +82,7 @@ jobs:
     env:
       exclude: ${{inputs.par && '-E "_ser"' || ''}}
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -107,7 +107,7 @@ jobs:
       run: ${{inputs.par && '-R "_cpu_np"' || ''}}
       exclude: ${{inputs.par && '"unit_tests|debug"' || '"^unit_tests$|debug"'}}
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -131,7 +131,7 @@ jobs:
     env:
       unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:
@@ -165,7 +165,7 @@ jobs:
       unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
       np: ${{inputs.par && '_np=2' || ''}}
     steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
       - uses: ./.github/actions/sanitize/restore
         id: restore
         with:

CHANGELOG

@@ -8,22 +8,6 @@
                                https://mfem.org
 
 
-Version 4.10 (development)
-==========================
-
-Discretization improvements
----------------------------
-- Replaced legacy simplex quadrature rules with symmetric positive-weight
-  rules for triangles (orders 0-25) and tetrahedra (orders 0-20). These
-  rules guarantee all-positive weights and interior quadrature points,
-  improving numerical stability. Higher orders fall back to Grundmann-Moller.
-  Triangle rules: Witherden & Vincent, Comput. Math. Appl. 69(10):1232-1241,
-  2015.
-  Tet rules (d=1-13): Witherden & Vincent (ibid).
-  Tet rules (d=14-20): Chuluunbaatar et al., Comput. Math. Appl. 124:89-97,
-  2022.
-
-
 Version 4.9.1 (development)
 ===========================
 

examples/ex12p.cpp

@@ -5,9 +5,9 @@
 // Sample runs:
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh
-//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 464 -n 10 -o 2 -elast
+//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 462 -n 10 -o 2 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-hex.mesh -s 3878
-//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 82
+//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 81
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh -s 3877 -o 2 -sys
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh -s 4544 -n 6 -o 3 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-quad-nurbs.mesh

examples/ex1p.cpp

@@ -66,6 +66,26 @@
 using namespace std;
 using namespace mfem;
 
+class VisServer
+{
+   const char *path{};
+   int port{-1};
+
+   struct Deleter { void operator()(FILE *); };
+   std::unique_ptr<FILE, Deleter> fglvis;
+
+   void Open();
+public:
+   VisServer(const char *path_, int port_ = 19916)
+      : path(path_), port(port_) { Open(); }
+
+   explicit VisServer() = default;
+
+   inline const char *GetPath() const { return path; }
+   inline int GetPort() const { return port; }
+   inline bool IsOpen() const { return (fglvis != nullptr); }
+};
+
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI and HYPRE.
@@ -82,6 +102,7 @@ int main(int argc, char *argv[])
    bool fa = false;
    const char *device_config = "cpu";
    bool visualization = true;
+   const char *visserver = "";
    bool algebraic_ceed = false;
 
    OptionsParser args(argc, argv);
@@ -106,6 +127,8 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&visserver, "-vs", "--visualization-server",
+                  "Path to GLVis binary to start own server.");
    args.Parse();
    if (!args.Good())
    {
@@ -301,6 +324,14 @@ int main(int argc, char *argv[])
    {
       char vishost[] = "localhost";
       int  visport   = 19916;
+
+      VisServer server;
+      if (strlen(visserver) > 0)
+      {
+         server = VisServer(visserver, visport);
+         visport = server.GetPort();
+      }
+
       socketstream sol_sock(vishost, visport);
       sol_sock << "parallel " << num_procs << " " << myid << "\n";
       sol_sock.precision(8);
@@ -315,3 +346,66 @@ int main(int argc, char *argv[])
 
    return 0;
 }
+
+void VisServer::Open()
+{
+   if (!path || strlen(path) <= 0) { return; }
+
+#ifdef MFEM_USE_MPI
+   if (Mpi::Root())
+#endif // MFEM_USE_MPI
+   {
+      // Open a log file
+      FILE *ftmp = popen("mktemp --tmpdir glvis-log.XXXXXX", "r");
+      constexpr size_t ssize = 256;
+      char stmp[ssize];
+      fgets(stmp, ssize, ftmp);
+      pclose(ftmp);
+      int len = strlen(stmp);
+      if (len > 0 && stmp[len-1] == '\n') { stmp[--len] = '\0'; }
+      std::cout << "Starting GLVis log in: " << stmp << std::endl;
+
+      // Start the server
+      std::stringstream ss;
+      ss << "trap '' SIGPIPE && " << path <<
+         " -no-pr -no-ex 2>&1 | tee \"" << stmp <<
+         "\" 2> /dev/null | grep -m 1 ^GLVIS_SERVER_PORT";
+      fglvis.reset(popen(ss.str().c_str(), "r"));
+
+      // Capture the port number
+      char line[ssize];
+      bool captured = false;
+      while (fgets(line, ssize, fglvis.get()))
+      {
+         if (strncmp(line, "GLVIS_SERVER_PORT=", 18) == 0)
+         {
+            sscanf(line, "GLVIS_SERVER_PORT=%d", &port);
+            captured = true;
+            break;
+         }
+      }
+      if (captured)
+      {
+         std::cout << "Started GLVis server at port: " << port << std::endl;
+      }
+      else
+      {
+         std::cerr << "GLVis server did not start normally." << std::endl;
+      }
+   }
+#ifdef MFEM_USE_MPI
+   MPI_Bcast(&port, 1, MPI_INT, 0, MPI_COMM_WORLD);
+#endif // MFEM_USE_MPI
+}
+
+void VisServer::Deleter::operator()(FILE *f)
+{
+   if (f != nullptr)
+   {
+      int ierr = pclose(f);
+      if (ierr != 0)
+      {
+         std::cerr << "GLVis server pclose() returns: " << ierr << std::endl;
+      }
+   }
+}

examples/ex22p.cpp

@@ -302,21 +302,15 @@ int main(int argc, char *argv[])
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
       sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
       sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
       sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                  << "window_title 'Exact: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
       sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                  << "window_title 'Exact: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
    }
 
    // 11. Set up the parallel sesquilinear form a(.,.) on the finite element
@@ -540,21 +534,15 @@ int main(int argc, char *argv[])
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
       sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
       sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
       sol_sock_r << "solution\n" << *pmesh << u.real()
                  << "window_title 'Solution: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
       sol_sock_i << "solution\n" << *pmesh << u.imag()
                  << "window_title 'Solution: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
    }
    if (visualization && exact_sol)
    {
@@ -563,21 +551,15 @@ int main(int argc, char *argv[])
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
       sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
       sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
       sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                  << "window_title 'Error: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
       sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                  << "window_title 'Error: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
    }
    if (visualization)
    {

examples/ex36.cpp

@@ -97,7 +97,13 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
-   args.ParseCheck();
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
 
    // 2. Read the mesh from the mesh file.
    const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -122,7 +128,7 @@ int main(int argc, char *argv[])
    *nodes /= scale;
 
    // 4. Define the necessary finite element spaces on the mesh.
-   H1Bubble_FECollection H1fec(order, order - 1, dim);
+   H1_FECollection H1fec(order+1, dim);
    FiniteElementSpace H1fes(&mesh, &H1fec);
 
    L2_FECollection L2fec(order-1, dim);

examples/ex36p.cpp

@@ -103,7 +103,19 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
-   args.ParseCheck();
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
 
    // 2. Read the mesh from the mesh file.
    const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -131,7 +143,7 @@ int main(int argc, char *argv[])
    mesh.Clear();
 
    // 4. Define the necessary finite element spaces on the mesh.
-   H1Bubble_FECollection H1fec(order, order - 1, dim);
+   H1_FECollection H1fec(order+1, dim);
    ParFiniteElementSpace H1fes(&pmesh, &H1fec);
 
    L2_FECollection L2fec(order-1, dim);

fem/CMakeLists.txt

@@ -73,7 +73,6 @@ set(SRCS
   fe/fe_base.cpp
   fe/fe_fixed_order.cpp
   fe/fe_h1.cpp
-  fe/fe_h1_bubble.cpp
   fe/fe_l2.cpp
   fe/fe_nd.cpp
   fe/fe_nurbs.cpp
@@ -134,7 +133,7 @@ set(SRCS
   tmop/assemble/diag2.cpp
   tmop/assemble/grad2_limit.cpp
   tmop/assemble/grad2.cpp
-  tmop/assemble/diag3_limit.cpp
+  tmop/assemble/diag3_limit.cpp  
   tmop/assemble/diag3.cpp
   tmop/assemble/grad3_limit.cpp
   tmop/assemble/grad3.cpp
@@ -222,7 +221,6 @@ set(HDRS
   fe/fe_base.hpp
   fe/fe_fixed_order.hpp
   fe/fe_h1.hpp
-  fe/fe_h1_bubble.hpp
   fe/fe_l2.hpp
   fe/fe_nd.hpp
   fe/fe_nurbs.hpp

fem/bilinearform.cpp

@@ -729,8 +729,7 @@ void BilinearForm::Assemble(int skip_zeros)
          tr = mesh -> GetBdrFaceTransformations (i);
          if (tr != NULL)
          {
-            mfem::DofTransformation doftrans;
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
             fe1 = fes -> GetFE (tr -> Elem1No);
             // The fe2 object is really a dummy and not used on the boundaries,
             // but we can't dereference a NULL pointer, and we don't want to
@@ -744,7 +743,6 @@ void BilinearForm::Assemble(int skip_zeros)
 
                boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
                                                               elemmat);
-               doftrans.TransformDual(elemmat);
                mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
             }
          }
@@ -1725,7 +1723,6 @@ void MixedBilinearForm::Assemble(int skip_zeros)
          }
       }
 
-      DofTransformation dom_dof_trans, ran_dof_trans;
       for (int i = 0; i < trial_fes -> GetNBE(); i++)
       {
          const int bdr_attr = mesh->GetBdrAttribute(i);
@@ -1734,8 +1731,8 @@ void MixedBilinearForm::Assemble(int skip_zeros)
          ftr = mesh -> GetBdrFaceTransformations (i);
          if (ftr != NULL)
          {
-            trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs, dom_dof_trans);
-            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs, ran_dof_trans);
+            trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
+            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
             trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
             test_fe1 = test_fes->GetFE(ftr->Elem1No);
             // The test_fe2 object is really a dummy and not used on the
@@ -1751,7 +1748,6 @@ void MixedBilinearForm::Assemble(int skip_zeros)
                boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
                                                            *test_fe2,
                                                            *ftr, elemmat);
-               TransformDual(ran_dof_trans, dom_dof_trans, elemmat);
                mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
             }
          }

fem/bilininteg.hpp

@@ -2710,7 +2710,7 @@ public:
 
 
 /** Integrator for $(-Q u, \nabla v)$ for Nedelec ($u$) and $H^1$ ($v$) elements.
-    This is equivalent to a weak divergence of the $H(curl)$ basis functions. */
+    This is equivalent to a weak divergence of the $H(curl$ basis functions. */
 class VectorFEWeakDivergenceIntegrator: public BilinearFormIntegrator
 {
 protected:

fem/fe.hpp

@@ -20,7 +20,6 @@
 #include "fe/fe_base.hpp"
 #include "fe/fe_fixed_order.hpp"
 #include "fe/fe_h1.hpp"
-#include "fe/fe_h1_bubble.hpp"
 #include "fe/fe_nd.hpp"
 #include "fe/fe_rt.hpp"
 #include "fe/fe_l2.hpp"

fem/fe/fe_h1_bubble.cpp

@@ -1,973 +0,0 @@
-// 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.
-
-// H1 Finite Element classes
-
-#include "fe_h1_bubble.hpp"
-
-namespace mfem
-{
-
-using namespace std;
-
-H1Bubble_TriangleElement::H1Bubble_TriangleElement(int p, int q, int btype)
-   : NodalFiniteElement(2, Geometry::TRIANGLE, 3*p + ((q+1)*(q+2))/2,
-                        max(p, 3 + q), FunctionSpace::Pk),
-     base_order(p), bubble_order(q)
-{
-   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
-   const real_t *cp2 = poly1d.ClosedPoints(
-                          q + 3, VerifyNodal(VerifyClosed(btype)));
-
-   const int n1d = max(p + 1, q + 1);
-   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
-
-#ifndef MFEM_THREAD_SAFE
-   shape_x.SetSize(n1d);
-   shape_y.SetSize(n1d);
-   shape_l.SetSize(n1d);
-   dshape_x.SetSize(n1d);
-   dshape_y.SetSize(n1d);
-   dshape_l.SetSize(n1d);
-   u.SetSize(npq);
-   du.SetSize(npq, dim);
-#endif
-
-   // vertices
-   Nodes.IntPoint(0).Set2(cp[0], cp[0]);
-   Nodes.IntPoint(1).Set2(cp[p], cp[0]);
-   Nodes.IntPoint(2).Set2(cp[0], cp[p]);
-
-   // edges
-   int o = 3;
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set2(cp[p-i], cp[i]);
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
-   }
-
-   // Interior P_{q+3} nodes
-   for (int j = 1; j < q + 3; j++)
-   {
-      for (int i = 1; i + j < q + 3; i++)
-      {
-         const real_t w = cp2[i] + cp2[j] + cp2[q+3-i-j];
-         Nodes.IntPoint(o++).Set2(cp2[i]/w, cp2[j]/w);
-      }
-   }
-
-#ifdef MFEM_THREAD_SAFE
-   Vector shape_x(n1d), shape_y(n1d), shape_l(n1d);
-#endif
-
-   DenseMatrix Tt(dof, npq);
-   for (int k = 0; k < dof; ++k)
-   {
-      const IntegrationPoint &ip = Nodes.IntPoint(k);
-      poly1d.CalcBasis(p, ip.x, shape_x);
-      poly1d.CalcBasis(p, ip.y, shape_y);
-      poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
-      o = 0;
-      for (int j = 0; j <= p; j++)
-      {
-         for (int i = 0; i + j <= p; i++)
-         {
-            Tt(k, o++) = shape_x[i]*shape_y[j]*shape_l[p-i-j];
-         }
-      }
-
-      poly1d.CalcBasis(q, ip.x, shape_x);
-      poly1d.CalcBasis(q, ip.y, shape_y);
-      poly1d.CalcBasis(q, 1. - ip.x - ip.y, shape_l);
-      const real_t b_T = ip.x * ip.y * (1 - ip.x - ip.y);
-      for (int j = 0; j <= q; j++)
-      {
-         for (int i = 0; i + j <= q; i++)
-         {
-            Tt(k, o++) = b_T*shape_x[i]*shape_y[j]*shape_l[q-i-j];
-         }
-      }
-   }
-
-   // Compute left inverse of T (given Tt = T^T).
-   DenseMatrix TtT(dof, dof);
-   MultAAt(Tt, TtT);
-
-   DenseMatrixInverse TtT_inv(TtT);
-   T_pinv.SetSize(dof, dof);
-   TtT_inv.Mult(Tt, T_pinv);
-}
-
-void H1Bubble_TriangleElement::CalcShape(const IntegrationPoint &ip,
-                                         Vector &shape) const
-{
-   const int p = base_order;
-   const int q = bubble_order;
-
-#ifdef MFEM_THREAD_SAFE
-   const int n1d = max(p + 1, q + 1);
-   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
-   Vector shape_x(n1d), shape_y(n1d), shape_l(n1d), u(npq);
-#endif
-
-   poly1d.CalcBasis(p, ip.x, shape_x);
-   poly1d.CalcBasis(p, ip.y, shape_y);
-   poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
-
-   int o = 0;
-   for (int j = 0; j <= p; j++)
-   {
-      for (int i = 0; i + j <= p; i++)
-      {
-         u(o++) = shape_x[i]*shape_y[j]*shape_l[p-i-j];
-      }
-   }
-
-   poly1d.CalcBasis(q, ip.x, shape_x);
-   poly1d.CalcBasis(q, ip.y, shape_y);
-   poly1d.CalcBasis(q, 1. - ip.x - ip.y, shape_l);
-   const real_t b_T = ip.x * ip.y * (1 - ip.x - ip.y);
-
-   for (int j = 0; j <= q; j++)
-   {
-      for (int i = 0; i + j <= q; i++)
-      {
-         u(o++) = b_T*shape_x[i]*shape_y[j]*shape_l[q-i-j];
-      }
-   }
-
-   T_pinv.Mult(u, shape);
-}
-
-void H1Bubble_TriangleElement::CalcDShape(const IntegrationPoint &ip,
-                                          DenseMatrix &dshape) const
-{
-   const int p = base_order;
-   const int q = bubble_order;
-
-#ifdef MFEM_THREAD_SAFE
-   const int n1d = max(p + 1, q + 1);
-   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
-   Vector  shape_x(n1d),  shape_y(n1d),  shape_l(n1d);
-   Vector dshape_x(n1d), dshape_y(n1d), dshape_l(n1d);
-   DenseMatrix du(npq, dim);
-#endif
-
-   const real_t lambda = 1.0 - ip.x - ip.y;
-
-   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
-   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
-   poly1d.CalcBasis(p, lambda, shape_l, dshape_l);
-
-   int o = 0;
-   for (int j = 0; j <= p; j++)
-   {
-      for (int i = 0; i + j <= p; i++)
-      {
-         int k = p - i - j;
-         du(o,0) = (dshape_x[i]*shape_l[k] - shape_x[i]*dshape_l[k])*shape_y[j];
-         du(o,1) = (dshape_y[j]* shape_l[k] - shape_y[j]*dshape_l[k])*shape_x[i];
-         o++;
-      }
-   }
-
-   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
-   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
-   poly1d.CalcBasis(q, lambda, shape_l, dshape_l);
-   const real_t b_T = ip.x * ip.y * lambda;
-   const real_t dxb_T = ip.y * (lambda - ip.x);
-   const real_t dyb_T = ip.x * (lambda - ip.y);
-
-   for (int j = 0; j <= q; j++)
-   {
-      for (int i = 0; i + j <= q; i++)
-      {
-         int k = q - i - j;
-         du(o,0) = shape_y[j]*(dxb_T*shape_x[i]*shape_l[k]
-                               + b_T*dshape_x[i]*shape_l[k]
-                               - b_T*shape_x[i]*dshape_l[k]);
-         du(o,1) = shape_x[i]*(dyb_T*shape_y[j]*shape_l[k]
-                               + b_T*dshape_y[j]*shape_l[k]
-                               - b_T*shape_y[j]*dshape_l[k]);
-         o++;
-      }
-   }
-
-   Mult(T_pinv, du, dshape);
-}
-
-H1Bubble_QuadrilateralElement::H1Bubble_QuadrilateralElement(
-   int p, int q, int btype)
-   : NodalFiniteElement(2, Geometry::SQUARE, 4*p + (q+1)*(q+1),
-                        max(p, 2 + q), FunctionSpace::Qk),
-     base_order(p), bubble_order(q)
-{
-   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
-   const real_t *cp2 = poly1d.ClosedPoints(
-                          q + 2, VerifyNodal(VerifyClosed(btype)));
-
-   const int n1d = max(p + 1, q + 1);
-   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
-
-#ifndef MFEM_THREAD_SAFE
-   shape_x.SetSize(n1d);
-   shape_y.SetSize(n1d);
-   dshape_x.SetSize(n1d);
-   dshape_y.SetSize(n1d);
-
-   u.SetSize(npq);
-   du.SetSize(npq, dim);
-#endif
-
-   // vertices
-   Nodes.IntPoint(0).Set2(cp[0], cp[0]);
-   Nodes.IntPoint(1).Set2(cp[p], cp[0]);
-   Nodes.IntPoint(2).Set2(cp[p], cp[p]);
-   Nodes.IntPoint(3).Set2(cp[0], cp[p]);
-
-   // edges
-   int o = 4;
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set2(cp[p], cp[i]);
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set2(cp[p-i], cp[p]);
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
-   }
-
-   // interior P_{q+2} nodes
-   for (int j = 1; j < q+2; j++)
-   {
-      for (int i = 1; i < q+2; i++)
-      {
-         Nodes.IntPoint(o++).Set2(cp2[i], cp2[j]);
-      }
-   }
-
-#ifdef MFEM_THREAD_SAFE
-   Vector shape_x(n1d), shape_y(n1d);
-#endif
-
-   DenseMatrix Tt(dof, npq);
-   for (int k = 0; k < dof; ++k)
-   {
-      const IntegrationPoint &ip = Nodes.IntPoint(k);
-      poly1d.CalcBasis(p, ip.x, shape_x);
-      poly1d.CalcBasis(p, ip.y, shape_y);
-
-      o = 0;
-      for (int j = 0; j <= p; j++)
-      {
-         for (int i = 0; i <= p; i++)
-         {
-            Tt(k, o++) = shape_x[i]*shape_y[j];
-         }
-      }
-
-      poly1d.CalcBasis(q, ip.x, shape_x);
-      poly1d.CalcBasis(q, ip.y, shape_y);
-      const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
-      for (int j = 0; j <= q; j++)
-      {
-         for (int i = 0; i <= q; i++)
-         {
-            Tt(k, o++) = b_T*shape_x[i]*shape_y[j];
-         }
-      }
-   }
-
-   // Compute left inverse of T (given Tt = T^T).
-   DenseMatrix TtT(dof, dof);
-   MultAAt(Tt, TtT);
-
-   DenseMatrixInverse TtT_inv(TtT);
-   T_pinv.SetSize(dof, dof);
-   TtT_inv.Mult(Tt, T_pinv);
-}
-
-void H1Bubble_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
-                                              Vector &shape) const
-{
-   const int p = base_order;
-   const int q = bubble_order;
-
-#ifdef MFEM_THREAD_SAFE
-   const int n1d = max(p + 1, q + 1);
-   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
-   Vector shape_x(n1d), shape_y(n1d), u(npq);
-#endif
-
-   poly1d.CalcBasis(p, ip.x, shape_x);
-   poly1d.CalcBasis(p, ip.y, shape_y);
-
-   int o = 0;
-   for (int j = 0; j <= p; j++)
-   {
-      for (int i = 0; i <= p; i++)
-      {
-         u(o++) = shape_x[i]*shape_y[j];
-      }
-   }
-
-   poly1d.CalcBasis(q, ip.x, shape_x);
-   poly1d.CalcBasis(q, ip.y, shape_y);
-   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
-
-   for (int j = 0; j <= q; j++)
-   {
-      for (int i = 0; i <= q; i++)
-      {
-         u(o++) = b_T*shape_x[i]*shape_y[j];
-      }
-   }
-
-   T_pinv.Mult(u, shape);
-}
-
-void H1Bubble_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
-                                               DenseMatrix &dshape) const
-{
-   const int p = base_order;
-   const int q = bubble_order;
-
-#ifdef MFEM_THREAD_SAFE
-   const int n1d = max(p + 1, q + 1);
-   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
-   Vector shape_x(n1d), shape_y(n1d), dshape_x(n1d), dshape_y(n1d);
-   DenseMatrix du(npq, dim);
-#endif
-
-   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
-   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
-
-   int o = 0;
-   for (int j = 0; j <= p; j++)
-   {
-      for (int i = 0; i <= p; i++)
-      {
-         du(o,0) = dshape_x[i]*shape_y[j];
-         du(o,1) = shape_x[i]*dshape_y[j];
-         o += 1;
-      }
-   }
-
-   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
-   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
-   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
-   const real_t dxb_T = (1.0 - 2*ip.x)*ip.y*(1.0 - ip.y);
-   const real_t dyb_T = ip.x*(1.0 - ip.x)*(1.0 - 2*ip.y);
-
-   for (int j = 0; j <= q; j++)
-   {
-      for (int i = 0; i <= q; i++)
-      {
-         du(o,0) = (dxb_T*shape_x[i] + b_T*dshape_x[i])*shape_y[j];
-         du(o,1) = (dyb_T*shape_y[j] + b_T*dshape_y[j])*shape_x[i];
-         o += 1;
-      }
-   }
-
-   Mult(T_pinv, du, dshape);
-}
-
-H1Bubble_TetrahedronElement::H1Bubble_TetrahedronElement(
-   int p, int q, int btype)
-   : NodalFiniteElement(3, Geometry::TETRAHEDRON,
-                        2*(p*p + 1) + ((q+1)*(q+2)*(q+3))/6,
-                        max(p, 4 + q), FunctionSpace::Pk),
-     base_order(p), bubble_order(q)
-{
-   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
-   const real_t *cp2 = poly1d.ClosedPoints(
-                          q + 4, VerifyNodal(VerifyClosed(btype)));
-
-   const int n1d = max(p+1, q+1);
-   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
-
-#ifndef MFEM_THREAD_SAFE
-   shape_x.SetSize(n1d);
-   shape_y.SetSize(n1d);
-   shape_z.SetSize(n1d);
-   shape_l.SetSize(n1d);
-   dshape_x.SetSize(n1d);
-   dshape_y.SetSize(n1d);
-   dshape_z.SetSize(n1d);
-   dshape_l.SetSize(n1d);
-   u.SetSize(npq);
-   du.SetSize(npq, dim);
-#else
-   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d);
-#endif
-
-   // vertices
-   Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
-   Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
-   Nodes.IntPoint(2).Set3(cp[0], cp[p], cp[0]);
-   Nodes.IntPoint(3).Set3(cp[0], cp[0], cp[p]);
-
-   // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
-   int o = 4;
-   for (int i = 1; i < p; i++)  // (0,1)
-   {
-      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
-   }
-   for (int i = 1; i < p; i++)  // (0,2)
-   {
-      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
-   }
-   for (int i = 1; i < p; i++)  // (0,3)
-   {
-      Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
-   }
-   for (int i = 1; i < p; i++)  // (1,2)
-   {
-      Nodes.IntPoint(o++).Set3(cp[p-i], cp[i], cp[0]);
-   }
-   for (int i = 1; i < p; i++)  // (1,3)
-   {
-      Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
-   }
-   for (int i = 1; i < p; i++)  // (2,3)
-   {
-      Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
-   }
-
-   // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i + j < p; i++)  // (1,2,3)
-      {
-         real_t w = cp[i] + cp[j] + cp[p-i-j];
-         Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, cp[i]/w, cp[j]/w);
-      }
-   }
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i + j < p; i++)  // (0,3,2)
-      {
-         real_t w = cp[i] + cp[j] + cp[p-i-j];
-         Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
-      }
-   }
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i + j < p; i++)  // (0,1,3)
-      {
-         real_t w = cp[i] + cp[j] + cp[p-i-j];
-         Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
-      }
-   }
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i + j < p; i++)  // (0,2,1)
-      {
-         real_t w = cp[i] + cp[j] + cp[p-i-j];
-         Nodes.IntPoint(o++).Set3(cp[j]/w, cp[i]/w, cp[0]);
-      }
-   }
-
-   // Interior P_{q+4} nodes
-   for (int k = 1; k < q + 4; k++)
-   {
-      for (int j = 1; j + k < q + 4; j++)
-      {
-         for (int i = 1; i + j + k < q + 4; i++)
-         {
-            real_t w = cp2[i] + cp2[j] + cp2[k] + cp2[q+4-i-j-k];
-            Nodes.IntPoint(o++).Set3(cp2[i]/w, cp2[j]/w, cp2[k]/w);
-         }
-      }
-   }
-
-   DenseMatrix Tt(dof, npq);
-   for (int m = 0; m < dof; ++m)
-   {
-      const IntegrationPoint &ip = Nodes.IntPoint(m);
-      poly1d.CalcBasis(p, ip.x, shape_x);
-      poly1d.CalcBasis(p, ip.y, shape_y);
-      poly1d.CalcBasis(p, ip.z, shape_z);
-      poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
-
-      o = 0;
-      for (int k = 0; k <= p; k++)
-      {
-         for (int j = 0; j + k <= p; j++)
-         {
-            for (int i = 0; i + j + k <= p; i++)
-            {
-               Tt(m, o++) = shape_x[i]*shape_y[j]*shape_z[k]*shape_l[p-i-j-k];
-            }
-         }
-      }
-
-      poly1d.CalcBasis(q, ip.x, shape_x);
-      poly1d.CalcBasis(q, ip.y, shape_y);
-      poly1d.CalcBasis(q, ip.z, shape_z);
-      poly1d.CalcBasis(q, 1. - ip.x - ip.y - ip.z, shape_l);
-      const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
-
-      for (int k = 0; k <= q; k++)
-      {
-         for (int j = 0; j + k <= q; j++)
-         {
-            for (int i = 0; i + j + k <= q; i++)
-            {
-               Tt(m, o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k]*shape_l[q-i-j-k];
-            }
-         }
-      }
-   }
-
-   // Compute left inverse of T (given Tt = T^T).
-   DenseMatrix TtT(dof, dof);
-   MultAAt(Tt, TtT);
-
-   DenseMatrixInverse TtT_inv(TtT);
-   T_pinv.SetSize(dof, dof);
-   TtT_inv.Mult(Tt, T_pinv);
-}
-
-void H1Bubble_TetrahedronElement::CalcShape(const IntegrationPoint &ip,
-                                            Vector &shape) const
-{
-   const int p = base_order;
-   const int q = bubble_order;
-
-#ifdef MFEM_THREAD_SAFE
-   const int n1d = max(p + 1, q + 1);
-   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
-   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d), u(npq);
-#endif
-
-   poly1d.CalcBasis(p, ip.x, shape_x);
-   poly1d.CalcBasis(p, ip.y, shape_y);
-   poly1d.CalcBasis(p, ip.z, shape_z);
-   poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
-
-   int o = 0;
-   for (int k = 0; k <= p; k++)
-   {
-      for (int j = 0; j + k <= p; j++)
-      {
-         for (int i = 0; i + j + k <= p; i++)
-         {
-            u[o++] = shape_x[i]*shape_y[j]*shape_z[k]*shape_l[p-i-j-k];
-         }
-      }
-   }
-
-   poly1d.CalcBasis(q, ip.x, shape_x);
-   poly1d.CalcBasis(q, ip.y, shape_y);
-   poly1d.CalcBasis(q, ip.z, shape_z);
-   poly1d.CalcBasis(q, 1. - ip.x - ip.y - ip.z, shape_l);
-   const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
-
-   for (int k = 0; k <= q; k++)
-   {
-      for (int j = 0; j + k <= q; j++)
-      {
-         for (int i = 0; i + j + k <= q; i++)
-         {
-            u(o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k]*shape_l[q-i-j-k];
-         }
-      }
-   }
-
-   T_pinv.Mult(u, shape);
-}
-
-void H1Bubble_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
-                                             DenseMatrix &dshape) const
-{
-   const int p = base_order;
-   const int q = bubble_order;
-
-#ifdef MFEM_THREAD_SAFE
-   const int n1d = max(p+1, q+1);
-   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
-
-   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d);
-   Vector dshape_x(n1d), dshape_y(n1d), dshape_z(n1d), dshape_l(n1d);
-   DenseMatrix du(npq, dim);
-#endif
-
-   const real_t lambda = 1.0 - ip.x - ip.y - ip.z;
-
-   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
-   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
-   poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
-   poly1d.CalcBasis(p, lambda, shape_l, dshape_l);
-
-   int o = 0;
-   for (int k = 0; k <= p; k++)
-   {
-      for (int j = 0; j + k <= p; j++)
-      {
-         for (int i = 0; i + j + k <= p; i++)
-         {
-            int l = p - i - j - k;
-            du(o,0) = (dshape_x[i]*shape_l[l] - shape_x[i]*dshape_l[l])
-                      *shape_y[j]*shape_z[k];
-            du(o,1) = (dshape_y[j]*shape_l[l] - shape_y[j]*dshape_l[l])
-                      *shape_x[i]*shape_z[k];
-            du(o,2) = (dshape_z[k]*shape_l[l] - shape_z[k]*dshape_l[l])
-                      *shape_x[i]*shape_y[j];
-            o++;
-         }
-      }
-   }
-
-   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
-   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
-   poly1d.CalcBasis(q, ip.z, shape_z, dshape_z);
-   poly1d.CalcBasis(q, lambda, shape_l, dshape_l);
-   const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
-   const real_t dxb_T = ip.y * ip.z * (lambda - ip.x);
-   const real_t dyb_T = ip.x * ip.z * (lambda - ip.y);
-   const real_t dzb_T = ip.x * ip.y * (lambda - ip.z);
-
-   for (int k = 0; k <= q; k++)
-   {
-      for (int j = 0; j + k <= q; j++)
-      {
-         for (int i = 0; i + j + k <= q; i++)
-         {
-            int l = q - i - j - k;
-            du(o,0) = shape_y[j]*shape_z[k]*(dxb_T*shape_x[i]*shape_l[l]
-                                             + b_T*dshape_x[i]*shape_l[l]
-                                             - b_T*shape_x[i]*dshape_l[l]);
-            du(o,1) = shape_x[i]*shape_z[k]*(dyb_T*shape_y[j]*shape_l[l]
-                                             + b_T*dshape_y[j]*shape_l[l]
-                                             - b_T*shape_y[j]*dshape_l[l]);
-            du(o,2) = shape_x[i]*shape_y[j]*(dzb_T*shape_z[k]*shape_l[l]
-                                             + b_T*dshape_z[k]*shape_l[l]
-                                             - b_T*shape_z[k]*dshape_l[l]);
-            o++;
-         }
-      }
-   }
-
-   Mult(T_pinv, du, dshape);
-}
-
-H1Bubble_HexahedronElement::H1Bubble_HexahedronElement(
-   int p, int q, int btype)
-   : NodalFiniteElement(3, Geometry::CUBE, (2 + 6*p*p) + (q+1)*(q+1)*(q+1),
-                        max(p, 2 + q), FunctionSpace::Qk),
-     base_order(p), bubble_order(q)
-{
-   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
-   const real_t *cp2 = poly1d.ClosedPoints(
-                          q + 2, VerifyNodal(VerifyClosed(btype)));
-
-   const int n1d = max(p + 1, q + 1);
-   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
-
-#ifndef MFEM_THREAD_SAFE
-   shape_x.SetSize(n1d);
-   shape_y.SetSize(n1d);
-   shape_z.SetSize(n1d);
-   dshape_x.SetSize(n1d);
-   dshape_y.SetSize(n1d);
-   dshape_z.SetSize(n1d);
-
-   u.SetSize(npq);
-   du.SetSize(npq, dim);
-#endif
-
-   // vertices
-   Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
-   Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
-   Nodes.IntPoint(2).Set3(cp[p], cp[p], cp[0]);
-   Nodes.IntPoint(3).Set3(cp[0], cp[p], cp[0]);
-
-   Nodes.IntPoint(4).Set3(cp[0], cp[0], cp[p]);
-   Nodes.IntPoint(5).Set3(cp[p], cp[0], cp[p]);
-   Nodes.IntPoint(6).Set3(cp[p], cp[p], cp[p]);
-   Nodes.IntPoint(7).Set3(cp[0], cp[p], cp[p]);
-
-   int o = 8;
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);   // (0,1)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[0]);   // (1,2)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[0]);   // (3,2)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);   // (0,3)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[p]);   // (4,5)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[p]);   // (5,6)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[p]);   // (7,6)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[p]);   // (4,7)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);   // (0,4)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[p], cp[0], cp[i]);   // (1,5)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[p], cp[p], cp[i]);   // (2,6)
-   }
-   for (int i = 1; i < p; i++)
-   {
-      Nodes.IntPoint(o++).Set3(cp[0], cp[p], cp[i]);   // (3,7)
-   }
-
-   // faces
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i < p; i++)
-      {
-         Nodes.IntPoint(o++).Set3(cp[i], cp[p-j], cp[0]);   // (3,2,1,0)
-      }
-   }
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i < p; i++)
-      {
-         Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[j]);   // (0,1,5,4)
-      }
-   }
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i < p; i++)
-      {
-         Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[j]);   // (1,2,6,5)
-      }
-   }
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i < p; i++)
-      {
-         Nodes.IntPoint(o++).Set3(cp[p-i], cp[p], cp[j]);   // (2,3,7,6)
-      }
-   }
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i < p; i++)
-      {
-         Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[j]);   // (3,0,4,7)
-      }
-   }
-   for (int j = 1; j < p; j++)
-   {
-      for (int i = 1; i < p; i++)
-      {
-         Nodes.IntPoint(o++).Set3(cp[i], cp[j], cp[p]);   // (4,5,6,7)
-      }
-   }
-
-   // interior P_{q+2} nodes
-   for (int k = 1; k < q+2; k++)
-   {
-      for (int j = 1; j < q+2; j++)
-      {
-         for (int i = 1; i < q+2; i++)
-         {
-            Nodes.IntPoint(o++).Set3(cp2[i], cp2[j], cp2[k]);
-         }
-      }
-   }
-
-#ifdef MFEM_THREAD_SAFE
-   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d);
-#endif
-
-   DenseMatrix Tt(dof, npq);
-   for (int m = 0; m < dof; ++m)
-   {
-      const IntegrationPoint &ip = Nodes.IntPoint(m);
-      poly1d.CalcBasis(p, ip.x, shape_x);
-      poly1d.CalcBasis(p, ip.y, shape_y);
-      poly1d.CalcBasis(p, ip.z, shape_z);
-
-      o = 0;
-      for (int k = 0; k <= p; k++)
-      {
-         for (int j = 0; j <= p; j++)
-         {
-            for (int i = 0; i <= p; i++)
-            {
-               Tt(m, o++) = shape_x[i]*shape_y[j]*shape_z[k];
-            }
-         }
-      }
-
-      poly1d.CalcBasis(q, ip.x, shape_x);
-      poly1d.CalcBasis(q, ip.y, shape_y);
-      poly1d.CalcBasis(q, ip.z, shape_z);
-      const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
-      for (int k = 0; k <= q; k++)
-      {
-         for (int j = 0; j <= q; j++)
-         {
-            for (int i = 0; i <= q; i++)
-            {
-               Tt(m, o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k];
-            }
-         }
-      }
-   }
-
-   // Compute left inverse of T (given Tt = T^T).
-   DenseMatrix TtT(dof, dof);
-   MultAAt(Tt, TtT);
-
-   DenseMatrixInverse TtT_inv(TtT);
-   T_pinv.SetSize(dof, dof);
-   TtT_inv.Mult(Tt, T_pinv);
-}
-
-void H1Bubble_HexahedronElement::CalcShape(const IntegrationPoint &ip,
-                                           Vector &shape) const
-{
-   const int p = base_order;
-   const int q = bubble_order;
-
-#ifdef MFEM_THREAD_SAFE
-   const int n1d = max(p + 1, q + 1);
-   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
-   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), u(npq);
-#endif
-
-   poly1d.CalcBasis(p, ip.x, shape_x);
-   poly1d.CalcBasis(p, ip.y, shape_y);
-   poly1d.CalcBasis(p, ip.z, shape_z);
-
-   int o = 0;
-   for (int k = 0; k <= p; k++)
-   {
-      for (int j = 0; j <= p; j++)
-      {
-         for (int i = 0; i <= p; i++)
-         {
-            u(o++) = shape_x[i]*shape_y[j]*shape_z[k];
-         }
-      }
-   }
-
-   poly1d.CalcBasis(q, ip.x, shape_x);
-   poly1d.CalcBasis(q, ip.y, shape_y);
-   poly1d.CalcBasis(q, ip.z, shape_z);
-   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
-
-   for (int k = 0; k <= q; k++)
-   {
-      for (int j = 0; j <= q; j++)
-      {
-         for (int i = 0; i <= q; i++)
-         {
-            u(o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k];
-         }
-      }
-   }
-
-   T_pinv.Mult(u, shape);
-}
-
-void H1Bubble_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
-                                            DenseMatrix &dshape) const
-{
-   const int p = base_order;
-   const int q = bubble_order;
-
-#ifdef MFEM_THREAD_SAFE
-   const int n1d = max(p + 1, q + 1);
-   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
-   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), dshape_x(n1d),
-          dshape_y(n1d), dshape_z(n1d);
-   DenseMatrix du(npq, dim);
-#endif
-
-   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
-   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
-   poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
-
-   int o = 0;
-   for (int k = 0; k <= p; k++)
-   {
-      for (int j = 0; j <= p; j++)
-      {
-         for (int i = 0; i <= p; i++)
-         {
-            du(o,0) = dshape_x[i]*shape_y[j]*shape_z[k];
-            du(o,1) = shape_x[i]*dshape_y[j]*shape_z[k];
-            du(o,2) = shape_x[i]*shape_y[j]*dshape_z[k];
-            o += 1;
-         }
-      }
-   }
-
-   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
-   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
-   poly1d.CalcBasis(q, ip.z, shape_z, dshape_z);
-   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
-   const real_t dxb_T = (1.0 - 2*ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
-   const real_t dyb_T = ip.x*(1.0 - ip.x)*(1.0 - 2*ip.y)*ip.z*(1.0 - ip.z);
-   const real_t dzb_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*(1.0 - 2*ip.z);
-
-   for (int k = 0; k <= q; k++)
-   {
-      for (int j = 0; j <= q; j++)
-      {
-         for (int i = 0; i <= q; i++)
-         {
-            du(o,0) = (dxb_T*shape_x[i] + b_T*dshape_x[i])*shape_y[j]*shape_z[k];
-            du(o,1) = (dyb_T*shape_y[j] + b_T*dshape_y[j])*shape_x[i]*shape_z[k];
-            du(o,2) = (dzb_T*shape_z[k] + b_T*dshape_z[k])*shape_x[i]*shape_y[j];
-            o += 1;
-         }
-      }
-   }
-
-   Mult(T_pinv, du, dshape);
-}
-
-}

fem/fe/fe_h1_bubble.hpp

@@ -1,109 +0,0 @@
-// 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.
-
-#ifndef MFEM_FE_H1_BUBBLE
-#define MFEM_FE_H1_BUBBLE
-
-#include "fe_base.hpp"
-
-namespace mfem
-{
-
-/// Arbitrary order H1 plus bubble elements in 2D on a triangle
-class H1Bubble_TriangleElement : public NodalFiniteElement
-{
-private:
-#ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
-   mutable DenseMatrix du;
-#endif
-   int base_order;
-   int bubble_order;
-   DenseMatrix T_pinv;
-
-public:
-   /// @brief Construct the triangular bubble element with degree-p polynomials,
-   /// enriched with cubic bubble times degree q polynomial.
-   H1Bubble_TriangleElement(int p, int q, int btype = BasisType::GaussLobatto);
-   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
-   void CalcDShape(const IntegrationPoint &ip,
-                   DenseMatrix &dshape) const override;
-};
-
-/// Arbitrary order H1 plus bubble elements in 2D on a quadrilateral
-class H1Bubble_QuadrilateralElement : public NodalFiniteElement
-{
-private:
-#ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, shape_y, dshape_x, dshape_y, u;
-   mutable DenseMatrix du;
-#endif
-   int base_order;
-   int bubble_order;
-   DenseMatrix T_pinv;
-
-public:
-   /// @brief Construct the quadrilateral bubble element with degree-p
-   /// polynomials, enriched with biquadratic bubble times degree q polynomial.
-   H1Bubble_QuadrilateralElement(
-      int p, int q, int btype = BasisType::GaussLobatto);
-   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
-   void CalcDShape(const IntegrationPoint &ip,
-                   DenseMatrix &dshape) const override;
-};
-
-/// Arbitrary order H1 plus bubble elements in 3D on a tetrahedron
-class H1Bubble_TetrahedronElement : public NodalFiniteElement
-{
-private:
-#ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, shape_y, shape_z, shape_l;
-   mutable Vector dshape_x, dshape_y, dshape_z, dshape_l, u;
-   mutable DenseMatrix du;
-#endif
-   int base_order;
-   int bubble_order;
-   DenseMatrix T_pinv;
-
-public:
-   /// @brief Construct the tetrahedral bubble element with degree-p
-   /// polynomials, enriched with quartic bubble times degree q polynomial.
-   H1Bubble_TetrahedronElement(int p, int q, int btype = BasisType::GaussLobatto);
-   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
-   void CalcDShape(const IntegrationPoint &ip,
-                   DenseMatrix &dshape) const override;
-};
-
-/// Arbitrary order H1 plus bubble elements in 3D on a hexahedron
-class H1Bubble_HexahedronElement : public NodalFiniteElement
-{
-private:
-#ifndef MFEM_THREAD_SAFE
-   mutable Vector shape_x, shape_y, shape_z;
-   mutable Vector dshape_x, dshape_y, dshape_z, u;
-   mutable DenseMatrix du;
-#endif
-   int base_order;
-   int bubble_order;
-   DenseMatrix T_pinv;
-
-public:
-   /// @brief Construct the hexahedral bubble element with degree-p polynomials,
-   /// enriched with triquadratic bubble times degree q polynomial.
-   H1Bubble_HexahedronElement(int p, int q, int btype = BasisType::GaussLobatto);
-   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
-   void CalcDShape(const IntegrationPoint &ip,
-                   DenseMatrix &dshape) const override;
-};
-
-} // namespace mfem
-
-#endif

fem/fe_coll.cpp

@@ -243,21 +243,11 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
    {
       fec = new H1Ser_FECollection(atoi(name + 10), atoi(name + 6));
    }
-   else if (!strncmp(name, "H1Bubble_", 9))
-   {
-      fec = new H1Bubble_FECollection(atoi(name + 13), atoi(name + 16),
-                                      atoi(name + 9));
-   }
    else if (!strncmp(name, "H1@", 3))
    {
       fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
                                 BasisType::GetType(name[3]));
    }
-   else if (!strncmp(name, "H1Bubble@", 9))
-   {
-      fec = new H1Bubble_FECollection(atoi(name + 15), atoi(name + 18),
-                                      atoi(name + 11), BasisType::GetType(name[9]));
-   }
    else if (!strncmp(name, "L2_T", 4))
       fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
                                 atoi(name + 4));
@@ -2132,173 +2122,6 @@ H1_FECollection::~H1_FECollection()
    }
 }
 
-static int GetBubbleSpaceOrder(int p, int q, int dim)
-{
-   switch (dim)
-   {
-      case 0: return 0;
-      case 1: return std::max(p, q + 2);
-      case 2: return std::max(p, q + 3);
-      case 3: return std::max(p, q + 4);
-   }
-   MFEM_ABORT("Unsupported dimension.");
-}
-
-H1Bubble_FECollection::H1Bubble_FECollection(const int p, const int q,
-                                             const int dim, const int btype)
-   : FiniteElementCollection(GetBubbleSpaceOrder(p, q, dim)),
-     dim(dim),
-     b_type(BasisType::Check(btype)),
-     h1_order(p),
-     bubble_order(q)
-{
-   MFEM_VERIFY(p >= 1, "H1Bubble_FECollection requires order >= 1.");
-   MFEM_VERIFY(dim >= 0 && dim <= 3, "Unsupported dimension.");
-
-   switch (btype)
-   {
-      case BasisType::GaussLobatto:
-      {
-         snprintf(fec_name, 32, "H1Bubble_%dD_P%d_P%d", dim, p, q);
-         break;
-      }
-      default:
-      {
-         const int pt_type = BasisType::GetQuadrature1D(btype);
-         MFEM_VERIFY(Quadrature1D::CheckClosed(pt_type) != Quadrature1D::Invalid,
-                     "unsupported BasisType: " << BasisType::Name(btype));
-         snprintf(fec_name, 32, "H1Bubble@%c_%dD_P%d_P%d",
-                  (int)BasisType::GetChar(btype), dim, p, q);
-      }
-   }
-
-   dofs[Geometry::POINT] = 1;
-   elements[Geometry::POINT] = make_unique<PointFiniteElement>();
-
-   if (dim >= 1)
-   {
-      dofs[Geometry::SEGMENT] = p - 1;
-      elements[Geometry::SEGMENT] = make_unique<H1_SegmentElement>(p, btype);
-   }
-
-   if (dim == 2)
-   {
-      dofs[Geometry::TRIANGLE] = ((q+1)*(q+2))/2;
-      dofs[Geometry::SQUARE] = (q+1)*(q+1);
-
-      elements[Geometry::TRIANGLE] =
-         make_unique<H1Bubble_TriangleElement>(p, q, btype);
-      elements[Geometry::SQUARE] =
-         make_unique<H1Bubble_QuadrilateralElement>(p, q, btype);
-   }
-
-   if (dim == 3)
-   {
-      dofs[Geometry::TRIANGLE] = ((p-1)*(p-2))/2;
-      dofs[Geometry::SQUARE] = (p-1)*(p-1);
-      dofs[Geometry::TETRAHEDRON] = ((q+1)*(q+2)*(q+3))/6;
-      dofs[Geometry::CUBE] = (q+1)*(q+1)*(q+1);
-
-      elements[Geometry::TRIANGLE] = make_unique<H1_TriangleElement>(p, btype);
-      elements[Geometry::SQUARE] = make_unique<H1_QuadrilateralElement>(p, btype);
-
-      elements[Geometry::TETRAHEDRON] =
-         make_unique<H1Bubble_TetrahedronElement>(p, q, btype);
-      elements[Geometry::CUBE] =
-         make_unique<H1Bubble_HexahedronElement>(p, q, btype);
-   }
-
-   // DOF orderings. Need only for lower-dimensional entities.
-   // Segment DOF orderings in 2D.
-   if (dim >= 2)
-   {
-      seg_dof_ord[0].resize(p - 1);
-      seg_dof_ord[1].resize(p - 1);
-      for (int i = 0; i < p - 1; i++)
-      {
-         seg_dof_ord[0][i] = i;
-         seg_dof_ord[1][i] = p - 2 - i;
-      }
-   }
-
-   // Face (triangle or quadrilateral) DOF orderings in 3D.
-   if (dim == 3)
-   {
-      const int n_tri_dof = dofs[Geometry::TRIANGLE];
-      for (int i = 0; i < 6; i++)
-      {
-         tri_dof_ord[i].resize(n_tri_dof);
-      }
-      // see Mesh::GetTriOrientation in mesh/mesh.cpp
-      const int pm1 = p - 1;
-      const int pm2 = p - 2;
-      for (int j = 0; j < pm2; j++)
-      {
-         for (int i = 0; i + j < pm2; i++)
-         {
-            int o = n_tri_dof - ((pm1 - j)*(pm2 - j))/2 + i;
-            int k = (p - 3) - j - i;
-            tri_dof_ord[0][o] = o;  // (0,1,2)
-            tri_dof_ord[1][o] = n_tri_dof - ((pm1-j)*(pm2-j))/2 + k;  // (1,0,2)
-            tri_dof_ord[2][o] = n_tri_dof - ((pm1-i)*(pm2-i))/2 + k;  // (2,0,1)
-            tri_dof_ord[3][o] = n_tri_dof - ((pm1-k)*(pm2-k))/2 + i;  // (2,1,0)
-            tri_dof_ord[4][o] = n_tri_dof - ((pm1-k)*(pm2-k))/2 + j;  // (1,2,0)
-            tri_dof_ord[5][o] = n_tri_dof - ((pm1-i)*(pm2-i))/2 + j;  // (0,2,1)
-         }
-      }
-
-      const int n_quad_dof = dofs[Geometry::SQUARE];
-      for (int i = 0; i < 8; i++)
-      {
-         quad_dof_ord[i].resize(n_quad_dof);
-      }
-      for (int j = 0; j < pm1; j++)
-      {
-         for (int i = 0; i < pm1; i++)
-         {
-            int o = i + j*pm1;
-            quad_dof_ord[0][o] = i + j*pm1;  // (0,1,2,3)
-            quad_dof_ord[1][o] = j + i*pm1;  // (0,3,2,1)
-            quad_dof_ord[2][o] = j + (pm2 - i)*pm1;  // (1,2,3,0)
-            quad_dof_ord[3][o] = (pm2 - i) + j*pm1;  // (1,0,3,2)
-            quad_dof_ord[4][o] = (pm2 - i) + (pm2 - j)*pm1;  // (2,3,0,1)
-            quad_dof_ord[5][o] = (pm2 - j) + (pm2 - i)*pm1;  // (2,1,0,3)
-            quad_dof_ord[6][o] = (pm2 - j) + i*pm1;  // (3,0,1,2)
-            quad_dof_ord[7][o] = i + (pm2 - j)*pm1;  // (3,2,1,0)
-         }
-      }
-   }
-}
-
-const FiniteElement *
-H1Bubble_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
-{
-   return elements[GeomType].get();
-}
-
-const int *H1Bubble_FECollection::DofOrderForOrientation(
-   Geometry::Type GeomType, int Or) const
-{
-   if (GeomType == Geometry::SEGMENT)
-   {
-      return (Or > 0) ? seg_dof_ord[0].data() : seg_dof_ord[1].data();
-   }
-   else if (GeomType == Geometry::TRIANGLE)
-   {
-      return tri_dof_ord[Or%6].data();
-   }
-   else if (GeomType == Geometry::SQUARE)
-   {
-      return quad_dof_ord[Or%8].data();
-   }
-   return nullptr;
-}
-
-FiniteElementCollection *H1Bubble_FECollection::GetTraceCollection() const
-{
-   return (dim < 0) ? NULL : new H1_Trace_FECollection(h1_order, dim, b_type);
-}
-
 
 H1_Trace_FECollection::H1_Trace_FECollection(const int p, const int dim,
                                              const int btype)

fem/fe_coll.hpp

@@ -111,8 +111,6 @@ public:
    | :------: | :---: | :---: | :-------: | :-----: | :---: |
    | H1_[DIM]_[ORDER] | H1 | * | 1 | VALUE | H1 nodal elements |
    | H1@[BTYPE]_[DIM]_[ORDER] | H1 | * | * | VALUE | H1 nodal elements |
-   | H1Bubble_[DIM]_[ORDER]_[BUBBLE_ORDER] | H1 | * | 1 | VALUE | H1 nodal elements enriched with bubble functions |
-   | H1Bubble@[BTYPE]_[DIM]_[ORDER]_[BUBBLE_ORDER] | H1 | * | 1 | VALUE | H1 nodal elements enriched with bubble functions |
    | H1Pos_[DIM]_[ORDER] | H1 | * | 2 | VALUE | H1 nodal elements |
    | H1Pos_Trace_[DIM]_[ORDER] | H^{1/2} | * | 2 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
    | H1_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
@@ -319,59 +317,6 @@ public:
    virtual ~H1_FECollection();
 };
 
-/// @brief Arbitrary order $H^1$-conforming (continuous) finite elements
-/// enriched with bubble functions.
-///
-/// The bubble space consists of the standard $P_p$ or $Q_p$ space, enriched
-/// with bubble functions, which are degree-$q$ polynomials times $b$, where $b$
-/// is the lowest-order bubble function.
-///
-/// The traces are the same as the standard $H^1$ traces.
-class H1Bubble_FECollection : public FiniteElementCollection
-{
-protected:
-   int dim;
-   int b_type;
-   int h1_order;
-   int bubble_order;
-
-   char fec_name[32];
-   std::array<int, Geometry::NumGeom> dofs{}; // zero initialize
-   std::array<std::unique_ptr<FiniteElement>, Geometry::NumGeom> elements;
-
-   std::array<std::vector<int>, 2> seg_dof_ord;
-   std::array<std::vector<int>, 6> tri_dof_ord;
-   std::array<std::vector<int>, 8> quad_dof_ord;
-   std::array<std::vector<int>, 24> tet_dof_ord;
-
-public:
-   /// Construct the $H^1$ bubble collection consisting of degree-$p$
-   /// polynomials enriched with the bubble function times degree-$q$
-   /// polynomials.
-   explicit H1Bubble_FECollection(const int p, const int q, const int dim = 3,
-                                  const int btype = BasisType::GaussLobatto);
-
-   const FiniteElement *
-   FiniteElementForGeometry(Geometry::Type GeomType) const override;
-
-   int DofForGeometry(Geometry::Type GeomType) const override
-   { return dofs[GeomType]; }
-
-   const int *DofOrderForOrientation(Geometry::Type GeomType,
-                                     int Or) const override;
-
-   const char *Name() const override { return fec_name; }
-
-   int GetContType() const override { return CONTINUOUS; }
-
-   int GetBasisType() const { return b_type; }
-
-   FiniteElementCollection *GetTraceCollection() const override;
-
-   FiniteElementCollection *Clone(int p) const override
-   { return new H1Bubble_FECollection(p, bubble_order, dim, b_type); }
-};
-
 /** @brief Arbitrary order H1-conforming (continuous) finite elements with
     positive basis functions. */
 class H1Pos_FECollection : public H1_FECollection

fem/fespace.cpp

@@ -3877,12 +3877,9 @@ const FiniteElement *FiniteElementSpace::GetFE(int i) const
    else
    {
 #ifdef MFEM_DEBUG
-      // Consistency check: fec->GetOrder() and FE->GetOrder() should return
-      // the same value (for standard, constant-order spaces). Skip this check
-      // even for constant-order bubble spaces, since the bubble functions on
-      // different geometries have different orders.
-      if (!IsVariableOrder() && FE->GetDim() > 0 &&
-          dynamic_cast<const H1Bubble_FECollection*>(fec) == nullptr)
+      // consistency check: fec->GetOrder() and FE->GetOrder() should return
+      // the same value (for standard, constant-order spaces)
+      if (!IsVariableOrder() && FE->GetDim() > 0)
       {
          MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
                      "internal error: " <<

fem/gridfunc.cpp

@@ -3137,29 +3137,52 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
 }
 
 void GridFunction::ProjectBdrCoefficientNormal(
-   Coefficient *coeff, VectorCoefficient *vcoeff, const Array<int> &bdr_attr)
+   VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
 {
-   if (fes->GetNBE() > 0)
-   {
-      // TODO: Replace this by GetTypicalBdrElement() once implemented
-      const FiniteElement *be = fes->GetBE(0);
-      MFEM_VERIFY(be->GetRangeType() == FiniteElement::SCALAR &&
-                  be->GetMapType() == FiniteElement::INTEGRAL, "Not an RT FE space!");
-   }
+#if 0
+   // implementation for the case when the face dofs are integrals of the
+   // normal component.
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   Array<int> dofs;
+   int dim = vcoeff.GetVDim();
+   Vector vc(dim), nor(dim), lvec, shape;
 
+   for (int i = 0; i < fes->GetNBE(); i++)
+   {
+      if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
+      {
+         continue;
+      }
+      fe = fes->GetBE(i);
+      T = fes->GetBdrElementTransformation(i);
+      int intorder = 2*fe->GetOrder(); // !!!
+      const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intorder);
+      int nd = fe->GetDof();
+      lvec.SetSize(nd);
+      shape.SetSize(nd);
+      lvec = 0.0;
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(j);
+         T->SetIntPoint(&ip);
+         vcoeff.Eval(vc, *T, ip);
+         CalcOrtho(T->Jacobian(), nor);
+         fe->CalcShape(ip, shape);
+         lvec.Add(ip.weight * (vc * nor), shape);
+      }
+      fes->GetBdrElementDofs(i, dofs);
+      SetSubVector(dofs, lvec);
+   }
+#else
    // implementation for the case when the face dofs are scaled point
    // values of the normal component.
    const FiniteElement *fe;
    ElementTransformation *T;
    Array<int> dofs;
-   Vector vc, nor, lvec;
+   int dim = vcoeff.GetVDim();
+   Vector vc(dim), nor(dim), lvec;
    DofTransformation doftrans;
-   if (vcoeff)
-   {
-      const int dim = vcoeff->GetVDim();
-      vc.SetSize(dim);
-      nor.SetSize(dim);
-   }
 
    for (int i = 0; i < fes->GetNBE(); i++)
    {
@@ -3175,22 +3198,15 @@ void GridFunction::ProjectBdrCoefficientNormal(
       {
          const IntegrationPoint &ip = ir.IntPoint(j);
          T->SetIntPoint(&ip);
-         if (coeff)
-         {
-            const real_t c = coeff->Eval(*T, ip);
-            lvec(j) = c * T->Weight();
-         }
-         else if (vcoeff)
-         {
-            vcoeff->Eval(vc, *T, ip);
-            CalcOrtho(T->Jacobian(), nor);
-            lvec(j) = (vc * nor);
-         }
+         vcoeff.Eval(vc, *T, ip);
+         CalcOrtho(T->Jacobian(), nor);
+         lvec(j) = (vc * nor);
       }
       fes->GetBdrElementDofs(i, dofs, doftrans);
       doftrans.TransformPrimal(lvec);
       SetSubVector(dofs, lvec);
    }
+#endif
 }
 
 void GridFunction::ProjectBdrCoefficientTangent(
@@ -4991,14 +5007,6 @@ real_t ExtrudeCoefficient::Eval(ElementTransformation &T,
    return sol_in.Eval(*T_in, ip);
 }
 
-void VectorExtrudeCoefficient::Eval(Vector &v, ElementTransformation &T,
-                                    const IntegrationPoint &ip)
-{
-   ElementTransformation *T_in =
-      mesh_in->GetElementTransformation(T.ElementNo / n);
-   T_in->SetIntPoint(&ip);
-   sol_in.Eval(v, *T_in, ip);
-}
 
 GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
                                     GridFunction *sol, const int ny)
@@ -5049,17 +5057,10 @@ GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
       return NULL;
    }
    FiniteElementSpace *solfes2d;
-   const int vdim = sol->FESpace()->GetVDim();
-   solfes2d = new FiniteElementSpace(mesh2d, solfec2d, vdim);
+   // assuming sol is scalar
+   solfes2d = new FiniteElementSpace(mesh2d, solfec2d);
    sol2d = new GridFunction(solfes2d);
    sol2d->MakeOwner(solfec2d);
-   if (vdim > 1)
-   {
-      VectorGridFunctionCoefficient vcsol(sol);
-      VectorExtrudeCoefficient vc2d(mesh, vcsol, ny);
-      sol2d->ProjectCoefficient(vc2d);
-   }
-   else
    {
       GridFunctionCoefficient csol(sol);
       ExtrudeCoefficient c2d(mesh, csol, ny);

fem/gridfunc.hpp

@@ -532,9 +532,6 @@ public:
    std::unique_ptr<GridFunction> ProlongateToMaxOrder() const;
 
 protected:
-   void ProjectBdrCoefficientNormal(Coefficient *coeff, VectorCoefficient *vcoeff,
-                                    const Array<int> &attr);
-
    /** @brief Accumulates (depending on @a type) the values of @a coeff at all
        shared vdofs and counts in how many zones each vdof appears. */
    void AccumulateAndCountZones(Coefficient &coeff, AvgType type,
@@ -659,26 +656,15 @@ public:
    virtual void ProjectBdrCoefficient(Coefficient *coeff[],
                                       const Array<int> &attr);
 
-   /** @brief Project the normal component of the given VectorCoefficient on
-       the boundary. */
-   /** Only boundary attributes that are marked in @a bdr_attr are
-       projected. Assumes RT-type vector finite element GridFunction. */
+   /** Project the normal component of the given VectorCoefficient on
+       the boundary. Only boundary attributes that are marked in
+       'bdr_attr' are projected. Assumes RT-type VectorFE GridFunction. */
    void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff,
-                                    const Array<int> &bdr_attr)
-   { ProjectBdrCoefficientNormal(NULL, &vcoeff, bdr_attr); }
-
-   /** @brief Project the given Coefficient in the normal direction on the
-       boundary. */
-   /** Only boundary attributes that are marked in @a bdr_attr are projected.
-       Assumes RT-type vector finite element GridFunction. */
-   void ProjectBdrCoefficientNormal(Coefficient &coeff,
-                                    const Array<int> &bdr_attr)
-   { ProjectBdrCoefficientNormal(&coeff, NULL, bdr_attr); }
+                                    const Array<int> &bdr_attr);
 
    /** @brief Project the tangential components of the given VectorCoefficient
-       on the boundary. */
-   /** Only boundary attributes that are marked in @a bdr_attr
-       are projected. Assumes ND-type vector finite element GridFunction. */
+       on the boundary. Only boundary attributes that are marked in @a bdr_attr
+       are projected. Assumes ND-type VectorFE GridFunction. */
    virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
                                              const Array<int> &bdr_attr);
 
@@ -1928,7 +1914,7 @@ real_t ComputeElementLpDistance(real_t p, int i,
                                 GridFunction& gf1, GridFunction& gf2);
 
 
-/// Class used for extruding a scalar coefficient
+/// Class used for extruding scalar GridFunctions
 class ExtrudeCoefficient : public Coefficient
 {
 private:
@@ -1936,52 +1922,13 @@ private:
    Mesh *mesh_in;
    Coefficient &sol_in;
 public:
-   /// Constructs an instance of VectorExtrudeCoefficient
-   /**
-    * @param m      1D mesh
-    * @param s      1D vector coefficient
-    * @param n_     number of transverse elements of the extruded mesh
-    */
    ExtrudeCoefficient(Mesh *m, Coefficient &s, int n_)
-      : n(n_), mesh_in(m), sol_in(s)
-   { MFEM_VERIFY(n > 0, "Number of transverse elements must be positive!"); }
-
+      : n(n_), mesh_in(m), sol_in(s) { }
    real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
-
    virtual ~ExtrudeCoefficient() { }
 };
 
-/// Class used for extruding a vector coefficient
-class VectorExtrudeCoefficient : public VectorCoefficient
-{
-private:
-   int n;
-   Mesh *mesh_in;
-   VectorCoefficient &sol_in;
-public:
-   /// Constructs an instance of VectorExtrudeCoefficient
-   /**
-    * @param m      1D mesh
-    * @param s      1D vector coefficient
-    * @param n_     number of transverse elements of the extruded mesh
-    */
-   VectorExtrudeCoefficient(Mesh *m, VectorCoefficient &s, int n_)
-      : VectorCoefficient(s.GetVDim()), n(n_), mesh_in(m), sol_in(s)
-   { MFEM_VERIFY(n > 0, "Number of transverse elements must be positive!"); }
-
-   void Eval(Vector &v, ElementTransformation &T,
-             const IntegrationPoint &ip) override;
-
-   virtual ~VectorExtrudeCoefficient() { }
-};
-
-/// Extrude a 1D GridFunction, after extruding the mesh with Extrude1D()
-/**
- * @param mesh      1D mesh
- * @param mesh2d    extruded mesh
- * @param sol       grid function
- * @param ny        number of transverse elements of the extruded mesh
- */
+/// Extrude a scalar 1D GridFunction, after extruding the mesh with Extrude1D.
 GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
                                     GridFunction *sol, const int ny);
 

fem/integ/bilininteg_hdiv_ea.cpp

@@ -197,21 +197,15 @@ static void EAHdivAssemble3D(const int NE,
       // Assemble (one row per thread)
       MFEM_FOREACH_THREAD(idx_i, x, NDOF)
       {
-         // NOTE: due to an llvm backend bug, usage of the modulus operator
-         // has been removed from this foreach section.
          const int ic = idx_i / NDOF_C;
-         const int idx_ii = idx_i - ic * NDOF_C; // idx_i % NDOF_C
+         const int idx_ii = idx_i % NDOF_C;
 
          const int nx_i = (ic == 0) ? D1D : D1D-1;
          const int ny_i = (ic == 1) ? D1D : D1D-1;
 
-         const int qx_i = idx_ii / nx_i;
-         const int ix = idx_ii - qx_i * nx_i; // idx_ii % nx_i
-
-         const int qy_i = qx_i / ny_i;
-         const int iy = qx_i - qy_i * ny_i; // (idx_ii / nx_i) % ny_i
-
-         const int iz = qy_i; // (idx_ii / nx_i) / ny_i
+         const int ix = idx_ii % nx_i;
+         const int iy = (idx_ii / nx_i) % ny_i;
+         const int iz = (idx_ii / nx_i) / ny_i;
 
          const real_t (&Bi1)[MQ1][MD1] = (ic == 0) ? r_Bc : r_Bo;
          const real_t (&Bi2)[MQ1][MD1] = (ic == 1) ? r_Bc : r_Bo;
@@ -220,18 +214,14 @@ static void EAHdivAssemble3D(const int NE,
          for (int idx_j = 0; idx_j < NDOF; ++idx_j)
          {
             const int jc = idx_j / NDOF_C;
-            const int idx_jj = idx_j - jc * NDOF_C; // idx_j % NDOF_C
+            const int idx_jj = idx_j % NDOF_C;
 
             const int nx_j = (jc == 0) ? D1D : D1D-1;
             const int ny_j = (jc == 1) ? D1D : D1D-1;
 
-            const int qx_j = idx_jj / nx_j;
-            const int jx = idx_jj - qx_j * nx_j; // idx_jj % nx_j
-
-            const int qy_j = qx_j / ny_j;
-            const int jy = qx_j - qy_j * ny_j; // (idx_jj / nx_j) % ny_j
-
-            const int jz = qy_j; // (idx_jj / nx_j) / ny_j
+            const int jx = idx_jj % nx_j;
+            const int jy = (idx_jj / nx_j) % ny_j;
+            const int jz = (idx_jj / nx_j) / ny_j;
 
             const real_t (&Bj1)[MQ1][MD1] = (jc == 0) ? r_Bc : r_Bo;
             const real_t (&Bj2)[MQ1][MD1] = (jc == 1) ? r_Bc : r_Bo;

fem/intrules.hpp

@@ -125,6 +125,18 @@ private:
    void AddTriPoints3b(const int off, const real_t b, const real_t weight)
    { AddTriPoints3(off, (1. - b)/2., b, weight); }
 
+   void AddTriPoints3R(const int off, const real_t a, const real_t b,
+                       const real_t c, const real_t weight)
+   {
+      IntPoint(off + 0).Set2w(a, b, weight);
+      IntPoint(off + 1).Set2w(c, a, weight);
+      IntPoint(off + 2).Set2w(b, c, weight);
+   }
+
+   void AddTriPoints3R(const int off, const real_t a, const real_t b,
+                       const real_t weight)
+   { AddTriPoints3R(off, a, b, 1. - a - b, weight); }
+
    void AddTriPoints6(const int off, const real_t a, const real_t b,
                       const real_t c, const real_t weight)
    {
@@ -171,6 +183,14 @@ private:
       AddTetPoints3(off + 1, a, 1. - 3.*a, weight);
    }
 
+   // given b, add the permutations of (a,a,a,b), where 3*a + b = 1
+   void AddTetPoints4b(const int off, const real_t b, const real_t weight)
+   {
+      const real_t a = (1. - b)/3.;
+      IntPoint(off).Set(a, a, a, weight);
+      AddTetPoints3(off + 1, a, b, weight);
+   }
+
    // add the permutations of (a,a,b,b), 2*(a + b) = 1
    void AddTetPoints6(const int off, const real_t a, const real_t weight)
    {
@@ -189,37 +209,14 @@ private:
       AddTetPoints6(off + 6, a, bc, cb, weight);
    }
 
-   // add all 24 permutations of (a,b,c,d) where a+b+c+d = 1, all distinct
-   void AddTetPoints24(const int off, const real_t a, const real_t b,
-                       const real_t c, const real_t weight)
+   // given (b,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
+   void AddTetPoints12bc(const int off, const real_t b, const real_t c,
+                         const real_t weight)
    {
-      const real_t d = 1. - a - b - c;
-      // all 24 permutations of 4 distinct barycentric coordinates
-      // permuting which coordinate goes to x, y, z (4th is 1-x-y-z)
-      IntPoint(off +  0).Set(a, b, c, weight);
-      IntPoint(off +  1).Set(a, b, d, weight);
-      IntPoint(off +  2).Set(a, c, b, weight);
-      IntPoint(off +  3).Set(a, c, d, weight);
-      IntPoint(off +  4).Set(a, d, b, weight);
-      IntPoint(off +  5).Set(a, d, c, weight);
-      IntPoint(off +  6).Set(b, a, c, weight);
-      IntPoint(off +  7).Set(b, a, d, weight);
-      IntPoint(off +  8).Set(b, c, a, weight);
-      IntPoint(off +  9).Set(b, c, d, weight);
-      IntPoint(off + 10).Set(b, d, a, weight);
-      IntPoint(off + 11).Set(b, d, c, weight);
-      IntPoint(off + 12).Set(c, a, b, weight);
-      IntPoint(off + 13).Set(c, a, d, weight);
-      IntPoint(off + 14).Set(c, b, a, weight);
-      IntPoint(off + 15).Set(c, b, d, weight);
-      IntPoint(off + 16).Set(c, d, a, weight);
-      IntPoint(off + 17).Set(c, d, b, weight);
-      IntPoint(off + 18).Set(d, a, b, weight);
-      IntPoint(off + 19).Set(d, a, c, weight);
-      IntPoint(off + 20).Set(d, b, a, weight);
-      IntPoint(off + 21).Set(d, b, c, weight);
-      IntPoint(off + 22).Set(d, c, a, weight);
-      IntPoint(off + 23).Set(d, c, b, weight);
+      const real_t a = (1. - b - c)/2.;
+      AddTetPoints3(off,     a, b, weight);
+      AddTetPoints3(off + 3, a, c, weight);
+      AddTetPoints6(off + 6, a, b, c, weight);
    }
 
 public:

fem/linearform.cpp

@@ -297,8 +297,7 @@ void LinearForm::Assemble()
          tr = mesh->GetBdrFaceTransformations(i);
          if (tr != NULL)
          {
-            mfem::DofTransformation doftrans;
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
             for (int k = 0; k < boundary_face_integs.Size(); k++)
             {
                if (boundary_face_integs_marker[k] &&
@@ -308,7 +307,6 @@ void LinearForm::Assemble()
                boundary_face_integs[k]->
                AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
                                       *tr, elemvect);
-               doftrans.TransformDual(elemvect);
                AddElementVector (vdofs, elemvect);
             }
          }

fem/prestriction.cpp

@@ -321,17 +321,12 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
    const int vd = vdim;
    const bool t = byvdim;
    const int threshold = ndofs;
-   const int nsdofs = pfes.GetFaceNbrVSize() / vd;
+   const int nsdofs = pfes.GetFaceNbrVSize();
    auto d_indices1 = scatter_indices1.Read();
    auto d_indices2 = scatter_indices2.Read();
    auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
-   const int ne_shared = nsdofs / elem_dofs;
-   const int nedof = elem_dofs;
-   // Note: the shape of face_nbr_data, as determined by
-   // ParFiniteElementSpace::ExchangeFaceNbrData, is (elem_dofs, vdim,
-   // ne_shared), independent of the ordering (byNODES or byVDIM) of the finite
-   // element space.
-   auto d_x_shared = Reshape(face_nbr_data.Read(), elem_dofs, vd, ne_shared);
+   auto d_x_shared = Reshape(face_nbr_data.Read(),
+                             t?vd:nsdofs, t?nsdofs:vd);
    auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
    mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
    {
@@ -351,9 +346,8 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
          }
          else if (idx2>=threshold) // shared boundary
          {
-            const int e_shared = (idx2 - threshold) / nedof;
-            const int i_shared = (idx2 - threshold) % nedof;
-            d_y(dof, c, 1, face) = d_x_shared(i_shared,c,e_shared);
+            d_y(dof, c, 1, face) = d_x_shared(t?c:(idx2-threshold),
+                                              t?(idx2-threshold):c);
          }
          else // true boundary
          {

fem/restriction.cpp

@@ -1398,17 +1398,20 @@ void L2FaceRestriction::PermuteAndSetSharedFaceDofsScatterIndices2(
    const int dim = fes.GetMesh()->Dimension();
    const int dof1d = fes.GetTypicalFE()->GetOrder()+1;
    fes.GetTypicalFE()->GetFaceMap(face_id2, face_map);
+   Array<int> face_nbr_dofs;
+   const ParFiniteElementSpace &pfes =
+      static_cast<const ParFiniteElementSpace&>(this->fes);
+   pfes.GetFaceNbrElementVDofs(elem_index, face_nbr_dofs);
 
    for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
    {
       const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
                                                orientation, dof1d, face_dof_elem1);
       const int volume_dof_elem2 = face_map[face_dof_elem2];
-      // Encode the volume DOF index and element index
-      const int global_dof_elem2 = elem_index*elem_dofs + volume_dof_elem2;
+      const int global_dof_elem2 = face_nbr_dofs[volume_dof_elem2];
       const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
       // Trick to differentiate dof location inter/shared
-      scatter_indices2[restriction_dof_elem2] = ndofs + global_dof_elem2;
+      scatter_indices2[restriction_dof_elem2] = ndofs+global_dof_elem2;
    }
 #endif
 }

general/arrays_by_name.hpp

@@ -278,18 +278,9 @@ void ArraysByName<T>::Load(std::istream &in)
          q1 = ArrayLine.find(' ');
          ArrayName = ArrayLine.substr(0,q1-1);
       }
-      if (q1+2 < ArrayLine.size())
-      {
-         // Read the remainder of the line which contains the array data
-         std::istringstream ArrayDataStream(ArrayLine.substr(q1+2,
-                                                             ArrayLine.size()));
-         data[ArrayName].Load(ArrayDataStream, 0);
-      }
-      else
-      {
-         // Read the array data starting on the next line
-         data[ArrayName].Load(in, 0);
-      }
+
+      // Ignore the remainder of the line which may contain explanatory comments
+      data[ArrayName].Load(in, 0);
    }
 
 }

general/device.cpp

@@ -726,16 +726,16 @@ std::string Device::GetUUID(const int device_id)
    MFEM_GPU_CHECK(cudaGetDeviceProperties(&prop, device_id));
    for (int i = 0; i < 16; ++i)
    {
-      const unsigned b = static_cast<unsigned char>(prop.uuid.bytes[i]);
-      res << std::setfill('0') << std::setw(2) << std::hex << b;
+      res << std::setfill('0') << std::setw(2) << std::hex
+          << static_cast<unsigned>(prop.uuid.bytes[i]);
    }
 #elif defined(MFEM_USE_HIP)
    hipUUID uuid;
    MFEM_GPU_CHECK(hipDeviceGetUuid(&uuid, device_id));
    for (int i = 0; i < 16; ++i)
    {
-      const unsigned b = static_cast<unsigned char>(uuid.bytes[i]);
-      res << std::setfill('0') << std::setw(2) << std::hex << b;
+      res << std::setfill('0') << std::setw(2) << std::hex
+          << static_cast<unsigned>(uuid.bytes[i]);
    }
 #endif
    return res.str();

linalg/hypre.cpp

@@ -317,9 +317,6 @@ void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)
 
 Vector * HypreParVector::GlobalVector() const
 {
-   MFEM_VERIFY(size > 0,
-               "GlobalVector method can only be called on vectors wherein each "
-               "process owns one or more entries");
    hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
    Vector *v = new Vector(hv->data, internal::to_int(hv->size));
    v->MakeDataOwner();

mesh/mesh.hpp

@@ -3206,22 +3206,10 @@ public:
 
 
 /// Extrude a 1D mesh
-/**
- * @param mesh      1D mesh
- * @param ny        number of transverse elements of the extruded mesh
- * @param sy        physical size in the direction of extrusion
- * @param closed    if false, only the original boundaries are extruded,
- *                  otherwise boundaries are generated all around the domain
- */
 Mesh *Extrude1D(Mesh *mesh, const int ny, const real_t sy,
                 const bool closed = false);
 
 /// Extrude a 2D mesh
-/**
- * @param mesh      2D mesh
- * @param nz        number of transverse elements of the extruded mesh
- * @param sz        physical size in the direction of extrusion
- */
 Mesh *Extrude2D(Mesh *mesh, const int nz, const real_t sz);
 
 /** @brief Constructs the smallest possible [0,1]^dim serial mesh that can be

mesh/mesh_readers.cpp

@@ -1516,15 +1516,12 @@ void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
 void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
 {
    string buff;
-   string version;
+   real_t version;
    int binary, dsize;
    input >> version >> binary >> dsize;
-   if (version != "2.2")
+   if (version < 2.2)
    {
-      MFEM_ABORT("Gmsh file version must be 2.2, found version "
-                 << version << ".\n"
-                 "To convert your mesh to the required format, use:\n"
-                 "  gmsh -format msh22 -save -o output.msh input.msh");
+      MFEM_ABORT("Gmsh file version < 2.2");
    }
    if (dsize != sizeof(double))
    {

tests/unit/fem/test_face_restriction.cpp

@@ -117,52 +117,3 @@ TEST_CASE("Vector FE Face Restriction", "[FaceRestriction]")
    gf2 -= gf;
    REQUIRE(gf2.Normlinf() == MFEM_Approx(0.0));
 }
-
-#ifdef MFEM_USE_MPI
-
-TEST_CASE("L2 Face Restriction", "[FaceRestriction][Parallel]")
-{
-   const int dim = GENERATE(2, 3);
-   constexpr int nx = 3;
-   constexpr int order = 2;
-   constexpr int vdim = 2;
-   const Ordering::Type ordering = GENERATE(Ordering::byNODES, Ordering::byVDIM);
-
-   Mesh serial_mesh = MakeCartesianMesh(nx, dim);
-   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
-
-   L2_FECollection fec(order, dim, BasisType::GaussLobatto);
-   ParFiniteElementSpace fes(&mesh, &fec, vdim, ordering);
-
-   auto *R = fes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,
-                                    FaceType::Interior);
-
-   Vector vals({1.0, 2.0});
-   VectorConstantCoefficient coeff(vals);
-
-   ParGridFunction gf(&fes);
-   gf.ProjectCoefficient(coeff);
-
-   Vector face_vec(R->Height());
-   R->Mult(gf, face_vec);
-
-   const int nf = mesh.GetNFbyType(FaceType::Interior);
-   const int face_dofs = fes.GetTypicalTraceElement()->GetDof();
-   auto h_face_vec = Reshape(face_vec.HostRead(), face_dofs, vdim, 2, nf);
-
-   for (int f = 0; f < nf; ++f)
-   {
-      for (int m = 0; m < 2; ++m)
-      {
-         for (int c = 0; c < vdim; ++c)
-         {
-            for (int i = 0; i < face_dofs; ++i)
-            {
-               REQUIRE(h_face_vec(i, c, m, f) == vals[c]);
-            }
-         }
-      }
-   }
-}
-
-#endif

tests/unit/fem/test_intrules.cpp

@@ -105,39 +105,39 @@ TEST_CASE("Integration rule order initialization", "[IntegrationRules]")
    SECTION("Segment rule constructed by accessing square rule")
    {
       auto &quad5_ir = intrules.Get(Geometry::SQUARE, 5);
-      REQUIRE(quad5_ir.GetOrder() >= 5);
+      REQUIRE(quad5_ir.GetOrder() == 5);
       // The segment integration rule of order 5 is lazy constructed when we get
       // the square integration rule of order 5. Make sure its order was
       // properly set:
       auto &line5_ir = intrules.Get(Geometry::SEGMENT, 5);
-      REQUIRE(line5_ir.GetOrder() >= 5);
+      REQUIRE(line5_ir.GetOrder() == 5);
    }
 
    SECTION("Segment rule constructed by accessing cube rule")
    {
       auto &hex7_ir = intrules.Get(Geometry::CUBE, 7);
-      REQUIRE(hex7_ir.GetOrder() >= 7);
+      REQUIRE(hex7_ir.GetOrder() == 7);
       // The segment integration rule of order 7 is lazy constructed when we get
       // the cube integration rule of order 7. Make sure its order was properly
       // set:
       auto &line7_ir = intrules.Get(Geometry::SEGMENT, 7);
-      REQUIRE(line7_ir.GetOrder() >= 7);
+      REQUIRE(line7_ir.GetOrder() == 7);
    }
 
    SECTION("Segment and triangle rules constructed by accessing prism rule")
    {
       auto &prism3_ir = intrules.Get(Geometry::PRISM, 3);
-      REQUIRE(prism3_ir.GetOrder() >= 3);
+      REQUIRE(prism3_ir.GetOrder() == 3);
       // The segment integration rule of order 3 is lazy constructed when we get
       // the prism integration rule of order 3. Make sure its order was properly
       // set:
       auto &line3_ir = intrules.Get(Geometry::SEGMENT, 3);
-      REQUIRE(line3_ir.GetOrder() >= 3);
+      REQUIRE(line3_ir.GetOrder() == 3);
       // The triangle integration rule of order 3 is lazy constructed when we
       // get the prism integration rule of order 3. Make sure its order was
       // properly set:
       auto &tri3_ir = intrules.Get(Geometry::TRIANGLE, 3);
-      REQUIRE(tri3_ir.GetOrder() >= 3);
+      REQUIRE(tri3_ir.GetOrder() == 3);
    }
 }
 
@@ -271,43 +271,3 @@ TEST_CASE("Simplex integration rules", "[SimplexRules]")
       }
    }
 }
-
-
-// Monomial exactness is tested by [SimplexRules] above, which now uses
-// positive-weight rules by default. The tests below verify properties
-// specific to the positive-weight rules: weight positivity, stability,
-// and interior point placement.
-
-TEST_CASE("Simplex rule positivity", "[IntegrationRules]")
-{
-   IntegrationRules rules;
-
-   SECTION("triangle rules have all positive weights for orders 0-25")
-   {
-      for (int order = 0; order <= 25; order++)
-      {
-         const IntegrationRule &ir = rules.Get(Geometry::TRIANGLE, order);
-         for (int i = 0; i < ir.GetNPoints(); i++)
-         {
-            INFO("order=" << order << ", point=" << i);
-            REQUIRE(ir.IntPoint(i).weight > 0.0);
-         }
-      }
-   }
-
-   SECTION("tet rules have all positive weights for orders 0-20")
-   {
-      for (int order = 0; order <= 20; order++)
-      {
-         const IntegrationRule &ir =
-            rules.Get(Geometry::TETRAHEDRON, order);
-         for (int i = 0; i < ir.GetNPoints(); i++)
-         {
-            INFO("order=" << order << ", point=" << i);
-            REQUIRE(ir.IntPoint(i).weight > 0.0);
-         }
-      }
-   }
-}
-
-

tests/unit/fem/test_project_bdr.cpp

@@ -25,201 +25,15 @@ void Func_3D_lin(const Vector &x, Vector &v)
    v[2] = -2.572 * x[0] + 1.321 * x[1] + 3.234 * x[2];
 }
 
-TEST_CASE("3D ProjectBdrCoefficientNormal Vector",
-          "[GridFunction]"
-          "[VectorGridFunctionCoefficient]")
-{
-   const int n = 1;
-   const int dim = 3;
-   const int order = 1;
-
-   const double tol = 1e-6;
-
-   for (int type = (int)Element::TETRAHEDRON;
-        type <= (int)Element::HEXAHEDRON; type++)
-   {
-      Mesh mesh = Mesh::MakeCartesian3D(
-                     n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
-
-      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
-
-      SECTION("3D GetVectorValue tests for element type " +
-              std::to_string(type))
-      {
-         RT_FECollection  rt_fec(order+1, dim);
-
-         FiniteElementSpace  rt_fespace(&mesh,  &rt_fec);
-
-         GridFunction  rt_x( &rt_fespace);
-
-         VectorGridFunctionCoefficient  rt_xCoef( &rt_x);
-
-         Array<int> bdr_marker(6);
-
-         Vector normal(dim);
-         Vector  f_val(dim);
-         Vector rt_val(dim);
-
-         for (int b = 1; b<=6; b++)
-         {
-            bdr_marker = 0;
-            bdr_marker[b-1] = 1;
-
-            rt_x = 0.0;
-            rt_x.ProjectBdrCoefficientNormal(funcCoef, bdr_marker);
-
-            for (int be = 0; be < mesh.GetNBE(); be++)
-            {
-               Element *e = mesh.GetBdrElement(be);
-               if (e->GetAttribute() != b) { continue; }
-
-               ElementTransformation *T = mesh.GetBdrElementTransformation(be);
-               const FiniteElement   *fe = rt_fespace.GetBE(be);
-               const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
-                                                        2*order + 2);
-
-               double rt_err = 0.0;
-
-               for (int j=0; j<ir.GetNPoints(); j++)
-               {
-                  const IntegrationPoint &ip = ir.IntPoint(j);
-                  T->SetIntPoint(&ip);
-
-                  CalcOrtho(T->Jacobian(), normal);
-
-                  funcCoef.Eval(f_val, *T, ip);
-                  rt_xCoef.Eval(rt_val, *T, ip);
-
-                  rt_val -= f_val;
-
-                  double rt_dist = rt_val * normal;
-
-                  rt_err += rt_dist;
-
-                  if (verbose_tests && rt_dist > tol)
-                  {
-                     mfem::out << be << ":" << j << " rt ("
-                               << f_val[0] << "," << f_val[1] << "," << f_val[2]
-                               << ") vs. ("
-                               << rt_val[0] << "," << rt_val[1] << ","
-                               << rt_val[2] << ") " << rt_dist << std::endl;
-                  }
-               }
-               rt_err  /= ir.GetNPoints();
-
-               REQUIRE( rt_err == MFEM_Approx(0.0));
-            }
-         }
-      }
-   }
-}
-
-TEST_CASE("3D ProjectBdrCoefficientNormal Scalar",
-          "[GridFunction]"
-          "[VectorGridFunctionCoefficient]")
-{
-   const int n = 1;
-   const int dim = 3;
-   const int order = 1;
-
-   const double tol = 1e-6;
-
-   const char bdrs_axis[] = {2, 1, 0, 1, 0, 2};
-   const char bdrs_sign[] = {-1, -1, +1, +1, -1, +1};
-
-   for (int type = (int)Element::TETRAHEDRON;
-        type <= (int)Element::HEXAHEDRON; type++)
-   {
-      Mesh mesh = Mesh::MakeCartesian3D(
-                     n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
-
-      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
-
-      SECTION("3D GetVectorValue tests for element type " +
-              std::to_string(type))
-      {
-         RT_FECollection  rt_fec(order+1, dim);
-
-         FiniteElementSpace  rt_fespace(&mesh,  &rt_fec);
-
-         GridFunction  rt_x( &rt_fespace);
-
-         VectorGridFunctionCoefficient  rt_xCoef( &rt_x);
-
-         Array<int> bdr_marker(6);
-
-         Vector normal(dim);
-         Vector  f_val(dim);
-         Vector rt_val(dim);
-
-         for (int b = 1; b<=6; b++)
-         {
-            bdr_marker = 0;
-            bdr_marker[b-1] = 1;
-
-            rt_x = 0.0;
-
-            normal = 0.;
-            normal(bdrs_axis[b-1]) = (bdrs_sign[b-1] > 0)?(+1.):(-1.);
-            VectorConstantCoefficient normCoef(normal);
-            InnerProductCoefficient prodCoef(funcCoef, normCoef);
-            rt_x.ProjectBdrCoefficientNormal(prodCoef, bdr_marker);
-
-            for (int be = 0; be < mesh.GetNBE(); be++)
-            {
-               Element *e = mesh.GetBdrElement(be);
-               if (e->GetAttribute() != b) { continue; }
-
-               ElementTransformation *T = mesh.GetBdrElementTransformation(be);
-               const FiniteElement   *fe = rt_fespace.GetBE(be);
-               const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
-                                                        2*order + 2);
-
-               double rt_err = 0.0;
-
-               for (int j=0; j<ir.GetNPoints(); j++)
-               {
-                  const IntegrationPoint &ip = ir.IntPoint(j);
-                  T->SetIntPoint(&ip);
-
-                  CalcOrtho(T->Jacobian(), normal);
-
-                  funcCoef.Eval(f_val, *T, ip);
-                  rt_xCoef.Eval(rt_val, *T, ip);
-
-                  rt_val -= f_val;
-
-                  double rt_dist = rt_val * normal;
-
-                  rt_err += rt_dist;
-
-                  if (verbose_tests && rt_dist > tol)
-                  {
-                     mfem::out << be << ":" << j << " rt ("
-                               << f_val[0] << "," << f_val[1] << "," << f_val[2]
-                               << ") vs. ("
-                               << rt_val[0] << "," << rt_val[1] << ","
-                               << rt_val[2] << ") " << rt_dist << std::endl;
-                  }
-               }
-               rt_err  /= ir.GetNPoints();
-
-               REQUIRE( rt_err == MFEM_Approx(0.0));
-            }
-         }
-      }
-   }
-}
-
 TEST_CASE("3D ProjectBdrCoefficientTangent",
           "[GridFunction]"
           "[VectorGridFunctionCoefficient]")
 {
-   const int n = 1;
-   const int dim = 3;
-   const int order = 1;
+   int n = 1;
+   int dim = 3;
+   int order = 1;
 
-   const double tol = 1e-6;
+   double tol = 1e-6;
 
    for (int type = (int)Element::TETRAHEDRON;
         type <= (int)Element::HEXAHEDRON; type++)

tests/unit/general/test_arrays_by_name.cpp

@@ -200,39 +200,3 @@ TEST_CASE("ArraysByName Sort/Unique Methods", "[ArraysByName]")
       }
    }
 }
-
-TEST_CASE("ArraysByName Print/Load Methods", "[ArraysByName]")
-{
-   ArraysByName<int> abn;
-
-   FillArraysByName(abn);
-
-   // Print object to string using default format
-   std::ostringstream oss1;
-   abn.Print(oss1);
-
-   // Load new object from printed output
-   ArraysByName<int> abn_load1;
-   std::istringstream iss1(oss1.str());
-   abn_load1.Load(iss1);
-   REQUIRE(abn == abn_load1);
-
-   // Print object to string using one line per array
-   std::ostringstream oss2;
-   oss2 << abn.Size() << '\n';
-   for (auto a : abn)
-   {
-      oss2 << '"' << a.first << "\" " << a.second.Size();
-      for (auto d : a.second)
-      {
-         oss2 << ' ' << d;
-      }
-      oss2 << '\n';
-   }
-
-   // Load new object from printed output
-   ArraysByName<int> abn_load2;
-   std::istringstream iss2(oss2.str());
-   abn_load2.Load(iss2);
-   REQUIRE(abn == abn_load2);
-}