merge with master and resolve conflicts

.github/workflows/stale.yml

@@ -29,47 +29,3 @@ jobs:
         operations-per-run: 500
         exempt-issue-labels: "bug,WIP,ready-for-review,in-review,in-next"
         exempt-pr-labels: "bug,WIP,ready-for-review,in-review,in-next"
-
-  # Stale action for PRs with "in-review" label.
-  stale-in-review-pr:
-
-    runs-on: ubuntu-latest
-    permissions:
-      issues: write
-      pull-requests: write
-      actions: write
-
-    steps:
-    - uses: actions/stale@v9
-      with:
-        repo-token: ${{ secrets.GITHUB_TOKEN }}
-        stale-pr-message: ':warning: This PR has been automatically marked as stale because it has not had any activity in the last 150 days. *If no activity occurs in the next 30 days, it will be automatically closed.* Thank you for your contributions.'
-        only-pr-labels: "in-review"
-        days-before-pr-stale: 150
-        days-before-pr-close: 30
-        days-before-issue-stale: -1
-        days-before-issue-close: -1
-        stale-pr-label: 'stale'
-        operations-per-run: 500
-
-  # Stale action for PRs with "WIP" label.
-  stale-wip-pr:
-
-    runs-on: ubuntu-latest
-    permissions:
-      issues: write
-      pull-requests: write
-      actions: write
-
-    steps:
-    - uses: actions/stale@v9
-      with:
-        repo-token: ${{ secrets.GITHUB_TOKEN }}
-        stale-pr-message: ':warning: This PR has been automatically marked as stale because it has not had any activity in the last 300 days. *If no activity occurs in the next 30 days, it will be automatically closed.* Thank you for your contributions.'
-        only-pr-labels: "WIP"
-        days-before-pr-stale: 300
-        days-before-pr-close: 30
-        days-before-issue-stale: -1
-        days-before-issue-close: -1
-        stale-pr-label: 'stale'
-        operations-per-run: 500

.gitignore

@@ -208,13 +208,10 @@ miniapps/electromagnetics/volta
 miniapps/electromagnetics/tesla
 miniapps/electromagnetics/maxwell
 miniapps/electromagnetics/joule
-miniapps/electromagnetics/lorentz
 miniapps/electromagnetics/Volta-AMR*
 miniapps/electromagnetics/Tesla-AMR*
 miniapps/electromagnetics/Maxwell-Parallel*
 miniapps/electromagnetics/Joule_[0-9]*
-miniapps/electromagnetics/Lorentz_[0-9]*
-miniapps/electromagnetics/Lorentz.dat
 
 miniapps/gslib/field-diff
 miniapps/gslib/field-interp

CHANGELOG

@@ -91,10 +91,6 @@ New and updated examples and miniapps
 - Added a new miniapp (tools/gridfunction-bounds) to compute piecewise linear
   bounds on a given high-order grid function.
 
-- Added a new miniapp (electromagnetics/lorentz) which computes the trajectory
-  of a charged particle, subject to Lorentz forces, in electrostatic and/or
-  magnetostatic fields as computed by the volta or tesla miniapps.
-
 API changes:
 -----------
 - mfem::internal::tensor and mfem::internal::dual have been moved to

CMakeLists.txt

@@ -278,11 +278,6 @@ if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
   endif()
 endif()
 
-# Umpire (must be included before hypre, so hypre can use it if needed)
-if (MFEM_USE_UMPIRE)
-   find_package(UMPIRE REQUIRED)
-endif()
-
 # MPI -> hypre; PETSc (optional)
 if (MFEM_USE_MPI)
   find_package(MPI REQUIRED)
@@ -507,6 +502,11 @@ if (MFEM_USE_RAJA)
    find_package(RAJA REQUIRED)
 endif()
 
+# UMPIRE
+if (MFEM_USE_UMPIRE)
+   find_package(UMPIRE REQUIRED)
+endif()
+
 # GOOGLE-BENCHMARK
 if (MFEM_USE_BENCHMARK)
    find_package(Benchmark REQUIRED)

config/cmake/modules/FindHYPRE.cmake

@@ -38,91 +38,51 @@ if (HYPRE_FOUND OR TARGET HYPRE)
 endif()
 
 if (HYPRE_FETCH OR FETCH_TPLS)
-  # Collect all HYPRE_ENABLE variables and pass them to hypre, assuming they are BOOL.
-  set(HYPRE_CMAKE_OPTIONS "")
-  get_cmake_property(all_vars VARIABLES)
-  foreach(var ${all_vars})
-    if(var MATCHES "^HYPRE_ENABLE")
-      list(APPEND HYPRE_CMAKE_OPTIONS "-D${var}:BOOL=${${var}}")
-    endif()
-  endforeach()
-
   set(HYPRE_FETCH_VERSION 2.33.0)
-  set(HYPRE_FETCH_TAG "v${HYPRE_FETCH_VERSION}" CACHE STRING "Tag, branch, or commit for HYPRE")
   add_library(HYPRE STATIC IMPORTED)
   # set options and associated dependencies
-  list(APPEND HYPRE_CMAKE_OPTIONS -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE})
+  set(CMAKE_OPTIONS)
+  list(APPEND CMAKE_OPTIONS -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE})
   if (MFEM_USE_CUDA)
-    list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_CUDA:BOOL=ON -DCMAKE_CUDA_ARCHITECTURES:STRING=${CMAKE_CUDA_ARCHITECTURES})
+    list(APPEND CMAKE_OPTIONS -DHYPRE_WITH_CUDA:BOOL=ON)
     find_package(CUDAToolkit REQUIRED)
     target_link_libraries(HYPRE INTERFACE CUDA::cusparse CUDA::curand CUDA::cublas)
   elseif (MFEM_USE_HIP)
-    list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_HIP:BOOL=ON)
+    list(APPEND CMAKE_OPTIONS -DHYPRE_WITH_HIP:BOOL=ON)
     find_package(rocsparse REQUIRED)
     find_package(rocrand REQUIRED)
     target_link_libraries(HYPRE INTERFACE rocsparse rocrand)
   endif()
-  if (MFEM_USE_CUDA OR MFEM_USE_HIP)
-    if (MFEM_USE_UMPIRE)
-      if (EXISTS ${umpire_DIR})
-        list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_UMPIRE:BOOL=ON -Dumpire_DIR:PATH=${umpire_DIR})
-      else()
-        message(FATAL_ERROR "MFEM_USE_UMPIRE=ON, however umpire_DIR isn't visible to HYPRE")
-      endif()
-    else()
-      list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_UMPIRE:BOOL=OFF)
-      message(WARNING
-"================================================================================
- Umpire is disabled while building HYPRE with GPU support.
- This is not recommended for performance reasons!
- Consider enabling Umpire with -DMFEM_USE_UMPIRE=ON and providing -DUMPIRE_DIR.
-================================================================================")
-    endif()
-  endif()
   if (MFEM_USE_SINGLE)
-    list(APPEND HYPRE_CMAKE_OPTIONS -DHYPRE_ENABLE_SINGLE:BOOL=ON)
+    list(APPEND CMAKE_OPTIONS -DHYPRE_ENABLE_SINGLE:BOOL=ON)
   endif()
   # define external project and create future include directory so it is present
   # to pass CMake checks at end of MFEM configuration step
-  message(STATUS "Will fetch HYPRE ${HYPRE_FETCH_TAG} to be built with ${HYPRE_CMAKE_OPTIONS}")
-  set(HYPRE_INSTALL ${CMAKE_BINARY_DIR}/fetch/hypre)
+  message(STATUS "Will fetch HYPRE ${HYPRE_FETCH_VERSION} to be built with ${CMAKE_OPTIONS}")
+  set(PREFIX ${CMAKE_BINARY_DIR}/fetch/hypre)
   include(ExternalProject)
   ExternalProject_Add(hypre
     GIT_REPOSITORY https://github.com/hypre-space/hypre.git
-    GIT_TAG ${HYPRE_FETCH_TAG}
+    GIT_TAG v${HYPRE_FETCH_VERSION}
     GIT_SHALLOW TRUE
-    GIT_PROGRESS TRUE
     UPDATE_DISCONNECTED TRUE
     SOURCE_SUBDIR src
-    PREFIX ${HYPRE_INSTALL}
-    BUILD_COMMAND ${CMAKE_COMMAND} --build . -- -j${CMAKE_BUILD_PARALLEL_LEVEL}
-    CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${HYPRE_INSTALL} -DCMAKE_INSTALL_LIBDIR:PATH=lib ${HYPRE_CMAKE_OPTIONS})
-  file(MAKE_DIRECTORY ${HYPRE_INSTALL}/include)
+    PREFIX ${PREFIX}
+    CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${PREFIX} -DCMAKE_INSTALL_LIBDIR:PATH=lib ${CMAKE_OPTIONS})
+  file(MAKE_DIRECTORY ${PREFIX}/include)
   # set imported library target properties
   add_dependencies(HYPRE hypre)
   set_target_properties(HYPRE PROPERTIES
-    IMPORTED_LOCATION ${HYPRE_INSTALL}/lib/libHYPRE.a
-    INTERFACE_INCLUDE_DIRECTORIES ${HYPRE_INSTALL}/include)
+    IMPORTED_LOCATION ${PREFIX}/lib/libHYPRE.a
+    INTERFACE_INCLUDE_DIRECTORIES ${PREFIX}/include)
   # convert HYPRE version to integer
-  if (HYPRE_FETCH_TAG MATCHES "^v?([0-9]+)\\.([0-9]+)\\.([0-9]+)$")
-    # Exact release tag X.Y.Z
-    string(REGEX MATCHALL "[0-9]+" HYPRE_SPLIT_VERSION "${HYPRE_FETCH_TAG}")
-  elseif (HYPRE_FETCH_VERSION MATCHES "([0-9]+)\\.([0-9]+)(\\.([0-9]+))?")
-    string(REGEX MATCHALL "[0-9]+" HYPRE_SPLIT_VERSION "${HYPRE_FETCH_VERSION}")
-  else (NOT DEFINED HYPRE_VERSION)
-    message(FATAL_ERROR "Unable to find HYPRE release version. Please provide it via -DHYPRE_VERSION")
-  endif()
-  if (HYPRE_SPLIT_VERSION AND NOT DEFINED HYPRE_VERSION)
-    list(GET HYPRE_SPLIT_VERSION 0 HYPRE_MAJOR_VERSION)
-    list(GET HYPRE_SPLIT_VERSION 1 HYPRE_MINOR_VERSION)
-    if (HYPRE_SPLIT_VERSION GREATER 2)
-      list(GET HYPRE_SPLIT_VERSION 2 HYPRE_PATCH_VERSION)
-    else()
-      set(HYPRE_PATCH_VERSION 0)
-    endif()
-    math(EXPR HYPRE_VERSION "10000*${HYPRE_MAJOR_VERSION} + 100*${HYPRE_MINOR_VERSION} + ${HYPRE_PATCH_VERSION}")
-    set(HYPRE_VERSION ${HYPRE_VERSION} CACHE STRING "HYPRE version." FORCE)
-  endif()
+  string(REGEX MATCHALL "[0-9]+" HYPRE_SPLIT_VERSION ${HYPRE_FETCH_VERSION})
+  list(GET HYPRE_SPLIT_VERSION 0 HYPRE_MAJOR_VERSION)
+  list(GET HYPRE_SPLIT_VERSION 1 HYPRE_MINOR_VERSION)
+  list(GET HYPRE_SPLIT_VERSION 2 HYPRE_PATCH_VERSION)
+  math(EXPR HYPRE_VERSION "10000*${HYPRE_MAJOR_VERSION} + 100*${HYPRE_MINOR_VERSION} + ${HYPRE_PATCH_VERSION}")
+  # set cache variables that would otherwise be set after mfem_find_package call
+  set(HYPRE_VERSION ${HYPRE_VERSION} CACHE STRING "HYPRE version." FORCE)
   return()
 endif()
 

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -932,14 +932,12 @@ function(mfem_export_mk_files)
   endif()
   set(MFEM_BUILD_TAG "${CMAKE_SYSTEM}")
   set(MFEM_PREFIX "${CMAKE_INSTALL_PREFIX}")
-  # For the next 4 variables, these are the values for the build-tree version of
+  # For the next 4 variable, these are the values for the build-tree version of
   # 'config.mk'
   set(MFEM_INC_DIR "${PROJECT_BINARY_DIR}")
   set(MFEM_LIB_DIR "${PROJECT_BINARY_DIR}")
   set(MFEM_TEST_MK "${PROJECT_SOURCE_DIR}/config/test.mk")
   set(MFEM_CONFIG_EXTRA "MFEM_BUILD_DIR ?= ${PROJECT_BINARY_DIR}")
-  # TODO: CUDA/HIP support:
-  set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
   set(MFEM_MPIEXEC ${MPIEXEC})
   if (NOT MFEM_MPIEXEC)
     set(MFEM_MPIEXEC "mpirun")

config/config.mk.in

@@ -88,7 +88,6 @@ MFEM_BUILD_TAG = @MFEM_BUILD_TAG@
 MFEM_PREFIX    = @MFEM_PREFIX@
 MFEM_INC_DIR   = @MFEM_INC_DIR@
 MFEM_LIB_DIR   = @MFEM_LIB_DIR@
-MFEM_XLINKER   = @MFEM_XLINKER@
 
 # Location of test.mk
 MFEM_TEST_MK = @MFEM_TEST_MK@

config/defaults.mk

