Merge remote-tracking branch 'origin/master' into complex-strumpack-dev

# Conflicts:
#	examples/ex11p.cpp
#	linalg/strumpack.cpp
#	linalg/strumpack.hpp

.appveyor.yml

@@ -43,9 +43,7 @@ before_build:
 build_script:
 - cmake --build build_parallel
 - cmake --build build_serial
-- cmake --build build_serial --target exec
 
 after_build:
 # - cmake --build build_parallel --target check
-- cmake --build build_serial --target RUN_TESTS
-
+- cmake --build build_serial --target check

.travis.yml

@@ -205,7 +205,6 @@ install:
         else
             echo "Reusing cached hypre-2.10.0b/";
         fi;
-        ln -s hypre-2.10.0b hypre;
      else
         echo "Serial build, not using hypre";
      fi

CHANGELOG

@@ -8,21 +8,22 @@
                                http://mfem.org
 
 
-Version 4.0-RC2, Apr 24, 2019
+Version 4.0-RC1, Apr 11, 2019
 =============================
 
 Requirements and Limitations
 ----------------------------
 - This is a release candidate for mfem-4.0.
-- Use at your own risk -- not everything will work and the API may change.
+- Use at your own risk -- not everything will work, the API may change.
 - We are looking for feedback from friendly users.
 - Unlike previous MFEM releases, this version requires a C++11 compiler.
 
 - GPU-related limitations:
-  * Hypre preconditioners are not yet available in GPU mode.
-  * Only constant coefficients are currently supported on GPUs.
+  * NVCC is not supported in the CMake build system yet.
+  * Element batching is currently ignored.
   * Full-assembly (on device), element assembly, and matrix-free bilinear forms
-    are not supported yet. Element batching is currently ignored.
+    are not supported yet.
+  * FunctionCoefficients do not currently work on GPUs.
   * Partial assembly kernels are not implemented yet for simplices.
 
 GPU support
@@ -144,10 +145,6 @@ New and improved solvers and preconditioners
 
 Miscellaneous
 -------------
-- In SparseMatrix added the option to perform MultTranspose() by matvec with
-  computed and stored transpose matrix. This is required for deterministic
-  results when using devices such as CUDA and OpenMP.
-
 - Added unit tests based on the Catch++ library.
 
 - Renamed the option MFEM_USE_OPENMP to MFEM_USE_LEGACY_OPENMP. This legacy

CMakeLists.txt

@@ -86,13 +86,6 @@ include("${CMAKE_CURRENT_SOURCE_DIR}/config/XSDKDefaults.cmake")
 
 # Enable languages.
 enable_language(CXX)
-if (MFEM_USE_CUDA)
-   # MFEM_USE_CUDA requires CMake 3.8 or newer (for direct CUDA support)
-   cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
-   enable_language(CUDA)
-   message(STATUS "Using CUDA architecture: ${CUDA_ARCH}")
-endif()
-
 if (XSDK_ENABLE_C)
    enable_language(C)
 endif()
@@ -254,7 +247,7 @@ endif()
 
 # Axom/Sidre
 if (MFEM_USE_SIDRE)
-  find_package(Axom REQUIRED Axom)
+  find_package(Axom REQUIRED Sidre SLIC axom_utils)
 endif()
 
 # PUMI
@@ -273,32 +266,6 @@ if (MFEM_USE_PUMI)
   endif()
 endif()
 
-# CUDA
-if (MFEM_USE_CUDA)
-   set(CMAKE_CUDA_STANDARD 11)
-   set(CMAKE_CUDA_STANDARD_REQUIRED ON)
-   set(CMAKE_CUDA_EXTENSIONS OFF)
-   set(CMAKE_CUDA_FLAGS "-arch=${CUDA_ARCH} --expt-extended-lambda"
-       CACHE STRING "CUDA flags set for MFEM" FORCE)
-   if (MFEM_USE_MPI)
-      set(CUDA_CCBIN_COMPILER ${MPI_CXX_COMPILER})
-   else()
-      set(CUDA_CCBIN_COMPILER ${CMAKE_CXX_COMPILER})
-   endif()
-   string(APPEND CMAKE_CUDA_FLAGS " -ccbin ${CUDA_CCBIN_COMPILER}")
-   set(MFEM_USE_MM YES CACHE BOOL "Enable MFEM's memory manager" FORCE)
-endif()
-
-# OCCA
-if (MFEM_USE_OCCA)
-   find_package(OCCA REQUIRED)
-endif()
-
-# RAJA
-if (MFEM_USE_RAJA)
-   find_package(RAJA REQUIRED)
-endif()
-
 # MFEM_TIMER_TYPE
 if (NOT DEFINED MFEM_TIMER_TYPE)
   if (APPLE)
@@ -324,7 +291,7 @@ endif()
 #    be before SuiteSparse.
 set(MFEM_TPLS MPI_CXX OPENMP BLAS LAPACK METIS HYPRE SuiteSparse SUNDIALS PETSC
     MESQUITE SuperLUDist STRUMPACK AXOM CONDUIT GECKO GNUTLS NETCDF MPFR PUMI
-    POSIXCLOCKS MFEMBacktrace ZLIB OCCA RAJA)
+    POSIXCLOCKS MFEMBacktrace ZLIB)
 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")
@@ -360,13 +327,6 @@ set(MFEM_SOURCE_DIRS general linalg mesh fem)
 foreach(DIR IN LISTS MFEM_SOURCE_DIRS)
   add_subdirectory(${DIR})
 endforeach()
-
-if (MFEM_USE_CUDA)
-   foreach(file IN LISTS SOURCES)
-      set_property(SOURCE ${file} PROPERTY LANGUAGE CUDA)
-   endforeach()
-endif()
-
 add_subdirectory(config)
 set(MASTER_HEADERS
   ${PROJECT_SOURCE_DIR}/mfem.hpp
@@ -377,11 +337,6 @@ set(CMAKE_INSTALL_RPATH_USE_LINK_PATH ON CACHE BOOL "")
 set(CMAKE_INSTALL_RPATH "${_lib_path}" CACHE PATH "")
 set(CMAKE_INSTALL_NAME_DIR "${_lib_path}" CACHE PATH "")
 
-set(MFEM_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR} CACHE PATH
-    "The MFEM source directory" FORCE)
-set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX} CACHE PATH
-    "The MFEM install directory" FORCE)
-
 # Declaring the library
 add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
@@ -479,12 +434,12 @@ endif()
 # Add 'check' target - quick test
 if (NOT MFEM_USE_MPI)
   add_custom_target(check
-    ${CMAKE_CTEST_COMMAND} -R \"^ex1_ser\" -C ${CMAKE_CFG_INTDIR}
+    ${CMAKE_CTEST_COMMAND} -R '^ex1_ser' -C ${CMAKE_CFG_INTDIR}
     USES_TERMINAL)
   add_dependencies(check ex1)
 else()
   add_custom_target(check
-    ${CMAKE_CTEST_COMMAND} -R \"^ex1p\" -C ${CMAKE_CFG_INTDIR}
+    ${CMAKE_CTEST_COMMAND} -R '^ex1p' -C ${CMAKE_CFG_INTDIR}
     USES_TERMINAL)
   add_dependencies(check ex1p)
 endif()
@@ -529,13 +484,6 @@ install(DIRECTORY ${MFEM_SOURCE_DIRS}
   DESTINATION ${INSTALL_INCLUDE_DIR}/mfem
   FILES_MATCHING PATTERN "*.hpp")
 
-# Install the okl files
-if (MFEM_USE_OCCA)
-  install(DIRECTORY ${MFEM_SOURCE_DIRS}
-    DESTINATION ${INSTALL_INCLUDE_DIR}/mfem
-    FILES_MATCHING PATTERN "*.okl")
-endif()
-
 # Install ${HEADERS}
 # ---
 # foreach (HDR ${HEADERS})

CONTRIBUTING.md

@@ -142,16 +142,6 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
   + [`HypreParMatrix`](http://mfem.github.io/doxygen/html/classmfem_1_1HypreParMatrix.html) and [`HypreParVector`](http://mfem.github.io/doxygen/html/classmfem_1_1HypreParVector.html)
   + [`HypreSolver`](http://mfem.github.io/doxygen/html/classmfem_1_1HypreSolver.html) and other [hypre classes](http://mfem.github.io/doxygen/html/hypre_8hpp.html)
 
-- GPU and multi-core CPU support is based on device kernels supporting different
-  backends (CUDA, OCCA, RAJA, OpenMP, etc.) and an internal lightweight
-  device/host memory manager.
-
-- The main device-relevant classes and sources are:
-  + [`Device`](http://mfem.github.io/doxygen/html/device_8hpp.html)
-  + [`MemoryManager`](http://mfem.github.io/doxygen/html/mem_manager_8hpp.html)
-  + the [`MFEM_FORALL`](http://mfem.github.io/doxygen/html/forall_8hpp.html) macro
-  + the [`cuda.hpp`](http://mfem.github.io/doxygen/html/cuda_8hpp.html) and [`occa.hpp`](http://mfem.github.io/doxygen/html/occa_8hpp.html) files
-
 - The `general/` directory contains C++ classes that serve as utilities for
   communication, error handling, arrays, (Boolean) tables, timing, etc.
 

INSTALL

@@ -13,17 +13,11 @@ of MFEM is a (modern) C++ compiler, such as g++. The parallel version of MFEM
 requires an MPI C++ compiler, as well as the following external libraries:
 
 - hypre (a library of high-performance preconditioners)
-  https://github.com/hypre-space/hypre
+  http://www.llnl.gov/CASC/hypre
 
 - METIS (a family of multilevel partitioning algorithms)
   http://glaros.dtc.umn.edu/gkhome/metis/metis/overview
 
-The hypre dependency can be downloaded as a tarball from GitHub or from the
-project webpage https://www.llnl.gov/casc/hypre. For example, the 2.16.0 release
-of hypre is available at
-
-   https://github.com/hypre-space/hypre/archive/v2.16.0.tar.gz
-
 The METIS dependency can be disabled but that is not generally recommended, see
 the option MFEM_USE_METIS.
 
@@ -54,7 +48,7 @@ following package managers:
 
 - Spack, https://github.com/spack/spack
 - OpenHPC, http://openhpc.community
-- Homebrew/Science, https://github.com/Homebrew/homebrew-science (deprecated)
+- Homebrew/Science, https://github.com/Homebrew/homebrew-science
 
 We also recommend downloading and building the MFEM-based GLVis visualization
 tool which can be used to visualize the meshes and solution in MFEM's examples
@@ -66,9 +60,9 @@ Serial build:
    make serial -j 4
 
 Parallel build:
-   (download hypre and METIS 4 from above URLs)
+   (download hypre 2.10.0b and METIS 4 from above URLs)
    (build METIS 4 in ../metis-4.0 relative to mfem/)
-   (build hypre in ../hypre relative to mfem/)
+   (build hypre 2.10.0b in ../hypre-2.10.0b relative to mfem/)
    make parallel -j 4
 
 CUDA build:
@@ -93,19 +87,13 @@ Serial build:
    make -j 4   (assuming "UNIX Makefiles" generator)
 
 Parallel build:
-   (download hypre and METIS 4 from above URLs)
+   (download hypre 2.10.0b and METIS 4 from above URLs)
    (build METIS 4 in ../metis-4.0 relative to mfem/)
-   (build hypre in ../hypre relative to mfem/)
+   (build hypre 2.10.0b in ../hypre-2.10.0b relative to mfem/)
    mkdir <mfem-build-dir> ; cd <mfem-build-dir>
    cmake <mfem-source-dir> -DMFEM_USE_MPI=YES
    make -j 4
 
-CUDA build:
-   (this build requires CMake 3.8 or newer)
-   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
-   cmake <mfem-source-dir> -DMFEM_USE_CUDA=YES
-   make -j 4
-
 Example codes (serial/parallel, depending on the build):
    make examples -j 4
 
@@ -290,7 +278,6 @@ MFEM_THREAD_SAFE = YES/NO
 
 MFEM_USE_LEGACY_OPENMP = YES/NO
    Enable (basic) experimental OpenMP support. Requires MFEM_THREAD_SAFE.
-   This option is deprecated.
 
 MFEM_USE_OPENMP = YES/NO
    Enable the OpenMP backend.
@@ -406,8 +393,7 @@ MFEM_USE_PUMI = YES/NO
 
 MFEM_USE_MM = YES/NO
    Enables support for the MFEM's memory manager (MM), which is required to
-   support devices with different memory spaces. This option is required when
-   CUDA support is enabled, i.e. when MFEM_USE_CUDA=YES.
+   support devices with different memory spaces.
 
 MFEM_USE_CUDA = YES/NO
    Enables support for CUDA devices in MFEM. CUDA is a parallel computing
@@ -420,15 +406,13 @@ MFEM_USE_CUDA = YES/NO
 MFEM_USE_RAJA = YES/NO
    Enable support for the RAJA performance portability layer in MFEM. RAJA
    provides a portable abstraction for loops, supporting different programming
-   model backends. When using RAJA built with CUDA support, CUDA support must be
-   also enabled in MFEM, i.e. MFEM_USE_CUDA=YES must be set.
+   model backends. When using the RAJA CUDA backend, MFEM_USE_MM is required.
 
 MFEM_USE_OCCA = YES/NO
    Enables support for the OCCA library in MFEM. OCCA is an open-source library
    which aims to make it easy to program different types of devices (e.g. CPU,
    GPU, FPGA) by providing an unified API for interacting with JIT-compiled
-   backends. In order to use the OCCA CUDA backend, CUDA support must be enabled
-   in MFEM as well, i.e. MFEM_USE_CUDA=YES must be set.
+   backends. When using the OCCA CUDA backend, MFEM_USE_MM is required.
 
 MFEM_BUILD_TAG = (any value)
    An optional tag to characterize the build. Exported to config/config.mk.
@@ -451,7 +435,7 @@ directory and use the string @MFEM_DIR@, e.g. HYPRE_OPT = -I@MFEM_DIR@/../hypre.
 The specific libraries and their options are:
 
 - HYPRE, required for the parallel build, i.e. when MFEM_USE_MPI = YES.
-  URL: https://github.com/hypre-space/hypre and https://www.llnl.gov/casc/hypre
+  URL: http://www.llnl.gov/CASC/hypre
   Options: HYPRE_OPT, HYPRE_LIB.
 
 - METIS, used when MFEM_USE_METIS = YES. If using METIS 5, set
@@ -661,8 +645,6 @@ Configuration variables (CMake)
 ===============================
 See the configuration file config/defaults.cmake for the default settings.
 
-Note: the option MFEM_USE_CUDA requires CMake version 3.8 or newer!
-
 Non-standard CMake variables for compilers:
    CXX    - If set, overwrite the auto-detected C++ compiler, serial build
    MPICXX - If set, overwrite the auto-detected MPI C++ compiler, parallel build
@@ -693,6 +675,9 @@ MFEM_USE_NETCDF
 MFEM_USE_MPFR
 MFEM_USE_GZSTREAM
 MFEM_USE_PUMI
+
+The following GNU make options are not supported with CMake yet:
+
 MFEM_USE_CUDA
 MFEM_USE_OCCA
 MFEM_USE_RAJA
@@ -743,8 +728,6 @@ The CMake build system adds auto-detection for the following packages/libraries:
  - LIBUNWIND
  - POSIXCLOCKS
  - PUMI
- - OCCA
- - RAJA
 
 The following built-in CMake packages are also used:
 

README

@@ -8,9 +8,9 @@
                                http://mfem.org
 
 MFEM is a modular parallel C++ library for finite element methods. Its goal is
-to enable high-performance scalable finite element discretization research and
-application development on a wide variety of platforms, ranging from laptops to
-supercomputers.
+to enable the research and development of scalable finite element discretization
+and solver algorithms through general finite element abstractions, accurate and
+flexible visualization, and tight integration with the hypre library.
 
 * For building instructions, see the file INSTALL, or type "make help".
 
@@ -39,24 +39,23 @@ conforming and non-conforming (AMR) adaptive refinement. Arbitrary element
 transformations, allowing for high-order mesh elements with curved boundaries,
 are also supported.
 
-When used as a "finite element to linear algebra translator", MFEM can take a
-problem described in terms of finite element-type objects, and produce the
-corresponding linear algebra vectors and fully or partially assembled operators,
-e.g. in the form of global sparse matrices or matrix-free operators. The library
-includes simple smoothers and Krylov solvers, such as PCG, MINRES and GMRES, as
-well as support for sequential sparse direct solvers from the SuiteSparse
+MFEM is commonly used as a "finite element to linear algebra translator", since
+it can take a problem described in terms of finite element-type objects, and
+produce the corresponding linear algebra vectors and sparse matrices. In order
+to facilitate this, MFEM uses compressed sparse row (CSR) sparse matrix storage
+and includes simple smoothers and Krylov solvers, such as PCG, MINRES and GMRES,
+as well as support for sequential sparse direct solvers from the SuiteSparse
 library. Nonlinear solvers (the Newton method), eigensolvers (LOBPCG), and
 several explicit and implicit Runge-Kutta time integrators are also available.
 
 MFEM supports MPI-based parallelism throughout the library, and can readily be
-used as a scalable unstructured finite element problem generator. As of version
-4.0, MFEM offers initial support for GPU acceleration, and programming models,
-such as CUDA, OCCA, RAJA and OpenMP. MFEM-based applications require minimal
-changes to switch from a serial to a high-performing MPI-parallel version of the
-code, where they can take advantage of the integrated linear solvers from the
-hypre library.  Comprehensive support for other external packages, e.g. PETSc
-and SUNDIALS is also included, giving access to many additional linear and
-nonlinear solvers, preconditioners, time integrators, etc.
+used as a scalable unstructured finite element problem generator. MFEM-based
+applications require minimal changes to transition from a serial to a
+high-performing parallel version of the code, where they can take advantage of
+the integrated scalable linear solvers from the hypre library. Comprehensive
+support for other external packages, e.g. PETSc and SUNDIALS is also included,
+giving access to many additional linear and nonlinear solvers, preconditioners,
+time integrators, etc.
 
 For examples of using MFEM, see the examples/ and miniapps/ directories, as well
 as the OpenGL visualization tool GLVis which is available at http://glvis.org.

config/XSDKDefaults.cmake

@@ -74,7 +74,7 @@
 
 IF (NOT COMMAND PRINT_VAR)
   FUNCTION(PRINT_VAR  VAR_NAME)
-    MESSAGE(STATUS "${VAR_NAME} = '${${VAR_NAME}}'")
+    MESSAGE("-- " "${VAR_NAME} = '${${VAR_NAME}}'")
   ENDFUNCTION()
 ENDIF()
 
@@ -166,14 +166,14 @@ IF (USE_XSDK_DEFAULTS)
     ENDIF()
     XSDK_HANDLE_LANG_DEFAULTS(Fortran  FC  "FFLAGS;FCFLAGS")
   ENDIF()
-
+  
   # Set XSDK defaults for other CMake variables
-
+  
   IF ("${BUILD_SHARED_LIBS}"  STREQUAL  "")
     MESSAGE("-- " "XSDK: Setting default BUILD_SHARED_LIBS=TRUE")
     SET(BUILD_SHARED_LIBS  TRUE  CACHE  BOOL  "Set by default in XSDK mode")
   ENDIF()
-
+  
   IF ("${CMAKE_BUILD_TYPE}"  STREQUAL  "")
     MESSAGE("-- " "XSDK: Setting default CMAKE_BUILD_TYPE=DEBUG")
     SET(CMAKE_BUILD_TYPE  DEBUG  CACHE  STRING  "Set by default in XSDK mode")
@@ -181,13 +181,6 @@ IF (USE_XSDK_DEFAULTS)
 
 ENDIF()
 
-
-##################################################################################
-#
-#               MFEM-specific additions: set TPL MFEM_USE_* defaults
-#
-##################################################################################
-
 IF (DEFINED TPL_ENABLE_MPI)
   SET(MFEM_USE_MPI ${TPL_ENABLE_MPI} CACHE BOOL "Enable MPI parallel build" FORCE)
 ENDIF()
@@ -259,15 +252,3 @@ ENDIF()
 IF (DEFINED TPL_ENABLE_PUMI)
   SET(MFEM_USE_PUMI ${TPL_ENABLE_PUMI} CACHE BOOL "Enable PUMI" FORCE)
 ENDIF()
-
-IF (DEFINED TPL_ENABLE_CUDA)
-  SET(MFEM_USE_CUDA ${TPL_ENABLE_CUDA} CACHE BOOL "Enable CUDA" FORCE)
-ENDIF()
-
-IF (DEFINED TPL_ENABLE_OCCA)
-  SET(MFEM_USE_OCCA ${TPL_ENABLE_OCCA} CACHE BOOL "Enable OCCA" FORCE)
-ENDIF()
-
-IF (DEFINED TPL_ENABLE_RAJA)
-  SET(MFEM_USE_RAJA ${TPL_ENABLE_RAJA} CACHE BOOL "Enable RAJA" FORCE)
-ENDIF()

config/cmake/MFEMConfig.cmake.in

@@ -41,10 +41,6 @@ set(MFEM_USE_MPFR @MFEM_USE_MPFR@)
 set(MFEM_USE_SIDRE @MFEM_USE_SIDRE@)
 set(MFEM_USE_CONDUIT @MFEM_USE_CONDUIT@)
 set(MFEM_USE_PUMI @MFEM_USE_PUMI@)
-set(MFEM_USE_MM @MFEM_USE_MM@)
-set(MFEM_USE_CUDA @MFEM_USE_CUDA@)
-set(MFEM_USE_OCCA @MFEM_USE_OCCA@)
-set(MFEM_USE_RAJA @MFEM_USE_RAJA@)
 
 set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
 set(MFEM_CXX_FLAGS "@CMAKE_CXX_FLAGS@")

config/cmake/config.hpp.in

@@ -30,12 +30,6 @@
 #define MFEM_VERSION_MINOR (((MFEM_VERSION)/100)%100)
 #define MFEM_VERSION_PATCH ((MFEM_VERSION)%100)
 
-// MFEM source directory.
-#define MFEM_SOURCE_DIR "@MFEM_SOURCE_DIR@"
-
-// MFEM install directory.
-#define MFEM_INSTALL_DIR "@MFEM_INSTALL_DIR@"
-
 // Description of the git commit used to build MFEM.
 #cmakedefine MFEM_GIT_STRING "@MFEM_GIT_STRING@"
 
@@ -110,19 +104,6 @@
 // Enable MFEM functionality based on the PUMI library
 #cmakedefine MFEM_USE_PUMI
 
-// Build the GPU/CUDA-enabled version of the MFEM library.
-// Requires a CUDA compiler (nvcc).
-#cmakedefine MFEM_USE_CUDA
-
-// Enable MFEM functionality based on the RAJA library
-#cmakedefine MFEM_USE_RAJA
-
-// Enable MFEM functionality based on the OCCA library
-#cmakedefine MFEM_USE_OCCA
-
-// Enable MFEM's internal Memory Manager (needed e.g. for MFEM_USE_CUDA)
-#cmakedefine MFEM_USE_MM
-
 // Which library functions to use in class StopWatch for measuring time.
 // For a list of the available options, see INSTALL.
 // If not defined, an option is selected automatically.

config/cmake/modules/FindAxom.cmake

@@ -18,4 +18,6 @@ include(MfemCmakeUtilities)
 # Note: components are enabled based on the find_package() parameters.
 mfem_find_package(Axom AXOM AXOM_DIR "include" "" "lib" ""
   "Paths to headers required by Axom." "Libraries required by Axom."
-  ADD_COMPONENT Axom "include" axom/config.hpp "lib" axom)
+  ADD_COMPONENT Sidre "include" sidre/sidre.hpp "lib" sidre
+  ADD_COMPONENT SLIC "include" slic/slic.hpp "lib" slic
+  ADD_COMPONENT axom_utils "include" axom_utils/Utilities.hpp "lib" axom_utils)

config/cmake/modules/FindOCCA.cmake

@@ -1,19 +0,0 @@
-# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
-# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
-# See file COPYRIGHT for details.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability see http://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the GNU Lesser General Public License (as published by the Free
-# Software Foundation) version 2.1 dated February 1999.
-
-# Defines the following variables:
-#   - OCCA_FOUND
-#   - OCCA_LIBRARIES
-#   - OCCA_INCLUDE_DIRS
-
-include(MfemCmakeUtilities)
-mfem_find_package(OCCA OCCA OCCA_DIR "include" "occa.hpp" "lib" "occa"
-  "Paths to headers required by OCCA." "Libraries required by OCCA.")

config/cmake/modules/FindRAJA.cmake

@@ -1,30 +0,0 @@
-# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
-# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
-# See file COPYRIGHT for details.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability see http://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the GNU Lesser General Public License (as published by the Free
-# Software Foundation) version 2.1 dated February 1999.
-
-# Defines the following variables:
-#   - RAJA_FOUND
-#   - RAJA_LIBRARIES
-#   - RAJA_INCLUDE_DIRS
-
-include(MfemCmakeUtilities)
-mfem_find_package(RAJA RAJA RAJA_DIR "include" "RAJA/RAJA.hpp" "lib" "RAJA"
-  "Paths to headers required by RAJA." "Libraries required by RAJA.")
-
-if (NOT RAJA_CONFIG_CMAKE)
-   set(RAJA_CONFIG_CMAKE "${RAJA_DIR}/share/raja/cmake/raja-config.cmake")
-endif()
-if (EXISTS "${RAJA_CONFIG_CMAKE}")
-   include("${RAJA_CONFIG_CMAKE}")
-   if (ENABLE_CUDA AND NOT MFEM_USE_CUDA)
-      message(FATAL_ERROR
-         "RAJA is built with CUDA: MFEM_USE_CUDA=YES is required")
-   endif()
-endif()

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -232,12 +232,10 @@ function(mfem_find_package Name Prefix DirVar IncSuffixes Header LibSuffixes
   # If we have the TPL_ versions of _INCLUDE_DIRS and _LIBRARIES then set the
   # standard ${Prefix} versions
   if (TPL_${Prefix}_INCLUDE_DIRS)
-    set(${Prefix}_INCLUDE_DIRS ${TPL_${Prefix}_INCLUDE_DIRS} CACHE STRING
-      "TPL_${Prefix}_INCLUDE_DIRS was found." FORCE)
+    set(${Prefix}_INCLUDE_DIRS ${TPL_${Prefix}_INCLUDE_DIRS} CACHE STRING "TPL_${Prefix}_INCLUDE_DIRS was found." FORCE)
   endif()
   if (TPL_${Prefix}_LIBRARIES)
-    set(${Prefix}_LIBRARIES ${TPL_${Prefix}_LIBRARIES} CACHE STRING
-      "TPL_${Prefix}_LIBRARIES was found." FORCE)
+    set(${Prefix}_LIBRARIES ${TPL_${Prefix}_LIBRARIES} CACHE STRING "TPL_${Prefix}_LIBRARIES was found." FORCE)
   endif()
 
   # Quick return
@@ -720,8 +718,7 @@ function(mfem_export_mk_files)
       MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
       MFEM_USE_SUPERLU MFEM_USE_STRUMPACK MFEM_USE_GECKO MFEM_USE_GNUTLS
       MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_MPFR MFEM_USE_SIDRE
-      MFEM_USE_CONDUIT MFEM_USE_PUMI MFEM_USE_MM MFEM_USE_CUDA MFEM_USE_OCCA
-      MFEM_USE_RAJA)
+      MFEM_USE_CONDUIT MFEM_USE_PUMI)
   foreach(var ${CONFIG_MK_BOOL_VARS})
     if (${var})
       set(${var} YES)
@@ -729,7 +726,6 @@ function(mfem_export_mk_files)
       set(${var} NO)
     endif()
   endforeach()
-  # TODO: Add support for MFEM_USE_CUDA=YES
   set(MFEM_CXX ${CMAKE_CXX_COMPILER})
   set(MFEM_CPPFLAGS "")
   string(STRIP "${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"

config/config.hpp

@@ -56,9 +56,4 @@
 #endif
 #endif // MFEM_USE_MPI not defined
 
-// CUDA requires the memory manager
-#if defined(MFEM_USE_CUDA) && !defined(MFEM_USE_MM)
-#error Building with CUDA (MFEM_USE_CUDA=YES) requires MFEM_USE_MM=YES
-#endif
-
 #endif // MFEM_CONFIG_HPP

config/defaults.cmake

@@ -42,10 +42,6 @@ option(MFEM_USE_MPFR "Enable MPFR usage." OFF)
 option(MFEM_USE_SIDRE "Enable Axom/Sidre usage" OFF)
 option(MFEM_USE_CONDUIT "Enable Conduit usage" OFF)
 option(MFEM_USE_PUMI "Enable PUMI" OFF)
-option(MFEM_USE_MM "Enable MFEM's memory manager" OFF)
-option(MFEM_USE_CUDA "Enable CUDA" OFF)
-option(MFEM_USE_OCCA "Enable OCCA" OFF)
-option(MFEM_USE_RAJA "Enable RAJA" OFF)
 
 set(MFEM_MPI_NP 4 CACHE STRING "Number of processes used for MPI tests")
 
@@ -63,16 +59,13 @@ option(MFEM_ENABLE_MINIAPPS "Build all of the miniapps" OFF)
 # set(CXX g++)
 # set(MPICXX mpicxx)
 
-# Set the target CUDA architecture
-set(CUDA_ARCH "sm_60" CACHE STRING "Target CUDA architecture.")
-
 set(MFEM_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 
 # The *_DIR paths below will be the first place searched for the corresponding
 # headers and library. If these fail, then standard cmake search is performed.
 # Note: if the variables are already in the cache, they are not overwritten.
 
-set(HYPRE_DIR "${MFEM_DIR}/../hypre/src/hypre" CACHE PATH
+set(HYPRE_DIR "${MFEM_DIR}/../hypre-2.10.0b/src/hypre" CACHE PATH
     "Path to the hypre library.")
 # If hypre was compiled to depend on BLAS and LAPACK:
 # set(HYPRE_REQUIRED_PACKAGES "BLAS" "LAPACK" CACHE STRING
@@ -161,9 +154,6 @@ set(Axom_REQUIRED_PACKAGES "Conduit/relay" CACHE STRING
 set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
     "Directory where PUMI is installed")
 
-set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
-set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
-
 set(BLAS_INCLUDE_DIRS "" CACHE STRING "Path to BLAS headers.")
 set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")
 set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")

config/defaults.mk

@@ -136,7 +136,7 @@ LIBUNWIND_OPT = -g
 LIBUNWIND_LIB = $(if $(NOTMAC),-lunwind -ldl,)
 
 # HYPRE library configuration (needed to build the parallel version)
-HYPRE_DIR = @MFEM_DIR@/../hypre/src/hypre
+HYPRE_DIR = @MFEM_DIR@/../hypre-2.10.0b/src/hypre
 HYPRE_OPT = -I$(HYPRE_DIR)/include
 HYPRE_LIB = -L$(HYPRE_DIR)/lib -lHYPRE
 
@@ -291,7 +291,7 @@ SIDRE_LIB = \
    -Wl,-rpath,$(SIDRE_DIR)/lib -L$(SIDRE_DIR)/lib \
    -Wl,-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
    -Wl,-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
-   -laxom -lconduit -lconduit_relay -lhdf5 $(ZLIB_LIB) -ldl
+   -lsidre -lslic -laxom_utils -lconduit -lconduit_relay -lhdf5 $(ZLIB_LIB) -ldl
 
 # PUMI
 # Note that PUMI_DIR is needed -- it is used to check for gmi_sim.h
@@ -300,17 +300,19 @@ PUMI_OPT = -I$(PUMI_DIR)/include
 PUMI_LIB = -L$(PUMI_DIR)/lib -lpumi -lcrv -lma -lmds -lapf -lpcu -lgmi -lparma\
    -llion -lmth -lapf_zoltan -lspr
 
-# CUDA library configuration (currently not needed)
+# CUDA library configuration. Since we compile and link with nvcc (when CUDA is
+# enabled) we only need to explicitly link with the CUDA driver, libcuda.*,
+# which is usually in a system path.
 CUDA_OPT =
-CUDA_LIB =
+CUDA_LIB = $(if $(NOTMAC),,-L/usr/local/cuda/lib) -lcuda
 
 # OCCA library configuration
-OCCA_DIR = @MFEM_DIR@/../occa
+OCCA_DIR ?= @MFEM_DIR@/../occa
 OCCA_OPT = -I$(OCCA_DIR)/include
 OCCA_LIB = $(XLINKER)-rpath,$(OCCA_DIR)/lib -L$(OCCA_DIR)/lib -locca
 
 # RAJA library configuration
-RAJA_DIR = @MFEM_DIR@/../raja
+RAJA_DIR ?= @MFEM_DIR@/../raja
 RAJA_OPT = -I$(RAJA_DIR)/include
 ifdef CUB_DIR
    RAJA_OPT += -I$(CUB_DIR)

config/sample-runs.sh

@@ -18,8 +18,6 @@ run_prefix=""
 run_vg="valgrind --leak-check=full --show-reachable=yes --track-origins=yes"
 run_suffix="-no-vis"
 skip_gen_meshes="yes"
-# filter-out device runs ("no") or non-device runs ("yes"):
-device_runs="no"
 cur_dir="${PWD}"
 mfem_dir="$(cd "$(dirname "$0")"/.. && pwd)"
 mfem_build_dir=""
@@ -150,11 +148,6 @@ function extract_sample_runs()
    if [ "$skip_gen_meshes" == "yes" ]; then
       runs=`printf "%s" "$runs" | grep -v ".* -m .*\.gen"`
    fi
-   if [ "$device_runs" == "yes" ]; then
-      runs=`printf "%s" "$runs" | grep ".* -d .*"`
-   else
-      runs=`printf "%s" "$runs" | grep -v ".* -d .*"`
-   fi
    IFS=$'\n'
    runs=(${runs})
    IFS="${old_IFS}"
@@ -176,9 +169,6 @@ function help_message()
       -g <dir> <pattern>
                   Specify explicitly a group (dir + file pattern) to run; This
                   option can be used multiple times to define multiple groups
-      -dev        configure only sample runs using devices.
-                  To test with a parallel build, the parallel (-p|-par) option
-                  should be set first on the command line.
       -v          Enable valgrind
       -o <dir>    [${output_dir:-"<empty>: output goes to stdout"}]
                   If not empty, save output to files inside <dir>
@@ -263,7 +253,7 @@ case "$1" in
    -h|-help)
       opt_help="yes"
       ;;
-   -p|-par)
+   -p|-parallel)
       mfem_config="MFEM_USE_MPI=YES MFEM_DEBUG=NO"
       ;;
    -g)
@@ -274,11 +264,6 @@ case "$1" in
       groups=("${groups[@]}" "${test_group}")
       shift 2
       ;;
-   -dev)
-       device_runs="yes"
-       mfem_config+=" MFEM_USE_CUDA=YES MFEM_USE_MM=YES \
-MFEM_USE_OCCA=YES MFEM_USE_RAJA=YES MFEM_USE_OPENMP=YES"
-      ;;
    -v)
       valgrind="yes"
       ;;
@@ -309,10 +294,6 @@ MFEM_USE_OCCA=YES MFEM_USE_RAJA=YES MFEM_USE_OPENMP=YES"
    -n)
       run_prefix="echo"
       ;;
