Try the 'lgtm' CodeQL queries

Windows CI METIS workaround

.github/workflows/build-container.yaml

@@ -1,4 +1,4 @@
-name: Build Deploy Container
+name: "Docker"
 
 on:
 

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

@@ -10,7 +10,7 @@
 # CONTRIBUTING.md for details.
 
 # In this CI section, we build different variants of mfem and run test on them.
-name: builds-and-tests
+name: "Tests"
 
 # Github actions can use the default "GITHUB_TOKEN". By default, this token
 # is set to have permissive access. However, this is not a good practice
@@ -47,17 +47,17 @@ jobs:
   builds-and-tests:
     strategy:
       matrix:
-        os: [ubuntu-20.04, macos-10.15, windows-2022]
+        os: [ubuntu-latest, macos-latest, windows-latest]
         target: [dbg, opt]
         mpi: [seq, par]
         build-system: [make, cmake]
         hypre-target: [int32]
         exclude:
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             build-system: cmake
-          - os: macos-10.15
+          - os: macos-latest
             build-system: cmake
-          - os: windows-2022
+          - os: windows-latest
             build-system: make
         # 'include' allows us to:
         # - Add a variable to all jobs without creating a new matrix dimension.
@@ -72,15 +72,15 @@ jobs:
             codecov: NO
           - target: opt
             codecov: YES
-          - os: windows-2022
+          - os: windows-latest
             codecov: NO
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             target: opt
             codecov: NO
             mpi: par
             build-system: cmake
             hypre-target: int32
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             target: opt
             codecov: NO
             mpi: par
@@ -112,35 +112,35 @@ jobs:
     # TODO: It would be nice to have only one step, e.g. with a dedicated
     # action, but I (@adrienbernede) don't see how at the moment.
     - name: get MPI (Linux)
-      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-20.04'
+      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-latest'
       run: |
         sudo apt-get install mpich libmpich-dev
         export MAKE_CXX_FLAG="MPICXX=mpic++"
 
     - name: get lcov (Linux)
-      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-20.04'
+      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-latest'
       run: |
         sudo apt-get install lcov
 
     - name: Set up Homebrew
-      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-10.15'
+      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
       uses: Homebrew/actions/setup-homebrew@c4aafe8c4620bf08883dd4679c374f11e73329d3
 
     - name: get MPI (MacOS)
-      if: matrix.mpi == 'par' && matrix.os == 'macos-10.15'
+      if: matrix.mpi == 'par' && matrix.os == 'macos-latest'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
         brew install openmpi
         export MAKE_CXX_FLAG="MPICXX=mpic++"
 
     - name: get MPI (MacOS)
-      if: matrix.codecov == 'YES' && matrix.os == 'macos-10.15'
+      if: matrix.codecov == 'YES' && matrix.os == 'macos-latest'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
         brew install lcov
 
     - name: get MPI (Windows)
-      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
       uses: mpi4py/setup-mpi@v1.0.3
 
     # Get Hypre through cache, or build it.
@@ -154,7 +154,7 @@ jobs:
         key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.2
 
     - name: get hypre
-      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-2022'
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-latest'
       uses: mfem/github-actions/build-hypre@v2.2
       with:
         archive: ${{ env.HYPRE_ARCHIVE }}
@@ -163,7 +163,7 @@ jobs:
         build-system: make
 
     - name: get hypre (Windows)
-      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-latest'
       uses: mfem/github-actions/build-hypre@v2.2
       with:
         archive: ${{ env.HYPRE_ARCHIVE }}
@@ -175,14 +175,14 @@ jobs:
     # Install will only run on cache miss.
     - name: cache metis
       id: metis-cache
-      if: matrix.mpi == 'par' && matrix.os != 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
       uses: actions/cache@v2
       with:
         path: ${{ env.METIS_TOP_DIR }}
         key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
 
     - name: install metis
-      if: matrix.mpi == 'par' && matrix.os != 'windows-2022' && steps.metis-cache.outputs.cache-hit != 'true'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-latest' && steps.metis-cache.outputs.cache-hit != 'true'
       uses: mfem/github-actions/build-metis@v2.2
       with:
         archive: ${{ env.METIS_ARCHIVE }}
@@ -196,20 +196,16 @@ jobs:
         key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
 
     - name: prepare binary cache location
-      if: matrix.os == 'windows-2022' && steps.vcpkg-cache.outputs.cache-hit != 'true'
+      if: matrix.os == 'windows-latest' && steps.vcpkg-cache.outputs.cache-hit != 'true'
       run: |
         mkdir -p vcpkg_cache
 
     - name: install metis (Windows)
-      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
       env:
         VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
       run: |
-        $PortFile = 'C:\vcpkg\ports\metis\portfile.cmake'
-        $OriginalURL = 'http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-${METIS_VERSION}.tar.gz'
-        $NewURL = 'https://github.com/mfem/tpls/raw/gh-pages/metis-5.1.0.tar.gz'
-        (Get-Content $PortFile).replace($OriginalURL, $NewURL) | Set-Content $PortFile
-        vcpkg install metis --triplet=x64-windows-static
+        vcpkg install metis-mfem --triplet=x64-windows-static --overlay-ports=${{ env.MFEM_TOP_DIR }}/config/vcpkg/ports
 
     # MFEM build and test
     - name: build
@@ -252,7 +248,7 @@ jobs:
       shell: bash
 
     - name: cmake unit tests (Ubuntu 20.04)
-      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-20.04'
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-latest'
       run: |
         CTEST_CONFIG="Release"
         [[ ${{ matrix.target }} == 'dbg' ]] && CTEST_CONFIG="Debug"
@@ -260,7 +256,7 @@ jobs:
       shell: bash
 
     - name: cmake tests
-      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-20.04'
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-latest'
       run: |
         CTEST_CONFIG="Release"
         cd ${{ env.MFEM_TOP_DIR }}/build && ctest --output-on-failure -C ${CTEST_CONFIG}

.github/workflows/codeql-analysis.yml

@@ -0,0 +1,71 @@
+# For most projects, this workflow file will not need changing; you simply need
+# to commit it to your repository.
+#
+# You may wish to alter this file to override the set of languages analyzed,
+# or to provide custom queries or build logic.
+#
+# ******** NOTE ********
+# We have attempted to detect the languages in your repository. Please check
+# the `language` matrix defined below to confirm you have the correct set of
+# supported CodeQL languages.
+#
+name: "Static Analysis"
+
+on:
+  push:
+    branches: [ "master", "next"]
+  pull_request:
+    # The branches below must be a subset of the branches above
+    branches: [ "master" ]
+
+jobs:
+  analyze:
+    name: Analyze
+    runs-on: ubuntu-latest
+    permissions:
+      actions: read
+      contents: read
+      security-events: write
+
+    strategy:
+      fail-fast: false
+      matrix:
+        language: [ 'cpp' ]
+        # CodeQL supports [ 'cpp', 'csharp', 'go', 'java', 'javascript', 'python', 'ruby' ]
+        # Learn more about CodeQL language support at https://aka.ms/codeql-docs/language-support
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v3
+
+    # Initializes the CodeQL tools for scanning.
+    - name: Initialize CodeQL
+      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.
+        # By default, queries listed here will override any specified in a config file.
+        # Prefix the list here with "+" to use these queries and those in the config file.
+
+        # Details on CodeQL's query packs refer to : https://docs.github.com/en/code-security/code-scanning/automatically-scanning-your-code-for-vulnerabilities-and-errors/configuring-code-scanning#using-queries-in-ql-packs
+        # queries: security-extended,security-and-quality
+        queries: lgtm
+
+
+    # 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@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
+
+    #   If the Autobuild fails above, remove it and uncomment the following three lines.
+    #   modify them (or add more) to build your code if your project, please refer to the EXAMPLE below for guidance.
+
+    # - run: |
+    #   echo "Run, Build Application using script"
+    #   ./location_of_script_within_repo/buildscript.sh
+
+    - name: Perform CodeQL Analysis
+      uses: github/codeql-action/analyze@v2

.github/workflows/mfem-analysis.yml

