add lshape example

Add AMD HIP support to CMake build system

CHANGELOG

@@ -35,9 +35,11 @@ Version 4.2.1 (development)
 - Implemented a filter method for the Navier miniapp to stabilize highly
   turbulent flows in direct numerical simulation.
 
+- Added HIP support to the CMake build system.
+
 - Added support for reading high-order Lagrange meshes in VTK format. Arbitrary-
   orders and all element types are supported. See the VTK blog for more info:
-  https://blog.kitware.com/wp-content/uploads/2018/09/Source_Issue_43.pdf
+  https://blog.kitware.com/wp-content/uploads/2018/09/Source_Issue_43.pdf.
 
 - Added support for reading VTK meshes in XML format.
 
@@ -66,6 +68,8 @@ Version 4.2.1 (development)
 - Implemented an adaptive linear solver tolerance option for NewtonSolver based
   on the algorithm of Eisenstat and Walker.
 
+- Added support for nonscalar coefficient with VectorDiffusionIntegrator.
+
 - Extending support for L2 basis functions using MapTypes VALUE and INTEGRAL in
   linear interpolators and GridFunction "GetValue" methods.
 

CMakeLists.txt

@@ -88,6 +88,9 @@ include("${CMAKE_CURRENT_SOURCE_DIR}/config/XSDKDefaults.cmake")
 # Enable languages.
 enable_language(CXX)
 if (MFEM_USE_CUDA)
+   if (MFEM_USE_HIP)
+      message(FATAL_ERROR " *** MFEM_USE_HIP cannot be combined with MFEM_USE_CUDA.")
+   endif()
    # MFEM_USE_CUDA requires CMake 3.8 or newer (for direct CUDA support)
    cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
    # Use ${CMAKE_CXX_COMPILER} as the cuda host compiler.
@@ -357,6 +360,15 @@ if (MFEM_USE_UMPIRE)
    find_package(UMPIRE REQUIRED)
 endif()
 
+# AMD HIP
+if (MFEM_USE_HIP)
+  find_package(HIP REQUIRED)
+  if (HIP_ARCH)
+    message(STATUS "Using HIP architecture: ${HIP_ARCH}")
+    list(APPEND HIP_HIPCC_FLAGS "--amdgpu-target=${HIP_ARCH}")
+  endif()
+endif()
+
 # ADIOS2 for parallel I/O
 if (MFEM_USE_ADIOS2)
    find_package(ADIOS2 REQUIRED)
@@ -432,9 +444,9 @@ foreach(DIR IN LISTS MFEM_SOURCE_DIRS)
 endforeach()
 
 if (MFEM_USE_CUDA)
-   foreach(file IN LISTS SOURCES)
-      set_property(SOURCE ${file} PROPERTY LANGUAGE CUDA)
-   endforeach()
+  set_source_files_properties(${SOURCES} PROPERTIES LANGUAGE CUDA)
+elseif(MFEM_USE_HIP)
+  set_source_files_properties(${SOURCES} PROPERTIES HIP_SOURCE_PROPERTY_FORMAT TRUE)
 endif()
 
 add_subdirectory(config)
@@ -453,7 +465,7 @@ set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX} CACHE PATH
     "The MFEM install directory" FORCE)
 
 # Declaring the library
-add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
+mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
 if (CMAKE_VERSION VERSION_GREATER 2.8.11)
   target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES})

config/cmake/config.hpp.in

@@ -132,6 +132,10 @@
 // Requires a CUDA compiler (nvcc).
 #cmakedefine MFEM_USE_CUDA
 
+// Build the HIP-enabled version of the MFEM library.
+// Requires a HIP compiler (hipcc).
+#cmakedefine MFEM_USE_HIP
+
 // Enable MFEM functionality based on the RAJA library
 #cmakedefine MFEM_USE_RAJA
 

config/cmake/modules/FindHIP.cmake