-   -*)
-      echo "unknown option: '$1'"
-      exit 1
-      ;;
    *=*)
       eval $1
       ;;

doc/CodeDocumentation.dox

@@ -35,10 +35,6 @@ namespace mfem {
  * - HypreParMatrix and HypreParVector
  * - HypreSolver and other \link hypre.hpp hypre classes\endlink
  *
- * <H3>Main GPU classes</H3>
- * - Device
- * - MemoryManager
- *
  * <H3>Example codes</H3>
  * - <a class="el" href="examples_2ex1_8cpp_source.html">Example 1</a>: nodal H1 FEM for the Laplace problem
  * - <a class="el" href="examples_2ex1p_8cpp_source.html">Example 1p</a>: parallel nodal H1 FEM for the Laplace problem

examples/CMakeLists.txt

@@ -27,6 +27,7 @@ list(APPEND ALL_EXE_SRCS
   ex18.cpp
   ex19.cpp
   ex20.cpp
+  ex21.cpp
   ex22.cpp
   )
 
@@ -52,6 +53,7 @@ if (MFEM_USE_MPI)
     ex18p.cpp
     ex19p.cpp
     ex20p.cpp
+    ex21p.cpp
     ex22p.cpp
     )
 endif()

examples/ex1.cpp

@@ -26,12 +26,12 @@
 //               ex1 -m ../data/mobius-strip.mesh -o -1 -sc
 //
 // Device sample runs:
-//               ex1 -pa -d cuda
-//               ex1 -pa -d raja-cuda
-//               ex1 -pa -d occa-cuda
-//               ex1 -pa -d raja-omp
-//               ex1 -pa -d occa-omp
-//               ex1 -m ../data/beam-hex.mesh -pa -d cuda
+//             > ex1 -pa -d cuda
+//             > ex1 -pa -d raja-cuda
+//             > ex1 -pa -d occa-cuda
+//             > ex1 -pa -d raja-omp
+//             > ex1 -pa -d occa-omp
+//             > ex1 -m ../data/beam-hex.mesh -pa -d cuda
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Laplace problem

examples/ex1p.cpp

@@ -26,9 +26,9 @@
 //               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
 //
 // Device sample runs:
-//               mpirun -np 4 ex1p -pa -d cuda
-//               mpirun -np 4 ex1p -pa -d occa-cuda
-//               mpirun -np 4 ex1p -pa -d raja-omp
+//             > mpirun -np 4 ex1p -pa -d cuda
+//             > mpirun -np 4 ex1p -pa -d occa-cuda
+//             > mpirun -np 4 ex1p -pa -d raja-omp
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Laplace problem

examples/ex21.cpp

@@ -0,0 +1,477 @@
+//                               MFEM Example 21
+//
+// Compile with: make ex21
+//
+// Sample runs:  ex21 -m ../data/inline-segment.mesh -o 3
+//               ex21 -m ../data/inline-tri.mesh -o 3
+//               ex21 -m ../data/inline-quad.mesh -o 3
+//               ex21 -m ../data/inline-quad.mesh -o 3 -p 1
+//               ex21 -m ../data/inline-quad.mesh -o 3 -p 2
+//               ex21 -m ../data/inline-tet.mesh -o 2
+//               ex21 -m ../data/inline-hex.mesh -o 2
+//               ex21 -m ../data/inline-hex.mesh -o 2 -p 1
+//               ex21 -m ../data/inline-hex.mesh -o 2 -p 2
+//               ex21 -m ../data/star.mesh -o 2 -sigma 10.0
+//
+// Description:  This example code demonstrates the use of MFEM to define and
+//               solve simple complex-valued linear systems.  We implement three
+//               variants of a damped harmonic oscillator:
+//
+//               1) A scalar H1 field
+//                  -Div(a Grad u) - omega^2 b u + i omega c u = 0
+//
+//               2) A vector H(Curl) field
+//                  Curl(a Curl u) - omega^2 b u + i omega c u = 0
+//
+//               3) A vector H(Div) field
+//                  -Grad(a Div u) - omega^2 b u + i omega c u = 0
+//
+//               In each case the field is driven by a forced oscillation, with
+//               angular frequency omega, imposed at the boundary or a portion
+//               of the boundary.
+//
+//               In electromagnetics the coefficients are typically named the
+//               permeability, mu = 1/a, permittivity, epsilon = b, and
+//               conductivity, sigma = c.  The user can specify these constants
+//               using either set of names.
+//
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+static double mu_ = 1.0;
+static double epsilon_ = 1.0;
+static double sigma_ = 20.0;
+static double omega_ = 10.0;
+
+double u0_real_exact(const Vector &);
+double u0_imag_exact(const Vector &);
+
+void u1_real_exact(const Vector &, Vector &);
+void u1_imag_exact(const Vector &, Vector &);
+
+void u2_real_exact(const Vector &, Vector &);
+void u2_imag_exact(const Vector &, Vector &);
+
+bool check_for_inline_mesh(const char * mesh_file);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int ref_levels = 0;
+   int order = 1;
+   int prob = 0;
+   double freq = -1.0;
+   double a_coef = 0.0;
+   bool visualization = 1;
+   bool herm_conv = true;
+   bool exact_sol = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&prob, "-p", "--problem-type",
+                  "Choose from 0: H_1, 1: H(Curl), or 2: H(Div) "
+                  "damped harmonic oscillator.");
+   args.AddOption(&a_coef, "-a", "--stiffness-coef",
+                  "Stiffness coefficient (spring constant or 1/mu).");
+   args.AddOption(&epsilon_, "-b", "--mass-coef",
+                  "Mass coefficient (or epsilon).");
+   args.AddOption(&sigma_, "-c", "--damping-coef",
+                  "Damping coefficient (or sigma).");
+   args.AddOption(&mu_, "-mu", "--permeability",
+                  "Permeability of free space (or 1/(spring constant)).");
+   args.AddOption(&epsilon_, "-eps", "--permittivity",
+                  "Permittivity of free space (or mass constant).");
+   args.AddOption(&sigma_, "-sigma", "--conductivity",
+                  "Conductivity (or damping constant).");
+   args.AddOption(&freq, "-f", "--frequency",
+                  "Frequency (in Hz).");
+   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
+                  "--no-hermitian", "Use convention for Hermitian operators.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   if ( a_coef != 0.0 )
+   {
+      mu_ = 1.0 / a_coef;
+   }
+   if ( freq > 0.0 )
+   {
+      omega_ = 2.0 * M_PI * freq;
+   }
+
+   exact_sol = check_for_inline_mesh(mesh_file);
+   if (exact_sol)
+   {
+      cout << "Identified an 'inline' mesh" << endl;
+   }
+
+   ComplexOperator::Convention conv =
+      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
+
+   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes
+   //    with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase resolution. In this example we do
+   //    'ref_levels' of uniform refinement where the user specifies
+   //    the number of levels with the '-r' option.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use continuous
+   //    Lagrange, Nedelec, or Raviart-Thomas finite elements of the specified
+   //    order.
+   if (dim == 1 && prob != 0 )
+   {
+      cout << "Switching to problem type 0, H1 basis functions, "
+           << "for 1 dimensional mesh." << endl;
+      prob = 0;
+   }
+
+   FiniteElementCollection *fec;
+   switch (prob)
+   {
+      case 0:  fec = new H1_FECollection(order, dim);      break;
+      case 1:  fec = new ND_FECollection(order, dim);      break;
+      case 2:  fec = new RT_FECollection(order - 1, dim);  break;
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
+        << endl;
+
+   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   //    In this example, the boundary conditions are defined based on the type
+   //    of mesh and the problem type.
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   if (mesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      if (exact_sol)
+      {
+         switch (prob)
+         {
+            case 0:   ess_bdr = 0; ess_bdr[0] = 1;  break;
+            default:  ess_bdr = 1; ess_bdr[2] = 0;  break;
+         }
+      }
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 6. Set up the linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system.
+   ComplexLinearForm b(fespace, conv);
+   b.Vector::operator=(0.0);
+
+   // 7. Define the solution vector u as a finite element grid function
+   //    corresponding to fespace. Initialize u with initial guess of 1+0i
+   //    or the exact solution if it is known.
+   ComplexGridFunction u(fespace);
+   ComplexGridFunction * u_exact = NULL;
+   if (exact_sol) { u_exact = new ComplexGridFunction(fespace); }
+
+   FunctionCoefficient u0_r(u0_real_exact);
+   FunctionCoefficient u0_i(u0_imag_exact);
+   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
+   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
+   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
+   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
+
+   ConstantCoefficient zeroCoef(0.0);
+   ConstantCoefficient oneCoef(1.0);
+
+   Vector zeroVec(dim); zeroVec = 0.0;
+   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
+   VectorConstantCoefficient zeroVecCoef(zeroVec);
+   VectorConstantCoefficient oneVecCoef(oneVec);
+
+   switch (prob)
+   {
+      case 0:
+         u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u0_r, u0_i); }
+         break;
+      case 1:
+         u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u1_r, u1_i); }
+         break;
+      case 2:
+         u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u2_r, u2_i); }
+         break;
+   }
+
+   if (visualization && exact_sol)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *mesh << u_exact->real()
+                 << "window_title 'Exact Real Part'" << flush;
+      sol_sock_i << "solution\n" << *mesh << u_exact->imag()
+                 << "window_title 'Exact Imaginary Part'" << flush;
+   }
+
+   // 8. Set up the sesquilinear form a(.,.) on the finite element
+   //    space corresponding to the damped harmonic oscillator operator
+   //    of the appropriate type:
+   //
+   //    0) A scalar H1 field
+   //       -Div(a Grad) - omega^2 b + i omega c
+   //
+   //    1) A vector H(Curl) field
+   //       Curl(a Curl) - omega^2 b + i omega c
+   //
+   //    2) A vector H(Div) field
+   //       -Grad(a Div) - omega^2 b + i omega c
+   //
+   ConstantCoefficient stiffnessCoef(1.0/mu_);
+   ConstantCoefficient massCoef(-omega_ * omega_ * epsilon_);
+   ConstantCoefficient lossCoef(omega_ * sigma_);
+   ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
+
+   SesquilinearForm *a = new SesquilinearForm(fespace, conv);
+   switch (prob)
+   {
+      case 0:
+         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new MassIntegrator(massCoef),
+                                new MassIntegrator(lossCoef));
+         break;
+      case 1:
+         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
+         break;
+      case 2:
+         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
+         break;
+   }
+
+   // 9. Assemble the bilinear form and the corresponding linear
+   //    system, applying any necessary transformations such as:
+   //    assembly, eliminating boundary conditions, applying conforming
+   //    constraints for non-conforming AMR, etc.
+   a->Assemble();
+
+   OperatorHandle A;
+   Vector B, U;
+
+   a->FormLinearSystem(ess_tdof_list, u, b, A, U, B);
+   u = 0.0;
+   U = 0.0;
+
+   {
+      ComplexSparseMatrix * Asp =
+         dynamic_cast<ComplexSparseMatrix*>(A.Ptr());
+
+      cout << "Size of linear system: "
+           << 2 * Asp->real().Width() << endl << endl;
+   }
+
+   // 10. Define and apply a GMRES solver for AU=B.
+   {
+      GMRESSolver gmres;
+      gmres.SetOperator(*A.Ptr());
+      gmres.SetRelTol(1e-12);
+      gmres.SetMaxIter(1000);
+      gmres.SetPrintLevel(1);
+      gmres.Mult(B, U);
+   }
+
+   // 11. Recover the solution as a finite element grid function and
+   //     compute the errors if the exact solution is known.
+   a->RecoverFEMSolution(U, b, u);
+
+   if (exact_sol)
+   {
+      double err_r = -1.0;
+      double err_i = -1.0;
+
+      switch (prob)
+      {
+         case 0:
+            err_r = u.real().ComputeL2Error(u0_r);
+            err_i = u.imag().ComputeL2Error(u0_i);
+            break;
+         case 1:
+            err_r = u.real().ComputeL2Error(u1_r);
+            err_i = u.imag().ComputeL2Error(u1_i);
+            break;
+         case 2:
+            err_r = u.real().ComputeL2Error(u2_r);
+            err_i = u.imag().ComputeL2Error(u2_i);
+            break;
+      }
+
+      cout << endl;
+      cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
+      cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
+      cout << endl;
+   }
+
+   // 12. Save the refined mesh and the solution. This output can be
+   //     viewed later using GLVis: "glvis -m mesh -g sol".
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+
+      ofstream sol_r_ofs("sol_r.gf");
+      ofstream sol_i_ofs("sol_i.gf");
+      sol_r_ofs.precision(8);
+      sol_i_ofs.precision(8);
+      u.real().Save(sol_r_ofs);
+      u.imag().Save(sol_i_ofs);
+   }
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *mesh << u.real()
+                 << "window_title 'Comp Real Part'" << flush;
+      sol_sock_i << "solution\n" << *mesh << u.imag()
+                 << "window_title 'Comp Imaginary Part'" << flush;
+   }
+   if (visualization && exact_sol)
+   {
+      *u_exact -= u;
+
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *mesh << u_exact->real()
+                 << "window_title 'Exact-Comp Real Part'" << flush;
+      sol_sock_i << "solution\n" << *mesh << u_exact->imag()
+                 << "window_title 'Exact-Comp Imaginary Part'" << flush;
+   }
+   if (visualization)
+   {
+      GridFunction u_t(fespace);
+      u_t = u.real();
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *mesh << u_t
+               << "window_title 'Harmonic Solution (t = 0.0 T)'"
+               << "pause\n" << flush;
+
+      cout << "GLVis visualization paused."
+           << " Press space (in the GLVis window) to resume it.\n";
+      int num_frames = 32;
+      int i = 0;
+      while (sol_sock)
+      {
+         double t = (double)(i % num_frames) / num_frames;
+         ostringstream oss;
+         oss << "Harmonic Solution (t = " << t << " T)";
+
+         add(cos( 2.0 * M_PI * t), u.real(),
+             sin(-2.0 * M_PI * t), u.imag(), u_t);
+         sol_sock << "solution\n" << *mesh << u_t
+                  << "window_title '" << oss.str() << "'" << flush;
+         i++;
+      }
+   }
+
+   // 14. Free the used memory.
+   delete a;
+   delete u_exact;
+   delete fespace;
+   delete fec;
+   delete mesh;
+
+   return 0;
+}
+
+bool check_for_inline_mesh(const char * mesh_file)
+{
+   string file(mesh_file);
+   size_t p0 = file.find_last_of("/");
+   string s0 = file.substr((p0==string::npos)?0:(p0+1),7);
+   return s0 == "inline-";
+}
+
+complex<double> u0_exact(const Vector &x)
+{
+   int dim = x.Size();
+   complex<double> i(0.0, 1.0);
+   complex<double> alpha = (epsilon_ * omega_ - i * sigma_);
+   complex<double> kappa = std::sqrt(mu_ * omega_* alpha);
+   return std::exp(-i * kappa * x[dim - 1]);
+}
+
+double u0_real_exact(const Vector &x)
+{
+   return u0_exact(x).real();
+}
+
+double u0_imag_exact(const Vector &x)
+{
+   return u0_exact(x).imag();
+}
+
+void u1_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
+}
+
+void u1_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
+}
+
+void u2_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
+}
+
+void u2_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
+}