@@ -9,7 +9,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-name: build-analysis
+name: "Build Analysis"
 
 permissions:
   actions: write
@@ -31,7 +31,7 @@ env:
 
 jobs:
   gitignore:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
 
     steps:
     - name: Cancel Previous Runs

.github/workflows/repo-check.yml

@@ -9,7 +9,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-name: repo-check
+name: "Checks"
 
 permissions:
   actions: write
@@ -28,7 +28,7 @@ on:
 
 jobs:
   file-headers-check:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -78,7 +78,7 @@ jobs:
         exit 1
 
   code-style:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -88,14 +88,14 @@ jobs:
 
     - name: get astyle
       run: |
-        sudo apt-get install astyle=3.1-1ubuntu2
+        sudo apt-get install astyle
 
     - name: style check
       run: |
         ./config/githooks/pre-push --style
 
   documentation:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -106,6 +106,8 @@ jobs:
     - name: get doxygen and graphviz
       run: |
         sudo apt-get install doxygen graphviz
+        cd doc
+        doxygen -u CodeDocumentation.conf.in 2>/dev/null
 
     - name: build documentation
       run: |
@@ -118,7 +120,7 @@ jobs:
       github.ref != 'refs/heads/master' &&
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     steps:
     - name: checkout mfem
       uses: actions/checkout@v2

config/vcpkg/ports/metis-mfem/disable-programs.patch

@@ -0,0 +1,8 @@
+--- a/CMakeLists.txt	Wed Dec 21 18:24:22 2016
++++ b/CMakeLists.txt	Wed Dec 21 18:24:26 2016
+@@ -20,4 +20,4 @@
+ # Recursively look for CMakeLists.txt in subdirs.
+ add_subdirectory("include")
+ add_subdirectory("libmetis")
+-add_subdirectory("programs")
++# add_subdirectory("programs")

config/vcpkg/ports/metis-mfem/enable-install.patch

@@ -0,0 +1,15 @@
+--- a/CMakeLists.txt	Sat Mar 30 17:24:45 2013
++++ b/CMakeLists.txt	Wed Dec 21 18:23:43 2016
+@@ -4,11 +4,7 @@
+ set(GKLIB_PATH "GKlib" CACHE PATH "path to GKlib")
+ set(SHARED FALSE CACHE BOOL "build a shared library")
+ 
+-if(MSVC)
+-  set(METIS_INSTALL FALSE)
+-else()
+-  set(METIS_INSTALL TRUE)
+-endif()
++set(METIS_INSTALL TRUE)
+ 
+ # Configure libmetis library.
+ if(SHARED)

config/vcpkg/ports/metis-mfem/fix-INT_MIN_define.patch

@@ -0,0 +1,34 @@
+diff --git a/include/metis.h b/include/metis.h
+index dc5406a..7732437 100644
+--- a/include/metis.h
++++ b/include/metis.h
+@@ -72,10 +72,14 @@ typedef __int64 int64_t;
+ #define PRId64       "I64d"
+ #define SCNd32       "ld"
+ #define SCNd64       "I64d"
++#ifdef _WIN32
++#include <stdint.h>
++#else
+ #define INT32_MIN    ((int32_t)_I32_MIN)
+ #define INT32_MAX    _I32_MAX
+ #define INT64_MIN    ((int64_t)_I64_MIN)
+ #define INT64_MAX    _I64_MAX
++#endif
+ #else
+ #include <inttypes.h>
+ #endif
+diff --git a/GKlib/gk_arch.h b/GKlib/gk_arch.h
+index 78b1431..7258763 100644
+--- a/GKlib/gk_arch.h
++++ b/GKlib/gk_arch.h
+@@ -32,8 +32,8 @@
+ 
+ 
+ #ifdef __MSC__ 
+-  #include "ms_stdint.h"
+-  #include "ms_inttypes.h"
++  #include <stdint.h>
++  #include <inttypes.h>
+   #include "ms_stat.h"
+ #else
+ #ifndef SUNOS

config/vcpkg/ports/metis-mfem/fix-gklib-vs14-math.patch

@@ -0,0 +1,11 @@
+--- a/GKlib/gk_arch.h	Wed Dec 21 18:34:18 2016
++++ b/GKlib/gk_arch.h	Wed Dec 21 18:30:49 2016
+@@ -58,7 +58,7 @@
+ #define PTRDIFF_MAX  INT64_MAX
+ #endif
+ 
+-#ifdef __MSC__
++#if defined(__MSC__) && (_MSC_VER < 1900)
+ /* MSC does not have rint() function */
+ #define rint(x) ((int)((x)+0.5))  
+ 

config/vcpkg/ports/metis-mfem/fix-linux-build-error.patch

@@ -0,0 +1,14 @@
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+index e94f050..b9613a7 100644
+--- a/CMakeLists.txt
++++ b/CMakeLists.txt
+@@ -1,7 +1,8 @@
+ cmake_minimum_required(VERSION 2.8)
+ project(METIS)
+ 
+-set(GKLIB_PATH "GKlib" CACHE PATH "path to GKlib")
++set(GKLIB_PATH "${CMAKE_SOURCE_DIR}/GKlib" CACHE PATH "path to GKlib")
++
+ set(SHARED FALSE CACHE BOOL "build a shared library")
+ 
+ set(METIS_INSTALL TRUE)

config/vcpkg/ports/metis-mfem/fix-metis-vs14-math.patch

@@ -0,0 +1,11 @@
+--- a/libmetis/metislib.h	Sat Mar 30 17:24:45 2013
++++ b/libmetis/metislib.h	Wed Dec 21 18:30:59 2016
+@@ -31,7 +31,7 @@
+ #include <proto.h>
+ 
+ 
+-#if defined(COMPILER_MSC)
++#if defined(COMPILER_MSC) && (_MSC_VER < 1900)
+ #if defined(rint)
+   #undef rint
+ #endif

config/vcpkg/ports/metis-mfem/fix-runtime-install-destination.patch

@@ -0,0 +1,10 @@
+--- a/libmetis/CMakeLists.txt	Sat Mar 30 17:24:45 2013
++++ b/libmetis/CMakeLists.txt	Wed Dec 21 17:41:37 2016
+@@ -11,6 +11,6 @@
+ if(METIS_INSTALL)
+   install(TARGETS metis
+     LIBRARY DESTINATION lib
+-    RUNTIME DESTINATION lib
++    RUNTIME DESTINATION bin
+     ARCHIVE DESTINATION lib)
+ endif()

config/vcpkg/ports/metis-mfem/install-metisConfig.patch

@@ -0,0 +1,44 @@
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+index b9613a7..e43ffee 100644
+--- a/CMakeLists.txt
++++ b/CMakeLists.txt
+@@ -22,3 +22,23 @@ include_directories(include)
+ add_subdirectory("include")
+ add_subdirectory("libmetis")
+ # add_subdirectory("programs")
++
++if(METIS_INSTALL)
++  set(PRJ_NAME metis)
++  set(PRJ_VER 5.1.0)
++  install(EXPORT metisTargets 
++    FILE ${PRJ_NAME}Targets.cmake
++    DESTINATION lib/cmake/${PRJ_NAME})
++  include(CMakePackageConfigHelpers)
++  write_basic_package_version_file(
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}ConfigVersion.cmake
++    VERSION ${PRJ_VER}
++    COMPATIBILITY SameMajorVersion)
++  file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}Config.cmake
++    "include(\${CMAKE_CURRENT_LIST_DIR}/${PRJ_NAME}Targets.cmake)")
++  install(FILES 
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}ConfigVersion.cmake
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}Config.cmake
++    DESTINATION lib/cmake/${PRJ_NAME})
++endif()
++
+diff --git a/libmetis/CMakeLists.txt b/libmetis/CMakeLists.txt
+index 7a5fc74..5a68cf0 100644
+--- a/libmetis/CMakeLists.txt
++++ b/libmetis/CMakeLists.txt
+@@ -9,8 +9,9 @@ if(UNIX)
+ endif()
+ 
+ if(METIS_INSTALL)
+-  install(TARGETS metis
++  install(TARGETS metis EXPORT metisTargets
+     LIBRARY DESTINATION lib
+     RUNTIME DESTINATION bin
+-    ARCHIVE DESTINATION lib)
++    ARCHIVE DESTINATION lib
++    INCLUDES DESTINATION include)
+ endif()