@@ -57,7 +57,7 @@ CUDA_DIR = $(or $(CUDA_HOME),$(patsubst %/,%,$(dir \
 CLANG_CUDA_FLAGS = -xcuda --cuda-path=$(CUDA_DIR) --cuda-gpu-arch=$(CUDA_ARCH)
 # flags for nvcc
 NVCC_FLAGS = -x=cu --expt-extended-lambda --expt-relaxed-constexpr \
- -arch=$(CUDA_ARCH) -isystem "$(CUDA_DIR)/include"
+ -arch=$(CUDA_ARCH)
 # Prefixes for passing flags to the host compiler and linker when using
 # CUDA_CXX=nvcc
 CUDA_XCOMPILER = -Xcompiler=

data/mfem.mesh

@@ -1,593 +0,0 @@
-MFEM mesh v1.0
-
-# Created by: Pointwise
-
-# MFEM Geometry Types:
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-# PRISM       = 6
-
-dimension
-2
-
-elements
-160
-1 3 1 164 163 0
-1 3 164 165 162 163
-1 3 2 166 164 1
-1 3 166 132 165 164
-1 3 3 167 166 2
-1 3 167 131 132 166
-1 3 4 168 167 3
-1 3 168 130 131 167
-1 3 5 169 168 4
-1 3 169 129 130 168
-1 3 6 170 169 5
-1 3 170 128 129 169
-1 3 171 172 170 6
-1 3 172 127 128 170
-1 3 124 125 172 171
-1 3 125 126 127 172
-1 3 162 165 173 161
-1 3 165 132 133 173
-1 3 161 173 174 160
-1 3 173 133 134 174
-1 3 160 174 175 159
-1 3 174 134 135 175
-1 3 6 7 176 171
-1 3 7 8 177 176
-1 3 171 176 123 124
-1 3 176 177 122 123
-1 3 159 175 178 158
-1 3 175 135 136 178
-1 3 158 178 179 157
-1 3 178 136 137 179
-1 3 157 179 180 156
-1 3 179 137 138 180
-1 3 122 177 181 121
-1 3 177 8 182 181
-1 3 8 9 183 182
-1 3 9 10 184 183
-1 3 10 11 185 184
-1 3 11 12 186 185
-1 3 12 13 187 186
-1 3 13 14 15 187
-1 3 121 181 119 120
-1 3 181 182 118 119
-1 3 182 183 117 118
-1 3 183 184 188 117
-1 3 184 185 109 188
-1 3 185 186 108 109
-1 3 186 187 189 108
-1 3 187 15 16 189
-1 3 109 110 190 188
-1 3 110 111 191 190
-1 3 111 112 113 191
-1 3 188 190 116 117
-1 3 190 191 115 116
-1 3 191 113 114 115
-1 3 189 192 107 108
-1 3 192 193 106 107
-1 3 193 194 105 106
-1 3 194 195 104 105
-1 3 195 196 103 104
-1 3 16 17 192 189
-1 3 17 18 193 192
-1 3 18 19 194 193
-1 3 19 20 195 194
-1 3 20 21 196 195
-1 3 97 98 197 96
-1 3 98 99 198 197
-1 3 99 100 199 198
-1 3 100 101 200 199
-1 3 101 102 201 200
-1 3 102 103 202 201
-1 3 103 196 203 202
-1 3 196 21 22 203
-1 3 96 197 204 95
-1 3 197 198 39 204
-1 3 198 199 38 39
-1 3 199 200 205 38
-1 3 200 201 32 205
-1 3 201 202 31 32
-1 3 202 203 206 31
-1 3 203 22 23 206
-1 3 32 33 207 205
-1 3 33 34 35 207
-1 3 205 207 37 38
-1 3 207 35 36 37
-1 3 39 40 208 204
-1 3 40 41 209 208
-1 3 41 42 210 209
-1 3 42 43 211 210
-1 3 43 44 212 211
-1 3 204 208 94 95
-1 3 208 209 93 94
-1 3 209 210 92 93
-1 3 210 211 91 92
-1 3 211 212 90 91
-1 3 90 212 213 89
-1 3 212 44 214 213
-1 3 44 45 215 214
-1 3 45 46 216 215
-1 3 46 47 217 216
-1 3 47 48 218 217
-1 3 48 49 219 218
-1 3 49 50 51 219
-1 3 89 213 87 88
-1 3 213 214 86 87
-1 3 214 215 85 86
-1 3 215 216 84 85
-1 3 216 217 83 84
-1 3 217 218 82 83
-1 3 218 219 220 82
-1 3 219 51 52 220
-1 3 53 221 220 52
-1 3 221 81 82 220
-1 3 54 222 221 53
-1 3 222 80 81 221
-1 3 55 223 222 54
-1 3 223 79 80 222
-1 3 26 27 224 25
-1 3 27 28 29 224
-1 3 25 224 225 24
-1 3 224 29 30 225
-1 3 24 225 206 23
-1 3 225 30 31 206
-1 3 154 155 226 153
-1 3 155 156 180 226
-1 3 153 226 227 152
-1 3 226 180 138 227
-1 3 152 227 228 151
-1 3 227 138 139 228
-1 3 151 228 229 150
-1 3 228 139 140 229
-1 3 150 229 230 149
-1 3 229 140 141 230
-1 3 149 230 231 148
-1 3 230 141 142 231
-1 3 148 231 232 147
-1 3 231 142 143 232
-1 3 147 232 145 146
-1 3 232 143 144 145
-1 3 56 233 223 55
-1 3 233 78 79 223
-1 3 57 234 233 56
-1 3 234 77 78 233
-1 3 58 235 234 57
-1 3 235 76 77 234
-1 3 61 236 59 60
-1 3 236 235 58 59
-1 3 62 237 236 61
-1 3 237 76 235 236
-1 3 63 238 237 62
-1 3 238 75 76 237
-1 3 64 239 238 63
-1 3 239 74 75 238
-1 3 65 240 239 64
-1 3 240 73 74 239
-1 3 66 241 240 65
-1 3 241 72 73 240
-1 3 67 242 241 66
-1 3 242 71 72 241
-1 3 68 69 242 67
-1 3 69 70 71 242
-
-boundary
-164
-3 1 0 1
-3 1 1 2
-3 1 2 3
-3 1 3 4
-3 1 4 5
-3 1 5 6
-3 1 6 7
-3 1 7 8
-3 1 8 9
-3 1 9 10
-3 1 10 11
-3 1 11 12
-3 1 12 13
-3 1 13 14
-3 1 16 17
-3 1 17 18
-3 1 18 19
-3 1 19 20
-3 1 20 21
-3 1 21 22
-3 1 22 23
-3 1 23 24
-3 1 24 25
-3 1 25 26
-3 1 26 27
-3 1 27 28
-3 1 28 29
-3 1 29 30
-3 1 30 31
-3 1 31 32
-3 1 32 33
-3 1 33 34
-3 1 34 35
-3 1 35 36
-3 1 36 37
-3 1 37 38
-3 1 38 39
-3 1 39 40
-3 1 40 41
-3 1 41 42
-3 1 42 43
-3 1 43 44
-3 1 49 50
-3 1 48 49
-3 1 47 48
-3 1 46 47
-3 1 45 46
-3 1 44 45
-3 1 52 53
-3 1 53 54
-3 1 54 55
-3 1 57 58
-3 1 56 57
-3 1 55 56
-3 1 60 61
-3 1 61 62
-3 1 62 63
-3 1 63 64
-3 1 64 65
-3 1 65 66
-3 1 66 67
-3 1 67 68
-3 1 75 76
-3 1 74 75
-3 1 73 74
-3 1 72 73
-3 1 71 72
-3 1 70 71
-3 1 76 77
-3 1 77 78
-3 1 78 79
-3 1 81 82
-3 1 80 81
-3 1 79 80
-3 1 82 83
-3 1 83 84
-3 1 84 85
-3 1 85 86
-3 1 86 87
-3 1 87 88
-3 1 94 95
-3 1 93 94
-3 1 92 93
-3 1 91 92
-3 1 90 91
-3 1 96 97
-3 1 95 96
-3 1 97 98
-3 1 98 99
-3 1 99 100
-3 1 100 101
-3 1 101 102
-3 1 102 103
-3 1 107 108
-3 1 106 107
-3 1 105 106
-3 1 104 105
-3 1 103 104
-3 1 108 109
-3 1 109 110
-3 1 110 111
-3 1 111 112
-3 1 112 113
-3 1 113 114
-3 1 114 115
-3 1 115 116
-3 1 116 117
-3 1 119 120
-3 1 118 119
-3 1 117 118
-3 1 131 132
-3 1 130 131
-3 1 129 130
-3 1 128 129
-3 1 127 128
-3 1 126 127
-3 1 132 133
-3 1 133 134
-3 1 134 135
-3 1 137 138
-3 1 136 137
-3 1 135 136
-3 1 138 139
-3 1 139 140
-3 1 140 141
-3 1 141 142
-3 1 142 143
-3 1 143 144
-3 1 147 148
-3 1 146 147
-3 1 153 154
-3 1 152 153
-3 1 151 152
-3 1 150 151
-3 1 149 150
-3 1 148 149
-3 1 156 157
-3 1 157 158
-3 1 158 159
-3 1 161 162
-3 1 160 161
-3 1 159 160
-2 1 69 70
-2 1 68 69
-3 1 88 89
-3 1 89 90
-3 1 121 122
-3 1 120 121
-3 1 123 124
-3 1 122 123
-3 1 125 126
-3 1 124 125
-1 1 144 145
-1 1 145 146
-3 1 15 16
-3 1 14 15
-3 1 50 51
-3 1 51 52
-3 1 59 60
-3 1 58 59
-3 1 154 155
-3 1 155 156
-3 1 163 0
-3 1 162 163
-
-vertices
-243
-2
-4 4
-4 3.5
-4 3
-4 2.5
-4 2
-4 1.5
-4 1
-4.5 1
-5 1
-5 1.5
-5 2
-5 2.5
-5 3
-5 3.5
-5 4
-5.500 4
-6 4
-6.500 4
-7 4
-7.5 4
-8 4
-8.5 4
-9 4
-9.5 4
-10 4
-10.5 4
-11 4
-11 3.5
-11 3
-10.5 3
-10 3
-9.5 3
-9.5 2.5
-10 2.5
-10.5 2.5
-10.5 2
-10.5 1.5
-10 1.5
-9.5 1.5
-9.5 1
-10 1
-10.5 1
-11 1
-11.5 1
-12 1
-12 1.5
-12 2
-12 2.5
-12 3
-12 3.5
-12 4
-12.5 4
-13 4
-13.333 3.75
-13.666 3.5
-14.000 3.25
-14.333 3.5
-14.666 3.75
-15.000 4
-15.500 4
-16.000 4
-16.000 3.5
-16.000 3
-16.000 2.5
-16.000 2
-16.000 1.5
-16.000 1
-16.000 0.5
-16.000 0
-15.500 0
-15.000 0
-15.000 0.5000000000000002
-15.000 1
-15.000 1.5
-15.000 2
-15.000 2.5
-15.000 3
-14.666 2.75
-14.333 2.5
-14.000 2.25
-13.666 2.5
-13.333 2.75
-13 3
-13 2.5
-13 2
-13 1.5
-13 1
-13 0.500
-13 0
-12.5 0
-12 0
-11.5 0
-11 0
-10.5 0
-10 0
-9.5 0
-9 0
-8.5 0
-8.5 0.5
-8.5 1
-8.5 1.5
-8.5 2
-8.5 2.5
-8.5 3
-8 3
-7.5 3
-7 3
-6.500 3
-6 3
-6 2.5
-6.5 2.5
-7 2.5
-7.5 2.5
-7.5 2
-7.5 1.5
-7.000 1.5
-6.5 1.5
-6 1.5
-6 1
-6 0.5
-6 0
-5.5 0
-5 0
-4.5 0
-4 0
-3.5 0
-3 0
-3 0.500
-3 1
-3 1.5
-3 2
-3 2.5
-3 3
-2.666 2.75
-2.333 2.5
-2.000 2.25
-1.666 2.5
-1.333 2.75
-1.000 3
-1.000 2.5
-1.000 2
-1.000 1.5
-1.000 1
-1.000 0.5000
-1.000 0
-0.5000 0
-0.0000 0
-0.0000 0.5
-0.0000 1
-0.0000 1.5
-0.0000 2
-0.0000 2.5
-0.0000 3
-0.0000 3.5
-0.0000 4
-0.5000 4
-1.000 4
-1.333 3.75
-1.666 3.5
-2.000 3.25
-2.333 3.5
-2.666 3.75
-3 4
-3.5 4
-3.5 3.5
-3 3.5
-3.5 3
-3.5 2.5
-3.5 2
-3.5 1.5
-3.5 1
-4 0.5
-3.5 0.5
-2.666 3.25
-2.333 3
-2.000 2.75
-4.5 0.5
-5 0.5
-1.666 3
-1.333 3.25
-1.000 3.5
-5.5 0.5
-5.500 1
-5.500 1.5
-5.500 2
-5.500 2.5
-5.500 3
-5.500 3.5
-6 2
-6 3.5
-6.5 2
-7 2
-6.5 3.5
-7 3.5
-7.5 3.5
-8 3.5
-8.5 3.5
-9 0.5
-9 1
-9 1.5
-9 2
-9 2.5
-9 3
-9 3.5
-9.5 0.5
-9.5 2
-9.5 3.5
-10 2
-10 0.5
-10.5 0.5
-11 0.5
-11.5 0.5
-12 0.5
-12.5 0.500
-12.5 1
-12.5 1.5
-12.5 2
-12.5 2.5
-12.5 3
-12.5 3.5
-13 3.5
-13.333 3.250
-13.666 3
-14.000 2.75
-10.5 3.5
-10 3.5
-0.500 3.5
-0.500 3
-0.500 2.5
-0.500 2
-0.500 1.5
-0.500 1
-0.500 0.5
-14.333 3
-14.666 3.25
-15.000 3.5
-15.500 3.5
-15.500 3
-15.500 2.5
-15.500 2
-15.500 1.5
-15.500 1
-15.500 0.5

doc/CodeDocumentation.dox

@@ -202,7 +202,6 @@ namespace mfem {
  * - <a class="el" href="tesla_8cpp_source.html">Tesla</a>: simple magnetostatics simulation code
  * - <a class="el" href="maxwell_8cpp_source.html">Maxwell</a>: simple transient full-wave electromagnetics simulation code
  * - <a class="el" href="joule_8cpp_source.html">Joule</a>: transient magnetics and Joule heating miniapp
- * - <a class="el" href="lorentz_8cpp_source.html">Lorentz</a>: simple particle tracking code based on the Lorentz force
  * - <a class="el" href="classmfem_1_1navier_1_1NavierSolver.html">Navier</a>: solve the transient incompressible Navier-Stokes equations
  * - <a class="el" href="mobius-strip_8cpp_source.html">Mobius Strip</a>: generate various Mobius strip-like meshes
  * - <a class="el" href="klein-bottle_8cpp_source.html">Klein Bottle</a>: generate three types of Klein bottle surfaces

fem/bilinearform.cpp

@@ -99,6 +99,7 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
    boundary_integs_marker = bf->boundary_integs_marker;
 
    interior_face_integs = bf->interior_face_integs;
+   interior_face_integs_marker = bf->interior_face_integs_marker;
 
    boundary_face_integs = bf->boundary_face_integs;
    boundary_face_integs_marker = bf->boundary_face_integs_marker;
@@ -254,9 +255,17 @@ void BilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
    boundary_integs_marker.Append(&bdr_marker);
 }
 
+void BilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi,
+                                             Array<int> *face_marker)
+{
+   interior_face_integs.Append (bfi);
+   interior_face_integs_marker.Append(face_marker);
+}
+
 void BilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi)
 {
    interior_face_integs.Append (bfi);
+   interior_face_integs_marker.Append(nullptr);
 }
 
 void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
@@ -685,6 +694,11 @@ void BilinearForm::Assemble(int skip_zeros)
             vdofs.Append (vdofs2);
             for (int k = 0; k < interior_face_integs.Size(); k++)
             {
+               // Skip the face if it's not in the integrator's attribute list.
+               if (interior_face_integs_marker[k] &&
+                   interior_face_integs_marker[k]->Find(tr->Attribute) == -1)
+               { continue; }
+
                interior_face_integs[k]->
                AssembleFaceMatrix(*fes->GetFE(tr->Elem1No),
                                   *fes->GetFE(tr->Elem2No),
@@ -1502,6 +1516,12 @@ void MixedBilinearForm::AddTraceFaceIntegrator (BilinearFormIntegrator * bfi)
    trace_face_integs.Append (bfi);
 }
 
+
+void MixedBilinearForm::AddFaceIntegrator (BilinearFormIntegrator *bfi)
+{
+   face_integs.Append(bfi);
+}
+
 void MixedBilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi)
 {
    boundary_trace_face_integs.Append(bfi);
@@ -1670,7 +1690,7 @@ void MixedBilinearForm::Assemble(int skip_zeros)
    if (boundary_face_integs.Size())
    {
       FaceElementTransformations *ftr;
-      Array<int> tr_vdofs2, te_vdofs2;
+      Array<int> trial_vdofs2, test_vdofs2;
       const FiniteElement *trial_fe1, *trial_fe2, *test_fe1, *test_fe2;
 
       // Which boundary attributes need to be processed?
@@ -1761,10 +1781,39 @@ void MixedBilinearForm::Assemble(int skip_zeros)
       }
    }
 
+   if (face_integs.Size())
+   {
+      FaceElementTransformations *ftr;
+      Array<int> trial_vdofs2, test_vdofs2;
+
+      int nfaces = mesh->GetNumFaces();
+      for (int i = 0; i < nfaces; i++)
+      {
+         ftr = mesh->GetFaceElementTransformations(i);
+         trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
+         test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
+         if (ftr->Elem2No >= 0)
+         {
+            trial_fes->GetElementVDofs(ftr->Elem2No, trial_vdofs2);
+            test_fes->GetElementVDofs(ftr->Elem2No, test_vdofs2);
+            trial_vdofs.Append(trial_vdofs2);
+            test_vdofs.Append(test_vdofs2);
+         }
+
+         for (int k = 0; k < face_integs.Size(); k++)
+         {
+            face_integs[k]->AssembleFaceMatrix(*trial_fes->GetFE(ftr->Elem1No),
+                                               *test_fes->GetFE(ftr->Elem1No),
+                                               *ftr, elemmat);
+            mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
+         }
+      }
+   }
+
    if (boundary_trace_face_integs.Size())
    {
       FaceElementTransformations *ftr;
-      Array<int> te_vdofs2;
+      Array<int> test_vdofs2;
       const FiniteElement *trial_face_fe, *test_fe1, *test_fe2;
 
       // Which boundary attributes need to be processed?
@@ -2355,6 +2404,8 @@ MixedBilinearForm::~MixedBilinearForm()
       { delete boundary_face_integs[i]; }
       for (i = 0; i < trace_face_integs.Size(); i++)
       { delete trace_face_integs[i]; }
+      for (i = 0; i < face_integs.Size(); i++)
+      { delete face_integs[i]; }
       for (i = 0; i < boundary_trace_face_integs.Size(); i++)
       { delete boundary_trace_face_integs[i]; }
    }

fem/bilinearform.hpp

@@ -114,6 +114,10 @@ protected:
 
    /// Set of interior face Integrators to be applied.
    Array<BilinearFormIntegrator*> interior_face_integs;
+   /// List of attributes for each integrator. The integrator is applie only on
+   /// faces that have such attributes. Corresponds to Mesh::GetFaceAttribute().
+   /// Note: it is a list; it's not a marker over all existing face attributes.
+   Array<Array<int>*> interior_face_integs_marker;  ///< Entries are not owned.
 
    /// Set of boundary face Integrators to be applied.
    Array<BilinearFormIntegrator*> boundary_face_integs;
@@ -278,6 +282,7 @@ public:
 
    /// Access all integrators added with AddInteriorFaceIntegrator().
    Array<BilinearFormIntegrator*> *GetFBFI() { return &interior_face_integs; }
+   Array<Array<int>*> *GetFBFI_Marker() { return &interior_face_integs_marker; }
 
    /// Access all integrators added with AddBdrFaceIntegrator().
    Array<BilinearFormIntegrator*> *GetBFBFI() { return &boundary_face_integs; }
@@ -424,6 +429,9 @@ public:
                               Array<int> &bdr_marker);
 
    /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
+   void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi,
+                                  Array<int> *face_marker);
+
    void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi);
 
    /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
@@ -796,6 +804,9 @@ protected:
    /// Trace face (skeleton) integrators.
    Array<BilinearFormIntegrator*> trace_face_integs;
 
+   /// face integrators.
+   Array<BilinearFormIntegrator*> face_integs;
+
    /// Boundary trace face (skeleton) integrators.
    Array<BilinearFormIntegrator*> boundary_trace_face_integs;
    /// Entries are not owned.
@@ -910,6 +921,9 @@ public:
    /// Adds a boundary integrator. Assumes ownership of @a bfi.
    void AddBoundaryIntegrator(BilinearFormIntegrator *bfi);
 
+   /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
+   void AddFaceIntegrator(BilinearFormIntegrator *bfi);
+
    /// Adds a boundary integrator. Assumes ownership of @a bfi.
    void AddBoundaryIntegrator(BilinearFormIntegrator * bfi,
                               Array<int> &bdr_marker);

fem/fe_coll.cpp

@@ -2456,7 +2456,8 @@ RT_FECollection::RT_FECollection(const int order, const int dim,
       const char *cb_name = BasisType::Name(cb_type); // this may abort
       MFEM_ABORT("unknown closed BasisType: " << cb_name);
    }
-   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid)
+   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid &&
+       ob_type != BasisType::IntegratedGLL)
    {
       const char *ob_name = BasisType::Name(ob_type); // this may abort
       MFEM_ABORT("unknown open BasisType: " << ob_name);
@@ -2783,7 +2784,8 @@ ND_FECollection::ND_FECollection(const int p, const int dim,
    int cp_type = BasisType::GetQuadrature1D(cb_type);
 
    // Error checking
-   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid)
+   if (Quadrature1D::CheckOpen(op_type) == Quadrature1D::Invalid &&
+       ob_type != BasisType::IntegratedGLL)
    {
       const char *ob_name = BasisType::Name(ob_type);
       MFEM_ABORT("Invalid open basis point type: " << ob_name);

fem/fespace.hpp

@@ -676,6 +676,20 @@ public:
    /// Copy assignment not supported
    FiniteElementSpace& operator=(const FiniteElementSpace&) = delete;
 
+   void ReplaceElemDofTable(const Table &new_elem_dof, int ndofs_new)
+   {
+      *elem_dof = new_elem_dof;
+      ndofs = ndofs_new;
+   }
+   void ReplaceBdrElemDofTable(const Table &new_bdrElem_dof)
+   {
+      *bdr_elem_dof = new_bdrElem_dof;
+   }
+   void ReplaceFaceDofTable(const Table &new_face_dof)
+   {
+      *face_dof = new_face_dof;
+   }
+
    /// Returns the mesh
    inline Mesh *GetMesh() const { return mesh; }
 

fem/intrules.cpp

@@ -947,7 +947,6 @@ int Quadrature1D::CheckOpen(int type)
       case OpenUniform:
       case ClosedUniform:
       case OpenHalfUniform:
-      case ClosedGL:
          return type; // all types can work as open
       default:
          return Invalid;

fem/linearform.cpp

@@ -105,6 +105,13 @@ void LinearForm::AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
    interior_face_integs.Append(lfi);
 }
 
+void LinearForm::AddTraceFaceIntegrator(LinearFormIntegrator *tfi,
+                                        Array<int> &attr_list)
+{
+   trace_face_integs.Append(tfi);
+   trace_face_integs_attributes.Append(&attr_list);
+}
+
 bool LinearForm::SupportsDevice() const
 {
    // return false for NURBS meshes, so we don’t convert it to non-NURBS
@@ -313,6 +320,39 @@ void LinearForm::Assemble()
       }
    }
 