examples/ex21p.cpp

@@ -0,0 +1,658 @@
+//                       MFEM Example 21 - Parallel Version
+//
+// Compile with: make ex21p
+//
+// Sample runs:  mpirun -np 4 ex21p -m ../data/inline-segment.mesh -o 3
+//               mpirun -np 4 ex21p -m ../data/inline-tri.mesh -o 3
+//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3
+//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3 -p 1
+//               mpirun -np 4 ex21p -m ../data/inline-quad.mesh -o 3 -p 2
+//               mpirun -np 4 ex21p -m ../data/inline-tet.mesh -o 2
+//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2
+//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2 -p 1
+//               mpirun -np 4 ex21p -m ../data/inline-hex.mesh -o 2 -p 2
+//               mpirun -np 4 ex21p -m ../data/star.mesh -o 2 -sigma 10.0
+//
+// Description:  This example code demonstrates the use of MFEM to define and
+//               solve simple complex-valued linear systems.  We implement three
+//               variants of a damped harmonic oscillator:
+//
+//               1) A scalar H1 field
+//                  -Div(a Grad u) - omega^2 b u + i omega c u = 0
+//
+//               2) A vector H(Curl) field
+//                  Curl(a Curl u) - omega^2 b u + i omega c u = 0
+//
+//               3) A vector H(Div) field
+//                  -Grad(a Div u) - omega^2 b u + i omega c u = 0
+//
+//               In each case the field is driven by a forced oscillation, with
+//               angular frequency omega, imposed at the boundary or a portion
+//               of the boundary.
+//
+//               In electromagnetics the coefficients are typically named the
+//               permeability, mu = 1/a, permittivity, epsilon = b, and
+//               conductivity, sigma = c.  The user can specify these constants
+//               using either set of names.
+//
+//#define MFEM_STRUMPACK_SRC
+
+#include <fstream>
+#include <iostream>
+#include "mfem.hpp"
+
+using namespace std;
+using namespace mfem;
+
+static double mu_ = 1.0;
+static double epsilon_ = 1.0;
+static double sigma_ = 20.0;
+static double omega_ = 10.0;
+
+double u0_real_exact(const Vector &);
+double u0_imag_exact(const Vector &);
+
+void u1_real_exact(const Vector &, Vector &);
+void u1_imag_exact(const Vector &, Vector &);
+
+void u2_real_exact(const Vector &, Vector &);
+void u2_imag_exact(const Vector &, Vector &);
+
+bool check_for_inline_mesh(const char * mesh_file);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm comm = MPI_COMM_WORLD;
+   MPI_Comm_size(comm, &num_procs);
+   MPI_Comm_rank(comm, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int ser_ref_levels = 1;
+   int par_ref_levels = 1;
+   int order = 1;
+   int prob = 0;
+   double freq = -1.0;
+   double a_coef = 0.0;
+   bool visualization = 1;
+   bool herm_conv = true;
+   bool exact_sol = true;
+#ifdef MFEM_USE_STRUMPACK
+   bool strumpack = false;
+#endif
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&prob, "-p", "--problem-type",
+                  "Choose from 0: H_1, 1: H(Curl), or 2: H(Div) "
+                  "damped harmonic oscillator.");
+   args.AddOption(&a_coef, "-a", "--stiffness-coef",
+                  "Stiffness coefficient (spring constant or 1/mu).");
+   args.AddOption(&epsilon_, "-b", "--mass-coef",
+                  "Mass coefficient (or epsilon).");
+   args.AddOption(&sigma_, "-c", "--damping-coef",
+                  "Damping coefficient (or sigma).");
+   args.AddOption(&mu_, "-mu", "--permeability",
+                  "Permeability of free space (or 1/(spring constant)).");
+   args.AddOption(&epsilon_, "-eps", "--permittivity",
+                  "Permittivity of free space (or mass constant).");
+   args.AddOption(&sigma_, "-sigma", "--conductivity",
+                  "Conductivity (or damping constant).");
+   args.AddOption(&freq, "-f", "--frequency",
+                  "Frequency (in Hz).");
+#ifdef MFEM_USE_STRUMPACK
+   args.AddOption(&strumpack, "-strumpack", "--strumpack-solver",
+                  "-no-strumpack", "--no-strumpack-solver",
+                  "Use STRUMPACK's double complex linear solver.");
+#endif
+   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
+                  "--no-hermitian", "Use convention for Hermitian operators.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   if ( a_coef != 0.0 )
+   {
+      mu_ = 1.0 / a_coef;
+   }
+   if ( freq > 0.0 )
+   {
+      omega_ = 2.0 * M_PI * freq;
+   }
+
+   exact_sol = check_for_inline_mesh(mesh_file);
+   if (myid == 0 && exact_sol)
+   {
+      cout << "Identified an 'inline' mesh" << endl;
+   }
+
+   ComplexOperator::Convention conv =
+      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   int dim = mesh->Dimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution.
+   for (int l = 0; l < ser_ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int l = 0; l < par_ref_levels; l++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define a parallel finite element space on the parallel
+   //    mesh. Here we use continuous Lagrange, Nedelec, or
+   //    Raviart-Thomas finite elements of the specified order.
+   if (dim == 1 && prob != 0 )
+   {
+      if (myid == 0)
+      {
+         cout << "Switching to problem type 0, H1 basis functions, "
+              << "for 1 dimensional mesh." << endl;
+      }
+      prob = 0;
+   }
+
+   FiniteElementCollection *fec;
+   switch (prob)
+   {
+      case 0:  fec = new H1_FECollection(order, dim);      break;
+      case 1:  fec = new ND_FECollection(order, dim);      break;
+      case 2:  fec = new RT_FECollection(order - 1, dim);  break;
+   }
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    based on the type of mesh and the problem type.
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      if (exact_sol)
+      {
+         switch (prob)
+         {
+            case 0:   ess_bdr = 0; ess_bdr[0] = 1;  break;
+            default:  ess_bdr = 1; ess_bdr[2] = 0;  break;
+         }
+      }
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 8. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system.
+   ParComplexLinearForm b(fespace, conv);
+   b.Vector::operator=(0.0);
+
+   // 9. Define the solution vector u as a parallel finite element
+   //    grid function corresponding to fespace. Initialize u with
+   //    initial guess of 1+0i or the exact solution if it is known.
+   ParComplexGridFunction u(fespace);
+   ParComplexGridFunction * u_exact = NULL;
+   if (exact_sol) { u_exact = new ParComplexGridFunction(fespace); }
+
+   FunctionCoefficient u0_r(u0_real_exact);
+   FunctionCoefficient u0_i(u0_imag_exact);
+   VectorFunctionCoefficient u1_r(dim, u1_real_exact);
+   VectorFunctionCoefficient u1_i(dim, u1_imag_exact);
+   VectorFunctionCoefficient u2_r(dim, u2_real_exact);
+   VectorFunctionCoefficient u2_i(dim, u2_imag_exact);
+
+   ConstantCoefficient zeroCoef(0.0);
+   ConstantCoefficient oneCoef(1.0);
+
+   Vector zeroVec(dim); zeroVec = 0.0;
+   Vector  oneVec(dim);  oneVec = 0.0; oneVec[(prob==2)?(dim-1):0] = 1.0;
+   VectorConstantCoefficient zeroVecCoef(zeroVec);
+   VectorConstantCoefficient oneVecCoef(oneVec);
+
+   switch (prob)
+   {
+      case 0:
+         u.ProjectBdrCoefficient(oneCoef, zeroCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u0_r, u0_i); }
+         break;
+      case 1:
+         u.ProjectBdrCoefficientTangent(oneVecCoef, zeroVecCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u1_r, u1_i); }
+         break;
+      case 2:
+         u.ProjectBdrCoefficientNormal(oneVecCoef, zeroVecCoef, ess_bdr);
+         if (exact_sol) { u_exact->ProjectCoefficient(u2_r, u2_i); }
+         break;
+   }
+
+   if (visualization && exact_sol)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
+                 << "window_title 'Exact Real Part'" << flush;
+      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
+                 << "window_title 'Exact Imaginary Part'" << flush;
+   }
+
+   // 10. Set up the parallel sesquilinear form a(.,.) on the finite element
+   //     space corresponding to the damped harmonic oscillator operator
+   //     of the appropriate type:
+   //
+   //     0) A scalar H1 field
+   //        -Div(a Grad) - omega^2 b + i omega c
+   //
+   //     1) A vector H(Curl) field
+   //        Curl(a Curl) - omega^2 b + i omega c
+   //
+   //     2) A vector H(Div) field
+   //        -Grad(a Div) - omega^2 b + i omega c
+   //
+   ConstantCoefficient stiffnessCoef(1.0/mu_);
+   ConstantCoefficient massCoef(-omega_ * omega_ * epsilon_);
+   ConstantCoefficient lossCoef(omega_ * sigma_);
+   ConstantCoefficient negMassCoef(omega_ * omega_ * epsilon_);
+
+   ParSesquilinearForm *a = new ParSesquilinearForm(fespace, conv);
+   switch (prob)
+   {
+      case 0:
+         a->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new MassIntegrator(massCoef),
+                                new MassIntegrator(lossCoef));
+         break;
+      case 1:
+         a->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
+         break;
+      case 2:
+         a->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef),
+                                NULL);
+         a->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef),
+                                new VectorFEMassIntegrator(lossCoef));
+         break;
+   }
+
+   // 10a. Set up the parallel bilinear form for the preconditioner
+   //      corresponding to the appropriate operator if the STRUMPACK solver
+   //      has not been selected.
+   //
+   //      0) A scalar H1 field
+   //         -Div(a Grad) - omega^2 b + omega c
+   //
+   //      1) A vector H(Curl) field
+   //         Curl(a Curl) + omega^2 b + omega c
+   //
+   //      2) A vector H(Div) field
+   //         -Grad(a Div) - omega^2 b + omega c
+   //
+   ParBilinearForm *pcOp = NULL;
+#ifdef MFEM_USE_STRUMPACK
+   if (!strumpack)
+#endif
+   {
+      pcOp = new ParBilinearForm(fespace);
+      switch (prob)
+      {
+         case 0:
+            pcOp->AddDomainIntegrator(new DiffusionIntegrator(stiffnessCoef));
+            pcOp->AddDomainIntegrator(new MassIntegrator(massCoef));
+            pcOp->AddDomainIntegrator(new MassIntegrator(lossCoef));
+            break;
+         case 1:
+            pcOp->AddDomainIntegrator(new CurlCurlIntegrator(stiffnessCoef));
+            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(negMassCoef));
+            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
+            break;
+         case 2:
+            pcOp->AddDomainIntegrator(new DivDivIntegrator(stiffnessCoef));
+            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(massCoef));
+            pcOp->AddDomainIntegrator(new VectorFEMassIntegrator(lossCoef));
+            break;
+      }
+   }
+
+   // 11. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, etc.
+   a->Assemble();
+   if (pcOp) { pcOp->Assemble(); }
+
+   OperatorHandle A;
+   Vector B, U;
+
+   a->FormLinearSystem(ess_tdof_list, u, b, A, U, B);
+   u = 0.0;
+   U = 0.0;
+
+   OperatorHandle PCOp;
+   if (pcOp) { pcOp->FormSystemMatrix(ess_tdof_list, PCOp); }
+
+   if (myid == 0)
+   {
+      ComplexHypreParMatrix * Ahyp =
+         dynamic_cast<ComplexHypreParMatrix*>(A.Ptr());
+
+      cout << "Size of linear system: "
+           << 2 * Ahyp->real().GetGlobalNumRows() << endl << endl;
+   }
+
+   // 12. Define and apply a parallel FGMRES solver for AU=B with a
+   //     block diagonal preconditioner based on the appropriate multigrid
+   //     preconditioner from hypre or simply use STRUMPACK.
+#ifdef MFEM_USE_STRUMPACK
+   if (!strumpack)
+#endif
+   {
+      Array<HYPRE_Int> blockTrueOffsets;
+      blockTrueOffsets.SetSize(3);
+      blockTrueOffsets[0] = 0;
+      blockTrueOffsets[1] = PCOp.Ptr()->Height();
+      blockTrueOffsets[2] = PCOp.Ptr()->Height();
+      blockTrueOffsets.PartialSum();
+
+      BlockDiagonalPreconditioner BDP(blockTrueOffsets);
+
+      Operator * pc_r = NULL;
+      Operator * pc_i = NULL;
+
+      switch (prob)
+      {
+         case 0:
+            pc_r =
+               new HypreBoomerAMG(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
+            pc_i = new ScaledOperator(pc_r,
+                                      (conv == ComplexOperator::HERMITIAN) ?
+                                      1.0:-1.0);
+            break;
+         case 1:
+            pc_r = new HypreAMS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
+                                fespace);
+            pc_i = new ScaledOperator(pc_r,
+                                      (conv == ComplexOperator::HERMITIAN) ?
+                                      1.0:-1.0);
+            break;
+         case 2:
+            if (dim == 2 )
+            {
+               pc_r = new HypreAMS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
+                                   fespace);
+            }
+            else
+            {
+               pc_r = new HypreADS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
+                                   fespace);
+            }
+            pc_i = new ScaledOperator(pc_r,
+                                      (conv == ComplexOperator::HERMITIAN) ?
+                                      1.0:-1.0);
+            break;
+      }
+      BDP.SetDiagonalBlock(0, pc_r);
+      BDP.SetDiagonalBlock(1, pc_i);
+      BDP.owns_blocks = 0;
+
+      FGMRESSolver fgmres(MPI_COMM_WORLD);
+      fgmres.SetPreconditioner(BDP);
+      fgmres.SetOperator(*A.Ptr());
+      fgmres.SetRelTol(1e-12);
+      fgmres.SetMaxIter(1000);
+      fgmres.SetPrintLevel(1);
+      fgmres.Mult(B, U);
+   }
+#ifdef MFEM_USE_STRUMPACK
+   else
+   {
+      ComplexHypreParMatrix * Ahyp =
+         dynamic_cast<ComplexHypreParMatrix*>(A.Ptr());
+
+      STRUMPACKRowLocCmplxMatrix A_strmp(Ahyp->real(), Ahyp->imag());
+
+      STRUMPACKCmplxSolver strmp(argc, argv, comm);
+
+      strmp.SetPrintFactorStatistics(true);
+      strmp.SetPrintSolveStatistics(true);
+      // strmp.SetKrylovSolver(strumpack::KrylovSolver::AUTO); // core dump
+      strmp.SetKrylovSolver(strumpack::KrylovSolver::DIRECT); // core dump
+      // strmp.SetKrylovSolver(strumpack::KrylovSolver::REFINE); // core dump
+      // strmp.SetKrylovSolver(strumpack::KrylovSolver::PREC_GMRES); // index out of range asserts from strumpack::DenseMatrix
+      // strmp.SetKrylovSolver(strumpack::KrylovSolver::GMRES); // WORKS
+      strmp.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+      strmp.SetOperator(A_strmp);
+      strmp.SetFromCommandLine();
+      strmp.Mult(B, U);
+   }
+#endif
+
+   // 13. Recover the parallel grid function corresponding to U. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(U, b, u);
+
+   if (exact_sol)
+   {
+      double err_r = -1.0;
+      double err_i = -1.0;
+
+      switch (prob)
+      {
+         case 0:
+            err_r = u.real().ComputeL2Error(u0_r);
+            err_i = u.imag().ComputeL2Error(u0_i);
+            break;
+         case 1:
+            err_r = u.real().ComputeL2Error(u1_r);
+            err_i = u.imag().ComputeL2Error(u1_i);
+            break;
+         case 2:
+            err_r = u.real().ComputeL2Error(u2_r);
+            err_i = u.imag().ComputeL2Error(u2_i);
+            break;
+      }
+
+      if ( myid == 0 )
+      {
+         cout << endl;
+         cout << "|| Re (u_h - u) ||_{L^2} = " << err_r << endl;
+         cout << "|| Im (u_h - u) ||_{L^2} = " << err_i << endl;
+         cout << endl;
+      }
+   }
+
+   // 14. Save the refined mesh and the solution in parallel. This output can be
+   //     viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_r_name, sol_i_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_r_name << "sol_r." << setfill('0') << setw(6) << myid;
+      sol_i_name << "sol_i." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_r_ofs(sol_r_name.str().c_str());
+      ofstream sol_i_ofs(sol_i_name.str().c_str());
+      sol_r_ofs.precision(8);
+      sol_i_ofs.precision(8);
+      u.real().Save(sol_r_ofs);
+      u.imag().Save(sol_i_ofs);
+   }
+
+   // 15. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *pmesh << u.real()
+                 << "window_title 'Comp Real Part'" << flush;
+      sol_sock_i << "solution\n" << *pmesh << u.imag()
+                 << "window_title 'Comp Imaginary Part'" << flush;
+   }
+   if (visualization && exact_sol)
+   {
+      *u_exact -= u;
+
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
+      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
+                 << "window_title 'Exact-Comp Real Part'" << flush;
+      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
+                 << "window_title 'Exact-Comp Imaginary Part'" << flush;
+   }
+   if (visualization)
+   {
+      ParGridFunction u_t(fespace);
+      u_t = u.real();
+      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 << u_t
+               << "window_title 'Harmonic Solution (t = 0.0 T)'"
+               << "pause\n" << flush;
+      if (myid == 0)
+         cout << "GLVis visualization paused."
+              << " Press space (in the GLVis window) to resume it.\n";
+      int num_frames = 32;
+      int i = 0;
+      while (sol_sock)
+      {
+         double t = (double)(i % num_frames) / num_frames;
+         ostringstream oss;
+         oss << "Harmonic Solution (t = " << t << " T)";
+
+         add(cos( 2.0 * M_PI * t), u.real(),
+             sin(-2.0 * M_PI * t), u.imag(), u_t);
+         sol_sock << "parallel " << num_procs << " " << myid << "\n";
+         sol_sock << "solution\n" << *pmesh << u_t
+                  << "window_title '" << oss.str() << "'" << flush;
+         i++;
+      }
+   }
+
+   // 16. Free the used memory.
+   delete a;
+   delete u_exact;
+   delete pcOp;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+bool check_for_inline_mesh(const char * mesh_file)
+{
+   string file(mesh_file);
+   size_t p0 = file.find_last_of("/");
+   string s0 = file.substr((p0==string::npos)?0:(p0+1),7);
+   return s0 == "inline-";
+}
+
+complex<double> u0_exact(const Vector &x)
+{
+   int dim = x.Size();
+   complex<double> i(0.0, 1.0);
+   complex<double> alpha = (epsilon_ * omega_ - i * sigma_);
+   complex<double> kappa = std::sqrt(mu_ * omega_* alpha);
+   return std::exp(-i * kappa * x[dim - 1]);
+}
+
+double u0_real_exact(const Vector &x)
+{
+   return u0_exact(x).real();
+}
+
+double u0_imag_exact(const Vector &x)
+{
+   return u0_exact(x).imag();
+}
+
+void u1_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_real_exact(x);
+}
+
+void u1_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[0] = u0_imag_exact(x);
+}
+
+void u2_real_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_real_exact(x);
+}
+
+void u2_imag_exact(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+   v.SetSize(dim); v = 0.0; v[dim-1] = u0_imag_exact(x);
+}

examples/ex3p_sp.cpp