@@ -0,0 +1,692 @@
+###############################################################################
+# FindHIP.cmake
+###############################################################################
+include(CheckCXXCompilerFlag)
+###############################################################################
+# SET: Variable defaults
+###############################################################################
+# User defined flags
+set(HIP_HIPCC_FLAGS "" CACHE STRING "Semicolon delimited flags for HIPCC")
+set(HIP_HCC_FLAGS "" CACHE STRING "Semicolon delimited flags for HCC")
+set(HIP_CLANG_FLAGS "" CACHE STRING "Semicolon delimited flags for CLANG")
+set(HIP_NVCC_FLAGS "" CACHE STRING "Semicolon delimted flags for NVCC")
+mark_as_advanced(HIP_HIPCC_FLAGS HIP_HCC_FLAGS HIP_CLANG_FLAGS HIP_NVCC_FLAGS)
+
+set(_hip_configuration_types ${CMAKE_CONFIGURATION_TYPES} ${CMAKE_BUILD_TYPE} Debug MinSizeRel Release RelWithDebInfo)
+list(REMOVE_DUPLICATES _hip_configuration_types)
+foreach(config ${_hip_configuration_types})
+    string(TOUPPER ${config} config_upper)
+    set(HIP_HIPCC_FLAGS_${config_upper} "" CACHE STRING "Semicolon delimited flags for HIPCC")
+    set(HIP_HCC_FLAGS_${config_upper} "" CACHE STRING "Semicolon delimited flags for HCC")
+    set(HIP_CLANG_FLAGS_${config_upper} "" CACHE STRING "Semicolon delimited flags for CLANG")
+    set(HIP_NVCC_FLAGS_${config_upper} "" CACHE STRING "Semicolon delimited flags for NVCC")
+    mark_as_advanced(HIP_HIPCC_FLAGS_${config_upper} HIP_HCC_FLAGS_${config_upper} HIP_CLANG_FLAGS_${config_upper} HIP_NVCC_FLAGS_${config_upper})
+endforeach()
+option(HIP_HOST_COMPILATION_CPP "Host code compilation mode" ON)
+option(HIP_VERBOSE_BUILD "Print out the commands run while compiling the HIP source file.  With the Makefile generator this defaults to VERBOSE variable specified on the command line, but can be forced on with this option." OFF)
+mark_as_advanced(HIP_HOST_COMPILATION_CPP)
+
+###############################################################################
+# FIND: HIP and associated helper binaries
+###############################################################################
+
+get_filename_component(_IMPORT_PREFIX "${CMAKE_CURRENT_LIST_DIR}/../" REALPATH)
+
+# HIP is supported on Linux only
+if(UNIX AND NOT APPLE AND NOT CYGWIN)
+    # Search for HIP installation
+    if(NOT HIP_ROOT_DIR)
+        # Search in user specified path first
+        find_path(
+            HIP_ROOT_DIR
+            NAMES bin/hipconfig
+            PATHS
+            "$ENV{ROCM_PATH}/hip"
+            ENV HIP_PATH
+            ${_IMPORT_PREFIX}
+            /opt/rocm/hip
+            DOC "HIP installed location"
+            NO_DEFAULT_PATH
+            )
+        if(NOT EXISTS ${HIP_ROOT_DIR})
+            if(HIP_FIND_REQUIRED)
+                message(FATAL_ERROR "Specify HIP_ROOT_DIR")
+            elseif(NOT HIP_FIND_QUIETLY)
+                message("HIP_ROOT_DIR not found or specified")
+            endif()
+        endif()
+        # And push it back to the cache
+        set(HIP_ROOT_DIR ${HIP_ROOT_DIR} CACHE PATH "HIP installed location" FORCE)
+    endif()
+
+    # Find HIPCC executable
+    find_program(
+        HIP_HIPCC_EXECUTABLE
+        NAMES hipcc
+        PATHS
+        "${HIP_ROOT_DIR}"
+        ENV ROCM_PATH
+        ENV HIP_PATH
+        /opt/rocm
+        /opt/rocm/hip
+        PATH_SUFFIXES bin
+        NO_DEFAULT_PATH
+        )
+    if(NOT HIP_HIPCC_EXECUTABLE)
+        # Now search in default paths
+        find_program(HIP_HIPCC_EXECUTABLE hipcc)
+    endif()
+    mark_as_advanced(HIP_HIPCC_EXECUTABLE)
+
+    # Find HIPCONFIG executable
+    find_program(
+        HIP_HIPCONFIG_EXECUTABLE
+        NAMES hipconfig
+        PATHS
+        "${HIP_ROOT_DIR}"
+        ENV ROCM_PATH
+        ENV HIP_PATH
+        /opt/rocm
+        /opt/rocm/hip
+        PATH_SUFFIXES bin
+        NO_DEFAULT_PATH
+        )
+    if(NOT HIP_HIPCONFIG_EXECUTABLE)
+        # Now search in default paths
+        find_program(HIP_HIPCONFIG_EXECUTABLE hipconfig)
+    endif()
+    mark_as_advanced(HIP_HIPCONFIG_EXECUTABLE)
+
+    # Find HIPCC_CMAKE_LINKER_HELPER executable
+    find_program(
+        HIP_HIPCC_CMAKE_LINKER_HELPER
+        NAMES hipcc_cmake_linker_helper
+        PATHS
+        "${HIP_ROOT_DIR}"
+        ENV ROCM_PATH
+        ENV HIP_PATH
+        /opt/rocm
+        /opt/rocm/hip
+        PATH_SUFFIXES bin
+        NO_DEFAULT_PATH
+        )
+    if(NOT HIP_HIPCC_CMAKE_LINKER_HELPER)
+        # Now search in default paths
+        find_program(HIP_HIPCC_CMAKE_LINKER_HELPER hipcc_cmake_linker_helper)
+    endif()
+    mark_as_advanced(HIP_HIPCC_CMAKE_LINKER_HELPER)
+
+    if(HIP_HIPCONFIG_EXECUTABLE AND NOT HIP_VERSION)
+        # Compute the version
+        execute_process(
+            COMMAND ${HIP_HIPCONFIG_EXECUTABLE} --version
+            OUTPUT_VARIABLE _hip_version
+            ERROR_VARIABLE _hip_error
+            OUTPUT_STRIP_TRAILING_WHITESPACE
+            ERROR_STRIP_TRAILING_WHITESPACE
+            )
+        if(NOT _hip_error)
+            set(HIP_VERSION ${_hip_version} CACHE STRING "Version of HIP as computed from hipcc")
+        else()
+            set(HIP_VERSION "0.0.0" CACHE STRING "Version of HIP as computed by FindHIP()")
+        endif()
+        mark_as_advanced(HIP_VERSION)
+    endif()
+    if(HIP_VERSION)
+        string(REPLACE "." ";" _hip_version_list "${HIP_VERSION}")
+        list(GET _hip_version_list 0 HIP_VERSION_MAJOR)
+        list(GET _hip_version_list 1 HIP_VERSION_MINOR)
+        list(GET _hip_version_list 2 HIP_VERSION_PATCH)
+        set(HIP_VERSION_STRING "${HIP_VERSION}")
+    endif()
+
+    if(HIP_HIPCONFIG_EXECUTABLE AND NOT HIP_PLATFORM)
+        # Compute the platform
+        execute_process(
+            COMMAND ${HIP_HIPCONFIG_EXECUTABLE} --platform
+            OUTPUT_VARIABLE _hip_platform
+            OUTPUT_STRIP_TRAILING_WHITESPACE
+            )
+        set(HIP_PLATFORM ${_hip_platform} CACHE STRING "HIP platform as computed by hipconfig")
+        mark_as_advanced(HIP_PLATFORM)
+    endif()
+
+    if(HIP_HIPCONFIG_EXECUTABLE AND NOT HIP_COMPILER)
+        # Compute the compiler
+        execute_process(
+            COMMAND ${HIP_HIPCONFIG_EXECUTABLE} --compiler
+            OUTPUT_VARIABLE _hip_compiler
+            OUTPUT_STRIP_TRAILING_WHITESPACE
+            )
+        set(HIP_COMPILER ${_hip_compiler} CACHE STRING "HIP compiler as computed by hipconfig")
+        mark_as_advanced(HIP_COMPILER)
+    endif()
+
+    if(HIP_HIPCONFIG_EXECUTABLE AND NOT HIP_RUNTIME)
+        # Compute the runtime
+        execute_process(
+            COMMAND ${HIP_HIPCONFIG_EXECUTABLE} --runtime
+            OUTPUT_VARIABLE _hip_runtime
+            OUTPUT_STRIP_TRAILING_WHITESPACE
+            )
+        set(HIP_RUNTIME ${_hip_runtime} CACHE STRING "HIP runtime as computed by hipconfig")
+        mark_as_advanced(HIP_RUNTIME)
+    endif()
+endif()
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(
+    HIP
+    REQUIRED_VARS
+    HIP_ROOT_DIR
+    HIP_HIPCC_EXECUTABLE
+    HIP_HIPCONFIG_EXECUTABLE
+    HIP_PLATFORM
+    HIP_COMPILER
+    HIP_RUNTIME
+    VERSION_VAR HIP_VERSION
+    )
+
+###############################################################################
+# Set HIP CMAKE Flags
+###############################################################################
+# Copy the invocation styles from CXX to HIP
+set(CMAKE_HIP_ARCHIVE_CREATE ${CMAKE_CXX_ARCHIVE_CREATE})
+set(CMAKE_HIP_ARCHIVE_APPEND ${CMAKE_CXX_ARCHIVE_APPEND})
+set(CMAKE_HIP_ARCHIVE_FINISH ${CMAKE_CXX_ARCHIVE_FINISH})
+set(CMAKE_SHARED_LIBRARY_SONAME_HIP_FLAG ${CMAKE_SHARED_LIBRARY_SONAME_CXX_FLAG})
+set(CMAKE_SHARED_LIBRARY_CREATE_HIP_FLAGS ${CMAKE_SHARED_LIBRARY_CREATE_CXX_FLAGS})
+set(CMAKE_SHARED_LIBRARY_HIP_FLAGS ${CMAKE_SHARED_LIBRARY_CXX_FLAGS})
+#set(CMAKE_SHARED_LIBRARY_LINK_HIP_FLAGS ${CMAKE_SHARED_LIBRARY_LINK_CXX_FLAGS})
+set(CMAKE_SHARED_LIBRARY_RUNTIME_HIP_FLAG ${CMAKE_SHARED_LIBRARY_RUNTIME_CXX_FLAG})
+set(CMAKE_SHARED_LIBRARY_RUNTIME_HIP_FLAG_SEP ${CMAKE_SHARED_LIBRARY_RUNTIME_CXX_FLAG_SEP})
+set(CMAKE_SHARED_LIBRARY_LINK_STATIC_HIP_FLAGS ${CMAKE_SHARED_LIBRARY_LINK_STATIC_CXX_FLAGS})
+set(CMAKE_SHARED_LIBRARY_LINK_DYNAMIC_HIP_FLAGS ${CMAKE_SHARED_LIBRARY_LINK_DYNAMIC_CXX_FLAGS})
+
+set(HIP_CLANG_PARALLEL_BUILD_COMPILE_OPTIONS "")
+set(HIP_CLANG_PARALLEL_BUILD_LINK_OPTIONS "")
+
+if("${HIP_COMPILER}" STREQUAL "nvcc")
+    # Set the CMake Flags to use the nvcc Compiler.
+    set(CMAKE_HIP_CREATE_SHARED_LIBRARY "${HIP_HIPCC_CMAKE_LINKER_HELPER} <CMAKE_SHARED_LIBRARY_CXX_FLAGS> <LANGUAGE_COMPILE_FLAGS> <LINK_FLAGS> <CMAKE_SHARED_LIBRARY_CREATE_CXX_FLAGS> <SONAME_FLAG><TARGET_SONAME> -o <TARGET> <OBJECTS> <LINK_LIBRARIES>")
+    set(CMAKE_HIP_CREATE_SHARED_MODULE "${HIP_HIPCC_CMAKE_LINKER_HELPER} <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> <SONAME_FLAG><TARGET_SONAME> -o <TARGET> <LINK_LIBRARIES> -shared" )
+    set(CMAKE_HIP_LINK_EXECUTABLE "${HIP_HIPCC_CMAKE_LINKER_HELPER} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
+elseif("${HIP_COMPILER}" STREQUAL "hcc")
+    # Set the CMake Flags to use the hcc Compiler.
+    set(CMAKE_HIP_CREATE_SHARED_LIBRARY "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HCC_HOME} <CMAKE_SHARED_LIBRARY_CXX_FLAGS> <LANGUAGE_COMPILE_FLAGS> <LINK_FLAGS> <CMAKE_SHARED_LIBRARY_CREATE_CXX_FLAGS> <SONAME_FLAG><TARGET_SONAME> -o <TARGET> <OBJECTS> <LINK_LIBRARIES>")
+    set(CMAKE_HIP_CREATE_SHARED_MODULE "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HCC_HOME} <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> <SONAME_FLAG><TARGET_SONAME> -o <TARGET> <LINK_LIBRARIES> -shared" )
+    set(CMAKE_HIP_LINK_EXECUTABLE "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HCC_HOME} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
+elseif("${HIP_COMPILER}" STREQUAL "clang")
+    #Number of parallel jobs by default is 1
+    if(NOT DEFINED HIP_CLANG_NUM_PARALLEL_JOBS)
+      set(HIP_CLANG_NUM_PARALLEL_JOBS 1)
+    endif()
+    #Add support for parallel build and link
+    if(${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang")
+      check_cxx_compiler_flag("-parallel-jobs=1" HIP_CLANG_SUPPORTS_PARALLEL_JOBS)
+    endif()
+    if(HIP_CLANG_NUM_PARALLEL_JOBS GREATER 1)
+      if(${HIP_CLANG_SUPPORTS_PARALLEL_JOBS})
+        set(HIP_CLANG_PARALLEL_BUILD_COMPILE_OPTIONS "-Wno-format-nonliteral -parallel-jobs=${HIP_CLANG_NUM_PARALLEL_JOBS}")
+        set(HIP_CLANG_PARALLEL_BUILD_LINK_OPTIONS "-parallel-jobs=${HIP_CLANG_NUM_PARALLEL_JOBS}")
+      else()
+        message("clang compiler doesn't support parallel jobs")
+      endif()
+    endif()
+
+    # Set the CMake Flags to use the HIP-Clang Compiler.
+    set(CMAKE_HIP_CREATE_SHARED_LIBRARY "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HIP_CLANG_PATH} ${HIP_CLANG_PARALLEL_BUILD_LINK_OPTIONS} <CMAKE_SHARED_LIBRARY_CXX_FLAGS> <LANGUAGE_COMPILE_FLAGS> <LINK_FLAGS> <CMAKE_SHARED_LIBRARY_CREATE_CXX_FLAGS> <SONAME_FLAG><TARGET_SONAME> -o <TARGET> <OBJECTS> <LINK_LIBRARIES>")
+    set(CMAKE_HIP_CREATE_SHARED_MODULE "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HIP_CLANG_PATH} ${HIP_CLANG_PARALLEL_BUILD_LINK_OPTIONS} <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> <SONAME_FLAG><TARGET_SONAME> -o <TARGET> <LINK_LIBRARIES> -shared" )
+    set(CMAKE_HIP_LINK_EXECUTABLE "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HIP_CLANG_PATH} ${HIP_CLANG_PARALLEL_BUILD_LINK_OPTIONS} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
+
+    if("${HIP_RUNTIME}" STREQUAL "rocclr")
+      if(TARGET host)
+        message(STATUS "host interface - found")
+        set(HIP_HOST_INTERFACE host)
+      endif()
+    endif()
+endif()
+
+###############################################################################
+# MACRO: Locate helper files
+###############################################################################
+macro(HIP_FIND_HELPER_FILE _name _extension)
+    set(_hip_full_name "${_name}.${_extension}")
+    get_filename_component(CMAKE_CURRENT_LIST_DIR "${CMAKE_CURRENT_LIST_FILE}" PATH)
+    set(HIP_${_name} "${CMAKE_CURRENT_LIST_DIR}/FindHIP/${_hip_full_name}")
+    if(NOT EXISTS "${HIP_${_name}}")
+        set(error_message "${_hip_full_name} not found in ${CMAKE_CURRENT_LIST_DIR}/FindHIP")
+        if(HIP_FIND_REQUIRED)
+            message(FATAL_ERROR "${error_message}")
+        else()
+            if(NOT HIP_FIND_QUIETLY)
+                message(STATUS "${error_message}")
+            endif()
+        endif()
+    endif()
+    # Set this variable as internal, so the user isn't bugged with it.
+    set(HIP_${_name} ${HIP_${_name}} CACHE INTERNAL "Location of ${_full_name}" FORCE)
+endmacro()
+
+###############################################################################
+hip_find_helper_file(run_make2cmake cmake)
+hip_find_helper_file(run_hipcc cmake)
+###############################################################################
+
+###############################################################################
+# MACRO: Reset compiler flags
+###############################################################################
+macro(HIP_RESET_FLAGS)
+    unset(HIP_HIPCC_FLAGS)
+    unset(HIP_HCC_FLAGS)
+    unset(HIP_CLANG_FLAGS)
+    unset(HIP_NVCC_FLAGS)
+    foreach(config ${_hip_configuration_types})
+        string(TOUPPER ${config} config_upper)
+        unset(HIP_HIPCC_FLAGS_${config_upper})
+        unset(HIP_HCC_FLAGS_${config_upper})
+        unset(HIP_CLANG_FLAGS_${config_upper})
+        unset(HIP_NVCC_FLAGS_${config_upper})
+    endforeach()
+endmacro()
+
+###############################################################################
+# MACRO: Separate the options from the sources
+###############################################################################
+macro(HIP_GET_SOURCES_AND_OPTIONS _sources _cmake_options _hipcc_options _hcc_options _clang_options _nvcc_options)
+    set(${_sources})
+    set(${_cmake_options})
+    set(${_hipcc_options})
+    set(${_hcc_options})
+    set(${_clang_options})
+    set(${_nvcc_options})
+    set(_hipcc_found_options FALSE)
+    set(_hcc_found_options FALSE)
+    set(_clang_found_options FALSE)
+    set(_nvcc_found_options FALSE)
+    foreach(arg ${ARGN})
+        if("x${arg}" STREQUAL "xHIPCC_OPTIONS")
+            set(_hipcc_found_options TRUE)
+            set(_hcc_found_options FALSE)
+            set(_clang_found_options FALSE)
+            set(_nvcc_found_options FALSE)
+        elseif("x${arg}" STREQUAL "xHCC_OPTIONS")
+            set(_hipcc_found_options FALSE)
+            set(_hcc_found_options TRUE)
+            set(_clang_found_options FALSE)
+            set(_nvcc_found_options FALSE)
+        elseif("x${arg}" STREQUAL "xCLANG_OPTIONS")
+            set(_hipcc_found_options FALSE)
+            set(_hcc_found_options FALSE)
+            set(_clang_found_options TRUE)
+            set(_nvcc_found_options FALSE)
+        elseif("x${arg}" STREQUAL "xNVCC_OPTIONS")
+            set(_hipcc_found_options FALSE)
+            set(_hcc_found_options FALSE)
+            set(_clang_found_options FALSE)
+            set(_nvcc_found_options TRUE)
+        elseif(
+                "x${arg}" STREQUAL "xEXCLUDE_FROM_ALL" OR
+                "x${arg}" STREQUAL "xSTATIC" OR
+                "x${arg}" STREQUAL "xSHARED" OR
+                "x${arg}" STREQUAL "xMODULE"
+                )
+            list(APPEND ${_cmake_options} ${arg})
+        else()
+            if(_hipcc_found_options)
+                list(APPEND ${_hipcc_options} ${arg})
+            elseif(_hcc_found_options)
+                list(APPEND ${_hcc_options} ${arg})
+            elseif(_clang_found_options)
+                list(APPEND ${_clang_options} ${arg})
+            elseif(_nvcc_found_options)
+                list(APPEND ${_nvcc_options} ${arg})
+            else()
+                # Assume this is a file
+                list(APPEND ${_sources} ${arg})
+            endif()
+        endif()
+    endforeach()
+endmacro()
+
+###############################################################################
+# MACRO: Add include directories to pass to the hipcc command
+###############################################################################
+set(HIP_HIPCC_INCLUDE_ARGS_USER "")
+macro(HIP_INCLUDE_DIRECTORIES)
+    foreach(dir ${ARGN})
+        list(APPEND HIP_HIPCC_INCLUDE_ARGS_USER $<$<BOOL:${dir}>:-I${dir}>)
+    endforeach()
+endmacro()
+
+###############################################################################
+# FUNCTION: Helper to avoid clashes of files with the same basename but different paths
+###############################################################################
+function(HIP_COMPUTE_BUILD_PATH path build_path)
+    # Convert to cmake style paths
+    file(TO_CMAKE_PATH "${path}" bpath)
+    if(IS_ABSOLUTE "${bpath}")
+        string(FIND "${bpath}" "${CMAKE_CURRENT_BINARY_DIR}" _binary_dir_pos)
+        if(_binary_dir_pos EQUAL 0)
+            file(RELATIVE_PATH bpath "${CMAKE_CURRENT_BINARY_DIR}" "${bpath}")
+        else()
+            file(RELATIVE_PATH bpath "${CMAKE_CURRENT_SOURCE_DIR}" "${bpath}")
+        endif()
+    endif()
+
+    # Remove leading /
+    string(REGEX REPLACE "^[/]+" "" bpath "${bpath}")
+    # Avoid absolute paths by removing ':'
+    string(REPLACE ":" "_" bpath "${bpath}")
+    # Avoid relative paths that go up the tree
+    string(REPLACE "../" "__/" bpath "${bpath}")
+    # Avoid spaces
+    string(REPLACE " " "_" bpath "${bpath}")
+    # Strip off the filename
+    get_filename_component(bpath "${bpath}" PATH)
+
+    set(${build_path} "${bpath}" PARENT_SCOPE)
+endfunction()
+
+###############################################################################
+# MACRO: Parse OPTIONS from ARGN & set variables prefixed by _option_prefix
+###############################################################################
+macro(HIP_PARSE_HIPCC_OPTIONS _option_prefix)
+    set(_hip_found_config)
+    foreach(arg ${ARGN})
+        # Determine if we are dealing with a per-configuration flag
+        foreach(config ${_hip_configuration_types})
+            string(TOUPPER ${config} config_upper)
+            if(arg STREQUAL "${config_upper}")
+                set(_hip_found_config _${arg})
+                # Clear arg to prevent it from being processed anymore
+                set(arg)
+            endif()
+        endforeach()
+        if(arg)
+            list(APPEND ${_option_prefix}${_hip_found_config} "${arg}")
+        endif()
+    endforeach()
+endmacro()
+
+###############################################################################
+# MACRO: Try and include dependency file if it exists
+###############################################################################
+macro(HIP_INCLUDE_HIPCC_DEPENDENCIES dependency_file)
+    set(HIP_HIPCC_DEPEND)
+    set(HIP_HIPCC_DEPEND_REGENERATE FALSE)
+
+    # Create the dependency file if it doesn't exist
+    if(NOT EXISTS ${dependency_file})
+        file(WRITE ${dependency_file} "# Generated by: FindHIP.cmake. Do not edit.\n")
+    endif()
+    # Include the dependency file
+    include(${dependency_file})
+
+    # Verify the existence of all the included files
+    if(HIP_HIPCC_DEPEND)
+        foreach(f ${HIP_HIPCC_DEPEND})
+            if(NOT EXISTS ${f})
+                # If they aren't there, regenerate the file again
+                set(HIP_HIPCC_DEPEND_REGENERATE TRUE)
+            endif()
+        endforeach()
+    else()
+        # No dependencies, so regenerate the file
+        set(HIP_HIPCC_DEPEND_REGENERATE TRUE)
+    endif()
+
+    # Regenerate the dependency file if needed
+    if(HIP_HIPCC_DEPEND_REGENERATE)
+        set(HIP_HIPCC_DEPEND ${dependency_file})
+        file(WRITE ${dependency_file} "# Generated by: FindHIP.cmake. Do not edit.\n")
+    endif()
+endmacro()
+
+###############################################################################
+# MACRO: Prepare cmake commands for the target
+###############################################################################
+macro(HIP_PREPARE_TARGET_COMMANDS _target _format _generated_files _source_files)
+    set(_hip_flags "")
+    string(TOUPPER "${CMAKE_BUILD_TYPE}" _hip_build_configuration)
+    if(HIP_HOST_COMPILATION_CPP)
+        set(HIP_C_OR_CXX CXX)
+    else()
+        set(HIP_C_OR_CXX C)
+    endif()
+    set(generated_extension ${CMAKE_${HIP_C_OR_CXX}_OUTPUT_EXTENSION})
+
+    # Initialize list of includes with those specified by the user. Append with
+    # ones specified to cmake directly.
+    set(HIP_HIPCC_INCLUDE_ARGS ${HIP_HIPCC_INCLUDE_ARGS_USER})
+
+    # Add the include directories
+    set(include_directories_generator "$<TARGET_PROPERTY:${_target},INCLUDE_DIRECTORIES>")
+    list(APPEND HIP_HIPCC_INCLUDE_ARGS "$<$<BOOL:${include_directories_generator}>:-I$<JOIN:${include_directories_generator}, -I>>")
+
+    get_directory_property(_hip_include_directories INCLUDE_DIRECTORIES)
+    list(REMOVE_DUPLICATES _hip_include_directories)
+    if(_hip_include_directories)
+        foreach(dir ${_hip_include_directories})
+            list(APPEND HIP_HIPCC_INCLUDE_ARGS $<$<BOOL:${dir}>:-I${dir}>)
+        endforeach()
+    endif()
+
+    HIP_GET_SOURCES_AND_OPTIONS(_hip_sources _hip_cmake_options _hipcc_options _hcc_options _clang_options _nvcc_options ${ARGN})
+    HIP_PARSE_HIPCC_OPTIONS(HIP_HIPCC_FLAGS ${_hipcc_options})
+    HIP_PARSE_HIPCC_OPTIONS(HIP_HCC_FLAGS ${_hcc_options})
+    HIP_PARSE_HIPCC_OPTIONS(HIP_CLANG_FLAGS ${_clang_options})
+    HIP_PARSE_HIPCC_OPTIONS(HIP_NVCC_FLAGS ${_nvcc_options})
+
+    # Add the compile definitions
+    set(compile_definition_generator "$<TARGET_PROPERTY:${_target},COMPILE_DEFINITIONS>")
+    list(APPEND HIP_HIPCC_FLAGS "$<$<BOOL:${compile_definition_generator}>:-D$<JOIN:${compile_definition_generator}, -D>>")
+
+    # Check if we are building shared library.
+    set(_hip_build_shared_libs FALSE)
+    list(FIND _hip_cmake_options SHARED _hip_found_SHARED)
+    list(FIND _hip_cmake_options MODULE _hip_found_MODULE)
+    if(_hip_found_SHARED GREATER -1 OR _hip_found_MODULE GREATER -1)
+        set(_hip_build_shared_libs TRUE)
+    endif()
+    list(FIND _hip_cmake_options STATIC _hip_found_STATIC)
+    if(_hip_found_STATIC GREATER -1)
+        set(_hip_build_shared_libs FALSE)
+    endif()
+
+    # If we are building a shared library, add extra flags to HIP_HIPCC_FLAGS
+    if(_hip_build_shared_libs)
+        list(APPEND HIP_HCC_FLAGS "-fPIC")
+        list(APPEND HIP_CLANG_FLAGS "-fPIC")
+        list(APPEND HIP_NVCC_FLAGS "--shared -Xcompiler '-fPIC'")
+    endif()
+
+    # Set host compiler
+    set(HIP_HOST_COMPILER "${CMAKE_${HIP_C_OR_CXX}_COMPILER}")
+
+    # Set compiler flags
+    set(_HIP_HOST_FLAGS "set(CMAKE_HOST_FLAGS ${CMAKE_${HIP_C_OR_CXX}_FLAGS})")
+    set(_HIP_HIPCC_FLAGS "set(HIP_HIPCC_FLAGS ${HIP_HIPCC_FLAGS})")
+    set(_HIP_HCC_FLAGS "set(HIP_HCC_FLAGS ${HIP_HCC_FLAGS})")
+    set(_HIP_CLANG_FLAGS "set(HIP_CLANG_FLAGS ${HIP_CLANG_FLAGS})")
+    set(_HIP_NVCC_FLAGS "set(HIP_NVCC_FLAGS ${HIP_NVCC_FLAGS})")
+    foreach(config ${_hip_configuration_types})
+        string(TOUPPER ${config} config_upper)
+        set(_HIP_HOST_FLAGS "${_HIP_HOST_FLAGS}\nset(CMAKE_HOST_FLAGS_${config_upper} ${CMAKE_${HIP_C_OR_CXX}_FLAGS_${config_upper}})")
+        set(_HIP_HIPCC_FLAGS "${_HIP_HIPCC_FLAGS}\nset(HIP_HIPCC_FLAGS_${config_upper} ${HIP_HIPCC_FLAGS_${config_upper}})")
+        set(_HIP_HCC_FLAGS "${_HIP_HCC_FLAGS}\nset(HIP_HCC_FLAGS_${config_upper} ${HIP_HCC_FLAGS_${config_upper}})")
+        set(_HIP_CLANG_FLAGS "${_HIP_CLANG_FLAGS}\nset(HIP_CLANG_FLAGS_${config_upper} ${HIP_CLANG_FLAGS_${config_upper}})")
+        set(_HIP_NVCC_FLAGS "${_HIP_NVCC_FLAGS}\nset(HIP_NVCC_FLAGS_${config_upper} ${HIP_NVCC_FLAGS_${config_upper}})")
+    endforeach()
+
+    # Reset the output variable
+    set(_hip_generated_files "")
+    set(_hip_source_files "")
+
+    # Iterate over all arguments and create custom commands for all source files
+    foreach(file ${ARGN})
+        # Ignore any file marked as a HEADER_FILE_ONLY
+        get_source_file_property(_is_header ${file} HEADER_FILE_ONLY)
+        # Allow per source file overrides of the format. Also allows compiling non .cu files.
+        get_source_file_property(_hip_source_format ${file} HIP_SOURCE_PROPERTY_FORMAT)
+        if((${file} MATCHES "\\.cu$" OR _hip_source_format) AND NOT _is_header)
+            set(host_flag FALSE)
+        else()
+            set(host_flag TRUE)
+        endif()
+
+        if(NOT host_flag)
+            # Determine output directory
+            HIP_COMPUTE_BUILD_PATH("${file}" hip_build_path)
+            set(hip_compile_output_dir "${CMAKE_CURRENT_BINARY_DIR}/CMakeFiles/${_target}.dir/${hip_build_path}")
+
+            get_filename_component(basename ${file} NAME)
+            set(generated_file_path "${hip_compile_output_dir}/${CMAKE_CFG_INTDIR}")
+            set(generated_file_basename "${_target}_generated_${basename}${generated_extension}")
+
+            # Set file names
+            set(generated_file "${generated_file_path}/${generated_file_basename}")
+            set(cmake_dependency_file "${hip_compile_output_dir}/${generated_file_basename}.depend")
+            set(custom_target_script_pregen "${hip_compile_output_dir}/${generated_file_basename}.cmake.pre-gen")
+            set(custom_target_script "${hip_compile_output_dir}/${generated_file_basename}.cmake")
+
+            # Set properties for object files
+            set_source_files_properties("${generated_file}"
+                PROPERTIES
+                EXTERNAL_OBJECT true # This is an object file not to be compiled, but only be linked
+                )
+
+            # Don't add CMAKE_CURRENT_SOURCE_DIR if the path is already an absolute path
+            get_filename_component(file_path "${file}" PATH)
+            if(IS_ABSOLUTE "${file_path}")
+                set(source_file "${file}")
+            else()
+                set(source_file "${CMAKE_CURRENT_SOURCE_DIR}/${file}")
+            endif()
+
+            # Bring in the dependencies
+            HIP_INCLUDE_HIPCC_DEPENDENCIES(${cmake_dependency_file})
+
+            # Configure the build script
+            configure_file("${HIP_run_hipcc}" "${custom_target_script_pregen}" @ONLY)
+            file(GENERATE
+                OUTPUT "${custom_target_script}"
+                INPUT "${custom_target_script_pregen}"
+                )
+            set(main_dep DEPENDS ${source_file})
+            if(CMAKE_GENERATOR MATCHES "Makefiles")
+                set(verbose_output "$(VERBOSE)")
+            elseif(HIP_VERBOSE_BUILD)
+                set(verbose_output ON)
+            else()
+                set(verbose_output OFF)
+            endif()
+
+            # Create up the comment string
+            file(RELATIVE_PATH generated_file_relative_path "${CMAKE_BINARY_DIR}" "${generated_file}")
+            set(hip_build_comment_string "Building HIPCC object ${generated_file_relative_path}")
+
+            # Build the generated file and dependency file
+            add_custom_command(
+                OUTPUT ${generated_file}
+                # These output files depend on the source_file and the contents of cmake_dependency_file
+                ${main_dep}
+                DEPENDS ${HIP_HIPCC_DEPEND}
+                DEPENDS ${custom_target_script}
+                # Make sure the output directory exists before trying to write to it.
+                COMMAND ${CMAKE_COMMAND} -E make_directory "${generated_file_path}"
+                COMMAND ${CMAKE_COMMAND} ARGS
+                -D verbose:BOOL=${verbose_output}
+                -D build_configuration:STRING=${_hip_build_configuration}
+                -D "generated_file:STRING=${generated_file}"
+                -P "${custom_target_script}"
+                WORKING_DIRECTORY "${hip_compile_output_dir}"
+                COMMENT "${hip_build_comment_string}"
+                )
+
+            # Make sure the build system knows the file is generated
+            set_source_files_properties(${generated_file} PROPERTIES GENERATED TRUE)
+            list(APPEND _hip_generated_files ${generated_file})
+            list(APPEND _hip_source_files ${file})
+        endif()
+    endforeach()
+
+    # Set the return parameter
+    set(${_generated_files} ${_hip_generated_files})
+    set(${_source_files} ${_hip_source_files})
+endmacro()
+
+###############################################################################
+# HIP_ADD_EXECUTABLE
+###############################################################################
+macro(HIP_ADD_EXECUTABLE hip_target)
+    # Separate the sources from the options
+    HIP_GET_SOURCES_AND_OPTIONS(_sources _cmake_options _hipcc_options _hcc_options _clang_options _nvcc_options ${ARGN})
+    HIP_PREPARE_TARGET_COMMANDS(${hip_target} OBJ _generated_files _source_files ${_sources} HIPCC_OPTIONS ${_hipcc_options} HCC_OPTIONS ${_hcc_options} CLANG_OPTIONS ${_clang_options} NVCC_OPTIONS ${_nvcc_options})
+    if(_source_files)
+        list(REMOVE_ITEM _sources ${_source_files})
+    endif()
+    if("${HIP_COMPILER}" STREQUAL "hcc")
+        if("x${HCC_HOME}" STREQUAL "x")
+            if (DEFINED ENV{ROCM_PATH})
+                set(HCC_HOME "$ENV{ROCM_PATH}/hcc")
+            elseif(DEFINED ENV{HIP_PATH})
+                set(HCC_HOME "$ENV{HIP_PATH}/../hcc")
+            else()
+                set(HCC_HOME "/opt/rocm/hcc")
+            endif()
+        endif()
+        set(CMAKE_HIP_LINK_EXECUTABLE "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HCC_HOME} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
+    elseif("${HIP_COMPILER}" STREQUAL "clang")
+        if("x${HIP_CLANG_PATH}" STREQUAL "x")
+            if(DEFINED ENV{HIP_CLANG_PATH})
+                set(HIP_CLANG_PATH $ENV{HIP_CLANG_PATH})
+            elseif(DEFINED ENV{ROCM_PATH})
+                set(HIP_CLANG_PATH "$ENV{ROCM_PATH}/llvm/bin")
+            elseif(DEFINED ENV{HIP_PATH})
+                set(HIP_CLANG_PATH "$ENV{HIP_PATH}/../llvm/bin")
+            else()
+                set(HIP_CLANG_PATH "/opt/rocm/llvm/bin")
+            endif()
+        endif()
+        set(CMAKE_HIP_LINK_EXECUTABLE "${HIP_HIPCC_CMAKE_LINKER_HELPER} ${HIP_CLANG_PATH} ${HIP_CLANG_PARALLEL_BUILD_LINK_OPTIONS} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
+    else()
+        set(CMAKE_HIP_LINK_EXECUTABLE "${HIP_HIPCC_CMAKE_LINKER_HELPER} <FLAGS> <CMAKE_CXX_LINK_FLAGS> <LINK_FLAGS> <OBJECTS> -o <TARGET> <LINK_LIBRARIES>")
+    endif()
+    if ("${_sources}" STREQUAL "")
+        add_executable(${hip_target} ${_cmake_options} ${_generated_files} "")
+    else()
+        add_executable(${hip_target} ${_cmake_options} ${_generated_files} ${_sources})
+    endif()
+    set_target_properties(${hip_target} PROPERTIES LINKER_LANGUAGE HIP)
+    # Link with host
+    if (HIP_HOST_INTERFACE)
+        # hip rt should be rocclr, compiler should be clang
+        target_link_libraries(${hip_target} ${HIP_HOST_INTERFACE})
+    endif()
+endmacro()
+
+###############################################################################
+# HIP_ADD_LIBRARY
+###############################################################################
+macro(HIP_ADD_LIBRARY hip_target)
+    # Separate the sources from the options
+    HIP_GET_SOURCES_AND_OPTIONS(_sources _cmake_options _hipcc_options _hcc_options _clang_options _nvcc_options ${ARGN})
+    HIP_PREPARE_TARGET_COMMANDS(${hip_target} OBJ _generated_files _source_files ${_sources} ${_cmake_options} HIPCC_OPTIONS ${_hipcc_options} HCC_OPTIONS ${_hcc_options} CLANG_OPTIONS ${_clang_options} NVCC_OPTIONS ${_nvcc_options})
+    if(_source_files)
+        list(REMOVE_ITEM _sources ${_source_files})
+    endif()
+    if ("${_sources}" STREQUAL "")
+        add_library(${hip_target} ${_cmake_options} ${_generated_files} "")
+    else()
+        add_library(${hip_target} ${_cmake_options} ${_generated_files} ${_sources})
+    endif()
+    set_target_properties(${hip_target} PROPERTIES LINKER_LANGUAGE ${HIP_C_OR_CXX})
+    # Link with host
+    if (HIP_HOST_INTERFACE)
+        # hip rt should be rocclr, compiler should be clang
+        target_link_libraries(${hip_target} ${HIP_HOST_INTERFACE})
+    endif()
+endmacro()
+
+# vim: ts=4:sw=4:expandtab:smartindent

config/cmake/modules/FindHIP/run_hipcc.cmake

@@ -0,0 +1,182 @@
+###############################################################################
+# Runs commands using HIPCC
+###############################################################################
+
+###############################################################################
+# This file runs the hipcc commands to produce the desired output file
+# along with the dependency file needed by CMake to compute dependencies.
+#
+# Input variables:
+#
+# verbose:BOOL=<>               OFF: Be as quiet as possible (default)
+#                               ON : Describe each step
+# build_configuration:STRING=<> Build configuration. Defaults to Debug.
+# generated_file:STRING=<>      File to generate. Mandatory argument.
+
+if(NOT build_configuration)
+    set(build_configuration Debug)
+endif()
+if(NOT generated_file)
+    message(FATAL_ERROR "You must specify generated_file on the command line")
+endif()
+
+# Set these up as variables to make reading the generated file easier
+set(HIP_HIPCC_EXECUTABLE "@HIP_HIPCC_EXECUTABLE@") # path
+set(HIP_HIPCONFIG_EXECUTABLE "@HIP_HIPCONFIG_EXECUTABLE@") #path
+set(HIP_HOST_COMPILER "@HIP_HOST_COMPILER@") # path
+set(CMAKE_COMMAND "@CMAKE_COMMAND@") # path
+set(HIP_run_make2cmake "@HIP_run_make2cmake@") # path
+set(HCC_HOME "@HCC_HOME@") #path
+set(HIP_CLANG_PATH "@HIP_CLANG_PATH@") #path
+set(HIP_CLANG_PARALLEL_BUILD_COMPILE_OPTIONS "@HIP_CLANG_PARALLEL_BUILD_COMPILE_OPTIONS@")
+
+@HIP_HOST_FLAGS@
+@_HIP_HIPCC_FLAGS@
+@_HIP_HCC_FLAGS@
+@_HIP_CLANG_FLAGS@
+@_HIP_NVCC_FLAGS@
+#Needed to bring the HIP_HIPCC_INCLUDE_ARGS variable in scope
+set(HIP_HIPCC_INCLUDE_ARGS @HIP_HIPCC_INCLUDE_ARGS@) # list
+
+set(cmake_dependency_file "@cmake_dependency_file@") # path
+set(source_file "@source_file@") # path
+set(host_flag "@host_flag@") # bool
+
+# Determine compiler and compiler flags
+execute_process(COMMAND ${HIP_HIPCONFIG_EXECUTABLE} --platform OUTPUT_VARIABLE HIP_PLATFORM OUTPUT_STRIP_TRAILING_WHITESPACE)
+execute_process(COMMAND ${HIP_HIPCONFIG_EXECUTABLE} --compiler OUTPUT_VARIABLE HIP_COMPILER OUTPUT_STRIP_TRAILING_WHITESPACE)
+execute_process(COMMAND ${HIP_HIPCONFIG_EXECUTABLE} --runtime OUTPUT_VARIABLE HIP_RUNTIME OUTPUT_STRIP_TRAILING_WHITESPACE)
+if(NOT host_flag)
+    set(__CC ${HIP_HIPCC_EXECUTABLE})
+    if("${HIP_PLATFORM}" STREQUAL "amd")
+        if("${HIP_COMPILER}" STREQUAL "hcc")
+            if(NOT "x${HCC_HOME}" STREQUAL "x")
+                set(ENV{HCC_HOME} ${HCC_HOME})
+            endif()
+            set(__CC_FLAGS ${HIP_HIPCC_FLAGS} ${HIP_HCC_FLAGS} ${HIP_HIPCC_FLAGS_${build_configuration}} ${HIP_HCC_FLAGS_${build_configuration}})
+        elseif("${HIP_COMPILER}" STREQUAL "clang")
+            if(NOT "x${HIP_CLANG_PATH}" STREQUAL "x")
+                set(ENV{HIP_CLANG_PATH} ${HIP_CLANG_PATH})
+            endif()
+            # Temporarily include HIP_HCC_FLAGS for HIP-Clang for PyTorch builds
+            set(__CC_FLAGS ${HIP_CLANG_PARALLEL_BUILD_COMPILE_OPTIONS} ${HIP_HIPCC_FLAGS} ${HIP_HCC_FLAGS} ${HIP_CLANG_FLAGS} ${HIP_HIPCC_FLAGS_${build_configuration}} ${HIP_HCC_FLAGS_${build_configuration}} ${HIP_CLANG_FLAGS_${build_configuration}})
+        endif()
+    else()
+        set(__CC_FLAGS ${HIP_HIPCC_FLAGS} ${HIP_NVCC_FLAGS} ${HIP_HIPCC_FLAGS_${build_configuration}} ${HIP_NVCC_FLAGS_${build_configuration}})
+    endif()
+else()
+    set(__CC ${HIP_HOST_COMPILER})
+    set(__CC_FLAGS ${CMAKE_HOST_FLAGS} ${CMAKE_HOST_FLAGS_${build_configuration}})
+endif()
+set(__CC_INCLUDES ${HIP_HIPCC_INCLUDE_ARGS})
+
+# hip_execute_process - Executes a command with optional command echo and status message.
+#   status     - Status message to print if verbose is true
+#   command    - COMMAND argument from the usual execute_process argument structure
+#   ARGN       - Remaining arguments are the command with arguments
+#   HIP_result - Return value from running the command
+macro(hip_execute_process status command)
+    set(_command ${command})
+    if(NOT "x${_command}" STREQUAL "xCOMMAND")
+        message(FATAL_ERROR "Malformed call to hip_execute_process.  Missing COMMAND as second argument. (command = ${command})")
+    endif()
+    if(verbose)
+        execute_process(COMMAND "${CMAKE_COMMAND}" -E echo -- ${status})
+        # Build command string to print
+        set(hip_execute_process_string)
+        foreach(arg ${ARGN})
+            # Escape quotes if any
+            string(REPLACE "\"" "\\\"" arg ${arg})
+            # Surround args with spaces with quotes
+            if(arg MATCHES " ")
+                list(APPEND hip_execute_process_string "\"${arg}\"")
+            else()
+                list(APPEND hip_execute_process_string ${arg})
+            endif()
+        endforeach()
+        # Echo the command
+        execute_process(COMMAND ${CMAKE_COMMAND} -E echo ${hip_execute_process_string})
+    endif()
+    # Run the command
+    execute_process(COMMAND ${ARGN} RESULT_VARIABLE HIP_result)
+endmacro()
+
+# Delete the target file
+hip_execute_process(
+    "Removing ${generated_file}"
+    COMMAND "${CMAKE_COMMAND}" -E remove "${generated_file}"
+    )
+
+# Generate the dependency file
+hip_execute_process(
+    "Generating dependency file: ${cmake_dependency_file}.pre"
+    COMMAND "${__CC}"
+    -M
+    "${source_file}"
+    -o "${cmake_dependency_file}.pre"
+    ${__CC_FLAGS}
+    ${__CC_INCLUDES}
+    )
+
+if(HIP_result)
+    message(FATAL_ERROR "Error generating ${generated_file}")
+endif()
+
+# Generate the cmake readable dependency file to a temp file
+hip_execute_process(
+    "Generating temporary cmake readable file: ${cmake_dependency_file}.tmp"
+    COMMAND "${CMAKE_COMMAND}"
+    -D "input_file:FILEPATH=${cmake_dependency_file}.pre"
+    -D "output_file:FILEPATH=${cmake_dependency_file}.tmp"
+    -D "verbose=${verbose}"
+    -P "${HIP_run_make2cmake}"
+    )
+
+if(HIP_result)
+    message(FATAL_ERROR "Error generating ${generated_file}")
+endif()
+
+# Copy the file if it is different
+hip_execute_process(
+    "Copy if different ${cmake_dependency_file}.tmp to ${cmake_dependency_file}"
+    COMMAND "${CMAKE_COMMAND}" -E copy_if_different "${cmake_dependency_file}.tmp" "${cmake_dependency_file}"
+    )
+
+if(HIP_result)
+    message(FATAL_ERROR "Error generating ${generated_file}")
+endif()
+
+# Delete the temporary file
+hip_execute_process(
+    "Removing ${cmake_dependency_file}.tmp and ${cmake_dependency_file}.pre"
+    COMMAND "${CMAKE_COMMAND}" -E remove "${cmake_dependency_file}.tmp" "${cmake_dependency_file}.pre"
+    )
+
+if(HIP_result)
+    message(FATAL_ERROR "Error generating ${generated_file}")
+endif()
+
+# Generate the output file
+hip_execute_process(
+    "Generating ${generated_file}"
+    COMMAND "${__CC}"
+    -c
+    "${source_file}"
+    -o "${generated_file}"
+    ${__CC_FLAGS}
+    ${__CC_INCLUDES}
+    )
+
+if(HIP_result)
+    # Make sure that we delete the output file
+    hip_execute_process(
+        "Removing ${generated_file}"
+        COMMAND "${CMAKE_COMMAND}" -E remove "${generated_file}"
+        )
+    message(FATAL_ERROR "Error generating file ${generated_file}")
+else()
+    if(verbose)
+        message("Generated ${generated_file} successfully.")
+    endif()
+endif()
+# vim: ts=4:sw=4:expandtab:smartindent

config/cmake/modules/FindHIP/run_make2cmake.cmake

@@ -0,0 +1,50 @@
+###############################################################################
+# Computes dependencies using HIPCC
+###############################################################################
+
+###############################################################################
+# This file converts dependency files generated using hipcc to a format that
+# cmake can understand.
+
+# Input variables:
+#
+# input_file:STRING=<> Dependency file to parse. Required argument
+# output_file:STRING=<> Output file to generate. Required argument
+
+if(NOT input_file OR NOT output_file)
+    message(FATAL_ERROR "You must specify input_file and output_file on the command line")
+endif()
+
+file(READ ${input_file} depend_text)
+
+if (NOT "${depend_text}" STREQUAL "")
+    string(REPLACE " /" "\n/" depend_text ${depend_text})
+    string(REGEX REPLACE "^.*:" "" depend_text ${depend_text})
+    string(REGEX REPLACE "[ \\\\]*\n" ";" depend_text ${depend_text})
+
+    set(dependency_list "")
+
+    foreach(file ${depend_text})
+        string(REGEX REPLACE "^ +" "" file ${file})
+        if(NOT EXISTS "${file}")
+            message(WARNING " Removing non-existent dependency file: ${file}")
+            set(file "")
+        endif()
+
+        if(NOT IS_DIRECTORY "${file}")
+            get_filename_component(file_absolute "${file}" ABSOLUTE)
+            list(APPEND dependency_list "${file_absolute}")
+        endif()
+    endforeach()
+endif()
+
+# Remove the duplicate entries and sort them.
+list(REMOVE_DUPLICATES dependency_list)
+list(SORT dependency_list)
+
+foreach(file ${dependency_list})
+    set(hip_hipcc_depend "${hip_hipcc_depend} \"${file}\"\n")
+endforeach()
+
+file(WRITE ${output_file} "# Generated by: FindHIP.cmake. Do not edit.\nSET(HIP_HIPCC_DEPEND\n ${hip_hipcc_depend})\n\n")
+# vim: ts=4:sw=4:expandtab:smartindent

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -43,6 +43,24 @@ function(convert_filenames_to_full_paths NAMES)
   set(${NAMES} ${tmp_names} PARENT_SCOPE)
 endfunction()
 
+# Wrapper for add_executable that calls the HIP wrapper if applicable
+macro(mfem_add_executable NAME)
+  if (MFEM_USE_HIP)
+    hip_add_executable(${NAME} ${ARGN})
+  else()
+    add_executable(${NAME} ${ARGN})
+  endif()
+endmacro()
+
+# Wrapper for add_library that calls the HIP wrapper if applicable
+macro(mfem_add_library NAME)
+  if (MFEM_USE_HIP)
+    hip_add_library(${NAME} ${ARGN})
+  else()
+    add_library(${NAME} ${ARGN})
+  endif()
+endmacro()
+
 # Simple shortcut to add_custom_target() with option to add the target to the
 # main target.
 function(add_mfem_target TARGET_NAME ADD_TO_ALL)
@@ -56,7 +74,7 @@ function(add_mfem_target TARGET_NAME ADD_TO_ALL)
 endfunction()
 
 # Add mfem examples
-function(add_mfem_examples EXE_SRCS)
+macro(add_mfem_examples EXE_SRCS)
   set(EXE_PREFIX "")
   set(EXE_PREREQUISITE "")
   set(EXE_NEEDED_BY "")
@@ -72,13 +90,15 @@ function(add_mfem_examples EXE_SRCS)
   foreach(SRC_FILE IN LISTS ${EXE_SRCS})
     # If CUDA is enabled, tag source files to be compiled with nvcc.
     if (MFEM_USE_CUDA)
-      set_property(SOURCE ${SRC_FILE} PROPERTY LANGUAGE CUDA)
+      set_source_files_properties(${SRC_FILE} PROPERTIES LANGUAGE CUDA)
+    elseif(MFEM_USE_HIP)
+      set_source_files_properties(${SRC_FILE} PROPERTIES HIP_SOURCE_PROPERTY_FORMAT TRUE)
     endif()
 
     get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
 
     string(REPLACE ".cpp" "" EXE_NAME "${EXE_PREFIX}${SRC_FILENAME}")
-    add_executable(${EXE_NAME} ${SRC_FILE})
+    mfem_add_executable(${EXE_NAME} ${SRC_FILE})
     add_dependencies(${MFEM_ALL_EXAMPLES_TARGET_NAME} ${EXE_NAME})
     if (EXE_NEEDED_BY)
       add_dependencies(${EXE_NEEDED_BY} ${EXE_NAME})
@@ -107,10 +127,10 @@ function(add_mfem_examples EXE_SRCS)
       endif()
     endif()
   endforeach(SRC_FILE)
-endfunction()
+endmacro()
 
 # A slightly more versatile function for adding miniapps to MFEM
-function(add_mfem_miniapp MFEM_EXE_NAME)
+macro(add_mfem_miniapp MFEM_EXE_NAME)
   # Parse the input arguments looking for the things we need
   set(POSSIBLE_ARGS "MAIN" "EXTRA_SOURCES" "EXTRA_HEADERS" "EXTRA_OPTIONS" "EXTRA_DEFINES" "LIBRARIES")
   set(CURRENT_ARG)
@@ -126,8 +146,7 @@ function(add_mfem_miniapp MFEM_EXE_NAME)
 
   # If CUDA is enabled, tag source files to be compiled with nvcc.
   if (MFEM_USE_CUDA)
-    set_property(SOURCE ${MAIN_LIST} ${EXTRA_SOURCES_LIST}
-      PROPERTY LANGUAGE CUDA)
+    set_source_files_properties(${MAIN_LIST} ${EXTRA_SOURCES_LIST} PROPERTIES LANGUAGE CUDA)
     if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.12.0)
       list(TRANSFORM EXTRA_OPTIONS_LIST PREPEND "-Xcompiler=")
     else()
@@ -137,11 +156,13 @@ function(add_mfem_miniapp MFEM_EXE_NAME)
       endforeach()
       set(EXTRA_OPTIONS_LIST ${LIST_})
     endif()
+  elseif(MFEM_USE_HIP)
+    set_source_files_properties(${MAIN_LIST} ${EXTRA_SOURCES_LIST} PROPERTIES HIP_SOURCE_PROPERTY_FORMAT TRUE)
   endif()
 
   # Actually add the executable
-  add_executable(${MFEM_EXE_NAME} ${MAIN_LIST}
-    ${EXTRA_SOURCES_LIST} ${EXTRA_HEADERS_LIST})
+  mfem_add_executable(${MFEM_EXE_NAME} ${MAIN_LIST}
+      ${EXTRA_SOURCES_LIST} ${EXTRA_HEADERS_LIST})
   add_dependencies(${MFEM_ALL_MINIAPPS_TARGET_NAME} ${MFEM_EXE_NAME})
   add_dependencies(${MFEM_EXE_NAME} ${MFEM_EXEC_PREREQUISITES_TARGET_NAME})
 
@@ -194,7 +215,7 @@ function(add_mfem_miniapp MFEM_EXE_NAME)
         LINK_FLAGS "${MPI_CXX_LINK_FLAGS}")
     endif()
   endif()
-endfunction()
+endmacro()
 
 
 # Auxiliary function, used in mfem_find_package().

config/defaults.cmake

@@ -76,6 +76,9 @@ option(MFEM_ENABLE_MINIAPPS "Build all of the miniapps" OFF)
 # Set the target CUDA architecture
 set(CUDA_ARCH "sm_60" CACHE STRING "Target CUDA architecture.")
 
+# Set the target HIP architecture
+set(HIP_ARCH "gfx900" CACHE STRING "Target HIP architecture.")
+
 set(MFEM_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 
 # The *_DIR paths below will be the first place searched for the corresponding

examples/CMakeLists.txt

@@ -34,7 +34,6 @@ list(APPEND ALL_EXE_SRCS
   ex25.cpp
   ex26.cpp
   ex27.cpp
-  ex91.cpp
   )
 
 if (MFEM_USE_MPI)

examples/drl_lshape.cpp

@@ -0,0 +1,340 @@
+// Compile with: make drl_shock_wave
+//
+// drl_shock_wave -o 2 -m ../data/inline-quad.mesh
+// for multi agent local, set the mesh to use 20x20 grid because that is what 
+// was used for training.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "multi_agent_local_refiner.hpp"
+
+#define MFEM_USE_RLLIB
+#ifdef MFEM_USE_RLLIB
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+
+using namespace std;
+using namespace mfem;
+
+#define alpha 2.0/3.0
+
+double exact_vel(const Vector &x)
+{
+   double xv = x(0), yv = x(1);
+   double rv = xv*xv + yv*yv;
+   if (rv > 0) { rv = pow(rv, 0.5); };
+   double theta = atan2(yv, xv);
+   if (theta < 0.0) { theta += 2*M_PI; }
+   return pow(rv, alpha)*sin(alpha*theta);
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "lshape.mesh";
+   int order = 2;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   int jobid = 0;
+   double error_threshold = 0.10;
+   double max_elem_error = 5.0e-3;
+   int refinement_levels = 2;
+
+
+#ifdef MFEM_USE_RLLIB
+   Py_Initialize();
+   import_array(); // numpy init
+#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(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&jobid, "-j", "--jobid",
+                  "slurb_jobid.");
+   args.AddOption(&error_threshold, "-err", "--err",
+                  "Total error fraction for zz or max_elem_error for policy.");
+   args.AddOption(&refinement_levels, "-r", "--ref",
+                 "Refinement levels");
+ 
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   mesh.SetCurvature(2);
+
+   // 4. Since a NURBS mesh can currently only be refined uniformly, we need to
+   //    convert it to a piecewise-polynomial curved mesh. First we refine the
+   //    NURBS mesh a bit more and then project the curvature to quadratic Nodes.
+   if (mesh.NURBSext)
+   {
+      for (int i = 0; i < 2; i++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.SetCurvature(2);
+   }
+   else {
+      // mesh.UniformRefinement();
+      //mesh.UniformRefinement();
+      for (int i = 0; i < refinement_levels; i++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.EnsureNCMesh();
+   }
+
+   // 5. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   // Create 0-order L2 gridfunction to hold errors
+   L2_FECollection fec0(0, dim);
+   FiniteElementSpace fes0(&mesh, &fec0);
+   GridFunction err(&fes0);
+
+   // 6. As in Example 1, we set up bilinear and linear forms corresponding to
+   //    the Laplace problem -\Delta u = 1. We don't assemble the discrete
+   //    problem yet, this will be done in the main loop.
+   BilinearForm a(&fespace);
+   if (pa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      a.SetDiagonalPolicy(Operator::DIAG_ONE);
+   }
+   LinearForm b(&fespace);
+
+   ConstantCoefficient rhs(0.0);
+   ConstantCoefficient one(1.0);
+   FunctionCoefficient exact(exact_vel);
+
+   BilinearFormIntegrator *integ = new DiffusionIntegrator(one);
+   a.AddDomainIntegrator(integ);
+   int int_order = 8;
+   int geom_type = mesh.GetElementBaseGeometry(0);
+   DomainLFIntegrator* dlfi = new DomainLFIntegrator(rhs);
+   dlfi->SetIntRule(&IntRules.Get(geom_type, int_order));
+   b.AddDomainIntegrator(dlfi);
+
+   // 7. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. All boundary attributes will be used for essential (Dirichlet) BC.
+   MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
+               "Boundary attributes required in the mesh.");
+   Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 9. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   socketstream err_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost, visport);
+      err_sock.open(vishost, visport);
+   }
+
+   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //     that uses the ComputeElementFlux method of the DiffusionIntegrator to
+   //     recover a smoothed flux (gradient) that is subtracted from the element
+   //     flux to get an error indicator. We need to supply the space for the
+   //     smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
+   FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
+   ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
+   //KellyErrorEstimator estimator2(*integ, x, flux_fespace);
+   //estimator.SetAnisotropic();
+
+   // 11. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+
+   bool zz = false;
+#if 0
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(error_threshold);
+   zz = true;
+#else
+   MAL_DRLRefiner refiner(x, error_threshold);
+#endif
+
+   bool derefine = false;
+   ThresholdDerefiner derefiner(estimator);
+   derefiner.SetThreshold(0.05);
+   derefiner.SetNCLimit(0);
+
+   string errorfilename;
+   errorfilename = to_string(jobid) + "_lshape_error.txt";
+   ofstream myfile;
+   myfile.open(errorfilename, ofstream::in | ofstream::out | ofstream::app);
+
+
+   // 12. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   const int max_dofs = 200000;
+   for (int it = 0; it < 6; it++)
+   {
+      int cdofs = fespace.GetTrueVSize();
+      cout << "\nAMR iteration " << it << endl;
+      cout << "Number of unknowns: " << cdofs << endl;
+
+      // 13. Assemble the right-hand side.
+      b.Assemble();
+
+      // 14. Set Dirichlet boundary values in the GridFunction x.
+      //     Determine the list of Dirichlet true DOFs in the linear system.
+      Array<int> ess_tdof_list;
+      x.ProjectBdrCoefficient(exact, ess_bdr);
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+      // 15. Assemble the stiffness matrix.
+      a.Assemble();
+      
+      // 16. Create the linear system: eliminate boundary conditions, constrain
+      //     hanging nodes and possibly apply other transformations. The system
+      //     will be solved for true (unconstrained) DOFs only.
+      OperatorPtr A;
+      Vector B, X;
+
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+
+      // 17. Solve the linear system A X = B.
+      if (!pa)
+      {
+#ifndef MFEM_USE_SUITESPARSE
+         // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+         GSSmoother M((SparseMatrix&)(*A));
+         PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
+#else
+         // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+         UMFPackSolver umf_solver;
+         umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+         umf_solver.SetOperator(*A);
+         umf_solver.Mult(B, X);
+#endif
+      }
+      else // Diagonal preconditioning in partial assembly mode.
+      {
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         PCG(*A, M, B, X, 3, 2000, 1e-12, 0.0);
+      }
+
+      // 18. After solving the linear system, reconstruct the solution as a
+      //     finite element GridFunction. Constrained nodes are interpolated
+      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
+      a.RecoverFEMSolution(X, b, x);
+
+      // Compute error against exact solution
+
+      x.ComputeElementL2Errors(exact, err);
+      int int_order = std::max(20 - it, 2*order+1);
+      double error;
+      error = err.Norml2();
+
+      if (derefine) {
+         myfile << error_threshold << " " << cdofs << " " << error << endl;
+      }
+      else {
+         myfile << -error_threshold << " " << cdofs << " " << error << endl;
+      }
+
+      // 19. Send solution by socket to the GLVis server.
+      if (visualization && sol_sock.good())
+      {
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << x << flush;
+      }
+      if (visualization && err_sock.good())
+      {
+         err_sock.precision(8);
+         err_sock << "solution\n" << mesh << err << flush;
+      }
+
+      if (cdofs > max_dofs)
+      {
+         cout << "Reached the maximum number of dofs. Stop." << endl;
+         break;
+      }
+
+      // 20. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+
+      refiner.Apply(mesh);
+      if (refiner.Stop())
+      {
+         cout << "Stopping criterion satisfied. Stop." << endl;
+         break;
+      }
+      fespace.Update();fes0.Update();x.Update();err.Update();
+      a.Update();b.Update();
+
+      if (derefine) {
+         derefiner.Apply(mesh);
+         fespace.Update();fes0.Update();x.Update();err.Update();
+         a.Update();b.Update();
+      }
+e
+      {
+         string solname = to_string(jobid) + "_lshape_amr" + to_string(it) + ".gf";
+         ofstream sol_ofs(solname);
+         x.Save(sol_ofs);
+      }
+
+      {
+         string meshname = to_string(jobid) + "_lshape_amr" + to_string(it) + ".mesh";
+         ofstream mesh_ofs(meshname);
+         mesh_ofs.precision(14);
+         mesh.Print(mesh_ofs);
+      }
+   }
+   myfile.close();
+
+
+   {
+      ofstream sol_ofs("lshape_amr.gf");
+      x.Save(sol_ofs);
+   }
+
+   {
+      ofstream mesh_ofs("lshape_amr.mesh");
+      mesh_ofs.precision(14);
+      mesh.Print(mesh_ofs);
+   }
+
+   return 0;
+}

examples/drl_shock_wave.cpp

@@ -0,0 +1,330 @@
+// Compile with: make drl_shock_wave
+//
+// drl_shock_wave -o 2 -m ../data/inline-quad.mesh
+// for multi agent local, set the mesh to use 20x20 grid because that is what 
+// was used for training.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "drl_shock_wave.hpp"
+#include "multi_agent_local_refiner.hpp"
+
+#define MFEM_USE_RLLIB
+#ifdef MFEM_USE_RLLIB
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+
+using namespace std;
+using namespace mfem;
+
+double bdry_function(const Vector& x)
+{
+   return 0.0; // default
+}
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/inline-quad.mesh";
+   int order = 2;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   int jobid = 0;
+
+#ifdef MFEM_USE_RLLIB
+   Py_Initialize();
+   import_array(); // numpy init
+#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(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&jobid, "-j", "--jobid",
+                  "slurb_jobid.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
+   //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
+   //    the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   mesh.SetCurvature(2);
+
+   // 4. Since a NURBS mesh can currently only be refined uniformly, we need to
+   //    convert it to a piecewise-polynomial curved mesh. First we refine the
+   //    NURBS mesh a bit more and then project the curvature to quadratic Nodes.
+   if (mesh.NURBSext)
+   {
+      for (int i = 0; i < 2; i++)
+      {
+         mesh.UniformRefinement();
+      }
+      mesh.SetCurvature(2);
+   }
+   else {
+      // mesh.UniformRefinement();
+      //mesh.UniformRefinement();
+      mesh.EnsureNCMesh();
+   }
+
+   // 5. Define a finite element space on the mesh. The polynomial order is
+   //    one (linear) by default, but this can be changed on the command line.
+   H1_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   // Create 0-order L2 gridfunction to hold errors
+   L2_FECollection fec0(0, dim);
+   FiniteElementSpace fes0(&mesh, &fec0);
+   GridFunction err(&fes0);
+
+   // 6. As in Example 1, we set up bilinear and linear forms corresponding to
+   //    the Laplace problem -\Delta u = 1. We don't assemble the discrete
+   //    problem yet, this will be done in the main loop.
+   BilinearForm a(&fespace);
+   if (pa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      a.SetDiagonalPolicy(Operator::DIAG_ONE);
+   }
+   LinearForm b(&fespace);
+
+   FunctionCoefficient bdry(bdry_function);
+   FunctionCoefficient rhs(layer2_laplace);
+   FunctionCoefficient exact(layer2_exsol);
+   ConstantCoefficient one(1.0);
+
+   BilinearFormIntegrator *integ = new DiffusionIntegrator(one);
+   a.AddDomainIntegrator(integ);
+   int int_order = 8;
+   int geom_type = mesh.GetElementBaseGeometry(0);
+   DomainLFIntegrator* dlfi = new DomainLFIntegrator(rhs);
+   dlfi->SetIntRule(&IntRules.Get(geom_type, int_order));
+   b.AddDomainIntegrator(dlfi);
+
+   // 7. The solution vector x and the associated finite element grid function
+   //    will be maintained over the AMR iterations. We initialize it to zero.
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. All boundary attributes will be used for essential (Dirichlet) BC.
+   MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
+               "Boundary attributes required in the mesh.");
+   Array<int> ess_bdr(mesh.bdr_attributes.Max());
+   ess_bdr = 1;
+
+   // 9. Connect to GLVis.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   socketstream err_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost, visport);
+      err_sock.open(vishost, visport);
+   }
+
+   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
+   //     that uses the ComputeElementFlux method of the DiffusionIntegrator to
+   //     recover a smoothed flux (gradient) that is subtracted from the element
+   //     flux to get an error indicator. We need to supply the space for the
+   //     smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
+   FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
+   ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
+   //KellyErrorEstimator estimator2(*integ, x, flux_fespace);
+   //estimator.SetAnisotropic();
+
+   // 11. A refiner selects and refines elements based on a refinement strategy.
+   //     The strategy here is to refine elements with errors larger than a
+   //     fraction of the maximum element error. Other strategies are possible.
+   //     The refiner will call the given error estimator.
+
+#if 0
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.10);
+#else
+   MAL_DRLRefiner refiner(x);
+   //DRLRefiner refiner(x);
+#endif
+  
+#if 0
+   x.ProjectCoefficient(exact);
+   refiner.Apply(mesh);
+   fespace.Update();
+   fes0.Update();
+   x.Update();
+   x.ProjectCoefficient(exact);
+   refiner.Apply(mesh);
+   fespace.Update();
+   fes0.Update();
+   x.Update();
+   {
+      ofstream sol_ofs("sol.gf");
+      x.Save(sol_ofs);
+   }
+   {
+      ofstream mesh_ofs("amr.mesh");
+      mesh_ofs.precision(14);
+      mesh.Print(mesh_ofs);
+   }
+   MFEM_ABORT(" ");
+#endif
+
+   // 12. The main AMR loop. In each iteration we solve the problem on the
+   //     current mesh, visualize the solution, and refine the mesh.
+   const int max_dofs = 20000;
+   for (int it = 0; it < 4; it++)
+   {
+      int cdofs = fespace.GetTrueVSize();
+      cout << "\nAMR iteration " << it << endl;
+      cout << "Number of unknowns: " << cdofs << endl;
+
+      // 13. Assemble the right-hand side.
+      b.Assemble();
+
+      // 14. Set Dirichlet boundary values in the GridFunction x.
+      //     Determine the list of Dirichlet true DOFs in the linear system.
+      Array<int> ess_tdof_list;
+      x.ProjectBdrCoefficient(exact, ess_bdr);
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+      // 15. Assemble the stiffness matrix.
+      a.Assemble();
+      
+      // 16. Create the linear system: eliminate boundary conditions, constrain
+      //     hanging nodes and possibly apply other transformations. The system
+      //     will be solved for true (unconstrained) DOFs only.
+      OperatorPtr A;
+      Vector B, X;
+
+      const int copy_interior = 1;
+      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+
+      // 17. Solve the linear system A X = B.
+      if (!pa)
+      {
+#ifndef MFEM_USE_SUITESPARSE
+         // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+         GSSmoother M((SparseMatrix&)(*A));
+         PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
+#else
+         // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+         UMFPackSolver umf_solver;
+         umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+         umf_solver.SetOperator(*A);
+         umf_solver.Mult(B, X);
+#endif
+      }
+      else // Diagonal preconditioning in partial assembly mode.
+      {
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         PCG(*A, M, B, X, 3, 2000, 1e-12, 0.0);
+      }
+
+      // 18. After solving the linear system, reconstruct the solution as a
+      //     finite element GridFunction. Constrained nodes are interpolated
+      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
+      a.RecoverFEMSolution(X, b, x);
+
+      // Compute error against exact solution
+
+      x.ComputeElementL2Errors(exact, err);
+
+      // 19. Send solution by socket to the GLVis server.
+      if (visualization && sol_sock.good())
+      {
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << x << flush;
+      }
+      if (visualization && err_sock.good())
+      {
+         err_sock.precision(8);
+         err_sock << "solution\n" << mesh << err << flush;
+      }
+
+      if (cdofs > max_dofs || it == 3)
+      {
+         cout << "Reached the maximum number of dofs. Stop." << endl;
+         break;
+      }
+
+      // 20. Call the refiner to modify the mesh. The refiner calls the error
+      //     estimator to obtain element errors, then it selects elements to be
+      //     refined and finally it modifies the mesh. The Stop() method can be
+      //     used to determine if a stopping criterion was met.
+
+      refiner.Apply(mesh);
+      if (refiner.Stop())
+      {
+         cout << "Stopping criterion satisfied. Stop." << endl;
+         break;
+      }
+
+      // 21. Update the space to reflect the new state of the mesh. Also,
+      //     interpolate the solution x so that it lies in the new space but
+      //     represents the same function. This saves solver iterations later
+      //     since we'll have a good initial guess of x in the next step.
+      //     Internally, FiniteElementSpace::Update() calculates an
+      //     interpolation matrix which is then used by GridFunction::Update().
+      fespace.Update();
+      fes0.Update();
+      x.Update();
+      err.Update();
+
+      {
+         string solname = to_string(jobid) + "_amrsol" + to_string(it) + ".gf";
+         ofstream sol_ofs(solname);
+         x.Save(sol_ofs);
+      }
+
+      {
+         string meshname = to_string(jobid) + "_amr" + to_string(it) + ".mesh";
+         ofstream mesh_ofs(meshname);
+         mesh_ofs.precision(14);
+         mesh.Print(mesh_ofs);
+      }
+
+      // 22. Inform also the bilinear and linear forms that the space has
+      //     changed.
+      a.Update();
+      b.Update();
+   }
+
+
+   {
+      ofstream sol_ofs("amrsol.gf");
+      x.Save(sol_ofs);
+   }
+
+   {
+      ofstream mesh_ofs("amr.mesh");
+      mesh_ofs.precision(14);
+      mesh.Print(mesh_ofs);
+   }
+
+   return 0;
+}