+   if (trace_face_integs.Size())
+   {
+      Mesh *mesh = fes->GetMesh();
+      FaceElementTransformations *tr;
+      const FiniteElement *fe_1, *fe_2;
+      Array<int> vdofs2;
+
+      const int nfaces = mesh->GetNumFaces();
+      for (int f = 0; f < nfaces; f++)
+      {
+         const int attr = mesh->GetFace(f)->GetAttribute();
+         tr   = mesh->GetFaceElementTransformations(f);
+         fe_1 = fes->GetFE(tr->Elem1No);
+         fes->GetElementVDofs(tr->Elem1No, vdofs);
+         if (tr->Elem2No >= 0)
+         {
+            fes->GetElementVDofs(tr->Elem2No, vdofs2);
+            vdofs.Append(vdofs2);
+            fe_2 = fes->GetFE(tr->Elem2No);
+         }
+         else
+         {
+            fe_2 = fe_1;
+         }
+
+         for (int k = 0; k < trace_face_integs.Size(); k++)
+         {
+            trace_face_integs[k]->AssembleRHSElementVect(*fe_1, *fe_2, *tr, elemvect);
+            AddElementVector(vdofs, elemvect);
+         }
+      }
+   }
+
    if (interior_face_integs.Size())
    {
       Mesh *mesh = fes->GetMesh();

fem/linearform.hpp

@@ -65,6 +65,11 @@ protected:
    /// Set of Internal Face Integrators to be applied.
    Array<LinearFormIntegrator*> interior_face_integs;
 
+   /// Set of trace (all faces - both interior and boundary) integrators.
+   Array<LinearFormIntegrator *> trace_face_integs;
+   Array<Array<int> *>
+   trace_face_integs_attributes; ///< Entries are not owned.
+
    /// The element ids where the centers of the delta functions lie
    Array<int> domain_delta_integs_elem_id;
 
@@ -151,6 +156,13 @@ public:
    /// Adds new Boundary Face Integrator. Assumes ownership of @a lfi.
    void AddBdrFaceIntegrator(LinearFormIntegrator *lfi);
 
+   /** @brief Add new Trace Face Integrator, restricted to the given face
+       attributes.
+       Assumes ownership of @a tfi. The array @a attr_list is stored
+       internally as a pointer to the given Array<int> object. */
+   void AddTraceFaceIntegrator(LinearFormIntegrator *tfi,
+                               Array<int> &attr_list);
+
    /** @brief Add new Boundary Face Integrator, restricted to the given boundary
        attributes.
 
@@ -182,6 +194,9 @@ public:
    /// Access all integrators added with AddBdrFaceIntegrator().
    Array<LinearFormIntegrator*> *GetFLFI() { return &boundary_face_integs; }
 
+   /// Access all integrators added with AddTraceFaceIntegrator().
+   Array<LinearFormIntegrator*> *GetTLFI() { return &trace_face_integs; }
+
    /// Access all integrators added with AddInteriorFaceIntegrator().
    Array<LinearFormIntegrator*> *GetIFLFI() { return &interior_face_integs; }
 

fem/lor/lor.cpp

@@ -61,6 +61,51 @@ void LORBase::AddIntegratorsAndMarkers(BilinearForm &a_from,
    }
 }
 
+void LORBase::AddIntegratorsAndMarkers2(BilinearForm &a_from,
+                                       BilinearForm &a_to,
+                                       GetIntegratorsFn get_integrators,
+                                       GetMarkersFn get_markers,
+                                       AddIntegratorMarkersFn add_integrator_marker,
+                                       AddIntegratorFn add_integrator,
+                                       const IntegrationRule *ir)
+{
+   Array<BilinearFormIntegrator*> *integrators = (a_from.*get_integrators)();
+   Array<Array<int>*> &markers = *(a_from.*get_markers)();
+
+   for (int i=0; i<integrators->Size(); ++i)
+   {
+      BilinearFormIntegrator *integrator = (*integrators)[i];
+      if (markers[i] != nullptr)
+      {
+         (a_to.*add_integrator_marker)(integrator, *markers[i]);
+      }
+      else
+      {
+         (a_to.*add_integrator)(integrator);
+      }
+      ir_map[integrator] = integrator->GetIntRule();
+      if (ir) { integrator->SetIntegrationRule(*ir); }
+   }
+}
+
+void LORBase::AddIntegratorsAndMarkers(BilinearForm &a_from,
+                                       BilinearForm &a_to,
+                                       GetIntegratorsFn get_integrators,
+                                       GetMarkersFn get_markers,
+                                       AddIntegratorMarkersPtrFn add_integrator,
+                                       const IntegrationRule *ir)
+{
+   Array<BilinearFormIntegrator*> *integrators = (a_from.*get_integrators)();
+   Array<Array<int>*> *markers = (a_from.*get_markers)();
+
+   for (int i=0; i<integrators->Size(); ++i)
+   {
+      (a_to.*add_integrator)((*integrators)[i], *markers[i]);
+      ir_map[(*integrators)[i]] = ((*integrators)[i])->GetIntegrationRule();
+      if (ir) { ((*integrators)[i])->SetIntegrationRule(*ir); }
+   }
+}
+
 void LORBase::ResetIntegrationRules(GetIntegratorsFn get_integrators)
 {
    Array<BilinearFormIntegrator*> *integrators = (a->*get_integrators)();

fem/lor/lor.hpp

@@ -27,6 +27,8 @@ private:
    using AddIntegratorFn = void (BilinearForm::*)(BilinearFormIntegrator*);
    using AddIntegratorMarkersFn =
       void (BilinearForm::*)(BilinearFormIntegrator*, Array<int>&);
+   using AddIntegratorMarkersPtrFn =
+      void (BilinearForm::*)(BilinearFormIntegrator*, Array<int>*);
 
    IntegrationRules irs;
    const IntegrationRule *ir_el, *ir_face;
@@ -52,6 +54,19 @@ private:
                                  AddIntegratorMarkersFn add_integrator_marker,
                                  AddIntegratorFn add_integrator,
                                  const IntegrationRule *ir);
+   void AddIntegratorsAndMarkers2(BilinearForm &a_from,
+                                 BilinearForm &a_to,
+                                 GetIntegratorsFn get_integrators,
+                                 GetMarkersFn get_markers,
+                                 AddIntegratorMarkersFn add_integrator_marker,
+                                 AddIntegratorFn add_integrator,
+                                 const IntegrationRule *ir);
+   void AddIntegratorsAndMarkers(BilinearForm &a_from,
+                                 BilinearForm &a_to,
+                                 GetIntegratorsFn get_integrators,
+                                 GetMarkersFn get_markers,
+                                 AddIntegratorMarkersPtrFn add_integrator,
+                                 const IntegrationRule *ir);
 
    /// Resets the integration rules of the integrators of @a a to their original
    /// values (after temporarily changing them for LOR assembly).

fem/nonlinearform.cpp

@@ -436,7 +436,7 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
    // In parallel, the result is in 'py' which is an alias for 'aux2'.
 }
 
-Operator &NonlinearForm::GetGradient(const Vector &x, bool finalize) const
+Operator &NonlinearForm::GetGradient(const Vector &x) const
 {
    if (ext)
    {
@@ -644,8 +644,6 @@ Operator &NonlinearForm::GetGradient(const Vector &x, bool finalize) const
       }
    }
 
-   if (!finalize) { return *Grad; }
-
    if (!Grad->Finalized())
    {
       Grad->Finalize(skip_zeros);
@@ -1205,14 +1203,7 @@ const BlockVector &BlockNonlinearForm::Prolongate(const BlockVector &bx) const
       aux1.Update(block_offsets);
       for (int s = 0; s < fes.Size(); s++)
       {
-         if (P[s])
-         {
-            P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
-         }
-         else
-         {
-            aux1.GetBlock(s) = bx.GetBlock(s);
-         }
+         P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
       }
       return aux1;
    }
@@ -1241,16 +1232,11 @@ void BlockNonlinearForm::Mult(const Vector &x, Vector &y) const
       {
          cP[s]->MultTranspose(pby.GetBlock(s), by.GetBlock(s));
       }
-      else if (needs_prolongation)
-      {
-         by.GetBlock(s) = pby.GetBlock(s);
-      }
       by.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
    }
 }
 
-void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx,
-                                                bool finalize) const
+void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
 {
    const int skip_zeros = 0;
    Array<Array<int> *> vdofs(fes.Size());
@@ -1504,7 +1490,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx,
       }
    }
 
-   if (finalize && !Grads(0,0)->Finalized())
+   if (!Grads(0,0)->Finalized())
    {
       for (int i=0; i<fes.Size(); ++i)
       {
@@ -1543,23 +1529,7 @@ Operator &BlockNonlinearForm::GetGradient(const Vector &x) const
          for (int s2 = 0; s2 < fes.Size(); ++s2)
          {
             delete cGrads(s1, s2);
-            if (cP[s1] && cP[s2])
-            {
-               cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
-            }
-            else if (cP[s1])
-            {
-               cGrads(s1, s2) = TransposeMult(*cP[s1], *Grads(s1, s2));
-            }
-            else if (cP[s2])
-            {
-               cGrads(s1, s2) = mfem::Mult(*Grads(s1, s2), *cP[s2]);
-            }
-            else
-            {
-               cGrads(s1, s2) = NULL;
-               continue;
-            }
+            cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
             mGrads(s1, s2) = cGrads(s1, s2);
          }
       }

fem/nonlinearform.hpp

@@ -217,12 +217,7 @@ public:
        In general, @a x may have non-homogeneous essential boundary values.
 
        The state @a x must be a true-dof vector. */
-   Operator &GetGradient(const Vector &x) const override { return GetGradient(x, true); }
-
-   /** @brief Compute the gradient Operator of the NonlinearForm corresponding
-       to the state @a x with optional finalization and elimintaion. */
-   /** @see GetGradient(const Vector &) */
-   Operator &GetGradient(const Vector &x, bool finalize) const;
+   Operator &GetGradient(const Vector &x) const override;
 
    /// Update the NonlinearForm to propagate updates of the associated FE space.
    /** After calling this method, the essential boundary conditions need to be
@@ -313,7 +308,7 @@ protected:
    void MultBlocked(const BlockVector &bx, BlockVector &by) const;
 
    /// Specialized version of GetGradient() for BlockVector
-   void ComputeGradientBlocked(const BlockVector &bx, bool finalize = true) const;
+   void ComputeGradientBlocked(const BlockVector &bx) const;
 
 public:
    /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().

fem/pbilinearform.cpp

@@ -151,15 +151,6 @@ void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
    }
 }
 
-HypreParMatrix *ParBilinearForm::ParallelAssembleInternalMatrix()
-{
-   if (p_mat.Ptr() == NULL)
-   {
-      ParallelAssemble(p_mat, mat);
-   }
-   return p_mat.As<HypreParMatrix>();
-}
-
 void ParBilinearForm::ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local)
 {
    A.Clear();
@@ -342,15 +333,6 @@ void ParBilinearForm
    A.EliminateRowsCols(dof_list, X, B);
 }
 
-void ParBilinearForm::ParallelEliminateEssentialBC(
-   const Array<int> &bdr_attr_is_ess, const HypreParVector &X, HypreParVector &B)
-{
-   Array<int> dof_list;
-   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
-
-   p_mat.As<HypreParMatrix>()->EliminateRowsCols(dof_list, X, B);
-}
-
 HypreParMatrix *ParBilinearForm::
 ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                              HypreParMatrix &A) const
@@ -362,26 +344,6 @@ ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
    return A.EliminateRowsCols(dof_list);
 }
 
-void ParBilinearForm::ParallelEliminateEssentialBC(const Array<int>
-                                                   &bdr_attr_is_ess)
-{
-   Array<int> tdofs_list;
-   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, tdofs_list);
-
-   ParallelEliminateTDofs(tdofs_list);
-}
-
-void ParBilinearForm::ParallelEliminateTDofs(const Array<int> &tdofs_list)
-{
-   p_mat_e.EliminateRowsCols(p_mat, tdofs_list);
-}
-
-void ParBilinearForm::ParallelEliminateTDofsInRHS(
-   const Array<int> &tdofs_list, const Vector &x, Vector &b)
-{
-   p_mat.EliminateBC(p_mat_e, tdofs_list, x, b);
-}
-
 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const real_t a)
 const
 {
@@ -523,7 +485,7 @@ void ParBilinearForm::FormLinearSystem(
       HypreParVector true_X(pfes), true_B(pfes);
       P.MultTranspose(b, true_B);
       R.Mult(x, true_X);
-      ParallelEliminateTDofsInRHS(ess_tdof_list, true_X, true_B);
+      p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
       R.MultTranspose(true_B, b);
       hybridization->ReduceRHS(true_B, B);
       X.SetSize(B.Size());
@@ -536,11 +498,17 @@ void ParBilinearForm::FormLinearSystem(
       B.SetSize(X.Size());
       P.MultTranspose(b, B);
       R.Mult(x, X);
-      ParallelEliminateTDofsInRHS(ess_tdof_list, X, B);
+      p_mat.EliminateBC(p_mat_e, ess_tdof_list, X, B);
       if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
    }
 }
 
+void ParBilinearForm::EliminateVDofsInRHS(
+   const Array<int> &vdofs, const Vector &x, Vector &b)
+{
+   p_mat.EliminateBC(p_mat_e, vdofs, x, b);
+}
+
 void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
                                        OperatorHandle &A)
 {
@@ -585,7 +553,7 @@ void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
          mat = NULL;
          delete mat_e;
          mat_e = NULL;
-         ParallelEliminateTDofs(ess_tdof_list);
+         p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
       }
       if (hybridization)
       {
@@ -647,180 +615,36 @@ void ParBilinearForm::Update(FiniteElementSpace *nfes)
    p_mat_e.Clear();
 }
 
-void ParMixedBilinearForm::pAllocMat()
-{
-   const int trial_nbr_size = trial_pfes->GetFaceNbrVSize();
-   const int test_nbr_size = test_pfes->GetFaceNbrVSize();
 
-   if (keep_nbr_block)
-   {
-      mat = new SparseMatrix(height + test_nbr_size, width + trial_nbr_size);
-   }
-   else
-   {
-      mat = new SparseMatrix(height, width + trial_nbr_size);
-   }
+HypreParMatrix *ParMixedBilinearForm::ParallelAssemble()
+{
+   // construct the block-diagonal matrix A
+   HypreParMatrix *A =
+      new HypreParMatrix(trial_pfes->GetComm(),
+                         test_pfes->GlobalVSize(),
+                         trial_pfes->GlobalVSize(),
+                         test_pfes->GetDofOffsets(),
+                         trial_pfes->GetDofOffsets(),
+                         mat);
+
+   HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
+                             trial_pfes->Dof_TrueDof_Matrix());
+
+   delete A;
+
+   return rap;
 }
 
-void ParMixedBilinearForm::AssembleSharedFaces(int skip_zeros)
+void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A)
 {
-   ParMesh *pmesh = trial_pfes->GetParMesh();
-   FaceElementTransformations *T;
-   Array<int> tr_vdofs1, tr_vdofs2, tr_vdofs_all;
-   Array<int> te_vdofs1, te_vdofs2, te_vdofs_all;
-   DenseMatrix elemmat;
-
-   int nfaces = pmesh->GetNSharedFaces();
-   for (int i = 0; i < nfaces; i++)
-   {
-      T = pmesh->GetSharedFaceTransformations(i);
-      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
-      trial_pfes->GetElementVDofs(T->Elem1No, tr_vdofs1);
-      test_pfes->GetElementVDofs(T->Elem1No, te_vdofs1);
-      trial_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, tr_vdofs2);
-      test_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, te_vdofs2);
-
-      tr_vdofs1.Copy(tr_vdofs_all);
-      for (int j = 0; j < tr_vdofs2.Size(); j++)
-      {
-         if (tr_vdofs2[j] >= 0)
-         {
-            tr_vdofs2[j] += width;
-         }
-         else
-         {
-            tr_vdofs2[j] -= width;
-         }
-      }
-      tr_vdofs_all.Append(tr_vdofs2);
-
-      if (keep_nbr_block)
-      {
-         te_vdofs1.Copy(te_vdofs_all);
-         for (int j = 0; j < te_vdofs2.Size(); j++)
-         {
-            if (te_vdofs2[j] >= 0)
-            {
-               te_vdofs2[j] += height;
-            }
-            else
-            {
-               te_vdofs2[j] -= height;
-            }
-         }
-         te_vdofs_all.Append(te_vdofs2);
-      }
-
-      for (int k = 0; k < interior_face_integs.Size(); k++)
-      {
-         interior_face_integs[k]->
-         AssembleFaceMatrix(*trial_pfes->GetFE(T->Elem1No),
-                            *test_pfes->GetFE(T->Elem1No),
-                            *trial_pfes->GetFaceNbrFE(Elem2NbrNo),
-                            *test_pfes->GetFaceNbrFE(Elem2NbrNo),
-                            *T, elemmat);
-         if (keep_nbr_block)
-         {
-            mat->AddSubMatrix(te_vdofs_all, tr_vdofs_all, elemmat, skip_zeros);
-         }
-         else
-         {
-            mat->AddSubMatrix(te_vdofs1, tr_vdofs_all, elemmat, skip_zeros);
-         }
-      }
-   }
-}
-
-void ParMixedBilinearForm::Assemble(int skip_zeros)
-{
-   if (interior_face_integs.Size())
-   {
-      trial_pfes->ExchangeFaceNbrData();
-      test_pfes->ExchangeFaceNbrData();
-      if (!ext && mat == NULL)
-      {
-         pAllocMat();
-      }
-   }
-
-   MixedBilinearForm::Assemble(skip_zeros);
-
-   if (!ext && interior_face_integs.Size() > 0)
-   {
-      AssembleSharedFaces(skip_zeros);
-   }
-}
-
-HypreParMatrix *ParMixedBilinearForm::ParallelAssembleInternalMatrix()
-{
-   if (p_mat.Ptr() == NULL)
-   {
-      ParallelAssemble(p_mat, mat);
-   }
-   return p_mat.As<HypreParMatrix>();
-}
-
-HypreParMatrix *ParMixedBilinearForm::ParallelAssemble(SparseMatrix *m)
-{
-   OperatorHandle Mh(Operator::Hypre_ParCSR);
-   ParallelAssemble(Mh, m);
-   Mh.SetOperatorOwner(false);
-   return Mh.As<HypreParMatrix>();
-}
-
-void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A,
-                                            SparseMatrix *A_local)
-{
-   A.Clear();
-
-   if (A_local == NULL) { return; }
-   MFEM_VERIFY(A_local->Finalized(), "the local matrix must be finalized");
-
-   OperatorHandle dA(A.Type()), hdA;
-
-   if (interior_face_integs.Size() == 0)
-   {
-      // construct the rectangular block-diagonal matrix dA
-      dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
-                                  test_pfes->GlobalVSize(),
-                                  trial_pfes->GlobalVSize(),
-                                  test_pfes->GetDofOffsets(),
-                                  trial_pfes->GetDofOffsets(),
-                                  A_local);
-   }
-   else
-   {
-      // handle the case when 'a' contains off-diagonal
-      const int lvrows = test_pfes->GetVSize();
-      const int lvcols = trial_pfes->GetVSize();
-      const HYPRE_BigInt *face_nbr_glob_lcol = trial_pfes->GetFaceNbrGlobalDofMap();
-      const HYPRE_BigInt lcol_offset = trial_pfes->GetMyDofOffset();
-
-      Array<HYPRE_BigInt> glob_J(A_local->NumNonZeroElems());
-      const int *J = A_local->GetJ();
-      for (int i = 0; i < glob_J.Size(); i++)
-      {
-         if (J[i] < lvcols)
-         {
-            glob_J[i] = J[i] + lcol_offset;
-         }
-         else
-         {
-            glob_J[i] = face_nbr_glob_lcol[J[i] - lvcols];
-         }
-      }
-
-      // TODO - construct dA directly in the A format
-      hdA.Reset(
-         new HypreParMatrix(trial_pfes->GetComm(), lvrows, test_pfes->GlobalVSize(),
-                            trial_pfes->GlobalVSize(), A_local->GetI(), glob_J,
-                            A_local->GetData(), test_pfes->GetDofOffsets(),
-                            trial_pfes->GetDofOffsets()));
-      // - hdA owns the new HypreParMatrix
-      // - the above constructor copies all input arrays
-      glob_J.DeleteAll();
-      dA.ConvertFrom(hdA);
-   }
+   // construct the rectangular block-diagonal matrix dA
+   OperatorHandle dA(A.Type());
+   dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
+                               test_pfes->GlobalVSize(),
+                               trial_pfes->GlobalVSize(),
+                               test_pfes->GetDofOffsets(),
+                               trial_pfes->GetDofOffsets(),
+                               mat);
 
    OperatorHandle P_test(A.Type()), P_trial(A.Type());
 
@@ -846,44 +670,6 @@ void ParMixedBilinearForm::TrueAddMult(const Vector &x, Vector &y,
    test_pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Yaux, 1.0, y);
 }
 
-void ParMixedBilinearForm::ParallelEliminateTrialEssentialBC(
-   const Array<int> &bdr_attr_is_ess)
-{
-   Array<int> trial_tdof_list;
-   trial_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, trial_tdof_list);
-
-   ParallelEliminateTrialTDofs(trial_tdof_list);
-}
-
-void ParMixedBilinearForm::ParallelEliminateTrialTDofs(
-   const Array<int> &trial_tdof_list)
-{
-   HypreParMatrix *temp = p_mat.As<HypreParMatrix>()->EliminateCols(
-                             trial_tdof_list);
-   p_mat_e.Reset(temp, true);
-}
-
-void ParMixedBilinearForm::ParallelEliminateTrialTDofsInRHS(
-   const Array<int> &trial_tdof_list, const Vector &x, Vector &b)
-{
-   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, x, 1.0, b);
-}
-
-void ParMixedBilinearForm::ParallelEliminateTestEssentialBC(
-   const Array<int> &bdr_attr_is_ess)
-{
-   Array<int> test_tdof_list;
-   test_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, test_tdof_list);
-
-   ParallelEliminateTestTDofs(test_tdof_list);
-}
-
-void ParMixedBilinearForm::ParallelEliminateTestTDofs(
-   const Array<int> &test_tdof_list)
-{
-   p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
-}
-
 void ParMixedBilinearForm::FormRectangularSystemMatrix(
    const Array<int>
    &trial_tdof_list,
@@ -904,8 +690,10 @@ void ParMixedBilinearForm::FormRectangularSystemMatrix(
       mat = NULL;
       delete mat_e;
       mat_e = NULL;
-      ParallelEliminateTrialTDofs(trial_tdof_list);
-      ParallelEliminateTestTDofs(test_tdof_list);
+      HypreParMatrix *temp =
+         p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
+      p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
+      p_mat_e.Reset(temp, true);
    }
 
    A = p_mat;
@@ -935,7 +723,7 @@ void ParMixedBilinearForm::FormRectangularLinearSystem(
    test_P->MultTranspose(b, B);
    trial_R->Mult(x, X);
 
-   ParallelEliminateTrialTDofsInRHS(trial_tdof_list, X, B);
+   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, B);
    B.SetSubVector(test_tdof_list, 0.0);
 }
 

fem/pbilinearform.hpp

@@ -73,7 +73,7 @@ public:
    /** When set to true and the ParBilinearForm has interior face integrators,
        the local SparseMatrix will include the rows (in addition to the columns)
        corresponding to face-neighbor dofs. The default behavior is to disregard
-       those rows. Must be called before the first Assemble() call. */
+       those rows. Must be called before the first Assemble call. */
    void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
 
    /** @brief Set the operator type id for the parallel matrix/operator when
@@ -101,14 +101,6 @@ public:
        diagonal for this case. */
    void AssembleDiagonal(Vector &diag) const override;
 
-   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
-   /** The returned matrix is the internal one, owned by the form. It is not
-       reassembled if it has been already constructed. If FormSystemMatrix()
-       has been called before, it is the system matrix with eliminated
-       essential DOFs, otherwise the parallel matrix is assembled here without
-       the elimination process. */
-   HypreParMatrix *ParallelAssembleInternalMatrix();
-
    /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
    /** The returned matrix has to be deleted by the caller. */
    HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
@@ -154,13 +146,6 @@ public:
                                      const HypreParVector &X,
                                      HypreParVector &B) const;
 
-   /// Eliminate essential boundary DOFs from the parallel system matrix.
-   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary. */
-   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
-                                     const HypreParVector &X,
-                                     HypreParVector &B);
-
    /// Eliminate essential boundary DOFs from a parallel assembled matrix @a A.
    /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
        the essential part of the boundary. The eliminated part is stored in a
@@ -172,12 +157,6 @@ public:
    HypreParMatrix *ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                                                 HypreParMatrix &A) const;
 
-   /// Eliminate essential boundary DOFs from the parallel system matrix.
-   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary. This method relies on
-       ParallelEliminateTDofs(const Array<int> &), see it for details. */
-   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess);
-
    /// Eliminate essential true DOFs from a parallel assembled matrix @a A.
    /** Given a list of essential true dofs and the parallel assembled matrix
        @a A, eliminate the true dofs from the matrix, storing the eliminated
@@ -190,28 +169,6 @@ public:
                                           HypreParMatrix &A) const
    { return A.EliminateRowsCols(tdofs_list); }
 
-   /// Eliminate essential true DOFs from the parallel system matrix.
-   /** Given a list of essential true dofs, eliminate the true dofs from
-       the parallel assembled system matrix, storing the eliminated part
-       internally. This method works in conjunction with
-       ParallelEliminateTDofsInRHS() and allows elimination of boundary
-       conditions in multiple right-hand sides. */
-   void ParallelEliminateTDofs(const Array<int> &tdofs_list);
-
-   /** @brief Use the stored eliminated part of the parallel system matrix for
-       elimination of boundary conditions in the r.h.s. */
-   /** Given a list of essential true dofs, eliminate the true dofs from the
-       right-hand side @a b using the solution vector @a x and the previously
-       stored eliminated part of the parallel assembled system matrix produced
-       by ParallelEliminateTDofs(const Array<int> &). */
-   void ParallelEliminateTDofsInRHS(const Array<int> &tdofs, const Vector &x,
-                                    Vector &b);
-
-   /// @deprecated Use ParallelEliminateTDofsInRHS() instead.
-   MFEM_DEPRECATED void EliminateVDofsInRHS(const Array<int> &vdofs,
-                                            const Vector &x, Vector &b)
-   { ParallelEliminateTDofsInRHS(vdofs, x, b); }
-
    /** @brief Compute @a y += @a a (P^t A P) @a x, where @a x and @a y are
        vectors on the true dofs. */
    void TrueAddMult(const Vector &x, Vector &y, const real_t a = 1.0) const;
@@ -281,6 +238,8 @@ public:
 
    void Update(FiniteElementSpace *nfes = NULL) override;
 
+   void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
+
    virtual ~ParBilinearForm() { }
 };
 
@@ -298,13 +257,6 @@ protected:
    /// Matrix and eliminated matrix
    OperatorHandle p_mat, p_mat_e;
 
-   bool keep_nbr_block;
-
-   // Allocate mat - called when (mat == NULL && fbfi.Size() > 0)
-   void pAllocMat();
-
-   void AssembleSharedFaces(int skip_zeros = 1);
-
 private:
    /// Copy construction is not supported; body is undefined.
    ParMixedBilinearForm(const ParMixedBilinearForm &);
@@ -324,7 +276,6 @@ public:
    {
       trial_pfes = trial_fes;
       test_pfes  = test_fes;
-      keep_nbr_block = false;
    }
 
    /** @brief Create a ParMixedBilinearForm on the given FiniteElementSpace%s
@@ -344,89 +295,15 @@ public:
    {
       trial_pfes = trial_fes;
       test_pfes  = test_fes;
-      keep_nbr_block = false;
    }
 
-   /** When set to true and the ParMixedBilinearForm has interior face
-       integrators, the local SparseMatrix will include the rows (in addition
-       to the columns) corresponding to face-neighbor dofs. The default
-       behavior is to disregard those rows. Must be called before the first
-       Assemble() call. */
-   void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
-
-   /// Assemble the local matrix
-   void Assemble(int skip_zeros = 1);
-
    /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
-   /** The returned matrix is the internal one, owned by the form. It is not
-       reassembled if it has been already constructed. If
-       FormRectangularSystemMatrix() has been called before, it is the system
-       matrix with eliminated essential DOFs, otherwise the parallel matrix is
-       assembled here without the elimination process. */
-   HypreParMatrix *ParallelAssembleInternalMatrix();
-
-   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
-   /** The returned matrix has to be deleted by the caller. */
-   HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
-
-   /** @brief Returns the eliminated matrix assembled on the true dofs, i.e.
-       P_test^t A_local P_trial. */
-   /** The returned matrix has to be deleted by the caller. */
-   HypreParMatrix *ParallelAssembleElim() { return ParallelAssemble(mat_e); }
-
-   /** @brief Return the matrix @a m assembled on the true dofs, i.e. P_test^t
-       A_local P_trial. */
-   /** The returned matrix has to be deleted by the caller. */
-   HypreParMatrix *ParallelAssemble(SparseMatrix *m);
+   HypreParMatrix *ParallelAssemble();
 
    /** @brief Returns the matrix assembled on the true dofs, i.e.
        @a A = P_test^t A_local P_trial, in the format (type id) specified by
        @a A. */