@@ -0,0 +1,334 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p
+//
+// Sample runs:  mpirun -np 4 ex3p -m ../data/star.mesh
+//               mpirun -np 4 ex3p -m ../data/square-disc.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/beam-tet.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/escher.mesh
+//               mpirun -np 4 ex3p -m ../data/escher.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/fichera.mesh
+//               mpirun -np 4 ex3p -m ../data/fichera-q2.vtk
+//               mpirun -np 4 ex3p -m ../data/fichera-q3.mesh
+//               mpirun -np 4 ex3p -m ../data/square-disc-nurbs.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex-nurbs.mesh
+//               mpirun -np 4 ex3p -m ../data/amr-quad.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/star-surf.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/mobius-strip.mesh -o 2 -f 0.1
+//               mpirun -np 4 ex3p -m ../data/klein-bottle.mesh -o 2 -f 0.1
+//
+// Description:  This example code solves a simple electromagnetic diffusion
+//               problem corresponding to the second order definite Maxwell
+//               equation curl curl E + E = f with boundary condition
+//               E x n = <given tangential field>. Here, we use a given exact
+//               solution E and compute the corresponding r.h.s. f.
+//               We discretize with Nedelec finite elements in 2D or 3D.
+//
+//               The example demonstrates the use of H(curl) finite element
+//               spaces with the curl-curl and the (vector finite element) mass
+//               bilinear form, as well as the computation of discretization
+//               error when the exact solution is known. Static condensation is
+//               also illustrated.
+//
+//               We recommend viewing examples 1-2 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution, E, and r.h.s., f. See below for implementation.
+void E_exact(const Vector &, Vector &);
+void f_exact(const Vector &, Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool visualization = 1;
+#ifdef MFEM_USE_STRUMPACK
+   bool use_strumpack = false;
+#endif
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+#ifdef MFEM_USE_STRUMPACK
+   args.AddOption(&use_strumpack, "-strumpack", "--strumpack-solver",
+                  "-no-strumpack", "--no-strumpack-solver",
+                  "Use STRUMPACK's double complex linear solver.");
+#endif
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   kappa = freq * M_PI;
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 1,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(100000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
+   //    meshes need to be reoriented before we can define high-order Nedelec
+   //    spaces on them.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 2;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new ND_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 7. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 8. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
+   //    basis functions in the finite element fespace.
+   VectorFunctionCoefficient f(sdim, f_exact);
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+   b->Assemble();
+
+   // 9. Define the solution vector x as a parallel finite element grid function
+   //    corresponding to fespace. Initialize x by projecting the exact
+   //    solution. Note that only values from the boundary edges will be used
+   //    when eliminating the non-homogeneous boundary condition to modify the
+   //    r.h.s. vector b.
+   ParGridFunction x(fespace);
+   VectorFunctionCoefficient E(sdim, E_exact);
+   x.ProjectCoefficient(E);
+
+   // 10. Set up the parallel bilinear form corresponding to the EM diffusion
+   //     operator curl muinv curl + sigma I, by adding the curl-curl and the
+   //     mass domain integrators.
+   Coefficient *muinv = new ConstantCoefficient(1.0);
+   Coefficient *sigma = new ConstantCoefficient(-1.0);
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
+
+   // 11. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   HypreParMatrix A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   if (myid == 0)
+   {
+      cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
+   }
+
+   StopWatch chrono;
+   chrono.Clear();
+   chrono.Start();
+
+#ifdef MFEM_USE_STRUMPACK
+   if (use_strumpack)
+   {
+      Operator * Arow = new STRUMPACKRowLocMatrix(A);
+
+      STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
+      strumpack->SetPrintFactorStatistics(true);
+      strumpack->SetPrintSolveStatistics(false);
+      strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+      strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+      // strumpack->SetMC64Job(strumpack::MC64Job::NONE);
+      // strumpack->SetSymmetricPattern(true);
+      strumpack->SetOperator(*Arow);
+      strumpack->SetFromCommandLine();
+      //Solver * precond = strumpack;
+
+      strumpack->Mult(B, X);
+
+      delete strumpack;
+      delete Arow;
+   }
+   else
+#endif
+   {
+      // 12. Define and apply a parallel PCG solver for AX=B with the AMS
+      //     preconditioner from hypre.
+      ParFiniteElementSpace *prec_fespace =
+         (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
+      HypreSolver *ams = new HypreAMS(A, prec_fespace);
+      HyprePCG *pcg = new HyprePCG(A);
+      pcg->SetTol(1e-12);
+      pcg->SetMaxIter(500);
+      pcg->SetPrintLevel(2);
+      pcg->SetPreconditioner(*ams);
+      pcg->Mult(B, X);
+
+      delete pcg;
+      delete ams;
+   }
+
+   chrono.Stop();
+   cout << "Solver time " << chrono.RealTime() << endl;
+
+   // 13. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 14. Compute and print the L^2 norm of the error.
+   {
+      double err = x.ComputeL2Error(E);
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 15. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      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 << x << flush;
+   }
+
+   // 17. Free the used memory.
+   delete a;
+   delete sigma;
+   delete muinv;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+
+void E_exact(const Vector &x, Vector &E)
+{
+   if (dim == 3)
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(2));
+      E(2) = sin(kappa * x(0));
+   }
+   else
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(0));
+      if (x.Size() == 3) { E(2) = 0.0; }
+   }
+}
+
+void f_exact(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
+      f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
+   }
+   else
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}

examples/ex6.cpp

@@ -16,9 +16,9 @@
 //               ex6 -m ../data/amr-quad.mesh
 //
 // Device sample runs:
-//               ex6 -pa -d cuda
-//               ex6 -pa -d occa-cuda
-//               ex6 -pa -d raja-omp
+//             > ex6 -pa -d cuda
+//             > ex6 -pa -d occa-cuda
+//             > ex6 -pa -d raja-omp
 //
 // Description:  This is a version of Example 1 with a simple adaptive mesh
 //               refinement loop. The problem being solved is again the Laplace

examples/ex6p.cpp

@@ -16,9 +16,9 @@
 //               mpirun -np 4 ex6p -m ../data/amr-quad.mesh
 //
 // Device sample runs:
-//               mpirun -np 4 ex6p -pa -d cuda
-//               mpirun -np 4 ex6p -pa -d occa-cuda
-//               mpirun -np 4 ex6p -pa -d raja-omp
+//             > mpirun -np 4 ex6p -pa -d cuda
+//             > mpirun -np 4 ex6p -pa -d occa-cuda
+//             > mpirun -np 4 ex6p -pa -d raja-omp
 //
 // Description:  This is a version of Example 1 with a simple adaptive mesh
 //               refinement loop. The problem being solved is again the Laplace

examples/makefile

@@ -22,9 +22,9 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES = ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 ex17\
-  ex18 ex19 ex20 ex22
+  ex18 ex19 ex20 ex21 ex22
 PAR_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p ex12p\
- ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex22p
+ ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p
 
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
@@ -118,7 +118,7 @@ clean-build:
 clean-exec:
 	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
 	@rm -rf Example5* Example9* Example15* Example16*
-	@rm -f sphere_refined.* sol.* sol_u.* sol_p.*
+	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
 	@rm -f ex9.mesh ex9-mesh.* ex9-init.* ex9-final.*
 	@rm -f deformed.* velocity.* elastic_energy.* mode_*
 	@rm -f ex16.mesh ex16-mesh.* ex16-init.* ex16-final.*

fem/bilinearform.cpp

@@ -588,18 +588,14 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
                                     Vector &b, OperatorHandle &A, Vector &X,
                                     Vector &B, int copy_interior)
 {
-   const SparseMatrix *P = fes->GetConformingProlongation();
-
    if (ext)
    {
-      if (P != NULL && assembly != AssemblyLevel::FULL && Device::IsEnabled())
-      {
-         P->BuildTranspose();
-      }
       ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
       return;
    }
 
+   const SparseMatrix *P = fes->GetConformingProlongation();
+
    FormSystemMatrix(ess_tdof_list, A);
 
    // Transform the system and perform the elimination in B, based on the

fem/complex_fem.cpp

@@ -0,0 +1,784 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#include "complex_fem.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+ComplexGridFunction::ComplexGridFunction(FiniteElementSpace *fes)
+   : Vector(2*(fes->GetVSize()))
+{
+   gfr_ = new GridFunction(fes, &data[0]);
+   gfi_ = new GridFunction(fes, &data[fes->GetVSize()]);
+}
+
+void
+ComplexGridFunction::Update()
+{
+   FiniteElementSpace * fes = gfr_->FESpace();
+
+   int vsize = fes->GetVSize();
+
+   const Operator *T = fes->GetUpdateOperator();
+   if (T)
+   {
+      // Update the individual GridFunction objects.  This will allocate
+      // new data arrays for each GridFunction.
+      gfr_->Update();
+      gfi_->Update();
+
+      // Our data array now contains old data as well as being the wrong size
+      // so reallocate it.
+      this->SetSize(2 * vsize);
+
+      // Create temporary vectors which point to the new data array
+      Vector gf_r(&data[0], vsize);
+      Vector gf_i(&data[vsize], vsize);
+
+      // Copy the updated GridFunctions into the new data array
+      gf_r = *gfr_;
+      gf_i = *gfi_;
+
+      // Replace the individual data arrays with pointers into the new data array
+      gfr_->NewDataAndSize(&data[0], vsize);
+      gfi_->NewDataAndSize(&data[vsize], vsize);
+   }
+   else
+   {
+      // The existing data will not be transferred to the new GridFunctions
+      // so delete it a allocate a new array
+      this->SetSize(2 * vsize);
+
+      // Point the individual GridFunctions to the new data array
+      gfr_->NewDataAndSize(&data[0], vsize);
+      gfi_->NewDataAndSize(&data[vsize], vsize);
+
+      // These updates will only set the proper 'sequence' value within
+      // the individual GridFunction objects because their sizes are
+      // already correct
+      gfr_->Update();
+      gfi_->Update();
+   }
+}
+
+void
+ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
+                                        Coefficient &imag_coeff)
+{
+   gfr_->ProjectCoefficient(real_coeff);
+   gfi_->ProjectCoefficient(imag_coeff);
+}
+
+void
+ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                        VectorCoefficient &imag_vcoeff)
+{
+   gfr_->ProjectCoefficient(real_vcoeff);
+   gfi_->ProjectCoefficient(imag_vcoeff);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                           Coefficient &imag_coeff,
+                                           Array<int> &attr)
+{
+   gfr_->ProjectBdrCoefficient(real_coeff, attr);
+   gfi_->ProjectBdrCoefficient(imag_coeff, attr);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
+                                                 VectorCoefficient &imag_vcoeff,
+                                                 Array<int> &attr)
+{
+   gfr_->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   gfi_->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
+}
+
+void
+ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                  &real_vcoeff,
+                                                  VectorCoefficient
+                                                  &imag_vcoeff,
+                                                  Array<int> &attr)
+{
+   gfr_->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   gfi_->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
+}
+
+
+ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *f,
+                                     ComplexOperator::Convention convention)
+   : Vector(2*(f->GetVSize())),
+     conv_(convention)
+{
+   lfr_ = new LinearForm(f, &data[0]);
+   lfi_ = new LinearForm(f, &data[f->GetVSize()]);
+}
+
+ComplexLinearForm::~ComplexLinearForm()
+{
+   delete lfr_;
+   delete lfi_;
+}
+
+void
+ComplexLinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi_real,
+                                       LinearFormIntegrator *lfi_imag)
+{
+   if ( lfi_real ) { lfr_->AddDomainIntegrator(lfi_real); }
+   if ( lfi_imag ) { lfi_->AddDomainIntegrator(lfi_imag); }
+}
+
+void
+ComplexLinearForm::Update()
+{
+   FiniteElementSpace *fes = lfr_->FESpace();
+
+   this->Update(fes);
+}
+
+void
+ComplexLinearForm::Update(FiniteElementSpace *fes)
+{
+   int vsize = fes->GetVSize();
+   SetSize(2 * vsize);
+
+   Vector lfr(&data[0], vsize);
+   Vector lfi(&data[vsize], vsize);
+
+   lfr_->Update(fes, lfr, 0);
+   lfi_->Update(fes, lfi, 0);
+}
+
+void
+ComplexLinearForm::Assemble()
+{
+   lfr_->Assemble();
+   lfi_->Assemble();
+   if (conv_ == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      *lfi_ *= -1.0;
+   }
+}
+
+complex<double>
+ComplexLinearForm::operator()(const ComplexGridFunction &gf) const
+{
+   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
+   return complex<double>((*lfr_)(gf.real()) - s * (*lfi_)(gf.imag()),
+                          (*lfr_)(gf.imag()) + s * (*lfi_)(gf.real()));
+}
+
+
+SesquilinearForm::SesquilinearForm(FiniteElementSpace *f,
+                                   ComplexOperator::Convention convention)
+   : conv_(convention),
+     blfr_(new BilinearForm(f)),
+     blfi_(new BilinearForm(f))
+{}
+
+SesquilinearForm::~SesquilinearForm()
+{
+   delete blfr_;
+   delete blfi_;
+}
+
+void SesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                                           BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { blfr_->AddDomainIntegrator(bfi_real); }
+   if (bfi_imag) { blfi_->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+SesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                        BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { blfr_->AddBoundaryIntegrator(bfi_real); }
+   if (bfi_imag) { blfi_->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+SesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                        BilinearFormIntegrator *bfi_imag,
+                                        Array<int> & bdr_marker)
+{
+   if (bfi_real) { blfr_->AddBoundaryIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { blfi_->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+SesquilinearForm::Assemble(int skip_zeros)
+{
+   blfr_->Assemble(skip_zeros);
+   blfi_->Assemble(skip_zeros);
+}
+
+void
+SesquilinearForm::Finalize(int skip_zeros)
+{
+   blfr_->Finalize(skip_zeros);
+   blfi_->Finalize(skip_zeros);
+}
+
+ComplexSparseMatrix *
+SesquilinearForm::AssembleCompSpMat()
+{
+   return new ComplexSparseMatrix(&blfr_->SpMat(),
+                                  &blfi_->SpMat(),
+                                  false, false, conv_);
+
+}
+
+void
+SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
+                                   Vector &x, Vector &b,
+                                   OperatorHandle &A,
+                                   Vector &X, Vector &B,
+                                   int ci)
+{
+   FiniteElementSpace * fes = blfr_->FESpace();
+
+   int vsize  = fes->GetVSize();
+   // int tvsize = pfes->GetTrueVSize();
+
+   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
+
+   // Allocate temporary vectors
+   Vector b_0(vsize);  b_0 = 0.0;
+   // Vector B_0(tvsize); B_0 = 0.0;
+
+   // Extract the real and imaginary parts of the input vectors
+   MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
+   Vector x_r(x.GetData(), vsize);
+   Vector x_i(&(x.GetData())[vsize], vsize);
+
+   MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
+   Vector b_r(b.GetData(), vsize);
+   Vector b_i(&(b.GetData())[vsize], vsize);
+   b_i *= s;
+   /*
+   X.SetSize(2 * tvsize);
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+   B.SetSize(2 * tvsize);
+   Vector B_r(B.GetData(), tvsize);
+   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   */
+   SparseMatrix * A_r = new SparseMatrix;
+   SparseMatrix * A_i = new SparseMatrix;
+   Vector X_0, B_0;
+
+   b_0 = b_r;
+   blfr_->FormLinearSystem(ess_tdof_list, x_r, b_r, *A_r, X_0, B_0, ci);
+
+   int tvsize = B_0.Size();
+   X.SetSize(2 * tvsize);
+   B.SetSize(2 * tvsize);
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+   Vector B_r(B.GetData(), tvsize);
+   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   X_r = X_0; B_r = B_0;
+
+   b_0 = 0.0;
+   blfi_->FormLinearSystem(ess_tdof_list, x_i, b_0, *A_i, X_0, B_0, false);
+   B_r -= B_0;
+
+   b_0 = b_i;
+   blfr_->FormLinearSystem(ess_tdof_list, x_i, b_0, *A_r, X_0, B_0, ci);
+   X_i = X_0; B_i = B_0;
+
+   b_0 = 0.0;
+   blfi_->FormLinearSystem(ess_tdof_list, x_r, b_0, *A_i, X_0, B_0, false);
+   B_i += B_0;
+
+   B_i *= s;
+   b_i *= s;
+
+   // A = A_r + i A_i
+   A.Clear();
+   ComplexSparseMatrix * A_sp =
+      new ComplexSparseMatrix(A_r, A_i, true, true, conv_);
+   A.Reset<ComplexSparseMatrix>(A_sp, true);
+}
+
+void
+SesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
+                                     Vector &x)
+{
+   FiniteElementSpace * fes = blfr_->FESpace();
+
+   const SparseMatrix *P = fes->GetConformingProlongation();
+
+   int vsize  = fes->GetVSize();
+   int tvsize = X.Size() / 2;
+
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+   Vector x_r(x.GetData(), vsize);
+   Vector x_i(&(x.GetData())[vsize], vsize);
+
+   if (!P)
+   {
+      x = X;
+   }
+   else
+   {
+      // Apply conforming prolongation
+      P->Mult(X_r, x_r);
+      P->Mult(X_i, x_i);
+   }
+}
+
+void
+SesquilinearForm::Update(FiniteElementSpace *nfes)
+{
+   if ( blfr_ ) { blfr_->Update(nfes); }
+   if ( blfi_ ) { blfi_->Update(nfes); }
+}
+
+
+#ifdef MFEM_USE_MPI
+
+ParComplexGridFunction::ParComplexGridFunction(ParFiniteElementSpace *pfes)
+   : Vector(2*(pfes->GetVSize()))
+{
+   pgfr_ = new ParGridFunction(pfes, &data[0]);
+   pgfi_ = new ParGridFunction(pfes, &data[pfes->GetVSize()]);
+}
+
+void
+ParComplexGridFunction::Update()
+{
+   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+
+   int vsize = pfes->GetVSize();
+
+   const Operator *T = pfes->GetUpdateOperator();
+   if (T)
+   {
+      // Update the individual GridFunction objects.  This will allocate
+      // new data arrays for each GridFunction.
+      pgfr_->Update();
+      pgfi_->Update();
+
+      // Our data array now contains old data as well as being the wrong size
+      // so reallocate it.
+      this->SetSize(2 * vsize);
+
+      // Create temporary vectors which point to the new data array
+      Vector gf_r(&data[0], vsize);
+      Vector gf_i(&data[vsize], vsize);
+
+      // Copy the updated GridFunctions into the new data array
+      gf_r = *pgfr_;
+      gf_i = *pgfi_;
+
+      // Replace the individual data arrays with pointers into the new data array
+      pgfr_->NewDataAndSize(&data[0], vsize);
+      pgfi_->NewDataAndSize(&data[vsize], vsize);
+   }
+   else
+   {
+      // The existing data will not be transferred to the new GridFunctions
+      // so delete it a allocate a new array
+      this->SetSize(2 * vsize);
+
+      // Point the individual GridFunctions to the new data array
+      pgfr_->NewDataAndSize(&data[0], vsize);
+      pgfi_->NewDataAndSize(&data[vsize], vsize);
+
+      // These updates will only set the proper 'sequence' value within
+      // the individual GridFunction objects because their sizes are
+      // already correct
+      pgfr_->Update();
+      pgfi_->Update();
+   }
+}
+
+void
+ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
+                                           Coefficient &imag_coeff)
+{
+   pgfr_->ProjectCoefficient(real_coeff);
+   pgfi_->ProjectCoefficient(imag_coeff);
+}
+
+void
+ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                           VectorCoefficient &imag_vcoeff)
+{
+   pgfr_->ProjectCoefficient(real_vcoeff);
+   pgfi_->ProjectCoefficient(imag_vcoeff);
+}
+
+void
+ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
+                                              Coefficient &imag_coeff,
+                                              Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficient(real_coeff, attr);
+   pgfi_->ProjectBdrCoefficient(imag_coeff, attr);
+}
+
+void
+ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
+                                                    &real_vcoeff,
+                                                    VectorCoefficient
+                                                    &imag_vcoeff,
+                                                    Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficientNormal(real_vcoeff, attr);
+   pgfi_->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
+}
+
+void
+ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
+                                                     &real_vcoeff,
+                                                     VectorCoefficient
+                                                     &imag_vcoeff,
+                                                     Array<int> &attr)
+{
+   pgfr_->ProjectBdrCoefficientTangent(real_vcoeff, attr);
+   pgfi_->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
+}
+
+void
+ParComplexGridFunction::Distribute(const Vector *tv)
+{
+   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+   HYPRE_Int size = pfes->GetTrueVSize();
+
+   double * tvd = tv->GetData();
+   Vector tvr(tvd, size);
+   Vector tvi(&tvd[size], size);
+
+   pgfr_->Distribute(tvr);
+   pgfi_->Distribute(tvi);
+}
+
+void
+ParComplexGridFunction::ParallelProject(Vector &tv) const
+{
+   ParFiniteElementSpace * pfes = pgfr_->ParFESpace();
+   HYPRE_Int size = pfes->GetTrueVSize();
+
+   double * tvd = tv.GetData();
+   Vector tvr(tvd, size);
+   Vector tvi(&tvd[size], size);
+
+   pgfr_->ParallelProject(tvr);
+   pgfi_->ParallelProject(tvi);
+}
+
+
+ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
+                                           ComplexOperator::Convention
+                                           convention)
+   : Vector(2*(pfes->GetVSize())),
+     conv_(convention)
+{
+   plfr_ = new ParLinearForm(pfes, &data[0]);
+   plfi_ = new ParLinearForm(pfes, &data[pfes->GetVSize()]);
+
+   HYPRE_Int * tdof_offsets = pfes->GetTrueDofOffsets();
+
+   int n = (HYPRE_AssumedPartitionCheck()) ? 2 : pfes->GetNRanks();
+   tdof_offsets_ = new HYPRE_Int[n+1];
+
+   for (int i=0; i<=n; i++)
+   {
+      tdof_offsets_[i] = 2 * tdof_offsets[i];
+   }
+}
+
+ParComplexLinearForm::~ParComplexLinearForm()
+{
+   delete plfr_;
+   delete plfi_;
+   delete [] tdof_offsets_;
+}
+
+void
+ParComplexLinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi_real,
+                                          LinearFormIntegrator *lfi_imag)
+{
+   if ( lfi_real ) { plfr_->AddDomainIntegrator(lfi_real); }
+   if ( lfi_imag ) { plfi_->AddDomainIntegrator(lfi_imag); }
+}
+
+void
+ParComplexLinearForm::Update(ParFiniteElementSpace *pf)
+{
+   ParFiniteElementSpace *pfes = (pf!=NULL)?pf:plfr_->ParFESpace();
+   int vsize = pfes->GetVSize();
+   SetSize(2 * vsize);
+
+   Vector plfr(&data[0], vsize);
+   Vector plfi(&data[vsize], vsize);
+
+   plfr_->Update(pfes, plfr, 0);
+   plfi_->Update(pfes, plfi, 0);
+}
+
+void
+ParComplexLinearForm::Assemble()
+{
+   plfr_->Assemble();
+   plfi_->Assemble();
+   if (conv_ == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      *plfi_ *= -1.0;
+   }
+}
+
+void
+ParComplexLinearForm::ParallelAssemble(Vector &tv)
+{
+   HYPRE_Int size = plfr_->ParFESpace()->GetTrueVSize();
+
+   double * tvd = tv.GetData();
+   Vector tvr(tvd, size);
+   Vector tvi(&tvd[size], size);
+
+   plfr_->ParallelAssemble(tvr);
+   plfi_->ParallelAssemble(tvi);
+}
+
+HypreParVector *
+ParComplexLinearForm::ParallelAssemble()
+{
+   const ParFiniteElementSpace * pfes = plfr_->ParFESpace();
+
+   HypreParVector * tv = new HypreParVector(pfes->GetComm(),
+                                            2*(pfes->GlobalTrueVSize()),
+                                            tdof_offsets_);
+
+   HYPRE_Int size = pfes->GetTrueVSize();
+
+   double * tvd = tv->GetData();
+   Vector tvr(tvd, size);
+   Vector tvi(&tvd[size], size);
+
+   plfr_->ParallelAssemble(tvr);
+   plfi_->ParallelAssemble(tvi);
+
+   return tv;
+}
+
+complex<double>
+ParComplexLinearForm::operator()(const ParComplexGridFunction &gf) const
+{
+   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
+   return complex<double>((*plfr_)(gf.real()) - s * (*plfi_)(gf.imag()),
+                          (*plfr_)(gf.imag()) + s * (*plfi_)(gf.real()));
+}
+
+
+ParSesquilinearForm::ParSesquilinearForm(ParFiniteElementSpace *pf,
+                                         ComplexOperator::Convention
+                                         convention)
+   : conv_(convention),
+     pblfr_(new ParBilinearForm(pf)),
+     pblfi_(new ParBilinearForm(pf))
+{}
+
+ParSesquilinearForm::~ParSesquilinearForm()
+{
+   delete pblfr_;
+   delete pblfi_;
+}
+
+void ParSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                                              BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr_->AddDomainIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi_->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+ParSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                           BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr_->AddBoundaryIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi_->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+ParSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                                           BilinearFormIntegrator *bfi_imag,
+                                           Array<int> & bdr_marker)
+{
+   if (bfi_real) { pblfr_->AddBoundaryIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { pblfi_->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+ParSesquilinearForm::Assemble(int skip_zeros)
+{
+   pblfr_->Assemble(skip_zeros);
+   pblfi_->Assemble(skip_zeros);
+}
+
+void
+ParSesquilinearForm::Finalize(int skip_zeros)
+{
+   pblfr_->Finalize(skip_zeros);
+   pblfi_->Finalize(skip_zeros);
+}
+
+ComplexHypreParMatrix *
+ParSesquilinearForm::ParallelAssemble()
+{
+   return new ComplexHypreParMatrix(pblfr_->ParallelAssemble(),
+                                    pblfi_->ParallelAssemble(),
+                                    true, true, conv_);
+
+}
+
+void
+ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
+                                      Vector &x, Vector &b,
+                                      OperatorHandle &A,
+                                      Vector &X, Vector &B,
+                                      int ci)
+{
+   ParFiniteElementSpace * pfes = pblfr_->ParFESpace();
+
+   int tvs = pfes->TrueVSize();
+   cout << "TrueVSize returns " << tvs << endl;
+   cout << "GetVSize returns " << pfes->GetVSize() << endl;
+
+   int vsize = x.Size() / 2;
+   // int vsize  = pfes->GetVSize();
+   // int tvsize = pfes->GetTrueVSize();
+
+   cout << "x.Size/2 returns " << vsize << endl;
+
+   double s = (conv_ == ComplexOperator::HERMITIAN)?1.0:-1.0;
+
+   // Allocate temporary vectors
+   Vector b_0(vsize);  b_0 = 0.0;
+   // Vector B_0(tvsize); B_0 = 0.0;
+
+   // Extract the real and imaginary parts of the input vectors
+   // MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
+   Vector x_r(x.GetData(), vsize);
+   Vector x_i(&(x.GetData())[vsize], vsize);
+
+   MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
+   Vector b_r(b.GetData(), vsize);
+   Vector b_i(&(b.GetData())[vsize], vsize);
+   b_i *= s;
+   /*
+   X.SetSize(2 * tvsize);
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+   B.SetSize(2 * tvsize);
+   Vector B_r(B.GetData(), tvsize);
+   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   */
+   OperatorHandle A_r, A_i;
+   Vector X_0, B_0;
+   cout << "pblfr fls 1" << endl << flush;
+   b_0 = b_r;
+   pblfr_->FormLinearSystem(ess_tdof_list, x_r, b_0, A_r, X_0, B_0, ci);
+
+   int tvsize = B_0.Size();
+   X.SetSize(2 * tvsize);
+   B.SetSize(2 * tvsize);
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+   Vector B_r(B.GetData(), tvsize);
+   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   X_r = X_0; B_r = B_0;
+   cout << "pblfi fls 1" << endl << flush;
+   b_0 = 0.0;
+   pblfi_->FormLinearSystem(ess_tdof_list, x_i, b_0, A_i, X_0, B_0, false);
+   B_r -= B_0;
+   cout << "pblfr fls 2" << endl << flush;
+   b_0 = b_i;
+   pblfr_->FormLinearSystem(ess_tdof_list, x_i, b_0, A_r, X_0, B_0, ci);
+   X_i = X_0; B_i = B_0;
+   cout << "pblfi fls 2" << endl << flush;
+   b_0 = 0.0;
+   pblfi_->FormLinearSystem(ess_tdof_list, x_r, b_0, A_i, X_0, B_0, false);
+   B_i += B_0;
+
+   B_i *= s;
+   b_i *= s;
+
+   // A = A_r + i A_i
+   A.Clear();
+   if ( A_r.Type() == Operator::Hypre_ParCSR &&
+        A_i.Type() == Operator::Hypre_ParCSR )
+   {
+      ComplexHypreParMatrix * A_hyp =
+         new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
+                                   A_i.As<HypreParMatrix>(),
+                                   A_r.OwnsOperator(),
+                                   A_i.OwnsOperator(),
+                                   conv_);
+      A.Reset<ComplexHypreParMatrix>(A_hyp, true);
+   }
+   else
+   {
+      ComplexOperator * A_op =
+         new ComplexOperator(A_r.As<Operator>(),
+                             A_i.As<Operator>(),
+                             A_r.OwnsOperator(),
+                             A_i.OwnsOperator(),
+                             conv_);
+      A.Reset<ComplexOperator>(A_op, true);
+   }
+}
+
+void
+ParSesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
+                                        Vector &x)
+{
+   ParFiniteElementSpace * pfes = pblfr_->ParFESpace();
+
+   const Operator &P = *pfes->GetProlongationMatrix();
+
+   int vsize  = pfes->GetVSize();
+   int tvsize = X.Size() / 2;
+
+   Vector X_r(X.GetData(), tvsize);
+   Vector X_i(&(X.GetData())[tvsize], tvsize);
+
+   Vector x_r(x.GetData(), vsize);
+   Vector x_i(&(x.GetData())[vsize], vsize);
+
+   // Apply conforming prolongation
+   P.Mult(X_r, x_r);
+   P.Mult(X_i, x_i);
+}
+
+void
+ParSesquilinearForm::Update(FiniteElementSpace *nfes)
+{
+   if ( pblfr_ ) { pblfr_->Update(nfes); }
+   if ( pblfi_ ) { pblfi_->Update(nfes); }
+}
+
+
+#endif // MFEM_USE_MPI
+
+}

