Adding unit test of GridFunction integrals in 1D

# Conflicts:
#	miniapps/common/CMakeLists.txt
#	miniapps/common/mesh_extras.cpp
#	miniapps/common/mesh_extras.hpp
#	miniapps/tools/makefile

.github/workflows/build-container.yaml

@@ -1,4 +1,4 @@
-name: "Docker"
+name: Build Deploy Container
 
 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: "Tests"
+name: builds-and-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-latest, macos-latest, windows-latest]
+        os: [ubuntu-20.04, macos-10.15, windows-2022]
         target: [dbg, opt]
         mpi: [seq, par]
         build-system: [make, cmake]
         hypre-target: [int32]
         exclude:
-          - os: ubuntu-latest
+          - os: ubuntu-20.04
             build-system: cmake
-          - os: macos-latest
+          - os: macos-10.15
             build-system: cmake
-          - os: windows-latest
+          - os: windows-2022
             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-latest
+          - os: windows-2022
             codecov: NO
-          - os: ubuntu-latest
+          - os: ubuntu-20.04
             target: opt
             codecov: NO
             mpi: par
             build-system: cmake
             hypre-target: int32
-          - os: ubuntu-latest
+          - os: ubuntu-20.04
             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-latest'
+      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-20.04'
       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-latest'
+      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-20.04'
       run: |
         sudo apt-get install lcov
 
     - name: Set up Homebrew
-      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
+      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-10.15'
       uses: Homebrew/actions/setup-homebrew@c4aafe8c4620bf08883dd4679c374f11e73329d3
 
     - name: get MPI (MacOS)
-      if: matrix.mpi == 'par' && matrix.os == 'macos-latest'
+      if: matrix.mpi == 'par' && matrix.os == 'macos-10.15'
       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-latest'
+      if: matrix.codecov == 'YES' && matrix.os == 'macos-10.15'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
         brew install lcov
 
     - name: get MPI (Windows)
-      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
+      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
       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-latest'
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-2022'
       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-latest'
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-2022'
       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-latest'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-2022'
       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-latest' && steps.metis-cache.outputs.cache-hit != 'true'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-2022' && steps.metis-cache.outputs.cache-hit != 'true'
       uses: mfem/github-actions/build-metis@v2.2
       with:
         archive: ${{ env.METIS_ARCHIVE }}
@@ -196,16 +196,20 @@ jobs:
         key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
 
     - name: prepare binary cache location
-      if: matrix.os == 'windows-latest' && steps.vcpkg-cache.outputs.cache-hit != 'true'
+      if: matrix.os == 'windows-2022' && steps.vcpkg-cache.outputs.cache-hit != 'true'
       run: |
         mkdir -p vcpkg_cache
 
     - name: install metis (Windows)
-      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
+      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
       env:
         VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
       run: |
-        vcpkg install metis-mfem --triplet=x64-windows-static --overlay-ports=${{ env.MFEM_TOP_DIR }}/config/vcpkg/ports
+        $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
 
     # MFEM build and test
     - name: build
@@ -248,7 +252,7 @@ jobs:
       shell: bash
 
     - name: cmake unit tests (Ubuntu 20.04)
-      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-latest'
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-20.04'
       run: |
         CTEST_CONFIG="Release"
         [[ ${{ matrix.target }} == 'dbg' ]] && CTEST_CONFIG="Debug"
@@ -256,7 +260,7 @@ jobs:
       shell: bash
 
     - name: cmake tests
-      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-latest'
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-20.04'
       run: |
         CTEST_CONFIG="Release"
         cd ${{ env.MFEM_TOP_DIR }}/build && ctest --output-on-failure -C ${CTEST_CONFIG}

.github/workflows/codeql-analysis.yml

@@ -1,71 +0,0 @@
-# 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-latest
+    runs-on: ubuntu-18.04
 
     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: "Checks"
+name: repo-check
 
 permissions:
   actions: write
@@ -28,7 +28,7 @@ on:
 
 jobs:
   file-headers-check:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-18.04
     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-latest
+    runs-on: ubuntu-18.04
     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
+        sudo apt-get install astyle=3.1-1ubuntu2
 
     - name: style check
       run: |
         ./config/githooks/pre-push --style
 
   documentation:
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-18.04
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -106,8 +106,6 @@ 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: |
@@ -120,7 +118,7 @@ jobs:
       github.ref != 'refs/heads/master' &&
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
-    runs-on: ubuntu-latest
+    runs-on: ubuntu-18.04
     steps:
     - name: checkout mfem
       uses: actions/checkout@v2

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

@@ -1,8 +0,0 @@
---- 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

@@ -1,15 +0,0 @@
---- 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

@@ -1,34 +0,0 @@
-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

@@ -1,11 +0,0 @@
---- 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

@@ -1,14 +0,0 @@
-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

@@ -1,11 +0,0 @@
---- 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

@@ -1,10 +0,0 @@
---- 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

@@ -1,44 +0,0 @@
-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

@@ -1,41 +0,0 @@
-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

@@ -1,7 +0,0 @@
-{
-  "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    = 100
+DOT_GRAPH_MAX_NODES    = 50
 
 # 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(os, buf, val(j), "\n", pv_data_format);
+            WriteBinaryOrASCII(out, 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(os, buf, vval(ii,jj), " ", pv_data_format);
+               WriteBinaryOrASCII(out, buf, vval(ii,jj), " ", pv_data_format);
             }
             if (pv_data_format == VTKFormat::ASCII) { os << '\n'; }
          }

fem/gridfunc.cpp

@@ -3653,6 +3653,278 @@ 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,6 +638,34 @@ 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,8 +467,7 @@ void LORDiscretization::FormLORSpace()
    mesh = new Mesh(Mesh::MakeRefined(mesh_ho, refinements, ref_type));
 
    fec = fes_ho.FEColl()->Clone(GetLOROrder());
-   const int vdim = fes_ho.GetVDim();
-   fes = new FiniteElementSpace(mesh, fec, vdim);
+   fes = new FiniteElementSpace(mesh, fec);
    SetupProlongationAndRestriction();
 }
 
@@ -512,8 +511,7 @@ void ParLORDiscretization::FormLORSpace()
    mesh = pmesh;
 
    fec = pfes_ho.FEColl()->Clone(GetLOROrder());
-   const int vdim = fes_ho.GetVDim();
-   ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec, vdim);
+   ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec);
    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 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.
+       @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
    */
    template <typename T_result_t, typename Q_t, typename q_t, typename asm_type>
    static inline MFEM_ALWAYS_INLINE

linalg/ode.cpp

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

linalg/ode.hpp

@@ -331,6 +331,7 @@ 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) pfem_extras.o
+   MINIOBJS = $(SEQ_MINIOBJS) pmesh_extras.o pfem_extras.o
 endif
 
 .SUFFIXES:

miniapps/common/mesh_extras.cpp

@@ -136,6 +136,324 @@ 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,6 +27,42 @@ 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

@@ -0,0 +1,109 @@
+// 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

@@ -0,0 +1,63 @@
+// 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,8 +25,10 @@ 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
-PAR_MINIAPPS =
+
+SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \
+	seq-moment
+PAR_MINIAPPS = par-moment
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -38,6 +40,7 @@ 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
@@ -60,6 +63,15 @@ 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

@@ -0,0 +1,181 @@
+#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

@@ -0,0 +1,144 @@
+#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,6 +67,7 @@ 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

@@ -0,0 +1,232 @@
+// 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