examples/drl_shock_wave.hpp

@@ -0,0 +1,136 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace mfem;
+using namespace std;
+
+const double alpha =  200.0; // standard params
+const double center = -0.05;
+const double radius = 0.7;
+template<typename T> T sqr(T x) { return x*x; }
+double layer2_exsol(Vector &p)
+{
+   double x = p(0), y = p(1);
+   double r = sqrt(sqr(x - center) + sqr(y - center));
+   return atan(alpha * (r - radius));
+}
+
+void layer2_exgrad(const Vector &p, Vector &grad)
+{
+   double x = p(0), y = p(1);
+   double r = sqrt(sqr(x - center) + sqr(y - center));
+   double u = r * (sqr(alpha) * sqr(r - radius) + 1);
+   grad(0) = alpha * (x - center) / u;
+   grad(1) = alpha * (y - center) / u;
+}
+
+double layer2_laplace(Vector &p)
+{
+   double x = p(0), y = p(1);
+   double r = sqr(y - center) + sqr(x - center);
+   double u = sqr(alpha) * sqr(sqrt(r) - radius) + 1;
+
+   return 2 * pow(alpha,3) * (sqrt(r) - radius) * sqr(y - center) / (r * sqr(u))
+          + alpha * sqr(y - center) / (pow(r, 1.5) * u)
+          - 2 * alpha / (sqrt(r) * u)
+          + 2 * pow(alpha,3) * (sqrt(r) - radius) * sqr(x - center) / (r * sqr(u))
+          + alpha * sqr(x - center) / (pow(r, 1.5) * u);
+}
+
+static double safeSqrt(double x)
+{
+   if (x < 0.0)
+      return -sqrt(-x);
+   else
+      return sqrt(x);
+}
+
+double CalculateH10Error(GridFunction *sol, VectorCoefficient *exgrad,
+                         Array<double> *elemError, Array<int> *elemRef,
+                         int intOrder)
+{
+   const FiniteElementSpace *fes = sol->FESpace();
+   Mesh* mesh = fes->GetMesh();
+
+   Vector e_grad, a_grad, el_dofs, q_grad;
+   DenseMatrix dshape, dshapet, Jinv;
+   Array<int> vdofs;
+   const FiniteElement *fe;
+   ElementTransformation *transf;
+
+   int dim = mesh->Dimension();
+   e_grad.SetSize(dim);
+   a_grad.SetSize(dim);
+   q_grad.SetSize(dim);
+   Jinv.SetSize(dim);
+
+   double error = 0.0;
+   if (elemError) elemError->SetSize(mesh->GetNE());
+   if (elemRef)   elemRef->SetSize(mesh->GetNE());
+
+   for (int i = 0; i < mesh->GetNE(); i++)
+   {
+      fe = fes->GetFE(i);
+      int fdof = fe->GetDof();
+      transf = mesh->GetElementTransformation(i);
+      el_dofs.SetSize(fdof);
+      dshape.SetSize(fdof, dim);
+      dshapet.SetSize(fdof, dim);
+
+      fes->GetElementVDofs(i, vdofs);
+      for (int k = 0; k < fdof; k++)
+         if (vdofs[k] >= 0)
+            el_dofs(k) =  (*sol)(vdofs[k]);
+         else
+            el_dofs(k) = -(*sol)(-1-vdofs[k]);
+
+      const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intOrder);
+
+      double el_err = 0.0, a_dxyz[3] = { 0, 0, 0 };
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(j);
+
+         transf->SetIntPoint(&ip);
+         CalcInverse(transf->Jacobian(), Jinv);
+         double w = ip.weight * transf->Weight();
+
+         exgrad->Eval(e_grad, *transf, ip);
+
+         fe->CalcDShape(ip, dshape);
+         Mult(dshape, Jinv, dshapet);
+         dshapet.MultTranspose(el_dofs, a_grad);
+
+         e_grad -= a_grad;
+         el_err += w * (e_grad * e_grad);
+
+         transf->Jacobian().MultTranspose(e_grad, q_grad);
+         for (int k = 0; k < dim; k++)
+         {
+            a_dxyz[k] += w * (q_grad[k] * q_grad[k]);
+         }
+      }
+
+      error += el_err;
+      if (elemError)
+         (*elemError)[i] = sqrt(fabs(el_err));
+
+      if (elemRef)
+      {
+         double sum = 0;
+         for (int k = 0; k < dim; k++)
+            sum += a_dxyz[k];
+
+         const double thresh = 0.2 * 3/dim;
+         int ref = 0;
+         for (int k = 0; k < dim; k++)
+            if (a_dxyz[k] / sum > thresh)
+               ref |= (1 << k);
+
+         (*elemRef)[i] = ref;
+      }
+   }
+
+   return safeSqrt(error);
+}