fem/complex_fem.hpp

@@ -0,0 +1,356 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_COMPLEX_FEM
+#define MFEM_COMPLEX_FEM
+
+#include "../linalg/complex_operator.hpp"
+#include "gridfunc.hpp"
+#include "linearform.hpp"
+#include "bilinearform.hpp"
+#ifdef MFEM_USE_MPI
+#include "pgridfunc.hpp"
+#include "plinearform.hpp"
+#include "pbilinearform.hpp"
+#endif
+#include <complex>
+
+namespace mfem
+{
+
+/// Class for complex-valued grid function - Vector with associated FE space.
+class ComplexGridFunction : public Vector
+{
+private:
+
+   GridFunction * gfr_;
+   GridFunction * gfi_;
+
+protected:
+   void Destroy() { delete gfr_; delete gfi_; }
+
+public:
+
+   /* @brief Construct a ComplexGridFunction associated with the
+      FiniteElementSpace @a *f. */
+   ComplexGridFunction(FiniteElementSpace *f);
+
+   void Update();
+
+   /// Assign constant values to the ComplexGridFunction data.
+   ComplexGridFunction &operator=(const std::complex<double> & value)
+   { *gfr_ = value.real(); *gfi_ = value.imag(); return *this; }
+
+   virtual void ProjectCoefficient(Coefficient &real_coeff,
+                                   Coefficient &imag_coeff);
+   virtual void ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                   VectorCoefficient &imag_vcoeff);
+
+   virtual void ProjectBdrCoefficient(Coefficient &real_coeff,
+                                      Coefficient &imag_coeff,
+                                      Array<int> &attr);
+   virtual void ProjectBdrCoefficientNormal(VectorCoefficient &real_coeff,
+                                            VectorCoefficient &imag_coeff,
+                                            Array<int> &attr);
+   virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff,
+                                             VectorCoefficient &imag_coeff,
+                                             Array<int> &attr);
+
+   FiniteElementSpace *FESpace() { return gfr_->FESpace(); }
+   const FiniteElementSpace *FESpace() const { return gfr_->FESpace(); }
+
+   GridFunction & real() { return *gfr_; }
+   GridFunction & imag() { return *gfi_; }
+   const GridFunction & real() const { return *gfr_; }
+   const GridFunction & imag() const { return *gfi_; }
+
+   /// Destroys grid function.
+   virtual ~ComplexGridFunction() { Destroy(); }
+
+};
+
+class ComplexLinearForm : public Vector
+{
+private:
+   ComplexOperator::Convention conv_;
+
+protected:
+   LinearForm * lfr_;
+   LinearForm * lfi_;
+
+   // HYPRE_Int * tdof_offsets_;
+
+public:
+
+   ComplexLinearForm(FiniteElementSpace *fes,
+                     ComplexOperator::Convention
+                     convention = ComplexOperator::HERMITIAN);
+
+   virtual ~ComplexLinearForm();
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(LinearFormIntegrator *lfi_real,
+                            LinearFormIntegrator *lfi_imag);
+
+   FiniteElementSpace *FESpace() const { return lfr_->FESpace(); }
+
+   LinearForm & real() { return *lfr_; }
+   LinearForm & imag() { return *lfi_; }
+   const LinearForm & real() const { return *lfr_; }
+   const LinearForm & imag() const { return *lfi_; }
+
+   void Update();
+   void Update(FiniteElementSpace *f);
+
+   /// Assembles the linear form i.e. sums over all domain/bdr integrators.
+   void Assemble();
+
+   std::complex<double> operator()(const ComplexGridFunction &gf) const;
+
+};
+
+// Class for sesquilinear form
+class SesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv_;
+
+   //protected:
+   BilinearForm *blfr_;
+   BilinearForm *blfi_;
+
+public:
+   SesquilinearForm(FiniteElementSpace *fes,
+                    ComplexOperator::Convention
+                    convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv_; }
+   void SetConvention(const ComplexOperator::Convention &
+                      convention) { conv_  = convention; }
+
+   BilinearForm & real() { return *blfr_; }
+   BilinearForm & imag() { return *blfi_; }
+   const BilinearForm & real() const { return *blfr_; }
+   const BilinearForm & imag() const { return *blfi_; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                            BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag,
+                              Array<int> &bdr_marker);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   /// Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
+   /** The returned matrix has to be deleted by the caller. */
+   ComplexSparseMatrix *AssembleCompSpMat();
+
+   /// Return the parallel FE space associated with the ParBilinearForm.
+   FiniteElementSpace *FESpace() const { return blfr_->FESpace(); }
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                         OperatorHandle &A, Vector &X, Vector &B,
+                         int copy_interior = 0);
+
+   /** Call this method after solving a linear system constructed using the
+       FormLinearSystem method to recover the solution as a ParGridFunction-size
+       vector in x. Use the same arguments as in the FormLinearSystem call. */
+   virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
+
+   virtual void Update(FiniteElementSpace *nfes = NULL);
+
+   virtual ~SesquilinearForm();
+};
+
+#ifdef MFEM_USE_MPI
+
+/// Class for complex-valued grid function - Vector with associated FE space.
+class ParComplexGridFunction : public Vector
+{
+private:
+
+   ParGridFunction * pgfr_;
+   ParGridFunction * pgfi_;
+
+protected:
+   void Destroy() { delete pgfr_; delete pgfi_; }
+
+public:
+
+   /* @brief Construct a ParComplexGridFunction associated with the
+      ParFiniteElementSpace @a *f. */
+   ParComplexGridFunction(ParFiniteElementSpace *pf);
+
+   void Update();
+
+   /// Assign constant values to the ParComplexGridFunction data.
+   ParComplexGridFunction &operator=(const std::complex<double> & value)
+   { *pgfr_ = value.real(); *pgfi_ = value.imag(); return *this; }
+
+   virtual void ProjectCoefficient(Coefficient &real_coeff,
+                                   Coefficient &imag_coeff);
+   virtual void ProjectCoefficient(VectorCoefficient &real_vcoeff,
+                                   VectorCoefficient &imag_vcoeff);
+
+   virtual void ProjectBdrCoefficient(Coefficient &real_coeff,
+                                      Coefficient &imag_coeff,
+                                      Array<int> &attr);
+   virtual void ProjectBdrCoefficientNormal(VectorCoefficient &real_coeff,
+                                            VectorCoefficient &imag_coeff,
+                                            Array<int> &attr);
+   virtual void ProjectBdrCoefficientTangent(VectorCoefficient &real_coeff,
+                                             VectorCoefficient &imag_coeff,
+                                             Array<int> &attr);
+
+   void Distribute(const Vector *tv);
+   void Distribute(const Vector &tv) { Distribute(&tv); }
+
+   /// Returns the vector restricted to the true dofs.
+   void ParallelProject(Vector &tv) const;
+
+   FiniteElementSpace *FESpace() { return pgfr_->FESpace(); }
+   const FiniteElementSpace *FESpace() const { return pgfr_->FESpace(); }
+
+   ParGridFunction & real() { return *pgfr_; }
+   ParGridFunction & imag() { return *pgfi_; }
+   const ParGridFunction & real() const { return *pgfr_; }
+   const ParGridFunction & imag() const { return *pgfi_; }
+
+   /// Destroys grid function.
+   virtual ~ParComplexGridFunction() { Destroy(); }
+
+};
+
+class ParComplexLinearForm : public Vector
+{
+private:
+   ComplexOperator::Convention conv_;
+
+protected:
+   ParLinearForm * plfr_;
+   ParLinearForm * plfi_;
+
+   HYPRE_Int * tdof_offsets_;
+
+public:
+
+   ParComplexLinearForm(ParFiniteElementSpace *pf,
+                        ComplexOperator::Convention
+                        convention = ComplexOperator::HERMITIAN);
+
+   virtual ~ParComplexLinearForm();
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(LinearFormIntegrator *lfi_real,
+                            LinearFormIntegrator *lfi_imag);
+
+   ParFiniteElementSpace *ParFESpace() const { return plfr_->ParFESpace(); }
+
+   ParLinearForm & real() { return *plfr_; }
+   ParLinearForm & imag() { return *plfi_; }
+   const ParLinearForm & real() const { return *plfr_; }
+   const ParLinearForm & imag() const { return *plfi_; }
+
+   void Update(ParFiniteElementSpace *pf = NULL);
+
+   /// Assembles the linear form i.e. sums over all domain/bdr integrators.
+   void Assemble();
+
+   /// Assemble the vector on the true dofs, i.e. P^t v.
+   void ParallelAssemble(Vector &tv);
+
+   /// Returns the vector assembled on the true dofs, i.e. P^t v.
+   HypreParVector *ParallelAssemble();
+
+   std::complex<double> operator()(const ParComplexGridFunction &gf) const;
+
+};
+
+// Class for parallel sesquilinear form
+class ParSesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv_;
+
+   //protected:
+   ParBilinearForm *pblfr_;
+   ParBilinearForm *pblfi_;
+
+public:
+   ParSesquilinearForm(ParFiniteElementSpace *pf,
+                       ComplexOperator::Convention
+                       convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv_; }
+   void SetConvention(const ComplexOperator::Convention &
+                      convention) { conv_  = convention; }
+
+   ParBilinearForm & real() { return *pblfr_; }
+   ParBilinearForm & imag() { return *pblfi_; }
+   const ParBilinearForm & real() const { return *pblfr_; }
+   const ParBilinearForm & imag() const { return *pblfi_; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                            BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator, restricted to specific boundary attributes.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag,
+                              Array<int> &bdr_marker);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   /// Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
+   /** The returned matrix has to be deleted by the caller. */
+   ComplexHypreParMatrix *ParallelAssemble();
+
+   /// Return the parallel FE space associated with the ParBilinearForm.
+   ParFiniteElementSpace *ParFESpace() const { return pblfr_->ParFESpace(); }
+
+   void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                         OperatorHandle &A, Vector &X, Vector &B,
+                         int copy_interior = 0);
+
+   /** Call this method after solving a linear system constructed using the
+       FormLinearSystem method to recover the solution as a ParGridFunction-size
+       vector in x. Use the same arguments as in the FormLinearSystem call. */
+   virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
+
+   virtual void Update(FiniteElementSpace *nfes = NULL);
+
+   virtual ~ParSesquilinearForm();
+};
+
+#endif // MFEM_USE_MPI
+
+}
+
+#endif // MFEM_COMPLEX_FEM

fem/fem.hpp

@@ -18,6 +18,7 @@
 #include "fe_coll.hpp"
 #include "eltrans.hpp"
 #include "coefficient.hpp"
+#include "complex_fem.hpp"
 #include "lininteg.hpp"
 #include "nonlininteg.hpp"
 #include "bilininteg.hpp"

fem/linearform.hpp

@@ -81,6 +81,15 @@ public:
        Update(FiniteElementSpace *, Vector &, int). */
    LinearForm() { fes = NULL; extern_lfs = 0; }
 
+   /// Construct a LinearForm using previously allocated array @a data.
+   /** The LinearForm does not assume ownership of @a data which is assumed to
+       be of size at least `f->GetVSize()`. Similar to the Vector constructor
+       for externally allocated array, the pointer @a data can be NULL. The data
+       array can be replaced later using the method SetData().
+    */
+   LinearForm(FiniteElementSpace *f, double *data) : Vector(data, f->GetVSize())
+   { fes = f; }
+
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
        have the same size.

fem/plinearform.hpp

@@ -45,6 +45,16 @@ public:
    /** The pointer @a pf is not owned by the newly constructed object. */
    ParLinearForm(ParFiniteElementSpace *pf) : LinearForm(pf) { pfes = pf; }
 
+   /// Construct a ParLinearForm using previously allocated array @a data.
+   /** The ParLinearForm does not assume ownership of @a data which is assumed
+       to be of size at least `pf->GetVSize()`. Similar to the LinearForm and
+       Vector constructors for externally allocated array, the pointer @a data
+       can be NULL. The data array can be replaced later using the method
+       SetData().
+    */
+   ParLinearForm(ParFiniteElementSpace *pf, double *data) :
+      LinearForm(pf, data), pfes(pf) { }
+
    /** @brief Create a ParLinearForm on the ParFiniteElementSpace @a *pf, using
        the same integrators as the ParLinearForm @a *plf.
 

fem/sidredatacollection.cpp

@@ -16,7 +16,9 @@
 
 #include "fem.hpp"
 
-#include <axom/sidre.hpp>
+#ifdef MFEM_USE_MPI
+#include <sidre/IOManager.hpp>
+#endif
 
 #include <string>
 #include <iomanip>      // for setw, setfill
@@ -202,10 +204,10 @@ SidreDataCollection::get_file_path(const std::string &filename) const
 
 axom::sidre::View *
 SidreDataCollection::AllocNamedBuffer(const std::string& buffer_name,
-                                      axom::sidre::IndexType sz,
+                                      axom::sidre::SidreLength sz,
                                       axom::sidre::TypeID type)
 {
-   sz = std::max(sz, sidre::IndexType(0));
+   sz = std::max(sz, sidre::SidreLength(0));
    sidre::Group *f = named_buffers_grp();
    sidre::View  *v = NULL;
 
@@ -823,7 +825,7 @@ void SidreDataCollection::Save(const std::string& filename,
 void SidreDataCollection::
 addScalarBasedGridFunction(const std::string &field_name, GridFunction *gf,
                            const std::string &buffer_name,
-                           axom::sidre::IndexType offset)
+                           axom::sidre::SidreLength offset)
 {
    sidre::Group* grp = m_bp_grp->getGroup("fields/" + field_name);
    MFEM_ASSERT(grp != NULL, "field " << field_name << " does not exist");
@@ -886,7 +888,7 @@ addScalarBasedGridFunction(const std::string &field_name, GridFunction *gf,
 void SidreDataCollection::
 addVectorBasedGridFunction(const std::string& field_name, GridFunction *gf,
                            const std::string &buffer_name,
-                           axom::sidre::IndexType offset)
+                           axom::sidre::SidreLength offset)
 {
    sidre::Group* grp = m_bp_grp->getGroup("fields/" + field_name);
    MFEM_ASSERT(grp != NULL, "field " << field_name << " does not exist");
@@ -1011,7 +1013,7 @@ DeregisterFieldInBPIndex(const std::string& field_name)
 void SidreDataCollection::RegisterField(const std::string &field_name,
                                         GridFunction *gf,
                                         const std::string &buffer_name,
-                                        axom::sidre::IndexType offset)
+                                        axom::sidre::SidreLength offset)
 {
    if ( field_name.empty() || buffer_name.empty() ||
         gf == NULL || gf->FESpace() == NULL )

fem/sidredatacollection.hpp

@@ -25,7 +25,7 @@
 #  pragma GCC diagnostic ignored "-Wpedantic"
 # endif
 #endif
-#include <axom/sidre.hpp>
+#include <sidre/sidre.hpp>
 #ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
 # pragma GCC diagnostic pop
 #endif
@@ -246,7 +246,7 @@ public:
     */
    void RegisterField(const std::string &field_name, GridFunction *gf,
                       const std::string &buffer_name,
-                      axom::sidre::IndexType offset);
+                      axom::sidre::SidreLength offset);
 
    /// Registers an attribute field in the Sidre DataStore
    /** The registration process is similar to that of RegisterField()
@@ -385,7 +385,7 @@ public:
     */
    axom::sidre::View *
    AllocNamedBuffer(const std::string& buffer_name,
-                    axom::sidre::IndexType sz,
+                    axom::sidre::SidreLength sz,
                     axom::sidre::TypeID type =
                        axom::sidre::DOUBLE_ID);
 
@@ -469,7 +469,7 @@ private:
    void addScalarBasedGridFunction(const std::string& field_name,
                                    GridFunction* gf,
                                    const std::string &buffer_name,
-                                   axom::sidre::IndexType offset);
+                                   axom::sidre::SidreLength offset);
 
    /**
     * \brief A private helper function to set up the views associated with the
@@ -483,7 +483,7 @@ private:
    void addVectorBasedGridFunction(const std::string& field_name,
                                    GridFunction* gf,
                                    const std::string &buffer_name,
-                                   axom::sidre::IndexType offset);
+                                   axom::sidre::SidreLength offset);
 
    /** @brief A private helper function to set up the Views associated with
        attribute field named @a field_name */

general/cuda.cpp

@@ -10,27 +10,17 @@
 // Software Foundation) version 2.1 dated February 1999.
 
 #include "cuda.hpp"
-#include "globals.hpp"
 
 namespace mfem
 {
 
-#ifdef MFEM_USE_CUDA
-void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
-                     const char *file, int line)
-{
-   mfem::err << "CUDA error: (" << expr << ") failed with error:\n --> "
-             << cudaGetErrorString(err)
-             << "\n ... in function: " << func
-             << "\n ... in file: " << file << ':' << line << '\n';
-   mfem_error();
-}
-#endif
-
 void* CuMemAlloc(void** dptr, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   MFEM_CUDA_CHECK(cudaMalloc(dptr, bytes));
+   if (CUDA_SUCCESS != ::cuMemAlloc((CUdeviceptr*)dptr, bytes))
+   {
+      mfem_error("Error in CuMemAlloc");
+   }
 #endif
    return *dptr;
 }
@@ -38,7 +28,10 @@ void* CuMemAlloc(void** dptr, size_t bytes)
 void* CuMemFree(void *dptr)
 {
 #ifdef MFEM_USE_CUDA
-   MFEM_CUDA_CHECK(cudaFree(dptr));
+   if (CUDA_SUCCESS != ::cuMemFree((CUdeviceptr)dptr))
+   {
+      mfem_error("Error in CuMemFree");
+   }
 #endif
    return dptr;
 }