config/vcpkg/ports/metis-mfem/portfile.cmake

@@ -0,0 +1,41 @@
+vcpkg_check_linkage(ONLY_STATIC_LIBRARY)
+set(OPTIONS -DSHARED=OFF)
+
+set(METIS_VERSION 5.1.0)
+
+vcpkg_download_distfile(ARCHIVE
+    URLS "https://github.com/mfem/tpls/raw/gh-pages/metis-${METIS_VERSION}.tar.gz"
+    FILENAME "metis-${METIS_VERSION}.tar.gz"
+    SHA512 deea47749d13bd06fbeaf98a53c6c0b61603ddc17a43dae81d72c8015576f6495fd83c11b0ef68d024879ed5415c14ebdbd87ce49c181bdac680573bea8bdb25
+)
+
+vcpkg_extract_source_archive_ex(
+    OUT_SOURCE_PATH SOURCE_PATH
+    ARCHIVE ${ARCHIVE}
+    REF ${METIS_VERSION}
+    PATCHES
+        enable-install.patch
+        disable-programs.patch
+        fix-runtime-install-destination.patch
+        fix-metis-vs14-math.patch
+        fix-gklib-vs14-math.patch
+        fix-linux-build-error.patch
+        install-metisConfig.patch
+        fix-INT_MIN_define.patch
+)
+
+vcpkg_configure_cmake(
+    SOURCE_PATH ${SOURCE_PATH}
+    PREFER_NINJA
+    OPTIONS ${OPTIONS}
+)
+
+vcpkg_install_cmake()
+vcpkg_copy_pdbs()
+vcpkg_fixup_cmake_targets(CONFIG_PATH lib/cmake/metis)
+
+file(REMOVE_RECURSE ${CURRENT_PACKAGES_DIR}/debug/include)
+
+# Handle copyright
+file(COPY ${SOURCE_PATH}/LICENSE.txt DESTINATION ${CURRENT_PACKAGES_DIR}/share/metis)
+file(INSTALL ${SOURCE_PATH}/LICENSE.txt DESTINATION ${CURRENT_PACKAGES_DIR}/share/${PORT} RENAME copyright)

config/vcpkg/ports/metis-mfem/vcpkg.json

@@ -0,0 +1,7 @@
+{
+  "name": "metis-mfem",
+  "version-string": "5.1.0",
+  "port-version": 0,
+  "description": "Serial Graph Partitioning and Fill-reducing Matrix Ordering",
+  "homepage": "https://glaros.dtc.umn.edu/gkhome/metis/metis/overview"
+}

doc/CodeDocumentation.conf.in

@@ -2349,7 +2349,7 @@ PLANTUML_INCLUDE_PATH  =
 # Minimum value: 0, maximum value: 10000, default value: 50.
 # This tag requires that the tag HAVE_DOT is set to YES.
 
-DOT_GRAPH_MAX_NODES    = 50
+DOT_GRAPH_MAX_NODES    = 100
 
 # The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the graphs
 # generated by dot. A depth value of 3 means that only nodes reachable from the

examples/ex1.cpp

@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               ex1 -pa -d cuda
-//               ex1 -fa -d cuda
+//             * ex1 -fa -d cuda
 //               ex1 -pa -d raja-cuda
 //             * ex1 -pa -d raja-hip
 //               ex1 -pa -d occa-cuda

examples/ex1p.cpp

@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               mpirun -np 4 ex1p -pa -d cuda
-//               mpirun -np 4 ex1p -fa -d cuda
+//             * mpirun -np 4 ex1p -fa -d cuda
 //               mpirun -np 4 ex1p -pa -d occa-cuda
 //               mpirun -np 4 ex1p -pa -d raja-omp
 //               mpirun -np 4 ex1p -pa -d ceed-cpu

fem/datacollection.cpp

@@ -1078,7 +1078,7 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
          it->second->GetValues(i, RefG->RefPts, val, pmat);
          for (int j = 0; j < val.Size(); j++)
          {
-            WriteBinaryOrASCII(out, buf, val(j), "\n", pv_data_format);
+            WriteBinaryOrASCII(os, buf, val(j), "\n", pv_data_format);
          }
       }
    }
@@ -1094,7 +1094,7 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
          {
             for (int ii = 0; ii < vval.Height(); ii++)
             {
-               WriteBinaryOrASCII(out, buf, vval(ii,jj), " ", pv_data_format);
+               WriteBinaryOrASCII(os, buf, vval(ii,jj), " ", pv_data_format);
             }
             if (pv_data_format == VTKFormat::ASCII) { os << '\n'; }
          }

fem/gridfunc.cpp

@@ -3653,278 +3653,6 @@ void GridFunction::ComputeElementLpErrors(const double p,
    }
 }
 