-   void ParallelAssemble(OperatorHandle &A) { ParallelAssemble(A, mat); }
-
-   /** Returns the eliminated matrix assembled on the true dofs, i.e.
-       @a A_elim = P^t A_elim_local P in the format (type id) specified by @a A.
-    */
-   void ParallelAssembleElim(OperatorHandle &A_elim)
-   { ParallelAssemble(A_elim, mat_e); }
-
-   /** Returns the matrix @a A_local assembled on the true dofs, i.e.
-       @a A = P_test^t A_local P_trial in the format (type id) specified by
-       @a A. */
-   void ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local);
-
-   /// Eliminate essential boundary trial DOFs from the parallel system matrix.
-   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary. This method relies on
-       ParallelEliminateTrialTDofs(const Array<int> &), see it for details. */
-   void ParallelEliminateTrialEssentialBC(const Array<int> &bdr_attr_is_ess);
-
-   /// Eliminate essential trial true DOFs from the parallel system matrix.
-   /** Given a list of essential trial true dofs, eliminate the trial true dofs
-       from the parallel assembled system matrix, storing the eliminated part
-       internally. This method works in conjunction with
-       ParallelEliminateTrialTDofsInRHS() and allows elimination of boundary
-       conditions in multiple right-hand sides. */
-   void ParallelEliminateTrialTDofs(const Array<int> &trial_tdof_list);
-
-   /** @brief Use the stored eliminated part of the parallel system matrix for
-       elimination of boundary conditions in the r.h.s. */
-   /** Given a list of essential trial true dofs, eliminate the trial true dofs
-       from the right-hand side @a B using the solution vector @a X and the
-       previously stored eliminated part of the parallel assembled system
-       matrix produced by ParallelEliminateTrialTDofs(const Array<int> &). */
-   void ParallelEliminateTrialTDofsInRHS(const Array<int> &trial_tdof_list,
-                                         const Vector &X, Vector &B);
-
-   /// Eliminate essential boundary test DOFs from the parallel system matrix.
-   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
-       the essential part of the boundary. */
-   void ParallelEliminateTestEssentialBC(const Array<int> &bdr_attr_is_ess);
-
-   /// Eliminate essential test true DOFs from the parallel system matrix.
-   /** Given a list of essential test true dofs, eliminate the test true dofs
-       from the parallel assembled system matrix. */
-   void ParallelEliminateTestTDofs(const Array<int> &test_tdof_list);
+   void ParallelAssemble(OperatorHandle &A);
 
    using MixedBilinearForm::FormRectangularSystemMatrix;
    using MixedBilinearForm::FormRectangularLinearSystem;

fem/pnonlinearform.cpp

@@ -105,59 +105,6 @@ const SparseMatrix &ParNonlinearForm::GetLocalGradient(const Vector &x) const
    return *Grad;
 }
 
-void ParNonlinearForm::GradientSharedFaces(const Vector &x,
-                                           int skip_zeros) const
-{
-   ParFiniteElementSpace *pfes = ParFESpace();
-   ParMesh *pmesh = pfes->GetParMesh();
-   FaceElementTransformations *T;
-   Array<int> vdofs1, vdofs2, vdofs_all;
-   DenseMatrix elemmat;
-   Vector el_x, nbr_x, face_x;
-   const Vector &px = Prolongate(x);
-
-   ParGridFunction pgf(pfes, const_cast<Vector&>(px), 0);
-   pgf.ExchangeFaceNbrData();
-
-   int nfaces = pmesh->GetNSharedFaces();
-   for (int i = 0; i < nfaces; i++)
-   {
-      T = pmesh->GetSharedFaceTransformations(i);
-      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
-
-      pfes->GetElementVDofs(T->Elem1No, vdofs1);
-      pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
-      face_x.SetSize(vdofs1.Size() + vdofs2.Size());
-
-      el_x.MakeRef(face_x, 0, vdofs1.Size());
-      pgf.GetSubVector(vdofs1, el_x);
-
-      nbr_x.MakeRef(face_x, vdofs1.Size(), vdofs2.Size());
-      pgf.FaceNbrData().GetSubVector(vdofs2, nbr_x);
-
-      vdofs1.Copy(vdofs_all);
-      for (int j = 0; j < vdofs2.Size(); j++)
-      {
-         if (vdofs2[j] >= 0)
-         {
-            vdofs2[j] += height;
-         }
-         else
-         {
-            vdofs2[j] -= height;
-         }
-      }
-      vdofs_all.Append(vdofs2);
-      for (int k = 0; k < fnfi.Size(); k++)
-      {
-         fnfi[k]->AssembleFaceGrad(*pfes->GetFE(T->Elem1No),
-                                   *pfes->GetFaceNbrFE(Elem2NbrNo),
-                                   *T, face_x, elemmat);
-         Grad->AddSubMatrix(vdofs1, vdofs_all, elemmat, skip_zeros);
-      }
-   }
-}
-
 Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 {
    if (NonlinearForm::ext) { return NonlinearForm::GetGradient(x); }
@@ -165,61 +112,19 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
    ParFiniteElementSpace *pfes = ParFESpace();
 
    pGrad.Clear();
-   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type()), hdA;
 
-   if (fnfi.Size())
+   NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
+
+   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
+
+   if (fnfi.Size() == 0)
    {
-      const int skip_zeros = 0;
-
-      pfes->ExchangeFaceNbrData();
-      if (Grad == NULL)
-      {
-         int nbr_size = pfes->GetFaceNbrVSize();
-         Grad = new SparseMatrix(pfes->GetVSize(), pfes->GetVSize() + nbr_size);
-      }
-
-      NonlinearForm::GetGradient(x, false); // (re)assemble Grad, no b.c.
-
-      GradientSharedFaces(x, skip_zeros);
-
-      Grad->Finalize(skip_zeros);
-
-      // handle the case when 'a' contains off-diagonal
-      int lvsize = pfes->GetVSize();
-      const HYPRE_BigInt *face_nbr_glob_ldof = pfes->GetFaceNbrGlobalDofMap();
-      HYPRE_BigInt ldof_offset = pfes->GetMyDofOffset();
-
-      Array<HYPRE_BigInt> glob_J(Grad->NumNonZeroElems());
-      int *J = Grad->GetJ();
-      for (int i = 0; i < glob_J.Size(); i++)
-      {
-         if (J[i] < lvsize)
-         {
-            glob_J[i] = J[i] + ldof_offset;
-         }
-         else
-         {
-            glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
-         }
-      }
-
-      // TODO - construct dA directly in the A format
-      hdA.Reset(
-         new HypreParMatrix(pfes->GetComm(), lvsize, pfes->GlobalVSize(),
-                            pfes->GlobalVSize(), Grad->GetI(), glob_J,
-                            Grad->GetData(), pfes->GetDofOffsets(),
-                            pfes->GetDofOffsets()));
-      // - hdA owns the new HypreParMatrix
-      // - the above constructor copies all input arrays
-      glob_J.DeleteAll();
-      dA.ConvertFrom(hdA);
+      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
+                             pfes->GetDofOffsets(), Grad);
    }
    else
    {
-      NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
-
-      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
-                             pfes->GetDofOffsets(), Grad);
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
    }
 
    // RAP the local gradient dA.
@@ -366,70 +271,7 @@ void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
 
    if (fnfi.Size() > 0)
    {
-      // Terms over shared interior faces in parallel.
-      ParMesh *pmesh = ParFESpace(0)->GetParMesh();
-      FaceElementTransformations *tr;
-
-      Array<Array<int> *>vdofs(fes.Size());
-      Array<Array<int> *>vdofs2(fes.Size());
-      Array<Vector *> el_x(fes.Size());
-      Array<const Vector *> el_x_const(fes.Size());
-      Array<Vector *> el_y(fes.Size());
-      Array<const FiniteElement *> fe(fes.Size());
-      Array<const FiniteElement *> fe2(fes.Size());
-      Array<ParGridFunction *> pgfs(fes.Size());
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         el_x_const[s] = el_x[s] = new Vector();
-         el_y[s] = new Vector();
-         vdofs[s] = new Array<int>;
-         vdofs2[s] = new Array<int>;
-         pgfs[s] = new ParGridFunction(const_cast<ParFiniteElementSpace*>(ParFESpace(s)),
-                                       xs.GetBlock(s));
-         pgfs[s]->ExchangeFaceNbrData();
-      }
-
-      const int n_shared_faces = pmesh->GetNSharedFaces();
-      for (int i = 0; i < n_shared_faces; i++)
-      {
-         tr = pmesh->GetSharedFaceTransformations(i, true);
-         int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
-
-         for (int s=0; s<fes.Size(); ++s)
-         {
-            const ParFiniteElementSpace *pfes = ParFESpace(s);
-            fe[s] = pfes->GetFE(tr->Elem1No);
-            fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
-
-            pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
-            pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
-
-            el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
-            xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
-            pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
-                                                el_x[s]->GetData() + vdofs[s]->Size());
-         }
-
-         for (int k = 0; k < fnfi.Size(); ++k)
-         {
-            fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
-
-            for (int s=0; s<fes.Size(); ++s)
-            {
-               if (el_y[s]->Size() == 0) { continue; }
-               ys.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
-            }
-         }
-      }
-
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         delete pgfs[s];
-         delete vdofs2[s];
-         delete vdofs[s];
-         delete el_y[s];
-         delete el_x[s];
-      }
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
    }
 
    for (int s=0; s<fes.Size(); ++s)
@@ -486,106 +328,6 @@ void ParBlockNonlinearForm::SetGradientType(Operator::Type tid)
    }
 }
 
-void ParBlockNonlinearForm::GradientSharedFaces(const BlockVector &xs,
-                                                int skip_zeros) const
-{
-   // Terms over shared interior faces in parallel.
-   ParMesh *pmesh = ParFESpace(0)->GetParMesh();
-   FaceElementTransformations *tr;
-
-   Array<Array<int> *>vdofs(fes.Size());
-   Array<Array<int> *>vdofs2(fes.Size());
-   Array<Array<int> *>vdofs_all(fes.Size());
-   Array<Vector *> el_x(fes.Size());
-   Array<const Vector *> el_x_const(fes.Size());
-   Array2D<DenseMatrix *> elmats(fes.Size(), fes.Size());
-   Array<const FiniteElement *> fe(fes.Size());
-   Array<const FiniteElement *> fe2(fes.Size());
-   Array<ParGridFunction *> pgfs(fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      el_x_const[s1] = el_x[s1] = new Vector();
-      vdofs[s1] = new Array<int>;
-      vdofs2[s1] = new Array<int>;
-      vdofs_all[s1] = new Array<int>;
-      pgfs[s1] = new ParGridFunction(
-         const_cast<ParFiniteElementSpace*>(ParFESpace(s1)),
-         const_cast<Vector&>(xs.GetBlock(s1)));
-      pgfs[s1]->ExchangeFaceNbrData();
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         elmats(s1,s2) = new DenseMatrix();
-      }
-   }
-
-   const int n_shared_faces = pmesh->GetNSharedFaces();
-   for (int i = 0; i < n_shared_faces; i++)
-   {
-      tr = pmesh->GetSharedFaceTransformations(i, true);
-      int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
-
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         const ParFiniteElementSpace *pfes = ParFESpace(s);
-         fe[s] = pfes->GetFE(tr->Elem1No);
-         fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
-
-         pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
-         pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
-
-         el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
-         xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
-         pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
-                                             el_x[s]->GetData() + vdofs[s]->Size());
-
-         vdofs[s]->Copy(*vdofs_all[s]);
-
-         const int lvsize = pfes->GetVSize();
-         for (int j = 0; j < vdofs2[s]->Size(); j++)
-         {
-            if ((*vdofs2[s])[j] >= 0)
-            {
-               (*vdofs2[s])[j] += lvsize;
-            }
-            else
-            {
-               (*vdofs2[s])[j] -= lvsize;
-            }
-         }
-         vdofs_all[s]->Append(*(vdofs2[s]));
-      }
-
-      for (int k = 0; k < fnfi.Size(); ++k)
-      {
-         fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
-
-         for (int s1=0; s1<fes.Size(); ++s1)
-         {
-            for (int s2=0; s2<fes.Size(); ++s2)
-            {
-               if (elmats(s1,s2)->Height() == 0) { continue; }
-               Grads(s1,s2)->AddSubMatrix(*vdofs[s1], *vdofs_all[s2],
-                                          *elmats(s1,s2), skip_zeros);
-            }
-         }
-      }
-   }
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      delete pgfs[s1];
-      delete vdofs_all[s1];
-      delete vdofs2[s1];
-      delete vdofs[s1];
-      delete el_x[s1];
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete elmats(s1,s2);
-      }
-   }
-}
-
 BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
 {
    if (pBlockGrad == NULL)
@@ -605,155 +347,49 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
       }
    }
 
-   // xs_true is not modified, so const_cast is okay
-   xs_true.Update(const_cast<Vector &>(x), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
+   GetLocalGradient(x); // gradients are stored in 'Grads'
 
    if (fnfi.Size() > 0)
    {
-      const int skip_zeros = 0;
-
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         const_cast<ParFiniteElementSpace*>(pfes[s])->ExchangeFaceNbrData();
-      }
-
-      for (int s1=0; s1<fes.Size(); ++s1)
-      {
-         for (int s2=0; s2<fes.Size(); ++s2)
-         {
-            if (Grads(s1,s2) == NULL)
-            {
-               int nbr_size = pfes[s2]->GetFaceNbrVSize();
-               Grads(s1,s2) = new SparseMatrix(pfes[s1]->GetVSize(),
-                                               pfes[s2]->GetVSize() + nbr_size);
-            }
-         }
-      }
-
-      // (re)assemble Grad without b.c. into 'Grads'
-      BlockNonlinearForm::ComputeGradientBlocked(xs, false);
-
-      GradientSharedFaces(xs, skip_zeros);
-
-      // finalize the gradients
-      for (int s1=0; s1<fes.Size(); ++s1)
-         for (int s2=0; s2<fes.Size(); ++s2)
-         {
-            Grads(s1,s2)->Finalize(skip_zeros);
-         }
-
-      for (int s1=0; s1<fes.Size(); ++s1)
-      {
-         for (int s2=0; s2<fes.Size(); ++s2)
-         {
-            OperatorHandle hdA;
-            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                           Ph(phBlockGrad(s1,s2)->Type()),
-                           Rh(phBlockGrad(s1,s2)->Type());
-
-            // handle the case when 'a' contains off-diagonal
-            int lvsize = pfes[s2]->GetVSize();
-            const HYPRE_BigInt *face_nbr_glob_ldof =
-               const_cast<ParFiniteElementSpace*>(pfes[s2])->GetFaceNbrGlobalDofMap();
-            HYPRE_BigInt ldof_offset = pfes[s2]->GetMyDofOffset();
-
-            Array<HYPRE_BigInt> glob_J(Grads(s1,s2)->NumNonZeroElems());
-            int *J = Grads(s1,s2)->GetJ();
-            for (int i = 0; i < glob_J.Size(); i++)
-            {
-               if (J[i] < lvsize)
-               {
-                  glob_J[i] = J[i] + ldof_offset;
-               }
-               else
-               {
-                  glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
-               }
-            }
-
-            // TODO - construct dA directly in the A format
-            hdA.Reset(
-               new HypreParMatrix(pfes[s2]->GetComm(), pfes[s1]->GetVSize(),
-                                  pfes[s1]->GlobalVSize(), pfes[s2]->GlobalVSize(),
-                                  Grads(s1,s2)->GetI(), glob_J, Grads(s1,s2)->GetData(),
-                                  pfes[s1]->GetDofOffsets(), pfes[s2]->GetDofOffsets()));
-            // - hdA owns the new HypreParMatrix
-            // - the above constructor copies all input arrays
-            glob_J.DeleteAll();
-            dA.ConvertFrom(hdA);
-
-            if (s1 == s2)
-            {
-               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-               OperatorHandle Ae;
-               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
-            }
-            else
-            {
-               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
-            }
-
-            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
-         }
-      }
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
    }
-   else
+
+   for (int s1=0; s1<fes.Size(); ++s1)
    {
-      // (re)assemble Grad without b.c. into 'Grads'
-      BlockNonlinearForm::ComputeGradientBlocked(xs);
-
-      for (int s1=0; s1<fes.Size(); ++s1)
+      for (int s2=0; s2<fes.Size(); ++s2)
       {
-         for (int s2=0; s2<fes.Size(); ++s2)
+         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                        Ph(phBlockGrad(s1,s2)->Type()),
+                        Rh(phBlockGrad(s1,s2)->Type());
+
+         if (s1 == s2)
          {
-            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                           Ph(phBlockGrad(s1,s2)->Type()),
-                           Rh(phBlockGrad(s1,s2)->Type());
+            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
+                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
+            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
 
-            if (s1 == s2)
-            {
-               dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                      pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-               OperatorHandle Ae;
-               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
-            }
-            else
-            {
-               dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                           pfes[s1]->GlobalVSize(),
-                                           pfes[s2]->GlobalVSize(),
-                                           pfes[s1]->GetDofOffsets(),
-                                           pfes[s2]->GetDofOffsets(),
-                                           Grads(s1,s2));
-               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
-            }
-
-            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+            OperatorHandle Ae;
+            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
          }
+         else
+         {
+            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
+                                        pfes[s1]->GlobalVSize(),
+                                        pfes[s2]->GlobalVSize(),
+                                        pfes[s1]->GetDofOffsets(),
+                                        pfes[s2]->GetDofOffsets(),
+                                        Grads(s1,s2));
+            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+         }
+
+         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
       }
    }
 

fem/pnonlinearform.hpp

@@ -29,8 +29,6 @@ protected:
    mutable ParGridFunction X, Y;
    mutable OperatorHandle pGrad;
 
-   void GradientSharedFaces(const Vector &x, int skip_zeros = 1) const;
-
 public:
    ParNonlinearForm(ParFiniteElementSpace *pf);
 
@@ -83,8 +81,6 @@ protected:
    mutable Array2D<OperatorHandle *> phBlockGrad;
    mutable BlockOperator *pBlockGrad;
 
-   void GradientSharedFaces(const BlockVector &xs, int skip_zeros) const;
-
 public:
    /// Computes the energy of the system
    real_t GetEnergy(const Vector &x) const override;

general/array.cpp

@@ -0,0 +1,214 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+
+// Abstract array data type
+
+#include "array.hpp"
+#include "../general/forall.hpp"
+#include <fstream>
+#include <type_traits>
+
+namespace mfem
+{
+
+template <class T>
+void Array<T>::Print(std::ostream &os, int width) const
+{
+   for (int i = 0; i < size; i++)
+   {
+      os << data[i];
+      if ( !((i+1) % width) || i+1 == size )
+      {
+         os << '\n';
+      }
+      else
+      {
+         os << " ";
+      }
+   }
+}
+
+template <class T>
+void Array<T>::Save(std::ostream &os, int fmt) const
+{
+   if (fmt == 0)
+   {
+      os << size << '\n';
+   }
+   for (int i = 0; i < size; i++)
+   {
+      os << operator[](i) << '\n';
+   }
+}
+
+template <class T>
+void Array<T>::Load(std::istream &in, int fmt)
+{
+   if (fmt == 0)
+   {
+      int new_size;
+      in >> new_size;
+      SetSize(new_size);
+   }
+   for (int i = 0; i < size; i++)
+   {
+      in >> operator[](i);
+   }
+}
+
+template <class T>
+T Array<T>::Max() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T max = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (max < operator[](i))
+      {
+         max = operator[](i);
+      }
+   }
+
+   return max;
+}
+
+template <class T>
+T Array<T>::Min() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T min = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (operator[](i) < min)
+      {
+         min = operator[](i);
+      }
+   }
+
+   return min;
+}
+
+// Partial Sum
+template <class T>
+void Array<T>::PartialSum()
+{
+   T sum = static_cast<T>(0);
+   for (int i = 0; i < size; i++)
+   {
+      sum+=operator[](i);
+      operator[](i) = sum;
+   }
+}
+
+template <class T>
+void Array<T>::Abs()
+{
+   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
+   const bool useDevice = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(useDevice);
+   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = std::abs(y[i]);
+   });
+}
+
+// Sum
+template <class T>
+T Array<T>::Sum() const
+{
+   T sum = static_cast<T>(0);
+   for (int i = 0; i < size; i++)
+   {
+      sum+=operator[](i);
+   }
+
+   return sum;
+}
+
+template <class T>
+int Array<T>::IsSorted() const
+{
+   T val_prev = operator[](0), val;
+   for (int i = 1; i < size; i++)
+   {
+      val=operator[](i);
+      if (val < val_prev)
+      {
+         return 0;
+      }
+      val_prev = val;
+   }
+
+   return 1;
+}
+
+template <class T>
+bool Array<T>::IsConstant() const
+{
+   if (size < 2) { return true; }
+   const T v0 = data[0];
+   for (int i = 1; i < size; i++)
+   {
+      if (data[i] != v0)
+      {
+         return false;
+      }
+   }
+
+   return true;
+}
+
+template <class T>
+void Array2D<T>::Load(const char *filename, int fmt)
+{
+   std::ifstream in;
+   in.open(filename, std::ifstream::in);
+   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
+   Load(in, fmt);
+   in.close();
+}
+
+template <class T>
+void Array2D<T>::Print(std::ostream &os, int width_)
+{
+   int height = this->NumRows();
+   int width  = this->NumCols();
+
+   for (int i = 0; i < height; i++)
+   {
+      os << "[row " << i << "]\n";
+      for (int j = 0; j < width; j++)
+      {
+         os << (*this)(i,j);
+         if ( (j+1) == width_ || (j+1) % width_ == 0 )
+         {
+            os << '\n';
+         }
+         else
+         {
+            os << ' ';
+         }
+      }
+   }
+}
+
+template class Array<char>;
+template class Array<int>;
+template class Array<long long>;
+template class Array<real_t>;
+template class Array2D<int>;
+template class Array2D<real_t>;
+
+} // namespace mfem

general/array.hpp

@@ -16,13 +16,9 @@
 #include "mem_manager.hpp"
 #include "device.hpp"
 #include "error.hpp"
-#include "forall.hpp"
 #include "globals.hpp"