examples/ex6.cpp

@@ -43,21 +43,159 @@
 #include <fstream>
 #include <iostream>
 
+#define MFEM_USE_RLLIB
+#ifdef MFEM_USE_RLLIB
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+
 using namespace std;
 using namespace mfem;
 
+int problem; // problem number, controls source term and bdry condition
+
+// Returns either rhs function f, or exact solution u.
+
+double rhs_function(const Vector &x, bool exact = false)
+{
+   int dim = x.Size();
+
+   if (0 == problem) {
+
+      // NIST: "Peak 2D" problem
+      // u = e^(-alpha*(x^2+y^2)), alpha = 1e3 or alpha = 1e5
+
+      double xc = 0.5;
+      double yc = 0.5;
+
+      double x0 = x(0)-xc;
+      double y0 = x(1)-yc;
+
+      double alpha = 1000;
+
+      double xx = x0*x0;
+      double yy = y0*y0;
+
+      if (exact) {
+
+         double u = exp(-alpha*(xx+yy));
+
+         return u;
+      }
+      else {
+
+         double f =
+            alpha*(4.*alpha*xx -2.0)*exp(-alpha*(xx+yy)) +
+            alpha*(4.*alpha*yy -2.0)*exp(-alpha*(xx+yy));
+
+         return -f;
+      }
+   }
+   if (1 == problem) {
+
+      // NIST "arctan circular wavefront" problem, w/ minor
+      // modifications.
+
+      double r0 = 0.25;
+      double a = 100.0;
+      double h = 1./3.;
+      double c = 0.5;
+
+      double x0 = 0.5;
+      double y0 = 0.5;
+
+      double dx = x(0)-x0;
+      double dy = x(1)-y0;
+      double dxdx = dx*dx;
+      double dydy = dy*dy;
+      double r = sqrt(dxdx+dydy);
+
+      if (exact) {
+         double u = c +h*atan(a*(r-r0));
+         return u;
+      }
+      else {
+         double aa = a*a;
+         double dr = r-r0;
+         double drdr = dr*dr;
+         double rr = r*r;
+         double t = 1+aa*drdr;
+         double tt = t*t;
+
+         double fx =
+            a/(t*r) -
+            a*dxdx/(t*rr*r) -
+            2*aa*a*dxdx*dr/(tt*rr);
+
+         double fy =
+            a/(t*r) -
+            a*dydy/(t*rr*r) -
+            2*aa*a*dydy*dr/(tt*rr);
+
+         return -h*(fx+fy);
+      }
+   }
+   if (2 == problem) {
+
+      // cross-shaped source
+
+      double x0 = x(0)-0.5;
+      double y0 = x(1)-0.5;
+
+      double w1 = 0.04;
+      double w2 = 0.20;
+
+      if (x(0) > 0.5-w1 && x(0) < 0.5+w1 &&
+          x(1) > 0.5-w2 && x(1) < 0.5+w2 ) {
+         return 100.0;
+      }
+      if (x(1) > 0.5-w1 && x(1) < 0.5+w1 &&
+          x(0) > 0.5-w2 && x(0) < 0.5+w2 ) {
+         return 100.0;
+      }
+      return 0.0;
+   }
+
+   return 1.0; // default
+}
+
+double exact_soln(const Vector& x)
+{
+   return rhs_function(x, true);
+}
+
+double bdry_function(const Vector& x)
+{
+   if (problem == 0) {
+      return 0.0;
+   }
+   if (problem == 1) {
+      return 1.0;
+   }
+
+   return 0.0; // default
+}
+
 int main(int argc, char *argv[])
 {
    // 1. Parse command-line options.
+   problem = 0;
    const char *mesh_file = "../data/star.mesh";
    int order = 1;
    bool pa = false;
    const char *device_config = "cpu";
    bool visualization = true;
 
+#ifdef MFEM_USE_RLLIB
+   Py_Initialize();
+   import_array(); // numpy init
+#endif
+
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
                   "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in rhs_function().");
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree).");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
@@ -98,12 +236,22 @@ int main(int argc, char *argv[])
       }
       mesh.SetCurvature(2);
    }
+   else {
+      // mesh.UniformRefinement();
+      //mesh.UniformRefinement();
+      mesh.EnsureNCMesh();
+   }
 
    // 5. Define a finite element space on the mesh. The polynomial order is
    //    one (linear) by default, but this can be changed on the command line.
    H1_FECollection fec(order, dim);
    FiniteElementSpace fespace(&mesh, &fec);
 
+   // Create 0-order L2 gridfunction to hold errors
+   L2_FECollection fec0(0, dim);
+   FiniteElementSpace fes0(&mesh, &fec0);
+   GridFunction err(&fes0);
+
    // 6. As in Example 1, we set up bilinear and linear forms corresponding to
    //    the Laplace problem -\Delta u = 1. We don't assemble the discrete
    //    problem yet, this will be done in the main loop.
@@ -115,12 +263,14 @@ int main(int argc, char *argv[])
    }
    LinearForm b(&fespace);
 
+   FunctionCoefficient rhs(rhs_function);
+   FunctionCoefficient exact(exact_soln);
+   FunctionCoefficient bdry(bdry_function);
    ConstantCoefficient one(1.0);
-   ConstantCoefficient zero(0.0);
 
    BilinearFormIntegrator *integ = new DiffusionIntegrator(one);
    a.AddDomainIntegrator(integ);
-   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.AddDomainIntegrator(new DomainLFIntegrator(rhs));
 
    // 7. The solution vector x and the associated finite element grid function
    //    will be maintained over the AMR iterations. We initialize it to zero.
@@ -137,9 +287,11 @@ int main(int argc, char *argv[])
    char vishost[] = "localhost";
    int  visport   = 19916;
    socketstream sol_sock;
+   socketstream err_sock;
    if (visualization)
    {
       sol_sock.open(vishost, visport);
+      err_sock.open(vishost, visport);
    }
 
    // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
@@ -149,18 +301,24 @@ int main(int argc, char *argv[])
    //     smoothed flux: an (H1)^sdim (i.e., vector-valued) space is used here.
    FiniteElementSpace flux_fespace(&mesh, &fec, sdim);
    ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
-   estimator.SetAnisotropic();
+   //KellyErrorEstimator estimator2(*integ, x, flux_fespace);
+   //estimator.SetAnisotropic();
 
    // 11. A refiner selects and refines elements based on a refinement strategy.
    //     The strategy here is to refine elements with errors larger than a
    //     fraction of the maximum element error. Other strategies are possible.
    //     The refiner will call the given error estimator.
+
+#if 1
    ThresholdRefiner refiner(estimator);
-   refiner.SetTotalErrorFraction(0.7);
+   refiner.SetTotalErrorFraction(0.10);
+#else
+   DRLRefiner refiner(x);
+#endif
 
    // 12. The main AMR loop. In each iteration we solve the problem on the
    //     current mesh, visualize the solution, and refine the mesh.
-   const int max_dofs = 50000;
+   const int max_dofs = 1000;
    for (int it = 0; ; it++)
    {
       int cdofs = fespace.GetTrueVSize();
@@ -173,7 +331,7 @@ int main(int argc, char *argv[])
       // 14. Set Dirichlet boundary values in the GridFunction x.
       //     Determine the list of Dirichlet true DOFs in the linear system.
       Array<int> ess_tdof_list;
-      x.ProjectBdrCoefficient(zero, ess_bdr);
+      x.ProjectBdrCoefficient(bdry, ess_bdr);
       fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
 
       // 15. Assemble the stiffness matrix.
@@ -214,14 +372,23 @@ int main(int argc, char *argv[])
       //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
       a.RecoverFEMSolution(X, b, x);
 
+      // Compute error against exact solution
+
+      x.ComputeElementL2Errors(exact, err);
+
       // 19. Send solution by socket to the GLVis server.
       if (visualization && sol_sock.good())
       {
          sol_sock.precision(8);
          sol_sock << "solution\n" << mesh << x << flush;
       }
+      if (visualization && err_sock.good())
+      {
+         err_sock.precision(8);
+         err_sock << "solution\n" << mesh << err << flush;
+      }
 
-      if (cdofs > max_dofs)
+      if (cdofs > max_dofs || it == 3)
       {
          cout << "Reached the maximum number of dofs. Stop." << endl;
          break;
@@ -231,6 +398,7 @@ int main(int argc, char *argv[])
       //     estimator to obtain element errors, then it selects elements to be
       //     refined and finally it modifies the mesh. The Stop() method can be
       //     used to determine if a stopping criterion was met.
+
       refiner.Apply(mesh);
       if (refiner.Stop())
       {
@@ -245,7 +413,9 @@ int main(int argc, char *argv[])
       //     Internally, FiniteElementSpace::Update() calculates an
       //     interpolation matrix which is then used by GridFunction::Update().
       fespace.Update();
+      fes0.Update();
       x.Update();
+      err.Update();
 
       // 22. Inform also the bilinear and linear forms that the space has
       //     changed.

examples/ex91.cpp

@@ -1,672 +0,0 @@
-#include"mfem.hpp"
-
-#include<memory>
-#include<iostream>
-#include<fstream>
-
-namespace mfem {
-
-
-class LinDiffQFunc
-{
-public:
-    LinDiffQFunc(mfem::Coefficient& dd, mfem::Coefficient& ll, double gg_,
-                 double pp0_, double pp1_):diff(dd),load(ll),gg(gg_),pp0(pp0_),pp1(pp1_)
-    {
-
-    }
-
-    double QEnergy(ElementTransformation &T,
-                  const IntegrationPoint &ip,
-                  mfem::Vector& param, mfem::Vector& uu)
-    {
-        double dd=diff.Eval(T,ip);
-        double ll=load.Eval(T,ip);
-
-        double rho0=param[0];
-        double rho1=param[1];
-        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
-
-        double rez = 0.5*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])*fd
-                        + 0.5*gg*uu[3]*uu[3] -uu[3]*ll;
-
-        return rez;
-    }
-
-    void QResidual(ElementTransformation &T,
-                   const IntegrationPoint &ip,
-                   mfem::Vector& param, mfem::Vector& uu, mfem::Vector& rr)
-    {
-        rr.SetSize(4);
-        double dd=diff.Eval(T,ip);
-        double ll=load.Eval(T,ip);
-
-        double rho0=param[0];
-        double rho1=param[1];
-
-        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
-
-        rr[0]=uu[0]*fd;
-        rr[1]=uu[1]*fd;
-        rr[2]=uu[2]*fd;
-        rr[3]=gg*uu[3]-ll;
-    }
-
-    void AQResidual(ElementTransformation &T,
-                    const IntegrationPoint &ip,
-                    mfem::Vector& param,
-                    mfem::Vector& uu, mfem::Vector& aa, mfem::Vector& rr)
-    {
-         rr.SetSize(2);
-         double dd=diff.Eval(T,ip);
-         double ll=load.Eval(T,ip);
-
-         double rho0=param[0];
-         double rho1=param[1];
-
-         double fd0=dd*pp0*std::pow(rho0,pp0-1.0)*std::pow(rho1,pp1);
-         double fd1=dd*std::pow(rho0,pp0)*pp1*std::pow(rho1,pp1-1.0);
-
-         rr[0] = (aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd0;
-         rr[1] = (aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd1;
-
-    }
-
-    void QGradResidual(ElementTransformation &T,
-                       const IntegrationPoint &ip,
-                       mfem::Vector& param, mfem::Vector& uu, mfem::DenseMatrix& hh)
-    {
-        hh.SetSize(4);
-        double dd=diff.Eval(T,ip);
-        //double ll=load.Eval(T,ip);
-
-
-        double rho0=param[0];
-        double rho1=param[1];
-
-        double fd=dd*std::pow(rho0,pp0)*std::pow(rho1,pp1);
-        hh=0.0;
-
-        hh(0,0)=fd;
-        hh(1,1)=fd;
-        hh(2,2)=fd;
-        hh(3,3)=gg;
-    }
-
-
-private:
-    mfem::Coefficient& diff;
-    mfem::Coefficient& load;
-    double gg;
-    double pp0;
-    double pp1;
-};
-
-
-class PrmBlockLSFEMDiffusion: public PrmBlockNonlinearFormIntegrator
-{
-public:
-    PrmBlockLSFEMDiffusion(LinDiffQFunc& qfun_)
-    {
-        qfunc=&qfun_;
-    }
-
-    /// Compute the local energy
-    virtual double GetElementEnergy(const Array<const FiniteElement *>&el,
-                                    const Array<const FiniteElement *>&pel,
-                                    ElementTransformation &Tr,
-                                    const Array<const Vector *>&elfun,
-                                    const Array<const Vector *>&pelfun)
-    {
-        int dof_u0 = el[0]->GetDof();
-        int dof_r0 = pel[0]->GetDof();
-        int dof_r1 = pel[1]->GetDof();
-
-        int dim = el[0]->GetDim();
-        int spaceDim = Tr.GetSpaceDim();
-        if (dim != spaceDim)
-        {
-            mfem::mfem_error(" PrmBlockLSFEMDiffusion::GetElementEnergy"
-                             " is not defined on manifold meshes");
-        }
-
-        //shape functions
-        Vector shu0(dof_u0);
-        Vector shr0(dof_r0);
-        Vector shr1(dof_r1);
-        DenseMatrix dsu0(dof_u0,dim);
-        DenseMatrix B(dof_u0, 4);
-        B=0.0;
-
-        double w;
-
-        Vector param(2); param=0.0;
-        Vector uu(4); uu=0.0;
-
-        double energy =0.0;
-
-        const IntegrationRule *ir = nullptr;
-        if(ir==nullptr){
-            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                                +pel[0]->GetOrder()+pel[1]->GetOrder();
-            ir=&IntRules.Get(Tr.GetGeometryType(),order);
-        }
-
-        for (int i = 0; i < ir->GetNPoints(); i++)
-        {
-            const IntegrationPoint &ip = ir->IntPoint(i);
-            Tr.SetIntPoint(&ip);
-            w=Tr.Weight();
-            w = ip.weight * w;
-
-            el[0]->CalcPhysDShape(Tr,dsu0);
-            el[0]->CalcPhysShape(Tr,shu0);
-            pel[0]->CalcPhysShape(Tr,shr0);
-            pel[1]->CalcPhysShape(Tr,shr1);
-
-            param[0]=shr0*(*pelfun[0]);
-            param[1]=shr1*(*pelfun[1]);
-
-            //set the matrix B
-            for(int jj=0;jj<dim;jj++)
-            {
-                B.SetCol(jj,dsu0.GetColumn(jj));
-            }
-            B.SetCol(3,shu0);
-            B.MultTranspose(*elfun[0],uu);
-            energy=energy+w * qfunc->QEnergy(Tr,ip,param,uu);
-        }
-        return energy;
-    }
-
-    /// Perform the local action of the BlockNonlinearFormIntegrator
-    virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
-                                       const Array<const FiniteElement *>&pel,
-                                       ElementTransformation &Tr,
-                                       const Array<const Vector *> &elfun,
-                                       const Array<const Vector *>&pelfun,
-                                       const Array<Vector *> &elvec)
-    {
-        int dof_u0 = el[0]->GetDof();
-        int dof_r0 = pel[0]->GetDof();
-        int dof_r1 = pel[1]->GetDof();
-
-        int dim = el[0]->GetDim();
-
-        elvec[0]->SetSize(dof_u0);
-        *elvec[0]=0.0;
-        int spaceDim = Tr.GetSpaceDim();
-        if (dim != spaceDim)
-        {
-            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
-                             " is not defined on manifold meshes");
-        }
-
-        //shape functions
-        Vector shu0(dof_u0);
-        Vector shr0(dof_r0);
-        Vector shr1(dof_r1);
-        DenseMatrix dsu0(dof_u0,dim);
-        DenseMatrix B(dof_u0, 4);
-        B=0.0;
-
-        double w;
-
-        Vector param(2); param=0.0;
-        Vector uu(4); uu=0.0;
-        Vector rr;
-        Vector lvec; lvec.SetSize(dof_u0);
-
-        const IntegrationRule *ir = nullptr;
-        if(ir==nullptr){
-            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                                +pel[0]->GetOrder()+pel[1]->GetOrder();
-            ir=&IntRules.Get(Tr.GetGeometryType(),order);
-        }
-
-        for (int i = 0; i < ir->GetNPoints(); i++)
-        {
-            const IntegrationPoint &ip = ir->IntPoint(i);
-            Tr.SetIntPoint(&ip);
-            w=Tr.Weight();
-            w = ip.weight * w;
-
-            el[0]->CalcPhysDShape(Tr,dsu0);
-            el[0]->CalcPhysShape(Tr,shu0);
-            pel[0]->CalcPhysShape(Tr,shr0);
-            pel[1]->CalcPhysShape(Tr,shr1);
-
-            param[0]=shr0*(*pelfun[0]);
-            param[1]=shr1*(*pelfun[1]);
-
-            //set the matrix B
-            for(int jj=0;jj<dim;jj++)
-            {
-                B.SetCol(jj,dsu0.GetColumn(jj));
-            }
-            B.SetCol(3,shu0);
-            B.MultTranspose(*elfun[0],uu);
-            qfunc->QResidual(Tr,ip,param, uu, rr);
-
-            B.Mult(rr,lvec);
-            elvec[0]->Add(w,lvec);
-        }
-
-    }
-
-    virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
-                                    const Array<const FiniteElement *> &el2,
-                                    const Array<const FiniteElement *> &pel1,
-                                    const Array<const FiniteElement *> &pel2,
-                                    FaceElementTransformations &Tr,
-                                    const Array<const Vector *> &elfun,
-                                    const Array<const Vector *>&pelfun,
-                                    const Array<Vector *> &elvect)
-    {
-
-    }
-
-    /// Assemble the local gradient matrix
-    virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
-                                     const Array<const FiniteElement *>&pel,
-                                     ElementTransformation &Tr,
-                                     const Array<const Vector *> &elfun,
-                                     const Array<const Vector *>&pelfun,
-                                     const Array2D<DenseMatrix *> &elmats)
-    {
-        int dof_u0 = el[0]->GetDof();
-        int dof_r0 = pel[0]->GetDof();
-        int dof_r1 = pel[1]->GetDof();
-
-        int dim = el[0]->GetDim();
-
-        //elmats[0]->Size(dof_u0, dof_u0);
-        //*elmats[0]=0.0;
-
-        DenseMatrix* K=elmats(0,0);
-        K->SetSize(dof_u0,dof_u0);
-        (*K)=0.0;
-
-        int spaceDim = Tr.GetSpaceDim();
-        if (dim != spaceDim)
-        {
-            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
-                             " is not defined on manifold meshes");
-        }
-
-        //shape functions
-        Vector shu0(dof_u0);
-        Vector shr0(dof_r0);
-        Vector shr1(dof_r1);
-        DenseMatrix dsu0(dof_u0,dim);
-        DenseMatrix B(dof_u0, 4);
-        DenseMatrix A(dof_u0, 4);
-        B=0.0;
-
-        double w;
-
-        Vector param(2); param=0.0;
-        Vector uu(4); uu=0.0;
-        DenseMatrix hh;
-        Vector lvec; lvec.SetSize(dof_u0);
-
-        const IntegrationRule *ir = nullptr;
-        if(ir==nullptr){
-            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                                +pel[0]->GetOrder()+pel[1]->GetOrder();
-            ir=&IntRules.Get(Tr.GetGeometryType(),order);
-        }
-
-        for (int i = 0; i < ir->GetNPoints(); i++)
-        {
-            const IntegrationPoint &ip = ir->IntPoint(i);
-            Tr.SetIntPoint(&ip);
-            w = Tr.Weight();
-            w = ip.weight * w;
-
-            el[0]->CalcPhysDShape(Tr,dsu0);
-            el[0]->CalcPhysShape(Tr,shu0);
-            pel[0]->CalcPhysShape(Tr,shr0);
-            pel[1]->CalcPhysShape(Tr,shr1);
-
-            param[0]=shr0*(*pelfun[0]);
-            param[1]=shr1*(*pelfun[1]);
-
-            //set the matrix B
-            for(int jj=0;jj<dim;jj++)
-            {
-                B.SetCol(jj,dsu0.GetColumn(jj));
-            }
-            B.SetCol(3,shu0);
-            B.MultTranspose(*elfun[0],uu);
-            qfunc->QGradResidual(Tr,ip,param,uu,hh);
-            Mult(B,hh,A);
-            AddMult_a_ABt(w,A,B,*K);
-        }
-    }
-
-    virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
-                                  const Array<const FiniteElement *>&el2,
-                                  const Array<const FiniteElement *> &pel1,
-                                  const Array<const FiniteElement *> &pel2,
-                                  FaceElementTransformations &Tr,
-                                  const Array<const Vector *> &elfun,
-                                  const Array<const Vector *>&pelfun,
-                                  const Array2D<DenseMatrix *> &elmats)
-    {
-
-    }
-
-    virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                          const Array<const FiniteElement *> &pel,
-                                          ElementTransformation &Tr,
-                                          const Array<const Vector *> &elfun,
-                                          const Array<const Vector *> &alfun,
-                                          const Array<const Vector *> &pelfun,
-                                          const Array<Vector *> &elvec)
-    {
-        int dof_u0 = el[0]->GetDof();
-        int dof_r0 = pel[0]->GetDof();
-        int dof_r1 = pel[1]->GetDof();
-
-        int dim = el[0]->GetDim();
-
-        Vector& e0 = *(elvec[0]);
-        Vector& e1 = *(elvec[1]);
-
-        e0.SetSize(dof_r0);
-        e0=0.0;
-        e1.SetSize(dof_r1);
-        e1=0.0;
-
-        int spaceDim = Tr.GetSpaceDim();
-        if (dim != spaceDim)
-        {
-            mfem::mfem_error(" PrmBlockLSFEMDiffusion::AssembleElementVector"
-                             " is not defined on manifold meshes");
-        }
-
-        //shape functions
-        Vector shu0(dof_u0);
-        Vector shr0(dof_r0);
-        Vector shr1(dof_r1);
-        DenseMatrix dsu0(dof_u0,dim);
-        DenseMatrix B(dof_u0, 4);
-        B=0.0;
-
-        double w;
-
-        Vector param(2); param=0.0;
-        Vector uu(4); uu=0.0;
-        Vector aa(4); aa=0.0;
-        Vector rr;
-        Vector lvec0; lvec0.SetSize(dof_r0);
-        Vector lvec1; lvec1.SetSize(dof_r1);
-
-        const IntegrationRule *ir = nullptr;
-        if(ir==nullptr){
-            int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                                +pel[0]->GetOrder()+pel[1]->GetOrder();
-            ir=&IntRules.Get(Tr.GetGeometryType(),order);
-        }
-
-        for (int i = 0; i < ir->GetNPoints(); i++)
-        {
-            const IntegrationPoint &ip = ir->IntPoint(i);
-            Tr.SetIntPoint(&ip);
-            w=Tr.Weight();
-            w = ip.weight * w;
-
-            el[0]->CalcPhysDShape(Tr,dsu0);
-            el[0]->CalcPhysShape(Tr,shu0);
-            pel[0]->CalcPhysShape(Tr,shr0);
-            pel[1]->CalcPhysShape(Tr,shr1);
-
-            param[0]=shr0*(*pelfun[0]);
-            param[1]=shr1*(*pelfun[1]);
-
-            //set the matrix B
-            for(int jj=0;jj<dim;jj++)
-            {
-                B.SetCol(jj,dsu0.GetColumn(jj));
-            }
-            B.SetCol(3,shu0);
-            B.MultTranspose(*elfun[0],uu);
-            B.MultTranspose(*alfun[0],aa);
-
-            qfunc->AQResidual(Tr, ip, param, uu, aa, rr);
-
-            lvec0=shr0;
-            lvec0*=rr[0];
-            lvec1=shr1;
-            lvec1*=rr[1];
-
-            e0.Add(w,lvec0);
-            e1.Add(w,lvec1);
-        }
-
-    }
-
-    virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
-                                       const Array<const FiniteElement *> &el2,
-                                       const Array<const FiniteElement *> &pel1,
-                                       const Array<const FiniteElement *> &pel2,
-                                       FaceElementTransformations &Tr,
-                                       const Array<const Vector *> &elfun,
-                                       const Array<const Vector *> &alfun,
-                                       const Array<const Vector *> &pelfun,
-                                       const Array<Vector *> &elvect)
-    {
-
-    }
-
-
-
-private:
-    LinDiffQFunc* qfunc;
-};
-
-}
-
-int main(int argc, char *argv[])
-{
-    const char *mesh_file = "../../data/beam-tet.mesh";
-    int ser_ref_levels = 1;
-    int order = 2;
-    bool visualization = true;
-    double newton_rel_tol = 1e-4;
-    double newton_abs_tol = 1e-6;
-    int newton_iter = 10;
-    int print_level = 0;
-
-    mfem::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(&order,
-                      "-o",
-                      "--order",
-                      "Order (degree) of the finite elements.");
-    args.AddOption(&visualization,
-                      "-vis",
-                      "--visualization",
-                      "-no-vis",
-                      "--no-visualization",
-                      "Enable or disable GLVis visualization.");
-    args.AddOption(&newton_rel_tol,
-                      "-rel",
-                      "--relative-tolerance",
-                      "Relative tolerance for the Newton solve.");
-    args.AddOption(&newton_abs_tol,
-                      "-abs",
-                      "--absolute-tolerance",
-                      "Absolute tolerance for the Newton solve.");
-    args.AddOption(&newton_iter,
-                      "-it",
-                      "--newton-iterations",
-                      "Maximum iterations for the Newton solve.");
-
-    args.Parse();
-    if (!args.Good())
-    {
-        args.PrintUsage(std::cout);
-        return 1;
-    }
-    args.PrintOptions(std::cout);
-    // 3. Read the (serial) mesh from the given mesh file on all processors.  We
-    //    can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
-    //    with the same code.
-    mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
-    int dim = mesh->Dimension();
-
-    // 4. Refine the mesh in serial to increase the resolution. In this example
-    //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
-    //    a command-line parameter.
-    for (int lev = 0; lev < ser_ref_levels; lev++)
-    {
-        mesh->UniformRefinement();
-    }
-
-    /// Define the q-function
-    mfem::ConstantCoefficient* dc=new mfem::ConstantCoefficient(1.0);
-    mfem::ConstantCoefficient* lc=new mfem::ConstantCoefficient(1.0);
-    mfem::LinDiffQFunc* qfun=new mfem::LinDiffQFunc(*dc,*lc,1.0,1.0,1.0);
-
-    mfem::H1_FECollection fec00(order, dim);
-    mfem::L2_FECollection fec01(order, dim);
-    mfem::FiniteElementSpace* bfes00=new mfem::FiniteElementSpace(mesh,&fec00,1,mfem::Ordering::byVDIM);
-    mfem::FiniteElementSpace* pfes00=new mfem::FiniteElementSpace(mesh,&fec00,1,mfem::Ordering::byVDIM);
-    mfem::FiniteElementSpace* pfes01=new mfem::FiniteElementSpace(mesh,&fec01,1,mfem::Ordering::byVDIM);
-
-    /// Define parametric nonlinear form
-    mfem::Array<mfem::FiniteElementSpace*> bfes;
-    mfem::Array<mfem::FiniteElementSpace*> pfes;
-
-    bfes.Append(bfes00);
-    pfes.Append(pfes00);
-    pfes.Append(pfes01);
-
-    mfem::PrmBlockNonlinearForm* nf=new mfem::PrmBlockNonlinearForm(bfes,pfes);
-    nf->AddDomainIntegrator(new mfem::PrmBlockLSFEMDiffusion(*qfun));
-
-    /// Define the grid functions
-    mfem::GridFunction* bgf00=new mfem::GridFunction(bfes00);
-    mfem::GridFunction* pgf00=new mfem::GridFunction(pfes00);
-    mfem::GridFunction* pgf01=new mfem::GridFunction(pfes01);
-    mfem::GridFunction* ggf00=new mfem::GridFunction(pfes00);
-    mfem::GridFunction* ggf01=new mfem::GridFunction(pfes01);
-
-    *bgf00=0.0;
-    *pgf00=1.0;
-    *pgf01=1.0;
-
-    mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets()); solbv=0.0;
-    mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets()); resbv=0.0;
-    mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets()); adjbv=0.0;
-    mfem::BlockVector prmbv; prmbv.Update(nf->PrmGetBlockTrueOffsets()); prmbv=1.0;
-    mfem::BlockVector grdbv; grdbv.Update(nf->PrmGetBlockTrueOffsets()); grdbv=0.0;
-
-    bgf00->SetFromTrueDofs(solbv.GetBlock(0));
-    pgf00->SetFromTrueDofs(prmbv.GetBlock(0));
-    pgf01->SetFromTrueDofs(prmbv.GetBlock(1));
-
-    nf->SetPrmFields(prmbv);
-    double energy = nf->GetEnergy(solbv);
-
-    nf->Mult(solbv,resbv);
-    std::cout<<"Norm res="<<resbv.Norml2()<<std::endl;
-
-    //mfem::Operator& K=nf->GetGradient(solbv);
-    std::cout<<"energy ="<< energy<<std::endl;
-
-    nf->SetStateFields(solbv);
-    nf->SetAdjointFields(adjbv);
-    nf->PrmMult(prmbv,grdbv);
-
-
-    //set the BC for the physics
-    mfem::Array<mfem::Array<int> *> ess_bdr;
-    mfem::Array<mfem::Vector*>      ess_rhs;
-    ess_bdr.Append(new mfem::Array<int>(mesh->bdr_attributes.Max()));
-    ess_rhs.Append(nullptr);
-    (*ess_bdr[0]) = 1;
-    nf->SetEssentialBC(ess_bdr,ess_rhs);
-
-    //define the solvers
-    mfem::UMFPackSolver* umfsolv=new mfem::UMFPackSolver();
-
-    mfem::GMRESSolver *gmres;
-    gmres = new mfem::GMRESSolver();
-    gmres->SetAbsTol(newton_abs_tol/10);
-    gmres->SetRelTol(newton_rel_tol/10);
-    gmres->SetMaxIter(100);
-    gmres->SetPrintLevel(print_level);
-    //gmres->SetPreconditioner(*prec);
-
-
-    mfem::NewtonSolver *ns;
-    ns = new mfem::NewtonSolver();
-    ns->iterative_mode = true;
-    ns->SetSolver(*gmres);
-    ns->SetOperator(*nf);
-    ns->SetPrintLevel(print_level);
-    ns->SetRelTol(newton_rel_tol);
-    ns->SetAbsTol(newton_abs_tol);
-    ns->SetMaxIter(newton_iter);
-
-    mfem::Vector b; //RHS is zero
-    solbv=0.0;
-    ns->Mult(b, solbv);
-
-
-    nf->SetStateFields(solbv);
-    nf->SetAdjointFields(solbv);
-    nf->PrmMult(prmbv,grdbv);
-
-
-
-
-    mfem::ParaViewDataCollection *dacol = new mfem::ParaViewDataCollection("Example91",
-                                                                      mesh);
-
-    ggf00->SetFromTrueDofs(grdbv.GetBlock(0));
-    ggf01->SetFromTrueDofs(grdbv.GetBlock(1));
-    pgf00->SetFromTrueDofs(solbv.GetBlock(0));
-
-    dacol->SetLevelsOfDetail(order);
-    dacol->RegisterField("sol", pgf00);
-    dacol->RegisterField("grad00", ggf00);
-    dacol->RegisterField("grad01", ggf01);
-
-    dacol->SetTime(1.0);
-    dacol->SetCycle(1);
-    dacol->Save();
-
-    delete dacol;
-
-    delete ns;
-    delete umfsolv;
-    delete gmres;
-    delete ess_bdr[0];
-
-    delete bgf00;
-    delete pgf00;
-    delete pgf01;
-    delete ggf00;
-    delete ggf01;
-
-
-    delete nf;
-    delete pfes01;
-    delete pfes00;
-    delete bfes00;
-
-    delete qfun;
-    delete lc;
-    delete dc;
-
-    delete mesh;
-
-}