-double GridFunction::ComputeIntegral(Coefficient &weight,
-                                     int ir_order)
-{
-   double val = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   Vector vals;
-
-   for (int i = 0; i < fes->GetNE(); i++)
-   {
-      fe = fes->GetFE(i);
-      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
-      T = fes->GetElementTransformation(i);
-      GetValues(*T, *ir, vals);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-
-         val += ip.weight * T->Weight() * vals[j] * weight.Eval(*T, ip);
-      }
-   }
-
-   return val;
-}
-
-double GridFunction::ComputeIntegral(const Array<int> &attr_marker,
-                                     Coefficient &weight,
-                                     int ir_order)
-{
-   double val = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   Vector vals;
-
-   for (int i = 0; i < fes->GetNE(); i++)
-   {
-      int attr = fes->GetAttribute(i);
-      if (attr_marker[attr-1] == 0) { continue; }
-
-      fe = fes->GetFE(i);
-      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
-      T = fes->GetElementTransformation(i);
-      GetValues(*T, *ir, vals);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-
-         val += ip.weight * T->Weight() * vals[j] * weight.Eval(*T, ip);
-      }
-   }
-
-   return val;
-}
-
-double GridFunction::ComputeIntegral(VectorCoefficient &v_weight,
-                                     int ir_order)
-{
-   double val = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   DenseMatrix vals, v_weights;
-
-   for (int i = 0; i < fes->GetNE(); i++)
-   {
-      fe = fes->GetFE(i);
-      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
-      T = fes->GetElementTransformation(i);
-      GetVectorValues(*T, *ir, vals);
-
-      v_weight.Eval(v_weights, *T, *ir);
-
-      // column-wise dot product of the vector field values (in vals) and the
-      // vector weights (in v_weight)
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-
-         double val_ip = 0.0;
-         for (int d = 0; d < vals.Height(); d++)
-         {
-            val_ip += vals(d,j)*v_weights(d,j);
-         }
-
-         val += ip.weight * T->Weight() * val_ip;
-      }
-   }
-
-   return val;
-}
-
-double GridFunction::ComputeIntegral(const Array<int> &attr_marker,
-                                     VectorCoefficient &v_weight,
-                                     int ir_order)
-{
-   double val = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   DenseMatrix vals, v_weights;
-
-   for (int i = 0; i < fes->GetNE(); i++)
-   {
-      int attr = fes->GetAttribute(i);
-      if (attr_marker[attr-1] == 0) { continue; }
-
-      fe = fes->GetFE(i);
-      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
-      T = fes->GetElementTransformation(i);
-      GetVectorValues(*T, *ir, vals);
-
-      v_weight.Eval(v_weights, *T, *ir);
-
-      // column-wise dot product of the vector field values (in vals) and the
-      // vector weights (in v_weight)
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-
-         double val_ip = 0.0;
-         for (int d = 0; d < vals.Height(); d++)
-         {
-            val_ip += vals(d,j)*v_weights(d,j);
-         }
-
-         val += ip.weight * T->Weight() * val_ip;
-      }
-   }
-
-   return val;
-}
-
-double GridFunction::ComputeBdrIntegral(Coefficient &weight,
-                                        int ir_order)
-{
-   double val = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   Vector vals;
-
-   for (int i = 0; i < fes->GetNBE(); i++)
-   {
-      fe = fes->GetBE(i);
-      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
-      T = fes->GetBdrElementTransformation(i);
-      GetValues(*T, *ir, vals);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-
-         val += ip.weight * T->Weight() * vals[j] * weight.Eval(*T, ip);
-      }
-   }
-
-   return val;
-}
-
-double GridFunction::ComputeBdrIntegral(const Array<int> &bdr_attr_marker,
-                                        Coefficient &weight,
-                                        int ir_order)
-{
-   double val = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   Vector vals;
-
-   for (int i = 0; i < fes->GetNBE(); i++)
-   {
-      int attr = fes->GetBdrAttribute(i);
-      if (bdr_attr_marker[attr-1] == 0) { continue; }
-
-      fe = fes->GetBE(i);
-      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
-      T = fes->GetBdrElementTransformation(i);
-      GetValues(*T, *ir, vals);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-
-         val += ip.weight * T->Weight() * vals[j] * weight.Eval(*T, ip);
-      }
-   }
-
-   return val;
-}
-
-double GridFunction::ComputeBdrIntegral(VectorCoefficient &v_weight,
-                                        int ir_order)
-{
-   double val = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   DenseMatrix vals, v_weights;
-
-   for (int i = 0; i < fes->GetNBE(); i++)
-   {
-      fe = fes->GetBE(i);
-      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
-      T = fes->GetBdrElementTransformation(i);
-      GetVectorValues(*T, *ir, vals);
-
-      v_weight.Eval(v_weights, *T, *ir);
-
-      // column-wise dot product of the vector field values (in vals) and the
-      // vector weights (in v_weight)
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-
-         double val_ip = 0.0;
-         for (int d = 0; d < vals.Height(); d++)
-         {
-            val_ip += vals(d,j)*v_weights(d,j);
-         }
-
-         val += ip.weight * T->Weight() * val_ip;
-      }
-   }
-
-   return val;
-}
-
-double GridFunction::ComputeBdrIntegral(const Array<int> &bdr_attr_marker,
-                                        VectorCoefficient &v_weight,
-                                        int ir_order)
-{
-   double val = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   DenseMatrix vals, v_weights;
-
-   for (int i = 0; i < fes->GetNBE(); i++)
-   {
-      int attr = fes->GetBdrAttribute(i);
-      if (bdr_attr_marker[attr-1] == 0) { continue; }
-
-      fe = fes->GetBE(i);
-      const IntegrationRule *ir = &(IntRules.Get(fe->GetGeomType(), ir_order));
-      T = fes->GetBdrElementTransformation(i);
-      GetVectorValues(*T, *ir, vals);
-
-      v_weight.Eval(v_weights, *T, *ir);
-
-      // column-wise dot product of the vector field values (in vals) and the
-      // vector weights (in v_weight)
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-
-         double val_ip = 0.0;
-         for (int d = 0; d < vals.Height(); d++)
-         {
-            val_ip += vals(d,j)*v_weights(d,j);
-         }
-
-         val += ip.weight * T->Weight() * val_ip;
-      }
-   }
-
-   return val;
-}
-
 GridFunction & GridFunction::operator=(double value)
 {
    Vector::operator=(value);

fem/gridfunc.hpp

@@ -638,34 +638,6 @@ public:
                                        ) const
    { ComputeElementLpErrors(infinity(), exsol, error, NULL, NULL, irs); }
 
-   virtual double ComputeIntegral(Coefficient &weight,
-                                  int ir_order);
-
-   virtual double ComputeIntegral(const Array<int> &attr_marker,
-                                  Coefficient &weight,
-                                  int ir_order);
-
-   virtual double ComputeIntegral(VectorCoefficient &v_weight,
-                                  int ir_order);
-
-   virtual double ComputeIntegral(const Array<int> &attr_marker,
-                                  VectorCoefficient &v_weight,
-                                  int ir_order);
-
-   virtual double ComputeBdrIntegral(Coefficient &weight,
-                                     int ir_order);
-
-   virtual double ComputeBdrIntegral(const Array<int> &bdr_attr_marker,
-                                     Coefficient &weight,
-                                     int ir_order);
-
-   virtual double ComputeBdrIntegral(VectorCoefficient &v_weight,
-                                     int ir_order);
-
-   virtual double ComputeBdrIntegral(const Array<int> &bdr_attr_marker,
-                                     VectorCoefficient &v_weight,
-                                     int ir_order);
-
    virtual void ComputeFlux(BilinearFormIntegrator &blfi,
                             GridFunction &flux,
                             bool wcoef = true, int subdomain = -1);

fem/lor/lor.cpp

@@ -467,7 +467,8 @@ void LORDiscretization::FormLORSpace()
    mesh = new Mesh(Mesh::MakeRefined(mesh_ho, refinements, ref_type));
 
    fec = fes_ho.FEColl()->Clone(GetLOROrder());
-   fes = new FiniteElementSpace(mesh, fec);
+   const int vdim = fes_ho.GetVDim();
+   fes = new FiniteElementSpace(mesh, fec, vdim);
    SetupProlongationAndRestriction();
 }
 
@@ -511,7 +512,8 @@ void ParLORDiscretization::FormLORSpace()
    mesh = pmesh;
 
    fec = pfes_ho.FEColl()->Clone(GetLOROrder());
-   ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec);
+   const int vdim = fes_ho.GetVDim();
+   ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec, vdim);
    fes = pfes;
    SetupProlongationAndRestriction();
 }

fem/tbilininteg.hpp

@@ -361,8 +361,8 @@ struct TDiffusionKernel<2,2,complex_t>
        @param F Jt [M x Dim x SDim x NE] - Jacobian transposed, data member in F
        @param Q CoefficientEval<>::Type
        @param q CoefficientEval<>::Type::result_t
-       @param A [M x Dim*(Dim+1)/2] partially assembled Dim x Dim symm. matrices
-       @param A [M x Dim x Dim]       partially assembled Dim x Dim matrices
+       @param A either [M x Dim*(Dim+1)/2] partially assembled Dim x Dim symm.
+       matrices, or [M x Dim x Dim] partially assembled Dim x Dim matrices.
    */
    template <typename T_result_t, typename Q_t, typename q_t, typename asm_type>
    static inline MFEM_ALWAYS_INLINE

linalg/ode.cpp

@@ -188,7 +188,6 @@ void ExplicitRKSolver::Step(Vector &x, double &t, double &dt)
 
 ExplicitRKSolver::~ExplicitRKSolver()
 {
-   delete [] k;
 }
 
 const double RK6Solver::a[] =

linalg/ode.hpp

@@ -331,7 +331,6 @@ public:
    ~AdamsMoultonSolver()
    {
       if (RKsolver) { delete RKsolver; }
-      delete [] k;
    };
 };
 

miniapps/common/makefile

@@ -49,7 +49,7 @@ SEQ_MINIOBJS = mesh_extras.o fem_extras.o
 ifeq ($(MFEM_USE_MPI),NO)
    MINIOBJS = $(SEQ_MINIOBJS)
 else
-   MINIOBJS = $(SEQ_MINIOBJS) pmesh_extras.o pfem_extras.o
+   MINIOBJS = $(SEQ_MINIOBJS) pfem_extras.o
 endif
 
 .SUFFIXES:

miniapps/common/mesh_extras.cpp

@@ -136,324 +136,6 @@ ElementMeshStream::ElementMeshStream(Element::Type e)
 
 }
 