@@ -46,15 +39,22 @@ void* CuMemFree(void *dptr)
 void* CuMemcpyHtoD(void* dst, const void* src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   MFEM_CUDA_CHECK(cudaMemcpy(dst, src, bytes, cudaMemcpyHostToDevice));
+   if (CUDA_SUCCESS != ::cuMemcpyHtoD((CUdeviceptr)dst, src, bytes))
+   {
+      mfem_error("Error in CuMemcpyHtoD");
+   }
 #endif
    return dst;
 }
 
-void* CuMemcpyHtoDAsync(void* dst, const void* src, size_t bytes)
+void* CuMemcpyHtoDAsync(void* dst, const void* src, size_t bytes, void *s)
 {
 #ifdef MFEM_USE_CUDA
-   MFEM_CUDA_CHECK(cudaMemcpyAsync(dst, src, bytes, cudaMemcpyHostToDevice));
+   if (CUDA_SUCCESS !=
+       ::cuMemcpyHtoDAsync((CUdeviceptr)dst, src, bytes, (CUstream)s))
+   {
+      mfem_error("Error in CuMemcpyHtoDAsync");
+   }
 #endif
    return dst;
 }
@@ -62,15 +62,24 @@ void* CuMemcpyHtoDAsync(void* dst, const void* src, size_t bytes)
 void* CuMemcpyDtoD(void* dst, void* src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   MFEM_CUDA_CHECK(cudaMemcpy(dst, src, bytes, cudaMemcpyDeviceToDevice));
+   if (CUDA_SUCCESS !=
+       ::cuMemcpyDtoD((CUdeviceptr)dst, (CUdeviceptr)src, bytes))
+   {
+      mfem_error("Error in CuMemcpyDtoD");
+   }
 #endif
    return dst;
 }
 
-void* CuMemcpyDtoDAsync(void* dst, void* src, size_t bytes)
+void* CuMemcpyDtoDAsync(void* dst, void* src, size_t bytes, void *s)
 {
 #ifdef MFEM_USE_CUDA
-   MFEM_CUDA_CHECK(cudaMemcpyAsync(dst, src, bytes, cudaMemcpyDeviceToDevice));
+   if (CUDA_SUCCESS !=
+       ::cuMemcpyDtoDAsync((CUdeviceptr)dst, (CUdeviceptr)src,
+                           bytes, (CUstream)s))
+   {
+      mfem_error("Error in CuMemcpyDtoDAsync");
+   }
 #endif
    return dst;
 }
@@ -78,7 +87,10 @@ void* CuMemcpyDtoDAsync(void* dst, void* src, size_t bytes)
 void* CuMemcpyDtoH(void *dst, void *src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA
-   MFEM_CUDA_CHECK(cudaMemcpy(dst, src, bytes, cudaMemcpyDeviceToHost));
+   if (CUDA_SUCCESS != ::cuMemcpyDtoH(dst, (CUdeviceptr)src, bytes))
+   {
+      mfem_error("Error in CuMemcpyDtoH");
+   }
 #endif
    return dst;
 }
@@ -86,7 +98,11 @@ void* CuMemcpyDtoH(void *dst, void *src, size_t bytes)
 void* CuMemcpyDtoHAsync(void* dst, void* src, size_t bytes, void *s)
 {
 #ifdef MFEM_USE_CUDA
-   MFEM_CUDA_CHECK(cudaMemcpyAsync(dst, src, bytes, cudaMemcpyDeviceToHost));
+   if (CUDA_SUCCESS !=
+       ::cuMemcpyDtoHAsync(dst, (CUdeviceptr)src, bytes, (CUstream)s))
+   {
+      mfem_error("Error in CuMemcpyDtoHAsync");
+   }
 #endif
    return dst;
 }

general/cuda.hpp

@@ -26,33 +26,89 @@
 #ifdef MFEM_USE_CUDA
 #define MFEM_ATTR_DEVICE __device__
 #define MFEM_ATTR_HOST_DEVICE __host__ __device__
-// Define a CUDA error check macro, MFEM_CUDA_CHECK(x), where x returns/is of
-// type 'cudaError_t'. This macro evaluates 'x' and raises an error if the
-// result is not cudaSuccess.
-#define MFEM_CUDA_CHECK(x) \
+// Define the CUDA debug macros:
+// - MFEM_CUDA_CHECK_DRV(x) where 'x' returns/is type 'CUresult'
+// - MFEM_CUDA_CHECK_RT(x)  where 'x' returns/is type 'cudaError_t'
+#ifdef MFEM_DEBUG
+#define MFEM_CUDA_CHECK_DRV(x) \
+   do \
+   { \
+      CUresult err = (x); \
+      if (err != CUDA_SUCCESS) \
+      { \
+         const char *error_string; \
+         cuGetErrorString(err, &error_string); \
+         _MFEM_MESSAGE("CUDA error: (" << #x \
+                       << ") failed with error:\n --> " \
+                       << error_string, 0); \
+      } \
+   } \
+   while (0)
+#define MFEM_CUDA_CHECK_RT(x) \
    do \
    { \
       cudaError_t err = (x); \
       if (err != cudaSuccess) \
       { \
-         mfem_cuda_error(err, #x, _MFEM_FUNC_NAME, __FILE__, __LINE__); \
+         _MFEM_MESSAGE("CUDA error: (" << #x \
+                       << ") failed with error:\n --> " \
+                       << cudaGetErrorString(err), 0); \
       } \
    } \
    while (0)
+#else
+#define MFEM_CUDA_CHECK_DRV(x) x
+#define MFEM_CUDA_CHECK_RT(x) x
+#endif
 #else // MFEM_USE_CUDA
 #define MFEM_ATTR_DEVICE
 #define MFEM_ATTR_HOST_DEVICE
+typedef int CUdevice;
+typedef int CUcontext;
+typedef void* CUstream;
 #endif // MFEM_USE_CUDA
 
 
 namespace mfem
 {
 
-#ifdef MFEM_USE_CUDA
-// Function used by the macro MFEM_CUDA_CHECK.
-void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
-                     const char *file, int line);
+// Define 'atomicAdd' function.
+#ifdef __CUDA_ARCH__
+#if __CUDA_ARCH__ < 600
+static __device__ inline double atomicAdd(double* address, double val)
+{
+   unsigned long long int* address_as_ull = (unsigned long long int*)address;
+   unsigned long long int old = *address_as_ull, assumed;
+   do
+   {
+      assumed = old;
+      old =
+         atomicCAS(address_as_ull, assumed,
+                   __double_as_longlong(val +
+                                        __longlong_as_double(assumed)));
+      // Note: uses integer comparison to avoid hang in case of NaN
+      // (since NaN != NaN)
+   }
+   while (assumed != old);
+   return __longlong_as_double(old);
+}
+#endif // __CUDA_ARCH__ < 600
+template<typename T> MFEM_ATTR_DEVICE
+inline T AtomicAdd(T volatile *address, T val)
+{
+   return atomicAdd((T *)address, val);
+}
+#else // __CUDA_ARCH__
+template<typename T> inline T AtomicAdd(T volatile *address, T val)
+{
+#ifdef MFEM_USE_OPENMP
+   #pragma omp atomic
 #endif
+   *address += val;
+   return *address;
+}
+#endif // __CUDA_ARCH__
+
 
 /// Allocates device memory
 void* CuMemAlloc(void **d_ptr, size_t bytes);
@@ -64,19 +120,20 @@ void* CuMemFree(void *d_ptr);
 void* CuMemcpyHtoD(void *d_dst, const void *h_src, size_t bytes);
 
 /// Copies memory from Host to Device
-void* CuMemcpyHtoDAsync(void *d_dst, const void *h_src, size_t bytes);
+void* CuMemcpyHtoDAsync(void *d_dst, const void *h_src,
+                        size_t bytes, void *stream);
 
 /// Copies memory from Device to Device
 void* CuMemcpyDtoD(void *d_dst, void *d_src, size_t bytes);
 
 /// Copies memory from Device to Device
-void* CuMemcpyDtoDAsync(void *d_dst, void *d_src, size_t bytes);
+void* CuMemcpyDtoDAsync(void *d_dst, void *d_src, size_t bytes, void *stream);
 
 /// Copies memory from Device to Host
 void* CuMemcpyDtoH(void *h_dst, void *d_src, size_t bytes);
 
 /// Copies memory from Device to Host
-void* CuMemcpyDtoHAsync(void *h_dst, void *d_src, size_t bytes);
+void* CuMemcpyDtoHAsync(void *h_dst, void *d_src, size_t bytes, void *stream);
 
 } // namespace mfem
 

general/device.cpp

@@ -24,6 +24,9 @@ namespace mfem
 namespace internal
 {
 
+CUstream *cuStream = NULL;
+static CUdevice cuDevice;
+static CUcontext cuContext;
 OccaDevice occaDevice;
 
 // Backends listed by priority, high to low:
@@ -97,9 +100,14 @@ void Device::Print(std::ostream &out)
 #ifdef MFEM_USE_CUDA
 static void DeviceSetup(const int dev, int &ngpu)
 {
-   MFEM_CUDA_CHECK(cudaGetDeviceCount(&ngpu));
-   MFEM_VERIFY(ngpu > 0, "No CUDA device found!");
-   MFEM_CUDA_CHECK(cudaSetDevice(dev));
+   cudaGetDeviceCount(&ngpu);
+   MFEM_VERIFY(ngpu>0, "No CUDA device found!");
+   cuInit(0);
+   cuDeviceGet(&internal::cuDevice, dev);
+   cuCtxCreate(&internal::cuContext, CU_CTX_SCHED_AUTO, internal::cuDevice);
+   internal::cuStream = new CUstream;
+   MFEM_VERIFY(internal::cuStream, "CUDA stream could not be created!");
+   cuStreamCreate(internal::cuStream, CU_STREAM_DEFAULT);
 }
 #endif
 
@@ -117,7 +125,7 @@ static void RajaDeviceSetup(const int dev, int &ngpu)
 #endif
 }
 
-static void OccaDeviceSetup(const int dev)
+static void OccaDeviceSetup(CUdevice cu_dev, CUcontext cu_ctx)
 {
 #ifdef MFEM_USE_OCCA
    const int cpu  = Device::Allows(Backend::OCCA_CPU);
@@ -130,8 +138,7 @@ static void OccaDeviceSetup(const int dev)
    if (cuda)
    {
 #if OCCA_CUDA_ENABLED
-      std::string mode("mode: 'CUDA', device_id : ");
-      internal::occaDevice.setup(mode.append(1,'0'+dev));
+      internal::occaDevice = occa::cuda::wrapDevice(cu_dev, cu_ctx);
 #else
       MFEM_ABORT("the OCCA CUDA backend requires OCCA built with CUDA!");
 #endif
@@ -190,10 +197,22 @@ void Device::Setup(const int device)
                "the OpenMP and RAJA OpenMP backends require MFEM built with"
                " MFEM_USE_OPENMP=YES");
 #endif
+   // The check for MFEM_USE_OCCA is in the function OccaDeviceSetup().
+
+   // We initialize CUDA and/or RAJA_CUDA first so OccaDeviceSetup() can reuse
+   // the same initialized cuDevice and cuContext objects when OCCA_CUDA is
+   // enabled.
    if (Allows(Backend::CUDA)) { CudaDeviceSetup(dev, ngpu); }
    if (Allows(Backend::RAJA_CUDA)) { RajaDeviceSetup(dev, ngpu); }
-   // The check for MFEM_USE_OCCA is in the function OccaDeviceSetup().
-   if (Allows(Backend::OCCA_MASK)) { OccaDeviceSetup(dev); }
+   if (Allows(Backend::OCCA_MASK))
+   {
+      OccaDeviceSetup(internal::cuDevice, internal::cuContext);
+   }
+}
+
+Device::~Device()
+{
+   delete internal::cuStream;
 }
 
 } // mfem

general/device.hpp

@@ -181,6 +181,8 @@ public:
        Backend::*_MASK, or combinations of those. */
    static inline bool Allows(unsigned long b_mask)
    { return Get().allowed_backends & b_mask; }
+
+   ~Device();
 };
 
 } // mfem

general/forall.hpp

@@ -22,9 +22,6 @@
 
 #ifdef MFEM_USE_RAJA
 #include "RAJA/RAJA.hpp"
-#if defined(RAJA_ENABLE_CUDA) && !defined(MFEM_USE_CUDA)
-#error When RAJA is built with CUDA, MFEM_USE_CUDA=YES is required
-#endif
 #endif
 
 namespace mfem
@@ -109,7 +106,8 @@ void CuWrap(const int N, DBODY &&d_body)
    if (N==0) { return; }
    const int GRID = (N+BLOCKS-1)/BLOCKS;
    CuKernel<<<GRID,BLOCKS>>>(N,d_body);
-   MFEM_CUDA_CHECK(cudaGetLastError());
+   const cudaError_t last = cudaGetLastError();
+   MFEM_VERIFY(last == cudaSuccess, cudaGetErrorString(last));
 }
 
 #else  // MFEM_USE_CUDA

general/mem_manager.cpp

@@ -312,6 +312,7 @@ void MemoryManager::Pull(const void *ptr, const std::size_t bytes)
    { mfem_error("Unknown pointer to pull from!"); }
 }
 
+namespace internal { extern CUstream *cuStream; }
 void* MemoryManager::Memcpy(void *dst, const void *src,
                             const std::size_t bytes, const bool async)
 {
@@ -321,7 +322,7 @@ void* MemoryManager::Memcpy(void *dst, const void *src,
    const bool run_on_host = !Device::Allows(Backend::DEVICE_MASK);
    if (run_on_host) { return std::memcpy(dst, src, bytes); }
    if (!async) { return CuMemcpyDtoD(d_dst, d_src, bytes); }
-   return CuMemcpyDtoDAsync(d_dst, d_src, bytes);
+   return CuMemcpyDtoDAsync(d_dst, d_src, bytes, internal::cuStream);
 }
 
 void MemoryManager::RegisterCheck(void *ptr)

general/occa.hpp

@@ -13,6 +13,7 @@
 #define MFEM_OCCA_HPP
 
 #include "../config/config.hpp"
+#include "cuda.hpp" // for CUdevice, CUcontext
 
 #ifdef MFEM_USE_OCCA
 #include <occa.hpp>

linalg/complex_operator.cpp

@@ -10,6 +10,8 @@
 // Software Foundation) version 2.1 dated February 1999.
 
 #include "complex_operator.hpp"
+#include <set>
+#include <map>
 
 namespace mfem
 {
@@ -39,6 +41,30 @@ ComplexOperator::~ComplexOperator()
    delete v_;
 }
 
+Operator & ComplexOperator::real()
+{
+   MFEM_ASSERT(Op_Real_, "ComplexOperator has no real part!");
+   return *Op_Real_;
+}
+
+Operator & ComplexOperator::imag()
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexOperator has no imaginary part!");
+   return *Op_Imag_;
+}
+
+const Operator & ComplexOperator::real() const
+{
+   MFEM_ASSERT(Op_Real_, "ComplexOperator has no real part!");
+   return *Op_Real_;
+}
+
+const Operator & ComplexOperator::imag() const
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexOperator has no imaginary part!");
+   return *Op_Imag_;
+}
+
 void ComplexOperator::Mult(const Vector &x, Vector &y) const
 {
    double * x_data = x.GetData();
@@ -120,6 +146,30 @@ void ComplexOperator::MultTranspose(const Vector &x_r, const Vector &x_i,
 }
 
 
+SparseMatrix & ComplexSparseMatrix::real()
+{
+   MFEM_ASSERT(Op_Real_, "ComplexSparseMatrix has no real part!");
+   return dynamic_cast<SparseMatrix &>(*Op_Real_);
+}
+
+SparseMatrix & ComplexSparseMatrix::imag()
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexSparseMatrix has no imaginary part!");
+   return dynamic_cast<SparseMatrix &>(*Op_Imag_);
+}
+
+const SparseMatrix & ComplexSparseMatrix::real() const
+{
+   MFEM_ASSERT(Op_Real_, "ComplexSparseMatrix has no real part!");
+   return dynamic_cast<const SparseMatrix &>(*Op_Real_);
+}
+
+const SparseMatrix & ComplexSparseMatrix::imag() const
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexSparseMatrix has no imaginary part!");
+   return dynamic_cast<const SparseMatrix &>(*Op_Imag_);
+}
+
 SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
 {
    SparseMatrix * A_r = dynamic_cast<SparseMatrix*>(Op_Real_);
@@ -184,4 +234,358 @@ SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
    return new SparseMatrix(I, J, D, this->Height(), this->Width());
 }
 
