Compare commits

..
Author SHA1 Message Date
Robert 1b31e07cbd Merge remote-tracking branch 'upstream/feature/partial_ess_bcs' into exaconstit-csr 2022-06-07 11:03:02 -07:00
Robert 63db5c481c fix build issue noted on github ci 2022-06-07 10:45:36 -07:00
Robert da100e4205 make style 2022-06-07 10:34:25 -07:00
Robert 1ae09de13d Merge remote-tracking branch 'origin/master' into feature/partial_ess_bcs 2022-06-07 10:14:43 -07:00
Robert cfcdffd0b1 Removed duplicated code to address reviewers comments
Removed the more or less duplicated FESpace::GetEssentialVDofs and moved the necessary book-keeping portions of things to GetEssentialTrueDofs
Additionally, I added an overwrite optional variable (defaulted to true) to GetEssentialVDofs in-order to allow one to overwrite the values of ess_vdofs if they'd like.
I still need to add some unit tests for the functionality though.
2022-06-07 10:10:31 -07:00
Robert 8216f862ce fix doxygen comments 2022-02-22 14:47:14 -08:00
Robert 6dabc0341e make style 2022-02-22 12:50:21 -08:00
Robert c328515f93 Changes from ExaConstit mfem branch that are useful for having varying partial Essential DOFs per boundary attribute 2022-02-22 12:42:03 -08:00
Robert d8a811cd0f Merge branch 'exaconstit-dev' into exaconstit-csr 2022-02-18 12:22:03 -08:00
Robert 0fb3ee7e98 Merge branch 'master' into exaconstit-dev 2022-02-18 12:21:37 -08:00
Robert Carson 444fa76832 Various fixes needed for BCs usage in exaconstit 2022-02-11 16:29:10 -08:00
Robert 5ebca1fc46 Merge fixes 2021-07-29 11:00:07 -07:00
Robert 6e2a238bfd Merge remote-tracking branch 'upstream/mfem-4.3-dev' into exaconstit-mfem-4.3-test 2021-07-29 10:31:45 -07:00
Robert 87d45dd5a6 Features needed for changing BCs in ExaConstit 2021-01-20 15:29:43 -08:00
Robert 77742567de Merge branch 'master' into exaconstit-dev 2020-12-17 12:45:58 -08:00
Robert 61840db326 Add Assemble Element Assembly hooks into NLFIntegrator class 2020-06-19 16:34:04 -07:00
Robert a365b639b9 Add AssembleDiagonalPA to NonlinearFormIntegrator class 2020-06-18 11:25:14 -07:00
Robert 8b4125f5d6 Merge branch 'master' into exaconstit-dev 2020-06-05 14:57:40 -07:00
Robert 1d3f056641 Merge branch 'master' into exaconstit-dev 2020-03-30 15:13:19 -07:00
Robert a32294af71 Merge branch 'master' into exaconstit-dev 2020-03-13 17:17:13 -07:00
Robert 589c1eeeea Merge branch 'master' into exaconstit-dev 2020-03-08 10:31:02 -07:00
rcarson3 cdd38a6419 Changes related CUDA runs for ExaConstit 2020-02-14 16:38:35 -08:00
rcarson3 d7a56a008a Merge branch 'master' into exaconstit-dev 2020-02-14 16:23:08 -08:00
rcarson3 5e2dccbbc9 NonlinearFormIntegrator addition to allow for Partial Assembly for GetGradient methods 2020-02-03 18:55:29 -08:00
rcarson3 3a296deaff Fix merge issues related IntegrationPoints.index
During the merge there were some leftover parts that used the old integration point ipID rather than the index.
2020-02-03 18:53:10 -08:00
rcarson3 640e35f419 Merge remote-tracking branch 'origin/master' into exaconstit-dev 2020-02-03 17:59:31 -08:00
rcarson3 df2e2fc516 Merge branch 'exaconstit-dev' of https://github.com/rcarson3/mfem into exaconstit-dev 2019-10-16 14:16:04 -07:00
rcarson3 beec7406cb Remove proof of concept xtal plasticity miniapp
The miniapp is no longer needed now that ExaConstit is a full fledge app available on LLNL's github page
2019-10-16 14:12:54 -07:00
rcarson3 a660b511df Merge branch 'master' into exaconstit-dev 2019-10-02 10:45:03 -07:00
Robert Carson ca8363c3e7 Bug fix in VectorFunctionRestrictedCoefficient class
Left in a small bug when updating the VectorFunctionRestrictedCoefficient class
2019-08-14 16:15:39 -07:00
rcarson3 20369095e7 Reduce number of variables passed into VectorFunctionRestrictedCoefficient's TDFunction
This was a class added by ExaConstit to MFEM's main library. It turns out that the function pointer used in that class requires fewer variables than were being used.
2019-08-02 15:30:26 -07:00
rcarson3 acb8766d23 Conduit fixes
Found a few buges. Conduit should now work no matter what. I ended up swapping back to using the MFEM Find_Conduit package to locate Conduit, since SPACK doesn't include all of the nice CMAKE config files.
2019-07-31 15:47:37 -07:00
rcarson3 abe80d02c2 Add support for Conduit post-processing files
Conduit provides a single binary data file that contains all of the outputted variables for each processor which greatly reduces the number of files created when using a large number of processors.
2019-07-31 14:07:26 -07:00
rcarson3 bf0e7c5d22 Merge local branch 'exaconstit-dev-4.0rc1' into exaconstit-dev 2019-07-31 13:47:47 -07:00
rcarson3 b5e5364b63 Better Initial Step Guess
The first initial step for the Exacmech model now constructs the tangent stiffness matrix by making the assumption that there is no velocity applied anywhere. This should result in the purely elastic stiffness matrix being returned. The stress during this step is still solved for with the applied velocity boundary condition. Overall, I have found this to lead to quicker solutions for the larger problem sizes.
2019-07-31 13:39:16 -07:00
rcarson3 f05f805cf9 Rectangular meshes autogenerated now supported
Rectangular meshes can now be specified in the option file and they'll be auto generated.
Also, the grain ID reordering is no longer required since the CA results now return things in the xyz ordering
2019-07-31 13:22:11 -07:00
Robert Carson 03213ce2f7 Minor bug fix to Vector.cpp
Forgot to include this bug fix in my updates.
2019-07-02 20:47:32 -07:00
Robert Carson 0d124d1aa9 Merge pull request #2 from rcarson3/exaconstit-dev-4.0rc1
MFEMv4.0 and initial c++ material model support
2019-07-02 15:49:46 -07:00
rcarson3 7ea8265a83 Bug fixes related to projection/data output
A bug was found in the QuadratureVectorFunction Eval method and was fixed. The projected quaternions are now returned normalized. The default precision for the Visit output has been increased to 12 decimal points. The projected stress is also now returned.
2019-07-02 15:46:16 -07:00
rcarson3 374f9a650c Added more data outputs for newer material models
If the newer C++ models are being used then additional data is outputted using the MFEM visit data format. Some of these formats can't be viewed in Visit itself.
2019-06-22 13:18:50 -07:00
rcarson3 e530f0f13b Bug fix for MTS like model
For the properties a total of 24 were needed instead of 23. I'd missed one while I was counting them all up.
2019-06-20 14:26:45 -07:00
rcarson3 be4dcbeb3d Bug fixes related to ExaCMech integration
Several bug fixes were needed to get at least the Voce ExaCMech model properly running in ExaConstit. The Kocks-Mecking model still hasn't been fully tested at this point.
2019-06-20 11:40:58 -07:00
rcarson3 51f6dbac7c Added plumbing for ExaCMech type models and makefile replaced with CMake plus reorganizing files/dirs
The initial plumbing for ExaCMech has been added. Also, the old makefile was getting tough to  keep up with. Therefore, I've made the move over to using CMake. The current implementation of this is rather ugly and really could use some help to not require the user to pass in where all the different libraries are located. Next, I haven't checked to see if moving over to this build system has broken any of the old code. Further testing is required.
2019-06-16 12:26:06 -07:00
rcarson3 a71548f0ba Merge branch 'master' into exaconstit-dev-4.0rc1 2019-06-12 15:41:53 -07:00
Robert Carson f8f48e1327 Merge pull request #1 from rcarson3/exaconstit-dev-4.0rc1
Exaconstit dev 4.0rc1
2019-04-25 12:38:13 -07:00
rcarson3 1fa911f522 Typo fix 2019-04-25 12:35:31 -07:00
rcarson3 2d06a2dda5 Bug fix related to updated inline mesh generation options
The newer inline mesh generation introduced in the pre-okina tag of MFEM changed how elements were ordered according to a space-filling order. This option was made default in-order to satisfy some issues related to non-conforming meshes. It turns out this caused some diverging material responses in some of our standard monotonic tension tests out to 1% strain when using symmetric BCs. We therefore need to tell the inline mesh generator that we want to use the old lexicographic ordering.
2019-04-25 12:16:28 -07:00
rcarson3 0bab3cb3dd Merge tag 'v4.0-rc1' into exaconstit-dev-4.0rc1
Release candidate #1 for version 4.0
2019-04-24 09:15:52 -07:00
rcarson3 d07f450843 Fix Potential Higher Order Bug
When using higher order elements the mesh wasn't initially promoted to the higher order element. Although, it was later on in the code when the original nodes are swapped out for beggining time step gridfunction. Next, it was noted that the mesh being printed out does not include the higher order mesh... I'm not exactly sure why this is the case.
2019-04-09 16:09:36 -07:00
rcarson3 976809a942 Get rid of compiler warnings and add hydrostatic stress as a visit output
Fixed the various compiler warnings that we were getting. The hydrostatic stress was also added as an output to the Visit files.
2019-04-01 13:09:04 -07:00
rcarson3 09fc9c6744 Modified Newton Raphson method added and some basic timers
A modified newton raphson method has been added in place of the old pure newton raphson. It's nothing fancy it just sets the relaxation factor for the next time step to 0.5 if it notes that the ratio in the norm of the current residual and the previous residual hasn't decreased by a factor of 10 or greater. The relaxation factor is set to 1.0 if it is converging fast enough. The main motivation behind this is to push the solution in the right direction when it starts to oscillate the actual answer. Next, a few basic timers have been set around the entire solution set to provide some very basic profiling for how long each time step takes to solve.
2019-02-26 09:04:06 -08:00
rcarson3 31cc114950 Refactor code related to UMAT in anticipation of writing GPGPU kernels for some of the UMAT code
Several changes were made to the code, but one of the biggest is the computation of all of the incremental and end time step deformation gradients all at once. Next, the local shape function gradiants associated with the initial mesh are saved, so we don't have to constantly swap between the beginning and end time step coordinates. From these changes, a few other functions have been touched. Overall, it looks like this new method does provide more stability when running certain types of xtal plasticity problems that caused the old method to fail. I have a feeling this was related to the end time step deformation gradient being calculated as F_end = F_incr F_beg which over time could lead to problems from floating point errors.
2019-02-14 16:11:20 -08:00
rcarson3 c33e377f75 Fix bug related to ProjectCoefficients for quadrature functions
The integration point ids were not initially set which was causing issues related to quadrature functions to grid functions.
2019-02-08 15:49:24 -08:00
rcarson3 3a6a3bdd35 Moved the NonlinearMechOperator class to its own file
The NonlinearMechOperator had previously been residing in the mechanics_driver file even though it was a fairly lengthy and extensive class. It has since been moved over to the mechanics_operator file, so it should now be easier to navigate and modify both the operator class and the functions/procedures related to the driver.
2019-02-06 12:27:02 -08:00
rcarson3 8559c826a6 Remove hyperelastic model and some minor bug fixes
The hyper elastic model was no longer supported when the code base was swapped over to a velocity based solution field. Therefore, it has been removed from the code base. Hyperelastic type models can still be used in this code through the UMAT interface. Next, a logical bug was introduced when portions of the checkMaterialArgs function were commented out in the last commit. This bug has since been fixed.
2019-02-06 11:47:54 -08:00
rcarson3 17b39bf850 New option file driver replacing old command line option driver
The old command line option driver has been completely replaced by a TOML file option driver. The new option file allowed us to move all of the various options over to a single class. It also allowed for the removal of several checks that were put into place related to the mesh and orientation types being used. Now only one type may be supplied at a time. This commit does rely on several C++11 features in order to work, and these are mainly found in the TOML parser.
2019-02-04 13:18:40 -08:00
rcarson3 e9d9ab94d8 Reduce runtime print output
Several functions were accidentally having all processors print out info. A small change has been made so now only 1 processor should be outputting data.
2019-01-25 14:52:20 -08:00
rcarson3 81baed1c1c Add a first step ramp solver for when solution doesn't converge
It was found with large meshes that the initial applied boundary conditions would not converge. However, the solution would converge if some reduced value of the initial applied boundary conditions is applied. Therefore, an initial step solver is added that ramps up the solution to the desired applied boundary conditions if it doesn't converge initially. After the first step, the old solver is used.
2019-01-24 16:04:27 -08:00
rcarson3 9d75ab3303 Preconditioner update for PCG
After a few discussions on mfem github issue #725, it was discussed that a few preconditioner options were leading to the bug noted in the PCG solver. The new changes seen here seemed to resolve these issues.
2019-01-22 13:36:55 -08:00
rcarson3 75039da7e4 Make GMRES default solver and bug fix to hexmesh BCs
The automatic hexmesh generated from the -hexmesh option was incorrectly relabeling the boundary element attributes over to the ExaConstit convention. This commit fixes this issue. Next, a bug appears to exist in the PCG iterative solver, and therefore the GMRES solver is being set to the default solver until the bug is found and fixed in the PCG solver.
2019-01-18 12:11:07 -08:00
Tzanio fd87a574a0 make style 2019-01-16 10:32:40 -08:00
rcarson3 4e8e1a1276 Merge branch 'master' into exaconstit-dev 2019-01-15 14:09:49 -08:00
rcarson3 69ce10a116 Bug fixes, preconditioner updates, and residual reformulations
A few bug fixes related to the strain calculations were found and fixed. A specific preconditioner is being used based on being used in other xtal plasticity codes that group members have used. The residual has been updated to explicitly solve for the divergence of the cauchy stress. The material stiffness matrix for the UMATs has also been updated so it now has something akin to the d sigma/d def_rate so we get out the divergence of cauchy stress for our [B]^t [C_stiff] [B] {Velocity} calculation.
2019-01-15 13:53:10 -08:00
rcarson3 07e296698c Rebase ExaConstit so it is now a velocity-based FEM code
Rewrote ExaConstit so that it solves for the velocity field instead of the displacement field. The advantages of this are that many crystal plasticity methods use quantities derived from the velocity gradient. Therefore, it's a lot easier to deal with these quantities when your solution field is the velocity field. One can also think of this as changing the code base from an Updated Total Lagrangian formulation to an Updated Lagrangian formulation.
2019-01-08 13:30:05 -08:00
rcarson3 6e814c7e1b Merge branch 'master' into exaconstit-dev 2018-11-16 09:40:00 -08:00
rcarson3 5b64fcc6ac Small bug fix related to how some UMAT variables were being initialized 2018-11-16 09:33:41 -08:00
rcarson3 e3c8fbc5c9 Bug fixes related to BC and reduced number of prints in driver file when running in parallel 2018-11-15 13:58:44 -08:00
rcarson3 bf1d5e8f23 Additional bug fixes, small refactoring, and added functionality 2018-11-12 14:22:09 -08:00
rcarson3 b566f22027 Bug fixes, added a few unit tests, and additional features added 2018-11-06 12:17:57 -08:00
rcarson3 655feea5e6 Merge latest changes from mfem/exaconstit-dev/branch 2018-10-05 13:24:38 -07:00
Steven Robert Wopschall 29239e34bf including bash script to run isochoric compression problem 2018-10-04 21:06:41 -07:00
Steven Robert Wopschall b0e94cfa70 functioning hyperelastic implementation for Robert to use to debug quadrature functions on the ExaModel 2018-10-04 20:59:00 -07:00
rcarson3 3ae3550282 Fixed the GetEssentialVDofs and GetEssentialTrueDofs issue related to examples not compiling 2018-09-21 09:22:34 -07:00
Carson a8dc62a6cf Merged changes from master branch into exaconstit-dev 2018-09-17 11:45:11 -07:00
Steven Robert Wopschall 7ebd03cb7f compiling version of the code with the ExaModel interface implementation 2018-09-05 15:43:05 -07:00
Steven Robert Wopschall b5332c0512 code is not debugged and tested, last commit works, this commit introduces incremental form, new model class, extensive work to umat interface, etc. 2018-08-09 16:10:21 -07:00
Steven Robert Wopschall 0c1103c6b5 added version printout and removed debug prints 2018-04-20 08:22:04 -07:00
Steven Robert Wopschall 2fba048e2b added quadrature functions for cp grain orientations, state variables, stress; added a boundary condition manager to more easily handle essential BCs; added plumbing for user defined material models including a UMAT interface; added a quadrature function coefficient wrapper; tested and debugged serial and parallel refinement and boundary conditions on isochoric cube compression problem using native NeoHookean hyperelastic model 2018-04-19 15:31:08 -07:00
Steven Robert Wopschall 3eed69f9db added a vector set function, a VectorFunctionCoefficient constructor, a GetElement QuadratureFunction routine and the mechanics coefficient and integrator files 2018-03-13 20:49:58 -07:00
Steven Robert Wopschall ed74c45d81 merging branch master into exavconstit-dev
Merge branch 'master' of https://github.com/mfem/mfem into exaconstit-dev
2018-01-25 09:33:02 -08:00
Steven R. Wopschall f36d0903cf fixed merge conflict in driver 2017-09-28 08:54:16 -07:00
Steven R. Wopschall 60ca5e54c1 slight modification to comments in driver and boundary ordering in cube-hex.mesh 2017-09-28 08:52:09 -07:00
Jamie A. Bramwell fa0eac49c6 Fixed run command in driver file 2017-09-27 15:28:56 -07:00
Steven R. Wopschall 97eee74bd6 added cube-hex mesh to be run with ExaConstit miniapp 2017-09-27 14:56:01 -07:00
Steven R. Wopschall d95afacf4b "made changes to gridfunc files"
Merge branch 'exaconstit-dev' of https://github.com/mfem/mfem into exaconstit-dev
2017-09-27 14:43:22 -07:00
Steven R. Wopschall b9e4a76487 modified gridfunc files to get branch up to date 2017-09-27 14:31:59 -07:00
Jamie A. Bramwell cad582afd2 Build fixes and cleanup 2017-09-27 13:15:53 -07:00
Jamie A. Bramwell 229073893d Moving over Steve Wopschall's initial implementation 2017-09-27 12:54:10 -07:00
237 changed files with 3897 additions and 11281 deletions
+2 -4
View File
@@ -1,4 +1,4 @@
name: "Docker"
name: Build Deploy Container
on:
@@ -20,7 +20,6 @@ on:
jobs:
build:
if: github.repository == 'mfem/mfem' # Don't run in forks
permissions:
packages: write
strategy:
@@ -28,8 +27,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
+4 -11
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 "")
@@ -563,9 +556,9 @@ endif()
message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")
target_include_directories(mfem
PUBLIC
${TPL_INCLUDE_DIRS}
$<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
${TPL_INCLUDE_DIRS})
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>)
set_target_properties(mfem PROPERTIES VERSION "${mfem_VERSION}")
set_target_properties(mfem PROPERTIES SOVERSION "${mfem_VERSION}")
-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"
}
+86
View File
@@ -0,0 +1,86 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
# TRIANGLE = 2
# SQUARE = 3
# TETRAHEDRON = 4
# CUBE = 5
#
dimension
3
elements
8
1 5 0 1 10 9 3 4 13 12
1 5 1 2 11 10 4 5 14 13
1 5 3 4 13 12 6 7 16 15
1 5 4 5 14 13 7 8 17 16
1 5 9 10 19 18 12 13 22 21
1 5 10 11 20 19 13 14 23 22
1 5 12 13 22 21 15 16 25 24
1 5 13 14 23 22 16 17 26 25
#
boundary
24
1 3 1 0 9 10
1 3 2 1 10 11
1 3 10 9 18 19
1 3 11 10 19 20
2 3 0 3 12 9
2 3 9 12 21 18
2 3 3 6 15 12
2 3 12 15 24 21
3 3 0 1 4 3
3 3 1 2 5 4
3 3 3 4 7 6
3 3 4 5 8 7
4 3 6 7 16 15
4 3 7 8 17 16
4 3 15 16 25 24
4 3 16 17 26 25
5 3 18 21 22 19
5 3 19 22 23 20
5 3 21 24 25 22
5 3 22 25 26 23
6 3 2 11 14 5
6 3 11 20 23 14
6 3 5 14 17 8
6 3 14 23 26 17
vertices
27
3
0.0 0.0 0.0
0.5 0.0 0.0
1.0 0.0 0.0
0.0 0.0 0.5
0.5 0.0 0.5
1.0 0.0 0.5
0.0 0.0 1.0
0.5 0.0 1.0
1.0 0.0 1.0
0.0 0.5 0.0
0.5 0.5 0.0
1.0 0.5 0.0
0.0 0.5 0.5
0.5 0.5 0.5
1.0 0.5 0.5
0.0 0.5 1.0
0.5 0.5 1.0
1.0 0.5 1.0
0.0 1.0 0.0
0.5 1.0 0.0
1.0 1.0 0.0
0.0 1.0 0.5
0.5 1.0 0.5
1.0 1.0 0.5
0.0 1.0 1.0
0.5 1.0 1.0
1.0 1.0 1.0
+84
View File
@@ -0,0 +1,84 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
# TRIANGLE = 2
# SQUARE = 3
# TETRAHEDRON = 4
# CUBE = 5
#
dimension
3
elements
8
1 5 0 1 4 3 9 10 13 12
1 5 1 2 5 4 10 11 14 13
1 5 9 10 13 12 18 19 22 21
1 5 10 11 14 13 19 20 23 22
1 5 3 4 7 6 12 13 16 15
1 5 4 5 8 7 13 14 17 16
1 5 12 13 16 15 21 22 25 24
1 5 13 14 17 16 22 23 26 25
boundary
24
1 3 0 1 10 9
1 3 1 2 11 10
1 3 9 10 19 18
1 3 10 11 20 19
3 3 0 3 4 1
3 3 1 4 5 2
3 3 3 6 7 4
3 3 4 7 8 5
3 3 18 19 22 21
3 3 19 20 23 22
3 3 21 22 25 24
3 3 22 23 26 25
3 3 2 5 14 11
3 3 11 14 23 20
3 3 5 8 17 14
3 3 14 17 26 23
3 3 0 9 12 3
3 3 9 18 21 12
3 3 3 12 15 6
3 3 12 21 24 15
2 3 6 15 16 7
2 3 7 16 17 8
2 3 15 24 25 16
2 3 16 25 26 17
vertices
27
3
0.0 0.0 0.0
0.5 0.0 0.0
1.0 0.0 0.0
0.0 0.0 0.5
0.5 0.0 0.5
1.0 0.0 0.5
0.0 0.0 1.0
0.5 0.0 1.0
1.0 0.0 1.0
0.0 0.5 0.0
0.5 0.5 0.0
1.0 0.5 0.0
0.0 0.5 0.5
0.5 0.5 0.5
1.0 0.5 0.5
0.0 0.5 1.0
0.5 0.5 1.0
1.0 0.5 1.0
0.0 1.0 0.0
0.5 1.0 0.0
1.0 1.0 0.0
0.0 1.0 0.5
0.5 1.0 0.5
1.0 1.0 0.5
0.0 1.0 1.0
0.5 1.0 1.0
1.0 1.0 1.0
+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());
}
+63 -275
View File
@@ -3,63 +3,34 @@
// Compile with: make ex33
//
// Sample runs: ex33 -m ../data/square-disc.mesh -alpha 0.33 -o 2
// ex33 -m ../data/square-disc.mesh -alpha 4.5 -o 3
// ex33 -m ../data/star.mesh -alpha 1.4 -o 3
// ex33 -m ../data/star.mesh -alpha 0.99 -o 3
// ex33 -m ../data/inline-quad.mesh -alpha 0.5 -o 3
// ex33 -m ../data/amr-quad.mesh -alpha 1.5 -o 3
// ex33 -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3
// ex33 -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4
// ex33 -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4
// ex33 -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5
//
// Verification runs:
// ex33 -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2
// ex33 -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2
// ex33 -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2
// ex33 -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1
//
// Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i)
// for all alpha.
//
// Description:
//
// In this example we solve the following fractional PDE with MFEM:
//
// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α,
// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α < 1,
//
// To solve this FPDE, we apply the operator ( - Δ )^(-N), where the integer
// N is given by floor(α). By doing so, we obtain
// To solve this FPDE, we rely on a rational approximation [2] of the normal
// linear operator A^{-α}, where A = - Δ (with associated homogeneous
// boundary conditions). Namely, we first approximate the operator
//
// ( - Δ )^(α-N) u = ( - Δ )^(-N) f in Ω, u = 0 on ∂Ω, 0 < α.
//
// We first compute the right hand side by solving the integer order PDE
//
// ( - Δ )^N g = f in Ω, g = ( - Δ )^k g = 0 on ∂Ω, k = 1,..,N-1
//
// The remaining FPDE is then given by
//
// ( - Δ )^(α-N) u = g in Ω, u = 0 on ∂Ω.
//
// We rely on a rational approximation [2] of the normal linear operator
// A^{-α + N}, where A = - Δ (with associated homogeneous boundary conditions)
// and (a-N) in (0,1). We approximate the operator
//
// A^{-α+N} ≈ Σ_{i=0}^M c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0,
// A^{-α} ≈ Σ_{i=0}^N c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0,
//
// where I is the L2-identity operator and the coefficients c_i and d_i
// are generated offline to a prescribed accuracy in a pre-processing step.
// We use the triple-A algorithm [1] to generate the rational approximation
// that this partial fractional expansion derives from. We then solve M+1
// that this partial fractional expansion derives from. We then solve N+1
// independent integer-order PDEs,
//
// A u_i + d_i u_i = c_i g in Ω, u_i = 0 on ∂Ω, i=0,...,M,
// A u_i + d_i u_i = c_i f in Ω, u_i = 0 on ∂Ω, i=0,...,N,
//
// using MFEM and sum u_i to arrive at an approximate solution of the FPDE
//
// u ≈ Σ_{i=0}^M u_i.
//
// (If alpha is an integer, we stop after the first PDE was solved.)
// u ≈ Σ_{i=0}^N u_i.
//
// References:
//
@@ -76,8 +47,6 @@
#include "mfem.hpp"
#include <fstream>
#include <iostream>
#include <math.h>
#include <string>
#include "ex33.hpp"
@@ -90,9 +59,8 @@ int main(int argc, char *argv[])
const char *mesh_file = "../data/star.mesh";
int order = 1;
int num_refs = 3;
double alpha = 0.5;
bool visualization = true;
bool verification = false;
double alpha = 0.5;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
@@ -107,9 +75,6 @@ int main(int argc, char *argv[])
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.AddOption(&verification, "-ver", "--verification", "-no-ver",
"--no-verification",
"Use sinusoidal function (f) for analytic comparison.");
args.Parse();
if (!args.Good())
{
@@ -119,31 +84,9 @@ int main(int argc, char *argv[])
args.PrintOptions(cout);
Array<double> coeffs, poles;
int progress_steps = 1;
// 2. Compute the rational expansion coefficients that define the
// integer-order PDEs.
const int power_of_laplace = floor(alpha);
double exponent_to_approximate = alpha - power_of_laplace;
bool integer_order = false;
// Check if alpha is an integer or not.
if (abs(exponent_to_approximate) > 1e-12)
{
mfem::out << "Approximating the fractional exponent "
<< exponent_to_approximate
<< endl;
ComputePartialFractionApproximation(exponent_to_approximate, coeffs,
poles);
// If the example is build without LAPACK, the exponent_to_approximate
// might be modified by the function call above.
alpha = exponent_to_approximate + power_of_laplace;
}
else
{
integer_order = true;
mfem::out << "Treating integer order PDE." << endl;
}
// 2. Compute the coefficients that define the integer-order PDEs.
ComputePartialFractionApproximation(alpha,coeffs,poles);
// 3. Read the mesh from the given mesh file.
Mesh mesh(mesh_file, 1, 1);
@@ -156,8 +99,8 @@ int main(int argc, char *argv[])
}
// 5. Define a finite element space on the mesh.
H1_FECollection fec(order, dim);
FiniteElementSpace fespace(&mesh, &fec);
FiniteElementCollection *fec = new H1_FECollection(order, dim);
FiniteElementSpace fespace(&mesh, fec);
cout << "Number of finite element unknowns: "
<< fespace.GetTrueVSize() << endl;
@@ -171,234 +114,79 @@ int main(int argc, char *argv[])
}
// 7. Define diffusion coefficient, load, and solution GridFunction.
auto func = [&alpha](const Vector &x)
{
double val = 1.0;
for (int i=0; i<x.Size(); i++)
{
val *= sin(M_PI*x(i));
}
return pow(x.Size()*pow(M_PI,2), alpha) * val;
};
FunctionCoefficient f(func);
ConstantCoefficient f(1.0);
ConstantCoefficient one(1.0);
GridFunction u(&fespace);
GridFunction x(&fespace);
GridFunction g(&fespace);
u = 0.0;
x = 0.0;
g = 0.0;
u = 0.;
// 8. Prepare for visualization.
char vishost[] = "localhost";
int visport = 19916;
// 9. Set up the linear form b(.) for integer-order PDE solves.
LinearForm b(&fespace);
if (verification)
socketstream xout, uout;
ostringstream oss_x, oss_u;
if (visualization)
{
// This statement is only relevant for the verification of the code. It
// uses a different f such that an analytic solution is known and easy
// to compare with the numerical one. The FPDE becomes:
// (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
// -> u(x,y) = sin(\pi x) sin(\pi y)
b.AddDomainIntegrator(new DomainLFIntegrator(f));
xout.open(vishost, visport);
xout.precision(8);
uout.open(vishost, visport);
uout.precision(8);
}
else
for (int i = 0; i < coeffs.Size(); i++)
{
b.AddDomainIntegrator(new DomainLFIntegrator(one));
}
b.Assemble();
// 9. Set up the linear form b(.) for integer-order PDE solve.
LinearForm b(&fespace);
ProductCoefficient cf(coeffs[i], f);
b.AddDomainIntegrator(new DomainLFIntegrator(cf));
b.Assemble();
// ------------------------------------------------------------------------
// 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f.
// ------------------------------------------------------------------------
// 10. Define GridFunction for integer-order PDE solve.
GridFunction x(&fespace);
x = 0.0;
if (power_of_laplace > 0)
{
// 10.1 Compute Stiffnes Matrix
BilinearForm k(&fespace);
k.AddDomainIntegrator(new DiffusionIntegrator(one));
k.Assemble();
// 11. Set up the bilinear form a(.,.) for integer-order PDE solve.
BilinearForm a(&fespace);
a.AddDomainIntegrator(new DiffusionIntegrator(one));
ConstantCoefficient c2(-poles[i]);
a.AddDomainIntegrator(new MassIntegrator(c2));
a.Assemble();
// 10.2 Compute Mass Matrix
BilinearForm m(&fespace);
m.AddDomainIntegrator(new MassIntegrator(one));
m.Assemble();
SparseMatrix mass;
Array<int> empty;
m.FormSystemMatrix(empty, mass);
// 10.3 Form the system of equations
// 12. Assemble the bilinear form and the corresponding linear system.
OperatorPtr A;
Vector B, X;
OperatorPtr Op;
k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B);
GSSmoother M((SparseMatrix&)(*Op));
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace
<< " ( f ) " << endl;
for (int i = 0; i < power_of_laplace; i++)
{
// 10.4 Solve the linear system Op X = B (N times).
PCG(*Op, M, B, X, 3, 300, 1e-12, 0.0);
// 13. Solve the linear system A X = B.
GSSmoother M((SparseMatrix&)(*A));
// 10.5 Visualize the solution g of -Δ ^ N g = f in the last step
if (i == power_of_laplace - 1)
{
// Needed for visualization and solution verification.
k.RecoverFEMSolution(X, b, g);
if (integer_order && verification)
{
// For an integer order PDE, g is also our solution u.
u+=g;
}
if (visualization)
{
socketstream fout;
ostringstream oss_f;
fout.open(vishost, visport);
fout.precision(8);
oss_f.str(""); oss_f.clear();
oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ "
<< power_of_laplace
<< " g = f";
fout << "solution\n" << mesh << g
<< "window_title '" << oss_f.str() << "'" << flush;
}
}
mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
<< " u = " << coeffs[i] << " f " << endl;
PCG(*A, M, B, X, 3, 200, 1e-12, 0.0);
// 10.6 Prepare for next iteration (primal / dual space)
mass.Mult(X, B);
X.SetSubVectorComplement(ess_tdof_list,0.0);
}
// 14. Recover the solution as a finite element grid function.
a.RecoverFEMSolution(X, b, x);
// 10.7 Extract solution for the next step. The b now corresponds to the
// function g in the PDE.
const SparseMatrix * R = fespace.GetRestrictionMatrix();
if (R)
{
R->MultTranspose(B,b);
}
else
{
b = B;
}
}
// 15. Accumulate integer-order PDE solutions.
u+=x;
// ------------------------------------------------------------------------
// 11. Solve the fractional PDE by solving M integer order PDEs and adding
// up the solutions.
// ------------------------------------------------------------------------
if (!integer_order)
{
// Setup visualization.
socketstream xout, uout;
ostringstream oss_x, oss_u;
// 16. Send the solutions by socket to a GLVis server.
if (visualization)
{
xout.open(vishost, visport);
xout.precision(8);
uout.open(vishost, visport);
uout.precision(8);
}
// Iterate over all expansion coefficient that contribute to the
// solution.
for (int i = 0; i < coeffs.Size(); i++)
{
mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
<< " u = " << coeffs[i] << " g " << endl;
oss_x.str(""); oss_x.clear();
oss_x << "Solution of PDE -Δ u + " << -poles[i]
<< " u = " << coeffs[i] << " f";
xout << "solution\n" << mesh << x
<< "window_title '" << oss_x.str() << "'" << flush;
// 11.1 Reset GridFunction for integer-order PDE solve.
x = 0.0;
// 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve.
BilinearForm a(&fespace);
a.AddDomainIntegrator(new DiffusionIntegrator(one));
ConstantCoefficient d_i(-poles[i]);
a.AddDomainIntegrator(new MassIntegrator(d_i));
a.Assemble();
// 11.3 Assemble the bilinear form and the corresponding linear system.
OperatorPtr A;
Vector B, X;
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
// 11.4 Solve the linear system A X = B.
GSSmoother M((SparseMatrix&)(*A));
PCG(*A, M, B, X, 3, 300, 1e-12, 0.0);
// 11.5 Recover the solution as a finite element grid function.
a.RecoverFEMSolution(X, b, x);
// 11.6 Accumulate integer-order PDE solutions.
x *= coeffs[i];
u += x;
// 11.7 Send fractional PDE solution to a GLVis server.
if (visualization)
{
oss_x.str(""); oss_x.clear();
oss_x << "Step " << progress_steps
<< ": Solution of PDE -Δ u + " << -poles[i]
<< " u = " << coeffs[i] << " g";
xout << "solution\n" << mesh << x
<< "window_title '" << oss_x.str() << "'" << flush;
oss_u.str(""); oss_u.clear();
oss_u << "Step " << progress_steps + 1
<< ": Solution of fractional PDE (-Δ)^" << alpha
<< " u = f";
uout << "solution\n" << mesh << u
<< "window_title '" << oss_u.str() << "'"
<< flush;
}
oss_u.str(""); oss_u.clear();
oss_u << "Solution of fractional PDE -Δ^" << alpha
<< " u = f";
uout << "solution\n" << mesh << u
<< "window_title '" << oss_u.str() << "'" << flush;
}
}
// ------------------------------------------------------------------------
// 12. (optional) Verify the solution.
// ------------------------------------------------------------------------
if (verification)
{
auto solution = [] (const Vector &x)
{
double val = 1.0;
for (int i=0; i<x.Size(); i++)
{
val *= sin(M_PI*x(i));
}
return val;
};
FunctionCoefficient sol(solution);
double l2_error = u.ComputeL2Error(sol);
string analytic_solution,expected_mesh;
switch (dim)
{
case 1:
analytic_solution = "sin(π x)";
expected_mesh = "inline_segment.mesh";
break;
case 2:
analytic_solution = "sin(π x) sin(π y)";
expected_mesh = "inline_quad.mesh";
break;
default:
analytic_solution = "sin(π x) sin(π y) sin(π z)";
expected_mesh = "inline_hex.mesh";
break;
}
mfem::out << "\n" << string(80,'=')
<< "\n\nSolution Verification in "<< dim << "D \n\n"
<< "Analytic solution : " << analytic_solution << "\n"
<< "Expected mesh : " << expected_mesh <<"\n"
<< "Your mesh : " << mesh_file << "\n"
<< "L2 error : " << l2_error << "\n\n"
<< string(80,'=') << endl;
}
// 17. Free the used memory.
delete fec;
return 0;
}
+4 -15
View File
@@ -32,7 +32,6 @@
#include "mfem.hpp"
#include <fstream>
#include <iostream>
#include <string>
using namespace std;
using namespace mfem;
@@ -250,13 +249,6 @@ void PartialFractionExpansion(double scale, Array<double> & poles,
coeffs.SetSize(psize);
coeffs = scale;
// Note: C p(z)/q(z) = Σ_i c_i / (z - p_i) results in an system of equations
// where the N unknowns are the coefficients c_i. After multiplying the
// system with q(z), the coefficients c_i can be computed analytically by
// choosing N values for z. Choosing z_j = = p_j diagonalizes the system and
// one can obtain an analytic form for the c_i coefficients. The result is
// implemented in the code block below.
for (int i=0; i<psize; i++)
{
double tmp_numer=1.0;
@@ -313,12 +305,9 @@ void ComputePartialFractionApproximation(double & alpha,
if (print_warning)
{
mfem::out
<< "\n" << string(80, '=')
<< "\nMFEM is compiled without LAPACK."
<< "\nUsing precomputed values for PartialFractionApproximation."
<< "\nOnly alpha = 0.33, 0.5, and 0.99 are available."
<< "\nThe default is alpha = 0.5.\n" << string(80, '=') << "\n"
<< endl;
<< "\nMFEM is compiled without LAPACK.\nUsing precomputed values for PartialFractionApproximation. \n"
<< "Only alpha = 0.33, 0.5, and 0.99 are available.\nThe default is alpha = 0.5."
<< std::endl;
}
const double eps = std::numeric_limits<double>::epsilon();
@@ -362,7 +351,7 @@ void ComputePartialFractionApproximation(double & alpha,
if (print_warning)
{
mfem::out << "=> Using precomputed values for alpha = "
mfem::out << "Using precomputed values for alpha = "
<< alpha << "\n" << std::endl;
}
+142 -293
View File
@@ -3,63 +3,34 @@
// Compile with: make ex33p
//
// Sample runs: mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 0.33 -o 2
// mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 4.5 -o 3
// mpirun -np 4 ex33p -m ../data/star.mesh -alpha 1.4 -o 3
// mpirun -np 4 ex33p -m ../data/star.mesh -alpha 0.99 -o 3
// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -alpha 0.5 -o 3
// mpirun -np 4 ex33p -m ../data/amr-quad.mesh -alpha 1.5 -o 3
// mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 -r 2
// mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4
// mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3
// mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4
// mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5
//
// Verification runs:
// mpirun -np 4 ex33p -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2
// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2
// mpirun -np 4 ex33p -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2
// mpirun -np 4 ex33p -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1
// Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i)
// for all alpha.
//
// Description:
//
// In this example we solve the following fractional PDE with MFEM:
//
// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α,
// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α < 1,
//
// To solve this FPDE, we apply the operator ( - Δ )^(-N), where the integer
// N is given by floor(α). By doing so, we obtain
// To solve this FPDE, we rely on a rational approximation [2] of the normal
// linear operator A^{-α}, where A = - Δ (with associated homogeneous
// boundary conditions). Namely, we first approximate the operator
//
// ( - Δ )^(α-N) u = ( - Δ )^(-N) f in Ω, u = 0 on ∂Ω, 0 < α.
//
// We first compute the right hand side by solving the integer order PDE
//
// ( - Δ )^N g = f in Ω, g = ( - Δ )^k g = 0 on ∂Ω, k = 1,..,N-1
//
// The remaining FPDE is then given by
//
// ( - Δ )^(α-N) u = g in Ω, u = 0 on ∂Ω.
//
// We rely on a rational approximation [2] of the normal linear operator
// A^{-α + N}, where A = - Δ (with associated homogeneous boundary conditions)
// and (a-N) in (0,1). We approximate the operator
//
// A^{-α+N} ≈ Σ_{i=0}^M c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0,
// A^{-α} ≈ Σ_{i=0}^N c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0,
//
// where I is the L2-identity operator and the coefficients c_i and d_i
// are generated offline to a prescribed accuracy in a pre-processing step.
// We use the triple-A algorithm [1] to generate the rational approximation
// that this partial fractional expansion derives from. We then solve M+1
// that this partial fractional expansion derives from. We then solve N+1
// independent integer-order PDEs,
//
// A u_i + d_i u_i = c_i g in Ω, u_i = 0 on ∂Ω, i=0,...,M,
// A u_i + d_i u_i = c_i f in Ω, u_i = 0 on ∂Ω, i=0,...,N,
//
// using MFEM and sum u_i to arrive at an approximate solution of the FPDE
//
// u ≈ Σ_{i=0}^M u_i.
//
// (If alpha is an integer, we stop after the first PDE was solved.)
// u ≈ Σ_{i=0}^N u_i.
//
// References:
//
@@ -76,8 +47,6 @@
#include "mfem.hpp"
#include <fstream>
#include <iostream>
#include <math.h>
#include <string>
#include "ex33.hpp"
@@ -96,9 +65,9 @@ int main(int argc, char *argv[])
const char *mesh_file = "../data/star.mesh";
int order = 1;
int num_refs = 3;
double alpha = 0.5;
bool visualization = true;
bool verification = false;
bool visualize_x = false;
double alpha = 0.5;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
@@ -110,12 +79,12 @@ int main(int argc, char *argv[])
"Number of uniform refinements");
args.AddOption(&alpha, "-alpha", "--alpha",
"Fractional exponent");
args.AddOption(&visualize_x, "-vis_x", "--visualize_x", "-no-vis_x",
"--no-visualization_x",
"Enable or disable GLVis visualization of each integer-order PDE solution.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.AddOption(&verification, "-ver", "--verification", "-no-ver",
"--no-verification",
"Use sinusoidal function (f) for analytic comparison.");
"Enable or disable GLVis visualization of the fractional PDE solution.");
args.Parse();
if (!args.Good())
{
@@ -128,51 +97,61 @@ int main(int argc, char *argv[])
}
Array<double> coeffs, poles;
int progress_steps = 1;
// 2. Compute the rational expansion coefficients that define the
// integer-order PDEs.
const int power_of_laplace = floor(alpha);
double exponent_to_approximate = alpha - power_of_laplace;
bool integer_order = false;
// Check if alpha is an integer or not.
if (abs(exponent_to_approximate) > 1e-12)
{
if (Mpi::Root())
{
mfem::out << "Approximating the fractional exponent "
<< exponent_to_approximate
<< endl;
}
ComputePartialFractionApproximation(exponent_to_approximate, coeffs,
poles);
// 2. Compute the coefficients that define the integer-order PDEs.
ComputePartialFractionApproximation(alpha,coeffs,poles);
// If the example is build without LAPACK, the exponent_to_approximate
// might be modified by the function call above.
alpha = exponent_to_approximate + power_of_laplace;
}
else
int num_par_solves;
int max_par_solves = max(1,num_procs/2);
for (num_par_solves=max_par_solves; num_par_solves>0; num_par_solves--)
{
integer_order = true;
if (Mpi::Root())
if (num_procs%num_par_solves==0 && num_par_solves<coeffs.Size())
{
mfem::out << "Treating integer order PDE." << endl;
break;
}
}
if (num_par_solves == 1) {num_par_solves = num_procs;}
// 3. Read the mesh from the given mesh file.
int solver_ranks = num_procs/num_par_solves;
// 3. Split the MPI communicator:
// row_comm is used for parallel partition of the mesh
// col_comm is used for independent integer-order solves
int row_color = myid / solver_ranks; // Determine color based on row
int col_color = myid % solver_ranks; // Determine color based on col
MPI_Comm row_comm, col_comm;
MPI_Comm_split(MPI_COMM_WORLD, row_color, myid, &row_comm);
MPI_Comm_split(MPI_COMM_WORLD, col_color, myid, &col_comm);
int row_rank, row_size, col_rank, col_size;
MPI_Comm_rank(row_comm, &row_rank);
MPI_Comm_size(row_comm, &row_size);
MPI_Comm_rank(col_comm, &col_rank);
MPI_Comm_size(col_comm, &col_size);
if (Mpi::Root())
{
mfem::out << "\nTotal number of MPI ranks = " << num_procs << endl;
mfem::out << "Number of independent parallel solves = " << col_size << endl;
mfem::out << "Number of MPI ranks within each solve = " << row_size
<<"\n" << endl;
}
// 4. Read the mesh from the given mesh file.
Mesh mesh(mesh_file, 1, 1);
int dim = mesh.Dimension();
// 4. Refine the mesh to increase the resolution.
// 5. Refine the mesh to increase the resolution.
for (int i = 0; i < num_refs; i++)
{
mesh.UniformRefinement();
}
ParMesh pmesh(MPI_COMM_WORLD, mesh);
ParMesh pmesh(row_comm, mesh);
mesh.Clear();
// 5. Define a finite element space on the mesh.
// 6. Define a finite element space on the mesh.
H1_FECollection fec(order, dim);
ParFiniteElementSpace fespace(&pmesh, &fec);
if (Mpi::Root())
@@ -181,7 +160,7 @@ int main(int argc, char *argv[])
<< fespace.GetTrueVSize() << endl;
}
// 6. Determine the list of true (i.e. conforming) essential boundary dofs.
// 7. Determine the list of true (i.e. conforming) essential boundary dofs.
Array<int> ess_tdof_list;
if (pmesh.bdr_attributes.Size())
{
@@ -190,250 +169,120 @@ int main(int argc, char *argv[])
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
// 7. Define diffusion coefficient, load, and solution GridFunction.
auto func = [&alpha](const Vector &x)
{
double val = 1.0;
for (int i=0; i<x.Size(); i++)
{
val *= sin(M_PI*x(i));
}
return pow(x.Size()*pow(M_PI,2), alpha) * val;
};
FunctionCoefficient f(func);
// 8. Define diffusion coefficient, load, and solution GridFunction.
ConstantCoefficient f(1.0);
ConstantCoefficient one(1.0);
ParGridFunction u(&fespace);
ParGridFunction x(&fespace);
ParGridFunction g(&fespace);
u = 0.0;
x = 0.0;
g = 0.0;
// 8. Prepare for visualization.
char vishost[] = "localhost";
int visport = 19916;
// 9. Set up the linear form b(.) for integer-order PDE solves.
ParLinearForm b(&fespace);
if (verification)
{
// This statement is only relevant for the verification of the code. It
// uses a different f such that an analytic solution is known and easy
// to compare with the numerical one. The FPDE becomes:
// (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
// -> u(x,y) = sin(\pi x) sin(\pi y)
b.AddDomainIntegrator(new DomainLFIntegrator(f));
}
else
{
b.AddDomainIntegrator(new DomainLFIntegrator(one));
}
b.AddDomainIntegrator(new DomainLFIntegrator(f));
b.Assemble();
// ------------------------------------------------------------------------
// 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f.
// ------------------------------------------------------------------------
if (power_of_laplace > 0)
int my_coeff_size = max(coeffs.Size()/col_size,1);
int ibeg = col_rank*my_coeff_size;
if (ibeg + 2*my_coeff_size > coeffs.Size())
{
// 10.1 Compute Stiffnes Matrix
ParBilinearForm k(&fespace);
k.AddDomainIntegrator(new DiffusionIntegrator(one));
k.Assemble();
my_coeff_size = coeffs.Size()-col_rank*my_coeff_size;
}
else if (ibeg > coeffs.Size() - 1)
{
my_coeff_size = 0;
}
// 10.2 Compute Mass Matrix
ParBilinearForm m(&fespace);
m.AddDomainIntegrator(new MassIntegrator(one));
m.Assemble();
HypreParMatrix mass;
Array<int> empty;
m.FormSystemMatrix(empty, mass);
int iend = ibeg+my_coeff_size;
// 10.3 Form the system of equations
for (int i = ibeg; i < iend; i++)
{
// 10. Reset GridFunction for integer-order PDE solve.
x = 0.0;
// 11. Set up the bilinear form a(.,.) for integer-order PDE solve.
ParBilinearForm a(&fespace);
a.AddDomainIntegrator(new DiffusionIntegrator(one));
ConstantCoefficient d_i(-poles[i]);
a.AddDomainIntegrator(new MassIntegrator(d_i));
a.Assemble();
// 12. Assemble the bilinear form and the corresponding linear system.
OperatorPtr A;
Vector B, X;
OperatorPtr Op;
k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B);
HypreBoomerAMG prec;
prec.SetPrintLevel(-1);
CGSolver cg(MPI_COMM_WORLD);
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
// 13. Solve the linear system A X = B.
HypreBoomerAMG * prec = new HypreBoomerAMG;
prec->SetPrintLevel(-1);
int print_level = (col_rank==0) ? 3 : 0;
if (Mpi::Root())
{
mfem::out << "\nMPI rank " << myid
<< ": Solving PDE -Δ u + " << -poles[i]
<< " u = " << coeffs[i] << " f " << endl;
}
CGSolver cg(row_comm);
cg.SetRelTol(1e-12);
cg.SetMaxIter(2000);
cg.SetPrintLevel(3);
cg.SetPreconditioner(prec);
cg.SetOperator(*Op);
cg.SetPrintLevel(print_level);
cg.SetPreconditioner(*prec);
cg.SetOperator(*A);
cg.Mult(B, X);
delete prec;
if (Mpi::Root())
// 14. Recover the solution as a finite element grid function.
a.RecoverFEMSolution(X, b, x);
// 15. Accumulate integer-order PDE solutions.
x *= coeffs[i];
u += x;
// 16. Send integer-order PDE solutions to a GLVis server.
if (visualize_x)
{
mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace
<< " ( f ) " << endl;
}
for (int i = 0; i < power_of_laplace; i++)
{
// 10.4 Solve the linear system Op X = B (N times).
cg.Mult(B, X);
// 10.5 Visualize the solution g of -Δ ^ N g = f in the last step
if (i == power_of_laplace - 1)
if (col_rank > 0 && i < iend-1)
{
// Needed for visualization and solution verification.
k.RecoverFEMSolution(X, b, g);
if (integer_order && verification)
{
// For an integer order PDE, g is also our solution u.
u+=g;
}
if (visualization)
{
socketstream fout;
ostringstream oss_f;
fout.open(vishost, visport);
fout.precision(8);
oss_f.str(""); oss_f.clear();
oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ "
<< power_of_laplace
<< " g = f";
fout << "parallel " << num_procs << " " << myid << "\n"
<< "solution\n" << pmesh << g
<< "window_title '" << oss_f.str() << "'" << flush;
}
MPI_Status status;
MPI_Recv(nullptr,0,MPI_INT, col_rank-1,0,col_comm,&status);
}
// 10.6 Prepare for next iteration (primal / dual space)
mass.Mult(X, B);
X.SetSubVectorComplement(ess_tdof_list,0.0);
}
// 10.7 Extract solution for the next step. The b now corresponds to the
// function g in the PDE.
const SparseMatrix* rm = fespace.GetRestrictionMatrix();
rm->MultTranspose(B, b);
}
// ------------------------------------------------------------------------
// 11. Solve the fractional PDE by solving M integer order PDEs and adding
// up the solutions.
// ------------------------------------------------------------------------
if (!integer_order)
{
// Setup visualization.
socketstream xout, uout;
ostringstream oss_x, oss_u;
if (visualization)
{
xout.open(vishost, visport);
char vishost[] = "localhost";
int visport = 19916;
socketstream xout(vishost, visport);
xout.precision(8);
uout.open(vishost, visport);
uout.precision(8);
}
// Iterate over all expansion coefficient that contribute to the
// solution.
for (int i = 0; i < coeffs.Size(); i++)
{
if (Mpi::Root())
ostringstream oss;
oss << "Solution of PDE -Δ u + " << -poles[i]
<< " u = " << coeffs[i] << " f" ;
xout << "parallel " << row_size << " " << row_rank << "\n";
xout << "solution\n" << pmesh << x
<< "window_title '" << oss.str() << "'" << flush;
if (col_rank < col_size-1)
{
mfem::out << "\nSolving PDE -Δ u + " << -poles[i]
<< " u = " << coeffs[i] << " g " << endl;
}
// 11.1 Reset GridFunction for integer-order PDE solve.
x = 0.0;
// 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve.
ParBilinearForm a(&fespace);
a.AddDomainIntegrator(new DiffusionIntegrator(one));
ConstantCoefficient d_i(-poles[i]);
a.AddDomainIntegrator(new MassIntegrator(d_i));
a.Assemble();
// 11.3 Assemble the bilinear form and the corresponding linear system.
OperatorPtr A;
Vector B, X;
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
// 11.4 Solve the linear system A X = B.
HypreBoomerAMG prec;
prec.SetPrintLevel(-1);
CGSolver cg(MPI_COMM_WORLD);
cg.SetRelTol(1e-12);
cg.SetMaxIter(2000);
cg.SetPrintLevel(3);
cg.SetPreconditioner(prec);
cg.SetOperator(*A);
cg.Mult(B, X);
// 11.5 Recover the solution as a finite element grid function.
a.RecoverFEMSolution(X, b, x);
// 11.6 Accumulate integer-order PDE solutions.
x *= coeffs[i];
u += x;
// 11.7 Send fractional PDE solution to a GLVis server.
if (visualization)
{
oss_x.str(""); oss_x.clear();
oss_x << "Step " << progress_steps
<< ": Solution of PDE -Δ u + " << -poles[i]
<< " u = " << coeffs[i] << " g";
xout << "parallel " << num_procs << " " << myid << "\n"
<< "solution\n" << pmesh << x
<< "window_title '" << oss_x.str() << "'" << flush;
oss_u.str(""); oss_u.clear();
oss_u << "Step " << progress_steps + 1
<< ": Solution of fractional PDE (-Δ)^" << alpha
<< " u = f";
uout << "parallel " << num_procs << " " << myid << "\n"
<< "solution\n" << pmesh << u
<< "window_title '" << oss_u.str() << "'"
<< flush;
MPI_Send(nullptr,0,MPI_INT,col_rank+1,0,col_comm);
}
}
}
// ------------------------------------------------------------------------
// 12. (optional) Verify the solution.
// ------------------------------------------------------------------------
if (verification)
// 17. Accumulate for the fractional PDE solution
MPI_Allreduce(MPI_IN_PLACE, u.GetData(), u.Size(),
MPI_DOUBLE, MPI_SUM,col_comm);
// 18. Send fractional PDE solution to a GLVis server.
if (visualization)
{
auto solution = [] (const Vector &x)
if (col_rank == 0)
{
double val = 1.0;
for (int i=0; i<x.Size(); i++)
{
val *= sin(M_PI*x(i));
}
return val;
};
FunctionCoefficient sol(solution);
double l2_error = u.ComputeL2Error(sol);
if (Mpi::Root())
{
string analytic_solution,expected_mesh;
switch (dim)
{
case 1:
analytic_solution = "sin(π x)";
expected_mesh = "inline_segment.mesh";
break;
case 2:
analytic_solution = "sin(π x) sin(π y)";
expected_mesh = "inline_quad.mesh";
break;
default:
analytic_solution = "sin(π x) sin(π y) sin(π z)";
expected_mesh = "inline_hex.mesh";
break;
}
mfem::out << "\n" << string(80,'=')
<< "\n\nSolution Verification in "<< dim << "D \n\n"
<< "Analytic solution : " << analytic_solution << "\n"
<< "Expected mesh : " << expected_mesh <<"\n"
<< "Your mesh : " << mesh_file << "\n"
<< "L2 error : " << l2_error << "\n\n"
<< string(80,'=') << endl;
char vishost[] = "localhost";
int visport = 19916;
socketstream uout(vishost, visport);
uout.precision(8);
ostringstream oss;
oss << "Solution of fractional PDE -Δ^" << alpha
<< " u = f" ;
uout << "parallel " << row_size << " " << row_rank << "\n";
uout << "solution\n" << pmesh << u
<< "window_title '" << oss.str() << "'" << flush;
}
}
+4 -40
View File
@@ -305,66 +305,38 @@ public:
/// Finalizes the matrix initialization.
virtual void Finalize(int skip_zeros = 1);
/** @brief Returns a const reference to the sparse matrix: \f$ M \f$
This will fail if HasSpMat() is false. */
/// Returns a const reference to the sparse matrix.
const SparseMatrix &SpMat() const
{
MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
return *mat;
}
/** @brief Returns a reference to the sparse matrix: \f$ M \f$
This will fail if HasSpMat() is false. */
/// Returns a reference to the sparse matrix: \f$ M \f$
SparseMatrix &SpMat()
{
MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
return *mat;
}
/** @brief Returns true if the sparse matrix is not null, false otherwise.
@sa SpMat(). */
bool HasSpMat()
{
return mat != nullptr;
}
/** @brief Nullifies the internal matrix \f$ M \f$ and returns a pointer
to it. Used for transfering ownership. */
SparseMatrix *LoseMat() { SparseMatrix *tmp = mat; mat = NULL; return tmp; }
/** @brief Returns a const reference to the sparse matrix of eliminated b.c.:
\f$ M_e \f$
This will fail if HasSpMatElim() is false. */
/// Returns a const reference to the sparse matrix of eliminated b.c.: \f$ M_e \f$
const SparseMatrix &SpMatElim() const
{
MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
return *mat_e;
}
/** @brief Returns a reference to the sparse matrix of eliminated b.c.:
\f$ M_e \f$
This will fail if HasSpMatElim() is false. */
/// Returns a reference to the sparse matrix of eliminated b.c.: \f$ M_e \f$
SparseMatrix &SpMatElim()
{
MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
return *mat_e;
}
/** @brief Returns true if the sparse matrix of eliminated b.c.s is not null,
false otherwise.
@sa SpMatElim(). */
bool HasSpMatElim()
{
return mat_e != nullptr;
}
/// Adds new Domain Integrator. Assumes ownership of @a bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi);
/// Adds new Domain Integrator restricted to certain elements specified by
@@ -438,14 +410,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 -12
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();
@@ -1506,7 +1497,6 @@ static void PAConvectionApply(const int dim,
{
switch ((D1D << 4 ) | Q1D)
{
case 0x22: return SmemPAConvectionApply3D<2,2>(NE,B,G,Bt,Gt,op,x,y);
case 0x23: return SmemPAConvectionApply3D<2,3>(NE,B,G,Bt,Gt,op,x,y);
case 0x24: return SmemPAConvectionApply3D<2,4>(NE,B,G,Bt,Gt,op,x,y);
case 0x26: return SmemPAConvectionApply3D<2,6>(NE,B,G,Bt,Gt,op,x,y);
@@ -1558,7 +1548,6 @@ static void PAConvectionApplyT(const int dim,
{
switch ((D1D << 4 ) | Q1D)
{
case 0x22: return SmemPAConvectionApplyT3D<2,2>(NE,B,G,Bt,Gt,op,x,y);
case 0x23: return SmemPAConvectionApplyT3D<2,3>(NE,B,G,Bt,Gt,op,x,y);
case 0x24: return SmemPAConvectionApplyT3D<2,4>(NE,B,G,Bt,Gt,op,x,y);
case 0x26: return SmemPAConvectionApplyT3D<2,6>(NE,B,G,Bt,Gt,op,x,y);
+1 -6
View File
@@ -136,9 +136,6 @@ static void PADGTraceSetup(const int dim,
void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
{
const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
Device::GetDeviceMemoryType() : pa_mt;
nf = fes.GetNFbyType(type);
if (nf==0) { return; }
// Assumes tensor-product elements
@@ -156,7 +153,7 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
geom = mesh->GetFaceGeometricFactors(
*ir,
FaceGeometricFactors::DETERMINANTS |
FaceGeometricFactors::NORMALS, type, mt);
FaceGeometricFactors::NORMALS, type);
maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
dofs1D = maps->ndof;
quad1D = maps->nqpt;
@@ -698,7 +695,6 @@ static void PADGTraceApply(const int dim,
{
switch ((D1D << 4 ) | Q1D)
{
case 0x22: return SmemPADGTraceApply3D<2,2,1>(NF,B,Bt,op,x,y);
case 0x23: return SmemPADGTraceApply3D<2,3,1>(NF,B,Bt,op,x,y);
case 0x34: return SmemPADGTraceApply3D<3,4,2>(NF,B,Bt,op,x,y);
case 0x45: return SmemPADGTraceApply3D<4,5,2>(NF,B,Bt,op,x,y);
@@ -1128,7 +1124,6 @@ static void PADGTraceApplyTranspose(const int dim,
{
switch ((D1D << 4 ) | Q1D)
{
case 0x22: return SmemPADGTraceApplyTranspose3D<2,2>(NF,B,Bt,op,x,y);
case 0x23: return SmemPADGTraceApplyTranspose3D<2,3>(NF,B,Bt,op,x,y);
case 0x34: return SmemPADGTraceApplyTranspose3D<3,4>(NF,B,Bt,op,x,y);
case 0x45: return SmemPADGTraceApplyTranspose3D<4,5>(NF,B,Bt,op,x,y);
+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];
+24
View File
@@ -517,6 +517,30 @@ void VectorRestrictedCoefficient::Eval(
}
}
void VectorFunctionRestrictedCoefficient::Eval(Vector &V,
ElementTransformation &T,
const IntegrationPoint &ip)
{
double x[3];
Vector transip(x, 3);
T.Transform(ip, transip);
V.SetSize(vdim);
if (active_attr[T.Attribute-1])
{
(*TDFunction)(T.Attribute, V);
}
else
{
V = 0.0;
}
if (Q)
{
V *= Q->Eval(T, ip, GetTime());
}
}
void PWMatrixCoefficient::InitMap(const Array<int> & attr,
const Array<MatrixCoefficient*> & coefs)
{
+38
View File
@@ -439,12 +439,19 @@ public:
/// Initialize the VectorCoefficient with vector dimension @a vd.
VectorCoefficient(int vd) { vdim = vd; time = 0.; }
// add constructor with no argument list (this is for the
// QuadratureVectorFunctionCoefficient extended class), srw
VectorCoefficient() { vdim = 0; time = 0.; }
/// Set the time for time dependent coefficients
virtual void SetTime(double t) { time = t; }
/// Get the time for time dependent coefficients
double GetTime() { return time; }
/// Set the dimension of the vector.
void SetVDim(int dim) { vdim = dim; }
/// Returns dimension of the vector.
int GetVDim() { return vdim; }
@@ -870,6 +877,37 @@ public:
const IntegrationRule &ir);
};
/// VectorFunctionCoefficient defined on a subset of domain or boundary attributes, srw
class VectorFunctionRestrictedCoefficient : public VectorCoefficient
{
private:
void (*TDFunction)(int, Vector &);
const Array<int> &active_attr;
const Array2D<double> &active_scale;
Coefficient *Q;
public:
/// Construct a time-dependent vector coefficient from a C-function
VectorFunctionRestrictedCoefficient(int dim,
void (*TDF)(int, Vector &),
Array<int> &attr,
Array2D<double> &scale,
Coefficient *q = NULL)
: VectorCoefficient(dim), Q(q), active_attr(attr), active_scale(scale)
{
TDFunction = TDF;
}
using VectorCoefficient::Eval;
virtual void Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip);
virtual ~VectorFunctionRestrictedCoefficient() { }
const Array<int> &GetActiveAttr() { return active_attr; }
const Array2D<double> &GetActiveScale() { return active_scale; }
};
typedef VectorCoefficient DiagonalMatrixCoefficient;
/// Base class for Matrix Coefficients that optionally depend on time and space.
+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'; }
}
+53 -3
View File
@@ -494,12 +494,21 @@ static void mark_dofs(const Array<int> &dofs, Array<int> &mark_array)
void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
Array<int> &ess_vdofs,
int component) const
int component,
bool overwrite) const
{
Array<int> vdofs, dofs;
ess_vdofs.SetSize(GetVSize());
ess_vdofs = 0;
if (overwrite)
{
ess_vdofs.SetSize(GetVSize());
ess_vdofs = 0;
}
else
{
MFEM_ASSERT(ess_vdofs.Size() == GetVSize(),
"ess_vdofs size is not equal to FESpaces GetVSize().");
}
for (int i = 0; i < GetNBE(); i++)
{
@@ -579,6 +588,47 @@ void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
MarkerToList(ess_tdofs, ess_tdof_list);
}
void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
Array<int> &ess_tdof_list,
const Array2D<bool> &component)
{
MFEM_ASSERT(component.NumCols() == vdim,
"Number of columns of component was not equal to FESpace vdim");
MFEM_ASSERT(component.NumRows() == bdr_attr_is_ess.Size(),
"Number of rows of component was not equal to bdr_attr_is_ess.Size()");
Array<int> ess_vdofs, ess_tdofs, bdr_attr_is_ess_single_comp;
bdr_attr_is_ess_single_comp.SetSize(bdr_attr_is_ess.Size());
for (int i = 0; i < vdim; i++)
{
// Only overwrite ess_vdofs on first iteration
// all other iterations we want to preserve values of
// ess_vdofs.
const bool overwrite = (i == 0) ? true : false;
bdr_attr_is_ess_single_comp = 0;
for (int j = 0; j < bdr_attr_is_ess.Size(); j++)
{
if (bdr_attr_is_ess[j] && component(j, i))
{
bdr_attr_is_ess_single_comp[j] = bdr_attr_is_ess[j];
}
}
GetEssentialVDofs(bdr_attr_is_ess_single_comp, ess_vdofs, i, overwrite);
}
const SparseMatrix *R = GetConformingRestriction();
if (!R)
{
ess_tdofs.MakeRef(ess_vdofs);
}
else
{
R->BooleanMult(ess_vdofs, ess_tdofs);
}
MarkerToList(ess_tdofs, ess_tdof_list);
}
void FiniteElementSpace::GetBoundaryTrueDofs(Array<int> &boundary_dofs,
int component)
{
+26 -7
View File
@@ -775,11 +775,16 @@ public:
/** @brief Mark degrees of freedom associated with boundary elements with
the specified boundary attributes (marked in 'bdr_attr_is_ess').
For spaces with 'vdim' > 1, the 'component' parameter can be used
to restricts the marked vDOFs to the specified component. */
to restricts the marked vDOFs to the specified component.
If overwrite is set to false then values in ess_vdofs are preserved
and not reset. However, the assumption here is that ess_vdofs is set to
the correct size already.*/
virtual void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
Array<int> &ess_vdofs,
int component = -1) const;
int component = -1,
bool overwrite = true) const;
/** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
boundary attributes marked in the array bdr_attr_is_ess.
@@ -789,6 +794,20 @@ public:
Array<int> &ess_tdof_list,
int component = -1);
/** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
boundary attributes marked in the array bdr_attr_is_ess.
For spaces with 'vdim' > 1, the 'component' array can be used
to restricts the marked tDOFs per boundary to the specified components.
If vdim > 1 then one can specify per boundary attribute which components
on a boundary are essential by assigning a value of true to its location
in the component array.
The component has dimensions number of boundary attributes x vdim. */
virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
Array<int> &ess_tdof_list,
const Array2D<bool> &component);
/** @brief Get a list of all boundary true dofs, @a boundary_dofs. For spaces
with 'vdim' > 1, the 'component' parameter can be used to restricts the
marked tDOFs to the specified component. Equivalent to
@@ -980,14 +999,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;
}
}
+98 -2
View File
@@ -2025,6 +2025,7 @@ void GridFunction::AccumulateAndCountZones(Coefficient &coeff,
// Local interpolation
Array<int> vdofs;
Vector vals;
this->HostReadWrite();
*this = 0.0;
HostReadWrite();
@@ -2066,6 +2067,7 @@ void GridFunction::AccumulateAndCountZones(VectorCoefficient &vcoeff,
// Local interpolation
Array<int> vdofs;
Vector vals;
this->HostReadWrite();
*this = 0.0;
HostReadWrite();
@@ -2672,6 +2674,100 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr)
#endif
}
void GridFunction::ProjectBdrCoefficient(VectorFunctionCoefficient &vfcoeff,
Array<int> &attr)
{
int i, j, fdof, d, ind, vdim;
Vector val;
const FiniteElement *fe;
ElementTransformation *transf;
Array<int> vdofs;
this->HostReadWrite();
vdim = fes->GetVDim();
// loop over boundary elements
for (i = 0; i < fes->GetNBE(); i++)
{
// if boundary attribute is 1 (Dirichlet)
if (attr[fes->GetBdrAttribute(i) - 1])
{
fe = fes->GetBE(i);
fdof = fe->GetDof();
transf = fes->GetBdrElementTransformation(i);
const IntegrationRule &ir = fe->GetNodes();
fes->GetBdrElementVDofs(i, vdofs);
// loop over dofs
for (j = 0; j < fdof; j++)
{
const IntegrationPoint &ip = ir.IntPoint(j);
transf->SetIntPoint(&ip);
vfcoeff.Eval(val, *transf, ip);
// loop over vector dimensions
for (d = 0; d < vdim; d++)
{
ind = vdofs[fdof*d+j];
if ( (ind = vdofs[fdof*d+j]) < 0 )
{
val(d) = -val(d), ind = -1-ind;
}
(*this)(ind) = val(d);
}
}
}
}
}
void GridFunction::ProjectBdrCoefficient(VectorFunctionRestrictedCoefficient &vfcoeff)
{
int i, j, fdof, d, ind, vdim;
Vector val;
const FiniteElement *fe;
ElementTransformation *transf;
Array<int> vdofs;
const Array<int> &active_attr = vfcoeff.GetActiveAttr();
const Array2D<double> &scale = vfcoeff.GetActiveScale();
this->HostReadWrite();
vdim = fes->GetVDim();
// loop over boundary elements
for (i = 0; i < fes->GetNBE(); i++) {
// if boundary attribute is 1 (Dirichlet)
if (active_attr[fes->GetBdrAttribute(i) - 1]) {
// instantiate a BC object
fe = fes->GetBE(i);
fdof = fe->GetDof();
transf = fes->GetBdrElementTransformation(i);
const IntegrationRule &ir = fe->GetNodes();
fes->GetBdrElementVDofs(i, vdofs);
// loop over dofs
for (j = 0; j < fdof; j++) {
const IntegrationPoint &ip = ir.IntPoint(j);
transf->SetIntPoint(&ip);
vfcoeff.Eval(val, *transf, ip);
// loop over vector dimensions
for (d = 0; d < vdim; d++) {
// check if the vector component (i.e. dof) is not constrained by a
// partial essential BC
if (scale(transf->Attribute - 1, d) > 0.0) {
ind = vdofs[fdof * d + j];
if ( (ind = vdofs[fdof * d + j]) < 0) {
val(d) = -val(d), ind = -1 - ind;
}
(*this)(ind) = val(d); // placing computed value in grid function
}
}
}
}
}
}
void GridFunction::ProjectBdrCoefficientNormal(
VectorCoefficient &vcoeff, Array<int> &bdr_attr)
{
@@ -4276,7 +4372,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 +4554,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);
+15 -2
View File
@@ -463,6 +463,20 @@ public:
will not be touched. */
virtual void ProjectBdrCoefficient(Coefficient *coeff[], Array<int> &attr);
/** add ProjectBdrCoefficient() taking a VectorFunctionCoefficient and projecting
onto the grid function, srw. */
void ProjectBdrCoefficient(VectorFunctionCoefficient &vfcoeff,
Array<int> &attr);
/** add ProjectBdrCoefficient() taking a VectorFunctionRestrictedCoefficient, srw
This is the same as the function above taking a VectorFunctionCoefficient, it
is just that the function coefficient is now restricted. This function
specifically projects values associated with ExaConstits essential BC
management and input. This routine does not, at this time, project a
VectorFunctionRestrictedCoefficient onto all grid function dofs.
This function should be renamed to reflect this. */
void ProjectBdrCoefficient(VectorFunctionRestrictedCoefficient &vfcoeff);
/** Project the normal component of the given VectorCoefficient on
the boundary. Only boundary attributes that are marked in
'bdr_attr' are projected. Assumes RT-type VectorFE GridFunction. */
@@ -724,7 +738,6 @@ public:
virtual ~GridFunction() { Destroy(); }
};
/** Overload operator<< for std::ostream and GridFunction; valid also for the
derived class ParGridFunction */
std::ostream &operator<<(std::ostream &out, const GridFunction &sol);
@@ -928,7 +941,6 @@ public:
/// Overload operator<< for std::ostream and QuadratureFunction.
std::ostream &operator<<(std::ostream &out, const QuadratureFunction &qf);
double ZZErrorEstimator(BilinearFormIntegrator &blfi,
GridFunction &u,
GridFunction &flux,
@@ -1042,6 +1054,7 @@ inline void QuadratureFunction::SetSpace(QuadratureSpace *qspace_,
inline void QuadratureFunction::GetElementValues(int idx, Vector &values)
{
// element offset is the number of quadrature points for that element
const int s_offset = qspace->element_offsets[idx];
const int sl_size = qspace->element_offsets[idx+1] - s_offset;
values.NewDataAndSize(data + vdim*s_offset, vdim*sl_size);
+7 -1
View File
@@ -129,7 +129,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;
}
+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
{
+16
View File
@@ -53,6 +53,22 @@ void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
}
}
void NonlinearForm::SetEssentialBC(const Array<int> &bdr_attr_is_ess,
const Array2D<bool> &bdr_component,
Vector *rhs)
{
// virtual call, works in parallel too
fes->GetEssentialTrueDofs(bdr_attr_is_ess, ess_tdof_list, bdr_component);
if (rhs)
{
for (int i = 0; i < ess_tdof_list.Size(); i++)
{
(*rhs)(ess_tdof_list[i]) = 0.0;
}
}
}
void NonlinearForm::SetEssentialVDofs(const Array<int> &ess_vdofs_list)
{
if (!P)
+20 -1
View File
@@ -146,8 +146,26 @@ public:
have zero entries at the essential true dofs. */
void SetEssentialBC(const Array<int> &bdr_attr_is_ess, Vector *rhs = NULL);
/** @brief Specify essential boundary conditions.
For spaces with 'vdim' > 1, the 'bdr_component' array can be used
to restricts the marked tDOFs per boundary to the specified components.
If vdim > 1 then one can specify per boundary attribute which components
on a boundary are essential by assigning a value of true to its location
in the bdr_component array.
The bdr_component has dimensions number of boundary attributes x vdim
This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
the result internally for use by other methods. If the @a rhs pointer is
not NULL, its essential true dofs will be set to zero. This makes it
"compatible" with the output vectors from the Mult() method which also
have zero entries at the essential true dofs. */
void SetEssentialBC(const Array<int> &bdr_attr_is_ess,
const Array2D<bool> &bdr_component,
Vector *rhs);
/// Specify essential boundary conditions.
/** Use either SetEssentialBC() or SetEssentialTrueDofs() if possible. */
/** @deprecated Use either SetEssentialBC() or SetEssentialTrueDofs(). */
void SetEssentialVDofs(const Array<int> &ess_vdofs_list);
/// Specify essential boundary conditions.
@@ -179,6 +197,7 @@ public:
vectors, i.e. their size must be fes->GetTrueVSize(). */
virtual void Mult(const Vector &x, Vector &y) const;
/** @brief Compute the gradient Operator of the NonlinearForm corresponding
to the state @a x. */
/** Any previously specified essential boundary conditions will be
+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
}
}
+23 -1
View File
@@ -80,6 +80,15 @@ void NonlinearFormIntegrator::AssembleElementVector(
" is not overloaded!");
}
void NonlinearFormIntegrator::AssembleElementVector(
const FiniteElement &el, ElementTransformation &Ttr_beg,
ElementTransformation &Ttr_end,
const Vector &elfun, Vector &elvect, const Vector &elvel)
{
mfem_error("NonlinearFormIntegrator::AssembleElementVector"
" is not overloaded!");
}
void NonlinearFormIntegrator::AssembleFaceVector(
const FiniteElement &el1, const FiniteElement &el2,
FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
@@ -105,6 +114,20 @@ void NonlinearFormIntegrator::AssembleFaceGrad(
" is not overloaded!");
}
void NonlinearFormIntegrator::AssembleGradPA(
const FiniteElementSpace &fes)
{
mfem_error("NonlinearFormIntegrator::AssembleGradPA"
" is not overloaded!");
}
void NonlinearFormIntegrator::AssembleEA(
const FiniteElementSpace &fes, Vector &emat)
{
mfem_error("NonlinearFormIntegrator::AssembleEA"
" is not overloaded!");
}
double NonlinearFormIntegrator::GetElementEnergy(
const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
{
@@ -166,7 +189,6 @@ double BlockNonlinearFormIntegrator::GetElementEnergy(
return 0.0;
}
double InverseHarmonicModel::EvalW(const DenseMatrix &J) const
{
Z.SetSize(J.Width());
+34 -10
View File
@@ -57,6 +57,16 @@ public:
ElementTransformation &Tr,
const Vector &elfun, Vector &elvect);
//This function really only be used whenever UMATs are being used,
//since UMATs require a velocity based solution vector,
//However, they depend on a lot of variables based on displacements
//to evolve their models. So, we need to be able to take
//derivatives using the beginning and the end time step mesh nodes.
//The other models should be using the one above this
virtual void AssembleElementVector(const FiniteElement &el,
ElementTransformation &Ttr_beg,
ElementTransformation &Ttr_end,
const Vector &elfun, Vector &elvect, const Vector &elvel);
/// @brief Perform the local action of the NonlinearFormIntegrator resulting
/// from a face integral term.
virtual void AssembleFaceVector(const FiniteElement &el1,
@@ -112,6 +122,14 @@ public:
This method can be called only after the method AssemblePA() has been
called. */
virtual void AddMultPA(const Vector &x, Vector &y) const;
/** @brief Prepare the integrator for partial assembly (PA) gradient
evaluations on the given FE space @a fes */
/** The result of the partial assembly is stored internally so that it can be
used later in the methods AddMultGradPA() and AssembleGradDiagonalPA().
The state Vector @a x is an E-vector. */
virtual void AssembleGradPA(const FiniteElementSpace &fes);
virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat);
/// Method for partially assembled gradient action.
/** All arguments are E-vectors. This method can be called only after the
@@ -188,17 +206,16 @@ public:
virtual ~BlockNonlinearFormIntegrator() { }
};
/// Abstract class for hyperelastic models
class HyperelasticModel
/// Abstract class for nonlinear models
class NonlinearModel
{
protected:
ElementTransformation *Ttr; /**< Reference-element to target-element
transformation. */
public:
HyperelasticModel() : Ttr(NULL) { }
virtual ~HyperelasticModel() { }
NonlinearModel() : Ttr(NULL) { }
virtual ~NonlinearModel() { }
/// A reference-element to target-element transformation that can be used to
/// evaluate Coefficient%s.
@@ -206,11 +223,6 @@ public:
point of interest. */
void SetTransformation(ElementTransformation &Ttr_) { Ttr = &Ttr_; }
/** @brief Evaluate the strain energy density function, W = W(Jpt).
@param[in] Jpt Represents the target->physical transformation
Jacobian matrix. */
virtual double EvalW(const DenseMatrix &Jpt) const = 0;
/** @brief Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
@param[in] Jpt Represents the target->physical transformation
Jacobian matrix.
@@ -235,6 +247,18 @@ public:
};
/// Abstract class for hyperelastic models
class HyperelasticModel : public NonlinearModel
{
public:
/** @brief Evaluate the strain energy density function, W = W(Jpt).
@param[in] Jpt Represents the target->physical transformation
Jacobian matrix. */
virtual double EvalW(const DenseMatrix &Jpt) const = 0;
};
/** Inverse-harmonic hyperelastic model with a strain energy density function
given by the formula: W(J) = (1/2) det(J) Tr((J J^t)^{-1}) where J is the
deformation gradient. */
+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.");

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