-#include "reducers.hpp"
-#include "scan.hpp"
 
 #include <iostream>
-#include <fstream>
 #include <cstdlib>
 #include <cstring>
 #include <algorithm>
@@ -139,8 +135,6 @@ public:
    /// Return the device flag of the Memory object used by the Array
    bool UseDevice() const { return data.UseDevice(); }
 
-   void UseDevice(bool use_dev) { data.UseDevice(use_dev); }
-
    /// Return true if the data will be deleted by the Array
    inline bool OwnsData() const { return data.OwnsHostPtr(); }
 
@@ -281,11 +275,11 @@ public:
 
    /** @brief Find the maximal element in the array, using the comparison
        operator `<` for class T. */
-   inline T Max() const;
+   T Max() const;
 
    /** @brief Find the minimal element in the array, using the comparison
        operator `<` for class T. */
-   inline T Min() const;
+   T Min() const;
 
    /// Sorts the array in ascending order. This requires operator< to be defined for T.
    void Sort() { std::sort((T*)data, data + size); }
@@ -303,22 +297,22 @@ public:
    }
 
    /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
-   inline int IsSorted() const;
+   int IsSorted() const;
 
    /// Does the Array have Size zero.
    bool IsEmpty() const { return Size() == 0; }
 
    /// Return true if all entries of the array are the same.
-   inline bool IsConstant() const;
+   bool IsConstant() const;
 
    /// Fill the entries of the array with the cumulative sum of the entries.
-   inline void PartialSum();
+   void PartialSum();
 
    /// Replace each entry of the array with its absolute value.
-   inline void Abs();
+   void Abs();
 
    /// Return the sum of all the array entries using the '+'' operator for class 'T'.
-   inline T Sum() const;
+   T Sum() const;
 
    /// Set all entries of the array to the provided constant.
    inline void operator=(const T &a);
@@ -803,14 +797,8 @@ template <typename T> template <typename CT>
 inline Array<T> &Array<T>::operator=(const Array<CT> &src)
 {
    SetSize(src.Size());
-
-   const bool use_dev = UseDevice() || src.UseDevice();
-   const auto x = src.Read(use_dev);
-   auto y = Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = x[i];
-   });
+   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
+   return *this;
 }
 
 template <class T>
@@ -1026,24 +1014,19 @@ template <class T>
 inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
 {
    sa.SetSize(sa_size);
-   const bool use_dev = UseDevice() || sa.UseDevice();
-   const auto x = Read(use_dev);
-   auto y = sa.Write(use_dev);
-   mfem::forall_switch(use_dev, sa_size, [=] MFEM_HOST_DEVICE (int i)
+   for (int i = 0; i < sa_size; i++)
    {
-      y[i] = x[offset + i];
-   });
+      sa[i] = (*this)[offset+i];
+   }
 }
 
 template <class T>
 inline void Array<T>::operator=(const T &a)
 {
-   const bool use_dev = UseDevice();
-   auto x = Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   for (int i = 0; i < size; i++)
    {
-      x[i] = a;
-   });
+      data[i] = a;
+   }
 }
 
 template <class T>
@@ -1052,153 +1035,6 @@ inline void Array<T>::Assign(const T *p)
    data.CopyFromHost(p, Size());
 }
 
-template <class T>
-inline void Array<T>::Print(std::ostream &os, int width) const
-{
-   for (int i = 0; i < size; i++)
-   {
-      os << data[i];
-      if ( !((i+1) % width) || i+1 == size )
-      {
-         os << '\n';
-      }
-      else
-      {
-         os << " ";
-      }
-   }
-}
-
-template <class T>
-inline void Array<T>::Save(std::ostream &os, int fmt) const
-{
-   if (fmt == 0)
-   {
-      os << size << '\n';
-   }
-   for (int i = 0; i < size; i++)
-   {
-      os << operator[](i) << '\n';
-   }
-}
-
-template <class T>
-void Array<T>::Load(std::istream &in, int fmt)
-{
-   if (fmt == 0)
-   {
-      int new_size;
-      in >> new_size;
-      SetSize(new_size);
-   }
-   for (int i = 0; i < size; i++)
-   {
-      in >> operator[](i);
-   }
-}
-
-template <class T>
-inline T Array<T>::Max() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T max = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (max < operator[](i))
-      {
-         max = operator[](i);
-      }
-   }
-
-   return max;
-}
-
-template <class T>
-inline T Array<T>::Min() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T min = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (operator[](i) < min)
-      {
-         min = operator[](i);
-      }
-   }
-
-   return min;
-}
-
-// Partial Sum
-template <class T>
-inline void Array<T>::PartialSum()
-{
-   auto data_ptr = ReadWrite(UseDevice());
-   InclusiveScan(UseDevice(), data_ptr, data_ptr, size);
-}
-
-template <class T>
-inline void Array<T>::Abs()
-{
-   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
-   const bool useDevice = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(useDevice);
-   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = std::abs(y[i]);
-   });
-}
-
-// Sum
-template <class T>
-inline T Array<T>::Sum() const
-{
-   T sum = static_cast<T>(0);
-   if (size > 0)
-   {
-      const auto m_data = Read(UseDevice());
-      reduce(size, sum, [=] MFEM_HOST_DEVICE(int i, T &r) { r += m_data[i]; },
-      /*  */ SumReducer<T> {}, UseDevice());
-   }
-   return sum;
-}
-
-template <class T>
-inline int Array<T>::IsSorted() const
-{
-   T val_prev = operator[](0), val;
-   for (int i = 1; i < size; i++)
-   {
-      val=operator[](i);
-      if (val < val_prev)
-      {
-         return 0;
-      }
-      val_prev = val;
-   }
-
-   return 1;
-}
-
-template <class T>
-inline bool Array<T>::IsConstant() const
-{
-   if (size < 2) { return true; }
-   const T v0 = data[0];
-   for (int i = 1; i < size; i++)
-   {
-      if (data[i] != v0)
-      {
-         return false;
-      }
-   }
-
-   return true;
-}
-
 
 template <class T>
 inline const T &Array2D<T>::operator()(int i, int j) const
@@ -1238,40 +1074,6 @@ inline T *Array2D<T>::operator[](int i)
    return &array1d[i*N];
 }
 
-template <class T>
-void Array2D<T>::Load(const char *filename, int fmt)
-{
-   std::ifstream in;
-   in.open(filename, std::ifstream::in);
-   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
-   Load(in, fmt);
-   in.close();
-}
-
-template <class T>
-void Array2D<T>::Print(std::ostream &os, int width_)
-{
-   int height = this->NumRows();
-   int width  = this->NumCols();
-
-   for (int i = 0; i < height; i++)
-   {
-      os << "[row " << i << "]\n";
-      for (int j = 0; j < width; j++)
-      {
-         os << (*this)(i,j);
-         if ( (j+1) == width_ || (j+1) % width_ == 0 )
-         {
-            os << '\n';
-         }
-         else
-         {
-            os << ' ';
-         }
-      }
-   }
-}
-
 
 template <class T>
 inline void Swap(Array2D<T> &a, Array2D<T> &b)

general/reducers.hpp

@@ -12,6 +12,7 @@
 #ifndef MFEM_REDUCERS_HPP
 #define MFEM_REDUCERS_HPP
 
+#include "array.hpp"
 #include "forall.hpp"
 
 #include <cmath>
@@ -513,33 +514,6 @@ template<class B, class R> struct reduction_kernel
       }
    }
 };
-
-template <class T>
-class ReductionWorkspace
-{
-   Memory<T> workspace;
-
-   static ReductionWorkspace &Instance()
-   {
-      static ReductionWorkspace instance;
-      return instance;
-   }
-
-   ~ReductionWorkspace() { workspace.Delete(); }
-
-public:
-   static T *Get(int num_blocks)
-   {
-      ReductionWorkspace &instance = Instance();
-      if (instance.workspace.Capacity() < num_blocks)
-      {
-         instance.workspace.Delete();
-         instance.workspace.New(num_blocks, MemoryType::HOST_PINNED);
-      }
-      return instance.workspace;
-   }
-};
-
 }
 
 /**
@@ -555,7 +529,8 @@ public:
  @tparam T value_type to operate on
  */
 template <class T, class B, class R>
-void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
+void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
+            Array<T> &workspace)
 {
    if (N == 0)
    {
@@ -592,7 +567,13 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
 
       red_type red{nullptr, std::forward<B>(body), reducer, N, items_per_thread};
       // allocate res to fit block_size entries
-      auto work = internal::ReductionWorkspace<T>::Get(nblocks);
+      auto mt = workspace.GetMemory().GetMemoryType();
+      if (mt != MemoryType::HOST_PINNED && mt != MemoryType::MANAGED)
+      {
+         mt = MemoryType::HOST_PINNED;
+      }
+      workspace.SetSize(nblocks, mt);
+      auto work = workspace.HostWrite();
       red.work = work;
       forall_2D(nblocks, block_size, 1, std::move(red));
       // wait for results

general/scan.hpp

@@ -28,37 +28,8 @@
 
 namespace mfem
 {
-
-namespace internal
-{
-class ScanWorkspace
-{
-   Memory<std::byte> workspace;
-   static ScanWorkspace &Instance()
-   {
-      static ScanWorkspace instance;
-      return instance;
-   }
-   ~ScanWorkspace() { workspace.Delete(); }
-public:
-   static std::byte *Get(int num_bytes)
-   {
-      ScanWorkspace &instance = Instance();
-      if (Size() < num_bytes)
-      {
-         instance.workspace.Delete();
-         instance.workspace.New(num_bytes);
-      }
-      return instance.workspace.Write(MemoryClass::DEVICE, Size());
-   }
-   static int Size()
-   {
-      return Instance().workspace.Capacity();
-   }
-};
-}
-
-/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, std::plus<>{})
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
+/// std::plus<>{})
 template <class InputIt, class OutputIt>
 void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 {
@@ -66,12 +37,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
-      if (bytes > 0)
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items);
+                       workspace.Write(), bytes, d_in, d_out, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -86,12 +57,11 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
       }
       // try allocating a larger buffer
       bytes = 0;
-      // get size of buffer
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
                         nullptr, bytes, d_in, d_out, num_items));
-      // resize buffer (in ScanWorkspace::Get) and try again
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items));
+                        workspace.Write(), bytes, d_in, d_out, num_items));
       return;
    }
 #endif
@@ -131,13 +101,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
-      if (bytes > 0)
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                       num_items);
+                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -154,9 +123,9 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, num_items));
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                        num_items));
+                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
       return;
    }
 #endif
@@ -195,13 +164,13 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
    if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
       if (bytes)
       {
          auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                       init_value, num_items);
+                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                       num_items);
 #if defined(MFEM_USE_CUDA)
          if (err == cudaSuccess)
          {
@@ -218,9 +187,10 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
       bytes = 0;
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
                         nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
+      workspace.SetSize(bytes);
       MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                        init_value, num_items));
+                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                        num_items));
       return;
    }
 #endif
@@ -243,7 +213,7 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 }
 
 /// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
-/// std::plus<>{})
+/// workspace, std::plus<>{})
 template <class InputIt, class OutputIt, class T>
 void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
                    T init_value)

general/table.cpp

@@ -193,6 +193,22 @@ void Table::ShiftUpI()
    I[0] = 0;
 }
 
+void Table::ReplaceConnection(int r, int c_old, int c_new)
+{
+   for (int j=I[r]; j<I[r+1]; j++)
+   {
+      if ( J[j] == c_old) { J[j] = c_new; }
+   }
+}
+
+void Table::RemoveRow(int r)
+{
+   for (int j=I[r]; j<I[r+1]; j++)
+   {
+      J[j] = -1;
+   }
+}
+
 void Table::SetSize(int dim, int connections_per_row)
 {
    SetDims (dim, dim * connections_per_row);

general/table.hpp

@@ -89,6 +89,9 @@ public:
    void AddConnections (int r, const int *c, int nc);
    void ShiftUpI();
 
+   void ReplaceConnection(int r, int c_old, int c_new);
+   void RemoveRow(int r);
+
    /// Set the size and the number of connections for the table.
    void SetSize(int dim, int connections_per_row);
 

linalg/batched/magma.cpp

@@ -167,7 +167,7 @@ void MagmaBatchedLinAlg::Invert(DenseTensor &A) const
    magma_int_t status;
 
    status = MFEM_MAGMA_PREFIX(getrf_batched)(
-               n, n, d_LU_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
+               n, n, d_A_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
                Magma::Queue());
    MFEM_VERIFY(status == MAGMA_SUCCESS, "");
 

linalg/hypre.cpp

@@ -561,8 +561,7 @@ void CopyMemory(Memory<T> &src, Memory<T> &dst, MemoryClass dst_mc,
     this function. In particular, @a dst should be empty or deleted before
     calling this function. */
 template <typename SrcT, typename DstT>
-void CopyConvertMemory(const Memory<SrcT> &src, MemoryClass dst_mc,
-                       Memory<DstT> &dst)
+void CopyConvertMemory(Memory<SrcT> &src, MemoryClass dst_mc, Memory<DstT> &dst)
 {
    auto capacity = src.Capacity();
    dst.New(capacity, GetMemoryType(dst_mc));
@@ -843,8 +842,8 @@ static int GetPartitioningArraySize(MPI_Comm comm)
 ///
 /// Both @a row and @a col are partitioning arrays, whose length is returned by
 /// GetPartitioningArraySize(), see @ref hypre_partitioning_descr.
-static bool RowAndColStartsAreEqual(MPI_Comm comm, const HYPRE_BigInt *rows,
-                                    const HYPRE_BigInt *cols)
+static bool RowAndColStartsAreEqual(MPI_Comm comm, HYPRE_BigInt *rows,
+                                    HYPRE_BigInt *cols)
 {
    const int part_size = GetPartitioningArraySize(comm);
    bool are_equal = true;
@@ -1132,7 +1131,7 @@ HypreParMatrix::HypreParMatrix(
 HypreParMatrix::HypreParMatrix(MPI_Comm comm,
                                HYPRE_BigInt *row_starts,
                                HYPRE_BigInt *col_starts,
-                               const SparseMatrix *sm_a)
+                               SparseMatrix *sm_a)
 {
    MFEM_ASSERT(sm_a != NULL, "invalid input");
    MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
@@ -1146,7 +1145,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,
 
    hypre_CSRMatrixSetDataOwner(csr_a,0);
    MemoryIJData mem_a;
-   CopyCSR(const_cast<SparseMatrix*>(sm_a), mem_a, csr_a, false);
+   CopyCSR(sm_a, mem_a, csr_a, false);
    hypre_CSRMatrixSetRownnz(csr_a);
 
    // NOTE: this call creates a matrix on host even when device support is
@@ -1308,11 +1307,10 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows,
                                HYPRE_BigInt glob_nrows,
                                HYPRE_BigInt glob_ncols,
-                               const int *I,
-                               const HYPRE_BigInt *J,
-                               const real_t *data,
-                               const HYPRE_BigInt *rows,
-                               const HYPRE_BigInt *cols)
+                               int *I, HYPRE_BigInt *J,
+                               real_t *data,
+                               HYPRE_BigInt *rows,
+                               HYPRE_BigInt *cols)
 {
    Init();
 
@@ -2329,8 +2327,8 @@ void HypreParMatrix::Threshold(real_t threshold)
    /* TODO: GenerateDiagAndOffd() uses an int array of size equal to the number
       of columns in csr_A_wo_z which is the global number of columns in A. This
       does not scale well. */
-   ierr += hypre_GenerateDiagAndOffd(csr_A_wo_z,parcsr_A_ptr,
-                                     col_start,col_end);
+   ierr += GenerateDiagAndOffd(csr_A_wo_z,parcsr_A_ptr,
+                               col_start,col_end);
 
    ierr += hypre_CSRMatrixDestroy(csr_A_wo_z);
 

linalg/hypre.hpp

@@ -25,18 +25,11 @@
 #define HYPRE_TIMING
 
 // hypre header files
-#if MFEM_HYPRE_VERSION < 30000
 #include <seq_mv.h>
 #include <temp_multivector.h>
-#else
-#include <_hypre_seq_mv.h>
-#include <_hypre_lobpcg_temp_multivector.h>
-#endif
 #include <_hypre_parcsr_mv.h>
 #include <_hypre_parcsr_ls.h>
 
-#include <HYPRE_parcsr_ls.h>
-
 #ifdef HYPRE_COMPLEX
 #error "MFEM does not work with HYPRE's complex numbers support"
 #endif
@@ -60,10 +53,6 @@
 #error "MFEM_USE_HIP=YES is required when HYPRE is built with HIP!"
 #endif
 
-#if MFEM_HYPRE_VERSION > 21500
-#define HYPRE_AssumedPartitionCheck() 1
-#endif
-
 namespace mfem
 {
 
@@ -565,7 +554,7 @@ public:
        partitioning arrays @a row_starts and @a col_starts. */
    HypreParMatrix(MPI_Comm comm, HYPRE_BigInt *row_starts,
                   HYPRE_BigInt *col_starts,
-                  const SparseMatrix *a); // constructor with 4 arguments, v2
+                  SparseMatrix *a); // constructor with 4 arguments, v2
 
    /// Creates boolean block-diagonal rectangular parallel matrix.
    /** The new HypreParMatrix does not take ownership of any of the input
@@ -594,9 +583,9 @@ public:
        arrays (so they can be deleted). See @ref hypre_partitioning_descr "here"
        for a description of the partitioning arrays @a rows and @a cols. */
    HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_BigInt glob_nrows,
-                  HYPRE_BigInt glob_ncols, const int *I, const HYPRE_BigInt *J,
-                  const real_t *data, const HYPRE_BigInt *rows,
-                  const HYPRE_BigInt *cols); // constructor with 9 arguments
+                  HYPRE_BigInt glob_ncols, int *I, HYPRE_BigInt *J,
+                  real_t *data, HYPRE_BigInt *rows,
+                  HYPRE_BigInt *cols); // constructor with 9 arguments
 
    /** @brief Copy constructor for a ParCSR matrix which creates a deep copy of
        structure and data from @a P. */

linalg/hypre_parcsr.cpp

@@ -1916,9 +1916,9 @@ hypre_ParCSRMatrixAdd(hypre_ParCSRMatrix *A,
       /* FIXME: GenerateDiagAndOffd() uses an int array of size equal to the
          number of columns in csr_C_temp which is the global number of columns
          in A and B. This does not scale well. */
-      ierr += hypre_GenerateDiagAndOffd(csr_C_temp, C,
-                                        hypre_ParCSRMatrixFirstColDiag(A),
-                                        hypre_ParCSRMatrixLastColDiag(A));
+      ierr += GenerateDiagAndOffd(csr_C_temp, C,
+                                  hypre_ParCSRMatrixFirstColDiag(A),
+                                  hypre_ParCSRMatrixLastColDiag(A));
 
       /* delete CSR version of C */
       ierr += hypre_CSRMatrixDestroy(csr_C_temp);

linalg/hypre_parcsr.hpp

@@ -21,10 +21,6 @@
 // hypre header files
 #include <_hypre_parcsr_mv.h>
 
-#if MFEM_HYPRE_VERSION < 30000
-#define hypre_GenerateDiagAndOffd GenerateDiagAndOffd
-#endif
-
 // Older hypre versions do not define HYPRE_BigInt and HYPRE_MPI_BIG_INT, so we
 // define them here for backward compatibility.
 #if MFEM_HYPRE_VERSION < 21600

linalg/sparsemat.cpp

@@ -46,12 +46,6 @@
 #define MFEM_GPUSPARSE_ALG HIPSPARSE_CSRMV_ALG1
 #endif // defined(MFEM_USE_CUDA)
 
-#if defined(MFEM_USE_SINGLE)
-#define MFEM_CUDA_or_HIP_REAL_T MFEM_CUDA_or_HIP(_R_32F)
-#elif defined(MFEM_USE_DOUBLE)
-#define MFEM_CUDA_or_HIP_REAL_T MFEM_CUDA_or_HIP(_R_64F)
-#endif
-
 namespace mfem
 {
 
@@ -63,10 +57,8 @@ int SparseMatrix::SparseMatrixCount = 0;
 /// @cond Suppress_Doxygen_warnings
 MFEM_cu_or_hip(sparseHandle_t) SparseMatrix::handle = nullptr;
 /// @endcond
-#ifndef MFEM_CUDA_1897_WORKAROUND
 size_t SparseMatrix::bufferSize = 0;
 void * SparseMatrix::dBuffer = nullptr;
-#endif
 #endif // MFEM_USE_CUDA_OR_HIP
 
 void SparseMatrix::InitGPUSparse()
@@ -472,67 +464,109 @@ void SparseMatrix::SortColumnIndices()
    }
 
 #ifdef MFEM_USE_CUDA_OR_HIP
-   if (Device::Allows(Backend::CUDA_MASK) || Device::Allows(Backend::HIP_MASK))
+   if ( Device::Allows( Backend::CUDA_MASK ))
    {
-      const int m = Height();
-      const int n = Width();
+#if defined(MFEM_USE_CUDA)
+      size_t pBufferSizeInBytes = 0;
+      void *pBuffer = NULL;
+
+      const int n = Height();
+      const int m = Width();
       const int nnzA = J.Capacity();
-      const int *d_ia = ReadI();
-      int *d_ja = ReadWriteJ();
+      real_t * d_a_sorted = ReadWriteData();
+      const int * d_ia = ReadI();
+      int * d_ja_sorted = ReadWriteJ();
+      csru2csrInfo_t sortInfoA;
 
-      // Get size of temporary buffer needed to sort the column indices,
-      // allocate the temporary buffer.
-      size_t pBufferSizeInBytes;
-      MFEM_cu_or_hip(sparseXcsrsort_bufferSizeExt)(handle, m, n, nnzA, d_ia,
-                                                   d_ja, &pBufferSizeInBytes);
-      void *pBuffer = MFEM_Cu_or_Hip(MemAlloc)(&pBuffer, pBufferSizeInBytes);
+      cusparseMatDescr_t matA_descr;
+      cusparseCreateMatDescr( &matA_descr );
+      cusparseSetMatIndexBase( matA_descr, CUSPARSE_INDEX_BASE_ZERO );
+      cusparseSetMatType( matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL );
 
-      // Create matrix descriptor, will have default values
-      // CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_MATRIX_TYPE_GENERAL.
-      MFEM_cu_or_hip(sparseMatDescr_t) matA_descr;
-      MFEM_cu_or_hip(sparseCreateMatDescr)(&matA_descr);
+      cusparseCreateCsru2csrInfo( &sortInfoA );
 
-      // Initialize permutation to identity
-      Array<int> P(nnzA);
-      int *d_P = P.Write();
-      mfem::forall(nnzA, [=] MFEM_HOST_DEVICE (int i) { d_P[i] = i; });
+#ifdef MFEM_USE_SINGLE
+      cusparseScsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
+                                       d_ja_sorted, sortInfoA,
+                                       &pBufferSizeInBytes);
+#elif defined MFEM_USE_DOUBLE
+      cusparseDcsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
+                                       d_ja_sorted, sortInfoA,
+                                       &pBufferSizeInBytes);
+#else
+      MFEM_ABORT("Floating point type undefined");
+#endif
 
-      // Sort the column indices. The array d_ja will now be sorted. The
-      // permutation required to sort the values will be returned in d_P.
-      MFEM_cu_or_hip(sparseXcsrsort)(handle, m, n, nnzA, matA_descr, d_ia, d_ja,
-                                     d_P, pBuffer);
+      CuMemAlloc( &pBuffer, pBufferSizeInBytes );
 
-      // Create a copy of the unsorted matrix values.
-      real_t *d_a = ReadWriteData();
-      void *d_a_unsorted = MFEM_Cu_or_Hip(MemAlloc)(&d_a_unsorted,
-                                                    nnzA * sizeof(real_t));
-      MFEM_Cu_or_Hip(MemcpyDtoD)(d_a_unsorted, d_a, nnzA * sizeof(real_t));
+#ifdef MFEM_USE_SINGLE
+      cusparseScsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
+                         d_ja_sorted, sortInfoA, pBuffer);
+#elif defined MFEM_USE_DOUBLE
+      cusparseDcsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
+                         d_ja_sorted, sortInfoA, pBuffer);
+#else
+      MFEM_ABORT("Floating point type undefined");
+#endif
 