examples/lshape_amr.mesh

@@ -0,0 +1,640 @@
+MFEM mesh v1.1
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+84
+1 3 0 21 53 24
+1 3 21 8 22 53
+1 3 53 22 18 23
+1 3 24 53 23 11
+1 3 8 25 54 22
+1 3 25 1 26 54
+1 3 54 26 9 27
+1 3 22 54 27 18
+1 3 18 65 92 68
+1 3 68 92 67 30
+1 3 92 66 55 67
+1 3 65 27 66 92
+1 3 27 9 28 55
+1 3 55 69 93 72
+1 3 69 28 70 93
+1 3 93 70 4 71
+1 3 72 93 71 29
+1 3 95 72 29 73
+1 3 74 95 73 10
+1 3 30 67 95 74
+1 3 67 55 72 95
+1 3 11 23 56 32
+1 3 23 18 30 56
+1 3 56 30 10 31
+1 3 32 56 31 3
+1 3 3 31 57 35
+1 3 31 10 33 57
+1 3 57 33 19 34
+1 3 35 57 34 14
+1 3 10 73 96 75
+1 3 75 96 77 33
+1 3 96 76 58 77
+1 3 73 29 76 96
+1 3 29 71 97 76
+1 3 71 94 108 98
+1 3 94 110 114 109
+1 3 110 4 111 114
+1 3 114 111 99 112
+1 3 109 114 112 108
+1 3 108 99 78 100
+1 3 98 108 100 97
+1 3 97 78 36 79
+1 3 76 97 79 58
+1 3 58 79 101 82
+1 3 79 36 80 101
+1 3 101 80 12 81
+1 3 82 101 81 37
+1 3 33 58 37 19
+1 3 19 37 59 40
+1 3 37 12 38 59
+1 3 59 38 7 39
+1 3 40 59 39 13
+1 3 14 34 60 42
+1 3 34 19 40 60
+1 3 60 40 13 41
+1 3 42 60 41 6
+1 3 4 103 113 99
+1 3 103 83 104 113
+1 3 113 104 102 105
+1 3 99 113 105 78
+1 3 78 102 85 36
+1 3 102 84 61 85
+1 3 83 43 84 102
+1 3 43 86 106 84
+1 3 86 15 87 106
+1 3 106 87 44 88
+1 3 84 106 88 61
+1 3 61 44 20 45
+1 3 36 61 45 12
+1 3 15 89 107 87
+1 3 87 107 91 44
+1 3 107 90 62 91
+1 3 89 46 90 107
+1 3 46 5 47 62
+1 3 62 47 16 48
+1 3 44 62 48 20
+1 3 20 48 63 51
+1 3 48 16 49 63
+1 3 63 49 2 50
+1 3 51 63 50 17
+1 3 12 45 64 38
+1 3 45 20 51 64
+1 3 64 51 17 52
+1 3 38 64 52 7
+
+boundary
+37
+1 1 0 21
+1 1 24 0
+1 1 21 8
+1 1 11 24
+1 1 8 25
+1 1 25 1
+1 1 1 26
+1 1 26 9
+1 1 9 28
+1 1 28 70
+1 1 70 4
+1 1 32 11
+1 1 3 32
+1 1 35 3
+1 1 14 35
+1 1 7 39
+1 1 39 13
+1 1 42 14
+1 1 13 41
+1 1 41 6
+1 1 6 42
+1 1 4 103
+1 1 103 83
+1 1 83 43
+1 1 43 86
+1 1 86 15
+1 1 15 89
+1 1 89 46
+1 1 46 5
+1 1 5 47
+1 1 47 16
+1 1 16 49
+1 1 49 2
+1 1 2 50
+1 1 50 17
+1 1 17 52
+1 1 52 7
+
+vertex_parents
+50
+65 18 27
+66 27 55
+67 30 55
+68 18 30
+69 28 55
+70 4 28
+71 4 29
+72 29 55
+73 10 29
+74 10 30
+75 10 33
+76 29 58
+77 33 58
+78 4 36
+79 36 58
+80 12 36
+81 12 37
+82 37 58
+83 4 43
+84 43 61
+85 36 61
+86 15 43
+87 15 44
+88 44 61
+89 15 46
+90 46 62
+91 44 62
+92 65 67
+93 69 71
+94 4 71
+95 67 73
+96 73 77
+97 71 79
+98 71 97
+99 4 78
+100 78 97
+101 79 81
+102 83 85
+103 4 83
+104 83 102
+105 78 102
+106 86 88
+107 89 91
+108 94 100
+109 94 108
+110 4 94
+111 4 99
+112 99 108
+113 103 105
+114 110 112
+
+coarse_elements
+12
+3 8 11 10 9
+3 13 14 15 16
+3 19 20 17 18
+3 29 32 31 30
+3 35 36 37 38
+3 34 88 39 40
+3 33 89 41 42
+3 43 44 45 46
+3 56 57 58 59
+3 92 62 61 60
+3 63 64 65 66
+3 69 72 71 70
+
+vertices
+115
+
+nodes
+FiniteElementSpace
+FiniteElementCollection: H1_2D_P2
+VDim: 2
+Ordering: 1
+
+-1 -1
+0 -1
+1 1
+-1 0
+0 0
+1 0
+-1 1
+0 1
+-0.5 -1
+0 -0.5
+-0.5 0
+-1 -0.5
+0 0.5
+-0.5 1
+-1 0.5
+0.5 0
+1 0.5
+0.5 1
+-0.5 -0.5
+-0.5 0.5
+0.5 0.5
+-0.75 -1
+-0.5 -0.75
+-0.75 -0.5
+-1 -0.75
+-0.25 -1
+0 -0.75
+-0.25 -0.5
+0 -0.25
+-0.25 0
+-0.5 -0.25
+-0.75 0
+-1 -0.25
+-0.5 0.25
+-0.75 0.5
+-1 0.25
+0 0.25
+-0.25 0.5
+0 0.75
+-0.25 1
+-0.5 0.75
+-0.75 1
+-1 0.75
+0.25 0
+0.5 0.25
+0.25 0.5
+0.75 0
+1 0.25
+0.75 0.5
+1 0.75
+0.75 1
+0.5 0.75
+0.25 1
+-0.75 -0.75
+-0.25 -0.75
+-0.25 -0.25
+-0.75 -0.25
+-0.75 0.25
+-0.25 0.25
+-0.25 0.75
+-0.75 0.75
+0.25 0.25
+0.75 0.25
+0.75 0.75
+0.25 0.75
+-0.375 -0.5
+-0.25 -0.375
+-0.375 -0.25
+-0.5 -0.375
+-0.125 -0.25
+0 -0.125
+-0.125 0
+-0.25 -0.125
+-0.375 0
+-0.5 -0.125
+-0.5 0.125
+-0.25 0.125
+-0.375 0.25
+0 0.125
+-0.125 0.25
+0 0.375
+-0.125 0.5
+-0.25 0.375
+0.125 0
+0.25 0.125
+0.125 0.25
+0.375 0
+0.5 0.125
+0.375 0.25
+0.625 0
+0.75 0.125
+0.625 0.25
+-0.375 -0.375
+-0.125 -0.125
+-0.0625 0
+-0.375 -0.125
+-0.375 0.125
+-0.125 0.125
+-0.125 0.0625
+0 0.0625
+-0.0625 0.125
+-0.125 0.375
+0.125 0.125
+0.0625 0
+0.125 0.0625
+0.0625 0.125
+0.375 0.125
+0.625 0.125
+-0.0625 0.0625
+-0.0625 0.03125
+-0.03125 0
+0 0.03125
+-0.03125 0.0625
+0.0625 0.0625
+-0.03125 0.03125
+-0.875 -1
+-0.75 -0.875
+-0.875 -0.75
+-1 -0.875
+-0.625 -1
+-0.5 -0.875
+-0.625 -0.75
+-0.5 -0.625
+-0.625 -0.5
+-0.75 -0.625
+-0.875 -0.5
+-1 -0.625
+-0.375 -1
+-0.25 -0.875
+-0.375 -0.75
+-0.125 -1
+0 -0.875
+-0.125 -0.75
+0 -0.625
+-0.125 -0.5
+-0.25 -0.625
+-0.375 -0.5
+-0.4375 -0.5
+-0.375 -0.4375
+-0.4375 -0.375
+-0.5 -0.4375
+-0.375 -0.3125
+-0.4375 -0.25
+-0.5 -0.3125
+-0.3125 -0.375
+-0.25 -0.3125
+-0.3125 -0.25
+-0.3125 -0.5
+-0.25 -0.4375
+0 -0.375
+-0.125 -0.25
+-0.25 -0.375
+-0.1875 -0.25
+-0.125 -0.1875
+-0.1875 -0.125
+-0.25 -0.1875
+-0.0625 -0.25
+0 -0.1875
+-0.0625 -0.125
+0 -0.0625
+-0.0625 0
+-0.125 -0.0625
+-0.1875 0
+-0.25 -0.0625
+-0.3125 -0.125
+-0.3125 0
+-0.375 -0.0625
+-0.4375 -0.125
+-0.4375 0
+-0.5 -0.0625
+-0.375 -0.1875
+-0.5 -0.1875
+-0.75 -0.375
+-0.875 -0.25
+-1 -0.375
+-0.5 -0.375
+-0.625 -0.25
+-0.5 -0.125
+-0.625 0
+-0.75 -0.125
+-0.875 0
+-1 -0.125
+-0.75 0.125
+-0.875 0.25
+-1 0.125
+-0.5 0.125
+-0.625 0.25
+-0.5 0.375
+-0.625 0.5
+-0.75 0.375
+-0.875 0.5
+-1 0.375
+-0.375 0.0625
+-0.4375 0.125
+-0.5 0.0625
+-0.375 0.1875
+-0.4375 0.25
+-0.5 0.1875
+-0.3125 0.125
+-0.25 0.1875
+-0.3125 0.25
+-0.25 0.0625
+-0.125 0.0625
+-0.1875 0.125
+-0.09375 0
+-0.0625 0.03125
+-0.09375 0.0625
+-0.125 0.03125
+-0.046875 0
+-0.03125 0.015625
+-0.046875 0.03125
+-0.0625 0.015625
+-0.015625 0
+0 0.015625
+-0.015625 0.03125
+0 0.046875
+-0.015625 0.0625
+-0.03125 0.046875
+-0.046875 0.0625
+-0.0625 0.046875
+-0.03125 0.0625
+0 0.09375
+-0.03125 0.125
+-0.0625 0.09375
+-0.09375 0.125
+-0.125 0.09375
+-0.0625 0.125
+0 0.1875
+-0.0625 0.25
+-0.125 0.1875
+-0.1875 0.25
+-0.125 0.3125
+-0.1875 0.375
+-0.25 0.3125
+0 0.3125
+-0.0625 0.375
+0 0.4375
+-0.0625 0.5
+-0.125 0.4375
+-0.1875 0.5
+-0.25 0.4375
+-0.375 0.25
+-0.25 0.375
+-0.375 0.5
+-0.25 0.625
+-0.375 0.75
+-0.5 0.625
+-0.125 0.5
+0 0.625
+-0.125 0.75
+0 0.875
+-0.125 1
+-0.25 0.875
+-0.375 1
+-0.5 0.875
+-0.75 0.625
+-0.875 0.75
+-1 0.625
+-0.625 0.75
+-0.625 1
+-0.75 0.875
+-0.875 1
+-1 0.875
+0.03125 0
+0.0625 0.03125
+0.03125 0.0625
+0 0.03125
+0.09375 0
+0.125 0.03125
+0.09375 0.0625
+0.125 0.09375
+0.09375 0.125
+0.0625 0.09375
+0.03125 0.125
+0.0625 0.125
+0.125 0.1875
+0.0625 0.25
+0.1875 0.125
+0.25 0.1875
+0.1875 0.25
+0.1875 0
+0.25 0.0625
+0.125 0.0625
+0.3125 0
+0.375 0.0625
+0.3125 0.125
+0.4375 0
+0.5 0.0625
+0.4375 0.125
+0.5 0.1875
+0.4375 0.25
+0.375 0.1875
+0.3125 0.25
+0.375 0.25
+0.5 0.375
+0.375 0.5
+0.25 0.375
+0.125 0.25
+0.125 0.5
+0 0.375
+0.5625 0
+0.625 0.0625
+0.5625 0.125
+0.625 0.1875
+0.5625 0.25
+0.6875 0.125
+0.75 0.1875
+0.6875 0.25
+0.6875 0
+0.75 0.0625
+0.875 0
+1 0.125
+0.875 0.25
+0.75 0.125
+1 0.375
+0.875 0.5
+0.75 0.375
+0.625 0.25
+0.625 0.5
+0.75 0.625
+0.625 0.75
+0.5 0.625
+1 0.625
+0.875 0.75
+1 0.875
+0.875 1
+0.75 0.875
+0.625 1
+0.5 0.875
+0.25 0.625
+0.125 0.75
+0.375 0.75
+0.375 1
+0.25 0.875
+0.125 1
+-0.875 -0.875
+-0.625 -0.875
+-0.625 -0.625
+-0.875 -0.625
+-0.375 -0.875
+-0.125 -0.875
+-0.125 -0.625
+-0.375 -0.625
+-0.4375 -0.4375
+-0.4375 -0.3125
+-0.3125 -0.3125
+-0.3125 -0.4375
+-0.125 -0.375
+-0.1875 -0.1875
+-0.0625 -0.1875
+-0.0625 -0.0625
+-0.1875 -0.0625
+-0.3125 -0.0625
+-0.4375 -0.0625
+-0.4375 -0.1875
+-0.3125 -0.1875
+-0.875 -0.375
+-0.625 -0.375
+-0.625 -0.125
+-0.875 -0.125
+-0.875 0.125
+-0.625 0.125
+-0.625 0.375
+-0.875 0.375
+-0.4375 0.0625
+-0.4375 0.1875
+-0.3125 0.1875
+-0.3125 0.0625
+-0.1875 0.0625
+-0.09375 0.03125
+-0.046875 0.015625
+-0.015625 0.015625
+-0.015625 0.046875
+-0.046875 0.046875
+-0.03125 0.09375
+-0.09375 0.09375
+-0.0625 0.1875
+-0.1875 0.1875
+-0.1875 0.3125
+-0.0625 0.3125
+-0.0625 0.4375
+-0.1875 0.4375
+-0.375 0.375
+-0.375 0.625
+-0.125 0.625
+-0.125 0.875
+-0.375 0.875
+-0.875 0.625
+-0.625 0.625
+-0.625 0.875
+-0.875 0.875
+0.03125 0.03125
+0.09375 0.03125
+0.09375 0.09375
+0.03125 0.09375
+0.0625 0.1875
+0.1875 0.1875
+0.1875 0.0625
+0.3125 0.0625
+0.4375 0.0625
+0.4375 0.1875
+0.3125 0.1875
+0.375 0.375
+0.125 0.375
+0.5625 0.0625
+0.5625 0.1875
+0.6875 0.1875
+0.6875 0.0625
+0.875 0.125
+0.875 0.375
+0.625 0.375
+0.625 0.625
+0.875 0.625
+0.875 0.875
+0.625 0.875
+0.125 0.625
+0.375 0.625
+0.375 0.875
+0.125 0.875