-double ComputeVolume(const Mesh &mesh, int ir_order)
-{
-   double vol = 0.0;
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNE(); i++)
-   {
-      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-
-         double w = T.Weight() * ip.weight;
-
-         vol += w;
-      }
-   }
-   return vol;
-}
-
-double ComputeVolume(const Mesh &mesh, const Array<int> &attr_marker,
-                     int ir_order)
-{
-   double vol = 0.0;
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNE(); i++)
-   {
-      int attr = mesh.GetAttribute(i);
-      if (attr_marker[attr-1] == 0) { continue; }
-
-      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-
-         double w = T.Weight() * ip.weight;
-
-         vol += w;
-      }
-   }
-   return vol;
-}
-
-double ComputeSurfaceArea(const Mesh &mesh, int ir_order)
-{
-   double area = 0.0;
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNBE(); i++)
-   {
-      // ElementTransformation *T = mesh.GetBdrElementTransformation(i);
-      const_cast<Mesh&>(mesh).GetBdrElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetBdrElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-
-         double w = T.Weight() * ip.weight;
-
-         area += w;
-      }
-   }
-   return area;
-}
-
-double ComputeSurfaceArea(const Mesh &mesh, const Array<int> &bdr_attr_marker,
-                          int ir_order)
-{
-   double area = 0.0;
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNBE(); i++)
-   {
-      int attr = mesh.GetBdrAttribute(i);
-      if (bdr_attr_marker[attr-1] == 0) { continue; }
-
-      const_cast<Mesh&>(mesh).GetBdrElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetBdrElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-
-         double w = T.Weight() * ip.weight;
-
-         area += w;
-      }
-   }
-   return area;
-}
-
-double ComputeZerothMoment(const Mesh &mesh, Coefficient &rho,
-                           int ir_order)
-{
-   double mom = 0.0;
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNE(); i++)
-   {
-      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-
-         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
-
-         mom += w;
-      }
-   }
-
-   return mom;
-}
-
-void ComputeElementZerothMoments(const Mesh &mesh, Coefficient &rho,
-                                 int ir_order, GridFunction &m)
-{
-   MFEM_ASSERT(m.Size() == mesh.GetNE(), "Invalid GridFunction.  "
-               "Must have a length equal to the number of mesh elements.");
-
-   FiniteElementSpace * fes = m.FESpace();
-
-   Array<int> vdofs;
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNE(); i++)
-   {
-      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      fes->GetElementVDofs(i, vdofs);
-
-      double mom = 0.0;
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-
-         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
-
-         mom += w;
-      }
-
-      m[vdofs[0]] = mom;
-   }
-}
-
-double ComputeFirstMoment(const Mesh &mesh, Coefficient &rho,
-                          int ir_order, Vector &mom)
-{
-   double mom0 = 0.0;
-
-   int sdim = mesh.SpaceDimension();
-
-   mom.SetSize(sdim);
-   mom = 0.0;
-
-   double x_data[3];
-   Vector x(x_data, sdim);
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNE(); i++)
-   {
-      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-         T.Transform(ip, x);
-
-         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
-
-         mom0 += w;
-         mom.Add(w, x);
-      }
-   }
-
-   return mom0;
-}
-
-double ComputeSecondMoment(const Mesh &mesh, Coefficient &rho,
-                           const Vector &center,
-                           int ir_order, DenseMatrix &mom)
-{
-   double mom0 = 0.0;
-
-   int sdim = mesh.SpaceDimension();
-
-   mom.SetSize(sdim);
-   mom = 0.0;
-
-   double x_data[3];
-   Vector x(x_data, sdim);
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNE(); i++)
-   {
-      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-         T.Transform(ip, x);
-         x.Add(-1.0, center);
-
-         double r2 = x * x;
-
-         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
-
-         mom0 += w;
-
-         for (int k=0; k<sdim; k++)
-         {
-            mom(k,k) += w * r2;
-            for (int l=k; l<sdim; l++)
-            {
-               mom(k,l) -= w * x[k] * x[l];
-            }
-         }
-      }
-   }
-
-   for (int k=0; k<sdim; k++)
-   {
-      for (int l=0; l<k; l++)
-      {
-         mom(k,l) = mom(l,k);
-      }
-   }
-
-   return mom0;
-}
-
-void ComputeElementCentersOfMass(const Mesh &mesh, Coefficient &rho,
-                                 int ir_order, GridFunction &c)
-{
-   int sdim = mesh.SpaceDimension();
-
-   MFEM_ASSERT(c.Size() == mesh.GetNE() * sdim,
-               "Invalid GridFunction.  Must have a length equal to the "
-               "number of mesh elements times the spatial dimension.");
-
-   FiniteElementSpace * fes = c.FESpace();
-
-   Array<int> vdofs;
-
-   double x_data[3];
-   Vector x(x_data, sdim);
-
-   c = 0.0;
-
-   IsoparametricTransformation T;
-
-   for (int i=0; i<mesh.GetNE(); i++)
-   {
-      const_cast<Mesh&>(mesh).GetElementTransformation(i, &T);
-      Geometry::Type geom = mesh.GetElementBaseGeometry(i);
-      const IntegrationRule *ir = &IntRules.Get(geom, ir_order);
-
-      fes->GetElementVDofs(i, vdofs);
-
-      double mom0 = 0.0;
-
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T.SetIntPoint(&ip);
-         T.Transform(ip, x);
-
-         double w = T.Weight() * ip.weight * rho.Eval(T, ip);
-
-         mom0 += w;
-
-         for (int k=0; k<sdim; k++)
-         {
-            c[vdofs[k]] += w * x[k];
-         }
-      }
-      for (int k=0; k<sdim; k++)
-      {
-         c[vdofs[k]] /= mom0;
-      }
-   }
-}
-
 void
 MergeMeshNodes(Mesh * mesh, int logging)
 {

miniapps/common/mesh_extras.hpp

@@ -27,42 +27,6 @@ public:
    ElementMeshStream(Element::Type e);
 };
 
-double ComputeVolume(const Mesh &mesh, int ir_order);
-double ComputeVolume(const Mesh &mesh, const Array<int> &attr_marker,
-                     int ir_order);
-double ComputeSurfaceArea(const Mesh &mesh, int ir_order);
-double ComputeSurfaceArea(const Mesh &mesh, const Array<int> &bdr_attr_marker,
-                          int ir_order);
-
-double ComputeZerothMoment(const Mesh &mesh, Coefficient &rho,
-                           int ir_order);
-double ComputeFirstMoment(const Mesh &mesh, Coefficient &rho,
-                          int ir_order, Vector &mom);
-double ComputeSecondMoment(const Mesh &mesh, Coefficient &rho,
-                           const Vector &center,
-                           int ir_order, DenseMatrix &mom);
-
-inline void ComputeNormalizedFirstMoment(const Mesh &mesh, Coefficient &rho,
-                                         int ir_order, Vector &mom)
-{
-   double mom0 = ComputeFirstMoment(mesh, rho, ir_order, mom);
-   mom /= mom0;
-}
-
-inline void ComputeNormalizedSecondMoment(const Mesh &mesh, Coefficient &rho,
-                                          const Vector &center,
-                                          int ir_order, DenseMatrix &mom)
-{
-   double mom0 = ComputeSecondMoment(mesh, rho, center, ir_order, mom);
-   mom *= 1.0 / mom0;
-}
-
-void ComputeElementZerothMoments(const Mesh &mesh, Coefficient &rho,
-                                 int ir_order, GridFunction &m);
-
-void ComputeElementCentersOfMass(const Mesh &mesh, Coefficient &rho,
-                                 int ir_order, GridFunction &c);
-
 /// Merges vertices which lie at the same location
 void MergeMeshNodes(Mesh * mesh, int logging);
 

miniapps/common/pmesh_extras.cpp