-      // Create the (input) dense vector with the unsorted values.
-      MFEM_cu_or_hip(sparseDnVecDescr_t) d_a_dense;
-      MFEM_cu_or_hip(sparseCreateDnVec)(&d_a_dense, nnzA, d_a_unsorted,
-                                        MFEM_CUDA_or_HIP_REAL_T);
-
-      // Create the (output) sparse vector that will have the sorted values.
-      MFEM_cu_or_hip(sparseSpVecDescr_t) d_a_sparse;
-      MFEM_cu_or_hip(sparseCreateSpVec)(&d_a_sparse, nnzA, nnzA, d_P, d_a,
-                                        MFEM_CU_or_HIP(SPARSE_INDEX_32I),
-                                        MFEM_CU_or_HIP(SPARSE_INDEX_BASE_ZERO),
-                                        MFEM_CUDA_or_HIP_REAL_T);
-
-      // Sort the matrix values using the permutation vector.
-      MFEM_cu_or_hip(sparseGather)(handle, d_a_dense, d_a_sparse);
-
-      // The above calls may be asynchronous, so we need to wait for them to
-      // finish before we can free memory.
+      // The above call is (at least in some cases) asynchronous, so we need to
+      // wait for it to finish before we can free device temporaries.
       MFEM_STREAM_SYNC;
 
-      MFEM_cu_or_hip(sparseDestroyDnVec)(d_a_dense);
-      MFEM_cu_or_hip(sparseDestroySpVec)(d_a_sparse);
-      MFEM_cu_or_hip(sparseDestroyMatDescr)(matA_descr);
+      cusparseDestroyCsru2csrInfo( sortInfoA );
+      cusparseDestroyMatDescr( matA_descr );
 
-      MFEM_Cu_or_Hip(MemFree)(d_a_unsorted);
-      MFEM_Cu_or_Hip(MemFree)(pBuffer);
+      CuMemFree( pBuffer );
+#endif
+   }
+   else if ( Device::Allows( Backend::HIP_MASK ))
+   {
+#if defined(MFEM_USE_HIP)
+      size_t pBufferSizeInBytes = 0;
+      void *pBuffer = NULL;
+      int *P = NULL;
+
+      const int n = Height();
+      const int m = Width();
+      const int nnzA = J.Capacity();
+      real_t * d_a_sorted = ReadWriteData();
+      const int * d_ia = ReadI();
+      int * d_ja_sorted = ReadWriteJ();
+
+      hipsparseMatDescr_t descrA;
+      hipsparseCreateMatDescr( &descrA );
+      // FIXME: There is not in-place version of csr sort in hipSPARSE currently, so we make
+      //        a temporary copy of the data for gthr, sort that, and then copy the sorted values
+      //        back to the array being returned. Where there is an in-place version available,
+      //        we should use it.
+      Array< real_t > a_tmp( nnzA );
+      real_t *d_a_tmp = a_tmp.Write();
+
+      hipsparseXcsrsort_bufferSizeExt(handle, n, m, nnzA, d_ia, d_ja_sorted,
+                                      &pBufferSizeInBytes);
+
+      HipMemAlloc( &pBuffer, pBufferSizeInBytes );
+      HipMemAlloc( (void**)&P, nnzA * sizeof(int) );
+
+      hipsparseCreateIdentityPermutation(handle, nnzA, P);
+      hipsparseXcsrsort(handle, n, m, nnzA, descrA, d_ia, d_ja_sorted, P, pBuffer);
+
+#if defined(MFEM_USE_SINGLE)
+      hipsparseSgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
+                     HIPSPARSE_INDEX_BASE_ZERO);
+#elif defined(MFEM_USE_DOUBLE)
+      hipsparseDgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
+                     HIPSPARSE_INDEX_BASE_ZERO);
+#else
+      MFEM_ABORT("Unsupported floating point type!");
+#endif
+
+      A.CopyFrom( a_tmp.GetMemory(), nnzA );
+      hipsparseDestroyMatDescr( descrA );
+
+      HipMemFree( pBuffer );
+      HipMemFree( P );
+#endif
    }
    else
 #endif // MFEM_USE_CUDA_OR_HIP
@@ -787,15 +821,27 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
             MFEM_CU_or_HIP(SPARSE_INDEX_32I),
             MFEM_CU_or_HIP(SPARSE_INDEX_32I),
             MFEM_CU_or_HIP(SPARSE_INDEX_BASE_ZERO),
-            MFEM_CUDA_or_HIP_REAL_T);
+#ifdef MFEM_USE_SINGLE
+            MFEM_CUDA_or_HIP(_R_32F));
+#else
+            MFEM_CUDA_or_HIP(_R_64F));
+#endif
 
          // Create handles for input/output vectors
          MFEM_cu_or_hip(sparseCreateDnVec)(&vecX_descr,
                                            x.Size(),
                                            const_cast<real_t *>(d_x),
-                                           MFEM_CUDA_or_HIP_REAL_T);
+#ifdef MFEM_USE_SINGLE
+                                           MFEM_CUDA_or_HIP(_R_32F));
+#else
+                                           MFEM_CUDA_or_HIP(_R_64F));
+#endif
          MFEM_cu_or_hip(sparseCreateDnVec)(&vecY_descr, y.Size(), d_y,
-                                           MFEM_CUDA_or_HIP_REAL_T);
+#ifdef MFEM_USE_SINGLE
+                                           MFEM_CUDA_or_HIP(_R_32F));
+#else
+                                           MFEM_CUDA_or_HIP(_R_64F));
+#endif
 #else
          cusparseCreateMatDescr(&matA_descr);
          cusparseSetMatIndexBase(matA_descr, CUSPARSE_INDEX_BASE_ZERO);
@@ -814,7 +860,11 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
          vecX_descr,
          &beta,
          vecY_descr,
-         MFEM_CUDA_or_HIP_REAL_T,
+#ifdef MFEM_USE_SINGLE
+         MFEM_CUDA_or_HIP(_R_32F),
+#else
+         MFEM_CUDA_or_HIP(_R_64F),
+#endif
          MFEM_GPUSPARSE_ALG,
          &newBufferSize);
 
@@ -841,7 +891,11 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
          vecX_descr,
          &beta,
          vecY_descr,
-         MFEM_CUDA_or_HIP_REAL_T,
+#ifdef MFEM_USE_SINGLE
+         MFEM_CUDA_or_HIP(_R_32F),
+#else
+         MFEM_CUDA_or_HIP(_R_64F),
+#endif
          MFEM_GPUSPARSE_ALG,
          dBuffer);
 #else
@@ -4318,14 +4372,6 @@ SparseMatrix::~SparseMatrix()
 #ifdef MFEM_USE_CUDA_OR_HIP
    if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-#ifdef MFEM_CUDA_1897_WORKAROUND
-      if (dBuffer)
-      {
-         MFEM_Cu_or_Hip(MemFree)(dBuffer);
-         dBuffer = nullptr;
-         bufferSize = 0;
-      }
-#endif
       if (SparseMatrixCount==1)
       {
          if (handle)
@@ -4333,14 +4379,12 @@ SparseMatrix::~SparseMatrix()
             MFEM_cu_or_hip(sparseDestroy)(handle);
             handle = nullptr;
          }
-#ifndef MFEM_CUDA_1897_WORKAROUND
          if (dBuffer)
          {
             MFEM_Cu_or_Hip(MemFree)(dBuffer);
             dBuffer = nullptr;
             bufferSize = 0;
          }
-#endif
       }
       SparseMatrixCount--;
    }

linalg/sparsemat.hpp

@@ -98,17 +98,9 @@ protected:
 #ifdef MFEM_USE_CUDA_OR_HIP
    // common for hipSPARSE and cuSPARSE
    static int SparseMatrixCount;
-   mutable bool initBuffers = false;
-
-#if defined(MFEM_USE_CUDA) && CUDA_VERSION >= 12300 && CUDA_VERSION < 12602
-   // Workaround for bug CUSPARSE-1897
-#define MFEM_CUDA_1897_WORKAROUND
-   mutable size_t bufferSize = 0;
-   mutable void *dBuffer = nullptr;
-#else
    static size_t bufferSize;
    static void *dBuffer;
-#endif
+   mutable bool initBuffers = false;
 
 #if defined(MFEM_USE_CUDA)
    cusparseStatus_t status;

linalg/vector.cpp

@@ -92,6 +92,18 @@ struct LpReducer
    }
 };
 
+static Array<real_t>& vector_workspace()
+{
+   static Array<real_t> instance;
+   return instance;
+}
+
+static Array<DevicePair<real_t, real_t>> &Lpvector_workspace()
+{
+   static Array<DevicePair<real_t, real_t>> instance;
+   return instance;
+}
+
 Vector::Vector(const Vector &v)
 {
    const int s = v.Size();
@@ -979,7 +991,7 @@ real_t Vector::Norml2() const
          }
       }
    },
-   L2Reducer{}, UseDevice());
+   L2Reducer{}, UseDevice(), Lpvector_workspace());
    // final answer
    return res.second * sqrt(res.first);
 }
@@ -994,7 +1006,7 @@ real_t Vector::Normlinf() const
    {
       r = fmax(r, fabs(m_data[i]));
    },
-   MaxReducer<real_t> {}, UseDevice());
+   MaxReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 
@@ -1008,7 +1020,7 @@ real_t Vector::Norml1() const
    {
       r += fabs(m_data[i]);
    },
-   SumReducer<real_t> {}, UseDevice());
+   SumReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 
@@ -1051,7 +1063,7 @@ real_t Vector::Normlp(real_t p) const
             }
          }
       },
-      LpReducer{p}, UseDevice());
+      LpReducer{p}, UseDevice(), Lpvector_workspace());
       // final answer
       return res.second * pow(res.first, 1.0 / p);
    } // end if p < infinity()
@@ -1084,7 +1096,7 @@ real_t Vector::operator*(const Vector &v) const
       {
          r += m_data[i] * v_data[i];
       },
-      SumReducer<real_t> {}, use_dev);
+      SumReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1155,7 +1167,7 @@ real_t Vector::Min() const
       {
          r = fmin(r, m_data[i]);
       },
-      MinReducer<real_t> {}, use_dev);
+      MinReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1201,7 +1213,7 @@ real_t Vector::Max() const
       {
          r = fmax(r, m_data[i]);
       },
-      MaxReducer<real_t> {}, use_dev);
+      MaxReducer<real_t> {}, use_dev, vector_workspace());
       return res;
    };
 
@@ -1236,7 +1248,7 @@ real_t Vector::Sum() const
    {
       r += m_data[i];
    },
-   SumReducer<real_t> {}, UseDevice());
+   SumReducer<real_t> {}, UseDevice(), vector_workspace());
    return res;
 }
 

makefile

@@ -377,7 +377,7 @@ MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
  MFEM_INC_DIR MFEM_TPLFLAGS MFEM_INCFLAGS MFEM_PICFLAG MFEM_FLAGS MFEM_LIB_DIR\
  MFEM_EXT_LIBS MFEM_LIBS MFEM_LIB_FILE MFEM_STATIC MFEM_SHARED MFEM_BUILD_TAG\
  MFEM_PREFIX MFEM_CONFIG_EXTRA MFEM_MPIEXEC MFEM_MPIEXEC_NP MFEM_MPI_NP\
- MFEM_TEST_MK MFEM_XLINKER
+ MFEM_TEST_MK
 
 # Config vars: values of the form @VAL@ are replaced by $(VAL) in config.mk
 MFEM_CPPFLAGS  ?= $(CPPFLAGS)
@@ -394,7 +394,6 @@ MFEM_BUILD_TAG ?= $(shell uname -snm)
 MFEM_PREFIX    ?= $(PREFIX)
 MFEM_INC_DIR   ?= $(if $(CONFIG_FILE_DEF),@MFEM_BUILD_DIR@,@MFEM_DIR@)
 MFEM_LIB_DIR   ?= $(if $(CONFIG_FILE_DEF),@MFEM_BUILD_DIR@,@MFEM_DIR@)
-MFEM_XLINKER   ?= $(XLINKER)
 MFEM_TEST_MK   ?= @MFEM_DIR@/config/test.mk
 # Use "\n" (interpreted by sed) to add a newline.
 MFEM_CONFIG_EXTRA ?= $(if $(CONFIG_FILE_DEF),MFEM_BUILD_DIR ?= @MFEM_DIR@,)

mesh/mesh.cpp

@@ -540,6 +540,7 @@ void Mesh::GetFaceTransformation(int FaceNo,
                                  IsoparametricTransformation *FTr) const
 {
    FTr->Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
+   FTr->Attribute = faces[FaceNo]->GetAttribute();
    FTr->ElementNo = FaceNo;
    FTr->ElementType = ElementTransformation::FACE;
    FTr->mesh = this;

mesh/mesh.hpp

@@ -2244,6 +2244,12 @@ public:
 
    /// @}
 
+   /// Return the attribute of element i.
+   int GetFaceAttribute(int i) const { return faces[i]->GetAttribute(); }
+
+   /// Set the attribute of face element i.
+   void SetFaceAttribute(int i, int attr) { faces[i]->SetAttribute(attr); }
+
    /// @name Methods related to mesh partitioning
    /// @{
 

mesh/mesh_readers.cpp

@@ -773,7 +773,7 @@ struct BufferReader : BufferReaderBase
          int header_entry_size = HeaderEntrySize();
          int nblocks = ReadHeaderEntry(header_buf);
          header_buf += header_entry_size;
-         std::vector<size_t> header(nblocks + 2);
+         std::vector<int> header(nblocks + 2);
          for (int i=0; i<nblocks+2; ++i)
          {
             header[i] = ReadHeaderEntry(header_buf);
@@ -792,7 +792,7 @@ struct BufferReader : BufferReaderBase
             dest_ptr += dest_len;
             source_ptr += source_len;
          }
-         MFEM_VERIFY(size_t(sizeof(F)*n) == (dest_ptr - dest_start),
+         MFEM_VERIFY(int(sizeof(F)*n) == (dest_ptr - dest_start),
                      "AppendedData: wrong data size");
          buf = uncompressed_data.data();
 #else

miniapps/autodiff/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/autodiff/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 

miniapps/common/makefile

@@ -20,7 +20,6 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 # Default target
 all: lib-common
 
-# Include defaults.mk to get the definition of BUILD_SOFLAGS
 DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
 include $(DEFAULTS_MK)
 
@@ -31,7 +30,13 @@ MFEM_LIB_FILE = mfem_is_not_built
 ifneq (clean,$(MAKECMDGOALS))
    -include $(CONFIG_MK)
 
-   XLINKER = $(MFEM_XLINKER)
+   ifeq ($(MFEM_USE_CUDA),YES)
+      XLINKER = $(CUDA_XLINKER)
+   else ifeq ($(MFEM_USE_HIP),YES)
+      XLINKER = $(HIP_XLINKER)
+   else
+      XLINKER = $(CXX_XLINKER)
+   endif
 
    BUILD_REAL_DIR = $(realpath .)
    BUILD_SOFLAGS := $(subst libmfem.,libmfem-common.,$(BUILD_SOFLAGS))

miniapps/diag-smoothers/makefile

@@ -15,9 +15,11 @@ MFEM_BUILD_DIR ?= ../..
 
 SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/diag-smoothers/,)
 CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
 
 MFEM_LIB_FILE = mfem_is_not_built
 
+-include $(DEFAULTS_MK)
 -include $(CONFIG_MK)
 
 DS_COMMON_SRC = ds-common.cpp
@@ -29,7 +31,8 @@ MINIAPPS = $(if $(MFEM_USE_MPI:NO=),$(PAR_MINIAPPS),)
 
 COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
 
 APP_DEPS = $(DS_COMMON_OBJ) $(MFEM_LIB_FILE) $(CONFIG_MK)
 APP_LIBS = $(COMMON_LIB) $(MFEM_LIBS)

miniapps/dpg/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/dpg/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -65,7 +69,8 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk

miniapps/electromagnetics/CMakeLists.txt

@@ -34,11 +34,6 @@ if (MFEM_USE_MPI)
    EXTRA_HEADERS maxwell_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
    LIBRARIES mfem-common)
 