examples/mal_rllib_eval.py

@@ -0,0 +1,188 @@
+
+import commentjson
+from ray.rllib.agents.registry import get_agent_class
+import amr_env
+import gym
+from gym import spaces
+import ray
+import ray.rllib.agents.ppo as ppo
+import tensorflow as tf
+
+import numpy as np
+
+
+import argparse
+import random
+from amr import models
+import pytest
+import os
+import copy
+from math import sqrt, nan, inf, isnan
+from ray.rllib.models import ModelCatalog
+from amr.models.cnn import CNNSmall
+
+# rllib requires you to give the policy an env with the same action
+# and observation spaces as used in training. The rest of it can be
+# "fake" if you provide your own observation data some other way.
+
+# USER INPUT - solution, 20x20 mesh, sine,tanh,steps,steps2, norm-diff reward with random threshold
+#866764
+#checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_5b99d_00000_0_2021-06-07_12-12-52/"
+#checkpoint_number = 900
+
+# USER INPUT - solution, 20x20 mesh, sine,tanh,steps,steps2, binary reward with random threshold
+#866762
+# checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_8d234_00000_0_2021-06-07_12-07-06/"
+# checkpoint_number = 900
+
+# solution, 20x20, steps2,sine,tanh,bumps,  binary with random
+# 880258 - fixed threshold 1.e-5
+checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_cc8cb_00000_0_2021-06-16_13-38-25"
+checkpoint_number = 1400
+
+# 880259 - random threshold [1.e-2, 1.e-6]
+checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_e1561_00000_0_2021-06-16_13-31-50"
+checkpoint_number = 900
+# with 1.e-3 - the second and third refinements are really good. still more than needed in first
+# with 1.e-2 - picks the right amount of elements.
+
+#880328 - fixed threshold 1.e-2
+#checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_3a8a0_00000_0_2021-06-16_16-18-58"
+#checkpoint_number = 300
+#refines 15 elements at first iteration.. not good with 300.
+
+# solution, 10x10, steps2,sine,tanh,bumps,  binary with random [1.e-2, 1.e-6]
+# slurm-880260.out
+#checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_cc8cb_00000_0_2021-06-16_13-38-25/"
+#checkpoint_number = 1500
+#1.e-3 -> right region but first iteration has too many
+#1.e-2 is almost perfect
+
+
+
+
+#this has the policy without error threshold
+#checkpoint_folder = "/p/lustre1/mittal3/local_deref/PPO/PPO_LocalAMR-v0_c6ca5_00000_0_2021-05-24_16-07-13/"
+
+# END OF USER INPUT
+
+# Read info from json file used for training.
+full_checkpoint_path = checkpoint_folder + '/checkpoint_' + str(checkpoint_number)+ '/checkpoint-' + str(checkpoint_number)
+path_env_config_file = checkpoint_folder + '/params.json'
+with open(path_env_config_file) as json_file:
+    env_trainer_config = commentjson.load(json_file)
+
+trainer_config = env_trainer_config
+trainer_config['env_config']['mesh_params']['nx'] = 1
+trainer_config['env_config']['mesh_params']['ny'] = 1
+local_sample = env_trainer_config['env_config']['local_sample']
+local_context = env_trainer_config['env_config']['local_context']
+reward_params = env_trainer_config['env_config']['reward_function_params']
+
+#set some default params
+observe_error = False
+observe_values = True
+observe_grads = False
+
+#get observing quantities
+observe_values = env_trainer_config['env_config']['observe_values']
+observe_depth = env_trainer_config['env_config']['observe_depth']
+observe_jacobian = env_trainer_config['env_config']['observe_jacobian']
+observe_ar = env_trainer_config['env_config']['observe_ar']
+observe_grads = env_trainer_config['env_config']['observe_grads']
+normalization = env_trainer_config['env_config']['normalization']
+
+#get reward_params
+reward_params = env_trainer_config['env_config']['reward_function_params']
+reward_params_name = reward_params['name']
+if reward_params_name == "random_penalized_norm_diff":
+    observe_error = True
+if reward_params_name == "random_binary":
+    observe_error = True
+
+class DummyEnv(gym.Env):
+
+    def __init__(self, config): # the config param is required by rllib
+
+        # image size is 42x42 (a size which uses CNN by default in rllib)
+        self.obsx = local_sample+2*local_context
+        self.obsy = self.obsx
+
+        # Either do nothing (0) or refine (1)
+        self.action_space = spaces.Discrete(2)
+
+        n_channels = 1
+        n_channels = observe_values + observe_depth + observe_grads
+
+        low = -np.inf
+        
+        high = np.inf
+        self.observation_space = spaces.Dict({
+                "scalar_info": spaces.Box(low=low, high=high, shape=(1 + observe_jacobian + observe_error, ), dtype=np.float32),
+                "obs_data":    spaces.Box(low=low, high=high,
+                shape=(self.obsx,
+                       self.obsy,
+                       n_channels), dtype=np.float32)
+                })
+
+        self.state = None
+        
+    def step(self, action):
+        pass
+    def reset(self):
+        pass
+    def render(self):
+        pass
+
+class Evaluator():
+
+    def __init__(self):
+
+        print("starting ray...")
+        ray.shutdown()
+        ray.init()
+        print("ray up...")
+
+        ModelCatalog.register_custom_model("cnn_small", CNNSmall)
+        trainer = ppo.PPOTrainer(config=trainer_config,env=DummyEnv)
+        trainer.restore(full_checkpoint_path)
+
+        trainer_config['evaluation_num_workers'] = 1
+        trainer_config['evaluation_interval'] = 0
+        trainer_config['num_workers'] = 0
+        trainer_config['num_envs_per_worker'] = 1
+
+        self.agent = ppo.PPOTrainer(config=trainer_config,env=DummyEnv)
+        self.agent.restore(full_checkpoint_path)
+        
+        self.env = DummyEnv({})
+
+    def eval(self,obso,scalar):
+        if observe_values and normalization:
+            obso -= np.mean(obso)
+        
+        obs = {
+            "obs_data" : obso,
+            "scalar_info" : scalar
+            }
+
+        pick = self.agent.compute_action(obs, explore=False)
+        return pick
+
+    def get_local_sample(self):
+        return local_sample
+
+    def get_local_context(self):
+        return local_context
+
+    def get_observe_error(self):
+        return observe_error
+
+    def get_observe_jacobian(self):
+        return observe_jacobian
+
+    def get_observe_values(self):
+        return observe_values
+
+    def get_observe_gradient(self):
+        return observe_grads

examples/multi_agent_local_refiner.hpp

@@ -0,0 +1,419 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace mfem;
+using namespace std;
+
+class MAL_DRLRefiner : public MeshOperator
+{
+protected:
+   GridFunction& u;
+
+   Array<Refinement> marked_elements;
+   long current_sequence;
+
+   int imgsz;
+   int local_sample;
+   int local_context;
+   bool observe_jacobian;
+   bool observe_error;
+   bool observe_values;
+   bool observe_gradient;
+
+   int nc_limit;
+
+   PyObject* eval_method;
+   PyObject* get_local_sample_method;
+   PyObject* get_local_context_method;
+   PyObject* get_observe_error_method;
+   PyObject* get_observe_jacobian_method;
+   PyObject* get_observe_values_method;
+   PyObject* get_observe_gradient_method;
+
+   FindPointsGSLIB *gslib;
+
+   /** @brief Apply the operator to the mesh.
+       @return STOP if a stopping criterion is satisfied or no elements were
+       marked for refinement; REFINED + CONTINUE otherwise. */
+   virtual int ApplyImpl(Mesh &mesh);
+
+public:
+
+   /// Construct a MAL_DRLRefiner that will operate on u.
+   MAL_DRLRefiner(GridFunction &u);
+
+   // default destructor (virtual)
+
+   /** @brief Set the maximum ratio of refinement levels of adjacent elements
+       (0 = unlimited). */
+   void SetNCLimit(int nc_limit)
+   {
+      MFEM_ASSERT(nc_limit >= 0, "Invalid NC limit");
+      this->nc_limit = nc_limit;
+   }
+
+   virtual void Reset();
+};
+
+MAL_DRLRefiner::MAL_DRLRefiner(GridFunction& u_) : u(u_)
+{
+   int ret = _import_array();
+   if (ret < 0) {
+      printf("problem with import_array\n");
+   }
+
+   PyRun_SimpleString("import sys");
+   PyRun_SimpleString("sys.path.append('.')");
+
+   // This is a workaround for something in tensorflow that dies without it.
+   PyRun_SimpleString("if not hasattr(sys, 'argv'):\n"
+                      "  sys.argv  = ['']");
+
+   PyObject* eval_mod = PyImport_ImportModule("mal_rllib_eval");
+   if (eval_mod == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   PyObject* eval_class = PyObject_GetAttrString(eval_mod, "Evaluator");
+   if (eval_class == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+   Py_DECREF(eval_mod);
+
+   PyObject* args = Py_BuildValue("()");
+   if (args == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   PyObject* eval_obj = PyEval_CallObject(eval_class, args);
+   if (eval_obj == NULL) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   eval_method = PyObject_GetAttrString(eval_obj, "eval");
+   if (eval_method == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   get_local_sample_method = PyObject_GetAttrString(eval_obj, "get_local_sample");
+   get_local_context_method = PyObject_GetAttrString(eval_obj, "get_local_context");
+   get_observe_error_method = PyObject_GetAttrString(eval_obj, "get_observe_error");
+   get_observe_jacobian_method = PyObject_GetAttrString(eval_obj, "get_observe_jacobian");
+   get_observe_values_method = PyObject_GetAttrString(eval_obj, "get_observe_values");
+   get_observe_gradient_method = PyObject_GetAttrString(eval_obj, "get_observe_gradient");
+
+   PyObject* local_sample_p = PyObject_CallFunctionObjArgs(
+                              get_local_sample_method, nullptr);
+   PyObject* local_context_p = PyObject_CallFunctionObjArgs(
+                              get_local_context_method, nullptr);
+   PyObject* observe_jacobian_p = PyObject_CallFunctionObjArgs(
+                              get_observe_jacobian_method, nullptr);
+   PyObject* observe_error_p = PyObject_CallFunctionObjArgs(
+                              get_observe_error_method, nullptr);
+   PyObject* observe_values_p = PyObject_CallFunctionObjArgs(
+                                get_observe_values_method, nullptr);
+   PyObject* observe_gradient_p = PyObject_CallFunctionObjArgs(
+                                  get_observe_gradient_method, nullptr);
+
+   PyArg_Parse(local_sample_p, "i", &local_sample);
+   PyArg_Parse(local_context_p, "i", &local_context);
+   int observe_jacobian_i, observe_error_i, observe_values_i, observe_gradient_i;
+   PyArg_Parse(observe_jacobian_p, "i", &observe_jacobian_i);
+   PyArg_Parse(observe_error_p, "i", &observe_error_i);
+   PyArg_Parse(observe_values_p, "i", &observe_values_i);
+   PyArg_Parse(observe_gradient_p, "i", &observe_gradient_i);
+
+   observe_jacobian = bool(observe_jacobian_i);
+   observe_error = bool(observe_error_i);
+   observe_values = bool(observe_values_i);
+   observe_gradient = bool(observe_gradient_i);
+
+#ifdef MFEM_USE_GSLIB
+   // setup gslib
+   gslib = new FindPointsGSLIB();
+   std::cout << " initialize findpts\n";
+#endif
+}
+
+int MAL_DRLRefiner::ApplyImpl(Mesh &mesh)
+{
+#ifdef MFEM_USE_GSLIB
+   // setup gslib
+   gslib->FreeData();
+   gslib->Setup(mesh);
+#endif
+   double u_min = u.Min();
+   double u_max = u.Max();
+   // u -= u_min;
+   // u /= (u_max-u_min); 
+
+   marked_elements.SetSize(0);
+   imgsz = local_sample + 2 * local_context;
+
+   GridFunction ugrad(u.FESpace());
+   GridFunction ugradmag(u.FESpace());
+
+   if (observe_gradient) {
+      const int s = ugrad.Size();
+
+      for (int d = 0; d < 2; d++)
+      {
+         u.GetDerivative(1, d, ugrad);
+         for (int i = 0; i < s; i++)
+         {
+            ugradmag(i) += pow(ugrad(i), 2.0);
+         }
+      }
+      for (int i = 0; i < s; i++)
+      {
+         ugradmag(i) = sqrt(ugradmag(i) + 1e-12);
+      }
+   }
+      
+   for (int k = 0; k < mesh.GetNE(); k++) {
+
+      // get scalar info
+      const int scalar_size = 1 + observe_error + observe_jacobian;
+      double* scalar_obs1 = new double[scalar_size];
+      scalar_obs1[0] = 1;
+      bool boundary = false;
+      Array<int> fcs, cor;
+      int e1, e2, inf1, inf2, ncf;
+      mesh.GetElementEdges(k, fcs, cor);
+      for (int f = 0; f < fcs.Size() && boundary == false; f++) {
+         mesh.GetFaceElements(fcs[f], &e1, &e2);
+         mesh.GetFaceInfos(fcs[f], &inf1, &inf2, &ncf);
+         if (e2 < 0 && inf2 < 0 && ncf == -1) {
+            boundary = true;
+         }
+      }
+      if (!boundary) {
+         scalar_obs1[0] = 2;
+      }
+
+
+      double error_threshold = 1.0e-2;
+      if (observe_jacobian && observe_error) {
+         scalar_obs1[1] = mesh.GetElementVolume(k);
+         scalar_obs1[2] = error_threshold;
+      }
+      else if (observe_jacobian) {
+         scalar_obs1[1] = mesh.GetElementVolume(k);
+      }
+      else if (observe_error) {
+         scalar_obs1[1] = error_threshold;
+      }
+
+      // assemble matrix of sample points, ref_space -> phys_space
+      ElementTransformation* trk = mesh.GetElementTransformation(k);
+      IntegrationPoint ipk;
+      Vector xk(2);
+      DenseMatrix m(2,imgsz*imgsz);
+      int c = 0;
+      double r_init = 0.001,
+             r_final = 1.0-r_init;
+      for (int j = 0; j < imgsz; ++j) {
+         for (int i = 0; i < imgsz; ++i) {
+            if (!boundary) {
+               ipk.y = (j - local_context + 0.5)/local_sample;
+               ipk.x = (i - local_context + 0.5)/local_sample;
+            }
+            else {
+               ipk.y = (r_init + j*(r_final-r_init))/(imgsz-1);
+               ipk.x = (r_init + i*(r_final-r_init))/(imgsz-1);
+            }
+
+            trk->Transform(ipk, xk);
+            m.SetCol(c++,xk);
+         }
+      } 
+
+      // phys_space -> elements, ips
+      Array<int> elems(imgsz*imgsz);
+      Array<IntegrationPoint> ips(imgsz*imgsz);
+      int n;
+      bool complete;
+      Vector ui;
+      double* obs1 = new double[imgsz*imgsz];
+#ifdef MFEM_USE_GSLIB
+      DenseMatrix mt(imgsz*imgsz, 2);
+      mt.Transpose(m);
+      Vector xyz(mt.GetData(), imgsz*imgsz*2);
+      ui.SetDataAndSize(obs1, imgsz*imgsz);
+      gslib->Interpolate(xyz, u, ui);
+#else
+      n = mesh.FindPoints(m, elems, ips, false);
+
+      // Build observation from GridFunction using elements, ips
+      n = 0;
+      complete = true;
+      for (int j = 0; j < imgsz; ++j) {
+         for (int i = 0; i < imgsz; ++i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs1[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               if (observe_gradient) {
+                  //std::cout << i << " " << j << " k10getgradmag\n";
+                  obs1[i*imgsz+j] = ugradmag.GetValue(el, ip);
+                  obs1[i*imgsz+j] = u.GetValue(el, ip);
+               }
+               else {
+                  obs1[i*imgsz+j] = u.GetValue(el, ip);
+               }
+            }
+            n++;
+         }
+      }
+#endif
+
+      // TODO: More efficient way to do the below mirroring
+      // invert i
+      double* obs2 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = 0; j < imgsz; ++j) {
+         for (int i = imgsz-1; i >= 0; --i) {
+            // int el = elems[n];
+            // if (el == -1) {
+            //    obs2[i*imgsz+j] = 0.0;
+            //    complete = false;
+            // }
+            // else {
+            //    IntegrationPoint& ip = ips[n];
+            //    obs2[i*imgsz+j] = u.GetValue(el, ip);
+            // }
+            obs2[i*imgsz+j] = obs1[n];
+            n++;
+         }
+      }
+
+      // invert j
+      double* obs3 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = imgsz-1; j >= 0; --j) {
+         for (int i = 0; i < imgsz; ++i) {
+            // int el = elems[n];
+            // if (el == -1) {
+            //    obs3[i*imgsz+j] = 0.0;
+            //    complete = false;
+            // }
+            // else {
+            //    IntegrationPoint& ip = ips[n];
+            //    obs3[i*imgsz+j] = u.GetValue(el, ip);
+            // }
+            obs3[i*imgsz+j] = obs1[n];
+            n++;
+         }
+      }
+
+      // invert i and j
+      double* obs4 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = imgsz-1; j >= 0; --j) {
+         for (int i = imgsz-1; i >= 0; --i) {
+            // int el = elems[n];
+            // if (el == -1) {
+            //    obs4[i*imgsz+j] = 0.0;
+            //    complete = false;
+            // }
+            // else {
+            //    IntegrationPoint& ip = ips[n];
+            //    obs4[i*imgsz+j] = u.GetValue(el, ip);
+            // }
+            obs4[i*imgsz+j] = obs1[n];
+            n++;
+         }
+      }
+
+
+      // apply policy: state -> action
+      bool refine = false;
+      if (complete) {
+
+         // convert to numpy array
+         npy_intp dims[3];
+         dims[0] = imgsz;
+         dims[1] = imgsz;
+         dims[2] = 1;
+
+         npy_intp scalar_dims[1];
+         scalar_dims[0] = scalar_size;
+
+         PyObject *pArray1 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs1));
+         PyObject *pArray2 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs2));
+         PyObject *pArray3 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs3));
+         PyObject *pArray4 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs4));   
+         PyObject *sArray1 = PyArray_SimpleNewFromData(
+            1, scalar_dims, NPY_DOUBLE, reinterpret_cast<void*>(scalar_obs1));
+         if (pArray1 == NULL) printf("pArray1 NULL!\n");
+
+         PyObject* action1 = PyObject_CallFunctionObjArgs(
+                     eval_method, pArray1, sArray1, nullptr);
+         PyObject* action2 = PyObject_CallFunctionObjArgs(
+                     eval_method, pArray2, sArray1, nullptr);
+         PyObject* action3 = PyObject_CallFunctionObjArgs(
+                     eval_method, pArray3, sArray1, nullptr);
+         PyObject* action4 = PyObject_CallFunctionObjArgs(
+                     eval_method, pArray4, sArray1, nullptr);
+
+         if (action1 == 0 || action2 == 0) {
+            PyErr_Print();
+            exit(1);
+         }
+
+         // parse integer return value
+         int action1_val;
+         PyArg_Parse(action1, "i", &action1_val);
+         int action2_val;
+         PyArg_Parse(action2, "i", &action2_val);
+         int action3_val;
+         PyArg_Parse(action3, "i", &action3_val);
+         int action4_val;
+         PyArg_Parse(action4, "i", &action4_val);
+         refine =
+            bool(action1_val) ||
+            bool(action2_val) ||
+            bool(action3_val) ||
+            bool(action4_val);
+      }
+
+      delete scalar_obs1;
+      delete obs1;
+      delete obs2;
+      delete obs3;
+      delete obs4;
+
+      if (refine) {
+         marked_elements.Append(Refinement(k));
+      }
+   }
+
+   long int num_marked_elements = mesh.ReduceInt(marked_elements.Size());
+   printf("marked %d elements\n",num_marked_elements);
+   if (num_marked_elements == 0) { return STOP; }
+
+   bool nonconforming = true;
+   mesh.GeneralRefinement(marked_elements, nonconforming, nc_limit);
+   return CONTINUE + REFINED;
+}
+
+void MAL_DRLRefiner::Reset()
+{
+   current_sequence = -1;
+}

examples/rllib_eval.py

@@ -0,0 +1,64 @@
+
+
+import gym
+from gym import spaces
+import ray
+import ray.rllib.agents.ppo as ppo
+import tensorflow as tf
+
+import numpy as np
+
+# rllib requires you to give the policy an env with the same action
+# and observation spaces as used in training. The rest of it can be
+# "fake" if you provide your own observation data some other way.
+
+class DummyEnv(gym.Env):
+
+    def __init__(self, config): # the config param is required by rllib
+
+        # image size is 42x42 (a size which uses CNN by default in rllib)
+        self.obsx = 42
+        self.obsy = 42
+
+        # Either do nothing (0) or refine (1)
+        self.action_space = spaces.Discrete(2)
+        self.observation_space = spaces.Box(-1.0, 2.0, shape=(self.obsx,self.obsy,1))
+
+        self.state = None
+        
+    def step(self, action):
+        pass
+    def reset(self):
+        pass
+    def render(self):
+        pass
+
+class Evaluator():
+
+    def __init__(self):
+
+        print("starting ray...")
+        ray.shutdown()
+        ray.init()
+        print("ray up...")
+
+        config = ppo.DEFAULT_CONFIG.copy()
+        config["log_level"] = "WARN"
+
+        # Create agent from checkpoint
+        self.agent = ppo.PPOTrainer(config,env=DummyEnv)
+        self.agent.restore("DRLRefinePolicy/checkpoint_210/checkpoint-210")
+
+        self.env = DummyEnv({})
+
+    def eval(self,obs):
+        pick = self.agent.compute_action(obs, explore=False)
+        return pick
+
+#evaluator = Evaluator()
+#obs = np.ones((42,42,1))
+#ref = evaluator.eval(obs)
+#evaluator.show_logits(obs)
+# evaluator.eval(np.random.rand(8))
+# evaluator.eval(np.random.rand(8))
+

fem/CMakeLists.txt

@@ -66,7 +66,6 @@ set(SRCS
   tmop_tools.cpp
   gslib.cpp
   transfer.cpp
-  prmnonlinearform.cpp
   )
 
 set(HDRS
@@ -113,7 +112,6 @@ set(HDRS
   tmop_tools.hpp
   gslib.hpp
   transfer.hpp
-  prmnonlinearform.hpp
   )
 
 if (MFEM_USE_SIDRE)
@@ -138,8 +136,7 @@ if (MFEM_USE_MPI)
     pgridfunc.cpp
     plinearform.cpp
     pnonlinearform.cpp
-    prestriction.cpp
-    pprmnonlinearform.cpp)
+    prestriction.cpp)
   # If this list (HDRS -> HEADERS) is used for install, we probably want the
   # headers added all the time.
   list(APPEND HDRS
@@ -148,8 +145,7 @@ if (MFEM_USE_MPI)
     pgridfunc.hpp
     plinearform.hpp
     pnonlinearform.hpp
-    prestriction.hpp
-    pprmnonlinearform.hpp)
+    prestriction.hpp)
 endif()
 
 convert_filenames_to_full_paths(SRCS)

fem/bilininteg.cpp