+#ifdef MFEM_USE_MPI
+
+ComplexHypreParMatrix::ComplexHypreParMatrix(HypreParMatrix * A_Real,
+                                             HypreParMatrix * A_Imag,
+                                             bool ownReal, bool ownImag,
+                                             Convention convention)
+   : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
+{
+   comm_ = (A_Real) ? A_Real->GetComm() :
+           ((A_Imag) ? A_Imag->GetComm() : MPI_COMM_WORLD);
+
+   MPI_Comm_rank(comm_, &myid_);
+   MPI_Comm_size(comm_, &nranks_);
+}
+
+HypreParMatrix & ComplexHypreParMatrix::real()
+{
+   MFEM_ASSERT(Op_Real_, "ComplexHypreParMatrix has no real part!");
+   return dynamic_cast<HypreParMatrix &>(*Op_Real_);
+}
+
+HypreParMatrix & ComplexHypreParMatrix::imag()
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexHypreParMatrix has no imaginary part!");
+   return dynamic_cast<HypreParMatrix &>(*Op_Imag_);
+}
+
+const HypreParMatrix & ComplexHypreParMatrix::real() const
+{
+   MFEM_ASSERT(Op_Real_, "ComplexHypreParMatrix has no real part!");
+   return dynamic_cast<const HypreParMatrix &>(*Op_Real_);
+}
+
+const HypreParMatrix & ComplexHypreParMatrix::imag() const
+{
+   MFEM_ASSERT(Op_Imag_, "ComplexHypreParMatrix has no imaginary part!");
+   return dynamic_cast<const HypreParMatrix &>(*Op_Imag_);
+}
+
+HypreParMatrix * ComplexHypreParMatrix::GetSystemMatrix() const
+{
+   HypreParMatrix * A_r = dynamic_cast<HypreParMatrix*>(Op_Real_);
+   HypreParMatrix * A_i = dynamic_cast<HypreParMatrix*>(Op_Imag_);
+
+   if ( A_r == NULL && A_i == NULL ) { return NULL; }
+
+   HYPRE_Int global_num_rows_r = (A_r) ? A_r->GetGlobalNumRows() : 0;
+   HYPRE_Int global_num_rows_i = (A_i) ? A_i->GetGlobalNumRows() : 0;
+   HYPRE_Int global_num_rows = std::max(global_num_rows_r, global_num_rows_i);
+
+   HYPRE_Int global_num_cols_r = (A_r) ? A_r->GetGlobalNumCols() : 0;
+   HYPRE_Int global_num_cols_i = (A_i) ? A_i->GetGlobalNumCols() : 0;
+   HYPRE_Int global_num_cols = std::max(global_num_cols_r, global_num_cols_i);
+
+   int row_starts_size = (HYPRE_AssumedPartitionCheck()) ? 2 : nranks_ + 1;
+   HYPRE_Int * row_starts = hypre_CTAlloc(HYPRE_Int, row_starts_size);
+   HYPRE_Int * col_starts = hypre_CTAlloc(HYPRE_Int, row_starts_size);
+
+   const HYPRE_Int * row_starts_z = (A_r) ? A_r->RowPart() :
+                                    ((A_i) ? A_i->RowPart() : NULL);
+   const HYPRE_Int * col_starts_z = (A_r) ? A_r->ColPart() :
+                                    ((A_i) ? A_i->ColPart() : NULL);
+
+   for (int i = 0; i < row_starts_size; i++)
+   {
+      row_starts[i] = 2 * row_starts_z[i];
+      col_starts[i] = 2 * col_starts_z[i];
+   }
+
+   SparseMatrix diag_r, diag_i, offd_r, offd_i;
+   HYPRE_Int * cmap_r, * cmap_i;
+
+   int nrows_r = 0, nrows_i = 0, ncols_r = 0, ncols_i = 0;
+   int ncols_offd_r = 0, ncols_offd_i = 0;
+   if (A_r)
+   {
+      A_r->GetDiag(diag_r);
+      A_r->GetOffd(offd_r, cmap_r);
+      nrows_r = diag_r.Height();
+      ncols_r = diag_r.Width();
+      ncols_offd_r = offd_r.Width();
+   }
+   if (A_i)
+   {
+      A_i->GetDiag(diag_i);
+      A_i->GetOffd(offd_i, cmap_i);
+      nrows_i = diag_i.Height();
+      ncols_i = diag_i.Width();
+      ncols_offd_i = offd_i.Width();
+   }
+   int nrows = std::max(nrows_r, nrows_i);
+   int ncols = std::max(ncols_r, ncols_i);
+
+   // Determine the unique set of off-diagonal columns global indices
+   std::set<int> cset;
+   for (int i=0; i<ncols_offd_r; i++)
+   {
+      cset.insert(cmap_r[i]);
+   }
+   for (int i=0; i<ncols_offd_i; i++)
+   {
+      cset.insert(cmap_i[i]);
+   }
+   int num_cols_offd = (int)cset.size();
+
+   // Exatract pointers to the various CSR arrays of the diagonal blocks
+   const int * diag_r_I = (A_r) ? diag_r.GetI() : NULL;
+   const int * diag_i_I = (A_i) ? diag_i.GetI() : NULL;
+
+   const int * diag_r_J = (A_r) ? diag_r.GetJ() : NULL;
+   const int * diag_i_J = (A_i) ? diag_i.GetJ() : NULL;
+
+   const double * diag_r_D = (A_r) ? diag_r.GetData() : NULL;
+   const double * diag_i_D = (A_i) ? diag_i.GetData() : NULL;
+
+   int diag_r_nnz = (diag_r_I) ? diag_r_I[nrows] : 0;
+   int diag_i_nnz = (diag_i_I) ? diag_i_I[nrows] : 0;
+   int diag_nnz = 2 * (diag_r_nnz + diag_i_nnz);
+
+   // Exatract pointers to the various CSR arrays of the off-diagonal blocks
+   const int * offd_r_I = (A_r) ? offd_r.GetI() : NULL;
+   const int * offd_i_I = (A_i) ? offd_i.GetI() : NULL;
+
+   const int * offd_r_J = (A_r) ? offd_r.GetJ() : NULL;
+   const int * offd_i_J = (A_i) ? offd_i.GetJ() : NULL;
+
+   const double * offd_r_D = (A_r) ? offd_r.GetData() : NULL;
+   const double * offd_i_D = (A_i) ? offd_i.GetData() : NULL;
+
+   int offd_r_nnz = (offd_r_I) ? offd_r_I[nrows] : 0;
+   int offd_i_nnz = (offd_i_I) ? offd_i_I[nrows] : 0;
+   int offd_nnz = 2 * (offd_r_nnz + offd_i_nnz);
+
+   // Allocate CSR arrays for the combined matrix
+   HYPRE_Int * diag_I = hypre_CTAlloc(HYPRE_Int, 2 * nrows + 1);
+   HYPRE_Int * diag_J = hypre_CTAlloc(HYPRE_Int, diag_nnz);
+   double    * diag_D = hypre_CTAlloc(double, diag_nnz);
+
+   HYPRE_Int * offd_I = hypre_CTAlloc(HYPRE_Int, 2 * nrows + 1);
+   HYPRE_Int * offd_J = hypre_CTAlloc(HYPRE_Int, offd_nnz);
+   double    * offd_D = hypre_CTAlloc(double, offd_nnz);
+   HYPRE_Int * cmap   = hypre_CTAlloc(HYPRE_Int, 2 * num_cols_offd);
+
+   // Fill the CSR arrays for the diagonal portion of the matrix
+   const double factor = (convention_ == HERMITIAN) ? 1.0 : -1.0;
+
+   diag_I[0] = 0;
+   diag_I[nrows] = diag_r_nnz + diag_i_nnz;
+   for (int i=0; i<nrows; i++)
+   {
+      diag_I[i + 1]         = ((diag_r_I)?diag_r_I[i+1]:0) +
+                              ((diag_i_I)?diag_i_I[i+1]:0);
+      diag_I[i + nrows + 1] = diag_I[i+1] + diag_r_nnz + diag_i_nnz;
+
+      if (diag_r_I)
+      {
+         for (int j=0; j<diag_r_I[i+1] - diag_r_I[i]; j++)
+         {
+            diag_J[diag_I[i] + j] = diag_r_J[diag_r_I[i] + j];
+            diag_D[diag_I[i] + j] = diag_r_D[diag_r_I[i] + j];
+
+            diag_J[diag_I[i+nrows] + j] =
+               diag_r_J[diag_r_I[i] + j] + ncols;
+            diag_D[diag_I[i+nrows] + j] =
+               factor * diag_r_D[diag_r_I[i] + j];
+         }
+      }
+      if (diag_i_I)
+      {
+         const int off_r = (diag_r_I)?(diag_r_I[i+1] - diag_r_I[i]):0;
+         for (int j=0; j<diag_i_I[i+1] - diag_i_I[i]; j++)
+         {
+            diag_J[diag_I[i] + off_r + j] =  diag_i_J[diag_i_I[i] + j] + ncols;
+            diag_D[diag_I[i] + off_r + j] = -diag_i_D[diag_i_I[i] + j];
+
+            diag_J[diag_I[i+nrows] + off_r + j] = diag_i_J[diag_i_I[i] + j];
+            diag_D[diag_I[i+nrows] + off_r + j] =
+               factor * diag_i_D[diag_i_I[i] + j];
+         }
+      }
+   }
+
+   // Determine the mappings describing the layout of off-diagonal columns
+   int num_recv_procs = 0;
+   HYPRE_Int * offd_col_start_stop = NULL;
+   this->getColStartStop(A_r, A_i, num_recv_procs, offd_col_start_stop);
+
+   std::set<int>::iterator sit;
+   std::map<int,int> cmapa, cmapb, cinvmap;
+   for (sit=cset.begin(); sit!=cset.end(); sit++)
+   {
+      int col_orig = *sit;
+      int col_2x2  = -1;
+      int col_size = 0;
+      for (int i=0; i<num_recv_procs; i++)
+      {
+         if (offd_col_start_stop[2*i] <= col_orig &&
+             col_orig < offd_col_start_stop[2*i+1])
+         {
+            col_2x2 = offd_col_start_stop[2*i] + col_orig;
+            col_size = offd_col_start_stop[2*i+1] - offd_col_start_stop[2*i];
+            break;
+         }
+      }
+      cmapa[*sit] = col_2x2;
+      cmapb[*sit] = col_2x2 + col_size;
+      cinvmap[col_2x2] = -1;
+      cinvmap[col_2x2 + col_size] = -1;
+   }
+   delete [] offd_col_start_stop;
+
+   std::map<int, int>::iterator mit;
+   int i = 0;
+   for (mit=cinvmap.begin(); mit!=cinvmap.end(); mit++, i++)
+   {
+      mit->second = i;
+      cmap[i] = mit->first;
+   }
+
+   // Fill the CSR arrays for the off-diagonal portion of the matrix
+   offd_I[0] = 0;
+   offd_I[nrows] = offd_r_nnz + offd_i_nnz;
+   for (int i=0; i<nrows; i++)
+   {
+      offd_I[i + 1]         = ((offd_r_I)?offd_r_I[i+1]:0) +
+                              ((offd_i_I)?offd_i_I[i+1]:0);
+      offd_I[i + nrows + 1] = offd_I[i+1] + offd_r_nnz + offd_i_nnz;
+
+      if (offd_r_I)
+      {
+         const int off_i = (offd_i_I)?(offd_i_I[i+1] - offd_i_I[i]):0;
+         for (int j=0; j<offd_r_I[i+1] - offd_r_I[i]; j++)
+         {
+            offd_J[offd_I[i] + j] =
+               cinvmap[cmapa[cmap_r[offd_r_J[offd_r_I[i] + j]]]];
+            offd_D[offd_I[i] + j] = offd_r_D[offd_r_I[i] + j];
+
+            offd_J[offd_I[i+nrows] + off_i + j] =
+               cinvmap[cmapb[cmap_r[offd_r_J[offd_r_I[i] + j]]]];
+            offd_D[offd_I[i+nrows] + off_i + j] =
+               factor * offd_r_D[offd_r_I[i] + j];
+         }
+      }
+      if (offd_i_I)
+      {
+         const int off_r = (offd_r_I)?(offd_r_I[i+1] - offd_r_I[i]):0;
+         for (int j=0; j<offd_i_I[i+1] - offd_i_I[i]; j++)
+         {
+            offd_J[offd_I[i] + off_r + j] =
+               cinvmap[cmapb[cmap_i[offd_i_J[offd_i_I[i] + j]]]];
+            offd_D[offd_I[i] + off_r + j] = -offd_i_D[offd_i_I[i] + j];
+
+            offd_J[offd_I[i+nrows] + j] =
+               cinvmap[cmapa[cmap_i[offd_i_J[offd_i_I[i] + j]]]];
+            offd_D[offd_I[i+nrows] + j] = factor * offd_i_D[offd_i_I[i] + j];
+         }
+      }
+   }
+
+   // Construct the combined matrix
+   HypreParMatrix * A = new HypreParMatrix(comm_,
+                                           2 * global_num_rows,
+                                           2 * global_num_cols,
+                                           row_starts, col_starts,
+                                           diag_I, diag_J, diag_D,
+                                           offd_I, offd_J, offd_D,
+                                           2 * num_cols_offd, cmap);
+
+   // Give the new matrix ownership of its interanl arrays
+   A->SetOwnerFlags(-1,-1,-1);
+   hypre_CSRMatrixSetDataOwner(((hypre_ParCSRMatrix*)(*A))->diag,1);
+   hypre_CSRMatrixSetDataOwner(((hypre_ParCSRMatrix*)(*A))->offd,1);
+   hypre_ParCSRMatrixSetRowStartsOwner((hypre_ParCSRMatrix*)(*A),1);
+   hypre_ParCSRMatrixSetColStartsOwner((hypre_ParCSRMatrix*)(*A),1);
+
+   return A;
+}
+
+void
+ComplexHypreParMatrix::getColStartStop(const HypreParMatrix * A_r,
+                                       const HypreParMatrix * A_i,
+                                       int & num_recv_procs,
+                                       HYPRE_Int *& offd_col_start_stop) const
+{
+   hypre_ParCSRCommPkg * comm_pkg_r =
+      (A_r) ? hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix*)(*A_r)) : NULL;
+   hypre_ParCSRCommPkg * comm_pkg_i =
+      (A_i) ? hypre_ParCSRMatrixCommPkg((hypre_ParCSRMatrix*)(*A_i)) : NULL;
+
+   std::set<HYPRE_Int> send_procs, recv_procs;
+   if ( comm_pkg_r )
+   {
+      for (HYPRE_Int i=0; i<comm_pkg_r->num_sends; i++)
+      {
+         send_procs.insert(comm_pkg_r->send_procs[i]);
+      }
+      for (HYPRE_Int i=0; i<comm_pkg_r->num_recvs; i++)
+      {
+         recv_procs.insert(comm_pkg_r->recv_procs[i]);
+      }
+   }
+   if ( comm_pkg_i )
+   {
+      for (HYPRE_Int i=0; i<comm_pkg_i->num_sends; i++)
+      {
+         send_procs.insert(comm_pkg_i->send_procs[i]);
+      }
+      for (HYPRE_Int i=0; i<comm_pkg_i->num_recvs; i++)
+      {
+         recv_procs.insert(comm_pkg_i->recv_procs[i]);
+      }
+   }
+
+   num_recv_procs = (int)recv_procs.size();
+
+   HYPRE_Int loc_start_stop[2];
+   offd_col_start_stop = new HYPRE_Int[2 * num_recv_procs];
+
+   const HYPRE_Int * row_part = (A_r) ? A_r->RowPart() :
+                                ((A_i) ? A_i->RowPart() : NULL);
+
+   int row_part_ind = (HYPRE_AssumedPartitionCheck()) ? 0 : myid_;
+   loc_start_stop[0] = row_part[row_part_ind];
+   loc_start_stop[1] = row_part[row_part_ind+1];
+
+   MPI_Request * req = new MPI_Request[send_procs.size()+recv_procs.size()];
+   MPI_Status * stat = new MPI_Status[send_procs.size()+recv_procs.size()];
+   int send_count = 0;
+   int recv_count = 0;
+   int tag = 0;
+
+   std::set<HYPRE_Int>::iterator sit;
+   for (sit=send_procs.begin(); sit!=send_procs.end(); sit++)
+   {
+      MPI_Isend(loc_start_stop, 2, HYPRE_MPI_INT,
+                *sit, tag, comm_, &req[send_count]);
+      send_count++;
+   }
+   for (sit=recv_procs.begin(); sit!=recv_procs.end(); sit++)
+   {
+      MPI_Irecv(&offd_col_start_stop[2*recv_count], 2, HYPRE_MPI_INT,
+                *sit, tag, comm_, &req[send_count+recv_count]);
+      recv_count++;
+   }
+
+   MPI_Waitall(send_count+recv_count, req, stat);
+
+   delete [] req;
+   delete [] stat;
+}
+
+
+#endif // MFEM_USE_MPI
+
 }

linalg/complex_operator.hpp

@@ -14,6 +14,9 @@
 
 #include "operator.hpp"
 #include "sparsemat.hpp"
+#ifdef MFEM_USE_MPI
+#include "hypre.hpp"
+#endif
 
 namespace mfem
 {
@@ -27,7 +30,8 @@ namespace mfem
     ComplexOperator allows one to choose a convention upon construction, which
     facilitates symmetry.
 
-    Matrix-vector products are then computed as:
+    If we let (y_r + i y_i) = (Op_r + i Op_i)(x_r + i x_i) then Matrix-vector
+    products are then computed as:
 
     1. When Convention::HERMITIAN is used (default)
     / y_r \   / Op_r -Op_i \ / x_r \
@@ -38,6 +42,8 @@ namespace mfem
     / y_r \   / Op_r -Op_i \ / x_r \
     |     | = |            | |     |
     \-y_i /   \-Op_i -Op_r / \ x_i /
+    In other words, Matrix-vector products with Convention::BLOCK_SYMMETRIC
+    compute the complex conjugate of Op*x.
 
     Either convention can be used with a given complex operator,
     however, each of them is best suited for certain classes of
@@ -82,9 +88,30 @@ public:
 
    virtual ~ComplexOperator();
 
+   /** @brief Check for existence of real or imaginary part of the operator
+
+       These methods do not check that the operators are non-zero but
+       only that the operators have been set.
+    */
+   bool hasRealPart() const { return Op_Real_ != NULL; }
+   bool hasImagPart() const { return Op_Imag_ != NULL; }
+
+   /** @brief Real or imaginary part accessor methods
+
+       The following accessor methods should only be called if the
+       requested part of the opertor is known to exist.  This
+       can be checked with hasRealPart() or hasImagPart().
+   */
+   virtual Operator & real();
+   virtual Operator & imag();
+   virtual const Operator & real() const;
+   virtual const Operator & imag() const;
+
    virtual void Mult(const Vector &x, Vector &y) const;
    virtual void MultTranspose(const Vector &x, Vector &y) const;
 
+   virtual Type GetType() const { return Complex_Operator; }
+
 protected:
    // Let this be hidden from the public interface since the implementation
    // depends on internal members
@@ -127,9 +154,74 @@ public:
       : ComplexOperator(A_Real, A_Imag, ownReal, ownImag, convention)
    {}
 
+   virtual SparseMatrix & real();
+   virtual SparseMatrix & imag();
+
+   virtual const SparseMatrix & real() const;
+   virtual const SparseMatrix & imag() const;
+
+   /** Combine the blocks making up this complex operator into a
+       single SparseMatrix.  The resulting matrix can be passed to
+       solvers which require access to the matrix entries themselves,
+       such as sparse direct solvers, rather than simply the action of
+       the opertor.  Note that this combined operator requires roughly
+       twice the memory of the block structured operator. */
    SparseMatrix * GetSystemMatrix() const;
+
+   virtual Type GetType() const { return MFEM_ComplexSparseMat; }
 };
 
+#ifdef MFEM_USE_MPI
+
+/** @brief Specialization of the ComplexOperator built from a pair of
+    HypreParMatrices.
+
+    The purpose of this specialization is to construct a single
+    HypreParMatrix object which is equivalent to the 2x2 block system
+    that the ComplexOperator mimics. The resulting HypreParMatrix can
+    then be passed along to solvers which require access to the CSR
+    matrix data such as SuperLU, STRUMPACK, or similar sparse linear
+    solvers.
+
+    See ComplexOperator documentation in operator.hpp for more information.
+ */
+class ComplexHypreParMatrix : public ComplexOperator
+{
+public:
+   ComplexHypreParMatrix(HypreParMatrix * A_Real, HypreParMatrix * A_Imag,
+                         bool ownReal, bool ownImag,
+                         Convention convention = HERMITIAN);
+
+   virtual HypreParMatrix & real();
+   virtual HypreParMatrix & imag();
+
+   virtual const HypreParMatrix & real() const;
+   virtual const HypreParMatrix & imag() const;
+
+   /** Combine the blocks making up this complex operator into a
+       single HypreParMatrix.  The resulting matrix can be passed to
+       solvers which require access to the matrix entries themselves,
+       such as sparse direct solvers or Hypre preconditioners, rather
+       than simply the action of the opertor.  Note that this combined
+       operator requires roughly twice the memory of the block
+       structured operator. */
+   HypreParMatrix * GetSystemMatrix() const;
+
+   virtual Type GetType() const { return Complex_Hypre_ParCSR; }
+
+private:
+   void getColStartStop(const HypreParMatrix * A_r,
+                        const HypreParMatrix * A_i,
+                        int & num_recv_procs,
+                        HYPRE_Int *& offd_col_start_stop) const;
+
+   MPI_Comm comm_;
+   int myid_;
+   int nranks_;
+};
+
+#endif // MFEM_USE_MPI
+
 }
 
-#endif
+#endif // MFEM_COMPLEX_OPERATOR

linalg/operator.hpp

@@ -124,14 +124,17 @@ public:
    enum Type
    {
       ANY_TYPE,         ///< ID for the base class Operator, i.e. any type.
-      MFEM_SPARSEMAT,   ///< ID for class SparseMatrix
+      MFEM_SPARSEMAT,   ///< ID for class SparseMatrix.
       Hypre_ParCSR,     ///< ID for class HypreParMatrix.
       PETSC_MATAIJ,     ///< ID for class PetscParMatrix, MATAIJ format.
       PETSC_MATIS,      ///< ID for class PetscParMatrix, MATIS format.
       PETSC_MATSHELL,   ///< ID for class PetscParMatrix, MATSHELL format.
       PETSC_MATNEST,    ///< ID for class PetscParMatrix, MATNEST format.
       PETSC_MATHYPRE,   ///< ID for class PetscParMatrix, MATHYPRE format.
-      PETSC_MATGENERIC  ///< ID for class PetscParMatrix, unspecified format.
+      PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
+      Complex_Operator, ///< ID for class ComplexOperator.
+      MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
+      Complex_Hypre_ParCSR   ///< ID for class ComplexHypreParMatrix.
    };
 
    /// Return the type ID of the Operator class.
@@ -302,6 +305,24 @@ public:
 };
 
 
+/// Scaled Operator B: x -> a A(x).
+class ScaledOperator : public Operator
+{
+private:
+   const Operator &A_;
+   double a_;
+
+public:
+   /// Create a scalar product operator related to A.
+   explicit ScaledOperator(const Operator *A, double a)
+      : Operator(A->Width(), A->Height()), A_(*A), a_(a) { }
+
+   /// Operator application
+   virtual void Mult(const Vector &x, Vector &y) const
+   { A_.Mult(x, y); y *= a_; }
+};
+
+
 /** @brief The transpose of a given operator. Switches the roles of the methods
     Mult() and MultTranspose(). */
 class TransposeOperator : public Operator

linalg/sparsemat.cpp

@@ -38,7 +38,6 @@ SparseMatrix::SparseMatrix(int nrows, int ncols)
      current_row(-1),
      ColPtrJ(NULL),
      ColPtrNode(NULL),
-     At(NULL),
      ownGraph(true),
      ownData(true),
      isSorted(false)
@@ -61,7 +60,6 @@ SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n)
      Rows(NULL),
      ColPtrJ(NULL),
      ColPtrNode(NULL),
-     At(NULL),
      ownGraph(true),
      ownData(true),
      isSorted(false)
@@ -80,7 +78,6 @@ SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n,
      Rows(NULL),
      ColPtrJ(NULL),
      ColPtrNode(NULL),
-     At(NULL),
      ownGraph(ownij),
      ownData(owna),
      isSorted(issorted)
@@ -106,7 +103,6 @@ SparseMatrix::SparseMatrix(int nrows, int ncols, int rowsize)
    , Rows(NULL)
    , ColPtrJ(NULL)
    , ColPtrNode(NULL)
-   , At(NULL)
    , ownGraph(true)
    , ownData(true)
    , isSorted(false)
@@ -187,7 +183,6 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
    current_row = -1;
    ColPtrJ = NULL;
    ColPtrNode = NULL;
-   At = NULL;
    isSorted = mat.isSorted;
 }
 
@@ -196,7 +191,6 @@ SparseMatrix::SparseMatrix(const Vector &v)
    , Rows(NULL)
    , ColPtrJ(NULL)
    , ColPtrNode(NULL)
-   , At(NULL)
    , ownGraph(true)
    , ownData(true)
    , isSorted(true)
@@ -251,7 +245,6 @@ void SparseMatrix::SetEmpty()
    current_row = -1;
    ColPtrJ = NULL;
    ColPtrNode = NULL;
-   At = NULL;
 #ifdef MFEM_USE_MEMALLOC
    NodesMem = NULL;
 #endif
@@ -340,7 +333,7 @@ void SparseMatrix::SetWidth(int newWidth)
       // Nothing to be done here
       return;
    }
-   else if (newWidth == -1)
+   else if ( newWidth == -1)
    {
       // Compute the actual width
       width = ActualWidth();
@@ -559,10 +552,12 @@ void SparseMatrix::Mult(const Vector &x, Vector &y) const
 
 void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
 {
-   MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
-               << ") must match matrix width (" << width << ")");
-   MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
-               << ") must match matrix height (" << height << ")");
+   MFEM_ASSERT(width == x.Size(),
+               "Input vector size (" << x.Size() << ") must match matrix width (" << width
+               << ")");
+   MFEM_ASSERT(height == y.Size(),
+               "Output vector size (" << y.Size() << ") must match matrix height (" << height
+               << ")");
 
    int i, j, end;
    double *Ap = A, *yp = y.GetData();
@@ -641,10 +636,12 @@ void SparseMatrix::MultTranspose(const Vector &x, Vector &y) const
 void SparseMatrix::AddMultTranspose(const Vector &x, Vector &y,
                                     const double a) const
 {
-   MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
-               << ") must match matrix height (" << height << ")");
-   MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()
-               << ") must match matrix width (" << width << ")");
+   MFEM_ASSERT(height == x.Size(),
+               "Input vector size (" << x.Size() << ") must match matrix height (" << height
+               << ")");
+   MFEM_ASSERT(width == y.Size(),
+               "Output vector size (" << y.Size() << ") must match matrix width (" << width
+               << ")");
 
    if (A == NULL)
    {
@@ -661,40 +658,23 @@ void SparseMatrix::AddMultTranspose(const Vector &x, Vector &y,
       }
       return;
    }
-
-   if (At)
+   // Prepare the lambda capture and get our pointers from the memory manager
+   const int d_height = height;
+   const DeviceArray d_I(I);
+   const DeviceArray d_J(J);
+   const DeviceVector d_A(A);
+   const DeviceVector d_x(x, x.Size());
+   DeviceVector d_y(y, y.Size());
+   MFEM_FORALL(i, d_height,
    {
-      At->AddMult(x, y, a);
-   }
-   else
-   {
-      MFEM_VERIFY(Device::IsDisabled(), "transpose action on device is not "
-                  "enabled; see BuildTranspose() for details.");
-      for (int i = 0; i < height; i++)
+      const double xi = a * d_x[i];
+      const int end = d_I[i+1];
+      for (int j = d_I[i]; j < end; j++)
       {
-         const double xi = a * x[i];
-         const int end = I[i+1];
-         for (int j = I[i]; j < end; j++)
-         {
-            const int Jj = J[j];
-            y[Jj] += A[j] * xi;
-         }
+         const int Jj = d_J[j];
+         AtomicAdd(&d_y[Jj], d_A[j] * xi);
       }
-   }
-}
-
-void SparseMatrix::BuildTranspose() const
-{
-   if (At == NULL)
-   {
-      At = Transpose(*this);
-   }
-}
-
-void SparseMatrix::ResetTranspose() const
-{
-   delete At;
-   At = NULL;
+   });
 }
 
 void SparseMatrix::PartMult(
@@ -2121,12 +2101,12 @@ void SparseMatrix::Set(const int i, const int j, const double A)
    if ((gi=i) < 0) { gi = -1-gi, s = -1; }
    else { s = 1; }
    MFEM_ASSERT(gi < height,
-               "Trying to set a row " << gi << " outside the matrix height "
+               "Trying to insert a row " << gi << " outside the matrix height "
                << height);
    if ((gj=j) < 0) { gj = -1-gj, t = -s; }
    else { t = s; }
    MFEM_ASSERT(gj < width,
-               "Trying to set a column " << gj << " outside the matrix width "
+               "Trying to insert a column " << gj << " outside the matrix width "
                << width);
    if (t < 0) { a = -a; }
    _Set_(gi, gj, a);
@@ -2162,7 +2142,7 @@ void SparseMatrix::SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
       if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
       else { s = 1; }
       MFEM_ASSERT(gi < height,
-                  "Trying to set a row " << gi << " outside the matrix height "
+                  "Trying to insert a row " << gi << " outside the matrix height "
                   << height);
       SetColPtr(gi);
       for (j = 0; j < cols.Size(); j++)
@@ -2175,7 +2155,7 @@ void SparseMatrix::SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to set a column " << gj << " outside the matrix width "
+                     "Trying to insert a column " << gj << " outside the matrix width "
                      << width);
          if (t < 0) { a = -a; }
          _Set_(gj, a);