- add_mfem_miniapp(lorentz
-   MAIN lorentz.cpp
-   EXTRA_HEADERS ${MFEM_MINIAPPS_COMMON_HEADERS}
-   LIBRARIES mfem-common)
-
   # Add the corresponding tests to the "test" target
   if (MFEM_ENABLE_TESTING)
     add_test(NAME tesla_np=4

miniapps/electromagnetics/lorentz.cpp

@@ -1,571 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-//
-//           -----------------------------------------------------
-//           Lorentz Miniapp:  Simple Lorentz Force Particle Mover
-//           -----------------------------------------------------
-//
-// This miniapp computes the trajectory of a single charged particle subject to
-// Lorentz forces.
-//
-//                           dp/dt = q (E + v x B)
-//
-// The method used is the explicit Boris algortihm which conserves phase space
-// volume for long term accuracy.
-//
-// The electric and magnetic fields are read from VisItDataCollection objects
-// such as those produced by the Volta and Tesla miniapps. It is notable that
-// these two fields do not need to be defined on the same mesh. Of course, the
-// particle trajectory can only be computed on the intersection of the two
-// domains. The starting point of the path must be chosen within in this
-// intersection and the trajectory will terminate when it leaves the
-// intersection or reaches a specified time duration.
-//
-// Note that the VisItDataCollection objects must have been stored using the
-// parallel format e.g. visit_dc.SetFormat(DataCollection::PARALLEL_FORMAT);.
-// Without this optional format specifier the vector field lookups will fail.
-//
-// Compile with: make lorentz
-//
-// Sample runs:
-//
-//   Free particle moving with constant velocity
-//      mpirun -np 4 lorentz -p0 '1 1 1'
-//
-//   Particle accelerating in a constant electric field
-//      mpirun -np 4 volta -m ../../data/inline-hex.mesh -dbcs '1 6' -dbcv '0 1'
-//      mpirun -np 4 lorentz -er Volta-AMR-Parallel -x0 '0.5 0.5 0.9' -p0 '1 0 0'
-//
-//   Particle accelerating in a constant magnetic field
-//      mpirun -np 4 tesla -m ../../data/inline-hex.mesh -ubbc '0 0 1'
-//      mpirun -np 4 lorentz -br Tesla-AMR-Parallel -x0 '0.1 0.5 0.1' -p0 '0 0.4 0.1' -tf 9
-//
-//   Magnetic mirror effect near a charged sphere and a bar magnet
-//      mpirun -np 4 volta -m ../../data/ball-nurbs.mesh -dbcs 1 -cs '0 0 0 0.1 2e-11' -rs 2 -maxit 4
-//      mpirun -np 4 tesla -m ../../data/fichera.mesh -maxit 4 -rs 3 -bm '-0.1 -0.1 -0.1 0.1 0.1 0.1 0.1 -1e10'
-//      mpirun -np 4 lorentz -er Volta-AMR-Parallel -ec 4 -br Tesla-AMR-Parallel -bc 4 -x0 '0.8 0 0' -p0 '-8 -4 4' -q -10 -tf 0.2 -dt 1e-3 -rf 1e-6
-//
-// This miniapp demonstrates the use of the ParMesh::FindPoints functionality
-// to evaluate field data from stored DataCollection objects.  While this
-// miniapp is far from a full particle-in-cell (PIC) code it does demonstrate
-// some of the building blocks that might be used to construct the particle
-// mover portion of a PIC code.
-
-#include "mfem.hpp"
-#include "../common/fem_extras.hpp"
-#include "../common/pfem_extras.hpp"
-#include "electromagnetics.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-using namespace mfem::common;
-using namespace mfem::electromagnetics;
-
-typedef DataCollection::FieldMapType fields_t;
-
-/// This class implements the Boris algorithm as described in the
-/// article `Why is Boris algorithm so good?` by H. Qin et al in
-/// Physics of Plasmas, Volume 20 Issue 8, August 2013,
-/// https://doi.org/10.1063/1.4818428.
-class BorisAlgorithm
-{
-private:
-   real_t charge_;
-   real_t mass_;
-
-   ParMesh         *E_pmesh_;
-   ParGridFunction *E_field_;
-
-   ParMesh         *B_pmesh_;
-   ParGridFunction *B_field_;
-
-   mutable Array<int>              elem_id_;
-   mutable Array<IntegrationPoint> ip_;
-
-   mutable Vector E_;
-   mutable Vector B_;
-   mutable Vector pxB_;
-   mutable Vector pm_;
-   mutable Vector pp_;
-
-   // Returns true if a usable V has been found. If @a pgf is NULL, V = 0 is
-   // returned as a default value.
-   bool GetValue(ParMesh *pmesh, ParGridFunction *pgf, Vector q, Vector &V)
-   {
-      DenseMatrix point(q.GetData(), 3, 1);
-
-      int pt_found =
-         (pmesh != NULL) ? pmesh->FindPoints(point, elem_id_, ip_, false) : -1;
-
-      // We have a mesh but the point was not found. The path must be outside
-      // the domain of interest.
-      if (pmesh != NULL && pt_found <= 0) { return false; }
-
-      int pt_root = -1;
-
-      if (pt_found > 0 && elem_id_[0] >= 0 && pgf != NULL)
-      {
-         pt_root = pmesh->GetMyRank();
-
-         pgf->GetVectorValue(elem_id_[0], ip_[0], V);
-      }
-      else
-      {
-         pt_root = 0;
-         V = 0.0;
-      }
-
-      // Determine processor which found the field point
-      int glb_pt_root = -1;
-      MPI_Allreduce(&pt_root, &glb_pt_root, 1,
-                    MPI_INT, MPI_MAX, MPI_COMM_WORLD);
-
-      // Send the field value to the root processor
-      if (pmesh != NULL && elem_id_[0] >= 0 && glb_pt_root != 0)
-      {
-         MPI_Send(V.GetData(), 3, MPITypeMap<real_t>::mpi_type,
-                  0, 1030, MPI_COMM_WORLD);
-      }
-
-      // Receive the field value on the root processor
-      if (Mpi::Root() && pmesh != NULL && glb_pt_root != 0)
-      {
-         MPI_Status status;
-         MPI_Recv(V.GetData(), 3, MPITypeMap<real_t>::mpi_type,
-                  glb_pt_root, 1030, MPI_COMM_WORLD, &status);
-      }
-      return true;
-   }
-
-public:
-   BorisAlgorithm(ParGridFunction *E_gf,
-                  ParGridFunction *B_gf,
-                  real_t charge, real_t mass)
-      : charge_(charge), mass_(mass),
-        E_field_(E_gf),
-        B_field_(B_gf),
-        E_(3), B_(3), pxB_(3), pm_(3), pp_(3)
-   {
-      E_pmesh_ = (E_field_) ? E_field_->ParFESpace()->GetParMesh() : NULL;
-      B_pmesh_ = (B_field_) ? B_field_->ParFESpace()->GetParMesh() : NULL;
-   }
-
-   bool Step(Vector &q, Vector &p, real_t &t, real_t &dt)
-   {
-      // Locate current point in each mesh, evaluate the fields, and collect
-      // field values on the root processor.
-      if (!GetValue(E_pmesh_, E_field_, q, E_)) { return false; }
-      if (!GetValue(B_pmesh_, B_field_, q, B_)) { return false; }
-
-      // Compute updated position and momentum using the Boris algorithm
-      if (Mpi::Root())
-      {
-         // Compute half of the contribution from q E
-         add(p, 0.5 * dt * charge_, E_, pm_);
-
-         // Compute the contributiobn from q p x B
-         const real_t B2 = B_ * B_;
-
-         // ... along pm x B
-         const real_t a1 = 4.0 * dt * charge_ * mass_;
-         pm_.cross3D(B_, pxB_);
-         pp_.Set(a1, pxB_);
-
-         // ... along pm
-         const real_t a2 = 4.0 * mass_ * mass_ -
-                           dt * dt * charge_ * charge_ * B2;
-         pp_.Add(a2, pm_);
-
-         // ... along B
-         const real_t a3 = 2.0 * dt * dt * charge_ * charge_ * (B_ * p);
-         pp_.Add(a3, B_);
-
-         // scale by common denominator
-         const real_t a4 = 4.0 * mass_ * mass_ +
-                           dt * dt * charge_ * charge_ * B2;
-         pp_ /= a4;
-
-         // Update the momentum
-         add(pp_, 0.5 * dt * charge_, E_, p);
-
-         // Update the position
-         q.Add(dt / mass_, p);
-      }
-
-      // Update the time
-      t += dt;
-
-      // Broadcast the updated position
-      MPI_Bcast(q.GetData(), 3, MPITypeMap<real_t>::mpi_type,
-                0, MPI_COMM_WORLD);
-
-      // Broadcast the updated momentum
-      MPI_Bcast(p.GetData(), 3, MPITypeMap<real_t>::mpi_type,
-                0, MPI_COMM_WORLD);
-
-      return true;
-   }
-};
-
-// Open the named VisItDataCollection and read the named field.
-// Returns pointers to the two new objects.
-int ReadGridFunction(const char * coll_name, const char * field_name,
-                     int pad_digits_cycle, int pad_digits_rank, int cycle,
-                     VisItDataCollection *&dc, ParGridFunction *& gf);
-
-// By default the initial position will be the center of the intersection
-// of the bounding boxes of the meshes containing the E and B fields.
-void SetInitialPosition(VisItDataCollection *E_dc,
-                        VisItDataCollection *B_dc,
-                        Vector &x_init);
-
-// Build a quadrilateral mesh approximating the trajectory as a
-// ribbon. One edge of the ribbon follows the trajectory of the
-// particle. The opposite edge is offset by the acceleration vector
-// (scaled by a constant called the r_factor).
-Mesh MakeTrajectoryMesh(int step, real_t m, real_t dt, real_t r_factor,
-                        const DenseMatrix &pos_data,
-                        const DenseMatrix &mom_data);
-
-// Prints the program's logo to the given output stream
-void display_banner(ostream & os);
-
-int main(int argc, char *argv[])
-{
-   Mpi::Init(argc, argv);
-   Hypre::Init();
-
-   if ( Mpi::Root() ) { display_banner(cout); }
-
-   const char *E_coll_name = "";
-   const char *E_field_name = "E";
-   int E_cycle = 10;
-   int E_pad_digits_cycle = 6;
-   int E_pad_digits_rank = 6;
-
-   const char *B_coll_name = "";
-   const char *B_field_name = "B";
-   int B_cycle = 10;
-   int B_pad_digits_cycle = 6;
-   int B_pad_digits_rank = 6;
-
-   real_t q = 1.0;
-   real_t m = 1.0;
-   real_t dt = 1e-2;
-   real_t t_init = 0.0;
-   real_t t_final = 1.0;
-   real_t r_factor = -1.0;
-   Vector x_init;
-   Vector p_init;
-   int visport = 19916;
-   bool visualization = true;
-   bool visit = true;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&E_coll_name, "-er", "--e-root-file",
-                  "Set the VisIt data collection E field root file prefix.");
-   args.AddOption(&E_field_name, "-ef", "--e-field-name",
-                  "Set the VisIt data collection E field name");
-   args.AddOption(&E_cycle, "-ec", "--e-cycle",
-                  "Set the E field cycle index to read.");
-   args.AddOption(&E_pad_digits_cycle, "-epdc", "--e-pad-digits-cycle",
-                  "Number of digits in E field cycle.");
-   args.AddOption(&E_pad_digits_rank, "-epdr", "--e-pad-digits-rank",
-                  "Number of digits in E field MPI rank.");
-   args.AddOption(&B_coll_name, "-br", "--b-root-file",
-                  "Set the VisIt data collection B field root file prefix.");
-   args.AddOption(&B_field_name, "-bf", "--b-field-name",
-                  "Set the VisIt data collection B field name");
-   args.AddOption(&B_cycle, "-bc", "--b-cycle",
-                  "Set the B field cycle index to read.");
-   args.AddOption(&B_pad_digits_cycle, "-bpdc", "--b-pad-digits-cycle",
-                  "Number of digits in B field cycle.");
-   args.AddOption(&B_pad_digits_rank, "-bpdr", "--b-pad-digits-rank",
-                  "Number of digits in B field MPI rank.");
-   args.AddOption(&q, "-q", "--charge",
-                  "Particle charge.");
-   args.AddOption(&m, "-m", "--mass",
-                  "Particle mass.");
-   args.AddOption(&dt, "-dt", "--time-step",
-                  "Time Step.");
-   args.AddOption(&t_init, "-ti", "--initial-time",
-                  "Initial Time.");
-   args.AddOption(&t_final, "-tf", "--final-time",
-                  "Final Time.");
-   args.AddOption(&x_init, "-x0", "--initial-position",
-                  "Initial position.");
-   args.AddOption(&p_init, "-p0", "--initial-momentum",
-                  "Initial momentum.");
-   args.AddOption(&r_factor, "-rf", "--ribbon-factor",
-                  "Scale factor for ribbon width (rf * (p1-p0) / (m * dt) "
-                  "where p0 and p1 are computed momenta).");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
-                  "Enable or disable VisIt visualization.");
-   args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (Mpi::Root())
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (r_factor <= 0.0)
-   {
-      r_factor = dt;
-   }
-   if (Mpi::Root())
-   {
-      args.PrintOptions(cout);
-   }
-
-   VisItDataCollection *E_dc = NULL;
-   ParGridFunction     *E_gf = NULL;
-
-   if (strcmp(E_coll_name, ""))
-   {
-      if (ReadGridFunction(E_coll_name, E_field_name, E_pad_digits_cycle,
-                           E_pad_digits_rank, E_cycle, E_dc, E_gf))
-      {
-         mfem::out << "Error loading E field" << endl;
-         return 1;
-      }
-   }
-
-   VisItDataCollection *B_dc = NULL;
-   ParGridFunction     *B_gf = NULL;
-
-   if (strcmp(B_coll_name, ""))
-   {
-      if (ReadGridFunction(B_coll_name, B_field_name, B_pad_digits_cycle,
-                           B_pad_digits_rank, B_cycle, B_dc, B_gf))
-      {
-         mfem::out << "Error loading B field" << endl;
-         return 1;
-      }
-   }
-
-   if (x_init.Size() < 3)
-   {
-      SetInitialPosition(E_dc, B_dc, x_init);
-   }
-   if (p_init.Size() < 3)
-   {
-      p_init.SetSize(3); p_init = 0.0;
-   }
-   if (Mpi::Root())
-   {
-      mfem::out << "Initial position: "; x_init.Print(mfem::out);
-      mfem::out << "Initial momentum: "; p_init.Print(mfem::out);
-   }
-
-   BorisAlgorithm boris(E_gf, B_gf, q, m);
-   Vector pos(x_init);
-   Vector mom(p_init);
-
-   ofstream ofs("Lorentz.dat");
-   ofs.precision(14);
-
-   int nsteps = 1 + (int)ceil((t_final - t_init) / dt);
-   DenseMatrix pos_data(3, nsteps);
-   DenseMatrix mom_data(3, nsteps + 1);
-   mom_data.SetCol(0, p_init);
-
-   if (Mpi::Root())
-   {
-      mfem::out << "Maximum number of steps: " << nsteps << endl;
-   }
-
-   int step = -1;
-   real_t t = t_init;
-   do
-   {
-      if (Mpi::Root())
-      {
-         ofs << t
-             << '\t' << pos[0] << '\t' << pos[1] << '\t' << pos[2]
-             << '\t' << mom[0] << '\t' << mom[1] << '\t' << mom[2]
-             << '\n';
-      }
-      step++;
-
-      pos_data.SetCol(step, pos);
-      mom_data.SetCol(step + 1, mom);
-   }
-   while (boris.Step(pos, mom, t, dt) && step < nsteps - 1);
-
-   if (Mpi::Root() && (visit || visualization))
-   {
-      Mesh trajectory = MakeTrajectoryMesh(step, m, dt, r_factor,
-                                           pos_data, mom_data);
-
-      L2_FECollection    fec_l2(0, 2);
-      FiniteElementSpace fes_l2(&trajectory, &fec_l2);
-      GridFunction traj_time(&fes_l2);
-      for (int i=0; i<step; i++)
-      {
-         traj_time[i] = dt * i;
-      }
-
-      if (visit)
-      {
-         VisItDataCollection visit_dc("Lorentz", &trajectory);
-         visit_dc.RegisterField("Time", &traj_time);
-         visit_dc.SetCycle(step);
-         visit_dc.SetTime(step * dt);
-         visit_dc.Save();
-      }
-
-      if (visualization)
-      {
-         socketstream traj_sock;
-         traj_sock.precision(8);
-
-         char vishost[] = "localhost";
-
-         int Wx = 0, Wy = 0; // window position
-         int Ww = 350, Wh = 350; // window size
-
-         VisualizeField(traj_sock, vishost, visport,
-                        traj_time, "Trajectory", Wx, Wy, Ww, Wh);
-      }
-   }
-   if (Mpi::Root())
-   {
-      mfem::out << "Number of steps taken: " << step << endl;
-   }
-
-   // Clean up
-   delete E_dc;
-   delete B_dc;
-}
-
-// Print the Lorentz ascii logo to the given ostream
-void display_banner(ostream & os)
-{
-   os << "   ____                                __          "
-      << endl
-      << "  |    |    ___________   ____   _____/  |_________"
-      << endl
-      << "  |    |   /  _ \\_  __ \\_/ __ \\ /    \\   __\\___   /"
-      << endl
-      << "  |    |__(  <_> )  | \\/\\  ___/|   |  \\  |  /    / "
-      << endl
-      << "  |_______ \\____/|__|    \\___  >___|  /__| /_____ \\"
-      << endl
-      << "          \\/                 \\/     \\/           \\/"
-      << endl << flush;
-}
-
-int ReadGridFunction(const char * coll_name, const char * field_name,
-                     int pad_digits_cycle, int pad_digits_rank, int cycle,
-                     VisItDataCollection *&dc, ParGridFunction *& gf)
-{
-   dc = new VisItDataCollection(MPI_COMM_WORLD, coll_name);
-   dc->SetPadDigitsCycle(pad_digits_cycle);
-   dc->SetPadDigitsRank(pad_digits_rank);
-   dc->Load(cycle);
-
-   if (dc->Error() != DataCollection::No_Error)
-   {
-      mfem::out << "Error loading VisIt data collection: "
-                << coll_name << endl;
-      return 1;
-   }
-
-   if (dc->GetMesh()->Dimension() < 3)
-   {
-      mfem::out << "Field must be defined on a three dimensional mesh"
-                << endl;
-      return 1;
-   }
-
-   if (dc->HasField(field_name))
-   {
-      gf = dc->GetParField(field_name);
-   }
-
-   return 0;
-}
-
-void SetInitialPosition(VisItDataCollection *E_dc,
-                        VisItDataCollection *B_dc,
-                        Vector &x_init)
-{
-   x_init.SetSize(3); x_init = 0.0;
-
-   if (E_dc != NULL || B_dc != NULL)
-   {
-      Vector E_p_min(3); E_p_min = -infinity();
-      Vector E_p_max(3); E_p_max = infinity();
-      if (E_dc != NULL)
-      {
-         ParMesh * E_pmesh = dynamic_cast<ParMesh*>(E_dc->GetMesh());
-         E_pmesh->GetBoundingBox(E_p_min, E_p_max);
-      }
-
-      Vector B_p_min(3); B_p_min = -infinity();
-      Vector B_p_max(3); B_p_max = infinity();
-      if (B_dc != NULL)
-      {
-         ParMesh *B_pmesh = dynamic_cast<ParMesh*>(B_dc->GetMesh());
-         B_pmesh->GetBoundingBox(B_p_min, B_p_max);
-      }
-
-      for (int d = 0; d<3; d++)
-      {
-         const real_t p_min = std::max(E_p_min[d], B_p_min[d]);
-         const real_t p_max = std::min(E_p_max[d], B_p_max[d]);
-         x_init[d] = 0.5 * (p_min + p_max);
-      }
-   }
-}
-
-Mesh MakeTrajectoryMesh(int step, real_t m, real_t dt, real_t r_factor,
-                        const DenseMatrix &pos_data,
-                        const DenseMatrix &mom_data)
-{
-   Mesh trajectory(2, 2 * (step + 1), step, 0, 3);
-
-   for (int i=0; i<=step; i++)
-   {
-      trajectory.AddVertex(pos_data(0,i), pos_data(1,i), pos_data(2,i));
-
-      real_t dpx = (mom_data(0, i + 1) - mom_data(0, i)) / (m * dt);
-      real_t dpy = (mom_data(1, i + 1) - mom_data(1, i)) / (m * dt);
-      real_t dpz = (mom_data(2, i + 1) - mom_data(2, i)) / (m * dt);
-
-      trajectory.AddVertex(pos_data(0,i) + r_factor * dpx,
-                           pos_data(1,i) + r_factor * dpy,
-                           pos_data(2,i) + r_factor * dpz);
-   }
-
-   int v[4];
-   for (int i=0; i<step; i++)
-   {
-      v[0] = 2 * i;
-      v[1] = 2 * (i + 1);
-      v[2] = 2 * (i + 1) + 1;
-      v[3] = 2 * i + 1;
-
-      trajectory.AddQuad(v);
-   }
-
-   trajectory.FinalizeQuadMesh(1);
-
-   return trajectory;
-}

miniapps/electromagnetics/makefile

@@ -17,11 +17,15 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/electromagnetics/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_MINIAPPS =
-PAR_MINIAPPS = volta tesla maxwell joule lorentz
+PAR_MINIAPPS = volta tesla maxwell joule
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -37,7 +41,8 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
 
 # Remove built-in rules
 %: %.cpp
@@ -51,10 +56,6 @@ all: $(MINIAPPS)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $@_solver.o $(COMMON_LIB) \
 	  $(MFEM_LIBS)
 
-lorentz: %: $(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)
-
 # Rules for compiling miniapp dependencies
 $(addsuffix _solver.o,$(MINIAPPS)): \
 %.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
@@ -85,7 +86,7 @@ include $(MFEM_TEST_MK)
 
 # Testing: Specific execution options
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-volta-test-par: volta-test-1 volta-test-2 volta-test-3
+volta-test-par: volta-test-1 volta-test-2
 volta-test-1: volta
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-maxit 2 -dbcs 1 -dbcg -ds '0.0 0.0 0.0 0.2 8.0')
@@ -93,29 +94,15 @@ volta-test-2: volta
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-maxit 2 -m ../../data/square-disc.mesh \
 	-dbcs '1 2 3 4 5 6 7 8' -dbcv '0 0 0 0 1 1 1 1')
-volta-test-3: volta
-	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-maxit 2 -m ../../data/inline-hex.mesh -dbcs '1 6' -dbcv '0 1')
-tesla-test-par: tesla-test-1 tesla-test-2
-tesla-test-1: tesla
+tesla-test-par: tesla
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-maxit 2 -cr '0 0 -0.2 0 0 0.2 0.2 0.4 1')
-tesla-test-2: tesla
-	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-maxit 2 -m ../../data/inline-hex.mesh -ubbc '0 0 1')
 maxwell-test-par: maxwell
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-abcs '-1' -dp '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5')
 joule-test-par: joule
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
 	-m cylinder-hex.mesh -p rod -tf 3)
-lorentz-test-par: lorentz-test-1 lorentz-test-2
-lorentz-test-1: lorentz volta-test-3
-	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-er Volta-AMR-Parallel -ec 2 -x0 '0.5 0.5 0.9' -p0 '1 0 0')
-lorentz-test-2: lorentz tesla-test-2
-	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-br Tesla-AMR-Parallel -bc 2 -x0 '0.1 0.5 0.1' -p0 '0 0.4 0.1' -tf 9)
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 
@@ -130,4 +117,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_* Lorentz*
+	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_*

miniapps/electromagnetics/tesla.cpp

@@ -253,7 +253,6 @@ int main(int argc, char *argv[])
 
    // Initialize VisIt visualization
    VisItDataCollection visit_dc("Tesla-AMR-Parallel", &pmesh);