@@ -2406,13 +2406,24 @@ void VectorDiffusionIntegrator::AssembleElementMatrix(
    const int dim = el.GetDim();
    const int dof = el.GetDof();
    const int sdim = Trans.GetSpaceDim();
-   const bool square = (dim == sdim);
-   double w;
 
-   elmat.SetSize(sdim * dof);
+   // If vdim is not set, set it to the space dimension;
+   vdim = (vdim <= 0) ? sdim : vdim;
+   const bool square = (dim == sdim);
+
+   if (VQ)
+   {
+      vcoeff.SetSize(vdim);
+   }
+   else if (MQ)
+   {
+      mcoeff.SetSize(vdim);
+   }
 
    dshape.SetSize(dof, dim);
    dshapedxt.SetSize(dof, sdim);
+
+   elmat.SetSize(vdim * dof);
    pelmat.SetSize(dof);
 
    const IntegrationRule *ir = IntRule;
@@ -2422,28 +2433,48 @@ void VectorDiffusionIntegrator::AssembleElementMatrix(
    }
 
    elmat = 0.0;
-   pelmat = 0.0;
 
    for (int i = 0; i < ir -> GetNPoints(); i++)
    {
+
       const IntegrationPoint &ip = ir->IntPoint(i);
-      el.CalcDShape (ip, dshape);
-      Trans.SetIntPoint (&ip);
-      w = Trans.Weight();
+      el.CalcDShape(ip, dshape);
+
+      Trans.SetIntPoint(&ip);
+      double w = Trans.Weight();
       w = ip.weight / (square ? w : w*w*w);
       // AdjugateJacobian = / adj(J),         if J is square
       //                    \ adj(J^t.J).J^t, otherwise
       Mult(dshape, Trans.AdjugateJacobian(), dshapedxt);
-      if (Q) {  w *= Q -> Eval (Trans, ip); }
-      AddMult_a_AAt(w, dshapedxt, pelmat);
-   }
-   for (int d = 0; d < sdim; d++)
-   {
-      for (int k = 0; k < dof; k++)
+
+      if (VQ)
       {
-         for (int l = 0; l < dof; l++)
+         VQ->Eval(vcoeff, Trans, ip);
+         for (int k = 0; k < vdim; ++k)
          {
-            elmat(dof*d+k, dof*d+l) = pelmat(k, l);
+            Mult_a_AAt(w*vcoeff(k), dshapedxt, pelmat);
+            elmat.AddMatrix(pelmat, dof*k, dof*k);
+         }
+      }
+      else if (MQ)
+      {
+         MQ->Eval(mcoeff, Trans, ip);
+         for (int i = 0; i < vdim; ++i)
+         {
+            for (int j = 0; j < vdim; ++j)
+            {
+               Mult_a_AAt(w*mcoeff(i,j), dshapedxt, pelmat);
+               elmat.AddMatrix(pelmat, dof*i, dof*j);
+            }
+         }
+      }
+      else
+      {
+         if (Q) { w *= Q->Eval(Trans, ip); }
+         Mult_a_AAt(w, dshapedxt, pelmat);
+         for (int k = 0; k < vdim; ++k)
+         {
+            elmat.AddMatrix(pelmat, dof*k, dof*k);
          }
       }
    }
@@ -2453,16 +2484,32 @@ void VectorDiffusionIntegrator::AssembleElementVector(
    const FiniteElement &el, ElementTransformation &Tr,
    const Vector &elfun, Vector &elvect)
 {
-   int dim = el.GetDim(); // assuming vector_dim == reference_dim
-   int dof = el.GetDof();
-   double w;
+   const int dim = el.GetDim();
+   const int dof = el.GetDof();
+   const int sdim = Tr.GetSpaceDim();
+
+   // If vdim is not set, set it to the space dimension;
+   vdim = (vdim <= 0) ? sdim : vdim;
+   const bool square = (dim == sdim);
+
+   if (VQ)
+   {
+      vcoeff.SetSize(vdim);
+   }
+   else if (MQ)
+   {
+      mcoeff.SetSize(vdim);
+   }
 
-   Jinv.SetSize(dim);
    dshape.SetSize(dof, dim);
-   pelmat.SetSize(dim);
-   gshape.SetSize(dim);
+   dshapedxt.SetSize(dof, dim);
+   // pelmat.SetSize(dim);
 
    elvect.SetSize(dim*dof);
+
+   // NOTE: DenseMatrix is in column-major order. This is consistent with
+   // vectors ordered byNODES. In the resulting DenseMatrix, each column
+   // corresponds to a particular vdim.
    DenseMatrix mat_in(elfun.GetData(), dof, dim);
    DenseMatrix mat_out(elvect.GetData(), dof, dim);
 
@@ -2476,22 +2523,50 @@ void VectorDiffusionIntegrator::AssembleElementVector(
    for (int i = 0; i < ir->GetNPoints(); i++)
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
-
-      Tr.SetIntPoint(&ip);
-      CalcAdjugate(Tr.Jacobian(), Jinv);
-      w = ip.weight / Tr.Weight();
-      if (Q)
-      {
-         w *= Q->Eval(Tr, ip);
-      }
-      MultAAt(Jinv, gshape);
-      gshape *= w;
-
       el.CalcDShape(ip, dshape);
 
-      MultAtB(mat_in, dshape, pelmat);
-      MultABt(pelmat, gshape, Jinv);
-      AddMultABt(dshape, Jinv, mat_out);
+      Tr.SetIntPoint(&ip);
+      double w = Tr.Weight();
+      w = ip.weight / (square ? w : w*w*w);
+      Mult(dshape, Tr.AdjugateJacobian(), dshapedxt);
+      MultAAt(dshapedxt, pelmat);
+
+      if (VQ)
+      {
+         VQ->Eval(vcoeff, Tr, ip);
+         for (int k = 0; k < vdim; ++k)
+         {
+            pelmat *= w*vcoeff(k);
+            const Vector vec_in(mat_in.GetColumn(k), dof);
+            Vector vec_out(mat_out.GetColumn(k), dof);
+            pelmat.AddMult(vec_in, vec_out);
+         }
+      }
+      else if (MQ)
+      {
+         MQ->Eval(mcoeff, Tr, ip);
+         for (int i = 0; i < vdim; ++i)
+         {
+            Vector vec_out(mat_out.GetColumn(i), dof);
+            for (int j = 0; j < vdim; ++j)
+            {
+               pelmat *= w*mcoeff(i,j);
+               const Vector vec_in(mat_in.GetColumn(j), dof);
+               pelmat.Mult(vec_in, vec_out);
+            }
+         }
+      }
+      else
+      {
+         if (Q) { w *= Q->Eval(Tr, ip); }
+         pelmat *= w;
+         for (int k = 0; k < vdim; ++k)
+         {
+            const Vector vec_in(mat_in.GetColumn(k), dof);
+            Vector vec_out(mat_out.GetColumn(k), dof);
+            pelmat.AddMult(vec_in, vec_out);
+         }
+      }
    }
 }
 

fem/bilininteg.hpp

@@ -2627,15 +2627,24 @@ public:
 };
 
 /** Integrator for
+
       (Q grad u, grad v) = sum_i (Q grad u_i, grad v_i) e_i e_i^T
-    for FE spaces defined by 'dim' copies of a scalar FE space. Where e_i
-    is the unit vector in the i-th direction.  The resulting local element
-    matrix is a block-diagonal matrix consisting of 'dim' copies of a scalar
-    diffusion matrix in each diagonal block. */
+
+    for vector FE spaces, where e_i is the unit vector in the i-th direction.
+    The resulting local element matrix is square, of size <tt> vdim*dof </tt>,
+    where \c vdim is the vector dimension space and \c dof is the local degrees
+    of freedom. The integrator is not aware of the true vector dimension and
+    must use \c VectorCoefficient, \c MatrixCoefficient, or a caller-specified
+    value to determine the vector space. For a scalar coefficient, the caller
+    may manually specify the vector dimension or the vector dimension is assumed
+    to be the spatial dimension (i.e. 2-dimension or 3-dimension).
+*/
 class VectorDiffusionIntegrator : public BilinearFormIntegrator
 {
 protected:
-   Coefficient *Q;
+   Coefficient *Q = NULL;
+   VectorCoefficient *VQ = NULL;
+   MatrixCoefficient *MQ = NULL;
 
    // PA extension
    const DofToQuad *maps;         ///< Not owned
@@ -2644,17 +2653,62 @@ protected:
    Vector pa_data;
 
    // CEED extension
-   CeedData* ceedDataPtr;
+   CeedData* ceedDataPtr = NULL;
 
 private:
    DenseMatrix dshape, dshapedxt, pelmat;
-   DenseMatrix Jinv, gshape;
+   int vdim = -1;
+   DenseMatrix mcoeff;
+   Vector vcoeff;
 
 public:
-   VectorDiffusionIntegrator()
-      : Q(NULL), ceedDataPtr(NULL) { }
+   VectorDiffusionIntegrator() { }
+
+   /** \brief Integrator with unit coefficient for caller-specified vector
+       dimension.
+
+       If the vector dimension does not match the true dimension of the space,
+       the resulting element matrix will be mathematically invalid. */
+   VectorDiffusionIntegrator(int vector_dimension)
+      : vdim(vector_dimension) { }
+
    VectorDiffusionIntegrator(Coefficient &q)
-      : Q(&q), ceedDataPtr(NULL) { }
+      : Q(&q) { }
+
+   /** \brief Integrator with scalar coefficient for caller-specified vector
+       dimension.
+
+       The element matrix is block-diagonal with \c vdim copies of the element
+       matrix integrated with the \c Coefficient.
+
+       If the vector dimension does not match the true dimension of the space,
+       the resulting element matrix will be mathematically invalid. */
+   VectorDiffusionIntegrator(Coefficient &q, int vector_dimension)
+      : Q(&q), vdim(vector_dimension) { }
+
+   /** \brief Integrator with \c VectorCoefficient. The vector dimension of the
+       \c FiniteElementSpace is assumed to be the same as the dimension of the
+       \c Vector.
+
+       The element matrix is block-diagonal and each block is integrated with
+       coefficient q_i.
+
+       If the vector dimension does not match the true dimension of the space,
+       the resulting element matrix will be mathematically invalid. */
+   VectorDiffusionIntegrator(VectorCoefficient &vq)
+      : VQ(&vq), vdim(vq.GetVDim()) { }
+
+   /** \brief Integrator with \c MatrixCoefficient. The vector dimension of the
+       \c FiniteElementSpace is assumed to be the same as the dimension of the
+       \c Matrix.
+
+       The element matrix is populated in each block. Each block is integrated
+       with coefficient q_ij.
+
+       If the vector dimension does not match the true dimension of the space,
+       the resulting element matrix will be mathematically invalid. */
+   VectorDiffusionIntegrator(MatrixCoefficient& mq)
+      : MQ(&mq), vdim(mq.GetVDim()) { }
 
    virtual ~VectorDiffusionIntegrator()
    {

fem/fem.hpp

@@ -41,7 +41,6 @@
 #include "transfer.hpp"
 #include "fespacehierarchy.hpp"
 #include "multigrid.hpp"
-#include "prmnonlinearform.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -49,7 +48,6 @@
 #include "plinearform.hpp"
 #include "pbilinearform.hpp"
 #include "pnonlinearform.hpp"
-#include "pprmnonlinearform.hpp"
 #endif
 
 #ifdef MFEM_USE_SIDRE

fem/nonlininteg.cpp

@@ -128,92 +128,6 @@ double BlockNonlinearFormIntegrator::GetElementEnergy(
    return 0.0;
 }
 
-double PrmBlockNonlinearFormIntegrator::GetElementEnergy(const Array<const FiniteElement *> &el,
-                                                         const Array<const FiniteElement *> &pel,
-                                                         ElementTransformation &Tr,
-                                                         const Array<const Vector *> &elfun,
-                                                         const Array<const Vector *> &pelfun)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::GetElementEnergy"
-               " is not overloaded!");
-    return 0.0;
-}
-
-void PrmBlockNonlinearFormIntegrator::AssembleFaceGrad(const Array<const FiniteElement *> &el1,
-                                                       const Array<const FiniteElement *> &el2,
-                                                       const Array<const FiniteElement *> &pel1,
-                                                       const Array<const FiniteElement *> &pel2,
-                                                       FaceElementTransformations &Tr,
-                                                       const Array<const Vector *> &elfun,
-                                                       const Array<const Vector *> &pelfun,
-                                                       const Array2D<DenseMatrix *> &elmats)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleFaceGrad"
-               " is not overloaded!");
-}
-
-void PrmBlockNonlinearFormIntegrator::AssembleElementGrad(const Array<const FiniteElement *> &el,
-                                                          const Array<const FiniteElement *> &pel,
-                                                          ElementTransformation &Tr,
-                                                          const Array<const Vector *> &elfun,
-                                                          const Array<const Vector *> &pelfun,
-                                                          const Array2D<DenseMatrix *> &elmats)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleElementGrad"
-               " is not overloaded!");
-}
-
-void PrmBlockNonlinearFormIntegrator::AssembleElementVector(const Array<const FiniteElement *> &el,
-                                                            const Array<const FiniteElement *> &pel,
-                                                            ElementTransformation &Tr,
-                                                            const Array<const Vector *> &elfun,
-                                                            const Array<const Vector *> &pelfun,
-                                                            const Array<Vector *> &elvec)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleElementVector"
-               " is not overloaded!");
-}
-
-void PrmBlockNonlinearFormIntegrator::AssembleFaceVector(const Array<const FiniteElement *> &el1,
-                                                         const Array<const FiniteElement *> &el2,
-                                                         const Array<const FiniteElement *> &pel1,
-                                                         const Array<const FiniteElement *> &pel2,
-                                                         FaceElementTransformations &Tr,
-                                                         const Array<const Vector *> &elfun,
-                                                         const Array<const Vector *> &pelfun,
-                                                         const Array<Vector *> &elvect)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssembleFaceVector"
-               " is not overloaded!");
-}
-
-
-void PrmBlockNonlinearFormIntegrator::AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                                               const Array<const FiniteElement *> &pel,
-                                                               ElementTransformation &Tr,
-                                                               const Array<const Vector *> &elfun,
-                                                               const Array<const Vector *> &alfun,
-                                                               const Array<const Vector *> &pelfun,
-                                                               const Array<Vector *> &elvec)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssemblePrmElementVector"
-               " is not overloaded!");
-}
-
-void PrmBlockNonlinearFormIntegrator::AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
-                                                         const Array<const FiniteElement *> &el2,
-                                                         const Array<const FiniteElement *> &pel1,
-                                                         const Array<const FiniteElement *> &pel2,
-                                                         FaceElementTransformations &Tr,
-                                                         const Array<const Vector *> &elfun,
-                                                         const Array<const Vector *> &alfun,
-                                                         const Array<const Vector *> &pelfun,
-                                                         const Array<Vector *> &elvect)
-{
-    mfem_error("PrmBlockNonlinearFormIntegrator::AssemblePrmFaceVector"
-               " is not overloaded!");
-}
-
 
 double InverseHarmonicModel::EvalW(const DenseMatrix &J) const
 {

fem/nonlininteg.hpp

@@ -130,80 +130,6 @@ public:
 };
 
 
-/** The abstract base class PrmBlockNonlinearFormIntegrator is
-    a generalization of the BlockNonlinearFormIntegrator class suitable
-    for block state and parameter vectors. */
-class PrmBlockNonlinearFormIntegrator
-{
-public:
-   /// Compute the local energy
-   virtual double GetElementEnergy(const Array<const FiniteElement *>&el,
-                                   const Array<const FiniteElement *>&pel,
-                                   ElementTransformation &Tr,
-                                   const Array<const Vector *>&elfun,
-                                   const Array<const Vector *>&pelfun);
-
-   /// Perform the local action of the BlockNonlinearFormIntegrator
-   virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
-                                      const Array<const FiniteElement *>&pel,
-                                      ElementTransformation &Tr,
-                                      const Array<const Vector *> &elfun,
-                                      const Array<const Vector *>&pelfun,
-                                      const Array<Vector *> &elvec);
-
-   virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
-                                   const Array<const FiniteElement *> &el2,
-                                   const Array<const FiniteElement *> &pel1,
-                                   const Array<const FiniteElement *> &pel2,
-                                   FaceElementTransformations &Tr,
-                                   const Array<const Vector *> &elfun,
-                                   const Array<const Vector *>&pelfun,
-                                   const Array<Vector *> &elvect);
-
-    /// Perform the local action on the parameters of the BlockNonlinearFormIntegrator
-    virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                       const Array<const FiniteElement *>&pel,
-                                       ElementTransformation &Tr,
-                                       const Array<const Vector *> &elfun,
-                                       const Array<const Vector *> &alfun,
-                                       const Array<const Vector *>&pelfun,
-                                       const Array<Vector *> &elvec);
-
-    virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
-                                    const Array<const FiniteElement *> &el2,
-                                    const Array<const FiniteElement *> &pel1,
-                                    const Array<const FiniteElement *> &pel2,
-                                    FaceElementTransformations &Tr,
-                                    const Array<const Vector *> &elfun,
-                                    const Array<const Vector *> &alfun,
-                                    const Array<const Vector *>&pelfun,
-                                    const Array<Vector *> &elvect);
-
-
-
-   /// Assemble the local gradient matrix
-   virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
-                                    const Array<const FiniteElement *>&pel,
-                                    ElementTransformation &Tr,
-                                    const Array<const Vector *> &elfun,
-                                    const Array<const Vector *>&pelfun,
-                                    const Array2D<DenseMatrix *> &elmats);
-
-   virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
-                                 const Array<const FiniteElement *>&el2,
-                                 const Array<const FiniteElement *> &pel1,
-                                 const Array<const FiniteElement *> &pel2,
-                                 FaceElementTransformations &Tr,
-                                 const Array<const Vector *> &elfun,
-                                 const Array<const Vector *>&pelfun,
-                                 const Array2D<DenseMatrix *> &elmats);
-
-
-   virtual ~PrmBlockNonlinearFormIntegrator() { }
-};
-
-
-
 /// Abstract class for hyperelastic models
 class HyperelasticModel
 {

fem/pprmnonlinearform.cpp

@@ -1,363 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_MPI
-
-#include "fem.hpp"
-
-
-namespace mfem
-{
-
-ParPrmBlockNonlinearForm::ParPrmBlockNonlinearForm(Array<ParFiniteElementSpace *> &pf,
-                                                   Array<ParFiniteElementSpace *> &ppf)
-    :PrmBlockNonlinearForm()
-{
-    pBlockGrad = nullptr;
-    SetParSpaces(pf,ppf);
-}
-
-void ParPrmBlockNonlinearForm::SetParSpaces(Array<ParFiniteElementSpace *> &pf,
-                                            Array<ParFiniteElementSpace *> &pprmf)
-{
-   delete pBlockGrad;
-   pBlockGrad = nullptr;
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-
-   Array<FiniteElementSpace *> serialSpaces(pf.Size());
-   Array<FiniteElementSpace *> prmserialSpaces(pprmf.Size());
-   for (int s=0; s<pf.Size(); s++)
-   {
-      serialSpaces[s] = (FiniteElementSpace *) pf[s];
-   }
-   for (int s=0; s<pprmf.Size(); s++)
-   {
-      prmserialSpaces[s] = (FiniteElementSpace *) pprmf[s];
-   }
-
-   SetSpaces(serialSpaces,prmserialSpaces);
-
-   phBlockGrad.SetSize(fes.Size(), fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2) = new OperatorHandle(Operator::Hypre_ParCSR);
-      }
-   }
-}
-
-ParFiniteElementSpace * ParPrmBlockNonlinearForm::ParFESpace(int k)
-{
-   return (ParFiniteElementSpace *)fes[k];
-}
-
-const ParFiniteElementSpace *ParPrmBlockNonlinearForm::ParFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)fes[k];
-}
-
-
-ParFiniteElementSpace * ParPrmBlockNonlinearForm::ParPrmFESpace(int k)
-{
-   return (ParFiniteElementSpace *)prmfes[k];
-}
-
-const ParFiniteElementSpace *ParPrmBlockNonlinearForm::ParPrmFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)prmfes[k];
-}
-
-// Here, rhs is a true dof vector
-void ParPrmBlockNonlinearForm::SetEssentialBC(const
-                                           Array<Array<int> *>&bdr_attr_is_ess,
-                                           Array<Vector *> &rhs)
-{
-   Array<Vector *> nullarray(fes.Size());
-   nullarray = NULL;
-
-   PrmBlockNonlinearForm::SetEssentialBC(bdr_attr_is_ess, nullarray);
-
-   for (int s = 0; s < fes.Size(); ++s)
-   {
-      if (rhs[s])
-      {
-         rhs[s]->SetSubVector(*ess_tdofs[s], 0.0);
-      }
-   }
-}
-
-void ParPrmBlockNonlinearForm::SetPrmEssentialBC(const
-                                           Array<Array<int> *>&bdr_attr_is_ess,
-                                           Array<Vector *> &rhs)
-{
-    Array<Vector *> nullarray(fes.Size());
-    nullarray = NULL;
-
-    PrmBlockNonlinearForm::SetPrmEssentialBC(bdr_attr_is_ess, nullarray);
-
-    for (int s = 0; s < prmfes.Size(); ++s)
-    {
-       if (rhs[s])
-       {
-          rhs[s]->SetSubVector(*prmess_tdofs[s], 0.0);
-       }
-    }
-}
-
-double ParPrmBlockNonlinearForm::GetEnergy(const Vector &x) const
-{
-   xs_true.Update(x.GetData(), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s = 0; s < fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   double enloc = PrmBlockNonlinearForm::GetEnergyBlocked(xs,xdv);
-   double englo = 0.0;
-
-   MPI_Allreduce(&enloc, &englo, 1, MPI_DOUBLE, MPI_SUM,
-                 ParFESpace(0)->GetComm());
-
-   return englo;
-}
-
-void ParPrmBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
-{
-   xs_true.Update(x.GetData(), block_trueOffsets);
-   ys_true.Update(y.GetData(), block_trueOffsets);
-   xs.Update(block_offsets);
-   ys.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   PrmBlockNonlinearForm::MultBlocked(xs, xdv, ys);
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->MultTranspose(
-         ys.GetBlock(s), ys_true.GetBlock(s));
-
-      ys_true.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
-   }
-}
-
-/// Block T-Vector to Block T-Vector
-void ParPrmBlockNonlinearForm::PrmMult(const Vector &x, Vector &y) const
-{
-    xs_true.Update(x.GetData(), prmblock_trueOffsets);
-    ys_true.Update(y.GetData(), prmblock_trueOffsets);
-    prmxs.Update(prmblock_offsets);
-    prmys.Update(prmblock_offsets);
-
-    for (int s=0; s<prmfes.Size(); ++s)
-    {
-       prmfes[s]->GetProlongationMatrix()->Mult(
-          xs_true.GetBlock(s), prmxs.GetBlock(s));
-    }
-
-    PrmBlockNonlinearForm::MultPrmBlocked(xsv,adv,xdv,prmys);
-
-    if (fnfi.Size() > 0)
-    {
-       MFEM_ABORT("TODO: assemble contributions from shared face terms");
-    }
-
-    for (int s=0; s<prmfes.Size(); ++s)
-    {
-       prmfes[s]->GetProlongationMatrix()->MultTranspose(
-          prmys.GetBlock(s), ys_true.GetBlock(s));
-
-       ys_true.GetBlock(s).SetSubVector(*prmess_tdofs[s], 0.0);
-    }
-
-}
-
-/// Return the local gradient matrix for the given true-dof vector x
-const BlockOperator & ParPrmBlockNonlinearForm::GetLocalGradient(
-   const Vector &x) const
-{
-   xs_true.Update(x.GetData(), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   PrmBlockNonlinearForm::ComputeGradientBlocked(xs,xdv); // (re)assemble Grad with b.c.
-
-   delete BlockGrad;
-   BlockGrad = new BlockOperator(block_offsets);
-
-   for (int i = 0; i < fes.Size(); ++i)
-   {
-      for (int j = 0; j < fes.Size(); ++j)
-      {
-         BlockGrad->SetBlock(i, j, Grads(i, j));
-      }
-   }
-   return *BlockGrad;
-}
-
-// Set the operator type id for the parallel gradient matrix/operator.
-void ParPrmBlockNonlinearForm::SetGradientType(Operator::Type tid)
-{
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2)->SetType(tid);
-      }
-   }
-}
-
-BlockOperator & ParPrmBlockNonlinearForm::GetGradient(const Vector &x) const
-{
-   if (pBlockGrad == NULL)
-   {
-      pBlockGrad = new BlockOperator(block_trueOffsets);
-   }
-
-   Array<const ParFiniteElementSpace *> pfes(fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      pfes[s1] = ParFESpace(s1);
-
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2)->Clear();
-      }
-   }
-
-   GetLocalGradient(x); // gradients are stored in 'Grads'
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                        Ph(phBlockGrad(s1,s2)->Type()),
-                        Rh(phBlockGrad(s1,s2)->Type());
-
-         if (s1 == s2)
-         {
-            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-            OperatorHandle Ae;
-            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
-         }
-         else
-         {
-            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                        pfes[s1]->GlobalVSize(),
-                                        pfes[s2]->GlobalVSize(),
-                                        pfes[s1]->GetDofOffsets(),
-                                        pfes[s2]->GetDofOffsets(),
-                                        Grads(s1,s2));
-            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
-         }
-
-         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
-      }
-   }
-
-   return *pBlockGrad;
-}
-
-ParPrmBlockNonlinearForm::~ParPrmBlockNonlinearForm()
-{
-   delete pBlockGrad;
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-}
-
-
-void ParPrmBlockNonlinearForm::SetStateFields(const Vector &xv) const
-{
-   xs_true.Update(xv.GetData(), block_trueOffsets);
-   xsv.Update(block_offsets);
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xsv.GetBlock(s));
-   }
-}
-
-
-void ParPrmBlockNonlinearForm::SetAdjointFields(const Vector &av) const
-{
-    xs_true.Update(av.GetData(), block_trueOffsets);
-    adv.Update(block_offsets);
-    for (int s=0; s<fes.Size(); ++s)
-    {
-       fes[s]->GetProlongationMatrix()->Mult(
-          xs_true.GetBlock(s), adv.GetBlock(s));
-    }
-}
-
-void ParPrmBlockNonlinearForm::SetPrmFields(const Vector &dv) const
-{
-    xs_true.Update(dv.GetData(),prmblock_trueOffsets);
-    xdv.Update(prmblock_offsets);
-    for (int s=0; s<prmfes.Size(); ++s)
-    {
-       prmfes[s]->GetProlongationMatrix()->Mult(
-          xs_true.GetBlock(s), xdv.GetBlock(s));
-    }
-}
-
-
-}
-
-#endif

fem/pprmnonlinearform.hpp