@@ -2197,7 +2177,7 @@ void SparseMatrix::SetSubMatrixTranspose(const Array<int> &rows,
       if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
       else { s = 1; }
       MFEM_ASSERT(gi < height,
-                  "Trying to set a row " << gi << " outside the matrix height "
+                  "Trying to insert a row " << gi << " outside the matrix height "
                   << height);
       SetColPtr(gi);
       for (j = 0; j < cols.Size(); j++)
@@ -2210,7 +2190,7 @@ void SparseMatrix::SetSubMatrixTranspose(const Array<int> &rows,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to set a column " << gj << " outside the matrix width "
+                     "Trying to insert a column " << gj << " outside the matrix width "
                      << width);
          if (t < 0) { a = -a; }
          _Set_(gj, a);
@@ -2230,7 +2210,7 @@ void SparseMatrix::GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
       if ((gi=rows[i]) < 0) { gi = -1-gi, s = -1; }
       else { s = 1; }
       MFEM_ASSERT(gi < height,
-                  "Trying to read a row " << gi << " outside the matrix height "
+                  "Trying to insert a row " << gi << " outside the matrix height "
                   << height);
       SetColPtr(gi);
       for (j = 0; j < cols.Size(); j++)
@@ -2238,7 +2218,7 @@ void SparseMatrix::GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to read a column " << gj << " outside the matrix width "
+                     "Trying to insert a column " << gj << " outside the matrix width "
                      << width);
          a = _Get_(gj);
          subm(i, j) = (t < 0) ? (-a) : (a);
@@ -2256,7 +2236,7 @@ bool SparseMatrix::RowIsEmpty(const int row) const
       gi = -1-gi;
    }
    MFEM_ASSERT(gi < height,
-               "Trying to query a row " << gi << " outside the matrix height "
+               "Trying to insert a row " << gi << " outside the matrix height "
                << height);
    if (Rows)
    {
@@ -2275,7 +2255,7 @@ int SparseMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
 
    if ((gi=row) < 0) { gi = -1-gi; }
    MFEM_ASSERT(gi < height,
-               "Trying to read a row " << gi << " outside the matrix height "
+               "Trying to insert a row " << gi << " outside the matrix height "
                << height);
    if (Rows)
    {
@@ -2302,7 +2282,7 @@ int SparseMatrix::GetRow(const int row, Array<int> &cols, Vector &srow) const
       j = I[gi];
       cols.MakeRef(J + j, I[gi+1]-j);
       srow.NewDataAndSize(A + j, cols.Size());
-      MFEM_ASSERT(row >= 0, "Row not valid: " << row << ", height: " << height);
+      MFEM_ASSERT(row >= 0, "Row not valid: " << row );
       return 1;
    }
 }
@@ -2316,7 +2296,7 @@ void SparseMatrix::SetRow(const int row, const Array<int> &cols,
    if ((gi=row) < 0) { gi = -1-gi, s = -1; }
    else { s = 1; }
    MFEM_ASSERT(gi < height,
-               "Trying to set a row " << gi << " outside the matrix height "
+               "Trying to insert a row " << gi << " outside the matrix height "
                << height);
 
    if (!Finalized())
@@ -2327,7 +2307,7 @@ void SparseMatrix::SetRow(const int row, const Array<int> &cols,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to set a column " << gj << " outside the matrix"
+                     "Trying to insert a column " << gj << " outside the matrix"
                      " width " << width);
          a = srow(j);
          if (t < 0) { a = -a; }
@@ -2345,7 +2325,7 @@ void SparseMatrix::SetRow(const int row, const Array<int> &cols,
          if ((gj=cols[j]) < 0) { gj = -1-gj, t = -s; }
          else { t = s; }
          MFEM_ASSERT(gj < width,
-                     "Trying to set a column " << gj << " outside the matrix"
+                     "Trying to insert a column " << gj << " outside the matrix"
                      " width " << width);
 
          J[i] = gj;
@@ -2791,10 +2771,9 @@ void SparseMatrix::Destroy()
       delete NodesMem;
    }
 #endif
-   delete At;
 }
 
-int SparseMatrix::ActualWidth() const
+int SparseMatrix::ActualWidth()
 {
    int awidth = 0;
    if (A)
@@ -2838,10 +2817,8 @@ SparseMatrix *Transpose (const SparseMatrix &A)
       "Finalize must be called before Transpose. Use TransposeRowMatrix instead");
 
    int i, j, end;
-   const int *A_i, *A_j;
-   int m, n, nnz, *At_i, *At_j;
-   const double *A_data;
-   double *At_data;
+   int m, n, nnz, *A_i, *A_j, *At_i, *At_j;
+   double *A_data, *At_data;
 
    m      = A.Height(); // number of rows of A
    n      = A.Width();  // number of columns of A
@@ -2968,10 +2945,8 @@ SparseMatrix *Mult (const SparseMatrix &A, const SparseMatrix &B,
                     SparseMatrix *OAB)
 {
    int nrowsA, ncolsA, nrowsB, ncolsB;
-   const int *A_i, *A_j, *B_i, *B_j;
-   int *C_i, *C_j, *B_marker;
-   const double *A_data, *B_data;
-   double *C_data;
+   int *A_i, *A_j, *B_i, *B_j, *C_i, *C_j, *B_marker;
+   double *A_data, *B_data, *C_data;
    int ia, ib, ic, ja, jb, num_nonzeros;
    int row_start, counter;
    double a_entry, b_entry;
@@ -3284,13 +3259,13 @@ SparseMatrix * Add(double a, const SparseMatrix & A, double b,
    int * C_j;
    double * C_data;
 
-   const int *A_i = A.GetI();
-   const int *A_j = A.GetJ();
-   const double *A_data = A.GetData();
+   int * A_i = A.GetI();
+   int * A_j = A.GetJ();
+   double * A_data = A.GetData();
 
-   const int *B_i = B.GetI();
-   const int *B_j = B.GetJ();
-   const double *B_data = B.GetData();
+   int * B_i = B.GetI();
+   int * B_j = B.GetJ();
+   double * B_data = B.GetData();
 
    int * marker = new int[ncols];
    std::fill(marker, marker+ncols, -1);
@@ -3518,7 +3493,6 @@ void SparseMatrix::Swap(SparseMatrix &other)
    mfem::Swap(current_row, other.current_row);
    mfem::Swap(ColPtrJ, other.ColPtrJ);
    mfem::Swap(ColPtrNode, other.ColPtrNode);
-   mfem::Swap(At, other.At);
 
 #ifdef MFEM_USE_MEMALLOC
    mfem::Swap(NodesMem, other.NodesMem);

linalg/sparsemat.hpp

@@ -64,9 +64,6 @@ protected:
    mutable int* ColPtrJ;
    mutable RowNode ** ColPtrNode;
 
-   /// Transpose of A. Owned. Used to perform MultTranspose() on devices.
-   mutable SparseMatrix *At;
-
 #ifdef MFEM_USE_MEMALLOC
    typedef MemAlloc <RowNode, 1024> RowNodeAlloc;
    RowNodeAlloc * NodesMem;
@@ -100,9 +97,6 @@ public:
 
    /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
        @a data is optionally transferred to the SparseMatrix. */
-   /** If the parameter @a data is NULL, then the internal #A array is allocated
-       by this constructor (initializing it with zeros and taking ownership,
-       regardless of the parameter @a owna). */
    SparseMatrix(int *i, int *j, double *data, int m, int n, bool ownij,
                 bool owna, bool issorted);
 
@@ -118,7 +112,7 @@ public:
        ownership. */
    SparseMatrix(const SparseMatrix &mat, bool copy_graph = true);
 
-   /// Create a SparseMatrix with diagonal @a v, i.e. A = Diag(v)
+   /// Create a SparseMatrix with diagonal v, i.e. A = Diag(v)
    SparseMatrix(const Vector & v);
 
 
@@ -140,35 +134,21 @@ public:
    /// Check if the SparseMatrix is empty.
    bool Empty() const { return (A == NULL) && (Rows == NULL); }
 
-   /// Return the array #I.
-   inline int *GetI() { return I; }
-   /// Return the array #I, const version.
-   inline const int *GetI() const { return I; }
-
-   /// Return the array #J.
-   inline int *GetJ() { return J; }
-   /// Return the array #J, const version.
-   inline const int *GetJ() const { return J; }
-
-   /// Return the element data, i.e. the array #A.
-   inline double *GetData() { return A; }
-   /// Return the element data, i.e. the array #A, const version.
-   inline const double *GetData() const { return A; }
-
-   /// Returns the number of elements in row @a i.
+   /// Return the array #I
+   inline int *GetI() const { return I; }
+   /// Return the array #J
+   inline int *GetJ() const { return J; }
+   /// Return element data, i.e. array #A
+   inline double *GetData() const { return A; }
+   /// Returns the number of elements in row @a i
    int RowSize(const int i) const;
-
-   /// Returns the maximum number of elements among all rows.
+   /// Returns the maximum number of elements among all rows
    int MaxRowSize() const;
-
-   /// Return a pointer to the column indices in a row.
+   /// Return a pointer to the column indices in a row
    int *GetRowColumns(const int row);
-   /// Return a pointer to the column indices in a row, const version.
    const int *GetRowColumns(const int row) const;
-
-   /// Return a pointer to the entries in a row.
+   /// Return a pointer to the entries in a row
    double *GetRowEntries(const int row);
-   /// Return a pointer to the entries in a row, const version.
    const double *GetRowEntries(const int row) const;
 
    /// Change the width of a SparseMatrix.
@@ -183,7 +163,7 @@ public:
 
    /// Returns the actual Width of the matrix.
    /*! This method can be called for matrices finalized or not. */
-   int ActualWidth() const;
+   int ActualWidth();
 
    /// Sort the column indices corresponding to each row.
    void SortColumnIndices();
@@ -226,45 +206,13 @@ public:
    void AddMultTranspose(const Vector &x, Vector &y,
                          const double a = 1.0) const;
 
-   /** @brief Build and store internally the transpose of this matrix which will
-       be used in the methods AddMultTranspose() and MultTranspose(). */
-   /** If this method has been called, the internal transpose matrix will be
-       used to perform the action of the transpose matrix in AddMultTranspose(),
-       and MultTranspose().
-
-       Warning: any changes in this matrix will invalidate the internal
-       transpose. To rebuild the transpose, call ResetTranspose() followed by a
-       call to this method. If the internal transpose is already built, this
-       method has no effect.
-
-       When any non-default backend is enabled, i.e. Device::IsEnabled() is
-       true, the methods AddMultTranspose(), and MultTranspose(), require the
-       internal transpose to be built. If that is not the case (i.e. the
-       internal transpose is not built), these methods will raise an error with
-       an appropriate message pointing to this method. When using the default
-       backend, calling this method is optional.
-
-       This method can only be used when the sparse matrix is finalized. */
-   void BuildTranspose() const;
-
-   /** Reset (destroy) the internal transpose matrix. See BuildTranspose() for
-       more details. */
-   void ResetTranspose() const;
-
    void PartMult(const Array<int> &rows, const Vector &x, Vector &y) const;
    void PartAddMult(const Array<int> &rows, const Vector &x, Vector &y,
                     const double a=1.0) const;
 
-   /// y = A * x, treating all entries as booleans (zero=false, nonzero=true).
-   /** The actual values stored in the data array, #A, are not used - this means
-       and that all entries in the sparsity pattern are considered to be true by
-       this method. */
+   /// y = A * x, but treat all elements as booleans (zero=false, nonzero=true).
    void BooleanMult(const Array<int> &x, Array<int> &y) const;
-
-   /// y = At * x, treating all entries as booleans (zero=false, nonzero=true).
-   /** The actual values stored in the data array, #A, are not used - this means
-       and that all entries in the sparsity pattern are considered to be true by
-       this method. */
+   /// y = At * x, but treat all elements as booleans (zero=false, nonzero=true).
    void BooleanMultTranspose(const Array<int> &x, Array<int> &y) const;
 
    /// Compute y^t A x
@@ -504,11 +452,11 @@ SparseMatrix *Transpose(const SparseMatrix &A);
 SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
                                              int useActualWidth);
 
-/// Matrix product A.B.
-/** If @a OAB is not NULL, we assume it has the structure of A.B and store the
-    result in @a OAB. If @a OAB is NULL, we create a new SparseMatrix to store
+/** Matrix product A.B.
+    If OAB is not NULL, we assume it has the structure
+    of A.B and store the result in OAB.
+    If OAB is NULL, we create a new SparseMatrix to store
     the result and return a pointer to it.
-
     All matrices must be finalized. */
 SparseMatrix *Mult(const SparseMatrix &A, const SparseMatrix &B,
                    SparseMatrix *OAB = NULL);
@@ -607,20 +555,16 @@ inline void SparseMatrix::SetColPtr(const int row) const
 inline void SparseMatrix::ClearColPtr() const
 {
    if (Rows)
-   {
       for (RowNode *node_p = Rows[current_row]; node_p != NULL;
            node_p = node_p->Prev)
       {
          ColPtrNode[node_p->Column] = NULL;
       }
-   }
    else
-   {
       for (int j = I[current_row], end = I[current_row+1]; j < end; j++)
       {
          ColPtrJ[J[j]] = -1;
       }
-   }
 }
 
 inline double &SparseMatrix::SearchRow(const int col)

linalg/vector.cpp

@@ -857,11 +857,11 @@ static double cuVectorMin(const int N, const double *X)
    const int bytes = min_sz*sizeof(double);
    static double *h_min = NULL;
    if (!h_min) { h_min = (double*)calloc(min_sz,sizeof(double)); }
-   static void *gdsr = NULL;
-   if (!gdsr) { MFEM_CUDA_CHECK(cudaMalloc(&gdsr, bytes)); }
+   static CUdeviceptr gdsr = (CUdeviceptr) NULL;
+   if (!gdsr) { ::cuMemAlloc(&gdsr,bytes); }
    cuKernelMin<<<gridSize,blockSize>>>(N, (double*)gdsr, x);
-   MFEM_CUDA_CHECK(cudaGetLastError());
-   MFEM_CUDA_CHECK(cudaMemcpy(h_min, gdsr, bytes, cudaMemcpyDeviceToHost));
+   MFEM_CUDA_CHECK_RT(cudaGetLastError());
+   ::cuMemcpy((CUdeviceptr)h_min,(CUdeviceptr)gdsr,bytes);
    double min = std::numeric_limits<double>::infinity();
    for (int i = 0; i < min_sz; i++) { min = fmin(min, h_min[i]); }
    return min;
@@ -909,19 +909,19 @@ static double cuVectorDot(const int N, const double *X, const double *Y)
       if (h_dot) { free(h_dot); }
       h_dot = (double*)calloc(dot_sz,sizeof(double));
    }
-   static void *gdsr = NULL;
+   static CUdeviceptr gdsr = (CUdeviceptr) NULL;
    if (!gdsr or dot_block_sz!=dot_sz)
    {
-      if (gdsr) { MFEM_CUDA_CHECK(cudaFree(gdsr)); }
-      MFEM_CUDA_CHECK(cudaMalloc(&gdsr,bytes));
+      if (gdsr) { MFEM_CUDA_CHECK_DRV(::cuMemFree(gdsr)); }
+      MFEM_CUDA_CHECK_DRV(::cuMemAlloc(&gdsr,bytes));
    }
    if (dot_block_sz!=dot_sz)
    {
       dot_block_sz = dot_sz;
    }
    cuKernelDot<<<gridSize,blockSize>>>(N, (double*)gdsr, x, y);
-   MFEM_CUDA_CHECK(cudaGetLastError());
-   MFEM_CUDA_CHECK(cudaMemcpy(h_dot, gdsr, bytes, cudaMemcpyDeviceToHost));
+   MFEM_CUDA_CHECK_RT(cudaGetLastError());
+   MFEM_CUDA_CHECK_DRV(::cuMemcpy((CUdeviceptr)h_dot,(CUdeviceptr)gdsr,bytes));
    double dot = 0.0;
    for (int i = 0; i < dot_sz; i++) { dot += h_dot[i]; }
    return dot;

makefile

@@ -211,9 +211,6 @@ ifneq ($(MFEM_USE_CUDA),YES)
    XCOMPILER = $(CXX_XCOMPILER)
    XLINKER   = $(CXX_XLINKER)
 else
-   ifneq ($(MFEM_USE_MM),YES)
-      $(error MFEM_USE_CUDA=YES requires MFEM_USE_MM=YES)
-   endif
    MFEM_CXX ?= $(CUDA_CXX)
    CXXFLAGS += $(CUDA_FLAGS) -ccbin $(CXX_OR_MPICXX)
    XCOMPILER = $(CUDA_XCOMPILER)
@@ -233,7 +230,7 @@ endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
- SUITESPARSE STRUMPACK GECKO GNUTLS NETCDF PETSC MPFR PUMI OCCA RAJA
+ SUITESPARSE STRUMPACK GECKO GNUTLS NETCDF PETSC MPFR PUMI CUDA OCCA RAJA
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 
 define mfem_check_dependency
@@ -249,7 +246,7 @@ ifeq ($(MAKECMDGOALS),config)
 endif
 
 # List of MFEM dependencies, processed below
-MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) LIBUNWIND OPENMP CUDA
+MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) LIBUNWIND OPENMP
 
 # List of deprecated MFEM dependencies, processed below
 MFEM_LEGACY_DEPENDENCIES = OPENMP
@@ -322,8 +319,8 @@ MFEM_TEST_MK   ?= @MFEM_DIR@/config/test.mk
 # Use "\n" (interpreted by sed) to add a newline.
 MFEM_CONFIG_EXTRA ?= $(if $(BUILD_DIR_DEF),MFEM_BUILD_DIR ?= @MFEM_DIR@,)
 
-MFEM_SOURCE_DIR  = $(MFEM_REAL_DIR)
-MFEM_INSTALL_DIR = $(abspath $(MFEM_PREFIX))
+MFEM_SOURCE_DIR  := $(MFEM_REAL_DIR)
+MFEM_INSTALL_DIR := $(BUILD_REAL_DIR)
 
 # If we have 'config' target, export variables used by config/makefile
 ifneq (,$(filter config,$(MAKECMDGOALS)))
@@ -347,11 +344,6 @@ ifneq (,$(filter install,$(MAKECMDGOALS)))
    MFEM_LIBS     = $(if $(shared),$(INSTALL_RPATH)) -L@MFEM_LIB_DIR@ -lmfem\
       @MFEM_EXT_LIBS@
    MFEM_LIB_FILE = @MFEM_LIB_DIR@/libmfem.$(if $(shared),$(SO_VER),a)
-   ifeq ($(MFEM_USE_OCCA),YES)
-      ifneq ($(MFEM_INSTALL_DIR),$(abspath $(PREFIX))
-         $(error OCCA is enabled: PREFIX must be set during configuration!)
-      endif
-   endif
    MFEM_PREFIX := $(abspath $(PREFIX))
    MFEM_INC_DIR = $(abspath $(PREFIX_INC))
    MFEM_LIB_DIR = $(abspath $(PREFIX_LIB))
@@ -366,7 +358,6 @@ DIRS = general linalg mesh fem
 SOURCE_FILES = $(foreach dir,$(DIRS),$(wildcard $(SRC)$(dir)/*.cpp))
 RELSRC_FILES = $(patsubst $(SRC)%,%,$(SOURCE_FILES))
 OBJECT_FILES = $(patsubst $(SRC)%,$(BLD)%,$(SOURCE_FILES:.cpp=.o))
-OKL_DIRS = fem
 
 .PHONY: lib all clean distclean install config status info deps serial parallel\
 	debug pdebug cuda pcuda cudebug pcudebug style check test unittest\
@@ -507,12 +498,7 @@ install: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT))
 # install remaining includes in each subdirectory
 	for dir in $(DIRS); do \
 	   mkdir -p $(PREFIX_INC)/mfem/$$dir && \
-	   $(INSTALL) -m 640 $(SRC)$$dir/*.hpp $(PREFIX_INC)/mfem/$$dir; \
-	done
-# install *.okl files
-	for dir in $(OKL_DIRS); do \
-	   mkdir -p $(PREFIX_INC)/mfem/$$dir && \
-	   $(INSTALL) -m 640 $(SRC)$$dir/*.okl $(PREFIX_INC)/mfem/$$dir; \
+	   $(INSTALL) -m 640 $(SRC)$$dir/*.hpp $(SRC)$$dir/*.okl $(PREFIX_INC)/mfem/$$dir; \
 	done
 # install config.mk in $(PREFIX_SHARE)
 	mkdir -p $(PREFIX_SHARE)

mesh/pmesh.cpp

@@ -3090,23 +3090,6 @@ void ParMesh::Rebalance()
                  " meshes.");
    }
 
-   // Make sure the Nodes use a ParFiniteElementSpace
-   if (Nodes && dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace()) == NULL)
-   {
-      ParFiniteElementSpace *pfes =
-         new ParFiniteElementSpace(*Nodes->FESpace(), *this);
-      ParGridFunction *new_nodes = new ParGridFunction(pfes);
-      *new_nodes = *Nodes;
-      if (Nodes->OwnFEC())
-      {
-         new_nodes->MakeOwner(Nodes->OwnFEC());
-         Nodes->MakeOwner(NULL); // takes away ownership of 'fec' and 'fes'
-         delete Nodes->FESpace();
-      }
-      delete Nodes;
-      Nodes = new_nodes;
-   }
-
    DeleteFaceNbrData();
 
    pncmesh->Rebalance();
@@ -3127,6 +3110,22 @@ void ParMesh::Rebalance()
    last_operation = Mesh::REBALANCE;
    sequence++;
 
+   // Make sure the Nodes use a ParFiniteElementSpace
+   if (Nodes && dynamic_cast<ParFiniteElementSpace*>(Nodes->FESpace()) == NULL)
+   {
+      ParFiniteElementSpace *pfes =
+         new ParFiniteElementSpace(*Nodes->FESpace(), *this);
+      ParGridFunction *new_nodes = new ParGridFunction(pfes);
+      *new_nodes = *Nodes;
+      if (Nodes->OwnFEC())
+      {
+         new_nodes->MakeOwner(Nodes->OwnFEC());
+         Nodes->MakeOwner(NULL); // takes away ownership of 'fec' and 'fes'
+         delete Nodes->FESpace();
+      }
+      delete Nodes;
+      Nodes = new_nodes;
+   }
    UpdateNodes();
 }
 

tests/unit/CMakeLists.txt

@@ -9,11 +9,11 @@
 # terms of the GNU Lesser General Public License (as published by the Free
 # Software Foundation) version 2.1 dated February 1999.
 
+# Include the build directory where mfem.hpp and mfem-performance.hpp are.
+include_directories(BEFORE ${PROJECT_BINARY_DIR})
 # Include the top mfem source directory - needed by some tests, e.g. to
 # #include "general/text.hpp".
 include_directories(BEFORE ${PROJECT_SOURCE_DIR})
-# Include the build directory where mfem.hpp and mfem-performance.hpp are.
-include_directories(BEFORE ${PROJECT_BINARY_DIR})
 # Include the source directory for the unit tests - catch.hpp is there.
 include_directories(BEFORE ${CMAKE_CURRENT_SOURCE_DIR})
 

tests/unit/fem/test_datacollection.cpp

@@ -12,13 +12,7 @@
 #include "mfem.hpp"
 #include "catch.hpp"
 #include <stdio.h>
-
-#ifndef _WIN32
-#include <unistd.h> // rmdir
-#else
-#include <direct.h> // _rmdir
-#define rmdir(dir) _rmdir(dir)
-#endif
+#include <unistd.h>  // rmdir
 
 using namespace mfem;
 

tests/unit/fem/test_inversetransform.cpp

@@ -125,7 +125,7 @@ TEST_CASE("InverseElementTransformation",
       REQUIRE( mesh_file.good() );
 
       const int npts = 100; // number of random points to test
-      const int min_found_pts = 93;
+      const int min_found_pts = 94;
       const int rand_seed = 189548;
       srand(rand_seed);