-   visit_dc.SetFormat(DataCollection::PARALLEL_FORMAT);
 
    if ( visit )
    {

miniapps/electromagnetics/volta.cpp

@@ -266,7 +266,6 @@ int main(int argc, char *argv[])
 
    // Initialize VisIt visualization
    VisItDataCollection visit_dc("Volta-AMR-Parallel", &pmesh);
-   visit_dc.SetFormat(DataCollection::PARALLEL_FORMAT);
 
    if ( visit )
    {

miniapps/gslib/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/gslib/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -38,7 +42,8 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk

miniapps/meshing/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/meshing/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -35,7 +39,8 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk

miniapps/mtop/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/mtop/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 

miniapps/navier/navier_cht.cpp

@@ -461,7 +461,7 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 
    Kmat.Mult(u, z);
    z.Neg(); // z = -z
-   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
+   K->EliminateVDofsInRHS(ess_tdof_list, u, z);
 
    M_solver.Mult(z, du_dt);
    du_dt.Print();
@@ -483,7 +483,7 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
    MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
    Kmat.Mult(u, z);
    z.Neg();
-   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
+   K->EliminateVDofsInRHS(ess_tdof_list, u, z);
 
    T_solver.Mult(z, du_dt);
    du_dt.SetSubVector(ess_tdof_list, 0.0);

miniapps/shifted/cutH1.cpp

@@ -0,0 +1,376 @@
+//            MFEM Example for cutting an H1 space along select faces.
+//
+// Compile with: make cutH1
+//
+// Sample runs:
+//   ./cutH1 -m ../../data/star.mesh -rs 1
+//
+
+#include "mfem.hpp"
+#include "../common/mfem-common.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace common;
+
+// Used for debugging the elem-to-dof tables when the elements' attributes
+// are associated with materials. Options for lvl:
+// 0 - only duplicated materials per DOF.
+// 1 - all materials per DOF .
+// 2 - full element/material output per DOF.
+void PrintDofElemTable(const Table &elem_dof, const ParMesh &pmesh,
+                       int lvl, bool boundary)
+{
+   Table dof_elem;
+   Transpose(elem_dof, dof_elem);
+
+   const int nrows = dof_elem.Size();
+   if (boundary == false)
+   {
+      std::cout << "--- Dof-to-Elem. Total elem DOFs: " << nrows << std::endl;
+   }
+   else
+   {
+      std::cout << "--- Dof-to-Bdr. Total bndry DOFs: " << nrows << std::endl;
+   }
+
+   Array<int> dof_elements;
+   for (int dof = 0; dof < nrows; dof++)
+   {
+      // Find the materials that share the current dof.
+      std::set<int> dof_materials;
+      dof_elem.GetRow(dof, dof_elements);
+      if (lvl == 2) { std::cout << "Elements for DOF " << dof << ": \n"; }
+      for (int e = 0; e < dof_elements.Size(); e++)
+      {
+         int mat_id;
+         if (boundary == false)
+         {
+            mat_id = pmesh.GetAttribute(dof_elements[e]);
+         }
+         else
+         {
+            int face_id = pmesh.GetBdrFace(dof_elements[e]);
+            int elem_id, tmp;
+            pmesh.GetFaceElements(face_id, &elem_id, &tmp);
+            mat_id = pmesh.GetAttribute(elem_id);
+         }
+
+         if (lvl == 2) { cout << dof_elements[e] << "(" << mat_id << ") "; }
+
+         dof_materials.insert(mat_id);
+      }
+      if (lvl == 2) { std::cout << std::endl; }
+
+      if (lvl == 2) { continue; }
+      if (lvl == 0 && dof_materials.size() < 2) { continue; }
+
+      std::cout << "Materials for DOF " << dof << ": " << std::endl;
+      for (auto it = dof_materials.cbegin(); it != dof_materials.cend(); it++)
+      { std::cout << *it << ' '; }
+      std::cout << std::endl;
+   }
+   std::cout << "--- End of Table" << std::endl;
+}
+
+void VisualizeL2(ParGridFunction &gf, int size, int x, int y)
+{
+   int myid = Mpi::WorldRank();
+   int num_procs = Mpi::WorldSize();
+
+   ParMesh *pmesh = gf.ParFESpace()->GetParMesh();
+   const int order = gf.ParFESpace()->GetOrder(0);
+   L2_FECollection fec(order, pmesh->Dimension());
+   ParFiniteElementSpace pfes(pmesh, &fec);
+   ParGridFunction gf_l2(&pfes);
+   gf_l2.ProjectGridFunction(gf);
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   socketstream sol_sock(vishost, visport);
+   sol_sock << "parallel " << num_procs << " " << myid << "\n";
+   sol_sock.precision(8);
+   sol_sock << "solution\n" << *pmesh << gf_l2;
+   sol_sock << "window_geometry " << x << " " << y << " "
+            << size << " " << size << "\n"
+            << "window_title '" << "Y" << "'\n"
+            << "keys mRjlc\n" << flush;
+}
+
+void cutH1Space(ParFiniteElementSpace &pfes, bool vis, bool print)
+{
+   ParMesh &pmesh = *pfes.GetParMesh();
+   ParGridFunction x_vis(&pfes);
+
+   // Duplicate DOFs on the material interface.
+   // That is, the DOF touches different element attributes.
+   const Table &elem_dof = pfes.GetElementToDofTable(),
+                &bdre_dof = pfes.GetBdrElementToDofTable();
+   Table dof_elem, dof_bdre;
+   Table new_elem_dof(elem_dof), new_bdre_dof(bdre_dof);
+   Transpose(elem_dof, dof_elem);
+   Transpose(bdre_dof, dof_bdre);
+   const int nrows = dof_elem.Size(), n_bdr_dofs = dof_bdre.Size();
+   int ndofs = nrows;
+   Array<int> dof_elements, dof_boundaries;
+   if (print)
+   {
+      PrintDofElemTable(elem_dof, pmesh, 0, false);
+      PrintDofElemTable(bdre_dof, pmesh, 2, true);
+   }
+   for (int dof = 0; dof < nrows; dof++)
+   {
+      // Check which materials share the current dof.
+      std::set<int> dof_materials;
+      dof_elem.GetRow(dof, dof_elements);
+      for (int e = 0; e < dof_elements.Size(); e++)
+      {
+         const int mat_id = pmesh.GetAttribute(dof_elements[e]);
+         dof_materials.insert(mat_id);
+      }
+      // Count the materials for the current DOF.
+      const int dof_mat_cnt = dof_materials.size();
+
+      // Duplicate the dof if it is shared between materials.
+      if (dof_mat_cnt > 1)
+      {
+         // The material with the lowest index keeps the old DOF id.
+         // All other materials duplicate the dof.
+         auto mat = dof_materials.cbegin();
+         mat++;
+         while (mat != dof_materials.cend())
+         {
+            // Replace in all elements with material mat.
+            const int new_dof_id = ndofs;
+            for (int e = 0; e < dof_elements.Size(); e++)
+            {
+               if (pmesh.GetAttribute(dof_elements[e]) == *mat)
+               {
+                  if (print)
+                  {
+                     std::cout << "Replacing DOF (for element) : "
+                               << dof << " -> " << new_dof_id
+                               << " in EL " << dof_elements[e] << std::endl;
+                  }
+                  new_elem_dof.ReplaceConnection(dof_elements[e],
+                                                 dof, new_dof_id);
+               }
+            }
+
+            // Replace in all boundary elements with material mat.
+            int dof_bdr_cnt = 0;
+            if (dof < n_bdr_dofs)
+            {
+               dof_bdre.GetRow(dof, dof_boundaries);
+               dof_bdr_cnt = dof_boundaries.Size();
+            }
+            for (int b = 0; b < dof_bdr_cnt; b++)
+            {
+               int face_id = pmesh.GetBdrFace(dof_boundaries[b]);
+               int elem_id, tmp;
+               pmesh.GetFaceElements(face_id, &elem_id, &tmp);
+               if (pmesh.GetAttribute(elem_id) == *mat)
+               {
+                  std::cout << "Replacing DOF (for boundary): "
+                            << dof << " -> " << new_dof_id
+                            << " in BE " << dof_boundaries[b] << std::endl;
+                  new_bdre_dof.ReplaceConnection(dof_boundaries[b],
+                                                 dof, new_dof_id);
+               }
+            }
+
+            // TODO go over faces (in face_dof) that have the replaced dof (the
+            // old id), and check if they have the higher el-attributes on
+            // noth sides. For such faces, the face_dof table should be updated
+            // with the  new_dof_id.
+            // These are faces that touch the interface at a point or an edge.
+
+            ndofs++;
+            mat++;
+         }
+      }
+
+      // Used only for visualization.
+      // Must be visualized before the space update.
+      x_vis(dof) = dof_mat_cnt;
+   }
+
+   // Send the solution by socket to a GLVis server.
+   if (vis)
+   {
+      int size = 500;
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      const int myid = pfes.GetMyRank(), num_procs = pfes.GetNRanks();
+
+      socketstream sol_sock_x(vishost, visport);
+      sol_sock_x << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_x.precision(8);
+      sol_sock_x << "solution\n" << pmesh << x_vis;
+      sol_sock_x << "window_geometry " << 0 << " " << 0 << " "
+                 << size << " " << size << "\n"
+                 << "window_title '" << "X" << "'\n"
+                 << "keys mRjlc\n" << flush;
+   }
+
+   if (print)
+   {
+      PrintDofElemTable(elem_dof, pmesh, 0, false);
+      PrintDofElemTable(new_elem_dof, pmesh, 0, false);
+   }
+
+   // Remove face dofs for cut faces.
+   const Table &face_dof = pfes.GetFaceToDofTable();
+   Table new_face_dof(face_dof);
+   for (int f = 0; f < pmesh.GetNumFaces(); f++)
+   {
+      auto *ftr = pmesh.GetFaceElementTransformations(f, 3);
+      if (ftr->Elem2No > 0 &&
+          pmesh.GetAttribute(ftr->Elem1No) != pmesh.GetAttribute(ftr->Elem2No))
+      {
+         if (print)
+         {
+            std::cout << ftr->Elem1No << " " << ftr->Elem2No << std::endl;
+            std::cout << pmesh.GetAttribute(ftr->Elem1No) << " "
+                      << pmesh.GetAttribute(ftr->Elem2No) << std::endl;
+            std::cout << "Removing face dofs for face " << f << std::endl;
+         }
+         new_face_dof.RemoveRow(f);
+      }
+   }
+   new_face_dof.Finalize();
+
+   // Cut the space.
+   pfes.ReplaceElemDofTable(new_elem_dof, ndofs);
+   pfes.ReplaceBdrElemDofTable(new_bdre_dof);
+   pfes.ReplaceFaceDofTable(new_face_dof);
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   Mpi::Init(argc, argv);
+   int myid = Mpi::WorldRank();
+   int num_procs = Mpi::WorldSize();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/inline-quad.mesh";
+   int rs_levels = 0;
+   int order = 2;
+   const char *device_config = "cpu";
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   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);
+      }
+      return 1;
+   }
+   if (myid == 0) { args.PrintOptions(cout); }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // Refine the mesh.
+   Mesh mesh(mesh_file, 1, 1);
+   const int dim = mesh.Dimension();
+   for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace pfes(&pmesh, &fec);
+
+   // Assign material indices to the element attributes.
+   const int NE = pmesh.GetNE();
+   for (int i = 0; i < NE; i++)
+   {
+      Vector center;
+      pmesh.GetElementCenter(i, center);
+      Element *el = pmesh.GetElement(i);
+      if (center(0) <= 0.5 || center(1) >= 0.5)
+      {
+         el->SetAttribute(0);
+      }
+      else { el->SetAttribute(1); }
+   }
+
+   cutH1Space(pfes, true, true);
+
+   // Set face_attribute = 77 to faces that are on the material interface.
+   // Remove face dofs for cut faces.
+   for (int f = 0; f < pmesh.GetNumFaces(); f++)
+   {
+      auto *ftr = pmesh.GetFaceElementTransformations(f, 3);
+      if (ftr->Elem2No > 0 &&
+          pmesh.GetAttribute(ftr->Elem1No) != pmesh.GetAttribute(ftr->Elem2No))
+      {
+         pmesh.SetFaceAttribute(f, 77);
+      }
+   }
+
+   // Simple Dirichlet BC.
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      pfes.FiniteElementSpace::GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // RHS.
+   ParLinearForm b(&pfes);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // LHS.
+   BilinearForm a(&pfes);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   Array<int> cut_face_attributes(1);
+   cut_face_attributes[0] = 77;
+   const double sigma = -1.0, kappa = -1.0;
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa),
+                               &cut_face_attributes);
+   a.Assemble();
+
+   // Form the system.
+   ParGridFunction u(&pfes);
+   u = 0.0;
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);
+
+   // Solve.
+   GSSmoother M((SparseMatrix&)(*A));
+   //PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+   CG(*A, B, X, 1, 1000, 1e-12, 0.0);
+   a.RecoverFEMSolution(X, b, u);
+
+   VisualizeL2(u, 500, 500, 0);
+
+   const double norm = u.Norml2();
+   std::cout << "Norm: " << norm << std::endl;
+   return 0;
+}

miniapps/shifted/distance.cpp

@@ -87,8 +87,6 @@
 //
 //   Problem 4: level set: Union of doughnut and swiss cheese shapes
 //     mpirun -np 4 distance -m ../../data/inline-hex.mesh -rs 3 -o 2 -t 1.0 -p 4
-//   Problem 5: point source in mfem mesh.
-//     mpirun -np 4 distance -m ../../data/mfem.mesh -p 5 -rs 3 -t 300.0
 
 #include <fstream>
 #include <iostream>
@@ -235,8 +233,7 @@ int main(int argc, char *argv[])
                   "1: Circle / sphere level set in 2D / 3D\n\t"
                   "2: 2D sine-looking level set\n\t"
                   "3: Gyroid level set in 2D or 3D\n\t"
-                  "4: Combo of a doughnut and swiss cheese shapes in 3D.\n\t"
-                  "5: Point source in MFEM mesh.");
+                  "4: Combo of a doughnut and swiss cheese shapes in 3D.");
    args.AddOption(&rs_levels, "-rs", "--refine-serial",
                   "Number of times to refine the mesh uniformly in serial.");
    args.AddOption(&order, "-o", "--order",
@@ -302,11 +299,6 @@ int main(int argc, char *argv[])
       ls_coeff = new FunctionCoefficient(doughnut_cheese);
       smooth_steps = 0;
    }
-   else if (problem == 5)
-   {
-      ls_coeff = new DeltaCoefficient(0.0, 0.0, 1000.0);
-      smooth_steps = 0;
-   }
    else { MFEM_ABORT("Unrecognized -problem option."); }
 
    const real_t dx = AvgElementSize(pmesh);
@@ -314,7 +306,7 @@ int main(int argc, char *argv[])
    if (solver_type == 0)
    {
       auto ds = new HeatDistanceSolver(t_param * dx * dx);
-      if (problem == 0 || problem == 5)
+      if (problem == 0)
       {
          ds->transform = false;
       }
@@ -342,7 +334,7 @@ int main(int argc, char *argv[])
    // Smooth-out Gibbs oscillations from the input level set. The smoothing
    // parameter here is specified to be mesh dependent with length scale dx.
    ParGridFunction filt_gf(&pfes_s);
-   if (problem != 0 && problem != 5)
+   if (problem != 0)
    {
       real_t filter_weight = dx;
       // The normalization-based solver needs a more diffused input.

miniapps/shifted/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/shifted/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -28,9 +32,11 @@ EXTRAPOLATE_SRC = extrapolate.cpp extrapolator.cpp marking.cpp
 EXTRAPOLATE_OBJ = $(EXTRAPOLATE_SRC:.cpp=.o)
 ALGOIM_SRC =  lsf_integral.cpp
 ALGOIM_OBJ = $(ALGOIM_SRC:.cpp=.o)
+CUTH1_SRC  = cutH1.cpp
+CUTH1_OBJ = $(CUTH1_SRC:.cpp=.o)
 
 SEQ_MINIAPPS = lsf_integral
-PAR_MINIAPPS = distance diffusion extrapolate
+PAR_MINIAPPS = distance diffusion extrapolate cutH1
 
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -42,7 +48,8 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
 
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk
@@ -70,6 +77,9 @@ diffusion: $(DIFFUSION_OBJ)
 extrapolate: $(EXTRAPOLATE_OBJ)
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(EXTRAPOLATE_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
 
+cutH1: $(CUTH1_OBJ)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(CUTH1_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
+
 # Rule for building lib-common
 lib-common:
 	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
@@ -94,7 +104,7 @@ $(MFEM_LIB_FILE):
 clean: clean-build clean-exec
 
 clean-build:
-	rm -f *.o *~ distance diffusion extrapolate lsf_integral
+	rm -f *.o *~ distance diffusion extrapolate lsf_integral cutH1
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:

miniapps/spde/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/spde/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 

miniapps/tools/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/tools/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -40,7 +44,8 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
 
 all: $(MINIAPPS)
 

miniapps/toys/makefile

@@ -17,6 +17,10 @@ SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/toys/,)
 CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
    $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
 
+# Include defaults.mk to get XLINKER
+DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
+include $(DEFAULTS_MK)
+
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
@@ -37,7 +41,8 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
 
 all: $(MINIAPPS)
 

tests/unit/general/test_reduction.cpp

@@ -22,6 +22,7 @@ using namespace mfem;
 
 TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -35,7 +36,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
       int res = 0;
       mfem::reduce(
       a.Size(), res, [=] MFEM_HOST_DEVICE(int i, int &r) { r += dptr[i]; },
-      SumReducer<int> {}, use_dev);
+      SumReducer<int> {}, use_dev, workspace);
       // correct for even-length summations
       int expected = (AsConst(a)[0] + AsConst(a)[a.Size() - 1]) * a.Size() / 2;
       CAPTURE(use_dev);
@@ -45,6 +46,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 {
+   Array<long long> workspace;
    Array<long long> a(64);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -62,7 +64,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev);
+         MultReducer<long long> {}, use_dev, workspace);
          long long expected = 0;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -74,7 +76,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev);
+         MultReducer<long long> {}, use_dev, workspace);
          long long expected = 21936950640377856;
          CAPTURE(use_dev);
          REQUIRE(res == expected);
@@ -84,6 +86,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 {
+   Array<unsigned> workspace;
    Array<unsigned> a(10);
    SECTION("{ Bit unset }")
    {
@@ -105,7 +108,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev);
+         BAndReducer<unsigned> {}, use_dev, workspace);
          CAPTURE(use_dev);
          REQUIRE(res == ((~1u) & ~(1u << unset_bit)));
          REQUIRE((res & (1u << unset_bit)) == 0);
@@ -129,7 +132,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
          mfem::reduce(
             a.Size(), res,
          [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev);
+         BAndReducer<unsigned> {}, use_dev, workspace);
          CAPTURE(use_dev);
          REQUIRE(res == (1u << set_bit));
       }
@@ -138,6 +141,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 {
+   Array<unsigned> workspace;
    Array<unsigned> a(0x210);
    a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -153,7 +157,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
       mfem::reduce(
          a.Size(), res,
       [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r |= dptr[i]; },
-      BOrReducer<unsigned> {}, use_dev);
+      BOrReducer<unsigned> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == 0x3ffu);
    }
@@ -161,6 +165,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -185,7 +190,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MinReducer<int> {}, use_dev);
+      MinReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == -10);
    }
@@ -193,6 +198,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 {
+   Array<int> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -217,7 +223,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
             r = dptr[i];
          }
       },
-      MaxReducer<int> {}, use_dev);
+      MaxReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res == 999 - 10);
    }
@@ -225,6 +231,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 {
+   Array<DevicePair<int, int>> workspace;
    Array<int> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -255,7 +262,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
             r.second = dptr[i];
          }
       },
-      MinMaxReducer<int> {}, use_dev);
+      MinMaxReducer<int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second == a.Size() - 11);
@@ -264,6 +271,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 {
+   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -289,7 +297,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMinReducer<double, int> {}, use_dev);
+      ArgMinReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == -10);
       REQUIRE(res.second >= 0);
@@ -300,6 +308,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 {
+   Array<DevicePair<double, int>> workspace;
    Array<double> a(1000);
 
    auto hptr = a.HostReadWrite();
@@ -328,7 +337,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
             r.second = i;
          }
       },
-      ArgMaxReducer<double, int> {}, use_dev);
+      ArgMaxReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.first == a.Size() - 11);
       REQUIRE(res.second >= 0);
@@ -339,6 +348,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 
 TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
 {
+   Array<MinMaxLocScalar<double, int>> workspace;
    Array<double> a(1000);
    auto hptr = a.HostReadWrite();
    for (int i = 0; i < a.Size(); ++i)
@@ -373,7 +383,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
             r.max_loc = i;
          }
       },
-      ArgMinMaxReducer<double, int> {}, use_dev);
+      ArgMinMaxReducer<double, int> {}, use_dev, workspace);
       CAPTURE(use_dev);
       REQUIRE(res.min_val == -10);
       REQUIRE(res.min_loc >= 0);

tests/unit/linalg/test_matrix_dense.cpp

@@ -373,7 +373,6 @@ TEST_CASE("Batched Linear Algebra",
    const int n_rhs = 2;
 
    DenseTensor A_batch(n, n, n_mat);
-   DenseTensor A_inv_batch(n, n, n_mat);
    Vector x_batch(n * n_rhs * n_mat), y_batch(n * n_rhs * n_mat);
    std::vector<DenseMatrix> As;
    std::vector<DenseMatrix> xs, ys;
@@ -405,7 +404,6 @@ TEST_CASE("Batched Linear Algebra",
       ys.back() = 0.0;
       AddMult_a(1.5, As.back(), xs.back(), ys.back());
       A_batch(i) = As.back();
-      A_inv_batch(i) = As.back();
    }
 
    // Test batched matrix-vector products
@@ -465,33 +463,6 @@ TEST_CASE("Batched Linear Algebra",
          }
       }
    }
-
-   // Test batched matrix inverse
-   BatchedLinAlg::Get(backend).Invert(A_inv_batch);
-   A_inv_batch.HostReadWrite();
-   Vector output_col(n);
-   Vector col;
-   for (int i = 0; i < n_mat; ++i)
-   {
-      DenseMatrix Ai_inv(A_inv_batch(i));
-      for (int j = 0; j < n; ++j)
-      {
-         output_col = 0.0;
-         As[i].GetColumnReference(j, col);
-         Ai_inv.Mult(col, output_col);
-         for (int k = 0; k < n; ++k)
-         {
-            if (j == k)
-            {
-               REQUIRE(output_col(k) == MFEM_Approx(1.0));
-            }
-            else
-            {
-               REQUIRE(output_col(k) == MFEM_Approx(0.0));
-            }
-         }
-      }
-   }
 }
 
 TEST_CASE("DenseTensor copy", "[DenseMatrix][DenseTensor]")

tests/unit/linalg/test_matrix_sparse.cpp

@@ -240,54 +240,4 @@ TEST_CASE("SparseMatrix printing", "[SparseMatrix]")
    }
 }
 
-TEST_CASE("SparseMatrix cuSPARSE Bug", "[SparseMatrix][GPU]")
-{
-   // This test case ensures that we have a functioning workaround for the bug
-   // CUSPARSE-1897. In versions of cuSPARSE before 12.8, the internal buffer
-   // used for cusparseSpMV must be the same when it is called with the same
-   // matrix.
-   //
-   // By default, MFEM uses one buffer, that is shared by all sparse matrices.
-   // In the code below, a buffer is created for A, then modified for B, then
-   // used again for A. Without the workaround, this fails with cuSPARSE version
-   // earlier than 12.8 (confirmed to fail with 12.4).
-
-   const int n = 100;
-   SparseMatrix A(n, n);
-   Vector d(n);
-   d.Randomize(1);
-   for (int i = 0; i < n; ++i)
-   {
-      A.Set(i, i, d[i]);
-   }
-   A.Finalize();
-
-   Vector x(n);
-   x = 1.0;
-
-   Vector y(n);
-   A.Mult(x, y);
-
-   {
-      SparseMatrix B(20, 20);
-      for (int i = 0; i < 20; ++i)
-      {
-         for (int j = 0; j < 20; ++j)
-         {
-            B.Set(i, j, 1.0);
-         }
-      }
-      B.Finalize();
-      Vector u(20);
-      u = 1.0;
-      Vector v(20);
-      B.Mult(u, v);
-   }
-
-   A.Mult(x, y);
-
-   y -= d;
-   REQUIRE(y.Normlinf() == MFEM_Approx(0.0));
-}
-
 } // namespace mfem