@@ -1,105 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_PPRMNONLINEARFORM
-#define MFEM_PPRMNONLINEARFORM
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_MPI
-
-#include "pgridfunc.hpp"
-#include "prmnonlinearform.hpp"
-
-namespace mfem
-{
-
-/** @brief A class representing a general parametric parallel block nonlinear operator
-    defined on the Cartesian product of multiple ParFiniteElementSpace%s. */
-/** The ParPrmBlockNonlinearForm takes as input, and returns as output, vectors on
-    the true dofs. */
-class ParPrmBlockNonlinearForm : public PrmBlockNonlinearForm
-{
-protected:
-   mutable BlockVector xs_true, ys_true;
-   mutable Array2D<OperatorHandle *> phBlockGrad;
-   mutable BlockOperator *pBlockGrad;
-
-public:
-   /// Computes the energy of the system
-   virtual double GetEnergy(const Vector &x) const;
-
-   /// Construct an empty ParPrmBlockNonlinearForm. Initialize with SetParSpaces().
-   ParPrmBlockNonlinearForm() : pBlockGrad(nullptr) { }
-
-   /** @brief Construct a ParPrmBlockNonlinearForm on the given set of
-       parametric and state ParFiniteElementSpace%s. */
-   ParPrmBlockNonlinearForm(Array<ParFiniteElementSpace *> &pf, Array<ParFiniteElementSpace *> &ppf );
-
-   /// Return the @a k-th parallel FE state space of the ParPrmBlockNonlinearForm.
-   ParFiniteElementSpace *ParFESpace(int k);
-   /** @brief Return the @a k-th parallel FE state space of the ParPrmBlockNonlinearForm
-       (const version). */
-   const ParFiniteElementSpace *ParFESpace(int k) const;
-
-   /// Return the @a k-th parallel FE parameters space of the ParPrmBlockNonlinearForm.
-   ParFiniteElementSpace *ParPrmFESpace(int k);
-   /** @brief Return the @a k-th parallel FE parameters space of the ParPrmBlockNonlinearForm
-       (const version). */
-   const ParFiniteElementSpace *ParPrmFESpace(int k) const;
-
-
-   /** @brief After a call to SetParSpaces(), the essential b.c. and the
-       gradient-type (if different from the default) must be set again. */
-   void SetParSpaces(Array<ParFiniteElementSpace *> &pf, Array<ParFiniteElementSpace *> &pprmf);
-
-   // Here, rhs is a true dof vector
-   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                                  Array<Vector *> &rhs);
-
-   // Here, rhs is a true dof vector
-   virtual void SetPrmEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                                  Array<Vector *> &rhs);
-
-
-   /// Block T-Vector to Block T-Vector
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   /// Block T-Vector to Block T-Vector
-   virtual void PrmMult(const Vector &x, Vector &y) const;
-
-   /// Return the local block gradient matrix for the given true-dof vector x
-   const BlockOperator &GetLocalGradient(const Vector &x) const;
-
-   virtual BlockOperator &GetGradient(const Vector &x) const;
-
-   /** @brief Set the operator type id for the blocks of the parallel gradient
-       matrix/operator. The default type is Operator::Hypre_ParCSR. */
-   void SetGradientType(Operator::Type tid);
-
-   /// Destructor.
-   virtual ~ParPrmBlockNonlinearForm();
-
-   /// Set the state fields
-   virtual void SetStateFields(const Vector &xv) const;
-
-   /// Set the adjoint fields
-   virtual void SetAdjointFields(const Vector &av) const;
-
-   /// Set the parameters/design fields
-   virtual void SetPrmFields(const Vector &dv) const;
-
-};
-
-}
-
-#endif
-#endif

fem/prmnonlinearform.hpp

@@ -1,231 +0,0 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_PRMNONLINEARFORM
-#define MFEM_PRMNONLINEARFORM
-
-#include "../config/config.hpp"
-#include "nonlininteg.hpp"
-#include "nonlinearform_ext.hpp"
-#include "bilinearform.hpp"
-#include "gridfunc.hpp"
-
-namespace mfem
-{
-
-/** @brief A class representing a general parametric block nonlinear operator
-    defined on the Cartesian product of multiple FiniteElementSpace%s. */
-class PrmBlockNonlinearForm : public Operator
-{
-protected:
-    /// FE spaces on which the form lives.
-    Array<FiniteElementSpace*> fes;
-
-    /// FE spaces for the parametric fields
-    Array<FiniteElementSpace*> prmfes;
-
-    int prmheight;
-    int prmwidth;
-
-
-    /// Set of Domain Integrators to be assembled (added).
-    Array<PrmBlockNonlinearFormIntegrator*> dnfi;
-
-    /// Set of interior face Integrators to be assembled (added).
-    Array<PrmBlockNonlinearFormIntegrator*> fnfi;
-
-    /// Set of Boundary Face Integrators to be assembled (added).
-    Array<PrmBlockNonlinearFormIntegrator*> bfnfi;
-    Array<Array<int>*>           bfnfi_marker;
-
-    /** Auxiliary block-vectors for wrapping input and output vectors or holding
-           GridFunction-like block-vector data (e.g. in parallel). */
-    mutable BlockVector xs, ys;
-    mutable BlockVector prmxs, prmys;
-
-
-    /** Auxiliary block-vectors for holding
-           GridFunction-like block-vector data (e.g. in parallel). */
-    mutable BlockVector xsv;
-
-    /** Auxiliary block-vectors for holding
-           GridFunction-like block-vector data for the parameter fields
-                                                (e.g. in parallel). */
-    mutable BlockVector xdv;
-    /** Auxiliary block-vectors for holding
-           GridFunction-like block-vector data for the adjoint fields
-                                                (e.g. in parallel). */
-    mutable BlockVector adv;
-
-
-    mutable Array2D<SparseMatrix*> Grads, cGrads;
-    mutable BlockOperator *BlockGrad;
-
-    // A list of the offsets
-    Array<int> block_offsets;
-    Array<int> block_trueOffsets;
-    // A list with the offsets for the parametric fields
-    Array<int> prmblock_offsets;
-    Array<int> prmblock_trueOffsets;
-
-    // Array of Arrays of tdofs for each space in 'fes'
-    Array<Array<int> *> ess_tdofs;
-
-    // Array of Arrays of tdofs for each space in 'prmfes'
-    Array<Array<int> *> prmess_tdofs;
-
-    /// Array of pointers to the prolongation matrix of fes, may be NULL
-    Array<const Operator *> P;
-
-    /// Array of pointers to the prolongation matrix of prmfes, may be NULL
-    Array<const Operator *> Pprm;
-
-    /// Array of results of dynamic-casting P to SparseMatrix pointer
-    Array<const SparseMatrix *> cP;
-
-    /// Array of results of dynamic-casting Pprm to SparseMatrix pointer
-    Array<const SparseMatrix *> cPprm;
-
-
-    /// Indicator if the Operator is part of a parallel run
-    bool is_serial = true;
-
-    /// Indicator if the Operator needs prolongation on assembly
-    bool needs_prolongation = false;
-
-    /// Indicator if the Operator needs prolongation on assembly
-    bool prmneeds_prolongation = false;
-
-
-    mutable BlockVector aux1, aux2;
-
-    mutable BlockVector prmaux1, prmaux2;
-
-    const BlockVector &Prolongate(const BlockVector &bx) const;
-
-    const BlockVector &PrmProlongate(const BlockVector &bx) const;
-
-    /// Specialized version of GetEnergy() for BlockVectors
-    //double GetEnergyBlocked(const BlockVector &bx) const;
-    double GetEnergyBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-
-    /// Specialized version of Mult() for BlockVector%s
-    /// Block L-Vector to Block L-Vector
-    void MultBlocked(const BlockVector &bx, const BlockVector &dx, BlockVector &by) const;
-
-    /// Specialized version of Mult() for BlockVector%s
-    /// Block L-Vector to Block L-Vector
-    /// bx - state vector, ax - adjoint vector, dx - parametric fields
-    /// dy = ax' d(residual(bx))/d(dx)
-    void MultPrmBlocked(const BlockVector &bx, const BlockVector & ax, const BlockVector &dx, BlockVector &dy) const;
-
-
-    /// Specialized version of GetGradient() for BlockVector
-    //void ComputeGradientBlocked(const BlockVector &bx) const;
-    void ComputeGradientBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-public:
-   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
-   PrmBlockNonlinearForm();
-
-   /// Construct a BlockNonlinearForm on the given set of FiniteElementSpace%s.
-   PrmBlockNonlinearForm(Array<FiniteElementSpace *> &f, Array<FiniteElementSpace *> &pf );
-
-   /// Return the @a k-th FE space of the PrmBlockNonlinearForm.
-   FiniteElementSpace *FESpace(int k) { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the PrmBlockNonlinearForm.
-   FiniteElementSpace *PrmFESpace(int k) { return prmfes[k]; }
-
-
-   /// Return the @a k-th FE space of the BlockNonlinearForm (const version).
-   const FiniteElementSpace *FESpace(int k) const { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the BlockNonlinearForm (const version).
-   const FiniteElementSpace *PrmFESpace(int k) const { return prmfes[k]; }
-
-   Array<PrmBlockNonlinearFormIntegrator*>& GetDNFI(){ return dnfi;}
-
-
-   /// (Re)initialize the PrmBlockNonlinearForm.
-   /** After a call to SetSpaces(), the essential b.c. must be set again. */
-   void SetSpaces(Array<FiniteElementSpace *> &f, Array<FiniteElementSpace *> &prmf);
-
-   /// Return the regular dof offsets.
-   const Array<int> &GetBlockOffsets() const { return block_offsets; }
-   /// Return the true-dof offsets.
-   const Array<int> &GetBlockTrueOffsets() const { return block_trueOffsets; }
-
-   /// Return the regular dof offsets for the parameters.
-   const Array<int> &PrmGetBlockOffsets() const { return prmblock_offsets; }
-   /// Return the true-dof offsets for the parameters.
-   const Array<int> &PrmGetBlockTrueOffsets() const { return prmblock_trueOffsets; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
-
-   /// Adds new Interior Face Integrator.
-   void AddInteriorFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
-
-   /// Adds new Boundary Face Integrator.
-   void AddBdrFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
-
-   /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
-       attributes. */
-   void AddBdrFaceIntegrator(PrmBlockNonlinearFormIntegrator *nlfi,
-                             Array<int> &bdr_marker);
-
-   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                               Array<Vector *> &rhs);
-
-   virtual void SetPrmEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                               Array<Vector *> &rhs);
-
-
-   virtual double GetEnergy(const Vector &x) const;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   virtual void Mult(const Vector &x, Vector &y) const;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   virtual void PrmMult(const Vector &x, Vector &t) const;
-
-   /// Method is only called in serial, the parallel version calls
-   /// GetGradientBlocked directly.
-   virtual Operator &GetGradient(const Vector &x) const;
-
-   /// Set the state fields
-   virtual void SetStateFields(const Vector &xv) const;
-
-   /// Set the adjoint fields
-   virtual void SetAdjointFields(const Vector &av) const;
-
-   /// Set the parameters/design fields
-   virtual void SetPrmFields(const Vector &dv) const;
-
-
-
-   /// Destructor.
-   virtual ~PrmBlockNonlinearForm();
-
-};
-
-}
-
-#endif
-
-

general/CMakeLists.txt

@@ -29,6 +29,7 @@ list(APPEND SRCS
   tic_toc.cpp
   tinyxml2.cpp
   version.cpp
+  hip.cpp
   )
 
 list(APPEND HDRS
@@ -59,6 +60,7 @@ list(APPEND HDRS
   tinyxml2.h
   text.hpp
   version.hpp
+  hip.hpp
   )
 
 if (MFEM_USE_MPI)

linalg/sparsemat.cpp

@@ -47,6 +47,19 @@ void SparseMatrix::InitCuSparse()
 #endif
 }
 
+void SparseMatrix::ClearCuSparse()
+{
+#ifdef MFEM_USE_CUDA
+   if (initBuffers)
+   {
+      cusparseDestroySpMat(matA_descr);
+      cusparseDestroyDnVec(vecX_descr);
+      cusparseDestroyDnVec(vecY_descr);
+      initBuffers = false;
+   }
+#endif
+}
+
 SparseMatrix::SparseMatrix(int nrows, int ncols)
    : AbstractSparseMatrix(nrows, (ncols >= 0) ? ncols : nrows),
      Rows(new RowNode *[nrows]),
@@ -282,15 +295,7 @@ void SparseMatrix::SetEmpty()
 #endif
    isSorted = false;
 
-#ifdef MFEM_USE_CUDA
-   if (initBuffers)
-   {
-      cusparseDestroySpMat(matA_descr);
-      cusparseDestroyDnVec(vecX_descr);
-      cusparseDestroyDnVec(vecY_descr);
-      initBuffers = false;
-   }
-#endif
+   ClearCuSparse();
 }
 
 int SparseMatrix::RowSize(const int i) const
@@ -3160,13 +3165,7 @@ void SparseMatrix::Destroy()
    delete At;
 
 #ifdef MFEM_USE_CUDA
-   if (initBuffers)
-   {
-      cusparseDestroySpMat(matA_descr);
-      cusparseDestroyDnVec(vecX_descr);
-      cusparseDestroyDnVec(vecY_descr);
-      initBuffers = false;
-   }
+   ClearCuSparse();
 #endif
 }
 

linalg/sparsemat.hpp

@@ -166,6 +166,10 @@ public:
    /// Clear the contents of the SparseMatrix.
    void Clear() { Destroy(); SetEmpty(); }
 
+   /** @brief Clear the CuSparse descriptors.
+       This must be called after releasing the device memory of A. */
+   void ClearCuSparse();
+
    /// Check if the SparseMatrix is empty.
    bool Empty() const { return (A == NULL) && (Rows == NULL); }
 

makefile

@@ -234,6 +234,9 @@ ifeq ($(MFEM_USE_CUDA),YES)
    endif
 endif
 
+MFEM_CXX += $(shell /usr/bin/python3-config --includes)
+#ALL_LIBS += $(shell /usr/bin/python3-config --ldflags)
+
 # HIP configuration
 ifeq ($(MFEM_USE_HIP),YES)
    ifeq ($(MFEM_USE_MPI),YES)
@@ -311,6 +314,8 @@ $(foreach dep,$(MFEM_DEPENDENCIES),$(eval $(call mfem_add_dependency,$(dep))))
 $(foreach dep,$(MFEM_LEGACY_DEPENDENCIES),$(eval $(call \
    mfem_add_legacy_dependency,$(dep))))
 
+ALL_LIBS += $(shell /usr/bin/python3-config --ldflags)
+
 # Timer option
 ifeq ($(MFEM_TIMER_TYPE),2)
    ALL_LIBS += $(POSIX_CLOCKS_LIB)

mesh/mesh_operators.cpp

@@ -145,6 +145,242 @@ void ThresholdRefiner::Reset()
    // marked_elements.SetSize(0); // not necessary
 }
 
+DRLRefiner::DRLRefiner(GridFunction& u_) : u(u_)
+{
+   int ret = _import_array();
+   if (ret < 0) {
+      printf("problem with import_array\n");
+   }
+
+   obs_x = 42;
+   obs_y = 42;
+
+   PyRun_SimpleString("import sys");
+   PyRun_SimpleString("sys.path.append('.')");
+
+   // This is a workaround for something in tensorflow that dies without it.
+   PyRun_SimpleString("if not hasattr(sys, 'argv'):\n"
+                      "  sys.argv  = ['']");
+
+   PyObject* eval_mod = PyImport_ImportModule("rllib_eval");
+   if (eval_mod == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   PyObject* eval_class = PyObject_GetAttrString(eval_mod, "Evaluator");
+   if (eval_class == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+   Py_DECREF(eval_mod);
+
+   PyObject* args = Py_BuildValue("()");
+   if (args == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   PyObject* eval_obj = PyEval_CallObject(eval_class, args);
+   if (eval_obj == NULL) {
+      PyErr_Print();
+      exit(1);
+   }
+
+   eval_method = PyObject_GetAttrString(eval_obj, "eval");
+   if (eval_method == 0) {
+      PyErr_Print();
+      exit(1);
+   }
+}
+
+int DRLRefiner::ApplyImpl(Mesh &mesh)
+{
+   marked_elements.SetSize(0);
+
+   const int imgsz = 42;
+
+   double ref_w = 2.0;
+   assert(obs_x == obs_y);
+   int ref_n = obs_x;
+   double ref_lo = 0.0 -ref_w;
+   double ref_hi = 1.0 +ref_w;
+   double ref_dx = (ref_hi -ref_lo)/ref_n;
+
+   for (int k = 0; k < mesh.GetNE(); k++) {
+
+      // assemble matrix of sample points, ref_space -> phys_space
+      ElementTransformation* trk = mesh.GetElementTransformation(k);
+      IntegrationPoint ipk;
+      Vector xk(2);
+      DenseMatrix m(2,imgsz*imgsz);
+      int c = 0;
+      for (int j = 0; j < obs_y; ++j) {
+         ipk.y = ref_lo +(j+0.5)*ref_dx;
+         for (int i = 0; i < obs_x; ++i) {
+            ipk.x = ref_lo +(i+0.5)*ref_dx;
+            trk->Transform(ipk, xk);
+            m.SetCol(c++,xk);
+         }
+      }
+
+      // phys_space -> elements, ips
+      Array<int> elems;
+      Array<IntegrationPoint> ips;
+      int n = mesh.FindPoints(m, elems, ips, false);
+
+      // Build observation from GridFunction using elements, ips
+      double* obs1 = new double[imgsz*imgsz];
+      n = 0;
+      bool complete = true;
+      for (int j = 0; j < obs_y; ++j) {
+         for (int i = 0; i < obs_x; ++i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs1[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               obs1[i*imgsz+j] = u.GetValue(el, ip);
+            }
+            n++;
+         }
+      }
+
+      // TODO: More efficient way to do the below mirroring
+
+      // invert i
+      double* obs2 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = 0; j < obs_y; ++j) {
+         for (int i = obs_x-1; i >= 0; --i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs2[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               obs2[i*imgsz+j] = u.GetValue(el, ip);
+            }
+            n++;
+         }
+      }
+
+      // invert j
+      double* obs3 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = obs_y-1; j >= 0; --j) {
+         for (int i = 0; i < obs_x; ++i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs3[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               obs3[i*imgsz+j] = u.GetValue(el, ip);
+            }
+            n++;
+         }
+      }
+
+      // invert i and j
+      double* obs4 = new double[imgsz*imgsz];
+      n = 0;
+      complete = true;
+      for (int j = obs_y-1; j >= 0; --j) {
+         for (int i = obs_x-1; i >= 0; --i) {
+            int el = elems[n];
+            if (el == -1) {
+               obs4[i*imgsz+j] = 0.0;
+               complete = false;
+            }
+            else {
+               IntegrationPoint& ip = ips[n];
+               obs4[i*imgsz+j] = u.GetValue(el, ip);
+            }
+            n++;
+         }
+      }
+
+      // apply policy: state -> action
+      bool refine = false;
+      if (complete) {
+
+         // convert to numpy array
+         npy_intp dims[3];
+         dims[0] = imgsz;
+         dims[1] = imgsz;
+         dims[2] = 1;
+
+         PyObject *pArray1 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs1));
+         if (pArray1 == NULL) printf("pArray1 NULL!\n");
+
+         PyObject *pArray2 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs2));
+         if (pArray2 == NULL) printf("pArray2 NULL!\n");
+
+         PyObject *pArray3 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs3));
+         if (pArray3 == NULL) printf("pArray3 NULL!\n");
+
+         PyObject *pArray4 = PyArray_SimpleNewFromData(
+            3, dims, NPY_DOUBLE, reinterpret_cast<void*>(obs4));
+         if (pArray4 == NULL) printf("pArray4 NULL!\n");
+
+         PyObject* action1 = PyObject_CallFunctionObjArgs(
+            eval_method, pArray1, nullptr);
+         PyObject* action2 = PyObject_CallFunctionObjArgs(
+            eval_method, pArray2, nullptr);
+         PyObject* action3 = PyObject_CallFunctionObjArgs(
+            eval_method, pArray3, nullptr);
+         PyObject* action4 = PyObject_CallFunctionObjArgs(
+            eval_method, pArray4, nullptr);
+
+         if (action1 == 0 || action2 == 0) {
+            PyErr_Print();
+            exit(1);
+         }
+
+         // parse integer return value
+         int action1_val;
+         PyArg_Parse(action1, "i", &action1_val);
+         int action2_val;
+         PyArg_Parse(action2, "i", &action2_val);
+         int action3_val;
+         PyArg_Parse(action3, "i", &action3_val);
+         int action4_val;
+         PyArg_Parse(action4, "i", &action4_val);
+         refine =
+            bool(action1_val) ||
+            bool(action2_val) ||
+            bool(action3_val) ||
+            bool(action4_val);
+      }
+
+      if (refine) {
+         marked_elements.Append(Refinement(k));
+      }
+   }
+
+   long int num_marked_elements = mesh.ReduceInt(marked_elements.Size());
+   printf("marked %d elements\n",num_marked_elements);
+   if (num_marked_elements == 0) { return STOP; }
+
+   bool nonconforming = true;
+   mesh.GeneralRefinement(marked_elements, nonconforming, nc_limit);
+   return CONTINUE + REFINED;
+}
+
+void DRLRefiner::Reset()
+{
+   current_sequence = -1;
+}
 
 int ThresholdDerefiner::ApplyImpl(Mesh &mesh)
 {

mesh/mesh_operators.hpp

@@ -19,6 +19,12 @@
 
 #include <limits>
 
+#define MFEM_USE_RLLIB
+#ifdef MFEM_USE_RLLIB
+#include <Python.h>
+#include "numpy/arrayobject.h"
+#endif
+
 namespace mfem
 {
 
@@ -254,6 +260,48 @@ public:
    virtual void Reset();
 };
 
+/** @brief Mesh refinement operator using a DRL policy from rllib.
+
+*/
+
+class DRLRefiner : public MeshOperator
+{
+protected:
+   GridFunction& u;
+
+   Array<Refinement> marked_elements;
+   long current_sequence;
+
+   int obs_x;
+   int obs_y;
+
+   int nc_limit;
+
+   PyObject* eval_method;
+
+   /** @brief Apply the operator to the mesh.
+       @return STOP if a stopping criterion is satisfied or no elements were
+       marked for refinement; REFINED + CONTINUE otherwise. */
+   virtual int ApplyImpl(Mesh &mesh);
+
+public:
+
+   /// Construct a DRLRefiner that will operate on u.
+   DRLRefiner(GridFunction &u);
+
+   // default destructor (virtual)
+
+   /** @brief Set the maximum ratio of refinement levels of adjacent elements
+       (0 = unlimited). */
+   void SetNCLimit(int nc_limit)
+   {
+      MFEM_ASSERT(nc_limit >= 0, "Invalid NC limit");
+      this->nc_limit = nc_limit;
+   }
+
+   virtual void Reset();
+};
+
 // TODO: BulkRefiner to refine a portion of the global error
 
 

miniapps/common/CMakeLists.txt

@@ -29,6 +29,9 @@ endif()
 if (MFEM_USE_CUDA)
   set_property(SOURCE ${SRCS} PROPERTY LANGUAGE CUDA)
 endif()
+if (MFEM_USE_HIP)
+  set_property(SOURCE ${SRCS} PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
+endif()
 
 convert_filenames_to_full_paths(SRCS)
 convert_filenames_to_full_paths(HDRS)
@@ -36,5 +39,5 @@ convert_filenames_to_full_paths(HDRS)
 set(MINIAPPS_COMMON_SOURCES ${MINIAPPS_COMMON_SOURCES} ${SRCS} PARENT_SCOPE)
 set(MINIAPPS_COMMON_HEADERS ${MINIAPPS_COMMON_HEADERS} ${HDRS} PARENT_SCOPE)
 
-add_library(mfem-common ${SRCS} ${HDRS})
+mfem_add_library(mfem-common ${SRCS} ${HDRS})
 target_link_libraries(mfem-common mfem)

tests/unit/CMakeLists.txt

@@ -68,9 +68,12 @@ set(UNIT_TESTS_SRCS
 if (MFEM_USE_CUDA)
    set_property(SOURCE ${UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
 endif()
+if (MFEM_USE_HIP)
+  set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS} PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
+endif()
 
 # All serial non-device unit tests are built into a single executable 'unit_tests'.
-add_executable(unit_tests unit_test_main.cpp ${UNIT_TESTS_SRCS})
+mfem_add_executable(unit_tests unit_test_main.cpp ${UNIT_TESTS_SRCS})
 # Unit tests need the ../../data directory.
 add_dependencies(unit_tests copy_data)
 target_link_libraries(unit_tests mfem)
@@ -104,13 +107,16 @@ set(SEDOV_TESTS_SRCS
 if (MFEM_USE_CUDA)
    set_property(SOURCE ${SEDOV_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
 endif()
+if (MFEM_USE_HIP)
+  set_property(SOURCE ${SEDOV_TESTS_SRCS} PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
+endif()
 
-add_executable(sedov_tests_cpu ${SEDOV_TESTS_SRCS})
+mfem_add_executable(sedov_tests_cpu ${SEDOV_TESTS_SRCS})
 target_compile_definitions(sedov_tests_cpu PUBLIC MFEM_SEDOV_TESTS=1)
 target_compile_definitions(sedov_tests_cpu PUBLIC MFEM_SEDOV_DEVICE="cpu")
 target_link_libraries(sedov_tests_cpu mfem)
 
-add_executable(sedov_tests_debug ${SEDOV_TESTS_SRCS})
+mfem_add_executable(sedov_tests_debug ${SEDOV_TESTS_SRCS})
 target_compile_definitions(sedov_tests_debug PUBLIC MFEM_SEDOV_TESTS=1)
 target_compile_definitions(sedov_tests_debug PUBLIC MFEM_SEDOV_DEVICE="debug")
 target_link_libraries(sedov_tests_debug mfem)