@@ -1,109 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-
-#include "pmesh_extras.hpp"
-#include "mesh_extras.hpp"
-
-#ifdef MFEM_USE_MPI
-
-using namespace std;
-
-namespace mfem
-{
-
-namespace common
-{
-
-double ComputeVolume(const ParMesh &pmesh, const Array<int> &attr_marker,
-                     int ir_order)
-{
-   double loc_vol = ComputeVolume(dynamic_cast<const Mesh&>(pmesh),
-                                  attr_marker, ir_order);
-   double glb_vol = 0.0;
-   MPI_Allreduce(&loc_vol, &glb_vol, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
-   return glb_vol;
-}
-
-double ComputeVolume(const ParMesh &pmesh, int ir_order)
-{
-   double loc_vol = ComputeVolume(dynamic_cast<const Mesh&>(pmesh), ir_order);
-   double glb_vol = 0.0;
-   MPI_Allreduce(&loc_vol, &glb_vol, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
-   return glb_vol;
-}
-
-double ComputeSurfaceArea(const ParMesh &pmesh, int ir_order)
-{
-   double loc_area = ComputeSurfaceArea(dynamic_cast<const Mesh&>(pmesh),
-                                        ir_order);
-   double glb_area = 0.0;
-   MPI_Allreduce(&loc_area, &glb_area, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
-   return glb_area;
-}
-
-double ComputeSurfaceArea(const ParMesh &pmesh,
-                          const Array<int> &bdr_attr_marker, int ir_order)
-{
-   double loc_area = ComputeSurfaceArea(dynamic_cast<const Mesh&>(pmesh),
-                                        bdr_attr_marker, ir_order);
-   double glb_area = 0.0;
-   MPI_Allreduce(&loc_area, &glb_area, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
-   return glb_area;
-}
-
-double ComputeZerothMoment(const ParMesh &pmesh, Coefficient &rho,
-                           int ir_order)
-{
-   double loc_mom = ComputeZerothMoment(dynamic_cast<const Mesh&>(pmesh),
-                                        rho, ir_order);
-   double glb_mom = 0.0;
-   MPI_Allreduce(&loc_mom, &glb_mom, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
-   return glb_mom;
-}
-
-double ComputeFirstMoment(const ParMesh &pmesh, Coefficient &rho,
-                          int ir_order, Vector &mom)
-{
-   int sdim = pmesh.SpaceDimension();
-   mom.SetSize(sdim);
-   double loc_mom_data[3];
-   Vector loc_mom(loc_mom_data, sdim);
-   double loc_mom0 = ComputeFirstMoment(dynamic_cast<const Mesh&>(pmesh),
-                                        rho, ir_order, loc_mom);
-   double glb_mom0 = 0.0;
-   MPI_Allreduce(&loc_mom0, &glb_mom0, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
-   MPI_Allreduce(loc_mom_data, mom.GetData(), sdim, MPI_DOUBLE, MPI_SUM,
-                 pmesh.GetComm());
-   return glb_mom0;
-}
-
-double ComputeSecondMoment(const ParMesh &pmesh, Coefficient &rho,
-                           const Vector &center,
-                           int ir_order, DenseMatrix &mom)
-{
-   int sdim = pmesh.SpaceDimension();
-   mom.SetSize(sdim);
-   double loc_mom_data[9];
-   DenseMatrix loc_mom(loc_mom_data, sdim, sdim);
-   double loc_mom0 = ComputeSecondMoment(dynamic_cast<const Mesh&>(pmesh),
-                                         rho, center, ir_order, loc_mom);
-   double glb_mom0 = 0.0;
-   MPI_Allreduce(&loc_mom0, &glb_mom0, 1, MPI_DOUBLE, MPI_SUM, pmesh.GetComm());
-   MPI_Allreduce(loc_mom_data, mom.GetData(), sdim * sdim, MPI_DOUBLE, MPI_SUM,
-                 pmesh.GetComm());
-   return glb_mom0;
-}
-
-} // namespace common
-
-} // namespace mfem
-
-#endif

miniapps/common/pmesh_extras.hpp

@@ -1,63 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-
-#ifndef MFEM_PMESH_EXTRAS
-#define MFEM_PMESH_EXTRAS
-
-#include "mfem.hpp"
-
-#ifdef MFEM_USE_MPI
-
-namespace mfem
-{
-
-namespace common
-{
-
-double ComputeVolume(const ParMesh &pmesh, int ir_order);
-double ComputeVolume(const ParMesh &mesh, const Array<int> &attr_marker,
-                     int ir_order);
-double ComputeSurfaceArea(const ParMesh &pmesh, int ir_order);
-double ComputeSurfaceArea(const ParMesh &mesh,
-                          const Array<int> &bdr_attr_marker,
-                          int ir_order);
-
-double ComputeZerothMoment(const ParMesh &pmesh, Coefficient &rho,
-                           int ir_order);
-double ComputeFirstMoment(const ParMesh &pmesh, Coefficient &rho,
-                          int ir_order, Vector &mom);
-double ComputeSecondMoment(const ParMesh &pmesh, Coefficient &rho,
-                           const Vector &center,
-                           int ir_order, DenseMatrix &mom);
-
-inline void ComputeNormalizedFirstMoment(const ParMesh &pmesh,
-                                         Coefficient &rho,
-                                         int ir_order, Vector &mom)
-{
-   double mom0 = ComputeFirstMoment(pmesh, rho, ir_order, mom);
-   mom /= mom0;
-}
-
-inline void ComputeNormalizedSecondMoment(const ParMesh &pmesh,
-                                          Coefficient &rho,
-                                          const Vector &center,
-                                          int ir_order, DenseMatrix &mom)
-{
-   double mom0 = ComputeSecondMoment(pmesh, rho, center, ir_order, mom);
-   mom *= 1.0 / mom0;
-}
-
-} // namespace common
-
-} // namespace mfem
-
-#endif // MFEM_USE_MPI
-#endif

miniapps/tools/makefile

@@ -25,10 +25,8 @@ include $(DEFAULTS_MK)
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-
-SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \
-	seq-moment
-PAR_MINIAPPS = par-moment
+SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer
+PAR_MINIAPPS =
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -40,7 +38,6 @@ endif
 .PHONY: all lib-common clean clean-build clean-exec
 .PRECIOUS: %.o
 
-#COMMON_O=../common/fem_extras.o ../common/mesh_extras.o ../common/pmesh_extras.o
 COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
@@ -63,15 +60,6 @@ display-basis: %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
 
-
-seq-moment: %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
-
-par-moment: %: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
-
 # Rule for building lib-common
 lib-common:
 	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common

miniapps/tools/par-moment.cpp

@@ -1,181 +0,0 @@
-#include "mfem.hpp"
-#include "../common/mesh_extras.hpp"
-#include "../common/pmesh_extras.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-using namespace mfem::common;
-
-static int prob_ = -1;
-
-double densityFunc(const Vector &x)
-{
-   switch (prob_)
-   {
-      case 0:
-         // Linear in the radius
-         return sqrt(x * x);
-      case 1:
-         // Hydrostatic equilibrium (1D)
-         return pow(cosh(sqrt(0.5 * (x * x))), -2.0);
-      case 2:
-         // Hydrostatic equilibrium (2D)
-         return pow(1.0 + 0.125 * (x * x), -2.0);
-      case 3:
-         // Hydrostatic equilibrium (3D), approx.
-         return pow(1.0 + (x * x) / 9.0, -1.5);
-      case 4:
-         // Off-center Gaussian
-      {
-         double s_data[3];
-         Vector s(s_data, x.Size());
-         s = 0.0; s[0] = -1.0; s += x;
-         return exp(-(s * s));
-      }
-   }
-   // Default to homogeneous
-   return 1.0;
-}
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../../data/ball-nurbs.mesh";
-   int ser_ref_levels = 2;
-   int par_ref_levels = 0;
-   int ir_order = 2;
-   bool visualization = true;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&ir_order, "-o", "--order",
-                  "Integration rule order.");
-   args.AddOption(&prob_, "-d", "--density",
-                  "Density profile:\n"
-                  "  0 - Linear increase with radius,\n"
-                  "  1 - Hydrostatic equilibrium in 1D,\n"
-                  "  2 - Hydrostatic equilibrium in 2D,\n"
-                  "  3 - Hydrostatic equilibrium in 3D (approximately),\n"
-                  "  4 - Gaussian centered at x = 1,\n"
-                  "  Default - homogeneous.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      MPI_Finalize();
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // 2. Read the mesh from the given mesh file. We can handle triangular,
-   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
-   //    the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-   int sdim = mesh->SpaceDimension();
-
-   // 3. Refine the mesh to increase the resolution. In this example we do
-   //    'ref_levels' of uniform refinement where 'ref_levels' is user defined.
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution. Once the
-   //    parallel mesh is defined, the serial mesh can be deleted.
-   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh->UniformRefinement();
-   }
-
-   FunctionCoefficient density(densityFunc);
-
-   L2_FECollection fec(0, dim);
-   ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
-   ParFiniteElementSpace *vfes = new ParFiniteElementSpace(pmesh, &fec, sdim);
-   ParGridFunction elemMass(fes);
-   ParGridFunction elemCent(vfes);
-
-   double vol = ComputeVolume(*pmesh, ir_order);
-   double area = ComputeSurfaceArea(*pmesh, ir_order);
-
-   Vector mom1(sdim);
-   Vector cent(sdim); cent = 0.0;
-   DenseMatrix mom2(sdim);
-
-   double mass  = ComputeZerothMoment(*pmesh, density, ir_order);
-   double mass1 = ComputeFirstMoment(*pmesh, density, ir_order, mom1);
-   double mass2 = ComputeSecondMoment(*pmesh, density, cent,
-                                      ir_order, mom2);
-
-   if (myid == 0)
-   {
-      cout << "Volume:       " << vol << endl;
-      cout << "Surface Area: " << area << endl;
-      cout << "Mass:         " << mass << endl;
-      cout << "Mass:         " << mass1 << endl;
-      cout << "Mass:         " << mass2 << endl;
-
-      cout << "First Moment:   "; mom1.Print(cout);
-      mom1 /= mass1;
-      cout << "Center of Mass: "; mom1.Print(cout);
-      cout << "Second Moment (moment of inertia):\n";
-      mom2.Print(cout);
-   }
-
-   ComputeElementZerothMoments(*pmesh, density, ir_order, elemMass);
-   ComputeElementCentersOfMass(*pmesh, density, ir_order, elemCent);
-
-   // 15. Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream m_sock(vishost, visport);
-      m_sock << "parallel " << num_procs << " " << myid << "\n";
-      m_sock.precision(8);
-      m_sock << "solution\n" << *pmesh << elemMass
-             << "window_title 'Element Masses'" << flush;
-
-      socketstream c_sock(vishost, visport);
-      c_sock << "parallel " << num_procs << " " << myid << "\n";
-      c_sock.precision(8);
-      c_sock << "solution\n" << *pmesh << elemCent
-             << "window_title 'Element Centers'"
-             << "keys vvv" << flush;
-   }
-
-   // 16. Free the used memory.
-   delete fes;
-   delete vfes;
-   delete pmesh;
-
-   MPI_Finalize();
-
-   return 0;
-}

miniapps/tools/seq-moment.cpp

@@ -1,144 +0,0 @@
-#include "mfem.hpp"
-#include "../common/mesh_extras.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-using namespace mfem::common;
-
-static int prob_ = -1;
-
-double densityFunc(const Vector &x)
-{
-   switch (prob_)
-   {
-      case 0:
-         // Linear in the radius
-         return sqrt(x * x);
-      case 1:
-         // Hydrostatic equilibrium (1D)
-         return pow(cosh(sqrt(0.5 * (x * x))), -2.0);
-      case 2:
-         // Hydrostatic equilibrium (2D)
-         return pow(1.0 + 0.125 * (x * x), -2.0);
-      case 3:
-         // Hydrostatic equilibrium (3D), approx.
-         return pow(1.0 + (x * x) / 9.0, -1.5);
-      case 4:
-         // Off-center Gaussian
-      {
-         double s_data[3];
-         Vector s(s_data, x.Size());
-         s = 0.0; s[0] = -1.0; s += x;
-         return exp(-(s * s));
-      }
-   }
-   // Default to homogeneous
-   return 1.0;
-}
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options.
-   const char *mesh_file = "../../data/ball-nurbs.mesh";
-   int ref_levels = 2;
-   int ir_order = 2;
-   bool visualization = true;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ref_levels, "-r", "--refine",
-                  "Number of times to refine the mesh uniformly.");
-   args.AddOption(&ir_order, "-o", "--order",
-                  "Integration rule order.");
-   args.AddOption(&prob_, "-d", "--density",
-                  "Density profile:\n"
-                  "  0 - Linear increase with radius,\n"
-                  "  1 - Hydrostatic equilibrium in 1D,\n"
-                  "  2 - Hydrostatic equilibrium in 2D,\n"
-                  "  3 - Hydrostatic equilibrium in 3D (approximately),\n"
-                  "  4 - Gaussian centered at x = 1,\n"
-                  "  Default - homogeneous.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   args.PrintOptions(cout);
-
-   // 2. Read the mesh from the given mesh file. We can handle triangular,
-   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
-   //    the same code.
-   Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-   int sdim = mesh.SpaceDimension();
-
-   // 3. Refine the mesh to increase the resolution. In this example we do
-   //    'ref_levels' of uniform refinement where 'ref_levels' is user defined.
-   {
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   FunctionCoefficient density(densityFunc);
-
-   L2_FECollection fec(0, dim);
-   FiniteElementSpace fes(&mesh, &fec);
-   FiniteElementSpace vfes(&mesh, &fec, sdim);
-   GridFunction elemMass(&fes);
-   GridFunction elemCent(&vfes);
-
-   double vol = ComputeVolume(mesh, ir_order);
-   double area = ComputeSurfaceArea(mesh, ir_order);
-
-   Vector mom1(sdim);
-   Vector cent(sdim); cent = 0.0;
-   DenseMatrix mom2(sdim);
-
-   double mass  = ComputeZerothMoment(mesh, density, ir_order);
-   double mass1 = ComputeFirstMoment(mesh, density, ir_order, mom1);
-   double mass2 = ComputeSecondMoment(mesh, density, cent,
-                                      ir_order, mom2);
-
-   cout << "Volume:       " << vol << endl;
-   cout << "Surface Area: " << area << endl;
-   cout << "Mass:         " << mass << endl;
-   cout << "Mass:         " << mass1 << endl;
-   cout << "Mass:         " << mass2 << endl;
-
-   cout << "First Moment:   "; mom1.Print(cout);
-   mom1 /= mass1;
-   cout << "Center of Mass: "; mom1.Print(cout);
-   cout << "Second Moment (moment of inertia):\n";
-   mom2.Print(cout);
-
-   ComputeElementZerothMoments(mesh, density, ir_order, elemMass);
-   ComputeElementCentersOfMass(mesh, density, ir_order, elemCent);
-
-   // 15. Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream m_sock(vishost, visport);
-      m_sock.precision(8);
-      m_sock << "solution\n" << mesh << elemMass
-             << "window_title 'Element Masses'" << flush;
-
-      socketstream c_sock(vishost, visport);
-      c_sock.precision(8);
-      c_sock << "solution\n" << mesh << elemCent
-             << "window_title 'Element Centers'"
-             << "keys vvv" << flush;
-   }
-
-   return 0;
-}

tests/unit/CMakeLists.txt

@@ -67,7 +67,6 @@ set(UNIT_TESTS_SRCS
   fem/test_fe.cpp
   fem/test_get_value.cpp
   fem/test_getderivative.cpp
-  fem/test_gridfunction_integrals.cpp
   fem/test_intrules.cpp
   fem/test_intruletypes.cpp
   fem/test_inversetransform.cpp

tests/unit/fem/test_gridfunction_integrals.cpp

@@ -1,232 +0,0 @@
-// Copyright (c) 2010-2022, 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.
-
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-
-using namespace mfem;
-
-namespace gridfunction_integrals
-{
-
-double func_1D_lin(const Vector &x)
-{
-   return x[0] - 0.5;
-}
-
-double wgt_1D_lin(const Vector &x)
-{
-   return 7.0 - 3.0 * x[0];
-}
-
-double func_1D_gaussian(const Vector &x)
-{
-   return sqrt(50.0/M_PI) * exp(-50.0 * pow(x[0] - 0.5, 2.0));
-}
-
-class FourierCosine1D : public Coefficient
-{
-private:
-   int n_;
-   double l_;
-   mutable Vector x_;
-
-public:
-   FourierCosine1D(double l) : n_(0), l_(l), x_(1) {}
-
-   void SetMode(int n) { n_ = n; }
-
-   double Eval(ElementTransformation &T, const IntegrationPoint &ip)
-   {
-      T.Transform(ip, x_);
-
-      return cos(2.0 * M_PI * (double)n_ * x_[0] / l_);
-   }
-};
-
-TEST_CASE("1D GridFunction::ComputeIntegral",
-          "[GridFunction]"
-          "[GridFunction::ComputeIntegral]")
-{
-   int n = 10;
-   int dim = 1;
-   int order = 1;
-
-   FunctionCoefficient funcCoef(func_1D_lin);
-   ConstantCoefficient wgt0Coef(1.0);
-   FunctionCoefficient wgt1Coef(wgt_1D_lin);
-
-   for (int type = (int)Element::SEGMENT;
-        type <= (int)Element::SEGMENT; type++)
-   {
-      Mesh mesh = Mesh::MakeCartesian1D(n, 2.0);
-
-      H1_FECollection h1_fec(order, dim);
-      DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
-                              FiniteElement::VALUE);
-      DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
-                              FiniteElement::INTEGRAL);
-
-      FiniteElementSpace h1_fespace(&mesh, &h1_fec);
-      FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
-      FiniteElementSpace dgi_fespace(&mesh, &dgi_fec);
-
-      GridFunction h1_x(&h1_fespace);
-      GridFunction dgv_x(&dgv_fespace);
-      GridFunction dgi_x(&dgi_fespace);
-
-      h1_x.ProjectCoefficient(funcCoef);
-      dgv_x.ProjectCoefficient(funcCoef);
-      dgi_x.ProjectCoefficient(funcCoef);
-
-      // First integrate with a weight of 1
-      double w0 = 1.0;
-
-      double w0_h1 = h1_x.ComputeIntegral(wgt0Coef, 0);
-      double w0_dgv = dgv_x.ComputeIntegral(wgt0Coef, 0);
-      double w0_dgi = dgi_x.ComputeIntegral(wgt0Coef, 0);
-
-      REQUIRE(w0_h1  == MFEM_Approx(w0));
-      REQUIRE(w0_dgv == MFEM_Approx(w0));
-      REQUIRE(w0_dgi == MFEM_Approx(w0));
-
-      // Integrate with a linear weight function
-      double w1 = 2.0;
-
-      double w1_h1  = h1_x.ComputeIntegral(wgt1Coef, 2);
-      double w1_dgv = dgv_x.ComputeIntegral(wgt1Coef, 2);
-      double w1_dgi = dgi_x.ComputeIntegral(wgt1Coef, 2);
-
-      REQUIRE(w1_h1  == MFEM_Approx(w1));
-      REQUIRE(w1_dgv == MFEM_Approx(w1));
-      REQUIRE(w1_dgi == MFEM_Approx(w1));
-   }
-}
-
-TEST_CASE("1D GridFunction::ComputeIntegral (Fourier)",
-          "[GridFunction]"
-          "[GridFunction::ComputeIntegral]")
-{
-   int n = 20;
-   int dim = 1;
-   int order = 3;
-   double l = 1.0;
-
-   FunctionCoefficient funcCoef(func_1D_gaussian);
-   FourierCosine1D wgtCoef(l);
-
-   for (int type = (int)Element::SEGMENT;
-        type <= (int)Element::SEGMENT; type++)
-   {
-      Mesh mesh = Mesh::MakeCartesian1D(n, l);
-
-      H1_FECollection h1_fec(order, dim);
-      DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
-                              FiniteElement::VALUE);
-      DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
-                              FiniteElement::INTEGRAL);
-
-      FiniteElementSpace h1_fespace(&mesh, &h1_fec);
-      FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
-      FiniteElementSpace dgi_fespace(&mesh, &dgi_fec);
-
-      GridFunction h1_x(&h1_fespace);
-      GridFunction dgv_x(&dgv_fespace);
-      GridFunction dgi_x(&dgi_fespace);
-
-      h1_x.ProjectCoefficient(funcCoef);
-      dgv_x.ProjectCoefficient(funcCoef);
-      dgi_x.ProjectCoefficient(funcCoef);
-
-      LinearForm h1_lf(&h1_fespace);
-      LinearForm dgv_lf(&dgv_fespace);
-      LinearForm dgi_lf(&dgi_fespace);
-
-      h1_lf.AddDomainIntegrator(new DomainLFIntegrator(wgtCoef));
-      dgv_lf.AddDomainIntegrator(new DomainLFIntegrator(wgtCoef));
-      dgi_lf.AddDomainIntegrator(new DomainLFIntegrator(wgtCoef));
-
-      // First integrate with a weight of 1
-      double w0 = 1.0;
-
-      h1_lf.Assemble();
-      dgi_lf.Assemble();
-      dgv_lf.Assemble();
-
-      double i1_w0_h1  = h1_lf(h1_x);
-      double i1_w0_dgv = dgv_lf(dgv_x);
-      double i1_w0_dgi = dgi_lf(dgi_x);
-
-      REQUIRE(i1_w0_h1  == MFEM_Approx(w0, 1e-6));
-      REQUIRE(i1_w0_dgv == MFEM_Approx(w0, 1e-6));
-      REQUIRE(i1_w0_dgi == MFEM_Approx(w0, 1e-6));
-
-      double i2_w0_h1  = h1_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-      double i2_w0_dgv = dgv_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-      double i2_w0_dgi = dgi_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-
-      REQUIRE(i2_w0_h1  == MFEM_Approx(i1_w0_h1));
-      REQUIRE(i2_w0_dgv == MFEM_Approx(i1_w0_dgv));
-      REQUIRE(i2_w0_dgi == MFEM_Approx(i1_w0_dgi));
-
-      // Integrate with a weight of cos(2 pi x)
-      wgtCoef.SetMode(1);
-
-      double w1 = -exp(-M_PI * M_PI / 50.0);
-
-      h1_lf.Assemble();
-      dgi_lf.Assemble();
-      dgv_lf.Assemble();
-
-      double i1_w1_h1  = h1_lf(h1_x);
-      double i1_w1_dgv = dgv_lf(dgv_x);
-      double i1_w1_dgi = dgi_lf(dgi_x);
-
-      REQUIRE(i1_w1_h1  == MFEM_Approx(w1, 1e-6));
-      REQUIRE(i1_w1_dgv == MFEM_Approx(w1, 1e-6));
-      REQUIRE(i1_w1_dgi == MFEM_Approx(w1, 1e-6));
-
-      double i2_w1_h1  = h1_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-      double i2_w1_dgv = dgv_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-      double i2_w1_dgi = dgi_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-
-      REQUIRE(i2_w1_h1  == MFEM_Approx(i1_w1_h1));
-      REQUIRE(i2_w1_dgv == MFEM_Approx(i1_w1_dgv));
-      REQUIRE(i2_w1_dgi == MFEM_Approx(i1_w1_dgi));
-
-      // Integrate with a weight of cos(4 pi x)
-      wgtCoef.SetMode(2);
-
-      double w2 = exp(-4.0 * M_PI * M_PI / 50.0);
-
-      h1_lf.Assemble();
-      dgi_lf.Assemble();
-      dgv_lf.Assemble();
-
-      double i1_w2_h1  = h1_lf(h1_x);
-      double i1_w2_dgv = dgv_lf(dgv_x);
-      double i1_w2_dgi = dgi_lf(dgi_x);
-
-      REQUIRE(i1_w2_h1  == MFEM_Approx(w2, 1e-6));
-      REQUIRE(i1_w2_dgv == MFEM_Approx(w2, 1e-6));
-      REQUIRE(i1_w2_dgi == MFEM_Approx(w2, 1e-6));
-
-      double i2_w2_h1  = h1_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-      double i2_w2_dgv = dgv_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-      double i2_w2_dgi = dgi_x.ComputeIntegral(wgtCoef, 2 * order + 1);
-
-      REQUIRE(i2_w2_h1  == MFEM_Approx(i1_w2_h1));
-      REQUIRE(i2_w2_dgv == MFEM_Approx(i1_w2_dgv));
-      REQUIRE(i2_w2_dgi == MFEM_Approx(i1_w2_dgi));
-   }
-}
-
-} // gridfunction_integrals