Merge branch 'master' into tribol-miniapp

This commit is contained in:
Eric B. Chin
2024-04-01 11:47:30 -07:00
635 changed files with 23693 additions and 15864 deletions
+1 -1
View File
@@ -65,7 +65,7 @@ jobs:
- name: GHCR Login
if: (github.event_name != 'pull_request')
uses: docker/login-action@v2
uses: docker/login-action@v3
with:
registry: ghcr.io
username: ${{ github.actor }}
+6
View File
@@ -27,6 +27,12 @@ jobs:
runs-on: ubuntu-latest
steps:
- name: Temporary workaround for sanitizer crashes
# See https://github.com/actions/runner-images/issues/9491
# The issue should be fixed in the next runner image for Ubuntu 22.04,
# see https://github.com/actions/runner-images/pull/9513
run: sudo sysctl vm.mmap_rnd_bits=28
- name: Cancel Previous Runs
uses: styfle/cancel-workflow-action@0.12.1
with:
+11 -1
View File
@@ -60,11 +60,18 @@ jobs:
./config/githooks/pre-push --release
continue-on-error: true
- name: math check
id: math
run: |
./config/githooks/pre-push --math
continue-on-error: true
- name: wrap-up
if: |
steps.copyright.outcome != 'success' ||
steps.license.outcome != 'success' ||
steps.release.outcome != 'success'
steps.release.outcome != 'success' ||
steps.math.outcome != 'success'
run: |
if [[ "${{ steps.copyright.outcome }}" != "success" ]]; then
echo "copyright check failed, unroll log for details"
@@ -75,6 +82,9 @@ jobs:
if [[ "${{ steps.release.outcome }}" != "success" ]]; then
echo "release check failed, unroll log for details"
fi
if [[ "${{ steps.math.outcome }}" != "success" ]]; then
echo "math check failed, unroll log for details"
fi
exit 1
code-style:
+1 -1
View File
@@ -87,7 +87,7 @@ report_baseline:
# We create an autotest-email.html file, because that's how we signal
# that there was an error / diff (temporary).
if [[ -f ${rundir}/${BASELINE_TEST}.err ]] || \
[[ -f ${rundir}/${BASELINE_TEST}-${SYS_TYPE}.diff ]]; then
[[ -f ${rundir}/${BASELINE_TEST}-${MACHINE_NAME}.diff ]]; then
cp ${rundir}/pipeline.txt ${rundir}/autotest-email.html
fi
msg="GitLab CI log for ${BASELINE_TEST} on ${MACHINE_NAME} ($(date +%Y-%m-%d))"
+3
View File
@@ -14,6 +14,9 @@
# locals
glob_err=${BASELINE_TEST}.err
base=${BASELINE_TEST}-${SYS_TYPE}
if [[ "${MACHINE_NAME}" == "quartz" ]]; then
base="${BASELINE_TEST}-${MACHINE_NAME}"
fi
base_diff=${base}.diff
base_patch=${base}.patch
base_out=${base}.out
+17
View File
@@ -11,6 +11,12 @@
Version 4.6.1 (development)
===========================
Meshing improvements
--------------------
- Introduced formulas for refinement of patches in NURBS meshes. Refinement by
arbitrary integer factors is also enabled, e.g. in the mesh-explorer miniapp.
NURBS coarsening and knot removal are also introduced.
Discretization improvements
---------------------------
- Introduced support for higher order non conformal Nedelec elements on
@@ -31,12 +37,20 @@ GPU support
New and updated examples and miniapps
-------------------------------------
- Added a new block solver in miniapp/solvers for the Darcy problem.
The new solver is based on a Bramble-Pasciak preconditioning. User can
use and implement their own preconditioner for the mass matrix.
- Added miniapp to demonstrate new elasticity integrator and unstructured element GPU support,
and a block diagonal preconditioner using low order refinement. Allows comparison with
currently existing legacy mode integrator. See miniapps/solvers/lor_elast.
Miscellaneous
-------------
- Added support for single and double precision, with corresponding hypre build.
Generalized the floating point type from `double` to `real_t`. For more
details see https://github.com/orgs/mfem/discussions/4207.
- The ReadCubit Genesis mesh importer has been rewritten to improve readability.
- Updated the Doxygen documentation style, which now requires Doxygen version
@@ -51,6 +65,9 @@ New and updated examples and miniapps
- Added a new mortar contact patch test miniapp using the Tribol interface
physics library (see https://github.com/LLNL/Tribol). See miniapps/tribol.
- RAJA backend will use seq_exec for serial loop execution when RAJA
v2023.06.00 and beyond is detected as loop_exec is deprecated.
Version 4.6, released on September 27, 2023
===========================================
+25 -8
View File
@@ -18,9 +18,10 @@ set(USER_CONFIG "${CMAKE_CURRENT_SOURCE_DIR}/config/user.cmake" CACHE PATH
"Path to optional user configuration file.")
# Require C++11 and disable compiler-specific extensions
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
set(CMAKE_CXX_STANDARD 11 CACHE STRING "C++ standard to use.")
set(CMAKE_CXX_STANDARD_REQUIRED ON CACHE BOOL
"Force the use of the chosen C++ standard.")
set(CMAKE_CXX_EXTENSIONS OFF CACHE BOOL "Enable C++ standard extensions.")
# Load user settings before the defaults - this way the defaults will not
# overwrite the user set options. If the user has not set all options, we still
@@ -94,7 +95,7 @@ if ((MFEM_USE_SUNDIALS OR
MFEM_USE_RAJA OR
MFEM_USE_UMPIRE) AND
("${CMAKE_CXX_STANDARD}" LESS "14"))
set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CXX_STANDARD 14 CACHE STRING "C++ standard to use." FORCE)
endif()
# Include xSDK default CMake file.
@@ -105,7 +106,8 @@ enable_language(CXX)
if (MINGW)
# MinGW GCC does not expose the functions jn/_jn, yn/_yn (used in Example
# 25/25p) unless we use '-std=gnu++11':
set(CMAKE_CXX_EXTENSIONS ON)
set(CMAKE_CXX_EXTENSIONS ON
CACHE BOOL "Enable C++ standard extensions." FORCE)
endif()
if (MFEM_USE_CUDA)
if (MFEM_USE_HIP)
@@ -116,9 +118,11 @@ if (MFEM_USE_CUDA)
set(CMAKE_CUDA_HOST_COMPILER ${CMAKE_CXX_COMPILER})
endif()
enable_language(CUDA)
set(CMAKE_CUDA_STANDARD ${CMAKE_CXX_STANDARD})
set(CMAKE_CUDA_STANDARD_REQUIRED ON)
set(CMAKE_CUDA_EXTENSIONS OFF)
set(CMAKE_CUDA_STANDARD ${CMAKE_CXX_STANDARD} CACHE STRING
"CUDA standard to use.")
set(CMAKE_CUDA_STANDARD_REQUIRED ON CACHE BOOL
"Force the use of the chosen CUDA standard.")
set(CMAKE_CUDA_EXTENSIONS OFF CACHE BOOL "Enable CUDA standard extensions.")
set(CUDA_FLAGS "--expt-extended-lambda")
if (CMAKE_VERSION VERSION_LESS 3.18.0)
set(CUDA_FLAGS "-arch=${CUDA_ARCH} ${CUDA_FLAGS}")
@@ -180,6 +184,19 @@ endif()
# Process configuration options
#-------------------------------------------------------------------------------
# MFEM_PRECISION -> MFEM_USE_SINGLE, MFEM_USE_DOUBLE
if (MFEM_PRECISION MATCHES "^(double|Double|DOUBLE)$")
set(MFEM_USE_SINGLE OFF)
set(MFEM_USE_DOUBLE ON)
elseif (MFEM_PRECISION MATCHES "^(single|Single|SINGLE)$")
set(MFEM_USE_SINGLE ON)
set(MFEM_USE_DOUBLE OFF)
else()
message(FATAL_ERROR " *** Invalid floating-point precision: "
"MFEM_PRECISION = ${MFEM_PRECISION}")
endif()
message(STATUS "Floating-point precision: MFEM_PRECISION = ${MFEM_PRECISION}")
# MFEM_DEBUG
if (CMAKE_BUILD_TYPE MATCHES "Debug|debug|DEBUG")
set(MFEM_DEBUG ON)
+5
View File
@@ -347,6 +347,9 @@ Before you can start, you need a GitHub account, here are a few suggestions:
- Code specifics
- All new public, protected, and private classes, methods, data members, and
functions have Doxygen-style documentation in source comments.
- Math formulas can be included in the Doxygen comments either with standard
LaTeX (`$..$` and `$$..$$`) for portions that need detailed explanation, or with
[Unicode](https://www.unicodeit.net/) or plain text when short or readable description is preferable.
- In addition to arguments and functionality, documentation should include the
current limitations of the code, any background information that is
implicitly assumed in the implementation, and the ownership and lifetime
@@ -360,6 +363,8 @@ Before you can start, you need a GitHub account, here are a few suggestions:
conflicted files in the commit message.
- All significant new features and changes should be documented in CHANGELOG.
- New examples and miniapps should have documentation on the MFEM webpage.
- The general floating-point type `real_t` should be used, rather than
`float` or `double`, except in special cases where only one is possible.
### Pull Requests
+10
View File
@@ -284,6 +284,15 @@ MFEM_USE_METIS = YES/NO
option in the library will be Cartesian partitioning with box meshes, and
thus most of the parallel examples and miniapps will fail.
MFEM_PRECISION = double/Double/DOUBLE/single/Single/SINGLE
Use single (float type) or double floating-point precision. In the
configuration header 'config/_config.hpp' this option is represented by
defining exactly one of the macros: MFEM_USE_DOUBLE, or MFEM_USE_SINGLE.
In the exported config files 'config.mk' and 'MFEMConfig.cmake', the option
is represented by the variables MFEM_USE_DOUBLE and MFEM_USE_SINGLE defined
as YES/NO (make) or ON/OFF (cmake). For more details see
https://github.com/orgs/mfem/discussions/4207
MFEM_DEBUG = YES/NO
Choose debug/optimized build. The debug build enables a number of messages
and consistency checks that may simplify bug-hunting.
@@ -968,6 +977,7 @@ The following options are equivalent to the GNU make options with the same name:
MFEM_USE_MPI
MFEM_USE_METIS - Set to ${MFEM_USE_MPI}, can be overwritten.
MFEM_PRECISION
MFEM_USE_LIBUNWIND
MFEM_USE_LAPACK
MFEM_THREAD_SAFE
+2
View File
@@ -18,6 +18,8 @@ set(MFEM_GIT_STRING "@MFEM_GIT_STRING@")
set(MFEM_USE_MPI @MFEM_USE_MPI@)
set(MFEM_USE_METIS @MFEM_USE_METIS@)
set(MFEM_USE_METIS_5 @MFEM_USE_METIS_5@)
set(MFEM_USE_DOUBLE @MFEM_USE_DOUBLE@)
set(MFEM_USE_SINGLE @MFEM_USE_SINGLE@)
set(MFEM_DEBUG @MFEM_DEBUG@)
set(MFEM_USE_EXCEPTIONS @MFEM_USE_EXCEPTIONS@)
set(MFEM_USE_ZLIB @MFEM_USE_ZLIB@)
+6
View File
@@ -46,6 +46,12 @@
// Requires an MPI compiler, and the libraries HYPRE and METIS.
#cmakedefine MFEM_USE_MPI
// Use double-precision floating point type
#cmakedefine MFEM_USE_DOUBLE
// Use single-precision floating point type
#cmakedefine MFEM_USE_SINGLE
// Enable debug checks in MFEM.
#cmakedefine MFEM_DEBUG
+3 -1
View File
@@ -79,7 +79,9 @@ if (HYPRE_FOUND AND HYPRE_USING_CUDA)
find_package(CUDAToolkit REQUIRED)
get_target_property(CUSPARSE_LIBRARIES CUDA::cusparse LOCATION)
get_target_property(CURAND_LIBRARIES CUDA::curand LOCATION)
list(APPEND HYPRE_LIBRARIES ${CUSPARSE_LIBRARIES} ${CURAND_LIBRARIES})
get_target_property(CUBLAS_LIBRARIES CUDA::cublas LOCATION)
list(APPEND HYPRE_LIBRARIES ${CUSPARSE_LIBRARIES} ${CURAND_LIBRARIES}
${CUBLAS_LIBRARIES})
set(HYPRE_LIBRARIES ${HYPRE_LIBRARIES} CACHE STRING
"HYPRE libraries + dependencies." FORCE)
message(STATUS "Updated HYPRE_LIBRARIES: ${HYPRE_LIBRARIES}")
+12 -11
View File
@@ -842,17 +842,18 @@ function(mfem_export_mk_files)
# Convert Boolean vars to YES/NO without writing the values to cache
set(CONFIG_MK_BOOL_VARS MFEM_USE_MPI MFEM_USE_METIS MFEM_USE_METIS_5
MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB MFEM_USE_LIBUNWIND
MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_LEGACY_OPENMP MFEM_USE_OPENMP
MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_SUITESPARSE
MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS MFEM_USE_STRUMPACK
MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_GNUTLS MFEM_USE_NETCDF
MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS
MFEM_USE_CONDUIT MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA
MFEM_USE_HIP MFEM_USE_RAJA MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER
MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO
MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG
MFEM_USE_TRIBOL MFEM_USE_MOONOLITH MFEM_USE_ALGOIM MFEM_USE_ENZYME)
MFEM_USE_SINGLE MFEM_USE_DOUBLE MFEM_DEBUG MFEM_USE_EXCEPTIONS
MFEM_USE_ZLIB MFEM_USE_LIBUNWIND MFEM_USE_LAPACK MFEM_THREAD_SAFE
MFEM_USE_LEGACY_OPENMP MFEM_USE_OPENMP MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS
MFEM_USE_SUITESPARSE MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS
MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_GNUTLS
MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE
MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB
MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA MFEM_USE_OCCA MFEM_USE_CEED
MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2
MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_ADFORWARD
MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL
MFEM_USE_MOONOLITH MFEM_USE_ALGOIM MFEM_USE_ENZYME)
foreach(var ${CONFIG_MK_BOOL_VARS})
if (${var})
set(${var} YES)
+55
View File
@@ -23,6 +23,61 @@
#include "_config.hpp"
#endif
namespace mfem
{
#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
#define MFEM_HOST_DEVICE __host__ __device__
#else
#define MFEM_HOST_DEVICE
#endif
// MFEM precision configuration
#if defined MFEM_USE_SINGLE && defined MFEM_USE_DOUBLE
#error "DOUBLE and SINGLE precision cannot both be specified"
#endif
#ifdef MFEM_USE_SINGLE
typedef float real_t;
#elif defined MFEM_USE_DOUBLE
typedef double real_t;
#else
#error "Either DOUBLE or SINGLE precision must be specified"
#endif
MFEM_HOST_DEVICE
constexpr real_t operator""_r(long double v)
{
return static_cast<real_t>(v);
}
MFEM_HOST_DEVICE
constexpr real_t operator""_r(unsigned long long v)
{
return static_cast<real_t>(v);
}
} // namespace mfem
// Return value for main function in examples that should be skipped by testing
// in some case. This return value prevents failures in testing.
#define MFEM_SKIP_RETURN_VALUE 242
// Request a global object to be instantiated for each thread in its TLS.
#define MFEM_THREAD_LOCAL thread_local
// MFEM_DEPRECATED macro to mark obsolete functions and methods
// see https://stackoverflow.com/questions/295120/c-mark-as-deprecated
#if defined(__GNUC__) || defined(__clang__)
#define MFEM_DEPRECATED __attribute__((deprecated))
#elif defined(_MSC_VER)
#define MFEM_DEPRECATED __declspec(deprecated)
#else
#pragma message("WARNING: You need to implement MFEM_DEPRECATED for this compiler")
#define MFEM_DEPRECATED
#endif
// Common configuration macros
#if (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) || defined(__clang__)
+6
View File
@@ -46,6 +46,12 @@
// Requires an MPI compiler, and the libraries HYPRE and METIS.
// #define MFEM_USE_MPI
// Use double-precision floating point type
// #define MFEM_USE_DOUBLE
// Use single-precision floating point type
// #define MFEM_USE_SINGLE
// Enable debug checks in MFEM.
// #define MFEM_DEBUG
+2
View File
@@ -18,6 +18,8 @@ MFEM_GIT_STRING = @MFEM_GIT_STRING@
MFEM_USE_MPI = @MFEM_USE_MPI@
MFEM_USE_METIS = @MFEM_USE_METIS@
MFEM_USE_METIS_5 = @MFEM_USE_METIS_5@
MFEM_USE_DOUBLE = @MFEM_USE_DOUBLE@
MFEM_USE_SINGLE = @MFEM_USE_SINGLE@
MFEM_DEBUG = @MFEM_DEBUG@
MFEM_USE_EXCEPTIONS = @MFEM_USE_EXCEPTIONS@
MFEM_USE_ZLIB = @MFEM_USE_ZLIB@
+2
View File
@@ -22,6 +22,8 @@ endif()
option(BUILD_SHARED_LIBS "Enable shared library build of MFEM" OFF)
option(MFEM_USE_MPI "Enable MPI parallel build" OFF)
option(MFEM_USE_METIS "Enable METIS usage" ${MFEM_USE_MPI})
set(MFEM_PRECISION "double" CACHE STRING
"Floating-point precision to use: single, or double")
option(MFEM_USE_EXCEPTIONS "Enable the use of exceptions" OFF)
option(MFEM_USE_ZLIB "Enable zlib for compressed data streams." OFF)
option(MFEM_USE_LIBUNWIND "Enable backtrace for errors." OFF)
+7 -1
View File
@@ -120,6 +120,7 @@ MFEM_MPI_NP = 4
MFEM_USE_MPI = NO
MFEM_USE_METIS = $(MFEM_USE_MPI)
MFEM_USE_METIS_5 = NO
MFEM_PRECISION = double
MFEM_DEBUG = NO
MFEM_USE_EXCEPTIONS = NO
MFEM_USE_ZLIB = NO
@@ -318,8 +319,13 @@ MPI_FORTRAN_LIB = -lmpifort
# MUMPS library configuration
MUMPS_DIR = @MFEM_DIR@/../MUMPS_5.5.0
MUMPS_OPT = -I$(MUMPS_DIR)/include
MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib -ldmumps\
MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib \
-lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
ifeq ($(MFEM_USE_SINGLE),YES)
MUMPS_LIB += -lsmumps
else
MUMPS_LIB += -ldmumps
endif
# STRUMPACK library configuration
STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build
+4
View File
@@ -9,7 +9,11 @@
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifdef MFEM_USE_SINGLE
#include "smumps_c.h"
#else
#include "dmumps_c.h"
#endif
#include <string>
#include <iostream>
#include <algorithm>
+17
View File
@@ -21,6 +21,7 @@ if [[ "${option}" == "--help" ]]; then
echo " --copyright"
echo " --license"
echo " --release"
echo " --math"
echo " --style"
echo " --history"
echo ""
@@ -63,6 +64,18 @@ if [[ "${option}" == "--release" || "${option}" == "" ]]; then
fi
fi
# math in doxygen check
math=true
if [[ "${option}" == "--math" || "${option}" == "" ]]; then
if grep '\\f' -R doc/CodeDocumentation.dox mfem.hpp config general linalg mesh fem examples miniapps | grep -v '\\frac' | grep -v fem/picojson.h | grep -v config/githooks/pre-push > matches.txt
then
echo "Please use $..$ and \$\$..\$\$ for LaTeX formulas in the following"
echo "comments instead of the Doxygen style \f$..\f$, \f[..\f], etc."
cat matches.txt
math=false
fi
fi
# wrap-up
code=0
if ! $copyright ; then
@@ -77,6 +90,10 @@ if ! $release ; then
echo "release check failed, unroll log for details"
code=1
fi
if ! $math ; then
echo "math in doxygen check failed, unroll log for details"
code=1
fi
# `code-style` is not just a check, it will actually reformat the code if
# necessary. This means that if one pushes while the repo is in dirty state
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+95
View File
@@ -0,0 +1,95 @@
MFEM NURBS mesh v1.1
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
# CUBE = 5
#
dimension
2
elements
2
1 3 0 1 4 5
2 3 1 2 3 4
boundary
6
3 1 0 1
3 1 4 5
3 1 1 2
3 1 3 4
1 1 5 0
2 1 2 3
edges
7
0 0 1
0 5 4
1 1 2
1 4 3
2 0 5
2 1 4
2 2 3
vertices
6
knotvectors
3
1 5 0 0 1 2 3 4 4
1 5 0 0 1 2 3 4 4
1 2 0 0 1 1
spacing
3
0 5 3 1 4 0 0 10.0
1 4 3 2 4 0 1 0.15 0.15
2 2 3 1 1 0 1 0.4
weights
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
FiniteElementSpace
FiniteElementCollection: NURBS1
VDim: 2
Ordering: 1
0 0
4 0
8 0
8 1
4 1
0 1
1 0
2 0
3 0
3 1
2 1
1 1
5 0
6 0
7 0
7 1
6 1
5 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,6 +1,6 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,6 +1,6 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,6 +1,6 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+381
View File
@@ -0,0 +1,381 @@
MFEM NURBS mesh v1.1
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
# CUBE = 5
#
dimension
3
elements
4
1 5 0 4 5 1 8 12 13 9
2 5 1 5 6 2 9 13 14 10
3 5 2 6 7 3 10 14 15 11
4 5 3 7 4 0 11 15 12 8
boundary
0
edges
32
0 0 1
0 1 2
0 2 3
0 3 0
0 4 5
0 5 6
0 6 7
0 7 4
1 0 4
1 1 5
1 2 6
1 3 7
0 8 9
0 9 10
0 10 11
0 11 8
0 12 13
0 13 14
0 14 15
0 15 12
1 8 12
1 9 13
1 10 14
1 11 15
2 0 8
2 1 9
2 2 10
2 3 11
2 4 12
2 5 13
2 6 14
2 7 15
vertices
16
knotvectors
3
2 3 0 0 0 1 1 1
2 3 0 0 0 1 1 1
2 5 0 0 0 1 1 2 2 2
spacing
2
0 5 3 1 1 0 1 10.0
1 5 3 1 1 0 1 10.0
weights
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
1
1
1
1
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
1
1
1
1
1
1
0.7071067811865475244
1
1
0.7071067811865475244
1
1
0.7071067811865475244
1
1
0.7071067811865475244
1
1
0.7071067811865475244
1
1
0.7071067811865475244
1
1
0.7071067811865475244
1
1
0.7071067811865475244
0.7071067811865475244
1
1
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.5
1
1
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.5
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.5
1
1
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.5
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.5
1
1
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.5
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.5
0.7071067811865475244
0.7071067811865475244
0.5
0.7071067811865475244
0.7071067811865475244
0.7071067811865475244
0.5
0.7071067811865475244
0.7071067811865475244
0.5
0.7071067811865475244
0.7071067811865475244
0.5
0.7071067811865475244
0.7071067811865475244
0.5
FiniteElementSpace
FiniteElementCollection: NURBS2
VDim: 3
Ordering: 1
2 0 0
0 2 0
-2 0 0
0 -2 0
3 0 0
0 3 0
-3 0 0
0 -3 0
10 0 18
10 2 20
10 0 22
10 -2 20
10 0 17
10 3 20
10 0 23
10 -3 20
2 2 0
-2 2 0
-2 -2 0
2 -2 0
3 3 0
-3 3 0
-3 -3 0
3 -3 0
2.5 0 0
0 2.5 0
-2.5 0 0
0 -2.5 0
10 2 18
10 2 22
10 -2 22
10 -2 18
10 3 17
10 3 23
10 -3 23
10 -3 17
10 0 17.5
10 2.5 20
10 0 22.5
10 -2.5 20
2 0 5
2 0 10
2 0 18
0 2 5
0 2 10
0 2 20
-2 0 5
-2 0 10
-2 0 22
0 -2 5
0 -2 10
0 -2 20
3 0 5
3 0 10
3 0 17
0 3 5
0 3 10
0 3 20
-3 0 5
-3 0 10
-3 0 23
0 -3 5
0 -3 10
0 -3 20
2.5 2.5 0
2.5 0 5
2.5 0 10
2.5 0 17.5
3 3 5
3 3 10
3 3 17
0 2.5 5
0 2.5 10
0 2.5 20
2 2 5
2 2 10
2 2 18
10 2.5 17.5
-2.5 2.5 0
-3 3 5
-3 3 10
-3 3 23
-2.5 0 5
-2.5 0 10
-2.5 0 22.5
-2 2 5
-2 2 10
-2 2 22
10 2.5 22.5
-2.5 -2.5 0
-3 -3 5
-3 -3 10
-3 -3 23
0 -2.5 5
0 -2.5 10
0 -2.5 20
-2 -2 5
-2 -2 10
-2 -2 22
10 -2.5 22.5
2.5 -2.5 0
3 -3 5
3 -3 10
3 -3 17
2 -2 5
2 -2 10
2 -2 18
10 -2.5 17.5
2.5 2.5 5
2.5 2.5 10
2.5 2.5 17.5
-2.5 2.5 5
-2.5 2.5 10
-2.5 2.5 22.5
-2.5 -2.5 5
-2.5 -2.5 10
-2.5 -2.5 22.5
2.5 -2.5 5
2.5 -2.5 10
2.5 -2.5 17.5
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+74
View File
@@ -0,0 +1,74 @@
MFEM NURBS mesh v1.1
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
# CUBE = 5
#
dimension
2
elements
1
1 3 0 1 2 3
boundary
4
1 1 0 1
1 1 2 3
1 1 3 0
1 1 1 2
edges
4
0 0 1
0 3 2
1 0 3
1 1 2
vertices
4
# Piecewise knotvectors
knotvectors
2
1 2 0 0 1 1
1 2 0 0 1 1
spacing
2
0 6 17 3
1 2 0
1 2
5 3 1 1 0 0
2 3 1 2 0 1
0.75 10.0 0.2
1 6 27 6
1 4 0
1 1 1 3
0 1 0 1
1 3 1 1 0 1
0 1 0 1
4 3 2 3 0 1
0.1 0.2 0.3
0.5
0.3 0.5
weights
1
1
1
1
FiniteElementSpace
FiniteElementCollection: NURBS1
VDim: 2
Ordering: 1
0 0
1 0
1 1
0 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM NURBS mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# SEGMENT = 1
# SQUARE = 3
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+1 -1
View File
@@ -1,7 +1,7 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
# MFEM Geometry Types (see fem/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
+11 -6
View File
@@ -1126,7 +1126,13 @@ IMAGE_PATH =
# need to set EXTENSION_MAPPING for the extension otherwise the files are not
# properly processed by doxygen.
INPUT_FILTER =
# Allow the use of standard LaTeX formulas in MFEM comments by replacing
# $$...$$ with \f[...\f] and $...$ with \f$...\f$ before running Doxygen.
# The regular expression ((?:.|\n)+?) is a lazy match for one or more
# characters, include newline. If present, Doxygen-style LaTeX commands
# such as \f$, \f[, etc., are left unchanged.
INPUT_FILTER = "perl -0777 -pe 's/\$\$((?:.|\n)+?)\$\$/\\f[\1\\f]/g; s/(?<!\\f)\$((?:.|\n)+?)\$/\\f\$\1\\f\$/g'"
# The FILTER_PATTERNS tag can be used to specify filters on a per file pattern
# basis. Doxygen will compare the file name with each pattern and apply the
@@ -1783,7 +1789,7 @@ USE_MATHJAX = YES
# The default value is: MathJax_2.
# This tag requires that the tag USE_MATHJAX is set to YES.
MATHJAX_VERSION = MathJax_2
MATHJAX_VERSION = MathJax_3
# When MathJax is enabled you can set the default output format to be used for
# the MathJax output. For more details about the output format see MathJax
@@ -1800,7 +1806,7 @@ MATHJAX_VERSION = MathJax_2
# The default value is: HTML-CSS.
# This tag requires that the tag USE_MATHJAX is set to YES.
MATHJAX_FORMAT = HTML-CSS
MATHJAX_FORMAT =
# When MathJax is enabled you need to specify the location relative to the HTML
# output directory using the MATHJAX_RELPATH option. The destination directory
@@ -1814,7 +1820,7 @@ MATHJAX_FORMAT = HTML-CSS
# - in case of MathJax version 3: https://cdn.jsdelivr.net/npm/mathjax@3
# This tag requires that the tag USE_MATHJAX is set to YES.
MATHJAX_RELPATH = http://cdn.mathjax.org/mathjax/latest
MATHJAX_RELPATH =
# The MATHJAX_EXTENSIONS tag can be used to specify one or more MathJax
# extension names that should be enabled during MathJax rendering. For example
@@ -1826,8 +1832,7 @@ MATHJAX_RELPATH = http://cdn.mathjax.org/mathjax/latest
# MATHJAX_EXTENSIONS = ams
# This tag requires that the tag USE_MATHJAX is set to YES.
MATHJAX_EXTENSIONS = TeX/AMSmath \
TeX/AMSsymbols
MATHJAX_EXTENSIONS = ams
# The MATHJAX_CODEFILE tag can be used to specify a file with javascript pieces
# of code that will be used on startup of the MathJax code. See the MathJax site
+12
View File
@@ -93,6 +93,15 @@ if (MFEM_USE_MPI)
)
endif()
# Examples that return MFEM_SKIP_RETURN_VALUE in some cases:
set(SKIP_TESTS)
if (HYPRE_USING_CUDA OR HYPRE_USING_HIP)
list(APPEND SKIP_TESTS ex19p.cpp ex28p.cpp)
endif()
if (MFEM_USE_SINGLE)
list(APPEND SKIP_TESTS ex33.cpp ex33p.cpp)
endif()
# Include the source directory where mfem.hpp and mfem-performance.hpp are.
include_directories(BEFORE ${PROJECT_BINARY_DIR})
@@ -102,6 +111,9 @@ add_mfem_examples(ALL_EXE_SRCS)
# Add a test for each example
if (MFEM_ENABLE_TESTING)
foreach(SRC_FILE ${ALL_EXE_SRCS})
if (SRC_FILE IN_LIST SKIP_TESTS)
continue()
endif()
get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
+41 -31
View File
@@ -62,7 +62,7 @@ protected:
BilinearForm M, S;
NonlinearForm H;
double viscosity;
real_t viscosity;
HyperelasticModel *model;
CGSolver M_solver; // Krylov solver for inverting the mass matrix M
@@ -84,16 +84,16 @@ protected:
public:
HyperelasticOperator(FiniteElementSpace &f, Array<int> &ess_bdr,
double visc, double mu, double K);
real_t visc, real_t mu, real_t K);
/// Compute the right-hand side of the ODE system.
virtual void Mult(const Vector &vx, Vector &dvx_dt) const;
/** Solve the Backward-Euler equation: k = f(x + dt*k, t), for the unknown k.
This is the only requirement for high-order SDIRK implicit integration.*/
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
virtual void ImplicitSolve(const real_t dt, const Vector &x, Vector &k);
double ElasticEnergy(const Vector &x) const;
double KineticEnergy(const Vector &v) const;
real_t ElasticEnergy(const Vector &x) const;
real_t KineticEnergy(const Vector &v) const;
void GetElasticEnergyDensity(const GridFunction &x, GridFunction &w) const;
virtual ~HyperelasticOperator();
@@ -109,7 +109,7 @@ private:
BilinearForm *M, *S;
NonlinearForm *H;
mutable SparseMatrix *Jacobian;
double dt;
real_t dt;
const Vector *v, *x;
mutable Vector w, z;
@@ -117,7 +117,7 @@ public:
ReducedSystemOperator(BilinearForm *M_, BilinearForm *S_, NonlinearForm *H_);
/// Set current dt, v, x values - needed to compute action and Jacobian.
void SetParameters(double dt_, const Vector *v_, const Vector *x_);
void SetParameters(real_t dt_, const Vector *v_, const Vector *x_);
/// Compute y = H(x + dt (v + dt k)) + M k + S (v + dt k).
virtual void Mult(const Vector &k, Vector &y) const;
@@ -141,7 +141,7 @@ private:
public:
ElasticEnergyCoefficient(HyperelasticModel &m, const GridFunction &x_)
: model(m), x(x_) { }
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
virtual ~ElasticEnergyCoefficient() { }
};
@@ -161,11 +161,11 @@ int main(int argc, char *argv[])
int ref_levels = 2;
int order = 2;
int ode_solver_type = 3;
double t_final = 300.0;
double dt = 3.0;
double visc = 1e-2;
double mu = 0.25;
double K = 5.0;
real_t t_final = 300.0;
real_t dt = 3.0;
real_t visc = 1e-2;
real_t mu = 0.25;
real_t K = 5.0;
bool visualization = true;
int vis_steps = 1;
@@ -309,13 +309,13 @@ int main(int argc, char *argv[])
<< " Press space (in the GLVis window) to resume it.\n";
}
double ee0 = oper.ElasticEnergy(x.GetTrueVector());
double ke0 = oper.KineticEnergy(v.GetTrueVector());
real_t ee0 = oper.ElasticEnergy(x.GetTrueVector());
real_t ke0 = oper.KineticEnergy(v.GetTrueVector());
cout << "initial elastic energy (EE) = " << ee0 << endl;
cout << "initial kinetic energy (KE) = " << ke0 << endl;
cout << "initial total energy (TE) = " << (ee0 + ke0) << endl;
double t = 0.0;
real_t t = 0.0;
oper.SetTime(t);
ode_solver->Init(oper);
@@ -324,7 +324,7 @@ int main(int argc, char *argv[])
bool last_step = false;
for (int ti = 1; !last_step; ti++)
{
double dt_real = min(dt, t_final - t);
real_t dt_real = min(dt, t_final - t);
ode_solver->Step(vx, t, dt_real);
@@ -332,8 +332,8 @@ int main(int argc, char *argv[])
if (last_step || (ti % vis_steps) == 0)
{
double ee = oper.ElasticEnergy(x.GetTrueVector());
double ke = oper.KineticEnergy(v.GetTrueVector());
real_t ee = oper.ElasticEnergy(x.GetTrueVector());
real_t ke = oper.KineticEnergy(v.GetTrueVector());
cout << "step " << ti << ", t = " << t << ", EE = " << ee << ", KE = "
<< ke << ", ΔTE = " << (ee+ke)-(ee0+ke0) << endl;
@@ -419,7 +419,7 @@ ReducedSystemOperator::ReducedSystemOperator(
dt(0.0), v(NULL), x(NULL), w(height), z(height)
{ }
void ReducedSystemOperator::SetParameters(double dt_, const Vector *v_,
void ReducedSystemOperator::SetParameters(real_t dt_, const Vector *v_,
const Vector *x_)
{
dt = dt_; v = v_; x = x_;
@@ -453,16 +453,26 @@ ReducedSystemOperator::~ReducedSystemOperator()
HyperelasticOperator::HyperelasticOperator(FiniteElementSpace &f,
Array<int> &ess_bdr, double visc,
double mu, double K)
: TimeDependentOperator(2*f.GetTrueVSize(), 0.0), fespace(f),
Array<int> &ess_bdr, real_t visc,
real_t mu, real_t K)
: TimeDependentOperator(2*f.GetTrueVSize(), (real_t) 0.0), fespace(f),
M(&fespace), S(&fespace), H(&fespace),
viscosity(visc), z(height/2)
{
const double rel_tol = 1e-8;
#if defined(MFEM_USE_DOUBLE)
const real_t rel_tol = 1e-8;
const real_t newton_abs_tol = 0.0;
#elif defined(MFEM_USE_SINGLE)
const real_t rel_tol = 1e-3;
const real_t newton_abs_tol = 1e-4;
#else
#error "Only single and double precision are supported!"
const real_t rel_tol = real_t(1);
const real_t newton_abs_tol = real_t(0);
#endif
const int skip_zero_entries = 0;
const double ref_density = 1.0; // density in the reference configuration
const real_t ref_density = 1.0; // density in the reference configuration
ConstantCoefficient rho0(ref_density);
M.AddDomainIntegrator(new VectorMassIntegrator(rho0));
M.Assemble(skip_zero_entries);
@@ -509,7 +519,7 @@ HyperelasticOperator::HyperelasticOperator(FiniteElementSpace &f,
newton_solver.SetOperator(*reduced_oper);
newton_solver.SetPrintLevel(1); // print Newton iterations
newton_solver.SetRelTol(rel_tol);
newton_solver.SetAbsTol(0.0);
newton_solver.SetAbsTol(newton_abs_tol);
newton_solver.SetMaxIter(10);
}
@@ -533,7 +543,7 @@ void HyperelasticOperator::Mult(const Vector &vx, Vector &dvx_dt) const
dx_dt = v;
}
void HyperelasticOperator::ImplicitSolve(const double dt,
void HyperelasticOperator::ImplicitSolve(const real_t dt,
const Vector &vx, Vector &dvx_dt)
{
int sc = height/2;
@@ -555,12 +565,12 @@ void HyperelasticOperator::ImplicitSolve(const double dt,
add(v, dt, dv_dt, dx_dt);
}
double HyperelasticOperator::ElasticEnergy(const Vector &x) const
real_t HyperelasticOperator::ElasticEnergy(const Vector &x) const
{
return H.GetEnergy(x);
}
double HyperelasticOperator::KineticEnergy(const Vector &v) const
real_t HyperelasticOperator::KineticEnergy(const Vector &v) const
{
return 0.5*M.InnerProduct(v, v);
}
@@ -581,7 +591,7 @@ HyperelasticOperator::~HyperelasticOperator()
}
double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
real_t ElasticEnergyCoefficient::Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{
model.SetTransformation(T);
@@ -601,7 +611,7 @@ void InitialDeformation(const Vector &x, Vector &y)
void InitialVelocity(const Vector &x, Vector &v)
{
const int dim = x.Size();
const double s = 0.1/64.;
const real_t s = 0.1/64.;
v = 0.0;
v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
+45 -35
View File
@@ -63,7 +63,7 @@ protected:
ParBilinearForm M, S;
ParNonlinearForm H;
double viscosity;
real_t viscosity;
HyperelasticModel *model;
HypreParMatrix *Mmat; // Mass matrix from ParallelAssemble()
@@ -86,16 +86,16 @@ protected:
public:
HyperelasticOperator(ParFiniteElementSpace &f, Array<int> &ess_bdr,
double visc, double mu, double K);
real_t visc, real_t mu, real_t K);
/// Compute the right-hand side of the ODE system.
virtual void Mult(const Vector &vx, Vector &dvx_dt) const;
/** Solve the Backward-Euler equation: k = f(x + dt*k, t), for the unknown k.
This is the only requirement for high-order SDIRK implicit integration.*/
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
virtual void ImplicitSolve(const real_t dt, const Vector &x, Vector &k);
double ElasticEnergy(const ParGridFunction &x) const;
double KineticEnergy(const ParGridFunction &v) const;
real_t ElasticEnergy(const ParGridFunction &x) const;
real_t KineticEnergy(const ParGridFunction &v) const;
void GetElasticEnergyDensity(const ParGridFunction &x,
ParGridFunction &w) const;
@@ -112,7 +112,7 @@ private:
ParBilinearForm *M, *S;
ParNonlinearForm *H;
mutable HypreParMatrix *Jacobian;
double dt;
real_t dt;
const Vector *v, *x;
mutable Vector w, z;
const Array<int> &ess_tdof_list;
@@ -122,7 +122,7 @@ public:
ParNonlinearForm *H_, const Array<int> &ess_tdof_list);
/// Set current dt, v, x values - needed to compute action and Jacobian.
void SetParameters(double dt_, const Vector *v_, const Vector *x_);
void SetParameters(real_t dt_, const Vector *v_, const Vector *x_);
/// Compute y = H(x + dt (v + dt k)) + M k + S (v + dt k).
virtual void Mult(const Vector &k, Vector &y) const;
@@ -146,7 +146,7 @@ private:
public:
ElasticEnergyCoefficient(HyperelasticModel &m, const ParGridFunction &x_)
: model(m), x(x_) { }
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
virtual ~ElasticEnergyCoefficient() { }
};
@@ -173,11 +173,11 @@ int main(int argc, char *argv[])
int par_ref_levels = 0;
int order = 2;
int ode_solver_type = 3;
double t_final = 300.0;
double dt = 3.0;
double visc = 1e-2;
double mu = 0.25;
double K = 5.0;
real_t t_final = 300.0;
real_t dt = 3.0;
real_t visc = 1e-2;
real_t mu = 0.25;
real_t K = 5.0;
bool adaptive_lin_rtol = true;
bool visualization = true;
int vis_steps = 1;
@@ -358,8 +358,8 @@ int main(int argc, char *argv[])
}
}
double ee0 = oper.ElasticEnergy(x_gf);
double ke0 = oper.KineticEnergy(v_gf);
real_t ee0 = oper.ElasticEnergy(x_gf);
real_t ke0 = oper.KineticEnergy(v_gf);
if (myid == 0)
{
cout << "initial elastic energy (EE) = " << ee0 << endl;
@@ -367,7 +367,7 @@ int main(int argc, char *argv[])
cout << "initial total energy (TE) = " << (ee0 + ke0) << endl;
}
double t = 0.0;
real_t t = 0.0;
oper.SetTime(t);
ode_solver->Init(oper);
@@ -376,7 +376,7 @@ int main(int argc, char *argv[])
bool last_step = false;
for (int ti = 1; !last_step; ti++)
{
double dt_real = min(dt, t_final - t);
real_t dt_real = min(dt, t_final - t);
ode_solver->Step(vx, t, dt_real);
@@ -386,8 +386,8 @@ int main(int argc, char *argv[])
{
v_gf.SetFromTrueVector(); x_gf.SetFromTrueVector();
double ee = oper.ElasticEnergy(x_gf);
double ke = oper.KineticEnergy(v_gf);
real_t ee = oper.ElasticEnergy(x_gf);
real_t ke = oper.KineticEnergy(v_gf);
if (myid == 0)
{
@@ -485,7 +485,7 @@ ReducedSystemOperator::ReducedSystemOperator(
ess_tdof_list(ess_tdof_list_)
{ }
void ReducedSystemOperator::SetParameters(double dt_, const Vector *v_,
void ReducedSystemOperator::SetParameters(real_t dt_, const Vector *v_,
const Vector *x_)
{
dt = dt_; v = v_; x = x_;
@@ -523,17 +523,27 @@ ReducedSystemOperator::~ReducedSystemOperator()
HyperelasticOperator::HyperelasticOperator(ParFiniteElementSpace &f,
Array<int> &ess_bdr, double visc,
double mu, double K)
: TimeDependentOperator(2*f.TrueVSize(), 0.0), fespace(f),
Array<int> &ess_bdr, real_t visc,
real_t mu, real_t K)
: TimeDependentOperator(2*f.TrueVSize(), (real_t) 0.0), fespace(f),
M(&fespace), S(&fespace), H(&fespace),
viscosity(visc), M_solver(f.GetComm()), newton_solver(f.GetComm()),
z(height/2)
{
const double rel_tol = 1e-8;
#if defined(MFEM_USE_DOUBLE)
const real_t rel_tol = 1e-8;
const real_t newton_abs_tol = 0.0;
#elif defined(MFEM_USE_SINGLE)
const real_t rel_tol = 1e-3;
const real_t newton_abs_tol = 1e-4;
#else
#error "Only single and double precision are supported!"
const real_t rel_tol = real_t(1);
const real_t newton_abs_tol = real_t(0);
#endif
const int skip_zero_entries = 0;
const double ref_density = 1.0; // density in the reference configuration
const real_t ref_density = 1.0; // density in the reference configuration
ConstantCoefficient rho0(ref_density);
M.AddDomainIntegrator(new VectorMassIntegrator(rho0));
M.Assemble(skip_zero_entries);
@@ -581,7 +591,7 @@ HyperelasticOperator::HyperelasticOperator(ParFiniteElementSpace &f,
newton_solver.SetOperator(*reduced_oper);
newton_solver.SetPrintLevel(1); // print Newton iterations
newton_solver.SetRelTol(rel_tol);
newton_solver.SetAbsTol(0.0);
newton_solver.SetAbsTol(newton_abs_tol);
newton_solver.SetAdaptiveLinRtol(2, 0.5, 0.9);
newton_solver.SetMaxIter(10);
}
@@ -607,7 +617,7 @@ void HyperelasticOperator::Mult(const Vector &vx, Vector &dvx_dt) const
dx_dt = v;
}
void HyperelasticOperator::ImplicitSolve(const double dt,
void HyperelasticOperator::ImplicitSolve(const real_t dt,
const Vector &vx, Vector &dvx_dt)
{
int sc = height/2;
@@ -629,17 +639,17 @@ void HyperelasticOperator::ImplicitSolve(const double dt,
add(v, dt, dv_dt, dx_dt);
}
double HyperelasticOperator::ElasticEnergy(const ParGridFunction &x) const
real_t HyperelasticOperator::ElasticEnergy(const ParGridFunction &x) const
{
return H.GetEnergy(x);
}
double HyperelasticOperator::KineticEnergy(const ParGridFunction &v) const
real_t HyperelasticOperator::KineticEnergy(const ParGridFunction &v) const
{
double loc_energy = 0.5*M.InnerProduct(v, v);
double energy;
MPI_Allreduce(&loc_energy, &energy, 1, MPI_DOUBLE, MPI_SUM,
fespace.GetComm());
real_t loc_energy = 0.5*M.InnerProduct(v, v);
real_t energy;
MPI_Allreduce(&loc_energy, &energy, 1, MPITypeMap<real_t>::mpi_type,
MPI_SUM, fespace.GetComm());
return energy;
}
@@ -660,7 +670,7 @@ HyperelasticOperator::~HyperelasticOperator()
}
double ElasticEnergyCoefficient::Eval(ElementTransformation &T,
real_t ElasticEnergyCoefficient::Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{
model.SetTransformation(T);
@@ -680,7 +690,7 @@ void InitialDeformation(const Vector &x, Vector &y)
void InitialVelocity(const Vector &x, Vector &v)
{
const int dim = x.Size();
const double s = 0.1/64.;
const real_t s = 0.1/64.;
v = 0.0;
v(dim-1) = s*x(0)*x(0)*(8.0-x(0));
+2 -2
View File
@@ -211,7 +211,7 @@ int main(int argc, char *argv[])
m->AddDomainIntegrator(new MassIntegrator(one));
m->Assemble();
// shift the eigenvalue corresponding to eliminated dofs to a large value
m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<real_t>::min());
m->Finalize();
HypreParMatrix *A = a->ParallelAssemble();
@@ -300,7 +300,7 @@ int main(int argc, char *argv[])
// 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
// parallel grid function to represent each of the eigenmodes returned by
// the solver.
Array<double> eigenvalues;
Array<real_t> eigenvalues;
lobpcg->Solve();
lobpcg->GetEigenvalues(eigenvalues);
ParGridFunction x(fespace);
+2 -2
View File
@@ -206,7 +206,7 @@ int main(int argc, char *argv[])
m->AddDomainIntegrator(new VectorMassIntegrator());
m->Assemble();
// shift the eigenvalue corresponding to eliminated dofs to a large value
m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<real_t>::min());
m->Finalize();
if (myid == 0)
{
@@ -247,7 +247,7 @@ int main(int argc, char *argv[])
// 10. Compute the eigenmodes and extract the array of eigenvalues. Define a
// parallel grid function to represent each of the eigenmodes returned by
// the solver.
Array<double> eigenvalues;
Array<real_t> eigenvalues;
lobpcg->Solve();
lobpcg->GetEigenvalues(eigenvalues);
ParGridFunction x(fespace);
+2 -2
View File
@@ -170,7 +170,7 @@ int main(int argc, char *argv[])
m->AddDomainIntegrator(new VectorFEMassIntegrator(one));
m->Assemble();
// shift the eigenvalue corresponding to eliminated dofs to a large value
m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<real_t>::min());
m->Finalize();
HypreParMatrix *A = a->ParallelAssemble();
@@ -198,7 +198,7 @@ int main(int argc, char *argv[])
// 10. Compute the eigenmodes and extract the array of eigenvalues. Define a
// parallel grid function to represent each of the eigenmodes returned by
// the solver.
Array<double> eigenvalues;
Array<real_t> eigenvalues;
ame->Solve();
ame->GetEigenvalues(eigenvalues);
ParGridFunction x(fespace);
+3 -3
View File
@@ -43,9 +43,9 @@ int main(int argc, char *argv[])
const char *mesh_file = "../data/star.mesh";
int ref_levels = -1;
int order = 1;
double sigma = -1.0;
double kappa = -1.0;
double eta = 0.0;
real_t sigma = -1.0;
real_t kappa = -1.0;
real_t eta = 0.0;
bool visualization = 1;
OptionsParser args(argc, argv);
+4 -4
View File
@@ -44,7 +44,7 @@ public:
pmesh(m),
pgf(f) {}
void MonitorSolution(int i, double norm, const Vector &x, bool final)
void MonitorSolution(int i, real_t norm, const Vector &x, bool final)
{
char vishost[] = "localhost";
int visport = 19916;
@@ -81,9 +81,9 @@ int main(int argc, char *argv[])
int ser_ref_levels = -1;
int par_ref_levels = 2;
int order = 1;
double sigma = -1.0;
double kappa = -1.0;
double eta = 0.0;
real_t sigma = -1.0;
real_t kappa = -1.0;
real_t eta = 0.0;
bool visualization = 1;
OptionsParser args(argc, argv);
+30 -30
View File
@@ -63,8 +63,8 @@ int problem;
int nfeatures;
// Prescribed time-dependent boundary and right-hand side functions.
double bdr_func(const Vector &pt, double t);
double rhs_func(const Vector &pt, double t);
real_t bdr_func(const Vector &pt, real_t t);
real_t rhs_func(const Vector &pt, real_t t);
// Update the finite element space, interpolate the solution and perform
// parallel load balancing.
@@ -79,9 +79,9 @@ int main(int argc, char *argv[])
nfeatures = 1;
const char *mesh_file = "../data/star-hilbert.mesh";
int order = 2;
double t_final = 1.0;
double max_elem_error = 5.0e-3;
double hysteresis = 0.15; // derefinement safety coefficient
real_t t_final = 1.0;
real_t max_elem_error = 5.0e-3;
real_t hysteresis = 0.15; // derefinement safety coefficient
int ref_levels = 0;
int nc_limit = 3; // maximum level of hanging nodes
bool visualization = true;
@@ -247,7 +247,7 @@ int main(int argc, char *argv[])
// refine the mesh as many times as necessary. Then we derefine any
// elements which have very small errors.
x = 0.0;
for (double time = 0.0; time < t_final + 1e-10; time += 0.01)
for (real_t time = 0.0; time < t_final + 1e-10; time += 0.01)
{
cout << "\nTime " << time << "\n\nRefinement:" << endl;
@@ -366,47 +366,47 @@ void UpdateProblem(Mesh &mesh, FiniteElementSpace &fespace,
}
const double alpha = 0.02;
const real_t alpha = 0.02;
// Spherical front with a Gaussian cross section and radius t
double front(double x, double y, double z, double t, int)
real_t front(real_t x, real_t y, real_t z, real_t t, int)
{
double r = sqrt(x*x + y*y + z*z);
real_t r = sqrt(x*x + y*y + z*z);
return exp(-0.5*pow((r - t)/alpha, 2));
}
double front_laplace(double x, double y, double z, double t, int dim)
real_t front_laplace(real_t x, real_t y, real_t z, real_t t, int dim)
{
double x2 = x*x, y2 = y*y, z2 = z*z, t2 = t*t;
double r = sqrt(x2 + y2 + z2);
double a2 = alpha*alpha, a4 = a2*a2;
real_t x2 = x*x, y2 = y*y, z2 = z*z, t2 = t*t;
real_t r = sqrt(x2 + y2 + z2);
real_t a2 = alpha*alpha, a4 = a2*a2;
return -exp(-0.5*pow((r - t)/alpha, 2)) / a4 *
(-2*t*(x2 + y2 + z2 - (dim-1)*a2/2)/r + x2 + y2 + z2 + t2 - dim*a2);
}
// Smooth spherical step function with radius t
double ball(double x, double y, double z, double t, int)
real_t ball(real_t x, real_t y, real_t z, real_t t, int)
{
double r = sqrt(x*x + y*y + z*z);
real_t r = sqrt(x*x + y*y + z*z);
return -atan(2*(r - t)/alpha);
}
double ball_laplace(double x, double y, double z, double t, int dim)
real_t ball_laplace(real_t x, real_t y, real_t z, real_t t, int dim)
{
double x2 = x*x, y2 = y*y, z2 = z*z, t2 = 4*t*t;
double r = sqrt(x2 + y2 + z2);
double a2 = alpha*alpha;
double den = pow(-a2 - 4*(x2 + y2 + z2 - 2*r*t) - t2, 2.0);
real_t x2 = x*x, y2 = y*y, z2 = z*z, t2 = 4*t*t;
real_t r = sqrt(x2 + y2 + z2);
real_t a2 = alpha*alpha;
real_t den = pow(-a2 - 4*(x2 + y2 + z2 - 2*r*t) - t2, 2.0);
return (dim == 2) ? 2*alpha*(a2 + t2 - 4*x2 - 4*y2)/r/den
/* */ : 4*alpha*(a2 + t2 - 4*r*t)/r/den;
}
// Composes several features into one function
template<typename F0, typename F1>
double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
real_t composite_func(const Vector &pt, real_t t, F0 f0, F1 f1)
{
int dim = pt.Size();
double x = pt(0), y = pt(1), z = 0.0;
real_t x = pt(0), y = pt(1), z = 0.0;
if (dim == 3) { z = pt(2); }
if (problem == 0)
@@ -417,11 +417,11 @@ double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
}
else
{
double sum = 0.0;
real_t sum = 0.0;
for (int i = 0; i < nfeatures; i++)
{
double x0 = 0.5*cos(2*M_PI * i / nfeatures);
double y0 = 0.5*sin(2*M_PI * i / nfeatures);
real_t x0 = 0.5*cos(2*M_PI * i / nfeatures);
real_t y0 = 0.5*sin(2*M_PI * i / nfeatures);
sum += f0(x - x0, y - y0, z, t, dim);
}
return sum;
@@ -429,11 +429,11 @@ double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
}
else
{
double sum = 0.0;
real_t sum = 0.0;
for (int i = 0; i < nfeatures; i++)
{
double x0 = 0.5*cos(2*M_PI * i / nfeatures + M_PI*t);
double y0 = 0.5*sin(2*M_PI * i / nfeatures + M_PI*t);
real_t x0 = 0.5*cos(2*M_PI * i / nfeatures + M_PI*t);
real_t y0 = 0.5*sin(2*M_PI * i / nfeatures + M_PI*t);
sum += f1(x - x0, y - y0, z, 0.25, dim);
}
return sum;
@@ -441,13 +441,13 @@ double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
}
// Exact solution, used for the Dirichlet BC.
double bdr_func(const Vector &pt, double t)
real_t bdr_func(const Vector &pt, real_t t)
{
return composite_func(pt, t, front, ball);
}
// Laplace of the exact solution, used for the right hand side.
double rhs_func(const Vector &pt, double t)
real_t rhs_func(const Vector &pt, real_t t)
{
return composite_func(pt, t, front_laplace, ball_laplace);
}
+30 -30
View File
@@ -68,8 +68,8 @@ int problem;
int nfeatures;
// Prescribed time-dependent boundary and right-hand side functions.
double bdr_func(const Vector &pt, double t);
double rhs_func(const Vector &pt, double t);
real_t bdr_func(const Vector &pt, real_t t);
real_t rhs_func(const Vector &pt, real_t t);
// Update the finite element space, interpolate the solution and perform
// parallel load balancing.
@@ -91,9 +91,9 @@ int main(int argc, char *argv[])
nfeatures = 1;
const char *mesh_file = "../data/star-hilbert.mesh";
int order = 2;
double t_final = 1.0;
double max_elem_error = 1.0e-4;
double hysteresis = 0.25; // derefinement safety coefficient
real_t t_final = 1.0;
real_t max_elem_error = 1.0e-4;
real_t hysteresis = 0.25; // derefinement safety coefficient
int ref_levels = 0;
int nc_limit = 3; // maximum level of hanging nodes
bool visualization = true;
@@ -282,7 +282,7 @@ int main(int argc, char *argv[])
// solve the problem on the current mesh, visualize the solution and
// refine the mesh as many times as necessary. Then we derefine any
// elements which have very small errors.
for (double time = 0.0; time < t_final + 1e-10; time += 0.01)
for (real_t time = 0.0; time < t_final + 1e-10; time += 0.01)
{
if (myid == 0)
{
@@ -427,47 +427,47 @@ void UpdateAndRebalance(ParMesh &pmesh, ParFiniteElementSpace &fespace,
}
const double alpha = 0.02;
const real_t alpha = 0.02;
// Spherical front with a Gaussian cross section and radius t
double front(double x, double y, double z, double t, int)
real_t front(real_t x, real_t y, real_t z, real_t t, int)
{
double r = sqrt(x*x + y*y + z*z);
real_t r = sqrt(x*x + y*y + z*z);
return exp(-0.5*pow((r - t)/alpha, 2));
}
double front_laplace(double x, double y, double z, double t, int dim)
real_t front_laplace(real_t x, real_t y, real_t z, real_t t, int dim)
{
double x2 = x*x, y2 = y*y, z2 = z*z, t2 = t*t;
double r = sqrt(x2 + y2 + z2);
double a2 = alpha*alpha, a4 = a2*a2;
real_t x2 = x*x, y2 = y*y, z2 = z*z, t2 = t*t;
real_t r = sqrt(x2 + y2 + z2);
real_t a2 = alpha*alpha, a4 = a2*a2;
return -exp(-0.5*pow((r - t)/alpha, 2)) / a4 *
(-2*t*(x2 + y2 + z2 - (dim-1)*a2/2)/r + x2 + y2 + z2 + t2 - dim*a2);
}
// Smooth spherical step function with radius t
double ball(double x, double y, double z, double t, int)
real_t ball(real_t x, real_t y, real_t z, real_t t, int)
{
double r = sqrt(x*x + y*y + z*z);
real_t r = sqrt(x*x + y*y + z*z);
return -atan(2*(r - t)/alpha);
}
double ball_laplace(double x, double y, double z, double t, int dim)
real_t ball_laplace(real_t x, real_t y, real_t z, real_t t, int dim)
{
double x2 = x*x, y2 = y*y, z2 = z*z, t2 = 4*t*t;
double r = sqrt(x2 + y2 + z2);
double a2 = alpha*alpha;
double den = pow(-a2 - 4*(x2 + y2 + z2 - 2*r*t) - t2, 2.0);
real_t x2 = x*x, y2 = y*y, z2 = z*z, t2 = 4*t*t;
real_t r = sqrt(x2 + y2 + z2);
real_t a2 = alpha*alpha;
real_t den = pow(-a2 - 4*(x2 + y2 + z2 - 2*r*t) - t2, 2.0);
return (dim == 2) ? 2*alpha*(a2 + t2 - 4*x2 - 4*y2)/r/den
/* */ : 4*alpha*(a2 + t2 - 4*r*t)/r/den;
}
// Composes several features into one function
template<typename F0, typename F1>
double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
real_t composite_func(const Vector &pt, real_t t, F0 f0, F1 f1)
{
int dim = pt.Size();
double x = pt(0), y = pt(1), z = 0.0;
real_t x = pt(0), y = pt(1), z = 0.0;
if (dim == 3) { z = pt(2); }
if (problem == 0)
@@ -478,11 +478,11 @@ double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
}
else
{
double sum = 0.0;
real_t sum = 0.0;
for (int i = 0; i < nfeatures; i++)
{
double x0 = 0.5*cos(2*M_PI * i / nfeatures);
double y0 = 0.5*sin(2*M_PI * i / nfeatures);
real_t x0 = 0.5*cos(2*M_PI * i / nfeatures);
real_t y0 = 0.5*sin(2*M_PI * i / nfeatures);
sum += f0(x - x0, y - y0, z, t, dim);
}
return sum;
@@ -490,11 +490,11 @@ double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
}
else
{
double sum = 0.0;
real_t sum = 0.0;
for (int i = 0; i < nfeatures; i++)
{
double x0 = 0.5*cos(2*M_PI * i / nfeatures + M_PI*t);
double y0 = 0.5*sin(2*M_PI * i / nfeatures + M_PI*t);
real_t x0 = 0.5*cos(2*M_PI * i / nfeatures + M_PI*t);
real_t y0 = 0.5*sin(2*M_PI * i / nfeatures + M_PI*t);
sum += f1(x - x0, y - y0, z, 0.25, dim);
}
return sum;
@@ -502,13 +502,13 @@ double composite_func(const Vector &pt, double t, F0 f0, F1 f1)
}
// Exact solution, used for the Dirichlet BC.
double bdr_func(const Vector &pt, double t)
real_t bdr_func(const Vector &pt, real_t t)
{
return composite_func(pt, t, front, ball);
}
// Laplace of the exact solution, used for the right hand side.
double rhs_func(const Vector &pt, double t)
real_t rhs_func(const Vector &pt, real_t t)
{
return composite_func(pt, t, front_laplace, ball_laplace);
}
+17 -17
View File
@@ -60,7 +60,7 @@ protected:
SparseMatrix Mmat, Kmat;
SparseMatrix *T; // T = M + dt K
double current_dt;
real_t current_dt;
CGSolver M_solver; // Krylov solver for inverting the mass matrix M
DSmoother M_prec; // Preconditioner for the mass matrix M
@@ -68,18 +68,18 @@ protected:
CGSolver T_solver; // Implicit solver for T = M + dt K
DSmoother T_prec; // Preconditioner for the implicit solver
double alpha, kappa;
real_t alpha, kappa;
mutable Vector z; // auxiliary vector
public:
ConductionOperator(FiniteElementSpace &f, double alpha, double kappa,
ConductionOperator(FiniteElementSpace &f, real_t alpha, real_t kappa,
const Vector &u);
virtual void Mult(const Vector &u, Vector &du_dt) const;
/** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
This is the only requirement for high-order SDIRK implicit integration.*/
virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
virtual void ImplicitSolve(const real_t dt, const Vector &u, Vector &k);
/// Update the diffusion BilinearForm K using the given true-dof vector `u`.
void SetParameters(const Vector &u);
@@ -87,7 +87,7 @@ public:
virtual ~ConductionOperator();
};
double InitialTemperature(const Vector &x);
real_t InitialTemperature(const Vector &x);
int main(int argc, char *argv[])
{
@@ -96,10 +96,10 @@ int main(int argc, char *argv[])
int ref_levels = 2;
int order = 2;
int ode_solver_type = 3;
double t_final = 0.5;
double dt = 1.0e-2;
double alpha = 1.0e-2;
double kappa = 0.5;
real_t t_final = 0.5;
real_t dt = 1.0e-2;
real_t alpha = 1.0e-2;
real_t kappa = 0.5;
bool visualization = true;
bool visit = false;
int vis_steps = 5;
@@ -246,7 +246,7 @@ int main(int argc, char *argv[])
// 8. Perform time-integration (looping over the time iterations, ti, with a
// time-step dt).
ode_solver->Init(oper);
double t = 0.0;
real_t t = 0.0;
bool last_step = false;
for (int ti = 1; !last_step; ti++)
@@ -293,12 +293,12 @@ int main(int argc, char *argv[])
return 0;
}
ConductionOperator::ConductionOperator(FiniteElementSpace &f, double al,
double kap, const Vector &u)
: TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
T(NULL), current_dt(0.0), z(height)
ConductionOperator::ConductionOperator(FiniteElementSpace &f, real_t al,
real_t kap, const Vector &u)
: TimeDependentOperator(f.GetTrueVSize(), (real_t) 0.0), fespace(f),
M(NULL), K(NULL), T(NULL), current_dt(0.0), z(height)
{
const double rel_tol = 1e-8;
const real_t rel_tol = 1e-8;
M = new BilinearForm(&fespace);
M->AddDomainIntegrator(new MassIntegrator());
@@ -336,7 +336,7 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
M_solver.Mult(z, du_dt);
}
void ConductionOperator::ImplicitSolve(const double dt,
void ConductionOperator::ImplicitSolve(const real_t dt,
const Vector &u, Vector &du_dt)
{
// Solve the equation:
@@ -382,7 +382,7 @@ ConductionOperator::~ConductionOperator()
delete K;
}
double InitialTemperature(const Vector &x)
real_t InitialTemperature(const Vector &x)
{
if (x.Norml2() < 0.5)
{
+17 -17
View File
@@ -62,7 +62,7 @@ protected:
HypreParMatrix Mmat;
HypreParMatrix Kmat;
HypreParMatrix *T; // T = M + dt K
double current_dt;
real_t current_dt;
CGSolver M_solver; // Krylov solver for inverting the mass matrix M
HypreSmoother M_prec; // Preconditioner for the mass matrix M
@@ -70,18 +70,18 @@ protected:
CGSolver T_solver; // Implicit solver for T = M + dt K
HypreSmoother T_prec; // Preconditioner for the implicit solver
double alpha, kappa;
real_t alpha, kappa;
mutable Vector z; // auxiliary vector
public:
ConductionOperator(ParFiniteElementSpace &f, double alpha, double kappa,
ConductionOperator(ParFiniteElementSpace &f, real_t alpha, real_t kappa,
const Vector &u);
virtual void Mult(const Vector &u, Vector &du_dt) const;
/** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
This is the only requirement for high-order SDIRK implicit integration.*/
virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
virtual void ImplicitSolve(const real_t dt, const Vector &u, Vector &k);
/// Update the diffusion BilinearForm K using the given true-dof vector `u`.
void SetParameters(const Vector &u);
@@ -89,7 +89,7 @@ public:
virtual ~ConductionOperator();
};
double InitialTemperature(const Vector &x);
real_t InitialTemperature(const Vector &x);
int main(int argc, char *argv[])
{
@@ -105,10 +105,10 @@ int main(int argc, char *argv[])
int par_ref_levels = 1;
int order = 2;
int ode_solver_type = 3;
double t_final = 0.5;
double dt = 1.0e-2;
double alpha = 1.0e-2;
double kappa = 0.5;
real_t t_final = 0.5;
real_t dt = 1.0e-2;
real_t alpha = 1.0e-2;
real_t kappa = 0.5;
bool visualization = true;
bool visit = false;
int vis_steps = 5;
@@ -313,7 +313,7 @@ int main(int argc, char *argv[])
// 10. Perform time-integration (looping over the time iterations, ti, with a
// time-step dt).
ode_solver->Init(oper);
double t = 0.0;
real_t t = 0.0;
bool last_step = false;
for (int ti = 1; !last_step; ti++)
@@ -382,13 +382,13 @@ int main(int argc, char *argv[])
return 0;
}
ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
double kap, const Vector &u)
: TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
T(NULL), current_dt(0.0),
ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, real_t al,
real_t kap, const Vector &u)
: TimeDependentOperator(f.GetTrueVSize(), (real_t) 0.0), fespace(f),
M(NULL), K(NULL), T(NULL), current_dt(0.0),
M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
{
const double rel_tol = 1e-8;
const real_t rel_tol = 1e-8;
M = new ParBilinearForm(&fespace);
M->AddDomainIntegrator(new MassIntegrator());
@@ -427,7 +427,7 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
M_solver.Mult(z, du_dt);
}
void ConductionOperator::ImplicitSolve(const double dt,
void ConductionOperator::ImplicitSolve(const real_t dt,
const Vector &u, Vector &du_dt)
{
// Solve the equation:
@@ -473,7 +473,7 @@ ConductionOperator::~ConductionOperator()
delete K;
}
double InitialTemperature(const Vector &x)
real_t InitialTemperature(const Vector &x)
{
if (x.Norml2() < 0.5)
{
+8 -8
View File
@@ -69,7 +69,7 @@ public:
void SetDisplacement(GridFunction &u_) { u = &u_; }
void SetComponent(int i, int j) { si = i; sj = j; }
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
};
// Simple GLVis visualization manager.
@@ -104,8 +104,8 @@ int main(int argc, char *argv[])
const char *mesh_file = "../data/beam-tri.mesh";
int ref_levels = -1;
int order = 1;
double alpha = -1.0;
double kappa = -1.0;
real_t alpha = -1.0;
real_t kappa = -1.0;
bool visualization = 1;
OptionsParser args(argc, argv);
@@ -245,7 +245,7 @@ int main(int argc, char *argv[])
// solve the system Ax=b with PCG for the symmetric formulation, or GMRES
// for the non-symmetric.
GSSmoother M(A);
const double rtol = 1e-6;
const real_t rtol = 1e-6;
if (alpha == -1.0)
{
PCG(A, M, B, X, 3, 5000, rtol*rtol, 0.0);
@@ -337,17 +337,17 @@ void InitDisplacement(const Vector &x, Vector &u)
}
double StressCoefficient::Eval(ElementTransformation &T,
real_t StressCoefficient::Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{
MFEM_ASSERT(u != NULL, "displacement field is not set");
double L = lambda.Eval(T, ip);
double M = mu.Eval(T, ip);
real_t L = lambda.Eval(T, ip);
real_t M = mu.Eval(T, ip);
u->GetVectorGradient(T, grad);
if (si == sj)
{
double div_u = grad.Trace();
real_t div_u = grad.Trace();
return L*div_u + 2*M*grad(si,si);
}
else
+8 -8
View File
@@ -69,7 +69,7 @@ public:
void SetDisplacement(GridFunction &u_) { u = &u_; }
void SetComponent(int i, int j) { si = i; sj = j; }
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip);
};
// Simple GLVis visualization manager.
@@ -108,8 +108,8 @@ int main(int argc, char *argv[])
int ser_ref_levels = -1;
int par_ref_levels = 1;
int order = 1;
double alpha = -1.0;
double kappa = -1.0;
real_t alpha = -1.0;
real_t kappa = -1.0;
bool amg_elast = false;
bool visualization = 1;
@@ -268,7 +268,7 @@ int main(int argc, char *argv[])
// 11. Define a simple symmetric Gauss-Seidel preconditioner and use it to
// solve the system Ax=b with PCG for the symmetric formulation, or GMRES
// for the non-symmetric.
const double rtol = 1e-6;
const real_t rtol = 1e-6;
HypreBoomerAMG amg(A);
if (amg_elast)
{
@@ -376,17 +376,17 @@ void InitDisplacement(const Vector &x, Vector &u)
}
double StressCoefficient::Eval(ElementTransformation &T,
real_t StressCoefficient::Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{
MFEM_ASSERT(u != NULL, "displacement field is not set");
double L = lambda.Eval(T, ip);
double M = mu.Eval(T, ip);
real_t L = lambda.Eval(T, ip);
real_t M = mu.Eval(T, ip);
u->GetVectorGradient(T, grad);
if (si == sj)
{
double div_u = grad.Trace();
real_t div_u = grad.Trace();
return L*div_u + 2*M*grad(si,si);
}
else
+10 -10
View File
@@ -52,11 +52,11 @@ int problem;
// Equation constant parameters.
const int num_equation = 4;
const double specific_heat_ratio = 1.4;
const double gas_constant = 1.0;
const real_t specific_heat_ratio = 1.4;
const real_t gas_constant = 1.0;
// Maximum characteristic speed (updated by integrators)
double max_char_speed;
real_t max_char_speed;
int main(int argc, char *argv[])
{
@@ -66,9 +66,9 @@ int main(int argc, char *argv[])
int ref_levels = 1;
int order = 3;
int ode_solver_type = 4;
double t_final = 2.0;
double dt = -0.01;
double cfl = 0.3;
real_t t_final = 2.0;
real_t dt = -0.01;
real_t cfl = 0.3;
bool visualization = true;
int vis_steps = 50;
@@ -228,7 +228,7 @@ int main(int argc, char *argv[])
}
// Determine the minimum element size.
double hmin = 0.0;
real_t hmin = 0.0;
if (cfl > 0)
{
hmin = mesh.GetElementSize(0, 1);
@@ -242,7 +242,7 @@ int main(int argc, char *argv[])
tic_toc.Clear();
tic_toc.Start();
double t = 0.0;
real_t t = 0.0;
euler.SetTime(t);
ode_solver->Init(euler);
@@ -260,7 +260,7 @@ int main(int argc, char *argv[])
bool done = false;
for (int ti = 0; !done; )
{
double dt_real = min(dt, t_final - t);
real_t dt_real = min(dt, t_final - t);
ode_solver->Step(sol, t, dt_real);
if (cfl > 0)
@@ -298,7 +298,7 @@ int main(int argc, char *argv[])
// 10. Compute the L2 solution error summed for all components.
if (t_final == 2.0)
{
const double error = sol.ComputeLpError(2, u0);
const real_t error = sol.ComputeLpError(2, u0);
cout << "Solution error: " << error << endl;
}

Some files were not shown because too many files have changed in this diff Show More