Compare commits

..
Author SHA1 Message Date
Tucker Babcock 0c6f51852b make style 2022-06-30 12:33:11 -06:00
Tucker Babcock 737883992d update linear form move ctor/assign to correctly move recently added LinearFormExtension member ext 2022-06-30 12:30:45 -06:00
Tucker Babcock d896f164fe Merge branch 'master' into move-nlf-dev 2022-06-30 12:27:48 -06:00
Tucker Babcock d60c544e0f merge master into move-nlf-dev 2022-03-24 10:11:05 -06:00
Tucker Babcock c5fbe6cb27 correctly handle moving the bilinear form integ markers 2022-01-19 10:06:55 -07:00
Tucker Babcock 981304f4f3 add element_matrices to bilinearform move constructor, I had forgotten it before 2022-01-17 15:14:17 -07:00
Tucker Babcock 196125ef2e make style 2022-01-15 20:08:31 -07:00
Tucker Babcock 57d196a963 reorder bilinear form constructor 2022-01-15 20:07:54 -07:00
Tucker Babcock e9b4e2534b first draft of bilinear form move semantics are complete, needs to be tested though 2022-01-15 19:18:33 -07:00
Tucker Babcock 96686455f6 Merge branch 'master' into move-nlf-dev 2022-01-15 18:27:24 -07:00
Tucker Babcock f363bf7468 working on adding move semantics to bilinear forms 2022-01-15 18:27:02 -07:00
Tucker Babcock bd6ac34b97 set assembly level of other nonlinear form to LEGACY in move assignment operator 2021-11-19 13:49:49 -07:00
Tucker Babcock 050bc5aacf update nonlinear form move assignment/ctors to use mfem::Swap for integrator and marker arrays, also add support for moving domain_integs_marker 2021-11-19 13:46:57 -07:00
Tucker Babcock 751e200d83 make style 2021-11-19 13:14:41 -07:00
Tucker Babcock d9fda017d4 merge master into move-nlf-dev 2021-11-19 13:13:36 -07:00
Tucker Babcock 2516c54386 add support for nonlinear form domain integrators restricted to certain attributes 2021-11-19 13:09:30 -07:00
Tucker Babcock 54e26b603a make style 2021-11-19 12:19:16 -07:00
Tucker Babcock d189d70da4 renamed nonlinear form integrator arrays to match bilinear form and linear form with more expressive names. Added method to add domain integrator restricted to certain element types for nonlinear form to match method for bilinear form. Need to update Mult, GetGradient, GetEnergy to use this and check for attributes 2021-11-19 12:18:36 -07:00
Tucker Babcock 7fd3ed3a94 make style 2021-11-09 15:23:59 -07:00
Tucker Babcock cce25fced8 added LinearFrom and ParLinearForm move ctors and assignment operators, but need to test 2021-11-09 15:20:21 -07:00
Tucker Babcock e25fb02eee add non-default move ctor/assign for ParNonlinearForm 2021-10-15 09:12:49 -06:00
Tucker Babcock c4f67c7dff explicitly defaulting ParNonlinearForm move ctor/assign 2021-10-14 16:40:21 -06:00
Tucker Babcock e5e92be833 make style 2021-10-14 15:26:30 -06:00
Tucker Babcock 88370b7fdc reverting all operator changes, returning *this in nlf move assignment operator 2021-10-14 15:20:11 -06:00
Tucker Babcock fc363ca778 added move ctors and assign operators to NonlinearForm classes, and explicitly delete copy ctor and assignment for NonlinearForm 2021-10-14 14:46:05 -06:00
189 changed files with 2502 additions and 9814 deletions
+2 -3
View File
@@ -1,4 +1,4 @@
name: "Docker"
name: Build Deploy Container
on:
@@ -28,8 +28,7 @@ jobs:
matrix:
# Dockerfiles to build, a matrix supports future expanded builds
container: [["config/docker/Dockerfile.base", "ghcr.io/mfem/mfem-ubuntu-base"],
["config/docker/Dockerfile", "ghcr.io/mfem/mfem-ubuntu"]]
container: [["config/docker/Dockerfile", "ghcr.io/mfem/mfem-ubuntu-base"]]
runs-on: ubuntu-latest
name: Build
+23 -23
View File
@@ -10,7 +10,7 @@
# CONTRIBUTING.md for details.
# In this CI section, we build different variants of mfem and run test on them.
name: "Tests"
name: builds-and-tests
# Github actions can use the default "GITHUB_TOKEN". By default, this token
# is set to have permissive access. However, this is not a good practice
@@ -47,17 +47,17 @@ jobs:
builds-and-tests:
strategy:
matrix:
os: [ubuntu-latest, macos-latest, windows-latest]
os: [ubuntu-20.04, macos-10.15, windows-2022]
target: [dbg, opt]
mpi: [seq, par]
build-system: [make, cmake]
hypre-target: [int32]
exclude:
- os: ubuntu-latest
- os: ubuntu-20.04
build-system: cmake
- os: macos-latest
- os: macos-10.15
build-system: cmake
- os: windows-latest
- os: windows-2022
build-system: make
# 'include' allows us to:
# - Add a variable to all jobs without creating a new matrix dimension.
@@ -72,15 +72,15 @@ jobs:
codecov: NO
- target: opt
codecov: YES
- os: windows-latest
- os: windows-2022
codecov: NO
- os: ubuntu-latest
- os: ubuntu-20.04
target: opt
codecov: NO
mpi: par
build-system: cmake
hypre-target: int32
- os: ubuntu-latest
- os: ubuntu-20.04
target: opt
codecov: NO
mpi: par
@@ -112,35 +112,35 @@ jobs:
# TODO: It would be nice to have only one step, e.g. with a dedicated
# action, but I (@adrienbernede) don't see how at the moment.
- name: get MPI (Linux)
if: matrix.mpi == 'par' && matrix.os == 'ubuntu-latest'
if: matrix.mpi == 'par' && matrix.os == 'ubuntu-20.04'
run: |
sudo apt-get install mpich libmpich-dev
export MAKE_CXX_FLAG="MPICXX=mpic++"
- name: get lcov (Linux)
if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-latest'
if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-20.04'
run: |
sudo apt-get install lcov
- name: Set up Homebrew
if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-10.15'
uses: Homebrew/actions/setup-homebrew@c4aafe8c4620bf08883dd4679c374f11e73329d3
- name: get MPI (MacOS)
if: matrix.mpi == 'par' && matrix.os == 'macos-latest'
if: matrix.mpi == 'par' && matrix.os == 'macos-10.15'
run: |
export HOMEBREW_NO_INSTALL_CLEANUP=1
brew install openmpi
export MAKE_CXX_FLAG="MPICXX=mpic++"
- name: get MPI (MacOS)
if: matrix.codecov == 'YES' && matrix.os == 'macos-latest'
if: matrix.codecov == 'YES' && matrix.os == 'macos-10.15'
run: |
export HOMEBREW_NO_INSTALL_CLEANUP=1
brew install lcov
- name: get MPI (Windows)
if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
uses: mpi4py/setup-mpi@v1.0.3
# Get Hypre through cache, or build it.
@@ -154,7 +154,7 @@ jobs:
key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.2
- name: get hypre
if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-latest'
if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-2022'
uses: mfem/github-actions/build-hypre@v2.2
with:
archive: ${{ env.HYPRE_ARCHIVE }}
@@ -163,7 +163,7 @@ jobs:
build-system: make
- name: get hypre (Windows)
if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-latest'
if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-2022'
uses: mfem/github-actions/build-hypre@v2.2
with:
archive: ${{ env.HYPRE_ARCHIVE }}
@@ -175,14 +175,14 @@ jobs:
# Install will only run on cache miss.
- name: cache metis
id: metis-cache
if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
if: matrix.mpi == 'par' && matrix.os != 'windows-2022'
uses: actions/cache@v2
with:
path: ${{ env.METIS_TOP_DIR }}
key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
- name: install metis
if: matrix.mpi == 'par' && matrix.os != 'windows-latest' && steps.metis-cache.outputs.cache-hit != 'true'
if: matrix.mpi == 'par' && matrix.os != 'windows-2022' && steps.metis-cache.outputs.cache-hit != 'true'
uses: mfem/github-actions/build-metis@v2.2
with:
archive: ${{ env.METIS_ARCHIVE }}
@@ -196,16 +196,16 @@ jobs:
key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
- name: prepare binary cache location
if: matrix.os == 'windows-latest' && steps.vcpkg-cache.outputs.cache-hit != 'true'
if: matrix.os == 'windows-2022' && steps.vcpkg-cache.outputs.cache-hit != 'true'
run: |
mkdir -p vcpkg_cache
- name: install metis (Windows)
if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
env:
VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
run: |
vcpkg install metis-mfem --triplet=x64-windows-static --overlay-ports=${{ env.MFEM_TOP_DIR }}/config/vcpkg/ports
vcpkg install metis --triplet=x64-windows-static
# MFEM build and test
- name: build
@@ -248,7 +248,7 @@ jobs:
shell: bash
- name: cmake unit tests (Ubuntu 20.04)
if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-latest'
if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-20.04'
run: |
CTEST_CONFIG="Release"
[[ ${{ matrix.target }} == 'dbg' ]] && CTEST_CONFIG="Debug"
@@ -256,7 +256,7 @@ jobs:
shell: bash
- name: cmake tests
if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-latest'
if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-20.04'
run: |
CTEST_CONFIG="Release"
cd ${{ env.MFEM_TOP_DIR }}/build && ctest --output-on-failure -C ${CTEST_CONFIG}
-71
View File
@@ -1,71 +0,0 @@
# For most projects, this workflow file will not need changing; you simply need
# to commit it to your repository.
#
# You may wish to alter this file to override the set of languages analyzed,
# or to provide custom queries or build logic.
#
# ******** NOTE ********
# We have attempted to detect the languages in your repository. Please check
# the `language` matrix defined below to confirm you have the correct set of
# supported CodeQL languages.
#
name: "Static Analysis"
on:
push:
branches: [ "master", "next"]
pull_request:
# The branches below must be a subset of the branches above
branches: [ "master" ]
jobs:
analyze:
name: Analyze
runs-on: ubuntu-latest
permissions:
actions: read
contents: read
security-events: write
strategy:
fail-fast: false
matrix:
language: [ 'cpp' ]
# CodeQL supports [ 'cpp', 'csharp', 'go', 'java', 'javascript', 'python', 'ruby' ]
# Learn more about CodeQL language support at https://aka.ms/codeql-docs/language-support
steps:
- name: Checkout repository
uses: actions/checkout@v3
# Initializes the CodeQL tools for scanning.
- name: Initialize CodeQL
uses: github/codeql-action/init@v2
with:
languages: ${{ matrix.language }}
# If you wish to specify custom queries, you can do so here or in a config file.
# By default, queries listed here will override any specified in a config file.
# Prefix the list here with "+" to use these queries and those in the config file.
# Details on CodeQL's query packs refer to : https://docs.github.com/en/code-security/code-scanning/automatically-scanning-your-code-for-vulnerabilities-and-errors/configuring-code-scanning#using-queries-in-ql-packs
# queries: security-extended,security-and-quality
queries: lgtm
# Autobuild attempts to build any compiled languages (C/C++, C#, or Java).
# If this step fails, then you should remove it and run the build manually (see below)
- name: Autobuild
uses: github/codeql-action/autobuild@v2
# ️ Command-line programs to run using the OS shell.
# 📚 See https://docs.github.com/en/actions/using-workflows/workflow-syntax-for-github-actions#jobsjob_idstepsrun
# If the Autobuild fails above, remove it and uncomment the following three lines.
# modify them (or add more) to build your code if your project, please refer to the EXAMPLE below for guidance.
# - run: |
# echo "Run, Build Application using script"
# ./location_of_script_within_repo/buildscript.sh
- name: Perform CodeQL Analysis
uses: github/codeql-action/analyze@v2
+2 -2
View File
@@ -9,7 +9,7 @@
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
name: "Build Analysis"
name: build-analysis
permissions:
actions: write
@@ -31,7 +31,7 @@ env:
jobs:
gitignore:
runs-on: ubuntu-latest
runs-on: ubuntu-18.04
steps:
- name: Cancel Previous Runs
+6 -8
View File
@@ -9,7 +9,7 @@
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
name: "Checks"
name: repo-check
permissions:
actions: write
@@ -28,7 +28,7 @@ on:
jobs:
file-headers-check:
runs-on: ubuntu-latest
runs-on: ubuntu-18.04
if: |
(github.event_name == 'push' ||
github.event.pull_request.head.repo.full_name != github.repository)
@@ -78,7 +78,7 @@ jobs:
exit 1
code-style:
runs-on: ubuntu-latest
runs-on: ubuntu-18.04
if: |
(github.event_name == 'push' ||
github.event.pull_request.head.repo.full_name != github.repository)
@@ -88,14 +88,14 @@ jobs:
- name: get astyle
run: |
sudo apt-get install astyle
sudo apt-get install astyle=3.1-1ubuntu2
- name: style check
run: |
./config/githooks/pre-push --style
documentation:
runs-on: ubuntu-latest
runs-on: ubuntu-18.04
if: |
(github.event_name == 'push' ||
github.event.pull_request.head.repo.full_name != github.repository)
@@ -106,8 +106,6 @@ jobs:
- name: get doxygen and graphviz
run: |
sudo apt-get install doxygen graphviz
cd doc
doxygen -u CodeDocumentation.conf.in 2>/dev/null
- name: build documentation
run: |
@@ -120,7 +118,7 @@ jobs:
github.ref != 'refs/heads/master' &&
(github.event_name == 'push' ||
github.event.pull_request.head.repo.full_name != github.repository)
runs-on: ubuntu-latest
runs-on: ubuntu-18.04
steps:
- name: checkout mfem
uses: actions/checkout@v2
-2
View File
@@ -307,8 +307,6 @@ miniapps/solvers/sol.*
miniapps/parelag/MultilevelHcurlHdivSolver
miniapps/parelag/*.mesh
miniapps/hooke/hooke
# Unit test binary and outputs
tests/unit/output_meshes
tests/unit/unit_tests
+16 -72
View File
@@ -10,29 +10,8 @@
Version 4.4.1 (development)
===========================
Meshing improvements
--------------------
- Added support for mixed meshes and pyramids in GSLIB-FindPoints.
Discretization improvements
---------------------------
- Added support for assembling low-order-refined matrices using a GPU-enabled
"batched" algorithm. The lor_solvers and plor_solvers now fully support GPU
acceleration.
- Added support for partial assembly and fully matrix-free operators on mixed
meshes (different element types and p-adaptivity) through libCEED, including
device acceleration, e.g. with NVIDIA and AMD GPUs. The p-adaptivity is
currently limited by MFEM capabilities, i.e. 2D serial meshes. All mixed
element topologies are supported in serial and parallel: segment, triangle,
square, tetrahedron, cube, prism, and pyramid.
- Added full assembly and device support for several LinearForm integrators:
* DomainLF: (f, v)
* VectorDomainLF: ((f1,...,fn), (v1,...,vn))
* DomainLFGrad: (f, grad(v))
* VectorDomainLFGrad: ((f1x,f1y,f1z,...,fnx,fny,fnz), grad(v1,...,vn))
- Added example for body-fitted volumetric and shape integration using the
Algoim library.
- Added WhiteGaussianNoiseDomainLFIntegrator: a LinearFormIntegrator class for
spatial Gaussian white noise.
@@ -40,53 +19,30 @@ Discretization improvements
- Added a new Zienkiewicz-Zhu patch recovery-based a posteriori error estimator.
See fem/estimators.hpp.
Linear and nonlinear solvers
----------------------------
New and updated examples and miniapps
-------------------------------------
- Added a new elasticity miniapp, Hooke, that showcases a low-level approach of
using MFEM to solve a nonlinear elasticity problem based on the fundamental
finite element operator decomposition. The miniapp also integrates with
automatic differentiation tools like a native dual number implementation or a
third party library such as Enzyme. See miniapps/elasticity for more details.
- Add a new example code, Example 33/33p, to demonstrate the solution of
spectral fractional PDEs with MFEM.
Integrations, testing and documentation
---------------------------------------
- Added a Dockerfile for a simple MFEM container, see config/docker/README.md.
- Added support for ParMoonolith, https://bitbucket.org/zulianp/par_moonolith,
which provides parallel non-conforming, non-matching, variational, volumetric
mesh information transfer. With ParMortarAssember, fields can be exchanged
between arbitrarily distributed and unrelated finite element meshes in a
variationally consistent way.
- Added support for the LLVM-based automatic differentiation tool Enzyme, see
https://github.com/EnzymeAD/Enzyme. Build system flags and a convenience
header are provided. The functionality and interaction are demonstrated in a
new miniapp in miniapps/elasticity.
- Added full assembly and device support for several LinearForm integrators:
* DomainLF: (f, v)
* VectorDomainLF: ((f1,...,fn), (v1,...,vn))
* DomainLFGrad: (f, grad(v))
* VectorDomainLFGrad: ((f1x,f1y,f1z,...,fnx,fny,fnz), grad(v1,...,vn))
- Added example for body-fitted volumetric and shape integration using the
Algoim library.
- Add a new example code, Example 33/33p, to demonstrate the solution of
spectral fractional PDEs with MFEM.
- Added a Dockerfile for a simple MFEM container, see config/docker/README.md.
- Added support for assembling low-order-refined matrices using a GPU-enabled
"batched" algorithm. The lor_solvers and plor_solvers now fully support GPU
acceleration.
- Added Windows 2022 CI testing with GitHub actions.
Miscellaneous
-------------
- Various other simplifications, extensions, and bugfixes in the code.
- Added boundary elimination with device support for `SparseMatrix` and
`HypreParMatrix`.
- When using `AssemblyLevel::FULL`, `FABilinearFormExtension::FormSystemMatrix`
outputs an `OperatorHandle` containing a `SparseMatrix` in serial, and an
`HypreParMatrix` in parallel (instead of a `ConstrainedOperator`).
- Added TMOP metrics for mesh untangling and worst-case quality improvement.
- Added support for mixed meshes and pyramids in GSLIB-FindPoints.
Version 4.4, released on March 21, 2022
=======================================
@@ -119,11 +75,6 @@ Meshing improvements
- Added a simpler interface to access mesh face information, see FaceInformation
and GetFaceInformation in the Mesh class.
- Added the method ParMesh::GetSerialMesh() that reconstructs a partitioned
parallel mesh on a given single rank. Also, added the method
ParMesh::PrintAsSerial() that saves the reconstructed serial mesh to a C++
stream on rank 0.
- Gmsh meshes where all elements have zero physical tag (the default Gmsh output
format if no physical groups are defined) are now successfully loaded, and
elements are reassigned attribute number 1.
@@ -226,13 +177,6 @@ Miscellaneous
- Fixed several MinGW build issues on Windows.
- In various places in the library, replace the use of 'long' with 'long long'
to better support Win64 builds where 'long' is 32-bit and 'long long' is
64-bit. On Linux and MacOS, both types are typically 64-bit.
- Update various "MemoryUsage" methods to return 'std::size_t' instead of 'long'
since the latter is 32-bit in Win64 builds.
- Added 'double' atomicAdd implementation for previous versions of CUDA.
- HypreParVector and Vector now support C++ move semantics, and the copy
+2 -9
View File
@@ -136,8 +136,6 @@ if (MFEM_USE_CUDA)
"CUDA flags set for MFEM" FORCE)
set(CUSPARSE_FOUND TRUE)
set(CUSPARSE_LIBRARIES "cusparse")
set(CUBLAS_FOUND TRUE)
set(CUSBLAS_LIBRARIES "cublas")
endif()
if (XSDK_ENABLE_C)
@@ -454,11 +452,6 @@ if (MFEM_USE_PARELAG)
find_package(PARELAG REQUIRED)
endif()
# Enzyme
if (MFEM_USE_ENZYME)
find_package(ENZYME REQUIRED)
endif()
# MFEM_TIMER_TYPE
if (NOT DEFINED MFEM_TIMER_TYPE)
if (APPLE)
@@ -485,8 +478,8 @@ endif()
set(MFEM_TPLS OPENMP HYPRE BLAS LAPACK SuperLUDist METIS SuiteSparse SUNDIALS
PETSC SLEPC MESQUITE MUMPS STRUMPACK AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
ADIOS2 CUBLAS CUSPARSE MKL_CPARDISO AMGX CALIPER CODIPACK BENCHMARK PARELAG
MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM ENZYME)
ADIOS2 CUSPARSE MKL_CPARDISO AMGX CALIPER CODIPACK BENCHMARK PARELAG
MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM)
# Add all *_FOUND libraries in the variable TPL_LIBRARIES.
set(TPL_LIBRARIES "")
-1
View File
@@ -131,7 +131,6 @@ The MFEM source code has the following structure:
│ ├── common
│ ├── electromagnetics
│ ├── gslib
│ ├── hooke
│ ├── meshing
│ ├── mtop
│ ├── navier
+2 -16
View File
@@ -558,14 +558,6 @@ MFEM_USE_PARELAG = YES/NO
use ParELAG. In fact, ParELAG is dependent on MFEM. Therefore, this option
currently only concerns the miniapps.
MFEM_USE_ENZYME = YES/NO
Enables automatic differentiation support through the LLVM plugin Enzyme.
This requires the compiler to be set to clang (>=14.0.0). We also advise to
use the link time optimization (LTO) plugin, to enable functions that you
define over multiple files (compilation units) and want to be differentiated
automatically, to work. This requires to also use LLVM/LLD for linking.
Recommended options are in config/defaults.mk.
MFEM_BUILD_TAG = (any value)
An optional tag to characterize the build. Exported to config/config.mk.
Can be used to identify the MFEM build from other makefiles.
@@ -768,6 +760,8 @@ The specific libraries and their options are:
Options: BLITZ_OPT, BLITZ_LIB
Versions: BLITZ = 1.0.2
- MKL CPardiso (optional), used when MFEM_USE_MKL_CPARDISO = YES.
URL: https://software.intel.com/content/www/us/en/develop/tools/math-kernel-library.html
Options: MKL_CPARDISO_OPT, MKL_CPARDISO_LIB.
@@ -844,12 +838,6 @@ The specific libraries and their options are:
URL: https://github.com/LLNL/parelag
Options: PARELAG_DIR, PARELAG_OPT, PARELAG_LIB.
- Enzyme, used when MFEM_USE_ENZYME = YES. Requires LLVM/Clang >= 14.0.0.
URL: https://github.com/EnzymeAD/Enzyme
Options: ENZYME_DIR, ENZYME_OPT, ENZYME_LIB.
Versions: Enzyme >= v0.0.33.
Building with CMake
===================
The MFEM build system consists of two steps: configuration and compilation.
@@ -988,7 +976,6 @@ MFEM_USE_CALIPER
MFEM_USE_FMS
MFEM_USE_BENCHMARK
MFEM_USE_PARELAG
MFEM_USE_ENZYME
The following options are CMake specific:
@@ -1048,7 +1035,6 @@ The CMake build system adds auto-detection for the following packages/libraries:
- FMS
- BENCHMARK
- ParELAG
- Enzyme
The following built-in CMake packages are also used:
-1
View File
@@ -61,7 +61,6 @@ set(MFEM_USE_CALIPER @MFEM_USE_CALIPER@)
set(MFEM_USE_ALGOIM @MFEM_USE_ALGOIM@)
set(MFEM_USE_BENCHMARK @MFEM_USE_BENCHMARK@)
set(MFEM_USE_PARELAG @MFEM_USE_PARELAG@)
set(MFEM_USE_ENZYME @MFEM_USE_ENZYME@)
set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
set(MFEM_CXX_FLAGS "@CMAKE_CXX_FLAGS@")
-3
View File
@@ -190,7 +190,4 @@
// Enable MFEM functionality based on the Google Benchmark library.
#cmakedefine MFEM_USE_BENCHMARK
// Enable Enzyme for AD
#cmakedefine MFEM_USE_ENZYME
#endif // MFEM_CONFIG_HEADER
-27
View File
@@ -1,27 +0,0 @@
# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
# LICENSE and NOTICE for details. LLNL-CODE-806117.
#
# This file is part of the MFEM library. For more information and source code
# availability visit https://mfem.org.
#
# MFEM is free software; you can redistribute it and/or modify it under the
# terms of the BSD-3 license. We welcome feedback and contributions, see file
# CONTRIBUTING.md for details.
message(STATUS "Looking for ENZYME ...")
message(STATUS " in ENZYME_DIR = ${ENZYME_DIR}")
# Make sure the directory and version combination works. Do nothing otherwise.
if(EXISTS "${ENZYME_DIR}/ClangEnzyme-${ENZYME_VERSION}.so")
message(STATUS "Found ENZYME: ${ENZYME_DIR}/ClangEnzyme-${ENZYME_VERSION}.so")
# Set ENZYME_FOUND
set(ENZYME_FOUND TRUE CACHE BOOL "ENZYME was found." FORCE)
# Set CXX flags to accomodate the Enzyme Clang plugin
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xclang -load -Xclang ${ENZYME_DIR}/ClangEnzyme-${ENZYME_VERSION}.so -mllvm -enzyme-loose-types=1")
set(MFEM_USE_ENZYME YES)
else()
endif()
@@ -894,7 +894,7 @@ function(mfem_export_mk_files)
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_MOONOLITH MFEM_USE_ALGOIM MFEM_USE_ENZYME)
MFEM_USE_MOONOLITH MFEM_USE_ALGOIM)
foreach(var ${CONFIG_MK_BOOL_VARS})
if (${var})
set(${var} YES)
-3
View File
@@ -195,7 +195,4 @@
// Enable functionality based on the Google Benchmark library.
// #define MFEM_USE_BENCHMARK
// Enable the Enzyme LLVM plugin
// #define MFEM_USE_ENZYME
#endif // MFEM_CONFIG_HEADER
-1
View File
@@ -63,7 +63,6 @@ MFEM_USE_ADFORWARD = @MFEM_USE_ADFORWARD@
MFEM_USE_CODIPACK = @MFEM_USE_CODIPACK@
MFEM_USE_BENCHMARK = @MFEM_USE_BENCHMARK@
MFEM_USE_PARELAG = @MFEM_USE_PARELAG@
MFEM_USE_ENZYME = @MFEM_USE_ENZYME@
# Compiler, compile options, and link options
MFEM_CXX = @MFEM_CXX@
-1
View File
@@ -64,7 +64,6 @@ option(MFEM_USE_ADFORWARD "Enable forward mode for AD" OFF)
option(MFEM_USE_CODIPACK "Enable automatic differentiation (AD) using CoDiPack" OFF)
option(MFEM_USE_BENCHMARK "Enable Google Benchmark" OFF)
option(MFEM_USE_PARELAG "Enable ParELAG" OFF)
option(MFEM_USE_ENZYME "Enable Enzyme" OFF)
# Optional overrides for autodetected MPIEXEC and MPIEXEC_NUMPROC_FLAG
# set(MFEM_MPIEXEC "mpirun" CACHE STRING "Command for running MPI tests")
+1 -21
View File
@@ -42,9 +42,6 @@ STATIC = YES
SHARED = NO
# CUDA configuration options
#
# If you set MFEM_USE_ENZYME=YES, CUDA_CXX has to be configured to use cuda with
# clang as its host compiler.
CUDA_CXX = nvcc
CUDA_ARCH = sm_60
CUDA_FLAGS = -x=cu --expt-extended-lambda -arch=$(CUDA_ARCH)
@@ -166,7 +163,6 @@ MFEM_USE_ADFORWARD = NO
MFEM_USE_CODIPACK = NO
MFEM_USE_BENCHMARK = NO
MFEM_USE_PARELAG = NO
MFEM_USE_ENZYME = NO
# MPI library compile and link flags
# These settings are used only when building MFEM with MPI + HIP
@@ -207,7 +203,7 @@ HYPRE_OPT = -I$(HYPRE_DIR)/include
HYPRE_LIB = -L$(HYPRE_DIR)/lib -lHYPRE
ifeq (YES,$(MFEM_USE_CUDA))
# This is only necessary when hypre is built with cuda:
HYPRE_LIB += -lcusparse -lcurand -lcublas
HYPRE_LIB += -lcusparse -lcurand
endif
ifeq (YES,$(MFEM_USE_HIP))
# This is only necessary when hypre is built with hip:
@@ -524,22 +520,6 @@ PARELAG_DIR = @MFEM_DIR@/../parelag
PARELAG_OPT = -I$(PARELAG_DIR)/src -I$(PARELAG_DIR)/build/src
PARELAG_LIB = -L$(PARELAG_DIR)/build/src -lParELAG
# Enzyme configuration
# If you want to enable automatic differentiation at compile time, use the
# options below, adapted to your configuration. To be more flexible, we
# recommend using the Enzyme plugin during link time optimization. One option is
# to add your options to the global compiler/linker flags like
#
# BASE_FLAGS += -flto
# CXX_XLINKER += -fuse-ld=lld -Wl,--lto-legacy-pass-manager\
# -Wl,-mllvm=-load=$(ENZYME_DIR)/LLDEnzyme-$(ENZYME_VERSION).so -Wl,
#
ENZYME_DIR ?= @MFEM_DIR@/../enzyme
ENZYME_VERSION ?= 14
ENZYME_OPT = -fno-experimental-new-pass-manager -Xclang -load -Xclang $(ENZYME_DIR)/ClangEnzyme-$(ENZYME_VERSION).so
ENZYME_LIB = ""
# If YES, enable some informational messages
VERBOSE = NO
+22 -19
View File
@@ -1,27 +1,30 @@
FROM ghcr.io/mfem/mfem-ubuntu-base:latest as builder
FROM ghcr.io/rse-ops/cuda-ubuntu-20.04:cuda-11.0.3
# docker build -t ghcr.io/mfem/mfem-ubuntu .
COPY ./config/docker/spack.yaml /opt/mfem-env/spack.yaml
RUN apt-get install -y python3 && \
cd /opt/mfem-env && \
. /opt/spack/share/spack/setup-env.sh && \
spack env activate . && \
spack env view regenerate
FROM ubuntu:22.04
COPY --from=builder /opt/view /opt/view
COPY --from=builder /opt/mfem-view /opt/mfem-view
# docker build -t ghcr.io/mfem/mfem-ubuntu-base .
RUN apt-get update && \
apt-get install -y unzip gfortran && \
spack compiler find && \
apt-get install -y libcurl4-openssl-dev libssl-dev
ENV PATH=$PATH:/opt/mfem-view/bin
ENV LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/opt/mfem-view/lib:/opt/mfem-view/lib64
ENV DEBIAN_FRONTEND=noninteractive
# /code is the working directory for code
WORKDIR /code
COPY . /code
# This is for a spack environment/view to install from there
WORKDIR /opt/mfem-env
RUN . /opt/spack/share/spack/setup-env.sh && \
spack env create -d . && \
echo " concretization: together" >> spack.yaml && \
spack env activate . && \
spack develop --path /code mfem@master+examples+miniapps && \
spack add mfem@master+examples+miniapps && \
spack install
# ensure mfem always on various paths
RUN cd /opt/mfem-env && \
spack env activate --sh -d . >> /etc/profile.d/z10_spack_environment.sh
# The user will see the view on shell into the container
WORKDIR /opt/mfem-view
ENTRYPOINT ["/bin/bash"]
WORKDIR /opt/mfem-env/.spack-env/view/
ENTRYPOINT ["/bin/bash", "--rcfile", "/etc/profile", "-l", "-c"]
-47
View File
@@ -1,47 +0,0 @@
FROM ghcr.io/rse-ops/cuda-ubuntu-20.04:cuda-11.0.3
# docker build -f Dockerfile.base -t ghcr.io/mfem/mfem-ubuntu-base .
RUN apt-get update && \
apt-get install -y unzip gfortran && \
spack compiler find && \
apt-get install -y libcurl4-openssl-dev libssl-dev
# /code is the working directory for code
WORKDIR /code
COPY . /code
# This is for a spack environment/view to install from there
RUN mkdir -p /opt/mfem-env \
&& (echo "spack:" \
&& echo " view:" \
&& echo " mfem:" \
&& echo " root: /opt/mfem-view" \
&& echo " link_type: copy" \
&& echo " packages:" \
&& echo " all:" \
&& echo " target:" \
&& echo " - x86_64_v3" \
&& echo " config:" \
&& echo " concretizer: clingo" \
&& echo " compiler:" \
&& echo " target:" \
&& echo " - x86_64_v3" \
&& echo " install_missing_compilers: true" \
&& echo " concretization: together") > /opt/mfem-env/spack.yaml
RUN cd /opt/mfem-env && \
. /opt/spack/share/spack/setup-env.sh && \
spack env activate . && \
spack develop --path /code mfem@master+examples+miniapps && \
spack add mfem@master+examples+miniapps # && \
# spack install
# ensure mfem always on various paths
#RUN cd /opt/mfem-env && \
# spack env activate --sh -d . >> /etc/profile.d/z10_spack_environment.sh
# Present the software install when we shell in
# The view is at /opt/mfem-env/.spack-env/view
#WORKDIR /opt/software
#ENTRYPOINT ["/bin/bash", "--rcfile", "/etc/profile", "-l", "-c"]
+7 -24
View File
@@ -1,8 +1,7 @@
# mfem Docker
We provide a [Dockerfile.base](Dockerfile.base) to build an ubuntu base image,
and a [Dockerfile](Dockerfile) to build a smaller one with a multi-stage build.
You can use this image for a demo of using mfem! 🎉️
We provide a [Dockerfile](Dockerfile) to build an ubuntu base image. You can use
this image for a demo of using mfem! 🎉️
Updated containers are built and deployed on merges to the main branch and releases.
If you want to request a build on demand, you can [manually run the workflow](https://docs.github.com/en/actions/managing-workflow-runs/manually-running-a-workflow) thanks to the workflow dispatch event.
@@ -15,33 +14,18 @@ is the [GitHub packages](https://github.com/features/packages) registry that sup
Docker images and other OCI artifacts. From the root of the repository:
```bash
$ docker build -f config/docker/Dockerfile -t ghcr.io/mfem/mfem-ubuntu .
$ docker build -f config/docker/Dockerfile.base -t ghcr.io/mfem/mfem-ubuntu-base .
$ docker build -f config/docker/Dockerfile -t ghcr.io/mfem/mfem-ubuntu-base .
```
### Shell Ubuntu
To shell into the container:
or this directory:
```bash
$ docker run -it ghcr.io/mfem/mfem-ubuntu
$ docker build -f Dockerfile -t ghcr.io/mfem/mfem-ubuntu-base ../../
```
This smaller image has a view where everything is installed.
### Shell
```bash
$ ls
bin etc include lib libexec sbin share var
```
- Examples are in share/mfem/examples
- Examples are in share/mfem/miniapps
You can read more about interaction with these examples and miniapps below.
### Shell Ubuntu Base
To shell into the container:
To shell into a container (here is an example with ubuntu):
```bash
$ docker run -it ghcr.io/mfem/mfem-ubuntu-base bash
@@ -144,4 +128,3 @@ $ docker run -it ghcr.io/mfem/mfem-ubuntu-base -v $PWD:/src bash
In the above, we can pretend your project is in the present working directory (PWD) and we are
binding to source. You can then use the mfem in the container for development, and if you
want to distribute your library or app in a container, you can use the mfem container as the base.
-11
View File
@@ -1,11 +0,0 @@
spack:
specs: [mfem@master+examples+miniapps]
view:
mfem:
root: /opt/mfem-view
link_type: copy
concretization: together
develop:
mfem:
path: /code
spec: mfem@master+examples+miniapps
@@ -1,8 +0,0 @@
--- a/CMakeLists.txt Wed Dec 21 18:24:22 2016
+++ b/CMakeLists.txt Wed Dec 21 18:24:26 2016
@@ -20,4 +20,4 @@
# Recursively look for CMakeLists.txt in subdirs.
add_subdirectory("include")
add_subdirectory("libmetis")
-add_subdirectory("programs")
+# add_subdirectory("programs")
@@ -1,15 +0,0 @@
--- a/CMakeLists.txt Sat Mar 30 17:24:45 2013
+++ b/CMakeLists.txt Wed Dec 21 18:23:43 2016
@@ -4,11 +4,7 @@
set(GKLIB_PATH "GKlib" CACHE PATH "path to GKlib")
set(SHARED FALSE CACHE BOOL "build a shared library")
-if(MSVC)
- set(METIS_INSTALL FALSE)
-else()
- set(METIS_INSTALL TRUE)
-endif()
+set(METIS_INSTALL TRUE)
# Configure libmetis library.
if(SHARED)
@@ -1,34 +0,0 @@
diff --git a/include/metis.h b/include/metis.h
index dc5406a..7732437 100644
--- a/include/metis.h
+++ b/include/metis.h
@@ -72,10 +72,14 @@ typedef __int64 int64_t;
#define PRId64 "I64d"
#define SCNd32 "ld"
#define SCNd64 "I64d"
+#ifdef _WIN32
+#include <stdint.h>
+#else
#define INT32_MIN ((int32_t)_I32_MIN)
#define INT32_MAX _I32_MAX
#define INT64_MIN ((int64_t)_I64_MIN)
#define INT64_MAX _I64_MAX
+#endif
#else
#include <inttypes.h>
#endif
diff --git a/GKlib/gk_arch.h b/GKlib/gk_arch.h
index 78b1431..7258763 100644
--- a/GKlib/gk_arch.h
+++ b/GKlib/gk_arch.h
@@ -32,8 +32,8 @@
#ifdef __MSC__
- #include "ms_stdint.h"
- #include "ms_inttypes.h"
+ #include <stdint.h>
+ #include <inttypes.h>
#include "ms_stat.h"
#else
#ifndef SUNOS
@@ -1,11 +0,0 @@
--- a/GKlib/gk_arch.h Wed Dec 21 18:34:18 2016
+++ b/GKlib/gk_arch.h Wed Dec 21 18:30:49 2016
@@ -58,7 +58,7 @@
#define PTRDIFF_MAX INT64_MAX
#endif
-#ifdef __MSC__
+#if defined(__MSC__) && (_MSC_VER < 1900)
/* MSC does not have rint() function */
#define rint(x) ((int)((x)+0.5))
@@ -1,14 +0,0 @@
diff --git a/CMakeLists.txt b/CMakeLists.txt
index e94f050..b9613a7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,7 +1,8 @@
cmake_minimum_required(VERSION 2.8)
project(METIS)
-set(GKLIB_PATH "GKlib" CACHE PATH "path to GKlib")
+set(GKLIB_PATH "${CMAKE_SOURCE_DIR}/GKlib" CACHE PATH "path to GKlib")
+
set(SHARED FALSE CACHE BOOL "build a shared library")
set(METIS_INSTALL TRUE)
@@ -1,11 +0,0 @@
--- a/libmetis/metislib.h Sat Mar 30 17:24:45 2013
+++ b/libmetis/metislib.h Wed Dec 21 18:30:59 2016
@@ -31,7 +31,7 @@
#include <proto.h>
-#if defined(COMPILER_MSC)
+#if defined(COMPILER_MSC) && (_MSC_VER < 1900)
#if defined(rint)
#undef rint
#endif
@@ -1,10 +0,0 @@
--- a/libmetis/CMakeLists.txt Sat Mar 30 17:24:45 2013
+++ b/libmetis/CMakeLists.txt Wed Dec 21 17:41:37 2016
@@ -11,6 +11,6 @@
if(METIS_INSTALL)
install(TARGETS metis
LIBRARY DESTINATION lib
- RUNTIME DESTINATION lib
+ RUNTIME DESTINATION bin
ARCHIVE DESTINATION lib)
endif()
@@ -1,44 +0,0 @@
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b9613a7..e43ffee 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -22,3 +22,23 @@ include_directories(include)
add_subdirectory("include")
add_subdirectory("libmetis")
# add_subdirectory("programs")
+
+if(METIS_INSTALL)
+ set(PRJ_NAME metis)
+ set(PRJ_VER 5.1.0)
+ install(EXPORT metisTargets
+ FILE ${PRJ_NAME}Targets.cmake
+ DESTINATION lib/cmake/${PRJ_NAME})
+ include(CMakePackageConfigHelpers)
+ write_basic_package_version_file(
+ ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}ConfigVersion.cmake
+ VERSION ${PRJ_VER}
+ COMPATIBILITY SameMajorVersion)
+ file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}Config.cmake
+ "include(\${CMAKE_CURRENT_LIST_DIR}/${PRJ_NAME}Targets.cmake)")
+ install(FILES
+ ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}ConfigVersion.cmake
+ ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}Config.cmake
+ DESTINATION lib/cmake/${PRJ_NAME})
+endif()
+
diff --git a/libmetis/CMakeLists.txt b/libmetis/CMakeLists.txt
index 7a5fc74..5a68cf0 100644
--- a/libmetis/CMakeLists.txt
+++ b/libmetis/CMakeLists.txt
@@ -9,8 +9,9 @@ if(UNIX)
endif()
if(METIS_INSTALL)
- install(TARGETS metis
+ install(TARGETS metis EXPORT metisTargets
LIBRARY DESTINATION lib
RUNTIME DESTINATION bin
- ARCHIVE DESTINATION lib)
+ ARCHIVE DESTINATION lib
+ INCLUDES DESTINATION include)
endif()
@@ -1,41 +0,0 @@
vcpkg_check_linkage(ONLY_STATIC_LIBRARY)
set(OPTIONS -DSHARED=OFF)
set(METIS_VERSION 5.1.0)
vcpkg_download_distfile(ARCHIVE
URLS "https://github.com/mfem/tpls/raw/gh-pages/metis-${METIS_VERSION}.tar.gz"
FILENAME "metis-${METIS_VERSION}.tar.gz"
SHA512 deea47749d13bd06fbeaf98a53c6c0b61603ddc17a43dae81d72c8015576f6495fd83c11b0ef68d024879ed5415c14ebdbd87ce49c181bdac680573bea8bdb25
)
vcpkg_extract_source_archive_ex(
OUT_SOURCE_PATH SOURCE_PATH
ARCHIVE ${ARCHIVE}
REF ${METIS_VERSION}
PATCHES
enable-install.patch
disable-programs.patch
fix-runtime-install-destination.patch
fix-metis-vs14-math.patch
fix-gklib-vs14-math.patch
fix-linux-build-error.patch
install-metisConfig.patch
fix-INT_MIN_define.patch
)
vcpkg_configure_cmake(
SOURCE_PATH ${SOURCE_PATH}
PREFER_NINJA
OPTIONS ${OPTIONS}
)
vcpkg_install_cmake()
vcpkg_copy_pdbs()
vcpkg_fixup_cmake_targets(CONFIG_PATH lib/cmake/metis)
file(REMOVE_RECURSE ${CURRENT_PACKAGES_DIR}/debug/include)
# Handle copyright
file(COPY ${SOURCE_PATH}/LICENSE.txt DESTINATION ${CURRENT_PACKAGES_DIR}/share/metis)
file(INSTALL ${SOURCE_PATH}/LICENSE.txt DESTINATION ${CURRENT_PACKAGES_DIR}/share/${PORT} RENAME copyright)
-7
View File
@@ -1,7 +0,0 @@
{
"name": "metis-mfem",
"version-string": "5.1.0",
"port-version": 0,
"description": "Serial Graph Partitioning and Fill-reducing Matrix Ordering",
"homepage": "https://glaros.dtc.umn.edu/gkhome/metis/metis/overview"
}
+1 -1
View File
@@ -2349,7 +2349,7 @@ PLANTUML_INCLUDE_PATH =
# Minimum value: 0, maximum value: 10000, default value: 50.
# This tag requires that the tag HAVE_DOT is set to YES.
DOT_GRAPH_MAX_NODES = 100
DOT_GRAPH_MAX_NODES = 50
# The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the graphs
# generated by dot. A depth value of 3 means that only nodes reachable from the
-9
View File
@@ -30,7 +30,6 @@
//
// Device sample runs:
// ex1 -pa -d cuda
// * ex1 -fa -d cuda
// ex1 -pa -d raja-cuda
// * ex1 -pa -d raja-hip
// ex1 -pa -d occa-cuda
@@ -38,13 +37,9 @@
// ex1 -pa -d occa-omp
// ex1 -pa -d ceed-cpu
// ex1 -pa -d ceed-cpu -o 4 -a
// ex1 -pa -d ceed-cpu -m ../data/square-mixed.mesh
// ex1 -pa -d ceed-cpu -m ../data/fichera-mixed.mesh
// * ex1 -pa -d ceed-cuda
// * ex1 -pa -d ceed-hip
// ex1 -pa -d ceed-cuda:/gpu/cuda/shared
// ex1 -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
// ex1 -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
// ex1 -m ../data/beam-hex.mesh -pa -d cuda
// ex1 -m ../data/beam-tet.mesh -pa -d ceed-cpu
// ex1 -m ../data/beam-tet.mesh -pa -d ceed-cuda:/gpu/cuda/ref
@@ -78,7 +73,6 @@ int main(int argc, char *argv[])
int order = 1;
bool static_cond = false;
bool pa = false;
bool fa = false;
const char *device_config = "cpu";
bool visualization = true;
bool algebraic_ceed = false;
@@ -93,8 +87,6 @@ int main(int argc, char *argv[])
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
"--no-partial-assembly", "Enable Partial Assembly.");
args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
"--no-full-assembly", "Enable Full Assembly.");
args.AddOption(&device_config, "-d", "--device",
"Device configuration string, see Device::Configure().");
#ifdef MFEM_USE_CEED
@@ -192,7 +184,6 @@ int main(int argc, char *argv[])
// domain integrator.
BilinearForm a(&fespace);
if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
if (fa) { a.SetAssemblyLevel(AssemblyLevel::FULL); }
a.AddDomainIntegrator(new DiffusionIntegrator(one));
// 10. Assemble the bilinear form and the corresponding linear system,
-9
View File
@@ -30,18 +30,13 @@
//
// Device sample runs:
// mpirun -np 4 ex1p -pa -d cuda
// * mpirun -np 4 ex1p -fa -d cuda
// mpirun -np 4 ex1p -pa -d occa-cuda
// mpirun -np 4 ex1p -pa -d raja-omp
// mpirun -np 4 ex1p -pa -d ceed-cpu
// mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
// mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/square-mixed.mesh
// mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/fichera-mixed.mesh
// * mpirun -np 4 ex1p -pa -d ceed-cuda
// * mpirun -np 4 ex1p -pa -d ceed-hip
// mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
// mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
// mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
// mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
//
// Description: This example code demonstrates the use of MFEM to define a
@@ -79,7 +74,6 @@ int main(int argc, char *argv[])
int order = 1;
bool static_cond = false;
bool pa = false;
bool fa = false;
const char *device_config = "cpu";
bool visualization = true;
bool algebraic_ceed = false;
@@ -94,8 +88,6 @@ int main(int argc, char *argv[])
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
"--no-partial-assembly", "Enable Partial Assembly.");
args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
"--no-full-assembly", "Enable Full Assembly.");
args.AddOption(&device_config, "-d", "--device",
"Device configuration string, see Device::Configure().");
#ifdef MFEM_USE_CEED
@@ -219,7 +211,6 @@ int main(int argc, char *argv[])
// Diffusion domain integrator.
ParBilinearForm a(&fespace);
if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
if (fa) { a.SetAssemblyLevel(AssemblyLevel::FULL); }
a.AddDomainIntegrator(new DiffusionIntegrator(one));
// 12. Assemble the parallel bilinear form and the corresponding linear
+1 -1
View File
@@ -182,7 +182,7 @@ int main(int argc, char *argv[])
}
for (int level = 0; level < order_refinements; ++level)
{
collections.Append(new H1_FECollection((int)std::pow(2, level+1), dim));
collections.Append(new H1_FECollection(std::pow(2, level+1), dim));
fespaces.AddOrderRefinedLevel(collections.Last());
}
+1 -1
View File
@@ -219,7 +219,7 @@ int main(int argc, char *argv[])
}
for (int level = 0; level < order_refinements; ++level)
{
collections.Append(new H1_FECollection((int)std::pow(2, level+1), dim));
collections.Append(new H1_FECollection(std::pow(2, level+1), dim));
fespaces->AddOrderRefinedLevel(collections.Last());
}
+133
View File
@@ -100,6 +100,7 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
// Copy the pointers to the integrators
domain_integs = bf->domain_integs;
domain_integs_marker = bf->domain_integs_marker;
boundary_integs = bf->boundary_integs;
boundary_integs_marker = bf->boundary_integs_marker;
@@ -112,6 +113,138 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
AllocMat();
}
BilinearForm::BilinearForm(BilinearForm &&other)
: Matrix(other.fes->GetVSize()), mat(other.mat), mat_e(other.mat_e),
fes(other.fes), assembly(other.assembly), batch(other.batch),
ext(other.ext), sequence(other.sequence), extern_bfs(other.extern_bfs),
element_matrices(other.element_matrices), static_cond(other.static_cond),
hybridization(other.hybridization), diag_policy(other.diag_policy),
precompute_sparsity(other.precompute_sparsity)
{
// We swap stored integrators and markers with the moved nonlinear form
mfem::Swap(domain_integs, other.domain_integs);
mfem::Swap(domain_integs_marker, other.domain_integs_marker);
mfem::Swap(boundary_integs, other.boundary_integs);
mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
mfem::Swap(interior_face_integs, other.interior_face_integs);
mfem::Swap(boundary_face_integs, other.boundary_face_integs);
mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
/// Leave the moved nonlinear form in a state as if it was just constructed
/// with fes
other.sequence = fes->GetSequence();
other.mat = nullptr;
other.mat_e = nullptr;
other.extern_bfs = 0;
other.element_matrices = nullptr;
other.static_cond = nullptr;
other.hybridization = nullptr;
other.precompute_sparsity = 0;
other.diag_policy = DIAG_KEEP;
other.assembly = AssemblyLevel::LEGACY;
other.batch = 1;
other.ext = nullptr;
}
BilinearForm& BilinearForm::operator=(BilinearForm &&other)
{
if (this != &other)
{
/// Cleanup current bilinear form first
delete mat_e;
delete mat;
delete element_matrices;
delete static_cond;
delete hybridization;
for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
for (int i = 0; i < boundary_integs.Size(); i++) { delete boundary_integs[i]; }
for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
delete ext;
/// Null out all our integs and set size of their arrays to zero
for (int k = 0; k < domain_integs.Size(); k++)
{
domain_integs[k] = nullptr;
}
domain_integs.SetSize(0);
for (int k = 0; k < boundary_integs.Size(); k++)
{
boundary_integs[k] = nullptr;
}
boundary_integs.SetSize(0);
for (int k = 0; k < interior_face_integs.Size(); k++)
{
interior_face_integs[k] = nullptr;
}
interior_face_integs.SetSize(0);
for (int k = 0; k < boundary_face_integs.Size(); ++k)
{
boundary_face_integs[k] = nullptr;
}
boundary_face_integs.SetSize(0);
/// Null out all our markers and set size of their arrays to zero
for (int k = 0; k < domain_integs_marker.Size(); ++k)
{
domain_integs_marker[k] = nullptr;
}
domain_integs_marker.SetSize(0);
for (int k = 0; k < boundary_integs_marker.Size(); ++k)
{
boundary_integs_marker[k] = nullptr;
}
boundary_integs_marker.SetSize(0);
for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
{
boundary_face_integs_marker[k] = nullptr;
}
boundary_face_integs_marker.SetSize(0);
/// Now steal data from other bilinear form leaving it in a state as if
/// it was just constructed with fes
Matrix::operator=(std::move(other));
mat = other.mat;
other.mat = nullptr;
mat_e = other.mat_e;
other.mat_e = nullptr;
fes = other.fes;
assembly = other.assembly;
other.assembly = AssemblyLevel::LEGACY;
batch = other.batch;
other.batch = 1;
ext = other.ext;
other.ext = nullptr;
sequence = other.sequence;
other.sequence = fes->GetSequence();
extern_bfs = other.extern_bfs;
other.extern_bfs = 0;
// Swap our empty integ and marker arrays with the moved bilinear form
mfem::Swap(domain_integs, other.domain_integs);
mfem::Swap(domain_integs_marker, other.domain_integs_marker);
mfem::Swap(boundary_integs, other.boundary_integs);
mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
mfem::Swap(interior_face_integs, other.interior_face_integs);
mfem::Swap(boundary_face_integs, other.boundary_face_integs);
mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
element_matrices = other.element_matrices;
other.element_matrices = nullptr;
static_cond = other.static_cond;
other.static_cond = nullptr;
hybridization = other.hybridization;
other.hybridization = nullptr;
diag_policy = other.diag_policy;
other.diag_policy = DIAG_KEEP;
precompute_sparsity = other.precompute_sparsity;
other.precompute_sparsity = 0;
}
return *this;
}
void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
{
if (ext)
+19 -20
View File
@@ -90,23 +90,23 @@ protected:
int extern_bfs;
/// Set of Domain Integrators to be applied.
Array<BilinearFormIntegrator*> domain_integs;
Array<BilinearFormIntegrator*> domain_integs; // owned
/// Element attribute marker (should be of length mesh->attributes.Max() or
/// 0 if mesh->attributes is empty)
/// Includes all by default.
/// 0 - ignore attribute
/// 1 - include attribute
Array<Array<int>*> domain_integs_marker;
Array<Array<int>*> domain_integs_marker; // not owned
/// Set of Boundary Integrators to be applied.
Array<BilinearFormIntegrator*> boundary_integs;
Array<BilinearFormIntegrator*> boundary_integs; // owned
Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned.
/// Set of interior face Integrators to be applied.
Array<BilinearFormIntegrator*> interior_face_integs;
Array<BilinearFormIntegrator*> interior_face_integs; // owned
/// Set of boundary face Integrators to be applied.
Array<BilinearFormIntegrator*> boundary_face_integs;
Array<BilinearFormIntegrator*> boundary_face_integs; // owned
Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned.
DenseMatrix elemmat;
@@ -141,13 +141,6 @@ protected:
ext = NULL;
}
private:
/// Copy construction is not supported; body is undefined.
BilinearForm(const BilinearForm &);
/// Copy assignment is not supported; body is undefined.
BilinearForm &operator=(const BilinearForm &);
public:
/// Creates bilinear form associated with FE space @a *f.
/** The pointer @a f is not owned by the newly constructed object. */
@@ -165,6 +158,20 @@ public:
#precompute_sparsity, see UsePrecomputedSparsity() for details. */
BilinearForm(FiniteElementSpace *f, BilinearForm *bf, int ps = 0);
/// Explicitly prohibit copy construction/assignment of BilinearForms
BilinearForm(const BilinearForm &) = delete;
BilinearForm &operator=(const BilinearForm &) = delete;
/// Move constructor for BilinearForm.
/** This constructor "steals" the owned data members from the @a other
BilinearForm. */
BilinearForm(BilinearForm &&other);
/// Move assignment operator for BilinearForm
/** This assignment first frees all owned data, then "steals" the owned data
members from the @a other BilinearForm. */
BilinearForm& operator=(BilinearForm &&other);
/// Get the size of the BilinearForm as a square matrix.
int Size() const { return height; }
@@ -438,14 +445,6 @@ public:
virtual const Operator *GetOutputRestriction() const
{ return GetRestriction(); }
/// @brief Compute serial RAP operator and store it in @a A as a SparseMatrix.
void SerialRAP(OperatorHandle &A)
{
MFEM_ASSERT(mat, "SerialRAP requires the SparseMatrix to be assembled.");
ConformingAssemble();
A.Reset(mat, false);
}
/** @brief Form the linear system A X = B, corresponding to this bilinear
form and the linear form @a b(.). */
/** This method applies any necessary transformations to the linear system
-52
View File
@@ -251,7 +251,6 @@ PABilinearFormExtension::PABilinearFormExtension(BilinearForm *form)
void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
{
if ( Device::Allows(Backend::CEED_MASK) ) { return; }
ElementDofOrdering ordering = UsesTensorBasis(*a->FESpace())?
ElementDofOrdering::LEXICOGRAPHIC:
ElementDofOrdering::NATIVE;
@@ -957,57 +956,6 @@ void FABilinearFormExtension::Assemble()
}
}
void FABilinearFormExtension::RAP(OperatorHandle &A)
{
#ifdef MFEM_USE_MPI
if ( auto pa = dynamic_cast<ParBilinearForm*>(a) )
{
pa->ParallelRAP(*pa->mat, A);
}
else
#endif
{
a->SerialRAP(A);
}
}
void FABilinearFormExtension::EliminateBC(const Array<int> &ess_dofs,
OperatorHandle &A)
{
#ifdef MFEM_USE_MPI
if ( dynamic_cast<ParBilinearForm*>(a) )
{
A.As<HypreParMatrix>()->EliminateBC(ess_dofs,
DiagonalPolicy::DIAG_ONE);
}
else
#endif
{
A.As<SparseMatrix>()->EliminateBC(ess_dofs,
DiagonalPolicy::DIAG_ONE);
}
}
void FABilinearFormExtension::FormSystemMatrix(const Array<int> &ess_dofs,
OperatorHandle &A)
{
RAP(A);
EliminateBC(ess_dofs, A);
}
void FABilinearFormExtension::FormLinearSystem(const Array<int> &ess_tdof_list,
Vector &x, Vector &b,
OperatorHandle &A,
Vector &X, Vector &B,
int copy_interior)
{
Operator *A_out;
Operator::FormLinearSystem(ess_tdof_list, x, b, A_out, X, B, copy_interior);
delete A_out;
FormSystemMatrix(ess_tdof_list, A);
}
void FABilinearFormExtension::DGMult(const Vector &x, Vector &y) const
{
#ifdef MFEM_USE_MPI
-9
View File
@@ -125,15 +125,6 @@ public:
FABilinearFormExtension(BilinearForm *form);
void Assemble();
void RAP(OperatorHandle &A);
/** @note Always does `DIAG_ONE` policy to be consistent with
`Operator::FormConstrainedSystemOperator`. */
void EliminateBC(const Array<int> &ess_dofs, OperatorHandle &A);
void FormSystemMatrix(const Array<int> &ess_tdof_list, OperatorHandle &A);
void FormLinearSystem(const Array<int> &ess_tdof_list,
Vector &x, Vector &b,
OperatorHandle &A, Vector &X, Vector &B,
int copy_interior = 0);
void Mult(const Vector &x, Vector &y) const;
void MultTranspose(const Vector &x, Vector &y) const;
+1 -2
View File
@@ -2737,8 +2737,7 @@ private:
public:
DivDivIntegrator() { Q = NULL; }
DivDivIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) :
BilinearFormIntegrator(ir), Q(&q) { }
DivDivIntegrator(Coefficient &q) : Q(&q) { }
virtual void AssembleElementMatrix(const FiniteElement &el,
ElementTransformation &Trans,
+1 -10
View File
@@ -30,16 +30,7 @@ void ConvectionIntegrator::AssembleMF(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedMFConvectionIntegrator(*this, fes, Q, alpha);
}
else
{
ceedOp = new ceed::MFConvectionIntegrator(fes, *ir, Q, alpha);
}
ceedOp = new ceed::MFConvectionIntegrator(fes, *ir, Q, alpha);
return;
}
MFEM_ABORT("Error: ConvectionIntegrator::AssembleMF only implemented with"
+1 -10
View File
@@ -1386,16 +1386,7 @@ void ConvectionIntegrator::AssemblePA(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedPAConvectionIntegrator(*this, fes, Q, alpha);
}
else
{
ceedOp = new ceed::PAConvectionIntegrator(fes, *ir, Q, alpha);
}
ceedOp = new ceed::PAConvectionIntegrator(fes, *ir, Q, alpha);
return;
}
const int dims = el.GetDim();
+1 -10
View File
@@ -33,16 +33,7 @@ void DiffusionIntegrator::AssembleMF(const FiniteElementSpace &fes)
MFEM_VERIFY(!VQ && !MQ,
"Only scalar coefficient supported for DiffusionIntegrator"
" with libCEED");
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedMFDiffusionIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::MFDiffusionIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::MFDiffusionIntegrator(fes, *ir, Q);
return;
}
MFEM_ABORT("Error: DiffusionIntegrator::AssembleMF only implemented with"
+4 -16
View File
@@ -271,21 +271,18 @@ void PADiffusionSetup3D(const int Q1D,
D(qx,qy,qz,1,e) = D12; // 1,2
D(qx,qy,qz,2,e) = w_detJ * (A11*R13 + A12*R23 + A13*R33); // 1,3
const double D21 = w_detJ * (A21*R11 + A22*R21 + A23*R31);
const double D22 = w_detJ * (A21*R12 + A22*R22 + A23*R32);
const double D23 = w_detJ * (A21*R13 + A22*R23 + A23*R33);
const double D33 = w_detJ * (A31*R13 + A32*R23 + A33*R33);
D(qx,qy,qz,3,e) = symmetric ? D22 : D21; // 2,2 or 2,1
D(qx,qy,qz,4,e) = symmetric ? D23 : D22; // 2,3 or 2,2
D(qx,qy,qz,5,e) = symmetric ? D33 : D23; // 3,3 or 2,3
if (symmetric)
if (!symmetric)
{
D(qx,qy,qz,3,e) = D22; // 2,2
}
else
{
D(qx,qy,qz,3,e) = w_detJ * (A21*R11 + A22*R21 + A23*R31); // 2,1
D(qx,qy,qz,6,e) = w_detJ * (A31*R11 + A32*R21 + A33*R31); // 3,1
D(qx,qy,qz,7,e) = w_detJ * (A31*R12 + A32*R22 + A33*R32); // 3,2
D(qx,qy,qz,8,e) = D33; // 3,3
@@ -368,16 +365,7 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
MFEM_VERIFY(!VQ && !MQ,
"Only scalar coefficient supported for DiffusionIntegrator"
" with libCEED");
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedPADiffusionIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
return;
}
const int dims = el.GetDim();
+31 -626
View File
@@ -24,20 +24,18 @@ namespace mfem
// PA H(div) Mass Assemble 2D kernel
void PAHdivSetup2D(const int Q1D,
const int coeffDim,
const int NE,
const Array<double> &w,
const Vector &j,
Vector &coeff_,
Vector &op)
{
const bool symmetric = (coeffDim != 4);
const int NQ = Q1D*Q1D;
auto W = w.Read();
auto J = Reshape(j.Read(), NQ, 2, 2, NE);
auto C = Reshape(coeff_.Read(), coeffDim, NQ, NE);
auto y = Reshape(op.Write(), NQ, symmetric ? 3 : 4, NE);
auto coeff = Reshape(coeff_.Read(), NQ, NE);
auto y = Reshape(op.Write(), NQ, 3, NE);
MFEM_FORALL(e, NE,
{
@@ -47,60 +45,28 @@ void PAHdivSetup2D(const int Q1D,
const double J21 = J(q,1,0,e);
const double J12 = J(q,0,1,e);
const double J22 = J(q,1,1,e);
const double c_detJ = W[q] / ((J11*J22)-(J21*J12));
// (1/detJ) J^T C J
if (coeffDim == 3 || coeffDim == 4) // Matrix coefficient
{
const double C11 = C(0,q,e);
const double C12 = C(1,q,e);
const double C21 = symmetric ? C12 : C(2,q,e);
const double C22 = symmetric ? C(2,q,e) : C(3,q,e);
const double R11 = C11*J11 + C12*J21;
const double R21 = C21*J11 + C22*J21;
const double R12 = C11*J12 + C12*J22;
const double R22 = C21*J12 + C22*J22;
y(q,0,e) = c_detJ * (J11*R11 + J21*R21); // 1,1
y(q,1,e) = c_detJ * (J11*R12 + J21*R22); // 1,2
if (symmetric)
{
y(q,2,e) = c_detJ * (J12*R12 + J22*R22); // 2,2
}
else
{
y(q,2,e) = c_detJ * (J12*R11 + J22*R21); // 2,1
y(q,3,e) = c_detJ * (J12*R12 + J22*R22); // 2,2
}
}
else // Vector or scalar coefficient
{
const double C1 = C(0,q,e);
const double C2 = (coeffDim == 2 ? C(1,q,e) : C1);
y(q,0,e) = c_detJ * (J11*C1*J11 + J21*C2*J21); // 1,1
y(q,1,e) = c_detJ * (J11*C1*J12 + J21*C2*J22); // 1,2
y(q,2,e) = c_detJ * (J12*C1*J12 + J22*C2*J22); // 2,2
}
const double c_detJ = W[q] * coeff(q, e) / ((J11*J22)-(J21*J12));
// (c/detJ) J^T J
y(q,0,e) = c_detJ * (J11*J11 + J21*J21); // 1,1
y(q,1,e) = c_detJ * (J11*J12 + J21*J22); // 1,2
y(q,2,e) = c_detJ * (J12*J12 + J22*J22); // 2,2
}
});
}
// PA H(div) Mass Assemble 3D kernel
void PAHdivSetup3D(const int Q1D,
const int coeffDim,
const int NE,
const Array<double> &w,
const Vector &j,
Vector &coeff_,
Vector &op)
{
const bool symmetric = (coeffDim != 9);
const int NQ = Q1D*Q1D*Q1D;
auto W = w.Read();
auto J = Reshape(j.Read(), NQ, 3, 3, NE);
auto C = Reshape(coeff_.Read(), coeffDim, NQ, NE);
auto y = Reshape(op.Write(), NQ, symmetric ? 6 : 9, NE);
auto coeff = Reshape(coeff_.Read(), NQ, NE);
auto y = Reshape(op.Write(), NQ, 6, NE);
MFEM_FORALL(e, NE,
{
@@ -118,58 +84,14 @@ void PAHdivSetup3D(const int Q1D,
const double detJ = J11 * (J22 * J33 - J32 * J23) -
/* */ J21 * (J12 * J33 - J32 * J13) +
/* */ J31 * (J12 * J23 - J22 * J13);
const double c_detJ = W[q] / detJ;
// (1/detJ) J^T C J
if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
{
double M[3][3];
M[0][0] = C(0, q, e);
M[0][1] = C(1, q, e);
M[0][2] = C(2, q, e);
M[1][0] = (!symmetric) ? C(3, q, e) : M[0][1];
M[1][1] = (!symmetric) ? C(4, q, e) : C(3, q, e);
M[1][2] = (!symmetric) ? C(5, q, e) : C(4, q, e);
M[2][0] = (!symmetric) ? C(6, q, e) : M[0][2];
M[2][1] = (!symmetric) ? C(7, q, e) : M[1][2];
M[2][2] = (!symmetric) ? C(8, q, e) : C(5, q, e);
int idx = 0;
for (int i=0; i<3; ++i)
for (int j = (symmetric ? i : 0); j<3; ++j)
{
y(q,idx,e) = 0.0;
for (int k=0; k<3; ++k)
{
double MJ_kj = 0.0;
for (int l=0; l<3; ++l)
{
MJ_kj += M[k][l] * J(q,l,j,e);
}
y(q,idx,e) += J(q,k,i,e) * MJ_kj;
}
y(q,idx,e) *= c_detJ;
idx++;
}
}
else // Vector or scalar coefficient version
{
int idx = 0;
for (int i=0; i<3; ++i)
for (int j=i; j<3; ++j)
{
y(q,idx,e) = 0.0;
for (int k=0; k<3; ++k)
{
y(q,idx,e) += J(q,k,i,e) * C(coeffDim == 3 ? k : 0, q, e) * J(q,k,j,e);
}
y(q,idx,e) *= c_detJ;
idx++;
}
}
const double c_detJ = W[q] * coeff(q, e) / detJ;
// (c/detJ) J^T J
y(q,0,e) = c_detJ * (J11*J11 + J21*J21 + J31*J31); // 1,1
y(q,1,e) = c_detJ * (J12*J11 + J22*J21 + J32*J31); // 2,1
y(q,2,e) = c_detJ * (J13*J11 + J23*J21 + J33*J31); // 3,1
y(q,3,e) = c_detJ * (J12*J12 + J22*J22 + J32*J32); // 2,2
y(q,4,e) = c_detJ * (J13*J12 + J23*J22 + J33*J32); // 3,2
y(q,5,e) = c_detJ * (J13*J13 + J23*J23 + J33*J33); // 3,3
}
});
}
@@ -177,7 +99,6 @@ void PAHdivSetup3D(const int Q1D,
void PAHdivMassApply2D(const int D1D,
const int Q1D,
const int NE,
const bool symmetric,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Array<double> &Bot_,
@@ -194,7 +115,7 @@ void PAHdivMassApply2D(const int D1D,
auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
auto Bct = Reshape(Bct_.Read(), D1D, Q1D);
auto op = Reshape(op_.Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
auto op = Reshape(op_.Read(), Q1D, Q1D, 3, NE);
auto x = Reshape(x_.Read(), 2*(D1D-1)*D1D, NE);
auto y = Reshape(y_.ReadWrite(), 2*(D1D-1)*D1D, NE);
@@ -257,12 +178,11 @@ void PAHdivMassApply2D(const int D1D,
{
const double O11 = op(qx,qy,0,e);
const double O12 = op(qx,qy,1,e);
const double O21 = symmetric ? O12 : op(qx,qy,2,e);
const double O22 = symmetric ? op(qx,qy,2,e) : op(qx,qy,3,e);
const double O22 = op(qx,qy,2,e);
const double massX = mass[qy][qx][0];
const double massY = mass[qy][qx][1];
mass[qy][qx][0] = (O11*massX)+(O12*massY);
mass[qy][qx][1] = (O21*massX)+(O22*massY);
mass[qy][qx][1] = (O12*massX)+(O22*massY);
}
}
@@ -305,179 +225,9 @@ void PAHdivMassApply2D(const int D1D,
}); // end of element loop
}
template<int T_D1D = 0, int T_Q1D = 0>
void SmemPAHdivMassApply2D(const int NE,
const bool symmetric,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Array<double> &Bot_,
const Array<double> &Bct_,
const Vector &op_,
const Vector &x_,
Vector &y_,
const int d1d = 0,
const int q1d = 0)
{
MFEM_CONTRACT_VAR(Bot_);
MFEM_CONTRACT_VAR(Bct_);
static constexpr int VDIM = 2;
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
const auto bo = Reshape(Bo_.Read(), Q1D, D1D-1);
const auto bc = Reshape(Bc_.Read(), Q1D, D1D);
const auto D = Reshape(op_.Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
const auto x = Reshape(x_.Read(), D1D*(D1D-1), VDIM, NE);
auto y = y_.ReadWrite();
MFEM_FORALL_3D(e, NE, Q1D, Q1D, VDIM,
{
const int tidz = MFEM_THREAD_ID(z);
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
constexpr int MQ1 = T_Q1D ? T_Q1D : HDIV_MAX_Q1D;
constexpr int MD1 = T_D1D ? T_D1D : HDIV_MAX_D1D;
constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
MFEM_SHARED double smo[MQ1*(MD1-1)];
DeviceMatrix Bo(smo, D1D-1, Q1D);
MFEM_SHARED double smc[MQ1*MD1];
DeviceMatrix Bc(smc, D1D, Q1D);
MFEM_SHARED double sm0[VDIM*MDQ*MDQ];
MFEM_SHARED double sm1[VDIM*MDQ*MDQ];
DeviceMatrix X(sm0, D1D*(D1D-1), VDIM);
DeviceCube QD(sm1, Q1D, D1D, VDIM);
DeviceCube QQ(sm0, Q1D, Q1D, VDIM);
// Load X, Bo and Bc into shared memory
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
MFEM_FOREACH_THREAD(dy,y,D1D)
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
if (qx < D1D && dy < (D1D-1)) { X(qx + dy*D1D,vd) = x(qx+dy*D1D,vd,e); }
if (tidz == 0)
{
if (dy < (D1D-1)) { Bo(dy,qx) = bo(qx,dy); }
Bc(dy,qx) = bc(qx,dy);
}
}
}
}
MFEM_SYNC_THREAD;
// Apply B operator
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int nx = (vd == 0) ? D1D : D1D-1;
const int ny = (vd == 1) ? D1D : D1D-1;
DeviceCube Xxy(X, nx, ny, VDIM);
DeviceMatrix Bx = (vd == 0) ? Bc : Bo;
MFEM_FOREACH_THREAD(dy,y,ny)
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
double dq = 0.0;
for (int dx = 0; dx < nx; ++dx)
{
dq += Xxy(dx,dy,vd) * Bx(dx,qx);
}
QD(qx,dy,vd) = dq;
}
}
}
MFEM_SYNC_THREAD;
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int ny = (vd == 1) ? D1D : D1D-1;
DeviceMatrix By = (vd == 1) ? Bc : Bo;
MFEM_FOREACH_THREAD(qy,y,Q1D)
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
double qq = 0.0;
for (int dy = 0; dy < ny; ++dy)
{
qq += QD(qx,dy,vd) * By(dy,qy);
}
QQ(qx,qy,vd) = qq;
}
}
}
MFEM_SYNC_THREAD;
// Apply D operator
if (tidz == 0)
{
MFEM_FOREACH_THREAD(qy,y,Q1D)
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
const double Qx = QQ(qx,qy,0);
const double Qy = QQ(qx,qy,1);
const double D11 = D(qx,qy,0,e);
const double D12 = D(qx,qy,1,e);
const double D21 = symmetric ? D12 : D(qx,qy,2,e);
const double D22 = symmetric ? D(qx,qy,2,e) : D(qx,qy,3,e);
QQ(qx,qy,0) = D11*Qx + D12*Qy;
QQ(qx,qy,1) = D21*Qx + D22*Qy;
}
}
}
MFEM_SYNC_THREAD;
// Apply Bt operator
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int nx = (vd == 0) ? D1D : D1D-1;
DeviceMatrix Btx = (vd == 0) ? Bc : Bo;
MFEM_FOREACH_THREAD(qy,y,Q1D)
{
MFEM_FOREACH_THREAD(dx,x,nx)
{
double qd = 0.0;
for (int qx = 0; qx < Q1D; ++qx)
{
qd += QQ(qx,qy,vd) * Btx(dx,qx);
}
QD(dx,qy,vd) = qd;
}
}
}
MFEM_SYNC_THREAD;
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int nx = (vd == 0) ? D1D : D1D-1;
const int ny = (vd == 1) ? D1D : D1D-1;
DeviceMatrix Bty = (vd == 1) ? Bc : Bo;
DeviceTensor<4> Yxy(y, nx, ny, VDIM, NE);
MFEM_FOREACH_THREAD(dy,y,ny)
{
MFEM_FOREACH_THREAD(dx,x,nx)
{
double dd = 0.0;
for (int qy = 0; qy < Q1D; ++qy)
{
dd += QD(dx,qy,vd) * Bty(dy,qy);
}
Yxy(dx,dy,vd,e) += dd;
}
}
}
MFEM_SYNC_THREAD;
});
}
void PAHdivMassAssembleDiagonal2D(const int D1D,
const int Q1D,
const int NE,
const bool symmetric,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Vector &op_,
@@ -488,7 +238,7 @@ void PAHdivMassAssembleDiagonal2D(const int D1D,
auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
auto op = Reshape(op_.Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
auto op = Reshape(op_.Read(), Q1D, Q1D, 3, NE);
auto diag = Reshape(diag_.ReadWrite(), 2*(D1D-1)*D1D, NE);
MFEM_FORALL(e, NE,
@@ -509,7 +259,7 @@ void PAHdivMassAssembleDiagonal2D(const int D1D,
for (int qy = 0; qy < Q1D; ++qy)
{
const double wy = (c == 1) ? Bc(qy,dy) : Bo(qy,dy);
mass[qx] += wy*wy*((c == 0) ? op(qx,qy,0,e) : op(qx,qy,symmetric ? 2 : 3,e));
mass[qx] += wy*wy*((c == 0) ? op(qx,qy,0,e) : op(qx,qy,2,e));
}
}
@@ -533,7 +283,6 @@ void PAHdivMassAssembleDiagonal2D(const int D1D,
void PAHdivMassAssembleDiagonal3D(const int D1D,
const int Q1D,
const int NE,
const bool symmetric,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Vector &op_,
@@ -545,7 +294,7 @@ void PAHdivMassAssembleDiagonal3D(const int D1D,
auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, 6, NE);
auto diag = Reshape(diag_.ReadWrite(), 3*(D1D-1)*(D1D-1)*D1D, NE);
MFEM_FORALL(e, NE,
@@ -558,8 +307,7 @@ void PAHdivMassAssembleDiagonal3D(const int D1D,
const int D1Dy = (c == 1) ? D1D : D1D - 1;
const int D1Dx = (c == 0) ? D1D : D1D - 1;
const int opc = (c == 0) ? 0 : ((c == 1) ? (symmetric ? 3 : 4) :
(symmetric ? 5 : 8));
const int opc = (c == 0) ? 0 : ((c == 1) ? 3 : 5);
double mass[HDIV_MAX_Q1D];
@@ -602,7 +350,6 @@ void PAHdivMassAssembleDiagonal3D(const int D1D,
void PAHdivMassApply3D(const int D1D,
const int Q1D,
const int NE,
const bool symmetric,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Array<double> &Bot_,
@@ -619,7 +366,7 @@ void PAHdivMassApply3D(const int D1D,
auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
auto Bct = Reshape(Bct_.Read(), D1D, Q1D);
auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
auto op = Reshape(op_.Read(), Q1D, Q1D, Q1D, 6, NE);
auto x = Reshape(x_.Read(), 3*(D1D-1)*(D1D-1)*D1D, NE);
auto y = Reshape(y_.ReadWrite(), 3*(D1D-1)*(D1D-1)*D1D, NE);
@@ -714,19 +461,15 @@ void PAHdivMassApply3D(const int D1D,
const double O11 = op(qx,qy,qz,0,e);
const double O12 = op(qx,qy,qz,1,e);
const double O13 = op(qx,qy,qz,2,e);
const double O21 = symmetric ? O12 : op(qx,qy,qz,3,e);
const double O22 = symmetric ? op(qx,qy,qz,3,e) : op(qx,qy,qz,4,e);
const double O23 = symmetric ? op(qx,qy,qz,4,e) : op(qx,qy,qz,5,e);
const double O31 = symmetric ? O13 : op(qx,qy,qz,6,e);
const double O32 = symmetric ? O23 : op(qx,qy,qz,7,e);
const double O33 = symmetric ? op(qx,qy,qz,5,e) : op(qx,qy,qz,8,e);
const double O22 = op(qx,qy,qz,3,e);
const double O23 = op(qx,qy,qz,4,e);
const double O33 = op(qx,qy,qz,5,e);
const double massX = mass[qz][qy][qx][0];
const double massY = mass[qz][qy][qx][1];
const double massZ = mass[qz][qy][qx][2];
mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
}
}
}
@@ -794,337 +537,6 @@ void PAHdivMassApply3D(const int D1D,
}); // end of element loop
}
template<int T_D1D = 0, int T_Q1D = 0>
void SmemPAHdivMassApply3D(const int NE,
const bool symmetric,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Array<double> &Bot_,
const Array<double> &Bct_,
const Vector &op_,
const Vector &x_,
Vector &y_,
const int d1d = 0,
const int q1d = 0)
{
MFEM_CONTRACT_VAR(Bot_);
MFEM_CONTRACT_VAR(Bct_);
static constexpr int VDIM = 3;
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
const auto bo = Reshape(Bo_.Read(), Q1D, D1D-1);
const auto bc = Reshape(Bc_.Read(), Q1D, D1D);
const auto D = Reshape(op_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
const auto x = Reshape(x_.Read(), D1D*(D1D-1)*(D1D-1), VDIM, NE);
auto y = y_.ReadWrite();
MFEM_FORALL_3D(e, NE, Q1D, Q1D, VDIM,
{
const int tidz = MFEM_THREAD_ID(z);
const int D1D = T_D1D ? T_D1D : d1d;
const int Q1D = T_Q1D ? T_Q1D : q1d;
constexpr int MQ1 = T_Q1D ? T_Q1D : HDIV_MAX_Q1D;
constexpr int MD1 = T_D1D ? T_D1D : HDIV_MAX_D1D;
constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
MFEM_SHARED double smo[MQ1*(MD1-1)];
DeviceMatrix Bo(smo, D1D-1, Q1D);
MFEM_SHARED double smc[MQ1*MD1];
DeviceMatrix Bc(smc, D1D, Q1D);
MFEM_SHARED double sm0[VDIM*MDQ*MDQ*MDQ];
MFEM_SHARED double sm1[VDIM*MDQ*MDQ*MDQ];
DeviceMatrix X(sm0, D1D*(D1D-1)*(D1D-1), VDIM);
DeviceTensor<4> QDD(sm1, Q1D, D1D, D1D, VDIM);
DeviceTensor<4> QQD(sm0, Q1D, Q1D, D1D, VDIM);
DeviceTensor<4> QQQ(sm1, Q1D, Q1D, Q1D, VDIM);
DeviceTensor<4> DQQ(sm0, D1D, Q1D, Q1D, VDIM);
DeviceTensor<4> DDQ(sm1, D1D, D1D, Q1D, VDIM);
// Load X into shared memory
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
MFEM_FOREACH_THREAD(dz,y,D1D-1)
{
MFEM_FOREACH_THREAD(dy,x,D1D-1)
{
MFEM_UNROLL(MD1)
for (int dx = 0; dx < D1D; ++dx)
{
X(dx+(dy+dz*(D1D-1))*D1D,vd) = x(dx+(dy+dz*(D1D-1))*D1D,vd,e);
}
}
}
}
// Load Bo and Bc into shared memory
if (tidz == 0)
{
MFEM_FOREACH_THREAD(d,y,D1D-1)
{
MFEM_FOREACH_THREAD(q,x,Q1D)
{
Bo(d,q) = bo(q,d);
}
}
MFEM_FOREACH_THREAD(d,y,D1D)
{
MFEM_FOREACH_THREAD(q,x,Q1D)
{
Bc(d,q) = bc(q,d);
}
}
}
MFEM_SYNC_THREAD;
// Apply B operator
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int nx = (vd == 0) ? D1D : D1D-1;
const int ny = (vd == 1) ? D1D : D1D-1;
const int nz = (vd == 2) ? D1D : D1D-1;
DeviceTensor<4> Xxyz(X, nx, ny, nz, VDIM);
DeviceMatrix Bx = (vd == 0) ? Bc : Bo;
MFEM_FOREACH_THREAD(dy,y,ny)
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
double u[D1D];
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
MFEM_UNROLL(MD1)
for (int dx = 0; dx < nx; ++dx)
{
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz)
{
u[dz] += Xxyz(dx,dy,dz,vd) * Bx(dx,qx);
}
}
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz) { QDD(qx,dy,dz,vd) = u[dz]; }
}
}
}
MFEM_SYNC_THREAD;
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int ny = (vd == 1) ? D1D : D1D-1;
const int nz = (vd == 2) ? D1D : D1D-1;
DeviceMatrix By = (vd == 1) ? Bc : Bo;
MFEM_FOREACH_THREAD(qy,y,Q1D)
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
double u[D1D];
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
MFEM_UNROLL(MD1)
for (int dy = 0; dy < ny; ++dy)
{
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz)
{
u[dz] += QDD(qx,dy,dz,vd) * By(dy,qy);
}
}
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz) { QQD(qx,qy,dz,vd) = u[dz]; }
}
}
}
MFEM_SYNC_THREAD;
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int nz = (vd == 2) ? D1D : D1D-1;
DeviceMatrix Bz = (vd == 2) ? Bc : Bo;
MFEM_FOREACH_THREAD(qy,y,Q1D)
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
double u[Q1D];
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz) { u[qz] = 0.0; }
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz)
{
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz)
{
u[qz] += QQD(qx,qy,dz,vd) * Bz(dz,qz);
}
}
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz) { QQQ(qx,qy,qz,vd) = u[qz]; }
}
}
}
MFEM_SYNC_THREAD;
// Apply D operator
if (tidz == 0)
{
MFEM_FOREACH_THREAD(qy,y,Q1D)
{
MFEM_FOREACH_THREAD(qx,x,Q1D)
{
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz)
{
const double Qx = QQQ(qx,qy,qz,0);
const double Qy = QQQ(qx,qy,qz,1);
const double Qz = QQQ(qx,qy,qz,2);
const double D11 = D(qx,qy,qz,0,e);
const double D12 = D(qx,qy,qz,1,e);
const double D13 = D(qx,qy,qz,2,e);
const double D21 = symmetric ? D12 : D(qx,qy,qz,3,e);
const double D22 = symmetric ? D(qx,qy,qz,3,e) : D(qx,qy,qz,4,e);
const double D23 = symmetric ? D(qx,qy,qz,4,e) : D(qx,qy,qz,5,e);
const double D31 = symmetric ? D13 : D(qx,qy,qz,6,e);
const double D32 = symmetric ? D23 : D(qx,qy,qz,7,e);
const double D33 = symmetric ? D(qx,qy,qz,5,e) : D(qx,qy,qz,8,e);
QQQ(qx,qy,qz,0) = D11*Qx + D12*Qy + D13*Qz;
QQQ(qx,qy,qz,1) = D21*Qx + D22*Qy + D23*Qz;
QQQ(qx,qy,qz,2) = D31*Qx + D32*Qy + D33*Qz;
}
}
}
}
MFEM_SYNC_THREAD;
// Apply Bt operator
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int nx = (vd == 0) ? D1D : D1D-1;
DeviceMatrix Btx = (vd == 0) ? Bc : Bo;
MFEM_FOREACH_THREAD(qy,y,Q1D)
{
MFEM_FOREACH_THREAD(dx,x,nx)
{
double u[Q1D];
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz) { u[qz] = 0.0; }
MFEM_UNROLL(MQ1)
for (int qx = 0; qx < Q1D; ++qx)
{
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz)
{
u[qz] += QQQ(qx,qy,qz,vd) * Btx(dx,qx);
}
}
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz) { DQQ(dx,qy,qz,vd) = u[qz]; }
}
}
}
MFEM_SYNC_THREAD;
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int nx = (vd == 0) ? D1D : D1D-1;
const int ny = (vd == 1) ? D1D : D1D-1;
DeviceMatrix Bty = (vd == 1) ? Bc : Bo;
MFEM_FOREACH_THREAD(dy,y,ny)
{
MFEM_FOREACH_THREAD(dx,x,nx)
{
double u[Q1D];
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz) { u[qz] = 0.0; }
MFEM_UNROLL(MQ1)
for (int qy = 0; qy < Q1D; ++qy)
{
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz)
{
u[qz] += DQQ(dx,qy,qz,vd) * Bty(dy,qy);
}
}
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz) { DDQ(dx,dy,qz,vd) = u[qz]; }
}
}
}
MFEM_SYNC_THREAD;
MFEM_FOREACH_THREAD(vd,z,VDIM)
{
const int nx = (vd == 0) ? D1D : D1D-1;
const int ny = (vd == 1) ? D1D : D1D-1;
const int nz = (vd == 2) ? D1D : D1D-1;
DeviceTensor<5> Yxyz(y, nx, ny, nz, VDIM, NE);
DeviceMatrix Btz = (vd == 2) ? Bc : Bo;
MFEM_FOREACH_THREAD(dy,y,ny)
{
MFEM_FOREACH_THREAD(dx,x,nx)
{
double u[D1D];
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz) { u[dz] = 0.0; }
MFEM_UNROLL(MQ1)
for (int qz = 0; qz < Q1D; ++qz)
{
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz)
{
u[dz] += DDQ(dx,dy,qz,vd) * Btz(dz,qz);
}
}
MFEM_UNROLL(MD1)
for (int dz = 0; dz < nz; ++dz) { Yxyz(dx,dy,dz,vd,e) += u[dz]; }
}
}
}
MFEM_SYNC_THREAD;
});
}
void PAHdivMassApply(const int dim,
const int D1D,
const int Q1D,
const int NE,
const bool symmetric,
const Array<double> &Bo,
const Array<double> &Bc,
const Array<double> &Bot,
const Array<double> &Bct,
const Vector &op,
const Vector &x,
Vector &y)
{
const int id = (D1D << 4) | Q1D;
if (dim == 2)
{
switch (id)
{
case 0x22: return SmemPAHdivMassApply2D<2,2>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
case 0x33: return SmemPAHdivMassApply2D<3,3>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
case 0x44: return SmemPAHdivMassApply2D<4,4>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
case 0x55: return SmemPAHdivMassApply2D<5,5>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
default: // fallback
return PAHdivMassApply2D(D1D,Q1D,NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
}
}
else if (dim == 3)
{
switch (id)
{
case 0x23: return SmemPAHdivMassApply3D<2,3>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
case 0x34: return SmemPAHdivMassApply3D<3,4>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
case 0x45: return SmemPAHdivMassApply3D<4,5>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
case 0x56: return SmemPAHdivMassApply3D<5,6>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
case 0x67: return SmemPAHdivMassApply3D<6,7>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
case 0x78: return SmemPAHdivMassApply3D<7,8>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
default: // fallback
return PAHdivMassApply3D(D1D,Q1D,NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
}
}
}
// PA H(div) div-div assemble 2D kernel
// NOTE: this is identical to PACurlCurlSetup3D
static void PADivDivSetup2D(const int Q1D,
@@ -1214,7 +626,7 @@ static void PADivDivApply2D(const int D1D,
{
double div[MAX_Q1D][MAX_Q1D];
// div[qy][qx] will be computed as du_x/dx + du_y/dy
// div[qy][qx] will be computed as du_x/dx + duy_/dy
for (int qy = 0; qy < Q1D; ++qy)
{
@@ -1797,13 +1209,6 @@ VectorFEDivergenceIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
}
}
if (test_el->GetMapType() == FiniteElement::INTEGRAL)
{
const GeometricFactors *geom =
mesh->GetGeometricFactors(*ir, GeometricFactors::DETERMINANTS);
coeff /= geom->detJ;
}
if (trial_el->GetDerivType() == mfem::FiniteElement::DIV && dim == 3)
{
PADivL2Setup3D(quad1D, ne, ir->GetWeights(), coeff, pa_data);
+1 -10
View File
@@ -31,16 +31,7 @@ void MassIntegrator::AssembleMF(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedMFMassIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::MFMassIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::MFMassIntegrator(fes, *ir, Q);
return;
}
MFEM_ABORT("Error: MassIntegrator::AssembleMF only implemented with"
+1 -10
View File
@@ -38,16 +38,7 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedPAMassIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
return;
}
int map_type = el.GetMapType();
+1 -10
View File
@@ -149,16 +149,7 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedPADiffusionIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::PADiffusionIntegrator(fes, *ir, Q);
return;
}
const int dims = el.GetDim();
+1 -13
View File
@@ -30,19 +30,7 @@ void VectorDiffusionIntegrator::AssembleMF(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
MFEM_VERIFY(!VQ && !MQ,
"Only scalar coefficient supported for DiffusionIntegrator"
" with libCEED");
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedMFDiffusionIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::MFDiffusionIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::MFDiffusionIntegrator(fes, *ir, Q);
return;
}
MFEM_ABORT("Error: VectorDiffusionIntegrator::AssembleMF only implemented"
+1 -10
View File
@@ -34,16 +34,7 @@ void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedPAMassIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::PAMassIntegrator(fes, *ir, Q);
return;
}
dim = mesh->Dimension();
+1 -10
View File
@@ -34,16 +34,7 @@ void VectorMassIntegrator::AssembleMF(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedMFMassIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::MFMassIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::MFMassIntegrator(fes, *ir, Q);
return;
}
MFEM_ABORT("Error: VectorMassIntegrator::AssembleMF only implemented with"
+83 -95
View File
@@ -11,7 +11,6 @@
#include "../general/forall.hpp"
#include "bilininteg.hpp"
#include "gridfunc.hpp"
namespace mfem
{
@@ -90,7 +89,6 @@ void SmemPAHcurlMassApply3D(const int D1D,
Vector &y);
void PAHdivSetup2D(const int Q1D,
const int coeffDim,
const int NE,
const Array<double> &w,
const Vector &j,
@@ -98,7 +96,6 @@ void PAHdivSetup2D(const int Q1D,
Vector &op);
void PAHdivSetup3D(const int Q1D,
const int coeffDim,
const int NE,
const Array<double> &w,
const Vector &j,
@@ -152,7 +149,6 @@ void PAHcurlH1ApplyTranspose3D(const int D1D,
void PAHdivMassAssembleDiagonal2D(const int D1D,
const int Q1D,
const int NE,
const bool symmetric,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Vector &op_,
@@ -161,24 +157,32 @@ void PAHdivMassAssembleDiagonal2D(const int D1D,
void PAHdivMassAssembleDiagonal3D(const int D1D,
const int Q1D,
const int NE,
const bool symmetric,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Vector &op_,
Vector &diag_);
void PAHdivMassApply(const int dim,
const int D1D,
const int Q1D,
const int NE,
const bool symmetric,
const Array<double> &Bo,
const Array<double> &Bc,
const Array<double> &Bot,
const Array<double> &Bct,
const Vector &op,
const Vector &x,
Vector &y);
void PAHdivMassApply2D(const int D1D,
const int Q1D,
const int NE,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Array<double> &Bot_,
const Array<double> &Bct_,
const Vector &op_,
const Vector &x_,
Vector &y_);
void PAHdivMassApply3D(const int D1D,
const int Q1D,
const int NE,
const Array<double> &Bo_,
const Array<double> &Bc_,
const Array<double> &Bot_,
const Array<double> &Bct_,
const Vector &op_,
const Vector &x_,
Vector &y_);
void PAHcurlL2Setup(const int NQ,
const int coeffDim,
@@ -814,80 +818,69 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
pa_data.SetSize((symmetric ? symmDims : MQfullDim) * nq * ne,
Device::GetMemoryType());
Vector coeff;
Vector coeff(coeffDim * ne * nq);
coeff = 1.0;
auto coeffh = Reshape(coeff.HostWrite(), coeffDim, nq, ne);
if (Q || DQ || MQ)
{
Vector DM(DQ ? coeffDim : 0);
DenseMatrix M;
DenseSymmetricMatrix SM;
auto *qf_c = dynamic_cast<QuadratureFunctionCoefficient*>(Q);
if (qf_c)
{
const QuadratureFunction &qf = qf_c->GetQuadFunction();
qf.Read();
coeff.MakeRef(const_cast<QuadratureFunction&>(qf), 0);
}
else
{
coeff.SetSize(coeffDim * ne * nq);
coeff = 1.0;
auto coeffh = Reshape(coeff.HostWrite(), coeffDim, nq, ne);
if (Q || DQ || MQ)
if (DQ)
{
Vector DM(DQ ? coeffDim : 0);
DenseMatrix M;
DenseSymmetricMatrix SM;
MFEM_VERIFY(coeffDim == dim, "");
}
if (SMQ)
{
MFEM_VERIFY(SMQ->GetSize() == dim, "");
SM.SetSize(dim);
}
else if (MQ)
{
MFEM_VERIFY(coeffDim == MQdim, "");
MFEM_VERIFY(MQ->GetHeight() == dim && MQ->GetWidth() == dim, "");
M.SetSize(dim);
}
if (DQ)
{
MFEM_VERIFY(coeffDim == dim, "");
}
if (SMQ)
{
MFEM_VERIFY(SMQ->GetSize() == dim, "");
SM.SetSize(dim);
}
else if (MQ)
{
MFEM_VERIFY(coeffDim == MQdim, "");
MFEM_VERIFY(MQ->GetHeight() == dim && MQ->GetWidth() == dim, "");
M.SetSize(dim);
}
for (int e=0; e<ne; ++e)
for (int e=0; e<ne; ++e)
{
ElementTransformation *tr = mesh->GetElementTransformation(e);
for (int p=0; p<nq; ++p)
{
ElementTransformation *tr = mesh->GetElementTransformation(e);
for (int p=0; p<nq; ++p)
if (SMQ)
{
if (SMQ)
{
SMQ->Eval(SM, *tr, ir->IntPoint(p));
int cnt = 0;
for (int i=0; i<dim; ++i)
for (int j=i; j<dim; ++j, ++cnt)
{
coeffh(cnt, p, e) = SM(i,j);
}
}
else if (MQ)
{
MQ->Eval(M, *tr, ir->IntPoint(p));
for (int i=0; i<dim; ++i)
for (int j=0; j<dim; ++j)
{
coeffh(j+(i*dim), p, e) = M(i,j);
}
}
else if (DQ)
{
DQ->Eval(DM, *tr, ir->IntPoint(p));
for (int i=0; i<coeffDim; ++i)
SMQ->Eval(SM, *tr, ir->IntPoint(p));
int cnt = 0;
for (int i=0; i<dim; ++i)
for (int j=i; j<dim; ++j, ++cnt)
{
coeffh(i, p, e) = DM[i];
coeffh(cnt, p, e) = SM(i,j);
}
}
else
}
else if (MQ)
{
MQ->Eval(M, *tr, ir->IntPoint(p));
for (int i=0; i<dim; ++i)
for (int j=0; j<dim; ++j)
{
coeffh(j+(i*dim), p, e) = M(i,j);
}
}
else if (DQ)
{
DQ->Eval(DM, *tr, ir->IntPoint(p));
for (int i=0; i<coeffDim; ++i)
{
coeffh(0, p, e) = Q->Eval(*tr, ir->IntPoint(p));
coeffh(i, p, e) = DM[i];
}
}
else
{
coeffh(0, p, e) = Q->Eval(*tr, ir->IntPoint(p));
}
}
}
}
@@ -904,12 +897,12 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
}
else if (trial_div && test_div && dim == 3)
{
PAHdivSetup3D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J,
PAHdivSetup3D(quad1D, ne, ir->GetWeights(), geom->J,
coeff, pa_data);
}
else if (trial_div && test_div && dim == 2)
{
PAHdivSetup2D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J,
PAHdivSetup2D(quad1D, ne, ir->GetWeights(), geom->J,
coeff, pa_data);
}
else if (((trial_curl && test_div) || (trial_div && test_curl)) &&
@@ -970,7 +963,7 @@ void VectorFEMassIntegrator::AssembleDiagonalPA(Vector& diag)
else if (trial_fetype == mfem::FiniteElement::DIV &&
test_fetype == trial_fetype)
{
PAHdivMassAssembleDiagonal3D(dofs1D, quad1D, ne, symmetric,
PAHdivMassAssembleDiagonal3D(dofs1D, quad1D, ne,
mapsO->B, mapsC->B, pa_data, diag);
}
else
@@ -978,7 +971,7 @@ void VectorFEMassIntegrator::AssembleDiagonalPA(Vector& diag)
MFEM_ABORT("Unknown kernel.");
}
}
else // 2D
else
{
if (trial_fetype == mfem::FiniteElement::CURL && test_fetype == trial_fetype)
{
@@ -988,7 +981,7 @@ void VectorFEMassIntegrator::AssembleDiagonalPA(Vector& diag)
else if (trial_fetype == mfem::FiniteElement::DIV &&
test_fetype == trial_fetype)
{
PAHdivMassAssembleDiagonal2D(dofs1D, quad1D, ne, symmetric,
PAHdivMassAssembleDiagonal2D(dofs1D, quad1D, ne,
mapsO->B, mapsC->B, pa_data, diag);
}
else
@@ -1041,8 +1034,8 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
}
else if (trial_div && test_div)
{
PAHdivMassApply(3, dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B, mapsO->Bt,
mapsC->Bt, pa_data, x, y);
PAHdivMassApply3D(dofs1D, quad1D, ne, mapsO->B, mapsC->B, mapsO->Bt,
mapsC->Bt, pa_data, x, y);
}
else if (trial_curl && test_div)
{
@@ -1063,7 +1056,7 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
MFEM_ABORT("Unknown kernel.");
}
}
else // 2D
else
{
if (trial_curl && test_curl)
{
@@ -1072,8 +1065,8 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
}
else if (trial_div && test_div)
{
PAHdivMassApply(2, dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B, mapsO->Bt,
mapsC->Bt, pa_data, x, y);
PAHdivMassApply2D(dofs1D, quad1D, ne, mapsO->B, mapsC->B, mapsO->Bt,
mapsC->Bt, pa_data, x, y);
}
else if ((trial_curl && test_div) || (trial_div && test_curl))
{
@@ -1118,11 +1111,6 @@ void VectorFEMassIntegrator::AddMultTransposePA(const Vector &x,
if (symmetricSpaces)
{
if (MQ && dynamic_cast<SymmetricMatrixCoefficient*>(MQ) == NULL)
{
MFEM_ABORT("VectorFEMassIntegrator transpose not implemented for asymmetric MatrixCoefficient");
}
this->AddMultPA(x, y);
}
}
@@ -62,20 +62,6 @@ PAConvectionIntegrator::PAConvectionIntegrator(
#endif
}
MixedPAConvectionIntegrator::MixedPAConvectionIntegrator(
const ConvectionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::VectorCoefficient *Q,
const double alpha)
{
#ifdef MFEM_USE_CEED
ConvectionOperatorInfo info(fes.GetMesh()->Dimension(), alpha);
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
MFConvectionIntegrator::MFConvectionIntegrator(
const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &irm,
@@ -91,20 +77,6 @@ MFConvectionIntegrator::MFConvectionIntegrator(
#endif
}
MixedMFConvectionIntegrator::MixedMFConvectionIntegrator(
const ConvectionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::VectorCoefficient *Q,
const double alpha)
{
#ifdef MFEM_USE_CEED
ConvectionOperatorInfo info(fes.GetMesh()->Dimension(), alpha);
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
} // namespace ceed
} // namespace mfem
+2 -21
View File
@@ -13,7 +13,6 @@
#define MFEM_LIBCEED_CONV_HPP
#include "../../interface/integrator.hpp"
#include "../../interface/mixed_integrator.hpp"
#include "../../../fespace.hpp"
namespace mfem
@@ -27,39 +26,21 @@ class PAConvectionIntegrator : public PAIntegrator
{
public:
PAConvectionIntegrator(const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
const mfem::IntegrationRule &irm,
mfem::VectorCoefficient *Q,
const double alpha);
};
class MixedPAConvectionIntegrator : public MixedIntegrator<PAIntegrator>
{
public:
MixedPAConvectionIntegrator(const ConvectionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::VectorCoefficient *Q,
const double alpha);
};
/// Represent a ConvectionIntegrator with AssemblyLevel::None using libCEED.
class MFConvectionIntegrator : public MFIntegrator
{
public:
MFConvectionIntegrator(const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
const mfem::IntegrationRule &irm,
mfem::VectorCoefficient *Q,
const double alpha);
};
class MixedMFConvectionIntegrator : public MixedIntegrator<MFIntegrator>
{
public:
MixedMFConvectionIntegrator(const ConvectionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::VectorCoefficient *Q,
const double alpha);
};
}
}
@@ -60,32 +60,6 @@ PADiffusionIntegrator::PADiffusionIntegrator(
#endif
}
MixedPADiffusionIntegrator::MixedPADiffusionIntegrator(
const DiffusionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
DiffusionOperatorInfo info(fes.GetMesh()->Dimension());
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
MixedPADiffusionIntegrator::MixedPADiffusionIntegrator(
const VectorDiffusionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
DiffusionOperatorInfo info(fes.GetMesh()->Dimension());
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
MFDiffusionIntegrator::MFDiffusionIntegrator(
const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &irm,
@@ -100,32 +74,6 @@ MFDiffusionIntegrator::MFDiffusionIntegrator(
#endif
}
MixedMFDiffusionIntegrator::MixedMFDiffusionIntegrator(
const DiffusionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
DiffusionOperatorInfo info(fes.GetMesh()->Dimension());
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
MixedMFDiffusionIntegrator::MixedMFDiffusionIntegrator(
const VectorDiffusionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
DiffusionOperatorInfo info(fes.GetMesh()->Dimension());
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
} // namespace ceed
} // namespace mfem
+2 -27
View File
@@ -13,7 +13,6 @@
#define MFEM_LIBCEED_DIFF_HPP
#include "../../interface/integrator.hpp"
#include "../../interface/mixed_integrator.hpp"
#include "../../../fespace.hpp"
namespace mfem
@@ -27,43 +26,19 @@ class PADiffusionIntegrator : public PAIntegrator
{
public:
PADiffusionIntegrator(const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
const mfem::IntegrationRule &irm,
mfem::Coefficient *Q);
};
class MixedPADiffusionIntegrator : public MixedIntegrator<PAIntegrator>
{
public:
MixedPADiffusionIntegrator(const DiffusionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
MixedPADiffusionIntegrator(const VectorDiffusionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
};
/// Represent a DiffusionIntegrator with AssemblyLevel::None using libCEED.
class MFDiffusionIntegrator : public MFIntegrator
{
public:
MFDiffusionIntegrator(const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
const mfem::IntegrationRule &irm,
mfem::Coefficient *Q);
};
class MixedMFDiffusionIntegrator : public MixedIntegrator<MFIntegrator>
{
public:
MixedMFDiffusionIntegrator(const DiffusionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
MixedMFDiffusionIntegrator(const VectorDiffusionIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
};
}
}
-48
View File
@@ -59,30 +59,6 @@ PAMassIntegrator::PAMassIntegrator(const mfem::FiniteElementSpace &fes,
#endif
}
MixedPAMassIntegrator::MixedPAMassIntegrator(const MassIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
MassOperatorInfo info;
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
MixedPAMassIntegrator::MixedPAMassIntegrator(const VectorMassIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
MassOperatorInfo info;
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
MFMassIntegrator::MFMassIntegrator(const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &irm,
mfem::Coefficient *Q)
@@ -96,30 +72,6 @@ MFMassIntegrator::MFMassIntegrator(const mfem::FiniteElementSpace &fes,
#endif
}
MixedMFMassIntegrator::MixedMFMassIntegrator(const MassIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
MassOperatorInfo info;
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
MixedMFMassIntegrator::MixedMFMassIntegrator(const VectorMassIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
MassOperatorInfo info;
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
} // namespace ceed
} // namespace mfem
+2 -27
View File
@@ -13,7 +13,6 @@
#define MFEM_LIBCEED_MASS_HPP
#include "../../interface/integrator.hpp"
#include "../../interface/mixed_integrator.hpp"
#include "../../../fespace.hpp"
namespace mfem
@@ -27,43 +26,19 @@ class PAMassIntegrator : public PAIntegrator
{
public:
PAMassIntegrator(const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
const mfem::IntegrationRule &irm,
mfem::Coefficient *Q);
};
class MixedPAMassIntegrator : public MixedIntegrator<PAIntegrator>
{
public:
MixedPAMassIntegrator(const MassIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
MixedPAMassIntegrator(const VectorMassIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
};
/// Represent a MassIntegrator with AssemblyLevel::None using libCEED.
class MFMassIntegrator : public MFIntegrator
{
public:
MFMassIntegrator(const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
const mfem::IntegrationRule &irm,
mfem::Coefficient *Q);
};
class MixedMFMassIntegrator : public MixedIntegrator<MFIntegrator>
{
public:
MixedMFMassIntegrator(const MassIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
MixedMFMassIntegrator(const VectorMassIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
};
}
}
@@ -60,19 +60,6 @@ PAVectorConvectionNLFIntegrator::PAVectorConvectionNLFIntegrator(
#endif
}
MixedPAVectorConvectionNLIntegrator::MixedPAVectorConvectionNLIntegrator(
const VectorConvectionNLFIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
NLConvectionOperatorInfo info(fes.GetMesh()->Dimension());
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
MFVectorConvectionNLFIntegrator::MFVectorConvectionNLFIntegrator(
const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &irm,
@@ -87,19 +74,6 @@ MFVectorConvectionNLFIntegrator::MFVectorConvectionNLFIntegrator(
#endif
}
MixedMFVectorConvectionNLIntegrator::MixedMFVectorConvectionNLIntegrator(
const VectorConvectionNLFIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q)
{
#ifdef MFEM_USE_CEED
NLConvectionOperatorInfo info(fes.GetMesh()->Dimension());
Assemble(integ, info, fes, Q);
#else
MFEM_ABORT("MFEM must be built with MFEM_USE_CEED=YES to use libCEED.");
#endif
}
} // namespace ceed
} // namespace mfem
@@ -13,7 +13,6 @@
#define MFEM_LIBCEED_NLCONV_HPP
#include "../../interface/integrator.hpp"
#include "../../interface/mixed_integrator.hpp"
#include "../../../fespace.hpp"
namespace mfem
@@ -32,15 +31,6 @@ public:
mfem::Coefficient *coeff);
};
class MixedPAVectorConvectionNLIntegrator : public MixedIntegrator<PAIntegrator>
{
public:
MixedPAVectorConvectionNLIntegrator(
const VectorConvectionNLFIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
};
/** Represent a VectorConvectionNLFIntegrator with AssemblyLevel::None
using libCEED. */
class MFVectorConvectionNLFIntegrator : public MFIntegrator
@@ -51,15 +41,6 @@ public:
mfem::Coefficient *coeff);
};
class MixedMFVectorConvectionNLIntegrator : public MixedIntegrator<MFIntegrator>
{
public:
MixedMFVectorConvectionNLIntegrator(
const VectorConvectionNLFIntegrator &integ,
const mfem::FiniteElementSpace &fes,
mfem::Coefficient *Q);
};
}
}
@@ -327,13 +327,13 @@ CEED_QFUNCTION(f_apply_conv_mf_const)(void *ctx, CeedInt Q,
const CeedScalar A33 = J11 * J22 - J12 * J21;
const CeedScalar w = qw[i] * coeff;
const CeedScalar qd00 = w * A11;
const CeedScalar qd10 = w * A21;
const CeedScalar qd20 = w * A31;
const CeedScalar qd01 = w * A12;
const CeedScalar qd01 = w * A21;
const CeedScalar qd02 = w * A31;
const CeedScalar qd10 = w * A12;
const CeedScalar qd11 = w * A22;
const CeedScalar qd21 = w * A32;
const CeedScalar qd02 = w * A13;
const CeedScalar qd12 = w * A23;
const CeedScalar qd12 = w * A32;
const CeedScalar qd20 = w * A13;
const CeedScalar qd21 = w * A23;
const CeedScalar qd22 = w * A33;
const CeedScalar u0 = u[i + Q * 0];
const CeedScalar u1 = u[i + Q * 1];
@@ -440,13 +440,13 @@ CEED_QFUNCTION(f_apply_conv_mf_quad)(void *ctx, CeedInt Q,
const CeedScalar A33 = J11 * J22 - J12 * J21;
const CeedScalar w = qw[i] * c[i];
const CeedScalar qd00 = w * A11;
const CeedScalar qd10 = w * A21;
const CeedScalar qd20 = w * A31;
const CeedScalar qd01 = w * A12;
const CeedScalar qd01 = w * A21;
const CeedScalar qd02 = w * A31;
const CeedScalar qd10 = w * A12;
const CeedScalar qd11 = w * A22;
const CeedScalar qd21 = w * A32;
const CeedScalar qd02 = w * A13;
const CeedScalar qd12 = w * A23;
const CeedScalar qd12 = w * A32;
const CeedScalar qd20 = w * A13;
const CeedScalar qd21 = w * A23;
const CeedScalar qd22 = w * A33;
const CeedScalar u0 = u[i + Q * 0];
const CeedScalar u1 = u[i + Q * 1];
+15 -37
View File
@@ -36,8 +36,6 @@ static CeedElemTopology GetCeedTopology(Geometry::Type geom)
return CEED_TOPOLOGY_HEX;
case Geometry::PRISM:
return CEED_TOPOLOGY_PRISM;
case Geometry::PYRAMID:
return CEED_TOPOLOGY_PYRAMID;
default:
MFEM_ABORT("This type of element is not supported");
return CEED_TOPOLOGY_PRISM; // Silence warning
@@ -45,11 +43,11 @@ static CeedElemTopology GetCeedTopology(Geometry::Type geom)
}
static void InitNonTensorBasis(const mfem::FiniteElementSpace &fes,
const mfem::FiniteElement &fe,
const mfem::IntegrationRule &ir,
Ceed ceed, CeedBasis *basis)
{
const mfem::DofToQuad &maps = fe.GetDofToQuad(ir, mfem::DofToQuad::FULL);
const mfem::DofToQuad &maps = fes.GetFE(0)->
GetDofToQuad(ir,mfem::DofToQuad::FULL);
mfem::Mesh *mesh = fes.GetMesh();
const int dim = mesh->Dimension();
const int ndofs = maps.ndof;
@@ -64,18 +62,18 @@ static void InitNonTensorBasis(const mfem::FiniteElementSpace &fes,
if (dim>2) { qX(2,i) = ip.z; }
qW(i) = ip.weight;
}
CeedBasisCreateH1(ceed, GetCeedTopology(fe.GetGeomType()),
CeedBasisCreateH1(ceed, GetCeedTopology(fes.GetFE(0)->GetGeomType()),
fes.GetVDim(), ndofs, nqpts,
maps.Bt.GetData(), maps.Gt.GetData(),
qX.GetData(), qW.GetData(), basis);
}
static void InitTensorBasis(const mfem::FiniteElementSpace &fes,
const mfem::FiniteElement &fe,
const mfem::IntegrationRule &ir,
Ceed ceed, CeedBasis *basis)
{
const mfem::DofToQuad &maps = fe.GetDofToQuad(ir, mfem::DofToQuad::TENSOR);
const mfem::DofToQuad &maps =
fes.GetFE(0)->GetDofToQuad(ir, mfem::DofToQuad::TENSOR);
mfem::Mesh *mesh = fes.GetMesh();
const int ndofs = maps.ndof;
const int nqpts = maps.nqpt;
@@ -98,30 +96,28 @@ static void InitTensorBasis(const mfem::FiniteElementSpace &fes,
qW.GetData(), basis);
}
static void InitBasisImpl(const FiniteElementSpace &fes,
const FiniteElement &fe,
const IntegrationRule &ir,
Ceed ceed, CeedBasis *basis)
void InitBasis(const FiniteElementSpace &fes,
const IntegrationRule &irm,
Ceed ceed, CeedBasis *basis)
{
// Check for FES -> basis, restriction in hash tables
const int P = fe.GetDof();
const int Q = ir.GetNPoints();
const mfem::FiniteElement *fe = fes.GetFE(0);
const int P = fe->GetDof();
const int Q = irm.GetNPoints();
const int ncomp = fes.GetVDim();
BasisKey basis_key(&fes, &ir, ncomp, P, Q);
BasisKey basis_key(&fes, &irm, ncomp, P, Q);
auto basis_itr = mfem::internal::ceed_basis_map.find(basis_key);
const bool tensor = dynamic_cast<const mfem::TensorBasisElement *>
(&fe) != nullptr;
// Init or retreive key values
if (basis_itr == mfem::internal::ceed_basis_map.end())
{
if ( tensor )
if (UsesTensorBasis(fes))
{
InitTensorBasis(fes, fe, ir, ceed, basis);
InitTensorBasis(fes, irm, ceed, basis);
}
else
{
InitNonTensorBasis(fes, fe, ir, ceed, basis);
InitNonTensorBasis(fes, irm, ceed, basis);
}
mfem::internal::ceed_basis_map[basis_key] = *basis;
}
@@ -131,24 +127,6 @@ static void InitBasisImpl(const FiniteElementSpace &fes,
}
}
void InitBasis(const FiniteElementSpace &fes,
const IntegrationRule &ir,
Ceed ceed, CeedBasis *basis)
{
const mfem::FiniteElement &fe = *fes.GetFE(0);
InitBasisImpl(fes, fe, ir, ceed, basis);
}
void InitBasisWithIndices(const FiniteElementSpace &fes,
const IntegrationRule &ir,
int nelem,
const int* indices,
Ceed ceed, CeedBasis *basis)
{
const mfem::FiniteElement &fe = *fes.GetFE(indices[0]);
InitBasisImpl(fes, fe, ir, ceed, basis);
}
#endif
} // namespace ceed
+3 -18
View File
@@ -22,32 +22,17 @@ namespace ceed
#ifdef MFEM_USE_CEED
/** @brief Initialize a CeedBasis for non-mixed meshes.
/** @brief Initialize a CeedBasis.
@param[in] fes Input finite element space.
@param[in] ir Input integration rule.
@param[in] irm Input integration rule.
@param[in] ceed Input Ceed object.
@param[out] basis The address of the initialized CeedBasis object.
*/
void InitBasis(const FiniteElementSpace &fes,
const IntegrationRule &ir,
const IntegrationRule &irm,
Ceed ceed, CeedBasis *basis);
/** @brief Initialize a CeedBasis for mixed meshes.
@param[in] fes The finite element space.
@param[in] ir is the integration rule for the operator.
@param[in] nelem The number of elements.
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`.
@param[in] ceed The Ceed object.
@param[out] basis The `CeedBasis` to initialize. */
void InitBasisWithIndices(const FiniteElementSpace &fes,
const IntegrationRule &ir,
int nelem,
const int* indices,
Ceed ceed, CeedBasis *basis);
#endif
} // namespace ceed
+3 -221
View File
@@ -14,7 +14,6 @@
#ifdef MFEM_USE_CEED
#include "../../../general/forall.hpp"
#include "../../../config/config.hpp"
#include "../../../linalg/vector.hpp"
#include "../../../linalg/dtensor.hpp"
@@ -78,14 +77,7 @@ struct QuadCoefficient : VariableCoefficient
/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
mfem::Coefficient @a Q, an mfem::Mesh @a mesh, and an mfem::IntegrationRule
@a ir.
@param[in] Q is the coefficient from the `Integrator`.
@param[in] mesh is the mesh.
@param[in] ir is the integration rule.
@param[out] coeff_ptr is the structure to store the coefficient for the
`CeedOperator`.
@param[out] ctx is the Context associated to the QFunction. */
@a ir. */
template <typename Context>
void InitCoefficient(mfem::Coefficient *Q, mfem::Mesh &mesh,
const mfem::IntegrationRule &ir,
@@ -151,15 +143,8 @@ void InitCoefficient(mfem::Coefficient *Q, mfem::Mesh &mesh,
/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
mfem::VectorCoefficient @a VQ, an mfem::Mesh @a mesh, and an
mfem::IntegrationRule @a ir.
@param[in] VQ is the vector coefficient from the `Integrator`.
@param[in] mesh is the mesh.
@param[in] ir is the integration rule.
@param[out] coeff_ptr is the structure to store the coefficient for the
`CeedOperator`.
@param[out] ctx is the Context associated to the QFunction. */
mfem::VectorCoefficient @a Q, an mfem::Mesh @a mesh, and an
mfem::IntegrationRule @a ir. */
template <typename Context>
void InitCoefficient(mfem::VectorCoefficient *VQ, mfem::Mesh &mesh,
const mfem::IntegrationRule &ir,
@@ -229,209 +214,6 @@ void InitCoefficient(mfem::VectorCoefficient *VQ, mfem::Mesh &mesh,
}
}
/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
mfem::Coefficient @a Q, an mfem::Mesh @a mesh, and an mfem::IntegrationRule
@a ir for the elements given by the indices @a indices.
@param[in] Q is the coefficient from the `Integrator`.
@param[in] mesh is the mesh.
@param[in] ir is the integration rule.
@param[in] nelem The number of elements.
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`.
@param[out] coeff_ptr is the structure to store the coefficient for the
`CeedOperator`.
@param[out] ctx is the Context associated to the QFunction. */
template <typename Context>
void InitCoefficientWithIndices(mfem::Coefficient *Q, mfem::Mesh &mesh,
const mfem::IntegrationRule &ir,
int nelem,
const int* indices,
Coefficient*& coeff_ptr, Context &ctx)
{
if ( Q == nullptr )
{
Coefficient *ceedCoeff = new Coefficient(1);
ctx.coeff = 1.0;
coeff_ptr = ceedCoeff;
}
else if (ConstantCoefficient *const_coeff =
dynamic_cast<ConstantCoefficient*>(Q))
{
Coefficient *ceedCoeff = new Coefficient(1);
ctx.coeff = const_coeff->constant;
coeff_ptr = ceedCoeff;
}
else if (GridFunctionCoefficient* gf_coeff =
dynamic_cast<GridFunctionCoefficient*>(Q))
{
GridCoefficient *ceedCoeff =
new GridCoefficient(*gf_coeff->GetGridFunction());
coeff_ptr = ceedCoeff;
}
else if (QuadratureFunctionCoefficient *cQ =
dynamic_cast<QuadratureFunctionCoefficient*>(Q))
{
QuadCoefficient *ceedCoeff = new QuadCoefficient(1);
const int ne = mesh.GetNE();
const int nq = ir.GetNPoints();
const mfem::QuadratureFunction &qFun = cQ->GetQuadFunction();
MFEM_VERIFY(qFun.Size() == nq * ne,
"Incompatible QuadratureFunction dimension \n");
MFEM_VERIFY(&ir == &qFun.GetSpace()->GetElementIntRule(0),
"IntegrationRule used within integrator and in"
" QuadratureFunction appear to be different");
ceedCoeff->coeff.SetSize(nq * nelem);
Memory<int> m_indices((int*)indices, nelem, false);
auto in = Reshape(qFun.Read(), nq, ne);
auto d_indices = Read(m_indices, nelem);
auto out = Reshape(ceedCoeff->coeff.Write(), nq, nelem);
MFEM_FORALL(i, nelem * nq,
{
const int q = i%nq;
const int sub_e = i/nq;
const int e = d_indices[sub_e];
out(q, sub_e) = in(q, e);
});
m_indices.DeleteDevice();
InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
coeff_ptr = ceedCoeff;
}
else
{
QuadCoefficient *ceedCoeff = new QuadCoefficient(1);
const int nq = ir.GetNPoints();
ceedCoeff->coeff.SetSize(nq * nelem);
auto C = Reshape(ceedCoeff->coeff.HostWrite(), nq, nelem);
for (int i = 0; i < nelem; ++i)
{
const int e = indices[i];
mfem::ElementTransformation &T = *mesh.GetElementTransformation(e);
for (int q = 0; q < nq; ++q)
{
C(q, i) = Q->Eval(T, ir.IntPoint(q));
}
}
InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
coeff_ptr = ceedCoeff;
}
}
/** @brief Initializes an mfem::ceed::Coefficient @a coeff_ptr from an
mfem::VectorCoefficient @a Q, an mfem::Mesh @a mesh, and an
mfem::IntegrationRule @a ir for the elements given by the indices @a indices.
@param[in] VQ is the vector coefficient from the `Integrator`.
@param[in] mesh is the mesh.
@param[in] ir is the integration rule.
@param[in] nelem The number of elements.
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`.
@param[out] coeff_ptr is the structure to store the coefficient for the
`CeedOperator`.
@param[out] ctx is the Context associated to the QFunction. */
template <typename Context>
void InitCoefficientWithIndices(mfem::VectorCoefficient *VQ, mfem::Mesh &mesh,
const mfem::IntegrationRule &ir,
int nelem,
const int* indices,
Coefficient *&coeff_ptr, Context &ctx)
{
if (VectorConstantCoefficient *const_coeff =
dynamic_cast<VectorConstantCoefficient*>(VQ))
{
const int vdim = const_coeff->GetVDim();
const mfem::Vector &val = const_coeff->GetVec();
Coefficient *ceedCoeff = new Coefficient(vdim);
for (int i = 0; i < vdim; i++)
{
ctx.coeff[i] = val[i];
}
coeff_ptr = ceedCoeff;
}
else if (VectorGridFunctionCoefficient* vgf_coeff =
dynamic_cast<VectorGridFunctionCoefficient*>(VQ))
{
GridCoefficient *ceedCoeff =
new GridCoefficient(*vgf_coeff->GetGridFunction());
coeff_ptr = ceedCoeff;
}
else if (VectorQuadratureFunctionCoefficient *cQ =
dynamic_cast<VectorQuadratureFunctionCoefficient*>(VQ))
{
QuadCoefficient *ceedCoeff = new QuadCoefficient(cQ->GetVDim());
const int dim = mesh.Dimension();
const int ne = mesh.GetNE();
const int nq = ir.GetNPoints();
const mfem::QuadratureFunction &qFun = cQ->GetQuadFunction();
MFEM_VERIFY(qFun.Size() == dim * nq * ne,
"Incompatible QuadratureFunction dimension \n");
MFEM_VERIFY(&ir == &qFun.GetSpace()->GetElementIntRule(0),
"IntegrationRule used within integrator and in"
" QuadratureFunction appear to be different");
ceedCoeff->coeff.SetSize(dim * nq * nelem);
Memory<int> m_indices((int*)indices, nelem, false);
auto in = Reshape(qFun.Read(), dim, nq, ne);
auto d_indices = Read(m_indices, nelem);
auto out = Reshape(ceedCoeff->coeff.Write(), dim, nq, nelem);
MFEM_FORALL(i, nelem * nq,
{
const int q = i%nq;
const int sub_e = i/nq;
const int e = d_indices[sub_e];
for (int d = 0; d < dim; d++)
{
out(d, q, sub_e) = in(d, q, e);
}
});
m_indices.DeleteDevice();
InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
coeff_ptr = ceedCoeff;
}
else
{
const int dim = mesh.Dimension();
QuadCoefficient *ceedCoeff = new QuadCoefficient(dim);
const int nq = ir.GetNPoints();
ceedCoeff->coeff.SetSize(dim * nq * nelem);
auto C = Reshape(ceedCoeff->coeff.HostWrite(), dim, nq, nelem);
mfem::DenseMatrix Q_ir;
for (int i = 0; i < nelem; ++i)
{
const int e = indices[i];
mfem::ElementTransformation &T = *mesh.GetElementTransformation(e);
VQ->Eval(Q_ir, T, ir);
for (int q = 0; q < nq; ++q)
{
for (int d = 0; d < dim; ++d)
{
C(d, q, i) = Q_ir(d, q);
}
}
}
InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
coeff_ptr = ceedCoeff;
}
}
template <typename Coeff, typename Context>
void InitCoefficient(Coeff *Q, mfem::Mesh &mesh,
const mfem::IntegrationRule &ir, int nelem,
const int* indices, Coefficient *&coeff_ptr, Context &ctx)
{
if (indices)
{
InitCoefficientWithIndices(Q, mesh, ir, nelem, indices, coeff_ptr, ctx);
}
else
{
InitCoefficient(Q, mesh, ir, coeff_ptr, ctx);
}
}
} // namespace ceed
} // namespace mfem
+83 -182
View File
@@ -18,7 +18,6 @@
#include "operator.hpp"
#include "coefficient.hpp"
#include "restriction.hpp"
#include "util.hpp"
#include "ceed.hpp"
namespace mfem
@@ -87,7 +86,6 @@ protected:
CeedQFunctionContext build_ctx;
CeedOperator build_oper;
public:
PAIntegrator()
: Operator(),
trial_basis(nullptr), test_basis(nullptr), mesh_basis(nullptr),
@@ -97,51 +95,23 @@ public:
qdata(nullptr), coeff(nullptr), build_ctx(nullptr), build_oper(nullptr)
{ }
/** @brief This method assembles the `PAIntegrator` with the given
`CeedOperatorInfo` @a info, an `mfem::FiniteElementSpace` @a fes, an
`mfem::IntegrationRule` @a ir, and `mfem::Coefficient` or
`mfem::VectorCoefficient` @a Q.
The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
and contain a `Context` type relevant to the qFunctions.
public:
/** This method assembles the PAIntegrator.
@param[in] info is the structure describing the CeedOperator to assemble.
@param[in] fes is the finite element space.
@param[in] ir is the integration rule for the operator.
@param[in] Q is the coefficient from the `Integrator`. */
@param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
the `CeedOperatorInfo` type is expected to inherit from
`OperatorInfo` and contain a `Context` type relevant to
the qFunctions.
@param[in] fes the `FiniteElementSpace` for the form,
@param[in] ir the `IntegrationRule` for the numerical integration,
@param[in] Q `Coefficient` or `VectorCoefficient`. */
template <typename CeedOperatorInfo, typename CoeffType>
void Assemble(CeedOperatorInfo &info,
const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
const mfem::IntegrationRule &irm,
CoeffType *Q)
{
Assemble(info, fes, ir, fes.GetNE(), nullptr, Q);
}
/** @brief This method assembles the `PAIntegrator` with the given
`CeedOperatorInfo` @a info, an `mfem::FiniteElementSpace` @a fes, an
`mfem::IntegrationRule` @a ir, and `mfem::Coefficient` or
`mfem::VectorCoefficient` @a Q for the elements given by the indices
@a indices.
The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
and contain a `Context` type relevant to the qFunctions.
@param[in] info is the structure describing the CeedOperator to assemble.
@param[in] fes is the finite element space.
@param[in] ir is the integration rule for the operator.
@param[in] nelem The number of elements.
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`. If `indices == nullptr`, assumes
that the `FiniteElementSpace` is not mixed.
@param[in] Q is the coefficient from the `Integrator`. */
template <typename CeedOperatorInfo, typename CoeffType>
void Assemble(CeedOperatorInfo &info,
const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
int nelem,
const int* indices,
CoeffType *Q)
{
Assemble(info, fes, fes, ir, nelem, indices, Q);
Assemble(info, fes, fes, irm, Q);
}
/** This method assembles the PAIntegrator for mixed forms.
@@ -158,40 +128,12 @@ public:
void Assemble(CeedOperatorInfo &info,
const mfem::FiniteElementSpace &trial_fes,
const mfem::FiniteElementSpace &test_fes,
const mfem::IntegrationRule &ir,
CoeffType *Q)
{
Assemble(info, trial_fes, test_fes, ir, trial_fes.GetNE(), nullptr, Q);
}
/** This method assembles the PAIntegrator for mixed forms on mixed meshes.
@param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
the `CeedOperatorInfo` type is expected to inherit from
`OperatorInfo` and contain a `Context` type relevant to
the qFunctions.
@param[in] trial_fes the trial `FiniteElementSpace` for the form,
@param[in] test_fes the test `FiniteElementSpace` for the form,
@param[in] ir the `IntegrationRule` for the numerical integration,
@param[in] nelem The number of elements,
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`. If `indices == nullptr`, assumes
that the `FiniteElementSpace` is not mixed,
@param[in] Q `Coefficient` or `VectorCoefficient`. */
template <typename CeedOperatorInfo, typename CoeffType>
void Assemble(CeedOperatorInfo &info,
const mfem::FiniteElementSpace &trial_fes,
const mfem::FiniteElementSpace &test_fes,
const mfem::IntegrationRule &ir,
int nelem,
const int* indices,
const mfem::IntegrationRule &irm,
CoeffType *Q)
{
Ceed ceed(internal::ceed);
mfem::Mesh &mesh = *trial_fes.GetMesh();
MFEM_VERIFY(!(!indices && mesh.GetNumGeometries(mesh.Dimension()) > 1),
"Use ceed::MixedIntegrator on mixed meshes.");
InitCoefficient(Q, mesh, ir, nelem, indices, coeff, info.ctx);
InitCoefficient(Q, mesh, irm, coeff, info.ctx);
bool const_coeff = coeff->IsConstant();
std::string build_func = const_coeff ? info.build_func_const
: info.build_func_quad;
@@ -203,6 +145,7 @@ public:
info.trial_op,
info.test_op
};
CeedInt nqpts, nelem = mesh.GetNE();
CeedInt dim = mesh.SpaceDimension();
CeedInt trial_vdim = trial_fes.GetVDim();
CeedInt test_vdim = test_fes.GetVDim();
@@ -210,23 +153,23 @@ public:
mesh.EnsureNodes();
if ( &trial_fes == &test_fes )
{
InitBasisAndRestriction(trial_fes, ir, nelem, indices,
ceed, &trial_basis, &trial_restr);
InitBasisAndRestriction(trial_fes, irm, ceed,
&trial_basis, &trial_restr);
test_basis = trial_basis;
test_restr = trial_restr;
}
else
{
InitBasisAndRestriction(trial_fes, ir, nelem, indices,
ceed, &trial_basis, &trial_restr);
InitBasisAndRestriction(test_fes, ir, nelem, indices,
ceed, &test_basis, &test_restr);
InitBasisAndRestriction(trial_fes, irm, ceed,
&trial_basis, &trial_restr);
InitBasisAndRestriction(test_fes, irm, ceed,
&test_basis, &test_restr);
}
const mfem::FiniteElementSpace *mesh_fes = mesh.GetNodalFESpace();
MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
InitBasisAndRestriction(*mesh_fes, ir, nelem, indices,
ceed, &mesh_basis, &mesh_restr);
InitBasisAndRestriction(*mesh_fes, irm, ceed, &mesh_basis,
&mesh_restr);
CeedInt trial_nqpts, test_nqpts;
CeedBasisGetNumQuadraturePoints(trial_basis, &trial_nqpts);
@@ -234,7 +177,7 @@ public:
MFEM_VERIFY(trial_nqpts == test_nqpts,
"Trial and test basis must have the same number of quadrature"
" points.");
CeedInt nqpts = trial_nqpts;
nqpts = trial_nqpts;
const int qdatasize = op.qdatasize;
InitStridedRestriction(*mesh_fes, nelem, nqpts, qdatasize,
@@ -278,10 +221,8 @@ public:
CeedOperatorCreate(ceed, build_qfunc, NULL, NULL, &build_oper);
if (GridCoefficient *gridCoeff = dynamic_cast<GridCoefficient*>(coeff))
{
InitBasisAndRestriction(*gridCoeff->gf.FESpace(), ir,
nelem, indices, ceed,
&gridCoeff->basis,
&gridCoeff->restr);
InitBasisAndRestriction(*gridCoeff->gf.FESpace(), irm, ceed,
&gridCoeff->basis, &gridCoeff->restr);
CeedOperatorSetField(build_oper, "coeff", gridCoeff->restr,
gridCoeff->basis, gridCoeff->coeffVector);
}
@@ -290,8 +231,7 @@ public:
{
const int ncomp = quadCoeff->ncomp;
CeedInt strides[3] = {ncomp, 1, ncomp*nqpts};
InitStridedRestriction(*mesh.GetNodalFESpace(),
nelem, nqpts, ncomp, strides,
InitStridedRestriction(*mesh_fes, nelem, nqpts, ncomp, strides,
&quadCoeff->restr);
CeedOperatorSetField(build_oper, "coeff", quadCoeff->restr,
CEED_BASIS_COLLOCATED, quadCoeff->coeffVector);
@@ -314,17 +254,22 @@ public:
switch (op.trial_op)
{
case EvalMode::None:
CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim, CEED_EVAL_NONE);
CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
CEED_EVAL_NONE);
break;
case EvalMode::Interp:
CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim, CEED_EVAL_INTERP);
CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
CEED_EVAL_INTERP);
break;
case EvalMode::Grad:
CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim, CEED_EVAL_GRAD);
CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
CEED_EVAL_GRAD);
break;
case EvalMode::InterpAndGrad:
CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim, CEED_EVAL_INTERP);
CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim, CEED_EVAL_GRAD);
CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
CEED_EVAL_INTERP);
CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
CEED_EVAL_GRAD);
break;
}
// qdata
@@ -333,17 +278,22 @@ public:
switch (op.test_op)
{
case EvalMode::None:
CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim, CEED_EVAL_NONE);
CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
CEED_EVAL_NONE);
break;
case EvalMode::Interp:
CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim, CEED_EVAL_INTERP);
CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
CEED_EVAL_INTERP);
break;
case EvalMode::Grad:
CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim, CEED_EVAL_GRAD);
CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
CEED_EVAL_GRAD);
break;
case EvalMode::InterpAndGrad:
CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim, CEED_EVAL_INTERP);
CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim, CEED_EVAL_GRAD);
CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
CEED_EVAL_INTERP);
CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
CEED_EVAL_GRAD);
break;
}
CeedQFunctionSetContext(apply_qfunc, build_ctx);
@@ -358,14 +308,18 @@ public:
CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
break;
case EvalMode::Interp:
CeedOperatorSetField(oper, "u", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
CEED_VECTOR_ACTIVE);
break;
case EvalMode::Grad:
CeedOperatorSetField(oper, "gu", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
CEED_VECTOR_ACTIVE);
break;
case EvalMode::InterpAndGrad:
CeedOperatorSetField(oper, "u", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "gu", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
CEED_VECTOR_ACTIVE);
break;
}
// qdata
@@ -379,14 +333,18 @@ public:
CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
break;
case EvalMode::Interp:
CeedOperatorSetField(oper, "v", test_restr, test_basis, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "v", test_restr, test_basis,
CEED_VECTOR_ACTIVE);
break;
case EvalMode::Grad:
CeedOperatorSetField(oper, "gv", test_restr, test_basis, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "gv", test_restr, test_basis,
CEED_VECTOR_ACTIVE);
break;
case EvalMode::InterpAndGrad:
CeedOperatorSetField(oper, "v", test_restr, test_basis, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "gv", test_restr, test_basis, CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "v", test_restr, test_basis,
CEED_VECTOR_ACTIVE);
CeedOperatorSetField(oper, "gv", test_restr, test_basis,
CEED_VECTOR_ACTIVE);
break;
}
@@ -444,7 +402,6 @@ protected:
Coefficient *coeff;
CeedQFunctionContext build_ctx;
public:
MFIntegrator()
: Operator(),
trial_basis(nullptr), test_basis(nullptr), mesh_basis(nullptr),
@@ -453,51 +410,23 @@ public:
apply_qfunc(nullptr), node_coords(nullptr),
qdata(nullptr), coeff(nullptr), build_ctx(nullptr) { }
/** @brief This method assembles the `MFIntegrator` with the given
`CeedOperatorInfo` @a info, an `mfem::FiniteElementSpace` @a fes, an
`mfem::IntegrationRule` @a ir, and `mfem::Coefficient` or
`mfem::VectorCoefficient` @a Q.
The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
and contain a `Context` type relevant to the qFunctions.
public:
/** This method assembles the MFIntegrator.
@param[in] info is the structure describing the CeedOperator to assemble.
@param[in] fes is the finite element space.
@param[in] ir is the integration rule for the operator.
@param[in] Q is the coefficient from the `Integrator`. */
@param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
the `CeedOperatorInfo` type is expected to inherit from
`OperatorInfo` and contain a `Context` type relevant to
the qFunctions.
@param[in] fes the `FiniteElementSpace` for the form,
@param[in] ir the `IntegrationRule` for the numerical integration,
@param[in] Q `Coefficient` or `VectorCoefficient`. */
template <typename CeedOperatorInfo, typename CoeffType>
void Assemble(CeedOperatorInfo &info,
const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
const mfem::IntegrationRule &irm,
CoeffType *Q)
{
Assemble(info, fes, ir, fes.GetNE(), nullptr, Q);
}
/** @brief This method assembles the `MFIntegrator` with the given
`CeedOperatorInfo` @a info, an `mfem::FiniteElementSpace` @a fes, an
`mfem::IntegrationRule` @a ir, and `mfem::Coefficient` or
`mfem::VectorCoefficient` @a Q for the elements given by the indices
@a indices.
The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
and contain a `Context` type relevant to the qFunctions.
@param[in] info is the structure describing the CeedOperator to assemble.
@param[in] fes is the finite element space.
@param[in] ir is the integration rule for the operator.
@param[in] nelem The number of elements.
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`. If `indices == nullptr`, assumes
that the `FiniteElementSpace` is not mixed.
@param[in] Q is the coefficient from the `Integrator`. */
template <typename CeedOperatorInfo, typename CoeffType>
void Assemble(CeedOperatorInfo &info,
const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
int nelem,
const int* indices,
CoeffType *Q)
{
Assemble(info, fes, fes, ir, nelem, indices, Q);
Assemble(info, fes, fes, irm, Q);
}
/** This method assembles the MFIntegrator for mixed forms.
@@ -514,40 +443,12 @@ public:
void Assemble(CeedOperatorInfo &info,
const mfem::FiniteElementSpace &trial_fes,
const mfem::FiniteElementSpace &test_fes,
const mfem::IntegrationRule &ir,
CoeffType *Q)
{
Assemble(info, trial_fes, test_fes, ir, trial_fes.GetNE(), nullptr, Q);
}
/** This method assembles the MFIntegrator for mixed forms.
@param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
the `CeedOperatorInfo` type is expected to inherit from
`OperatorInfo` and contain a `Context` type relevant to
the qFunctions.
@param[in] trial_fes the trial `FiniteElementSpace` for the form,
@param[in] test_fes the test `FiniteElementSpace` for the form,
@param[in] ir the `IntegrationRule` for the numerical integration,
@param[in] nelem The number of elements,
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`. If `indices == nullptr`, assumes
that the `FiniteElementSpace` is not mixed,
@param[in] Q `Coefficient` or `VectorCoefficient`. */
template <typename CeedOperatorInfo, typename CoeffType>
void Assemble(CeedOperatorInfo &info,
const mfem::FiniteElementSpace &trial_fes,
const mfem::FiniteElementSpace &test_fes,
const mfem::IntegrationRule &ir,
int nelem,
const int* indices,
const mfem::IntegrationRule &irm,
CoeffType *Q)
{
Ceed ceed(internal::ceed);
Mesh &mesh = *trial_fes.GetMesh();
MFEM_VERIFY(!(!indices && mesh.GetNumGeometries(mesh.Dimension()) > 1),
"Use ceed::MixedIntegrator on mixed meshes.");
InitCoefficient(Q, mesh, ir, nelem, indices, coeff, info.ctx);
InitCoefficient(Q, mesh, irm, coeff, info.ctx);
bool const_coeff = coeff->IsConstant();
std::string apply_func = const_coeff ? info.apply_func_mf_const
: info.apply_func_mf_quad;
@@ -558,7 +459,7 @@ public:
info.trial_op,
info.test_op
};
CeedInt nqpts, nelem = mesh.GetNE();
CeedInt dim = mesh.SpaceDimension();
CeedInt trial_vdim = trial_fes.GetVDim();
CeedInt test_vdim = test_fes.GetVDim();
@@ -566,22 +467,22 @@ public:
mesh.EnsureNodes();
if ( &trial_fes == &test_fes )
{
InitBasisAndRestriction(trial_fes, ir, nelem, indices, ceed,
InitBasisAndRestriction(trial_fes, irm, ceed,
&trial_basis, &trial_restr);
test_basis = trial_basis;
test_restr = trial_restr;
}
else
{
InitBasisAndRestriction(trial_fes, ir, nelem, indices, ceed,
InitBasisAndRestriction(trial_fes, irm, ceed,
&trial_basis, &trial_restr);
InitBasisAndRestriction(test_fes, ir, nelem, indices, ceed,
InitBasisAndRestriction(test_fes, irm, ceed,
&test_basis, &test_restr);
}
const mfem::FiniteElementSpace *mesh_fes = mesh.GetNodalFESpace();
MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
InitBasisAndRestriction(*mesh_fes, ir, nelem, indices, ceed, &mesh_basis,
InitBasisAndRestriction(*mesh_fes, irm, ceed, &mesh_basis,
&mesh_restr);
CeedInt trial_nqpts, test_nqpts;
@@ -590,7 +491,7 @@ public:
MFEM_VERIFY(trial_nqpts == test_nqpts,
"Trial and test basis must have the same number of quadrature"
" points.");
CeedInt nqpts = trial_nqpts;
nqpts = trial_nqpts;
InitVector(*mesh.GetNodes(), node_coords);
@@ -671,8 +572,8 @@ public:
// coefficient
if (GridCoefficient *gridCoeff = dynamic_cast<GridCoefficient*>(coeff))
{
InitBasisAndRestriction(*gridCoeff->gf.FESpace(), ir, nelem, indices,
ceed, &gridCoeff->basis, &gridCoeff->restr);
InitBasisAndRestriction(*gridCoeff->gf.FESpace(), irm, ceed,
&gridCoeff->basis, &gridCoeff->restr);
CeedOperatorSetField(oper, "coeff", gridCoeff->restr,
gridCoeff->basis, gridCoeff->coeffVector);
}
-2
View File
@@ -22,8 +22,6 @@
#include "coefficient.hpp"
// PA or MF Operator using libCEED.
#include "integrator.hpp"
// PA Operator supporting mixed finite element spaces.
#include "mixed_integrator.hpp"
// Utility functions
#include "util.hpp"
// Wrapper to include <ceed.h>
-126
View File
@@ -1,126 +0,0 @@
// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#ifndef MFEM_LIBCEED_MIXED_INTEGRATOR
#define MFEM_LIBCEED_MIXED_INTEGRATOR
#include "ceed.hpp"
#include "integrator.hpp"
#include <unordered_map>
namespace mfem
{
namespace ceed
{
/** @brief This class wraps a `ceed::PAIntegrator` or `ceed::MFIntegrator` to
support mixed finite element spaces. */
template <typename CeedInteg>
class MixedIntegrator : public ceed::Operator
{
#ifdef MFEM_USE_CEED
using ElementKey = std::pair<int, int>; //< Element::Type, Order >
struct key_hash
{
std::size_t operator()(const ElementKey& k) const
{
return k.first + 2 * k.second;
}
};
using ElementsMap = std::unordered_map<const ElementKey, int*, key_hash>;
std::vector<CeedInteg*> sub_ops;
public:
template <typename Integrator, typename CeedOperatorInfo, typename CoeffType>
void Assemble(const Integrator &integ,
CeedOperatorInfo &info,
const mfem::FiniteElementSpace &fes,
CoeffType *Q)
{
ElementsMap count;
ElementsMap element_indices;
ElementsMap offsets;
// Count the number of elements of each type
for (int i = 0; i < fes.GetNE(); i++)
{
ElementKey key(fes.GetElementType(i), fes.GetElementOrder(i));
auto value = count.find(key);
if (value == count.end())
{
count[key] = new int(1);
}
else
{
(*value->second)++;
}
}
// Initialization of the arrays
for ( const auto& value : count )
{
element_indices[value.first] = new int[*value.second];
offsets[value.first] = new int(0);
}
// Populates the indices arrays for each element type
for (int i = 0; i < fes.GetNE(); i++)
{
ElementKey key(fes.GetElementType(i), fes.GetElementOrder(i));
int &offset = *(offsets[key]);
int* indices_array = element_indices[key];
indices_array[offset] = i;
offset++;
}
// Create composite CeedOperator
CeedCompositeOperatorCreate(internal::ceed, &oper);
// Create each sub-CeedOperator
sub_ops.reserve(element_indices.size());
for (const auto& value : element_indices)
{
const int* indices = value.second;
const int first_index = indices[0];
const mfem::FiniteElement &el = *fes.GetFE(first_index);
auto &T = *fes.GetMesh()->GetElementTransformation(first_index);
MFEM_ASSERT(!integ.GetIntegrationRule(),
"Mixed mesh integrators should not have an"
" IntegrationRule.");
const IntegrationRule &ir = GetRule(integ, el, el, T);
auto sub_op = new CeedInteg();
int nelem = *count[value.first];
sub_op->Assemble(info, fes, ir, nelem, indices, Q);
sub_ops.push_back(sub_op);
CeedCompositeOperatorAddSub(oper, sub_op->GetCeedOperator());
}
const int ndofs = fes.GetVDim() * fes.GetNDofs();
CeedVectorCreate(internal::ceed, ndofs, &u);
CeedVectorCreate(internal::ceed, ndofs, &v);
}
virtual ~MixedIntegrator()
{
for (auto sub_op : sub_ops)
{
delete sub_op;
}
}
#endif
};
} // namespace ceed
} // namespace mfem
#endif // MFEM_LIBCEED_MIXED_INTEGRATOR
+56 -195
View File
@@ -20,8 +20,8 @@ namespace ceed
#ifdef MFEM_USE_CEED
static void InitNativeRestr(const mfem::FiniteElementSpace &fes,
Ceed ceed, CeedElemRestriction *restr)
static void InitNonTensorRestriction(const mfem::FiniteElementSpace &fes,
Ceed ceed, CeedElemRestriction *restr)
{
const mfem::FiniteElement *fe = fes.GetFE(0);
const int P = fe->GetDof();
@@ -31,173 +31,77 @@ static void InitNativeRestr(const mfem::FiniteElementSpace &fes,
const mfem::TensorBasisElement * tfe =
dynamic_cast<const mfem::TensorBasisElement *>(fe);
const int stride = compstride == 1 ? fes.GetVDim() : 1;
const mfem::Array<int>& dof_map = tfe->GetDofMap();
for (int i = 0; i < fes.GetNE(); i++)
if (tfe) // Lexicographic ordering using dof_map
{
const int el_offset = P * i;
for (int j = 0; j < P; j++)
const mfem::Array<int>& dof_map = tfe->GetDofMap();
for (int i = 0; i < fes.GetNE(); i++)
{
tp_el_dof[j+el_offset] = stride*el_dof.GetJ()[dof_map[j]+el_offset];
const int el_offset = P * i;
for (int j = 0; j < P; j++)
{
tp_el_dof[j+el_offset] = stride*el_dof.GetJ()[dof_map[j]+el_offset];
}
}
}
else // Native ordering
{
for (int e = 0; e < fes.GetNE(); e++)
{
for (int i = 0; i < P; i++)
{
tp_el_dof[i + e*P] = stride*el_dof.GetJ()[i + e*P];
}
}
}
CeedElemRestrictionCreate(ceed, fes.GetNE(), P, fes.GetVDim(),
compstride, (fes.GetVDim())*(fes.GetNDofs()),
CEED_MEM_HOST, CEED_COPY_VALUES,
tp_el_dof.GetData(), restr);
}
static void InitLexicoRestr(const mfem::FiniteElementSpace &fes,
Ceed ceed, CeedElemRestriction *restr)
// TODO fuse Tensor and NonTensor Restriction
void InitTensorRestriction(const mfem::FiniteElementSpace &fes,
Ceed ceed, CeedElemRestriction *restr)
{
const mfem::FiniteElement *fe = fes.GetFE(0);
const int P = fe->GetDof();
const mfem::TensorBasisElement * tfe =
dynamic_cast<const mfem::TensorBasisElement *>(fe);
MFEM_VERIFY(tfe, "invalid FE");
const mfem::Array<int>& dof_map = tfe->GetDofMap();
CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
const mfem::Table &el_dof = fes.GetElementToDofTable();
mfem::Array<int> tp_el_dof(el_dof.Size_of_connections());
const int dof = fe->GetDof();
const int stride = compstride == 1 ? fes.GetVDim() : 1;
for (int e = 0; e < fes.GetNE(); e++)
if (dof_map.Size()>0)
{
for (int i = 0; i < P; i++)
for (int i = 0; i < fes.GetNE(); i++)
{
tp_el_dof[i + e*P] = stride*el_dof.GetJ()[i + e*P];
const int el_offset = dof * i;
for (int j = 0; j < dof; j++)
{
tp_el_dof[j+el_offset] = stride*el_dof.GetJ()[dof_map[j]+el_offset];
}
}
}
CeedElemRestrictionCreate(ceed, fes.GetNE(), P, fes.GetVDim(),
else // dof_map.Size == 0, means dof_map[j]==j;
{
for (int i = 0; i < fes.GetNE(); i++)
{
const int el_offset = dof * i;
for (int j = 0; j < dof; j++)
{
tp_el_dof[j+el_offset] = stride*el_dof.GetJ()[j+el_offset];
}
}
}
CeedElemRestrictionCreate(ceed, fes.GetNE(), dof, fes.GetVDim(),
compstride, (fes.GetVDim())*(fes.GetNDofs()),
CEED_MEM_HOST, CEED_COPY_VALUES,
tp_el_dof.GetData(), restr);
}
static void InitRestrictionImpl(const mfem::FiniteElementSpace &fes,
Ceed ceed, CeedElemRestriction *restr)
{
const mfem::FiniteElement *fe = fes.GetFE(0);
const mfem::TensorBasisElement * tfe =
dynamic_cast<const mfem::TensorBasisElement *>(fe);
if ( tfe && tfe->GetDofMap().Size()>0 ) // Native ordering using dof_map
{
InitNativeRestr(fes, ceed, restr);
}
else // Lexicographic ordering
{
InitLexicoRestr(fes, ceed, restr);
}
}
static void InitNativeRestrWithIndices(
const mfem::FiniteElementSpace &fes,
int nelem,
const int* indices,
Ceed ceed, CeedElemRestriction *restr)
{
const mfem::FiniteElement *fe = fes.GetFE(indices[0]);
const int P = fe->GetDof();
CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
mfem::Array<int> tp_el_dof(nelem*P);
const mfem::TensorBasisElement * tfe =
dynamic_cast<const mfem::TensorBasisElement *>(fe);
Array<int> dofs;
const int stride = compstride == 1 ? fes.GetVDim() : 1;
const mfem::Array<int>& dof_map = tfe->GetDofMap();
for (int i = 0; i < nelem; i++)
{
const int elem_index = indices[i];
fes.GetElementDofs(elem_index, dofs);
const int el_offset = P * i;
for (int j = 0; j < P; j++)
{
tp_el_dof[j + el_offset] = stride*dofs[dof_map[j]];
}
}
CeedElemRestrictionCreate(ceed, nelem, P, fes.GetVDim(),
compstride, (fes.GetVDim())*(fes.GetNDofs()),
CEED_MEM_HOST, CEED_COPY_VALUES,
tp_el_dof.GetData(), restr);
}
static void InitLexicoRestrWithIndices(
const mfem::FiniteElementSpace &fes,
int nelem,
const int* indices,
Ceed ceed, CeedElemRestriction *restr)
{
const mfem::FiniteElement *fe = fes.GetFE(indices[0]);
const int P = fe->GetDof();
CeedInt compstride = fes.GetOrdering()==Ordering::byVDIM ? 1 : fes.GetNDofs();
mfem::Array<int> tp_el_dof(nelem*P);
Array<int> dofs;
const int stride = compstride == 1 ? fes.GetVDim() : 1;
for (int i = 0; i < nelem; i++)
{
const int elem_index = indices[i];
fes.GetElementDofs(elem_index, dofs);
const int el_offset = P * i;
for (int j = 0; j < P; j++)
{
tp_el_dof[j + el_offset] = stride*dofs[j];
}
}
CeedElemRestrictionCreate(ceed, nelem, P, fes.GetVDim(),
compstride, (fes.GetVDim())*(fes.GetNDofs()),
CEED_MEM_HOST, CEED_COPY_VALUES,
tp_el_dof.GetData(), restr);
}
static void InitRestrictionWithIndicesImpl(
const mfem::FiniteElementSpace &fes,
int nelem,
const int* indices,
Ceed ceed, CeedElemRestriction *restr)
{
const mfem::FiniteElement *fe = fes.GetFE(indices[0]);
const mfem::TensorBasisElement * tfe =
dynamic_cast<const mfem::TensorBasisElement *>(fe);
if ( tfe && tfe->GetDofMap().Size()>0 ) // Native ordering using dof_map
{
InitNativeRestrWithIndices(fes, nelem, indices, ceed, restr);
}
else // Lexicographic ordering
{
InitLexicoRestrWithIndices(fes, nelem, indices, ceed, restr);
}
}
static void InitCoeffRestrictionWithIndicesImpl(
const mfem::FiniteElementSpace &fes,
int nelem,
const int* indices,
int nquads,
int ncomp,
Ceed ceed,
CeedElemRestriction *restr)
{
mfem::Array<int> tp_el_dof(nelem*nquads);
const int stride_quad = ncomp;
const int stride_elem = ncomp*nquads;
// TODO generalize to support different #quads
for (int i = 0; i < nelem; i++)
{
const int elem_index = indices[i];
const int el_offset = elem_index * stride_elem;
for (int j = 0; j < nquads; j++)
{
tp_el_dof[j + nquads * i] = j * stride_quad + el_offset;
}
}
CeedElemRestrictionCreate(ceed, nelem, nquads, ncomp, 1,
ncomp*fes.GetNE()*nquads,
CEED_MEM_HOST, CEED_COPY_VALUES,
tp_el_dof.GetData(), restr);
}
void InitStridedRestriction(const mfem::FiniteElementSpace &fes,
CeedInt nelem, CeedInt nqpts, CeedInt qdatasize,
const CeedInt *strides,
@@ -235,57 +139,14 @@ void InitRestriction(const FiniteElementSpace &fes,
// Init or retreive key values
if (restr_itr == mfem::internal::ceed_restr_map.end())
{
InitRestrictionImpl(fes, ceed, restr);
mfem::internal::ceed_restr_map[restr_key] = *restr;
}
else
{
*restr = restr_itr->second;
}
}
void InitRestrictionWithIndices(const FiniteElementSpace &fes,
int nelem,
const int* indices,
Ceed ceed,
CeedElemRestriction *restr)
{
// Check for FES -> basis, restriction in hash tables
const mfem::FiniteElement *fe = fes.GetFE(indices[0]);
const int P = fe->GetDof();
const int ncomp = fes.GetVDim();
RestrKey restr_key(&fes, nelem, P, ncomp, restr_type::Standard);
auto restr_itr = mfem::internal::ceed_restr_map.find(restr_key);
// Init or retreive key values
if (restr_itr == mfem::internal::ceed_restr_map.end())
{
InitRestrictionWithIndicesImpl(fes, nelem, indices, ceed, restr);
mfem::internal::ceed_restr_map[restr_key] = *restr;
}
else
{
*restr = restr_itr->second;
}
}
void InitCoeffRestrictionWithIndices(const FiniteElementSpace &fes,
int nelem,
const int* indices,
int nquads,
int ncomp,
Ceed ceed,
CeedElemRestriction *restr)
{
// Check for FES -> basis, restriction in hash tables
RestrKey restr_key(&fes, nelem, nquads, ncomp, restr_type::Coeff);
auto restr_itr = mfem::internal::ceed_restr_map.find(restr_key);
// Init or retreive key values
if (restr_itr == mfem::internal::ceed_restr_map.end())
{
InitCoeffRestrictionWithIndicesImpl(fes, nelem, indices, nquads, ncomp,
ceed, restr);
if (UsesTensorBasis(fes))
{
InitTensorRestriction(fes, ceed, restr);
}
else
{
InitNonTensorRestriction(fes, ceed, restr);
}
mfem::internal::ceed_restr_map[restr_key] = *restr;
}
else
+23 -49
View File
@@ -21,63 +21,37 @@ namespace ceed
{
#ifdef MFEM_USE_CEED
/** @brief Initialize a CeedElemRestriction for non-mixed meshes.
@param[in] fes Input finite element space.
@param[in] ceed Input Ceed object.
@param[out] restr The address of the initialized CeedElemRestriction object.
*/
void InitRestriction(const FiniteElementSpace &fes,
Ceed ceed,
CeedElemRestriction *restr);
/** @brief Initialize a CeedElemRestriction for mixed meshes.
@param[in] fes The finite element space.
@param[in] ceed The Ceed object.
@param[in] nelem The number of elements.
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`.
@param[out] restr The `CeedElemRestriction` to initialize. */
void InitRestrictionWithIndices(const FiniteElementSpace &fes,
int nelem,
const int* indices,
Ceed ceed,
CeedElemRestriction *restr);
/** @brief Initialize a strided CeedElemRestriction
@param[in] nelem is the number of elements.
@param[in] nqpts is the total number of quadrature points.
@param[in] qdatasize is the number of data per quadrature point.
@param[in] strides Array for strides between [nodes, components, elements].
/// @brief Initialize a strided CeedElemRestriction
/** @a nelem is the number of elements,
@a nqpts is the total number of quadrature points
@a qdatasize is the number of data per quadrature point
@a strides Array for strides between [nodes, components, elements].
Data for node i, component j, element k can be found in the L-vector at
index i*strides[0] + j*strides[1] + k*strides[2]. CEED_STRIDES_BACKEND may
be used with vectors created by a Ceed backend.
@param[out] restr The `CeedElemRestriction` to initialize. */
be used with vectors created by a Ceed backend. */
void InitStridedRestriction(const mfem::FiniteElementSpace &fes,
CeedInt nelem, CeedInt nqpts, CeedInt qdatasize,
const CeedInt *strides,
CeedElemRestriction *restr);
/** @brief Initialize a CeedElemRestriction for a mfem::Coefficient on a mixed
mesh.
/** @brief Initialize a CeedElemRestriction.
*
* @param[in] fes Input finite element space.
* @param[in] ceed Input Ceed object.
@param[out] restr The address of the initialized CeedElemRestriction object.
*/
void InitRestriction(const FiniteElementSpace &fes,
Ceed ceed,
CeedElemRestriction *restr);
@param[in] fes The finite element space.
@param[in] nelem is the number of elements.
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`.
@param[in] nquads is the total number of quadrature points
@param[in] ncomp is the number of data per quadrature point
@param[in] ceed The Ceed object.
@param[out] restr The `CeedElemRestriction` to initialize. */
void InitCoeffRestrictionWithIndices(const FiniteElementSpace &fes,
int nelem,
const int* indices,
int nquads,
int ncomp,
Ceed ceed,
CeedElemRestriction *restr);
/** @brief Initialize a CeedElemRestriction.
*
* @param[in] fes Input finite element space.
* @param[in] ceed Input Ceed object.
@param[out] restr The address of the initialized CeedElemRestriction object.
*/
void InitTensorRestriction(const FiniteElementSpace &fes,
Ceed ceed, CeedElemRestriction *restr);
#endif
-88
View File
@@ -99,34 +99,6 @@ void InitBasisAndRestriction(const FiniteElementSpace &fes,
InitRestriction(fes, ceed, restr);
}
void InitBasisAndRestrictionWithIndices(const FiniteElementSpace &fes,
const IntegrationRule &irm,
int nelem,
const int* indices,
Ceed ceed, CeedBasis *basis,
CeedElemRestriction *restr)
{
InitBasisWithIndices(fes, irm, nelem, indices, ceed, basis);
InitRestrictionWithIndices(fes, nelem, indices, ceed, restr);
}
void InitBasisAndRestriction(const FiniteElementSpace &fes,
const IntegrationRule &irm,
int nelem,
const int* indices,
Ceed ceed, CeedBasis *basis,
CeedElemRestriction *restr)
{
if (indices)
{
InitBasisAndRestrictionWithIndices(fes,irm,nelem,indices,ceed,basis,restr);
}
else
{
InitBasisAndRestriction(fes,irm,ceed,basis,restr);
}
}
// Assumes a tensor-product operator with one active field
int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
{
@@ -186,66 +158,6 @@ int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
return 0;
}
template <>
const IntegrationRule & GetRule<MassIntegrator>(
const MassIntegrator &integ,
const FiniteElement &trial_fe,
const FiniteElement &test_fe,
ElementTransformation &trans)
{
return MassIntegrator::GetRule(trial_fe, test_fe, trans);
}
template <>
const IntegrationRule & GetRule<VectorMassIntegrator>(
const VectorMassIntegrator &integ,
const FiniteElement &trial_fe,
const FiniteElement &test_fe,
ElementTransformation &trans)
{
return MassIntegrator::GetRule(trial_fe, test_fe, trans);
}
template <>
const IntegrationRule & GetRule<ConvectionIntegrator>(
const ConvectionIntegrator &integ,
const FiniteElement &trial_fe,
const FiniteElement &test_fe,
ElementTransformation &trans)
{
return ConvectionIntegrator::GetRule(trial_fe, test_fe, trans);
}
template <>
const IntegrationRule & GetRule<VectorConvectionNLFIntegrator>(
const VectorConvectionNLFIntegrator &integ,
const FiniteElement &trial_fe,
const FiniteElement &test_fe,
ElementTransformation &trans)
{
return VectorConvectionNLFIntegrator::GetRule(trial_fe, trans);
}
template <>
const IntegrationRule & GetRule<DiffusionIntegrator>(
const DiffusionIntegrator &integ,
const FiniteElement &trial_fe,
const FiniteElement &test_fe,
ElementTransformation &trans)
{
return DiffusionIntegrator::GetRule(trial_fe, test_fe);
}
template <>
const IntegrationRule & GetRule<VectorDiffusionIntegrator>(
const VectorDiffusionIntegrator &integ,
const FiniteElement &trial_fe,
const FiniteElement &test_fe,
ElementTransformation &trans)
{
return DiffusionIntegrator::GetRule(trial_fe, test_fe);
}
std::string ceed_path;
const std::string &GetCeedPath()
+3 -43
View File
@@ -26,9 +26,7 @@
namespace mfem
{
class FiniteElement;
class FiniteElementSpace;
class ElementTransformation;
class IntegrationRule;
class Vector;
@@ -57,51 +55,15 @@ void RemoveBasisAndRestriction(const mfem::FiniteElementSpace *fes);
/// Initialize a CeedVector from an mfem::Vector
void InitVector(const mfem::Vector &v, CeedVector &cv);
/** @brief Initialize a CeedBasis and a CeedElemRestriction based on an
mfem::FiniteElementSpace @a fes, and an mfem::IntegrationRule @a ir.
@param[in] fes The finite element space.
@param[in] ir The integration rule.
@param[in] ceed The Ceed object.
@param[out] basis The `CeedBasis` to initialize.
@param[out] restr The `CeedElemRestriction` to initialize.
@warning Only for non-mixed finite element spaces. */
/** Initialize a CeedBasis and a CeedElemRestriction based on an
mfem::FiniteElementSpace @a fes, and an mfem::IntegrationRule @a ir. */
void InitBasisAndRestriction(const mfem::FiniteElementSpace &fes,
const mfem::IntegrationRule &ir,
Ceed ceed, CeedBasis *basis,
CeedElemRestriction *restr);
/** @brief Initialize a CeedBasis and a CeedElemRestriction based on an
mfem::FiniteElementSpace @a fes, and an mfem::IntegrationRule @a ir,
and a list of @a nelem elements of indices @a indices.
@param[in] fes The finite element space.
@param[in] ir The integration rule.
@param[in] nelem The number of elements.
@param[in] indices The indices of the elements of same type in the
`FiniteElementSpace`. If `indices == nullptr`, assumes
that the `FiniteElementSpace` is not mixed.
@param[in] ceed The Ceed object.
@param[out] basis The `CeedBasis` to initialize.
@param[out] restr The `CeedElemRestriction` to initialize. */
void InitBasisAndRestriction(const FiniteElementSpace &fes,
const IntegrationRule &ir,
int nelem,
const int* indices,
Ceed ceed, CeedBasis *basis,
CeedElemRestriction *restr);
int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field);
template <typename Integrator>
const IntegrationRule & GetRule(
const Integrator &integ,
const FiniteElement &trial_fe,
const FiniteElement &test_fe,
ElementTransformation &Trans);
/// Return the path to the libCEED q-function headers.
const std::string &GetCeedPath();
@@ -125,7 +87,7 @@ struct BasisHash
};
using BasisMap = std::unordered_map<const BasisKey, CeedBasis, BasisHash>;
enum restr_type {Standard, Strided, Coeff};
enum restr_type {Standard, Strided};
// Hash table for CeedElemRestriction
using RestrKey =
@@ -155,8 +117,6 @@ namespace internal
{
#ifdef MFEM_USE_CEED
/** @warning These maps have a tendency to create bugs when adding new "types"
of CeedBasis and CeedElemRestriction. */
extern ceed::BasisMap ceed_basis_map;
extern ceed::RestrMap ceed_restr_map;
#endif
+1 -1
View File
@@ -633,7 +633,7 @@ AlgebraicSpaceHierarchy::AlgebraicSpaceHierarchy(FiniteElementSpace &fes)
current_order = order;
Ceed ceed = internal::ceed;
InitRestriction(fes, ceed, &fine_er);
InitTensorRestriction(fes, ceed, &fine_er);
CeedElemRestriction er = fine_er;
int dim = fes.GetMesh()->Dimension();
+1
View File
@@ -134,6 +134,7 @@ public:
}
~AlgebraicSpaceHierarchy()
{
CeedElemRestrictionDestroy(&fine_er);
for (int i=0; i<R_tr.Size(); ++i)
{
delete R_tr[i];
+2 -2
View File
@@ -1078,7 +1078,7 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
it->second->GetValues(i, RefG->RefPts, val, pmat);
for (int j = 0; j < val.Size(); j++)
{
WriteBinaryOrASCII(os, buf, val(j), "\n", pv_data_format);
WriteBinaryOrASCII(out, buf, val(j), "\n", pv_data_format);
}
}
}
@@ -1094,7 +1094,7 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
{
for (int ii = 0; ii < vval.Height(); ii++)
{
WriteBinaryOrASCII(os, buf, vval(ii,jj), " ", pv_data_format);
WriteBinaryOrASCII(out, buf, vval(ii,jj), " ", pv_data_format);
}
if (pv_data_format == VTKFormat::ASCII) { os << '\n'; }
}
+5 -5
View File
@@ -980,14 +980,14 @@ public:
void Save(std::ostream &out) const;
};
/// @brief Return true if the mesh contains only one topology and the elements are tensor elements.
inline bool UsesTensorBasis(const FiniteElementSpace& fes)
{
Mesh & mesh = *fes.GetMesh();
const bool mixed = mesh.GetNumGeometries(mesh.Dimension()) > 1;
// TODO: mixed meshes: return true if there is at least one tensor-product
// Geometry in the global mesh and the FE collection returns a
// TensorBasisElement for that Geometry?
// Potential issue: empty local mesh --> no element 0.
return !mixed &&
dynamic_cast<const mfem::TensorBasisElement *>(fes.GetFE(0))!=nullptr;
return dynamic_cast<const mfem::TensorBasisElement *>(fes.GetFE(0))!=nullptr;
}
}
+2 -2
View File
@@ -4276,7 +4276,7 @@ void TensorProductLegendre(int dim, // input
poly1d.CalcLegendre(order, x3, poly_z);
}
int basis_dimension = static_cast<int>(pow(order+1,dim));
int basis_dimension = pow(order+1,dim);
poly.SetSize(basis_dimension);
switch (dim)
{
@@ -4458,7 +4458,7 @@ double LSZZErrorEstimator(BilinearFormIntegrator &blfi, // input
const int patch_order = max(ufes->GetElementOrder(el1),
ufes->GetElementOrder(el2));
int num_basis_functions = static_cast<int>(pow(patch_order+1,dim));
int num_basis_functions = pow(patch_order+1,dim);
int flux_order = 2*patch_order + 1;
DenseMatrix A(num_basis_functions);
Array<double> b(sdim * num_basis_functions);
+127 -2
View File
@@ -26,15 +26,134 @@ LinearForm::LinearForm(FiniteElementSpace *f, LinearForm *lf)
ext = nullptr;
extern_lfs = 1;
// Copy the pointers to the integrators
// Copy the pointers to the integrators and markers
domain_integs = lf->domain_integs;
domain_integs_marker = lf->domain_integs_marker;
domain_delta_integs = lf->domain_delta_integs;
boundary_integs = lf->boundary_integs;
boundary_integs_marker = lf->boundary_integs_marker;
boundary_face_integs = lf->boundary_face_integs;
boundary_face_integs_marker = lf->boundary_face_integs_marker;
interior_face_integs = lf->interior_face_integs;
}
LinearForm::LinearForm(LinearForm &&other)
: Vector(std::move(other)), fes(other.fes), ext(other.ext),
extern_lfs(other.extern_lfs)
{
// Linear forms are stored on the device
UseDevice(true);
// We swap stored integrators and markers with the moved LinearForm
mfem::Swap(domain_integs, other.domain_integs);
mfem::Swap(domain_integs_marker, other.domain_integs_marker);
mfem::Swap(domain_delta_integs, other.domain_delta_integs);
mfem::Swap(boundary_integs, other.boundary_integs);
mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
mfem::Swap(boundary_face_integs, other.boundary_face_integs);
mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
mfem::Swap(interior_face_integs, other.interior_face_integs);
other.fes = nullptr;
other.ext = nullptr;
// moved LinearForm now set to owns its integrators, though none are stored
// since we swapped our newly constructed empty Arrays with its Arrays
other.extern_lfs = 0;
}
LinearForm& LinearForm::operator=(LinearForm &&other)
{
if (this != &other)
{
Vector::operator=(std::move(other));
/// If we own our integrators, delete them
if (extern_lfs == 0)
{
for (int k = 0; k < domain_delta_integs.Size(); k++)
{ delete domain_delta_integs[k]; }
for (int k = 0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
for (int k = 0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
for (int k = 0; k < boundary_face_integs.Size(); k++)
{ delete boundary_face_integs[k]; }
for (int k = 0; k < interior_face_integs.Size(); k++)
{ delete interior_face_integs[k]; }
}
/// Null out all integs
for (int k = 0; k < domain_delta_integs.Size(); k++)
{
domain_delta_integs[k] = nullptr;
}
domain_delta_integs.SetSize(0);
for (int k = 0; k < domain_integs.Size(); k++)
{
domain_integs[k] = nullptr;
}
domain_integs.SetSize(0);
for (int k = 0; k < boundary_integs.Size(); k++)
{
boundary_integs[k] = nullptr;
}
boundary_integs.SetSize(0);
for (int k = 0; k < boundary_face_integs.Size(); k++)
{
boundary_face_integs[k] = nullptr;
}
boundary_face_integs.SetSize(0);
for (int k = 0; k < interior_face_integs.Size(); k++)
{
interior_face_integs[k] = nullptr;
}
interior_face_integs.SetSize(0);
/// Null out all markers
for (int k = 0; k < domain_integs_marker.Size(); ++k)
{
domain_integs_marker[k] = nullptr;
}
domain_integs_marker.SetSize(0);
for (int k = 0; k < boundary_integs_marker.Size(); ++k)
{
boundary_integs_marker[k] = nullptr;
}
boundary_integs_marker.SetSize(0);
for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
{
boundary_face_integs_marker[k] = nullptr;
}
boundary_face_integs_marker.SetSize(0);
Vector::operator=(std::move(other));
// swap stored integrators and markers with the moved LinearForm
mfem::Swap(domain_integs, other.domain_integs);
mfem::Swap(domain_integs_marker, other.domain_integs_marker);
mfem::Swap(domain_delta_integs, other.domain_delta_integs);
mfem::Swap(boundary_integs, other.boundary_integs);
mfem::Swap(boundary_integs_marker, other.boundary_integs_marker);
mfem::Swap(boundary_face_integs, other.boundary_face_integs);
mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
mfem::Swap(interior_face_integs, other.interior_face_integs);
fes = other.fes;
other.fes = nullptr;
ext = other.ext;
other.ext = nullptr;
extern_lfs = other.extern_lfs;
// moved LinearForm now set to owns its integrators, though none are stored
// since we swapped our empty Arrays with its Arrays
other.extern_lfs = 0;
}
return *this;
}
void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi)
@@ -129,7 +248,13 @@ bool LinearForm::SupportsDevice()
if (mesh_dim == 1 || mesh_dim != mesh.SpaceDimension()) { return false; }
// tensor-product finite element space only
if (!UsesTensorBasis(*fes)) { return false; }
// with point values preserving scalar fields
for (int e = 0; e < fes->GetNE(); ++e)
{
const FiniteElement *fe = fes->GetFE(e);
if (fe->GetMapType() != FiniteElement::VALUE) { return false; }
if (!dynamic_cast<const TensorBasisElement*>(fe)) { return false; }
}
return true;
}
+13 -4
View File
@@ -74,10 +74,6 @@ protected:
/// Force (re)computation of delta locations.
void ResetDeltaLocations() { domain_delta_integs_elem_id.SetSize(0); }
private:
/// Copy construction is not supported; body is undefined.
LinearForm(const LinearForm &);
public:
/// Creates linear form associated with FE space @a *f.
/** The pointer @a f is not owned by the newly constructed object. */
@@ -107,6 +103,9 @@ public:
LinearForm(FiniteElementSpace *f, double *data) : Vector(data, f->GetVSize())
{ fes = f; ext = nullptr; extern_lfs = 0; }
/// Explicitly prohibit copy construction of LinearForm
LinearForm(const LinearForm &other) = delete;
/// Copy assignment. Only the data of the base class Vector is copied.
/** It is assumed that this object and @a rhs use FiniteElementSpace%s that
have the same size.
@@ -116,6 +115,16 @@ public:
LinearForm &operator=(const LinearForm &rhs)
{ return operator=((const Vector &)rhs); }
/// Move constructor for LinearForm.
/** This constructor "steals" the owned data members from the @a other
LinearForm. */
LinearForm(LinearForm &&other);
/// Move assignment operator for LinearForm
/** This assignment first frees all owned data, then "steals" the owned data
members from the @a other LinearForm. */
LinearForm& operator=(LinearForm &&other);
/// (DEPRECATED) Return the FE space associated with the LinearForm.
/** @deprecated Use FESpace() instead. */
MFEM_DEPRECATED FiniteElementSpace *GetFES() { return fes; }
+22 -41
View File
@@ -18,9 +18,8 @@ namespace mfem
template<int T_D1D = 0, int T_Q1D = 0> static
void DLFEvalAssemble2D(const int vdim, const int ne, const int d, const int q,
const int map_type, const int *markers, const double *b,
const double *j, const double *weights,
const Vector &coeff, double *y)
const int *markers, const double *b, const double *j,
const double *weights, const Vector &coeff, double *y)
{
const auto F = coeff.Read();
const auto M = Reshape(markers, ne);
@@ -55,19 +54,11 @@ void DLFEvalAssemble2D(const int vdim, const int ne, const int d, const int q,
{
MFEM_FOREACH_THREAD(y,y,q)
{
double detJ;
if (map_type == FiniteElement::VALUE)
{
const double J11 = J(x,y,0,0,e);
const double J21 = J(x,y,1,0,e);
const double J12 = J(x,y,0,1,e);
const double J22 = J(x,y,1,1,e);
detJ = J11 * J22 - J21 * J12;
}
else
{
detJ = 1.0;
}
const double J11 = J(x,y,0,0,e);
const double J21 = J(x,y,1,0,e);
const double J12 = J(x,y,0,1,e);
const double J22 = J(x,y,1,1,e);
const double detJ = J11 * J22 - J21 * J12;
const double coeff_val = cst ? cst_val : C(c,x,y,e);
QQ(y,x) = W(x,y) * coeff_val * detJ;
}
@@ -99,9 +90,8 @@ void DLFEvalAssemble2D(const int vdim, const int ne, const int d, const int q,
template<int T_D1D = 0, int T_Q1D = 0> static
void DLFEvalAssemble3D(const int vdim, const int ne, const int d, const int q,
const int map_type, const int *markers, const double *b,
const double *j, const double *weights,
const Vector &coeff, double *y)
const int *markers, const double *b, const double *j,
const double *weights, const Vector &coeff, double *y)
{
const auto F = coeff.Read();
const auto M = Reshape(markers, ne);
@@ -138,26 +128,18 @@ void DLFEvalAssemble3D(const int vdim, const int ne, const int d, const int q,
{
for (int z = 0; z < q; ++z)
{
double detJ;
if (map_type == FiniteElement::VALUE)
{
const double J11 = J(x,y,z,0,0,e);
const double J21 = J(x,y,z,1,0,e);
const double J31 = J(x,y,z,2,0,e);
const double J12 = J(x,y,z,0,1,e);
const double J22 = J(x,y,z,1,1,e);
const double J32 = J(x,y,z,2,1,e);
const double J13 = J(x,y,z,0,2,e);
const double J23 = J(x,y,z,1,2,e);
const double J33 = J(x,y,z,2,2,e);
detJ = J11 * (J22 * J33 - J32 * J23) -
/* */ J21 * (J12 * J33 - J32 * J13) +
/* */ J31 * (J12 * J23 - J22 * J13);
}
else
{
detJ = 1.0;
}
const double J11 = J(x,y,z,0,0,e);
const double J21 = J(x,y,z,1,0,e);
const double J31 = J(x,y,z,2,0,e);
const double J12 = J(x,y,z,0,1,e);
const double J22 = J(x,y,z,1,1,e);
const double J32 = J(x,y,z,2,1,e);
const double J13 = J(x,y,z,0,2,e);
const double J23 = J(x,y,z,1,2,e);
const double J33 = J(x,y,z,2,2,e);
const double detJ = J11 * (J22 * J33 - J32 * J23) -
/* */ J21 * (J12 * J33 - J32 * J13) +
/* */ J31 * (J12 * J23 - J22 * J13);
const double coeff_val = cst_coeff ? cst_val : C(c,x,y,z,e);
QQQ(z,y,x) = W(x,y,z) * coeff_val * detJ;
}
@@ -224,7 +206,6 @@ static void DLFEvalAssemble(const FiniteElementSpace &fes,
const int d = maps.ndof, q = maps.nqpt;
constexpr int flags = GeometricFactors::JACOBIANS;
const GeometricFactors *geom = mesh->GetGeometricFactors(*ir, flags, mt);
const int map_type = fes.GetFE(0)->GetMapType();
decltype(&DLFEvalAssemble2D<>) ker =
dim == 2 ? DLFEvalAssemble2D<> : DLFEvalAssemble3D<>;
@@ -263,7 +244,7 @@ static void DLFEvalAssemble(const FiniteElementSpace &fes,
const double *J = geom->J.Read();
const double *W = ir->GetWeights().Read();
double *Y = y.ReadWrite();
ker(vdim, ne, d, q, map_type, M, B, J, W, coeff, Y);
ker(vdim, ne, d, q, M, B, J, W, coeff, Y);
}
void DomainLFIntegrator::AssembleDevice(const FiniteElementSpace &fes,
+2 -4
View File
@@ -467,8 +467,7 @@ void LORDiscretization::FormLORSpace()
mesh = new Mesh(Mesh::MakeRefined(mesh_ho, refinements, ref_type));
fec = fes_ho.FEColl()->Clone(GetLOROrder());
const int vdim = fes_ho.GetVDim();
fes = new FiniteElementSpace(mesh, fec, vdim);
fes = new FiniteElementSpace(mesh, fec);
SetupProlongationAndRestriction();
}
@@ -512,8 +511,7 @@ void ParLORDiscretization::FormLORSpace()
mesh = pmesh;
fec = pfes_ho.FEColl()->Clone(GetLOROrder());
const int vdim = fes_ho.GetVDim();
ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec, vdim);
ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec);
fes = pfes;
SetupProlongationAndRestriction();
}
+6 -6
View File
@@ -39,7 +39,7 @@ void BatchedLOR_AMS::Form2DEdgeToVertex_ND(Array<int> &edge2vert)
{
const int o = order;
const int op1 = o + 1;
const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
const int nedge = dim*o*pow(op1, dim-1);
edge2vert.SetSize(2*nedge);
auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
@@ -73,7 +73,7 @@ void BatchedLOR_AMS::Form2DEdgeToVertex_RT(Array<int> &edge2vert)
{
const int o = order;
const int op1 = o + 1;
const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
const int nedge = dim*o*pow(op1, dim-1);
edge2vert.SetSize(2*nedge);
auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
@@ -106,7 +106,7 @@ void BatchedLOR_AMS::Form3DEdgeToVertex(Array<int> &edge2vert)
{
const int o = order;
const int op1 = o + 1;
const int nedge = static_cast<int>(dim*o*pow(op1, dim-1));
const int nedge = dim*o*pow(op1, dim-1);
edge2vert.SetSize(2*nedge);
auto e2v = Reshape(edge2vert.HostWrite(), 2, nedge);
@@ -170,8 +170,8 @@ void BatchedLOR_AMS::FormGradientMatrix()
MFEM_VERIFY(R_v != NULL && R_e != NULL, "");
const int nel_ho = edge_fes.GetNE();
const int nedge_per_el = static_cast<int>(dim*order*pow(order + 1, dim - 1));
const int nvert_per_el = static_cast<int>(pow(order + 1, dim));
const int nedge_per_el = dim*order*pow(order + 1, dim - 1);
const int nvert_per_el = pow(order + 1, dim);
const auto offsets_e = R_e->Offsets().Read();
const auto indices_e = R_e->Indices().Read();
@@ -274,7 +274,7 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
const int nel_ho = vert_fes.GetNE();
const int ndp1 = order + 1;
const int ndof_per_el = static_cast<int>(pow(ndp1, dim));
const int ndof_per_el = pow(ndp1, dim);
const int sdim = dim;
const int ntdofs = R->Height();
+246 -11
View File
@@ -13,7 +13,6 @@
#include "../../fem/quadinterpolator.hpp"
#include "../../general/forall.hpp"
#include <climits>
#include "../pbilinearform.hpp"
// Specializations
#include "lor_h1.hpp"
@@ -46,6 +45,30 @@ bool HasIntegrators(BilinearForm &a)
return false;
}
#ifdef MFEM_USE_MPI
void HypreStealOwnership(HypreParMatrix &A_hyp, SparseMatrix &A_diag)
{
#ifndef HYPRE_BIGINT
bool own_i = A_hyp.GetDiagMemoryI().OwnsHostPtr();
bool own_j = A_hyp.GetDiagMemoryJ().OwnsHostPtr();
MFEM_CONTRACT_VAR(own_j);
MFEM_ASSERT(own_i == own_j, "Inconsistent ownership");
if (!own_i)
{
std::swap(A_diag.GetMemoryI(), A_hyp.GetDiagMemoryI());
std::swap(A_diag.GetMemoryJ(), A_hyp.GetDiagMemoryJ());
}
#endif
if (!A_hyp.GetDiagMemoryData().OwnsHostPtr())
{
std::swap(A_diag.GetMemoryData(), A_hyp.GetDiagMemoryData());
}
A_hyp.SetOwnerFlags(3, A_hyp.OwnsOffd(), A_hyp.OwnsColMap());
}
#endif
bool BatchedLORAssembly::FormIsSupported(BilinearForm &a)
{
const FiniteElementCollection *fec = a.FESpace()->FEColl();
@@ -80,7 +103,7 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
const int nel_ho = mesh_ho.GetNE();
const int order = fes_ho.GetMaxElementOrder();
const int nd1d = order + 1;
const int ndof_per_el = static_cast<int>(pow(nd1d, dim));
const int ndof_per_el = pow(nd1d, dim);
const GridFunction *nodal_gf = mesh_ho.GetNodes();
const FiniteElementSpace *nodal_fes = nodal_gf->FESpace();
@@ -284,7 +307,7 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
const bool plus = si_E >= 0;
const int i_E = plus ? si_E : -1 - si_E;
i_elts[e_i] = i_E/ndof_per_el;
const int i_Bi = i_E % ndof_per_el;
const double i_Bi = i_E%ndof_per_el;
i_B[e_i] = plus ? i_Bi : -1 - i_Bi; // encode with sign
}
for (int j=0; j<nnz_per_row; ++j)
@@ -314,7 +337,7 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
const bool plus = sj_E >= 0;
const int j_E = plus ? sj_E : -1 - sj_E;
j_elts[e_j] = j_E/ndof_per_el;
const int j_Bj = j_E % ndof_per_el;
const double j_Bj = j_E%ndof_per_el;
j_B[e_j] = plus ? j_Bj : -1 - j_Bj; // encode with sign
}
const int min_e = GetMinElt(i_elts, i_ne, j_elts, j_ne);
@@ -460,13 +483,199 @@ void BatchedLORAssembly::ParAssemble(
OperatorHandle A_local;
AssembleWithoutBC(a, A_local);
ParBilinearForm *pa =
dynamic_cast<ParBilinearForm*>(&a);
ParFiniteElementSpace *pfes_ho =
dynamic_cast<ParFiniteElementSpace*>(&fes_ho);
MFEM_VERIFY(pfes_ho != nullptr,
"ParAssemble must be called with ParFiniteElementSpace");
pa->ParallelRAP(*A_local.As<SparseMatrix>(), A, true);
// Create a block diagonal parallel matrix
OperatorHandle A_diag(Operator::Hypre_ParCSR);
A_diag.MakeSquareBlockDiag(pfes_ho->GetComm(),
pfes_ho->GlobalVSize(),
pfes_ho->GetDofOffsets(),
A_local.As<SparseMatrix>());
A.As<HypreParMatrix>()->EliminateBC(ess_dofs,
Operator::DiagonalPolicy::DIAG_ONE);
// Parallel matrix assembly using P^t A P (if needed)
if (IsIdentityProlongation(pfes_ho->GetProlongationMatrix()))
{
A_diag.SetOperatorOwner(false);
A.Reset(A_diag.Ptr());
HypreStealOwnership(*A.As<HypreParMatrix>(), *A_local.As<SparseMatrix>());
}
else
{
OperatorHandle P(Operator::Hypre_ParCSR);
P.ConvertFrom(pfes_ho->Dof_TrueDof_Matrix());
A.MakePtAP(A_diag, P);
}
// Eliminate the boundary conditions
HypreParMatrix *A_mat = A.As<HypreParMatrix>();
hypre_ParCSRMatrix *A_hypre = *A_mat;
A_mat->HypreReadWrite();
hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A_hypre);
hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A_hypre);
HYPRE_Int diag_nrows = hypre_CSRMatrixNumRows(diag);
HYPRE_Int offd_ncols = hypre_CSRMatrixNumCols(offd);
const int n_ess_dofs = ess_dofs.Size();
const auto ess_dofs_d = ess_dofs.GetMemory().Read(
GetHypreMemoryClass(), n_ess_dofs);
// Start communication to figure out which columns need to be eliminated in
// the off-diagonal block
hypre_ParCSRCommHandle *comm_handle;
HYPRE_Int *int_buf_data, *eliminate_row, *eliminate_col;
{
eliminate_row = mfem_hypre_CTAlloc_host(HYPRE_Int, diag_nrows);
eliminate_col = mfem_hypre_CTAlloc_host(HYPRE_Int, offd_ncols);
// Make sure A has a communication package
hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A_hypre);
if (!comm_pkg)
{
hypre_MatvecCommPkgCreate(A_hypre);
comm_pkg = hypre_ParCSRMatrixCommPkg(A_hypre);
}
// Which of the local rows are to be eliminated?
for (int i = 0; i < diag_nrows; i++)
{
eliminate_row[i] = 0;
}
ess_dofs.HostRead();
for (int i = 0; i < n_ess_dofs; i++)
{
eliminate_row[ess_dofs[i]] = 1;
}
// Use a matvec communication pattern to find (in eliminate_col) which of
// the local offd columns are to be eliminated
HYPRE_Int num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
int_buf_data = mfem_hypre_CTAlloc_host(
HYPRE_Int,
hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
int index = 0;
for (int i = 0; i < num_sends; i++)
{
int start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
for (int j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
{
int k = hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j);
int_buf_data[index++] = eliminate_row[k];
}
}
comm_handle = hypre_ParCSRCommHandleCreate(
11, comm_pkg, int_buf_data, eliminate_col);
}
// Eliminate rows and columns in the diagonal block
{
const auto I = diag->i;
const auto J = diag->j;
auto data = diag->data;
MFEM_HYPRE_FORALL(i, n_ess_dofs,
{
const int idof = ess_dofs_d[i];
for (int j=I[idof]; j<I[idof+1]; ++j)
{
const int jdof = J[j];
if (jdof == idof)
{
// Set eliminate diagonal equal to identity
data[j] = 1.0;
}
else
{
data[j] = 0.0;
for (int k=I[jdof]; k<I[jdof+1]; ++k)
{
if (J[k] == idof)
{
data[k] = 0.0;
break;
}
}
}
}
});
}
// Eliminate rows in the off-diagonal block
{
const auto I = offd->i;
auto data = offd->data;
MFEM_HYPRE_FORALL(i, n_ess_dofs,
{
const int idof = ess_dofs_d[i];
for (int j=I[idof]; j<I[idof+1]; ++j)
{
data[j] = 0.0;
}
});
}
// Wait for MPI communication to finish
Array<HYPRE_Int> cols_to_eliminate;
{
hypre_ParCSRCommHandleDestroy(comm_handle);
// set the array cols_to_eliminate
int ncols_to_eliminate = 0;
for (int i = 0; i < offd_ncols; i++)
{
if (eliminate_col[i]) { ncols_to_eliminate++; }
}
cols_to_eliminate.SetSize(ncols_to_eliminate);
cols_to_eliminate = 0.0;
ncols_to_eliminate = 0;
for (int i = 0; i < offd_ncols; i++)
{
if (eliminate_col[i])
{
cols_to_eliminate[ncols_to_eliminate++] = i;
}
}
mfem_hypre_TFree_host(int_buf_data);
mfem_hypre_TFree_host(eliminate_row);
mfem_hypre_TFree_host(eliminate_col);
}
// Eliminate columns in the off-diagonal block
{
const int ncols_to_eliminate = cols_to_eliminate.Size();
const int nrows_offd = hypre_CSRMatrixNumRows(offd);
const auto cols = cols_to_eliminate.GetMemory().Read(
GetHypreMemoryClass(), ncols_to_eliminate);
const auto I = offd->i;
const auto J = offd->j;
auto data = offd->data;
// Note: could also try a different strategy, looping over nnz in the
// matrix and then doing a binary search in ncols_to_eliminate to see if
// the column should be eliminated.
MFEM_HYPRE_FORALL(idx, ncols_to_eliminate,
{
const int j = cols[idx];
for (int i=0; i<nrows_offd; ++i)
{
for (int jj=I[i]; jj<I[i+1]; ++jj)
{
if (J[jj] == j)
{
data[jj] = 0.0;
break;
}
}
}
});
}
}
#endif
@@ -483,8 +692,34 @@ void BatchedLORAssembly::Assemble(
AssembleWithoutBC(a, A);
SparseMatrix *A_mat = A.As<SparseMatrix>();
A_mat->EliminateBC(ess_dofs,
Operator::DiagonalPolicy::DIAG_KEEP);
// Eliminate essential DOFs (BCs) from the matrix (what we do here is
// equivalent to DiagonalPolicy::DIAG_KEEP).
const int n_ess_dofs = ess_dofs.Size();
const auto ess_dofs_d = ess_dofs.Read();
const auto I = A_mat->ReadI();
const auto J = A_mat->ReadJ();
auto dA = A_mat->ReadWriteData();
MFEM_FORALL(i, n_ess_dofs,
{
const int idof = ess_dofs_d[i];
for (int j=I[idof]; j<I[idof+1]; ++j)
{
const int jdof = J[j];
if (jdof != idof)
{
dA[j] = 0.0;
for (int k=I[jdof]; k<I[jdof+1]; ++k)
{
if (J[k] == idof)
{
dA[k] = 0.0;
break;
}
}
}
}
});
}
BatchedLORAssembly::BatchedLORAssembly(FiniteElementSpace &fes_ho_)
+16
View File
@@ -143,6 +143,22 @@ static T *GetIntegrator(BilinearForm &a)
return nullptr;
}
#ifdef MFEM_USE_MPI
/// @brief Make @a A_hyp steal ownership of its diagonal part @a A_diag.
///
/// If @a A_hyp does not own I and J, then they are aliases pointing to the I
/// and J arrays in @a A_diag. In that case, this function swaps the memory
/// objects. Similarly for the data array.
///
/// After this function is called, @a A_hyp will own all of the arrays of its
/// diagonal part.
///
/// @note I and J can only be aliases when HYPRE_BIGINT is disabled.
void HypreStealOwnership(HypreParMatrix &A_hyp, SparseMatrix &A_diag);
#endif
/// Abstract base class for the batched LOR assembly kernels.
class BatchedLORKernel
{
+254 -85
View File
@@ -15,6 +15,92 @@
namespace mfem
{
NonlinearForm::NonlinearForm(NonlinearForm &&other)
: Operator(other.fes->GetTrueVSize()), assembly(other.assembly),
ext(other.ext), fes(other.fes), Grad(other.Grad), cGrad(other.cGrad),
sequence(other.fes->GetSequence()), P(other.fes->GetProlongationMatrix()),
cP(dynamic_cast<const SparseMatrix*>(P))
{
// We swap stored integrators and markers with the moved nonlinear form
mfem::Swap(domain_integs, other.domain_integs);
mfem::Swap(domain_integs_marker, other.domain_integs_marker);
mfem::Swap(interior_face_integs, other.interior_face_integs);
mfem::Swap(boundary_face_integs, other.boundary_face_integs);
mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
/// Leave the moved nonlinear form in a state as if it was just constructed
/// with fes
other.ext = nullptr;
other.cGrad = nullptr;
other.Grad = nullptr;
other.assembly = AssemblyLevel::LEGACY;
}
NonlinearForm& NonlinearForm::operator=(NonlinearForm &&other)
{
if (this != &other)
{
/// Cleanup current nonlinear form first
delete cGrad;
delete Grad;
for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
delete ext;
/// Null out all our integs and set size of their arrays to zero
for (int k = 0; k < domain_integs.Size(); k++)
{
domain_integs[k] = nullptr;
}
domain_integs.SetSize(0);
for (int k = 0; k < boundary_face_integs.Size(); k++)
{
boundary_face_integs[k] = nullptr;
}
boundary_face_integs.SetSize(0);
for (int k = 0; k < interior_face_integs.Size(); k++)
{
interior_face_integs[k] = nullptr;
}
interior_face_integs.SetSize(0);
/// Null out all our markers and set size of their arrays to zero
for (int k = 0; k < domain_integs_marker.Size(); ++k)
{
domain_integs_marker[k] = nullptr;
}
domain_integs_marker.SetSize(0);
for (int k = 0; k < boundary_face_integs_marker.Size(); ++k)
{
boundary_face_integs_marker[k] = nullptr;
}
boundary_face_integs_marker.SetSize(0);
/// Now steal data from other nonlinear form leaving it in a state as if
/// it was just constructed with fes
Operator::operator=(std::move(other));
assembly = other.assembly;
other.assembly = AssemblyLevel::LEGACY;
Grad = other.Grad;
other.Grad = nullptr;
cGrad = other.cGrad;
other.cGrad = nullptr;
// Swap our empty integ and marker arrays with the moved nonlinear form
mfem::Swap(domain_integs, other.domain_integs);
mfem::Swap(domain_integs_marker, other.domain_integs_marker);
mfem::Swap(interior_face_integs, other.interior_face_integs);
mfem::Swap(boundary_face_integs, other.boundary_face_integs);
mfem::Swap(boundary_face_integs_marker, other.boundary_face_integs_marker);
ext = other.ext;
other.ext = nullptr;
}
return *this;
}
void NonlinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
{
if (ext)
@@ -38,6 +124,38 @@ void NonlinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
}
}
void NonlinearForm::AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
{
domain_integs.Append(nlfi);
domain_integs_marker.Append(nullptr); // null marker means apply everywhere
}
void NonlinearForm::AddDomainIntegrator(NonlinearFormIntegrator *nlfi,
Array<int> &elem_marker)
{
domain_integs.Append(nlfi);
domain_integs_marker.Append(&elem_marker);
}
void NonlinearForm::AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi)
{
interior_face_integs.Append(nlfi);
}
void NonlinearForm::AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi)
{
boundary_face_integs.Append(nlfi);
// null marker means apply everywhere
boundary_face_integs_marker.Append(nullptr);
}
void NonlinearForm::AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi,
Array<int> &bdr_marker)
{
boundary_face_integs.Append(nlfi);
boundary_face_integs_marker.Append(&bdr_marker);
}
void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
Vector *rhs)
{
@@ -87,8 +205,10 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
{
if (ext)
{
MFEM_VERIFY(!fnfi.Size(), "Interior faces terms not yet implemented!");
MFEM_VERIFY(!bfnfi.Size(), "Boundary face terms not yet implemented!");
MFEM_VERIFY(!interior_face_integs.Size(),
"Interior faces terms not yet implemented!");
MFEM_VERIFY(!boundary_face_integs.Size(),
"Boundary face terms not yet implemented!");
return ext->GetGridFunctionEnergy(x);
}
@@ -97,30 +217,47 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
const FiniteElement *fe;
ElementTransformation *T;
DofTransformation *doftrans;
Mesh *mesh = fes->GetMesh();
double energy = 0.0;
if (dnfi.Size())
if (domain_integs.Size())
{
for (int k = 0; k < domain_integs.Size(); k++)
{
if (domain_integs_marker[k] != nullptr)
{
MFEM_VERIFY(mesh->attributes.Size() ==
domain_integs_marker[k]->Size(),
"invalid element marker for domain integrator #"
<< k << ", counting from zero");
}
}
for (int i = 0; i < fes->GetNE(); i++)
{
int elem_attr = mesh->GetAttribute(i);
fe = fes->GetFE(i);
doftrans = fes->GetElementVDofs(i, vdofs);
T = fes->GetElementTransformation(i);
x.GetSubVector(vdofs, el_x);
if (doftrans) {doftrans->InvTransformPrimal(el_x); }
for (int k = 0; k < dnfi.Size(); k++)
for (int k = 0; k < domain_integs.Size(); k++)
{
energy += dnfi[k]->GetElementEnergy(*fe, *T, el_x);
if (domain_integs_marker[k] == nullptr ||
(*(domain_integs_marker[k]))[elem_attr-1] == 1)
{
energy += domain_integs[k]->GetElementEnergy(*fe, *T, el_x);
}
}
}
}
if (fnfi.Size())
if (interior_face_integs.Size())
{
MFEM_ABORT("TODO: add energy contribution from interior face terms");
}
if (bfnfi.Size())
if (boundary_face_integs.Size())
{
MFEM_ABORT("TODO: add energy contribution from boundary face terms");
}
@@ -173,25 +310,41 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
py = 0.0;
if (dnfi.Size())
if (domain_integs.Size())
{
for (int k = 0; k < domain_integs.Size(); k++)
{
if (domain_integs_marker[k] != nullptr)
{
MFEM_VERIFY(mesh->attributes.Size() ==
domain_integs_marker[k]->Size(),
"invalid element marker for domain integrator #"
<< k << ", counting from zero");
}
}
for (int i = 0; i < fes->GetNE(); i++)
{
int elem_attr = mesh->GetAttribute(i);
fe = fes->GetFE(i);
doftrans = fes->GetElementVDofs(i, vdofs);
T = fes->GetElementTransformation(i);
px.GetSubVector(vdofs, el_x);
if (doftrans) {doftrans->InvTransformPrimal(el_x); }
for (int k = 0; k < dnfi.Size(); k++)
for (int k = 0; k < domain_integs.Size(); k++)
{
dnfi[k]->AssembleElementVector(*fe, *T, el_x, el_y);
if (doftrans) {doftrans->TransformDual(el_y); }
py.AddElementVector(vdofs, el_y);
if (domain_integs_marker[k] == nullptr ||
(*(domain_integs_marker[k]))[elem_attr-1] == 1)
{
domain_integs[k]->AssembleElementVector(*fe, *T, el_x, el_y);
if (doftrans) {doftrans->TransformDual(el_y); }
py.AddElementVector(vdofs, el_y);
}
}
}
}
if (fnfi.Size())
if (interior_face_integs.Size())
{
FaceElementTransformations *tr;
const FiniteElement *fe1, *fe2;
@@ -211,16 +364,16 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
fe1 = fes->GetFE(tr->Elem1No);
fe2 = fes->GetFE(tr->Elem2No);
for (int k = 0; k < fnfi.Size(); k++)
for (int k = 0; k < interior_face_integs.Size(); k++)
{
fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
interior_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
py.AddElementVector(vdofs, el_y);
}
}
}
}
if (bfnfi.Size())
if (boundary_face_integs.Size())
{
FaceElementTransformations *tr;
const FiniteElement *fe1, *fe2;
@@ -229,14 +382,14 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < bfnfi.Size(); k++)
for (int k = 0; k < boundary_face_integs.Size(); k++)
{
if (bfnfi_marker[k] == NULL)
if (boundary_face_integs_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *bfnfi_marker[k];
Array<int> &bdr_marker = *boundary_face_integs_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary face integrator #"
<< k << ", counting from zero");
@@ -262,12 +415,12 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
// but we can't dereference a NULL pointer, and we don't want to
// actually make a fake element.
fe2 = fe1;
for (int k = 0; k < bfnfi.Size(); k++)
for (int k = 0; k < boundary_face_integs.Size(); k++)
{
if (bfnfi_marker[k] &&
(*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
if (boundary_face_integs_marker[k] &&
(*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
bfnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
boundary_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
py.AddElementVector(vdofs, el_y);
}
}
@@ -320,26 +473,42 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
*Grad = 0.0;
}
if (dnfi.Size())
if (domain_integs.Size())
{
for (int k = 0; k < domain_integs.Size(); k++)
{
if (domain_integs_marker[k] != nullptr)
{
MFEM_VERIFY(mesh->attributes.Size() ==
domain_integs_marker[k]->Size(),
"invalid element marker for domain integrator #"
<< k << ", counting from zero");
}
}
for (int i = 0; i < fes->GetNE(); i++)
{
int elem_attr = fes->GetMesh()->GetAttribute(i);
fe = fes->GetFE(i);
doftrans = fes->GetElementVDofs(i, vdofs);
T = fes->GetElementTransformation(i);
px.GetSubVector(vdofs, el_x);
if (doftrans) {doftrans->InvTransformPrimal(el_x); }
for (int k = 0; k < dnfi.Size(); k++)
for (int k = 0; k < domain_integs.Size(); k++)
{
dnfi[k]->AssembleElementGrad(*fe, *T, el_x, elmat);
if (doftrans) { doftrans->TransformDual(elmat); }
Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
// Grad->AddSubMatrix(vdofs, vdofs, elmat, 1);
if (domain_integs_marker[k] == nullptr ||
(*(domain_integs_marker[k]))[elem_attr-1] == 1)
{
domain_integs[k]->AssembleElementGrad(*fe, *T, el_x, elmat);
if (doftrans) { doftrans->TransformDual(elmat); }
Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
// Grad->AddSubMatrix(vdofs, vdofs, elmat, 1);
}
}
}
}
if (fnfi.Size())
if (interior_face_integs.Size())
{
FaceElementTransformations *tr;
const FiniteElement *fe1, *fe2;
@@ -359,16 +528,16 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
fe1 = fes->GetFE(tr->Elem1No);
fe2 = fes->GetFE(tr->Elem2No);
for (int k = 0; k < fnfi.Size(); k++)
for (int k = 0; k < interior_face_integs.Size(); k++)
{
fnfi[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
interior_face_integs[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
}
}
}
}
if (bfnfi.Size())
if (boundary_face_integs.Size())
{
FaceElementTransformations *tr;
const FiniteElement *fe1, *fe2;
@@ -377,14 +546,14 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < bfnfi.Size(); k++)
for (int k = 0; k < boundary_face_integs.Size(); k++)
{
if (bfnfi_marker[k] == NULL)
if (boundary_face_integs_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *bfnfi_marker[k];
Array<int> &bdr_marker = *boundary_face_integs_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary face integrator #"
<< k << ", counting from zero");
@@ -410,12 +579,12 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
// but we can't dereference a NULL pointer, and we don't want to
// actually make a fake element.
fe2 = fe1;
for (int k = 0; k < bfnfi.Size(); k++)
for (int k = 0; k < boundary_face_integs.Size(); k++)
{
if (bfnfi_marker[k] &&
(*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
if (boundary_face_integs_marker[k] &&
(*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
bfnfi[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
boundary_face_integs[k]->AssembleFaceGrad(*fe1, *fe2, *tr, el_x, elmat);
Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
}
}
@@ -471,9 +640,9 @@ NonlinearForm::~NonlinearForm()
{
delete cGrad;
delete Grad;
for (int i = 0; i < dnfi.Size(); i++) { delete dnfi[i]; }
for (int i = 0; i < fnfi.Size(); i++) { delete fnfi[i]; }
for (int i = 0; i < bfnfi.Size(); i++) { delete bfnfi[i]; }
for (int i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
for (int i = 0; i < interior_face_integs.Size(); i++) { delete interior_face_integs[i]; }
for (int i = 0; i < boundary_face_integs.Size(); i++) { delete boundary_face_integs[i]; }
delete ext;
}
@@ -564,8 +733,8 @@ BlockNonlinearForm::BlockNonlinearForm(Array<FiniteElementSpace *> &f) :
void BlockNonlinearForm::AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nfi,
Array<int> &bdr_attr_marker)
{
bfnfi.Append(nfi);
bfnfi_marker.Append(&bdr_attr_marker);
boundary_face_integs.Append(nfi);
boundary_face_integs_marker.Append(&bdr_attr_marker);
}
void BlockNonlinearForm::SetEssentialBC(
@@ -600,7 +769,7 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
vdofs[i] = new Array<int>;
}
if (dnfi.Size())
if (domain_integs.Size())
for (int i = 0; i < fes[0]->GetNE(); ++i)
{
T = fes[0]->GetElementTransformation(i);
@@ -612,9 +781,9 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
if (doftrans) {doftrans->InvTransformPrimal(*el_x[s]); }
}
for (int k = 0; k < dnfi.Size(); ++k)
for (int k = 0; k < domain_integs.Size(); ++k)
{
energy += dnfi[k]->GetElementEnergy(fe, *T, el_x_const);
energy += domain_integs[k]->GetElementEnergy(fe, *T, el_x_const);
}
}
@@ -625,12 +794,12 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
delete vdofs[i];
}
if (fnfi.Size())
if (interior_face_integs.Size())
{
MFEM_ABORT("TODO: add energy contribution from interior face terms");
}
if (bfnfi.Size())
if (boundary_face_integs.Size())
{
MFEM_ABORT("TODO: add energy contribution from boundary face terms");
}
@@ -668,7 +837,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
vdofs2[s] = new Array<int>;
}
if (dnfi.Size())
if (domain_integs.Size())
{
for (int i = 0; i < fes[0]->GetNE(); ++i)
{
@@ -681,10 +850,10 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
if (doftrans[s]) {doftrans[s]->InvTransformPrimal(*el_x[s]); }
}
for (int k = 0; k < dnfi.Size(); ++k)
for (int k = 0; k < domain_integs.Size(); ++k)
{
dnfi[k]->AssembleElementVector(fe, *T,
el_x_const, el_y);
domain_integs[k]->AssembleElementVector(fe, *T,
el_x_const, el_y);
for (int s=0; s<fes.Size(); ++s)
{
@@ -696,7 +865,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
}
}
if (fnfi.Size())
if (interior_face_integs.Size())
{
Mesh *mesh = fes[0]->GetMesh();
FaceElementTransformations *tr;
@@ -719,10 +888,10 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
}
for (int k = 0; k < fnfi.Size(); ++k)
for (int k = 0; k < interior_face_integs.Size(); ++k)
{
fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
interior_face_integs[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
for (int s=0; s<fes.Size(); ++s)
{
@@ -734,7 +903,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
}
}
if (bfnfi.Size())
if (boundary_face_integs.Size())
{
Mesh *mesh = fes[0]->GetMesh();
FaceElementTransformations *tr;
@@ -742,14 +911,14 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < bfnfi.Size(); ++k)
for (int k = 0; k < boundary_face_integs.Size(); ++k)
{
if (bfnfi_marker[k] == NULL)
if (boundary_face_integs_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *bfnfi_marker[k];
Array<int> &bdr_marker = *boundary_face_integs_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary face integrator #"
<< k << ", counting from zero");
@@ -776,12 +945,12 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
bx.GetBlock(s).GetSubVector(*(vdofs[s]), *el_x[s]);
}
for (int k = 0; k < bfnfi.Size(); ++k)
for (int k = 0; k < boundary_face_integs.Size(); ++k)
{
if (bfnfi_marker[k] &&
(*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
if (boundary_face_integs_marker[k] &&
(*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
bfnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
boundary_face_integs[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
for (int s=0; s<fes.Size(); ++s)
{
@@ -886,7 +1055,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
}
}
if (dnfi.Size())
if (domain_integs.Size())
{
for (int i = 0; i < fes[0]->GetNE(); ++i)
{
@@ -899,9 +1068,9 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
if (doftrans[s]) {doftrans[s]->InvTransformPrimal(*el_x[s]); }
}
for (int k = 0; k < dnfi.Size(); ++k)
for (int k = 0; k < domain_integs.Size(); ++k)
{
dnfi[k]->AssembleElementGrad(fe, *T, el_x_const, elmats);
domain_integs[k]->AssembleElementGrad(fe, *T, el_x_const, elmats);
for (int j=0; j<fes.Size(); ++j)
{
@@ -920,7 +1089,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
}
}
if (fnfi.Size())
if (interior_face_integs.Size())
{
FaceElementTransformations *tr;
Mesh *mesh = fes[0]->GetMesh();
@@ -941,9 +1110,9 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
}
for (int k = 0; k < fnfi.Size(); ++k)
for (int k = 0; k < interior_face_integs.Size(); ++k)
{
fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
interior_face_integs[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
for (int j=0; j<fes.Size(); ++j)
{
for (int l=0; l<fes.Size(); ++l)
@@ -957,7 +1126,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
}
}
if (bfnfi.Size())
if (boundary_face_integs.Size())
{
FaceElementTransformations *tr;
Mesh *mesh = fes[0]->GetMesh();
@@ -966,14 +1135,14 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < bfnfi.Size(); ++k)
for (int k = 0; k < boundary_face_integs.Size(); ++k)
{
if (bfnfi_marker[k] == NULL)
if (boundary_face_integs_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *bfnfi_marker[k];
Array<int> &bdr_marker = *boundary_face_integs_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary face integrator #"
<< k << ", counting from zero");
@@ -1000,11 +1169,11 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
bx.GetBlock(s).GetSubVector(*vdofs[s], *el_x[s]);
}
for (int k = 0; k < bfnfi.Size(); ++k)
for (int k = 0; k < boundary_face_integs.Size(); ++k)
{
if (bfnfi_marker[k] &&
(*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
bfnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
if (boundary_face_integs_marker[k] &&
(*boundary_face_integs_marker[k])[bdr_attr-1] == 0) { continue; }
boundary_face_integs[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
for (int l=0; l<fes.Size(); ++l)
{
for (int j=0; j<fes.Size(); ++j)
@@ -1109,19 +1278,19 @@ BlockNonlinearForm::~BlockNonlinearForm()
delete ess_tdofs[i];
}
for (int i = 0; i < dnfi.Size(); ++i)
for (int i = 0; i < domain_integs.Size(); ++i)
{
delete dnfi[i];
delete domain_integs[i];
}
for (int i = 0; i < fnfi.Size(); ++i)
for (int i = 0; i < interior_face_integs.Size(); ++i)
{
delete fnfi[i];
delete interior_face_integs[i];
}
for (int i = 0; i < bfnfi.Size(); ++i)
for (int i = 0; i < boundary_face_integs.Size(); ++i)
{
delete bfnfi[i];
delete boundary_face_integs[i];
}
}
+54 -28
View File
@@ -36,14 +36,19 @@ protected:
FiniteElementSpace *fes; // not owned
/// Set of Domain Integrators to be assembled (added).
Array<NonlinearFormIntegrator*> dnfi; // owned
Array<NonlinearFormIntegrator*> domain_integs; // owned
/// Element attribute marker (should be of length mesh->attributes)
/// Includes all by default.
/// 0 - ignore attribute
/// 1 - include attribute
Array<Array<int>*> domain_integs_marker; // not owned
/// Set of interior face Integrators to be assembled (added).
Array<NonlinearFormIntegrator*> fnfi; // owned
Array<NonlinearFormIntegrator*> interior_face_integs; // owned
/// Set of boundary face Integrators to be assembled (added).
Array<NonlinearFormIntegrator*> bfnfi; // owned
Array<Array<int>*> bfnfi_marker; // not owned
Array<NonlinearFormIntegrator*> boundary_face_integs; // owned
Array<Array<int>*> boundary_face_integs_marker; // not owned
mutable SparseMatrix *Grad, *cGrad; // owned
/// Gradient Operator when not assembled as a matrix.
@@ -77,6 +82,20 @@ public:
cP(dynamic_cast<const SparseMatrix*>(P))
{ }
/// Explicitly prohibit copy construction/assignment of NonlinearForms
NonlinearForm(const NonlinearForm &other) = delete;
NonlinearForm& operator=(const NonlinearForm &other) = delete;
/// Move constructor for NonlinearForm.
/** This constructor "steals" the owned data members from the @a other
NonlinearForm. */
NonlinearForm(NonlinearForm &&other);
/// Move assignment operator for NonlinearForm
/** This assignment first frees all owned data, then "steals" the owned data
members from the @a other NonlinearForm. */
NonlinearForm& operator=(NonlinearForm &&other);
/// Set the desired assembly level. The default is AssemblyLevel::LEGACY.
/** For nonlinear operators, the "matrix" assembly levels usually do not make
sense, so only LEGACY, NONE (matrix-free) and PARTIAL are supported.
@@ -106,37 +125,44 @@ public:
FiniteElementSpace *FESpace() { return fes; }
const FiniteElementSpace *FESpace() const { return fes; }
/// Adds new Domain Integrator.
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
{ dnfi.Append(nlfi); }
/** @brief Adds new domain integrator.
@note Assumes ownership of @a nlfi. */
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi);
/** @brief Adds new domain integrator restricted to certain elements
specified by the @a elem_marker.
@note Assumes ownership of @a nlfi. The array @a elem_marker is stored
internally as a pointer to the given Array<int> object. */
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi,
Array<int> &elem_marker);
/// Access all integrators added with AddDomainIntegrator().
Array<NonlinearFormIntegrator*> *GetDNFI() { return &dnfi; }
const Array<NonlinearFormIntegrator*> *GetDNFI() const { return &dnfi; }
Array<NonlinearFormIntegrator*> *GetDNFI() { return &domain_integs; }
const Array<NonlinearFormIntegrator*> *GetDNFI() const { return &domain_integs; }
/// Adds new Interior Face Integrator.
void AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi)
{ fnfi.Append(nlfi); }
void AddInteriorFaceIntegrator(NonlinearFormIntegrator *nlfi);
/** @brief Access all interior face integrators added with
AddInteriorFaceIntegrator(). */
const Array<NonlinearFormIntegrator*> &GetInteriorFaceIntegrators() const
{ return fnfi; }
{ return interior_face_integs; }
/// Adds new Boundary Face Integrator.
void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi)
{ bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
/** @brief Adds new boundary face integrator.
@note Assumes ownership of @a nlfi. */
void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi);
/** @brief Adds new Boundary Face Integrator, restricted to specific boundary
attributes. */
void AddBdrFaceIntegrator(NonlinearFormIntegrator *nfi,
Array<int> &bdr_marker)
{ bfnfi.Append(nfi); bfnfi_marker.Append(&bdr_marker); }
/** @brief Adds new boundary face integrator restricted to specific boundary
attributes by the @a bdr_marker.
@note Assumes ownership of @a nlfi. The array @a bdr_marker is stored
internally as a pointer to the given Array<int> object. */
void AddBdrFaceIntegrator(NonlinearFormIntegrator *nlfi,
Array<int> &bdr_marker);
/** @brief Access all boundary face integrators added with
AddBdrFaceIntegrator(). */
const Array<NonlinearFormIntegrator*> &GetBdrFaceIntegrators() const
{ return bfnfi; }
{ return boundary_face_integs; }
/// Specify essential boundary conditions.
/** This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
@@ -226,14 +252,14 @@ protected:
Array<FiniteElementSpace*> fes;
/// Set of Domain Integrators to be assembled (added).
Array<BlockNonlinearFormIntegrator*> dnfi;
Array<BlockNonlinearFormIntegrator*> domain_integs;
/// Set of interior face Integrators to be assembled (added).
Array<BlockNonlinearFormIntegrator*> fnfi;
Array<BlockNonlinearFormIntegrator*> interior_face_integs;
/// Set of Boundary Face Integrators to be assembled (added).
Array<BlockNonlinearFormIntegrator*> bfnfi;
Array<Array<int>*> bfnfi_marker;
Array<BlockNonlinearFormIntegrator*> boundary_face_integs;
Array<Array<int>*> boundary_face_integs_marker;
/** Auxiliary block-vectors for wrapping input and output vectors or holding
GridFunction-like block-vector data (e.g. in parallel). */
@@ -298,15 +324,15 @@ public:
/// Adds new Domain Integrator.
void AddDomainIntegrator(BlockNonlinearFormIntegrator *nlfi)
{ dnfi.Append(nlfi); }
{ domain_integs.Append(nlfi); }
/// Adds new Interior Face Integrator.
void AddInteriorFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
{ fnfi.Append(nlfi); }
{ interior_face_integs.Append(nlfi); }
/// Adds new Boundary Face Integrator.
void AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
{ bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
{ boundary_face_integs.Append(nlfi); boundary_face_integs_marker.Append(NULL); }
/** @brief Adds new Boundary Face Integrator, restricted to specific boundary
attributes. */
+10 -17
View File
@@ -25,16 +25,12 @@ PANonlinearFormExtension::PANonlinearFormExtension(const NonlinearForm *nlf):
NonlinearFormExtension(nlf),
fes(*nlf->FESpace()),
dnfi(*nlf->GetDNFI()),
elemR(nullptr),
elemR(fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC)),
Grad(*this)
{
if (!DeviceCanUseCeed())
{
elemR = fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
// TODO: optimize for the case when 'elemR' is identity
xe.SetSize(elemR->Height(), Device::GetMemoryType());
ye.SetSize(elemR->Height(), Device::GetMemoryType());
}
// TODO: optimize for the case when 'elemR' is identity
xe.SetSize(elemR->Height(), Device::GetMemoryType());
ye.SetSize(elemR->Height(), Device::GetMemoryType());
ye.UseDevice(true);
}
@@ -139,16 +135,13 @@ void PANonlinearFormExtension::Gradient::Update()
MFNonlinearFormExtension::MFNonlinearFormExtension(const NonlinearForm *form):
NonlinearFormExtension(form), fes(*form->FESpace())
{
if (!DeviceCanUseCeed())
const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
elem_restrict_lex = fes.GetElementRestriction(ordering);
if (elem_restrict_lex) // replace with a check for not identity
{
const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
elem_restrict_lex = fes.GetElementRestriction(ordering);
if (elem_restrict_lex) // replace with a check for not identity
{
localX.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
localY.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
localY.UseDevice(true); // ensure 'localY = 0.0' is done on device
}
localX.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
localY.SetSize(elem_restrict_lex->Height(), Device::GetMemoryType());
localY.UseDevice(true); // ensure 'localY = 0.0' is done on device
}
}
+1 -10
View File
@@ -28,16 +28,7 @@ void VectorConvectionNLFIntegrator::AssemblePA(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedPAVectorConvectionNLIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::PAVectorConvectionNLFIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::PAVectorConvectionNLFIntegrator(fes, *ir, Q);
return;
}
dim = mesh->Dimension();
+1 -10
View File
@@ -28,16 +28,7 @@ void VectorConvectionNLFIntegrator::AssembleMF(const FiniteElementSpace &fes)
if (DeviceCanUseCeed())
{
delete ceedOp;
const bool mixed = mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
fes.IsVariableOrder();
if (mixed)
{
ceedOp = new ceed::MixedMFVectorConvectionNLIntegrator(*this, fes, Q);
}
else
{
ceedOp = new ceed::MFVectorConvectionNLFIntegrator(fes, *ir, Q);
}
ceedOp = new ceed::MFVectorConvectionNLFIntegrator(fes, *ir, Q);
return;
}
MFEM_ABORT("Not yet implemented.");
+49 -30
View File
@@ -19,6 +19,55 @@
namespace mfem
{
ParBilinearForm::ParBilinearForm(ParBilinearForm &&other)
: BilinearForm(std::move(other)), pfes(other.pfes),
Xaux(other.pfes, other.Xaux.GetData()), Yaux(other.pfes, other.Yaux.GetData()),
Ytmp(std::move(other.Ytmp)), p_mat(other.p_mat), p_mat_e(other.p_mat_e),
keep_nbr_block(other.keep_nbr_block)
{
other.Xaux.MakeRef(other.pfes, nullptr);
other.Yaux.MakeRef(other.pfes, nullptr);
p_mat.SetOperatorOwner();
other.p_mat.SetOperatorOwner(false);
other.p_mat.SetType(Operator::Hypre_ParCSR);
p_mat_e.SetOperatorOwner();
other.p_mat_e.SetOperatorOwner(false);
other.p_mat_e.SetType(Operator::Hypre_ParCSR);
other.keep_nbr_block = false;
}
ParBilinearForm& ParBilinearForm::operator=(ParBilinearForm &&other)
{
if (this != &other)
{
BilinearForm::operator=(std::move(other));
pfes = other.pfes;
Xaux.MakeRef(other.pfes, other.Xaux.GetData());
other.Xaux.MakeRef(other.pfes, nullptr);
Yaux.MakeRef(other.pfes, other.Yaux.GetData());
other.Yaux.MakeRef(other.pfes, nullptr);
Ytmp = std::move(other.Ytmp);
p_mat = other.p_mat;
p_mat.SetOperatorOwner();
other.p_mat.SetOperatorOwner(false);
other.p_mat.SetType(Operator::Hypre_ParCSR);
p_mat_e = other.p_mat_e;
p_mat_e.SetOperatorOwner();
other.p_mat_e.SetOperatorOwner(false);
other.p_mat_e.SetType(Operator::Hypre_ParCSR);
keep_nbr_block = other.keep_nbr_block;
other.keep_nbr_block = false;
}
return *this;
}
void ParBilinearForm::pAllocMat()
{
int nbr_size = pfes->GetFaceNbrVSize();
@@ -121,36 +170,6 @@ void ParBilinearForm::pAllocMat()
dof_dof.LoseData();
}
void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
bool steal_loc_A)
{
ParFiniteElementSpace &pfespace = *ParFESpace();
// Create a block diagonal parallel matrix
OperatorHandle A_diag(Operator::Hypre_ParCSR);
A_diag.MakeSquareBlockDiag(pfespace.GetComm(),
pfespace.GlobalVSize(),
pfespace.GetDofOffsets(),
&loc_A);
// Parallel matrix assembly using P^t A P (if needed)
if (IsIdentityProlongation(pfespace.GetProlongationMatrix()))
{
A_diag.SetOperatorOwner(false);
A.Reset(A_diag.As<HypreParMatrix>());
if (steal_loc_A)
{
HypreStealOwnership(*A.As<HypreParMatrix>(), loc_A);
}
}
else
{
OperatorHandle P(Operator::Hypre_ParCSR);
P.ConvertFrom(pfespace.Dof_TrueDof_Matrix());
A.MakePtAP(A_diag, P);
}
}
void ParBilinearForm::ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local)
{
A.Clear();
+7 -18
View File
@@ -44,13 +44,6 @@ protected:
void AssembleSharedFaces(int skip_zeros = 1);
private:
/// Copy construction is not supported; body is undefined.
ParBilinearForm(const ParBilinearForm &);
/// Copy assignment is not supported; body is undefined.
ParBilinearForm &operator=(const ParBilinearForm &);
public:
/// Creates parallel bilinear form associated with the FE space @a *pf.
/** The pointer @a pf is not owned by the newly constructed object. */
@@ -71,6 +64,13 @@ public:
p_mat(Operator::Hypre_ParCSR), p_mat_e(Operator::Hypre_ParCSR)
{ keep_nbr_block = false; }
/// Explicitly prohibit copy construction/assignment of ParBilinearForm
ParBilinearForm(const ParBilinearForm &) = delete;
ParBilinearForm &operator=(const ParBilinearForm &) = delete;
ParBilinearForm(ParBilinearForm &&other);
ParBilinearForm& operator=(ParBilinearForm &&other);
/** When set to true and the ParBilinearForm has interior face integrators,
the local SparseMatrix will include the rows (in addition to the columns)
corresponding to face-neighbor dofs. The default behavior is to disregard
@@ -114,17 +114,6 @@ public:
/** The returned matrix has to be deleted by the caller. */
HypreParMatrix *ParallelAssemble(SparseMatrix *m);
/** @brief Compute parallel RAP operator and store it in @a A as a HypreParMatrix.
@param[in] loc_A The rank-local `SparseMatrix`.
@param[out] A The `OperatorHandle` containing the global `HypreParMatrix`.
@param[in] steal_loc_A Have the `HypreParMatrix` in @a A take ownership of
the memory objects in @a loc_A.
*/
void ParallelRAP(SparseMatrix &loc_A,
OperatorHandle &A,
bool steal_loc_A = false);
/** @brief Returns the matrix assembled on the true dofs, i.e.
@a A = P^t A_local P, in the format (type id) specified by @a A. */
void ParallelAssemble(OperatorHandle &A) { ParallelAssemble(A, mat); }
+7 -15
View File
@@ -194,12 +194,12 @@ void ParFiniteElementSpace::Construct()
void ParFiniteElementSpace::PrintPartitionStats()
{
long long ltdofs = ltdof_size;
long long min_ltdofs, max_ltdofs, sum_ltdofs;
long ltdofs = ltdof_size;
long min_ltdofs, max_ltdofs, sum_ltdofs;
MPI_Reduce(&ltdofs, &min_ltdofs, 1, MPI_LONG_LONG, MPI_MIN, 0, MyComm);
MPI_Reduce(&ltdofs, &max_ltdofs, 1, MPI_LONG_LONG, MPI_MAX, 0, MyComm);
MPI_Reduce(&ltdofs, &sum_ltdofs, 1, MPI_LONG_LONG, MPI_SUM, 0, MyComm);
MPI_Reduce(&ltdofs, &min_ltdofs, 1, MPI_LONG, MPI_MIN, 0, MyComm);
MPI_Reduce(&ltdofs, &max_ltdofs, 1, MPI_LONG, MPI_MAX, 0, MyComm);
MPI_Reduce(&ltdofs, &sum_ltdofs, 1, MPI_LONG, MPI_SUM, 0, MyComm);
if (MyRank == 0)
{
@@ -219,14 +219,14 @@ void ParFiniteElementSpace::PrintPartitionStats()
for (int i = 1; i < NRanks; i++)
{
MPI_Status status;
MPI_Recv(&ltdofs, 1, MPI_LONG_LONG, i, 123, MyComm, &status);
MPI_Recv(&ltdofs, 1, MPI_LONG, i, 123, MyComm, &status);
mfem::out << " " << ltdofs;
}
mfem::out << "\n";
}
else
{
MPI_Send(&ltdofs, 1, MPI_LONG_LONG, 0, 123, MyComm);
MPI_Send(&ltdofs, 1, MPI_LONG, 0, 123, MyComm);
}
}
}
@@ -959,10 +959,6 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
SparseMatrix Pdiag;
P->GetDiag(Pdiag);
R = Transpose(Pdiag);
// The following call ensures that the action of the transpose of P is
// performed fast when HYPRE is built for GPUs.
P->EnsureMultTranspose();
}
HypreParMatrix *ParFiniteElementSpace::GetPartialConformingInterpolation()
@@ -2628,10 +2624,6 @@ int ParFiniteElementSpace
{
*P_ = MakeVDimHypreMatrix(pmatrix, ndofs, num_true_dofs,
dof_offs, tdof_offs);
// The following call ensures that the action of the transpose of *P_ is
// performed fast when HYPRE is built for GPUs.
(*P_)->EnsureMultTranspose();
}
// clean up possible remaining messages in the queue to avoid receiving
+11
View File
@@ -18,6 +18,17 @@
namespace mfem
{
ParLinearForm& ParLinearForm::operator=(ParLinearForm &&other)
{
if (this != &other)
{
LinearForm::operator=(std::move(other));
pfes = other.pfes;
other.pfes = nullptr;
}
return *this;
}
void ParLinearForm::Update(ParFiniteElementSpace *pf)
{
if (pf) { pfes = pf; }
+15 -4
View File
@@ -28,10 +28,6 @@ class ParLinearForm : public LinearForm
protected:
ParFiniteElementSpace *pfes; ///< Points to the same object as #fes
private:
/// Copy construction is not supported; body is undefined.
ParLinearForm(const ParLinearForm &);
public:
/** @brief Create an empty ParLinearForm without an associated
ParFiniteElementSpace.
@@ -64,6 +60,9 @@ public:
ParLinearForm(ParFiniteElementSpace *pf, ParLinearForm * plf)
: LinearForm(pf, plf) { pfes = pf; }
/// Explicitly prohibit copy construction of ParLinearForm
ParLinearForm(const ParLinearForm &other) = delete;
/// Copy assignment. Only the data of the base class Vector is copied.
/** It is assumed that this object and @a rhs use ParFiniteElementSpace%s
that have the same size.
@@ -73,6 +72,18 @@ public:
ParLinearForm &operator=(const ParLinearForm &rhs)
{ return operator=((const Vector &)rhs); }
/// Move constructor for ParLinearForm.
/** This constructor "steals" the owned data members from the @a other
ParLinearForm. */
ParLinearForm(ParLinearForm &&other)
: LinearForm(std::move(other)), pfes(other.pfes)
{ other.pfes = nullptr; }
/// Move assignment operator for ParLinearForm
/** This assignment first frees all owned data, then "steals" the owned data
members from the @a other ParLinearForm. */
ParLinearForm& operator=(ParLinearForm &&other);
ParFiniteElementSpace *ParFESpace() const { return pfes; }
/// Update the object according to the given new FE space @a *pf.
+39 -7
View File
@@ -27,13 +27,45 @@ ParNonlinearForm::ParNonlinearForm(ParFiniteElementSpace *pf)
MFEM_VERIFY(!Serial(), "internal MFEM error");
}
ParNonlinearForm::ParNonlinearForm(ParNonlinearForm &&other)
: NonlinearForm(std::move(other)),
X((ParFiniteElementSpace *)other.fes, other.X.GetData()),
Y((ParFiniteElementSpace *)other.fes, other.Y.GetData()), pGrad(other.pGrad)
{
other.X.MakeRef(other.fes, nullptr);
other.Y.MakeRef(other.fes, nullptr);
pGrad.SetOperatorOwner();
other.pGrad.SetOperatorOwner(false);
other.pGrad.SetType(Operator::Hypre_ParCSR);
}
ParNonlinearForm& ParNonlinearForm::operator=(ParNonlinearForm &&other)
{
if (this != &other)
{
NonlinearForm::operator=(std::move(other));
X.MakeRef(other.fes, other.X.GetData());
other.X.MakeRef(other.fes, nullptr);
Y.MakeRef(other.fes, other.Y.GetData());
other.Y.MakeRef(other.fes, nullptr);
pGrad = other.pGrad;
pGrad.SetOperatorOwner();
other.pGrad.SetOperatorOwner(false);
other.pGrad.SetType(Operator::Hypre_ParCSR);
}
return *this;
}
double ParNonlinearForm::GetParGridFunctionEnergy(const Vector &x) const
{
double loc_energy, glob_energy;
loc_energy = GetGridFunctionEnergy(x);
if (fnfi.Size())
if (interior_face_integs.Size())
{
MFEM_ABORT("TODO: add energy contribution from shared faces");
}
@@ -48,7 +80,7 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
{
NonlinearForm::Mult(x, y); // x --(P)--> aux1 --(A_local)--> aux2
if (fnfi.Size())
if (interior_face_integs.Size())
{
MFEM_VERIFY(!NonlinearForm::ext, "Not implemented (extensions + faces");
// Terms over shared interior faces in parallel.
@@ -78,9 +110,9 @@ void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
X.GetSubVector(vdofs1, el_x.GetData());
X.FaceNbrData().GetSubVector(vdofs2, el_x.GetData() + vdofs1.Size());
for (int k = 0; k < fnfi.Size(); k++)
for (int k = 0; k < interior_face_integs.Size(); k++)
{
fnfi[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
interior_face_integs[k]->AssembleFaceVector(*fe1, *fe2, *tr, el_x, el_y);
aux2.AddElementVector(vdofs1, el_y.GetData());
}
}
@@ -116,7 +148,7 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
if (fnfi.Size() == 0)
if (interior_face_integs.Size() == 0)
{
dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
pfes->GetDofOffsets(), Grad);
@@ -252,7 +284,7 @@ void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const
BlockNonlinearForm::MultBlocked(xs, ys);
if (fnfi.Size() > 0)
if (interior_face_integs.Size() > 0)
{
MFEM_ABORT("TODO: assemble contributions from shared face terms");
}
@@ -332,7 +364,7 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
GetLocalGradient(x); // gradients are stored in 'Grads'
if (fnfi.Size() > 0)
if (interior_face_integs.Size() > 0)
{
MFEM_ABORT("TODO: assemble contributions from shared face terms");
}
+7
View File
@@ -32,6 +32,13 @@ protected:
public:
ParNonlinearForm(ParFiniteElementSpace *pf);
/// Explicitly prohibit copy construction/assignment of ParNonlinearForm
ParNonlinearForm(const ParNonlinearForm &other) = delete;
ParNonlinearForm& operator=(const ParNonlinearForm &other) = delete;
ParNonlinearForm(ParNonlinearForm &&other);
ParNonlinearForm& operator=(ParNonlinearForm &&other);
ParFiniteElementSpace *ParFESpace() const
{ return (ParFiniteElementSpace *)fes; }
+2 -2
View File
@@ -361,8 +361,8 @@ struct TDiffusionKernel<2,2,complex_t>
@param F Jt [M x Dim x SDim x NE] - Jacobian transposed, data member in F
@param Q CoefficientEval<>::Type
@param q CoefficientEval<>::Type::result_t
@param A either [M x Dim*(Dim+1)/2] partially assembled Dim x Dim symm.
matrices, or [M x Dim x Dim] partially assembled Dim x Dim matrices.
@param A [M x Dim*(Dim+1)/2] partially assembled Dim x Dim symm. matrices
@param A [M x Dim x Dim] partially assembled Dim x Dim matrices
*/
template <typename T_result_t, typename Q_t, typename q_t, typename asm_type>
static inline MFEM_ALWAYS_INLINE
+12 -213
View File
@@ -58,39 +58,6 @@ void TMOP_Combo_QualityMetric::AssembleH(const DenseMatrix &Jpt,
}
}
double TMOP_WorstCaseUntangleOptimizer_Metric::EvalW(const DenseMatrix &Jpt)
const
{
double metric_tilde = EvalWBarrier(Jpt);
double metric = metric_tilde;
if (wctype == WorstCaseType::PMean)
{
metric = std::pow(metric_tilde, exponent);
}
else if (wctype == WorstCaseType::Beta)
{
double beta = max_muT+muT_ep;
metric = metric_tilde/(beta-metric_tilde);
}
return metric;
}
double TMOP_WorstCaseUntangleOptimizer_Metric::EvalWBarrier(
const DenseMatrix &Jpt) const
{
double denominator = 1.0;
if (btype == BarrierType::Shifted)
{
denominator = 2.0*(Jpt.Det()-std::min(alpha*min_detT-detT_ep, 0.0));
}
else if (btype == BarrierType::Pseudo)
{
double detT = Jpt.Det();
denominator = detT + std::sqrt(detT*detT + detT_ep*detT_ep);
}
return tmop_metric.EvalW(Jpt)/denominator;
}
double TMOP_Metric_001::EvalW(const DenseMatrix &Jpt) const
{
ie.SetJacobian(Jpt.GetData());
@@ -254,30 +221,6 @@ void TMOP_Metric_002::AssembleH(const DenseMatrix &Jpt,
ie.Assemble_ddI1b(0.5*weight, A.GetData());
}
double TMOP_Metric_004::EvalW(const DenseMatrix &Jpt) const
{
ie.SetJacobian(Jpt.GetData());
return ie.Get_I1() - 2.0*ie.Get_I2b();
}
void TMOP_Metric_004::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
{
ie.SetJacobian(Jpt.GetData());
Add(1.0, ie.Get_dI1(), -2.0, ie.Get_dI2b(), P);
}
void TMOP_Metric_004::AssembleH(const DenseMatrix &Jpt,
const DenseMatrix &DS,
const double weight,
DenseMatrix &A) const
{
ie.SetJacobian(Jpt.GetData());
ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
ie.Assemble_ddI1(weight, A.GetData());
ie.Assemble_ddI2b(-2.0*weight, A.GetData());
}
double TMOP_Metric_007::EvalW(const DenseMatrix &Jpt) const
{
// mu_7 = |J-J^{-t}|^2 = |J|^2 + |J^{-1}|^2 - 4
@@ -3682,18 +3625,20 @@ void TMOP_Integrator::ComputeFDh(const Vector &x, const FiniteElementSpace &fes)
{
if (!fdflag) { return; }
ComputeMinJac(x, fes);
#ifdef MFEM_USE_MPI
const ParFiniteElementSpace *pfes =
dynamic_cast<const ParFiniteElementSpace *>(&fes);
if (pfes)
{
double min_jac_all;
MPI_Allreduce(&dx, &min_jac_all, 1, MPI_DOUBLE, MPI_MIN, pfes->GetComm());
dx = min_jac_all;
}
#endif
}
#ifdef MFEM_USE_MPI
void TMOP_Integrator::ComputeFDh(const Vector &x,
const ParFiniteElementSpace &pfes)
{
if (!fdflag) { return; }
ComputeMinJac(x, pfes);
double min_jac_all;
MPI_Allreduce(&dx, &min_jac_all, 1, MPI_DOUBLE, MPI_MIN, pfes.GetComm());
dx = min_jac_all;
}
#endif
void TMOP_Integrator::EnableFiniteDifferences(const GridFunction &x)
{
fdflag = true;
@@ -3722,152 +3667,6 @@ void TMOP_Integrator::EnableFiniteDifferences(const ParGridFunction &x)
}
#endif
double TMOP_Integrator::ComputeMinDetT(const Vector &x,
const FiniteElementSpace &fes)
{
double min_detT = std::numeric_limits<double>::infinity();
const int NE = fes.GetMesh()->GetNE();
const int dim = fes.GetMesh()->Dimension();
Array<int> xdofs;
Jpr.SetSize(dim);
Jpt.SetSize(dim);
Jrt.SetSize(dim);
for (int i = 0; i < NE; i++)
{
const FiniteElement *fe = fes.GetFE(i);
const IntegrationRule &ir = EnergyIntegrationRule(*fe);
const int dof = fe->GetDof(), nsp = ir.GetNPoints();
DSh.SetSize(dof, dim);
PMatI.SetSize(dof, dim);
Vector posV(dof * dim);
PMatI.UseExternalData(posV.GetData(), dof, dim);
fes.GetElementVDofs(i, xdofs);
x.GetSubVector(xdofs, posV);
DenseTensor Jtr(dim, dim, ir.GetNPoints());
targetC->ComputeElementTargets(i, *fe, ir, posV, Jtr);
for (int q = 0; q < nsp; q++)
{
const IntegrationPoint &ip = ir.IntPoint(q);
const DenseMatrix &Jtr_q = Jtr(q);
CalcInverse(Jtr_q, Jrt);
fe->CalcDShape(ip, DSh);
MultAtB(PMatI, DSh, Jpr);
Mult(Jpr, Jrt, Jpt);
double detT = Jpt.Det();
min_detT = std::min(min_detT, detT);
}
}
return min_detT;
}
double TMOP_Integrator::ComputeUntanglerMaxMuBarrier(const Vector &x,
const FiniteElementSpace &fes)
{
double max_muT = -std::numeric_limits<double>::infinity();
const int NE = fes.GetMesh()->GetNE();
const int dim = fes.GetMesh()->Dimension();
Array<int> xdofs;
Jpr.SetSize(dim);
Jpt.SetSize(dim);
Jrt.SetSize(dim);
TMOP_WorstCaseUntangleOptimizer_Metric *wcuo =
dynamic_cast<TMOP_WorstCaseUntangleOptimizer_Metric *>(metric);
if (!wcuo || wcuo->GetWorstCaseType() !=
TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::Beta)
{
return 0.0;
}
for (int i = 0; i < NE; i++)
{
const FiniteElement *fe = fes.GetFE(i);
const IntegrationRule &ir = EnergyIntegrationRule(*fe);
const int dof = fe->GetDof(), nsp = ir.GetNPoints();
Jpr.SetSize(dim);
Jrt.SetSize(dim);
Jpt.SetSize(dim);
DSh.SetSize(dof, dim);
PMatI.SetSize(dof, dim);
Vector posV(dof * dim);
PMatI.UseExternalData(posV.GetData(), dof, dim);
fes.GetElementVDofs(i, xdofs);
x.GetSubVector(xdofs, posV);
DenseTensor Jtr(dim, dim, ir.GetNPoints());
targetC->ComputeElementTargets(i, *fe, ir, posV, Jtr);
for (int q = 0; q < nsp; q++)
{
const IntegrationPoint &ip = ir.IntPoint(q);
const DenseMatrix &Jtr_q = Jtr(q);
CalcInverse(Jtr_q, Jrt);
fe->CalcDShape(ip, DSh);
MultAtB(PMatI, DSh, Jpr);
Mult(Jpr, Jrt, Jpt);
double metric_val = 0.0;
if (wcuo)
{
wcuo->SetTargetJacobian(Jtr_q);
metric_val = wcuo->EvalWBarrier(Jpt);
}
max_muT = std::max(max_muT, metric_val);
}
}
return max_muT;
}
void TMOP_Integrator::ComputeUntangleMetricQuantiles(const Vector &x,
const FiniteElementSpace &fes)
{
TMOP_WorstCaseUntangleOptimizer_Metric *wcuo =
dynamic_cast<TMOP_WorstCaseUntangleOptimizer_Metric *>(metric);
if (!wcuo) { return; }
#ifdef MFEM_USE_MPI
const ParFiniteElementSpace *pfes =
dynamic_cast<const ParFiniteElementSpace *>(&fes);
#endif
if (wcuo && wcuo->GetBarrierType() ==
TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::Shifted)
{
double min_detT = ComputeMinDetT(x, fes);
double min_detT_all = min_detT;
#ifdef MFEM_USE_MPI
if (pfes)
{
MPI_Allreduce(&min_detT, &min_detT_all, 1, MPI_DOUBLE, MPI_MIN,
pfes->GetComm());
}
#endif
if (wcuo) { wcuo->SetMinDetT(min_detT_all); }
}
double max_muT = ComputeUntanglerMaxMuBarrier(x, fes);
double max_muT_all = max_muT;
#ifdef MFEM_USE_MPI
if (pfes)
{
MPI_Allreduce(&max_muT, &max_muT_all, 1, MPI_DOUBLE, MPI_MAX,
pfes->GetComm());
}
#endif
wcuo->SetMaxMuT(max_muT_all);
}
void TMOPComboIntegrator::EnableLimiting(const GridFunction &n0,
const GridFunction &dist,
Coefficient &w0,
+4 -138
View File
@@ -104,89 +104,6 @@ public:
const double weight, DenseMatrix &A) const;
};
/// Simultaneous Untangler + Worst Case Improvement Metric
/// Uses a base metric mu and is defined as:
/// mu_tilde = mu_hat, when WorstCaseType = None,
/// = mu_hat/(beta - mu_hat), when WorstCaseType = Beta,
/// = mu_hat^p, when WorstCaseType = PMean,
/// where beta = max(mu_hat) + muT_ep,
/// and mu_hat = (mu/2phi(tau,ep)) where
/// 2phi(tau,ep) = 1, when when BarrierType = None,
/// = 2*(tau - min(alpha*min(tau)-detT_ep,0)), when BarrierType = Shifted
/// = tau^2 + sqrt(tau^2 + ep^2), when BarrierType = Pseuso
/// where tau = det(T), and max(mu_hat) and min(tau) are computed over the
/// entire mesh.
/// Ultimately, this metric can be used for mesh untangling with the BarrierType
/// option and for worst case quality improvement with the WorstCaseType option.
class TMOP_WorstCaseUntangleOptimizer_Metric : public TMOP_QualityMetric
{
public:
enum class BarrierType
{
None,
Shifted,
Pseudo
};
enum class WorstCaseType
{
None,
Beta,
PMean
};
protected:
TMOP_QualityMetric &tmop_metric; // non-barrier metric to use
double min_detT; // minimum Jacobian in the mesh
double max_muT; // max mu_k/phi(tau,ep) in the mesh
int exponent; // used for p-mean metrics
double alpha; // scaling factor for min(det(T))
double detT_ep; // small constant subtracted from min(detT)
double muT_ep; // small constant added to muT term
BarrierType btype;
WorstCaseType wctype;
public:
TMOP_WorstCaseUntangleOptimizer_Metric(TMOP_QualityMetric &tmop_metric_,
int exponent_ = 1,
double alpha_ = 1.5,
double detT_ep_ = 0.0001,
double muT_ep_ = 0.0001,
BarrierType btype_ = BarrierType::None,
WorstCaseType wctype_ = WorstCaseType::None) :
tmop_metric(tmop_metric_), exponent(exponent_), alpha(alpha_),
detT_ep(detT_ep_), muT_ep(muT_ep_), btype(btype_), wctype(wctype_)
{
MFEM_VERIFY(wctype == WorstCaseType::None,
"Worst-case optimization has not been fully developed!");
if (btype != BarrierType::None)
{
const int m_id = tmop_metric.Id();
MFEM_VERIFY(m_id == 4 || m_id == 14 || m_id == 66,
"Incorrect input barrier metric -- must be 4 / 14 / 66");
}
}
virtual double EvalW(const DenseMatrix &Jpt) const;
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
{ MFEM_ABORT("Not implemented"); }
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
const double weight, DenseMatrix &A) const
{ MFEM_ABORT("Not implemented"); }
// Compute mu_hat.
virtual double EvalWBarrier(const DenseMatrix &Jpt) const;
virtual void SetMinDetT(double min_detT_) { min_detT = min_detT_; }
virtual void SetMaxMuT(double max_muT_) { max_muT = max_muT_; }
virtual BarrierType GetBarrierType() { return btype; }
virtual WorstCaseType GetWorstCaseType() { return wctype; }
};
/// 2D non-barrier metric without a type.
class TMOP_Metric_001 : public TMOP_QualityMetric
{
@@ -283,24 +200,6 @@ public:
virtual int Id() const { return 2; }
};
/// 2D non-barrier shape (S) metric.
class TMOP_Metric_004 : public TMOP_QualityMetric
{
protected:
mutable InvariantsEvaluator2D<double> ie;
public:
// W = |J|^2 - 2*det(J)
virtual double EvalW(const DenseMatrix &Jpt) const;
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const;
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
const double weight, DenseMatrix &A) const;
virtual int Id() const { return 4; }
};
/// 2D barrier Shape+Size (VS) metric (not polyconvex).
class TMOP_Metric_007 : public TMOP_QualityMetric
{
@@ -437,29 +336,6 @@ public:
const double weight, DenseMatrix &A) const;
};
/// 2D non-barrier Shape+Size (VS) metric.
class TMOP_Metric_066 : public TMOP_Combo_QualityMetric
{
protected:
mutable InvariantsEvaluator2D<double> ie;
double gamma;
TMOP_QualityMetric *sh_metric, *sz_metric;
public:
TMOP_Metric_066(double gamma_) : gamma(gamma_),
sh_metric(new TMOP_Metric_004),
sz_metric(new TMOP_Metric_055)
{
// (1-gamma) mu_4 + gamma mu_55
AddQualityMetric(sh_metric, 1.-gamma_);
AddQualityMetric(sz_metric, gamma_);
}
virtual int Id() const { return 66; }
double GetGamma() const { return gamma; }
virtual ~TMOP_Metric_066() { delete sh_metric; delete sz_metric; }
};
/// 2D barrier size (V) metric (polyconvex).
class TMOP_Metric_077 : public TMOP_QualityMetric
{
@@ -611,7 +487,7 @@ protected:
mutable InvariantsEvaluator3D<double> ie;
public:
// W = |J|^2 / (3 * det(J)^(2/3)) - 1.
// W = |J|^2 / 3 * det(J)^(-2/3) - 1.
virtual double EvalW(const DenseMatrix &Jpt) const;
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const;
@@ -1585,6 +1461,9 @@ protected:
/** @brief Determines the perturbation, h, for FD-based approximation. */
void ComputeFDh(const Vector &x, const FiniteElementSpace &fes);
#ifdef MFEM_USE_MPI
void ComputeFDh(const Vector &x, const ParFiniteElementSpace &pfes);
#endif
void ComputeMinJac(const Vector &x, const FiniteElementSpace &fes);
void UpdateAfterMeshPositionChange(const Vector &new_x);
@@ -1644,13 +1523,6 @@ protected:
void AssemblePA_Limiting();
void ComputeAllElementTargets(const Vector &xe = Vector()) const;
// Compute Min(Det(Jpt)) in the mesh, does not reduce over MPI.
double ComputeMinDetT(const Vector &x, const FiniteElementSpace &fes);
// Compute Max(mu_hat) for the TMOP_WorstCaseUntangleOptimizer_Metric,
// does not reduce over MPI.
double ComputeUntanglerMaxMuBarrier(const Vector &x,
const FiniteElementSpace &fes);
public:
/** @param[in] m TMOP_QualityMetric for r-adaptivity (not owned).
@param[in] tc Target-matrix construction algorithm to use (not owned).
@@ -1843,12 +1715,6 @@ public:
/// Get the surface fitting weight.
double GetSurfaceFittingWeight();
/// Computes quantiles needed for UntangleMetrics. Note that in parallel,
/// the ParFiniteElementSpace must be passed as argument for consistency
/// across MPI ranks.
void ComputeUntangleMetricQuantiles(const Vector &x,
const FiniteElementSpace &fes);
};
class TMOPComboIntegrator : public NonlinearFormIntegrator

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