Compare commits

..
95 changed files with 2437 additions and 2584 deletions
+31 -42
View File
@@ -47,36 +47,26 @@ jobs:
strategy:
matrix:
os: [ubuntu-18.04, macos-10.15]
target: [dbg, opt]
mpi: [seq, par]
target: [debug, optim]
mpi: [sequential, parallel]
build-system: [make]
hypre-target: [int32]
# 'include' allows us to:
# - Add a variable to all jobs without creating a new matrix dimension.
# Codecov is defined that way.
# - Add a new combination.
# 'build-system: cmake' and 'hypre-target: int64'
# 'include' allows us to
# - add a variable without creating a new matrix dimension.
# - add a new combination ('build-system: cmake' case here)
#
# note: we will gather coverage info for any non-debug run except the
# CMake build.
include:
- target: dbg
- target: debug
codecov: NO
- target: opt
- target: optim
codecov: YES
- os: ubuntu-18.04
target: opt
target: optim
codecov: NO
mpi: par
mpi: parallel
build-system: cmake
hypre-target: int32
- os: ubuntu-18.04
target: opt
codecov: NO
mpi: par
build-system: make
hypre-target: int64
name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}
name: ${{ matrix.os }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.build-system }}
runs-on: ${{ matrix.os }}
@@ -102,7 +92,7 @@ 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-18.04'
if: matrix.mpi == 'parallel' && matrix.os == 'ubuntu-18.04'
run: |
sudo apt-get install mpich libmpich-dev
export MAKE_CXX_FLAG="MPICXX=mpic++"
@@ -113,11 +103,11 @@ jobs:
sudo apt-get install lcov
- name: Set up Homebrew
if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-10.15'
if: ( matrix.mpi == 'parallel' || 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-10.15'
if: matrix.mpi == 'parallel' && matrix.os == 'macos-10.15'
run: |
export HOMEBREW_NO_INSTALL_CLEANUP=1
brew install openmpi
@@ -133,40 +123,39 @@ jobs:
# Install will only run on cache miss.
- name: cache hypre
id: hypre-cache
if: matrix.mpi == 'par'
if: matrix.mpi == 'parallel'
uses: actions/cache@v2
with:
path: ${{ env.HYPRE_TOP_DIR }}
key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.0
key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2
- name: get hypre
if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true'
uses: mfem/github-actions/build-hypre@v2.0
if: matrix.mpi == 'parallel' && steps.hypre-cache.outputs.cache-hit != 'true'
uses: mfem/github-actions/build-hypre@v1.0
with:
archive: ${{ env.HYPRE_ARCHIVE }}
dir: ${{ env.HYPRE_TOP_DIR }}
target: ${{ matrix.hypre-target }}
hypre-archive: ${{ env.HYPRE_ARCHIVE }}
hypre-dir: ${{ env.HYPRE_TOP_DIR }}
# Get Metis through cache, or build it.
# Install will only run on cache miss.
- name: cache metis
id: metis-cache
if: matrix.mpi == 'par'
if: matrix.mpi == 'parallel'
uses: actions/cache@v2
with:
path: ${{ env.METIS_TOP_DIR }}
key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.0
key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2
- name: install metis
if: matrix.mpi == 'par' && steps.metis-cache.outputs.cache-hit != 'true'
uses: mfem/github-actions/build-metis@v2.0
if: matrix.mpi == 'parallel' && steps.metis-cache.outputs.cache-hit != 'true'
uses: mfem/github-actions/build-metis@v1.0
with:
archive: ${{ env.METIS_ARCHIVE }}
dir: ${{ env.METIS_TOP_DIR }}
metis-archive: ${{ env.METIS_ARCHIVE }}
metis-dir: ${{ env.METIS_TOP_DIR }}
# MFEM build and test
- name: build
uses: mfem/github-actions/build-mfem@v2.0
uses: mfem/github-actions/build-mfem@v1.0
with:
os: ${{ matrix.os }}
target: ${{ matrix.target }}
@@ -179,17 +168,17 @@ jobs:
# Run checks (and only checks) on debug targets
- name: checks
if: matrix.build-system == 'make' && matrix.target == 'dbg'
if: matrix.build-system == 'make' && matrix.target == 'debug'
run: |
cd ${{ env.MFEM_TOP_DIR }} && make check
- name: unit tests
if: matrix.build-system == 'make' && matrix.target == 'opt'
if: matrix.build-system == 'make' && matrix.target == 'optim'
run: |
cd ${{ env.MFEM_TOP_DIR }} && make unittest
- name: tests
if: matrix.build-system == 'make' && matrix.target == 'opt'
if: matrix.build-system == 'make' && matrix.target == 'optim'
run: |
cd ${{ env.MFEM_TOP_DIR }} && make test
@@ -201,8 +190,8 @@ jobs:
# Code coverage (process and upload reports)
- name: codecov
if: matrix.codecov == 'YES'
uses: mfem/github-actions/upload-coverage@v2.0
uses: mfem/github-actions/upload-coverage@v1.0
with:
name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}
name: ${{ matrix.os }}-${{ matrix.mpi }}
project_dir: ${{ env.MFEM_TOP_DIR }}
directories: "fem general linalg mesh"
+9 -10
View File
@@ -53,33 +53,32 @@ jobs:
uses: actions/cache@v2
with:
path: ${{ env.HYPRE_TOP_DIR }}
key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2.0
key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2
- name: Get Hypre
if: steps.hypre-cache.outputs.cache-hit != 'true'
uses: mfem/github-actions/build-hypre@v2.0
uses: mfem/github-actions/build-hypre@master
with:
archive: ${{ env.HYPRE_ARCHIVE }}
dir: ${{ env.HYPRE_TOP_DIR }}
target: int32
hypre-archive: ${{ env.HYPRE_ARCHIVE }}
hypre-dir: ${{ env.HYPRE_TOP_DIR }}
- name: Cache Metis Install
id: metis-cache
uses: actions/cache@v2
with:
path: ${{ env.METIS_TOP_DIR }}
key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.0
key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2
- name: Install Metis
if: steps.metis-cache.outputs.cache-hit != 'true'
uses: mfem/github-actions/build-metis@v2.0
uses: mfem/github-actions/build-metis@master
with:
archive: ${{ env.METIS_ARCHIVE }}
dir: ${{ env.METIS_TOP_DIR }}
metis-archive: ${{ env.METIS_ARCHIVE }}
metis-dir: ${{ env.METIS_TOP_DIR }}
# MFEM build and test
- name: build-mfem
uses: mfem/github-actions/build-mfem@v2.0
uses: mfem/github-actions/build-mfem@master
with:
os: ${{ runner.os }}
target: optim
+33 -6
View File
@@ -28,24 +28,49 @@ jobs:
access_token: ${{ github.token }}
- name: checkout mfem
uses: actions/checkout@v2
with:
path: mfem
- name: copyright check
id: copyright
run: |
./config/githooks/pre-push --copyright
cd mfem
if git grep -l "^\(#\|//\).*\(\-2020\|\ 2010,\)" > matches.txt
then
echo "Please update the following files to Copyright (c) 2010-2021:"
cat matches.txt
exit 1
else
echo "No outdated copyright found."
fi
continue-on-error: true
- name: license check
id: license
run: |
./config/githooks/pre-push --license
cd mfem
if git grep -li "^\(#\|//\).*GNU\ Lesser\ General\ Public\ License" > matches.txt
then
echo "Please update the following files to the BSD-3 license:"
cat matches.txt
exit 1
else
echo "No GNU GPL license found."
fi
continue-on-error: true
- name: release check
id: release
run: |
./config/githooks/pre-push --release
cd mfem
if git grep -l "^\(#\|//\).*LLNL\-CODE\-443211" > matches.txt
then
echo "Please update the following files to LLNL-CODE-806117:"
cat matches.txt
exit 1
else
echo "No outdated release number found."
fi
continue-on-error: true
- name: wrap-up
@@ -75,7 +100,8 @@ jobs:
- name: style check
run: |
./config/githooks/pre-push --style
cd tests/scripts
./runtest code-style
documentation:
runs-on: ubuntu-18.04
@@ -107,4 +133,5 @@ jobs:
run: |
git fetch origin master:master
git checkout -b gh-actions-branch-history
./config/githooks/pre-push --history
cd tests/scripts
./runtest branch-history
-1
View File
@@ -26,7 +26,6 @@ CMakeFiles/
config/_config.hpp
config/config.mk
config/sample-runs-build.log
config/user.mk
doc/CodeDocumentation.conf
doc/CodeDocumentation.html
doc/CodeDocumentation
+8 -27
View File
@@ -48,58 +48,39 @@ variables:
TPLS_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/tpls.git
TESTS_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/tests.git
AUTOTEST_REPO: ssh://git@mybitbucket.llnl.gov:7999/mfem/autotest.git
MFEM_DATA_REPO: https://github.com/mfem/data.git
ARTIFACTS_DIR: artifacts
# The pipeline is divided into stages. Usually, these are also synchronization
# points, however, we use "needs" keyword to express the DAG of jobs for more
# efficiency.
# - We use setup and setup_baseline phases to download content outside of mfem
# directory.
# - We use setup phase to download content outside of mfem directory.
# - Allocate/Release is where quartz resources are allocated/released once for all.
# - Build and Test is where we build and MFEM for multiple toolchains.
# - Baseline_checks gathers baseline-type test suites execution
# - Baseline_publish, only available on master, allows to update baseline
# results
stages:
- setup
- q_allocate_resources
- q_build_and_test
- q_release_resources
- l_build_and_test
- c_build_and_test
- setup_baseline
- setup
- baseline_check
- baseline_to_autotest
- baseline_publish
# setup clones the mfem/data repo in ${BUILD_ROOT}. The build_and_test script
# then symlinks the repo to the parent directory of the MFEM source directory.
# Unit tests that depend on the mfem/data repo will then detect that this
# directory is present and be enabled.
# The setup job in setup stage don't rely on MFEM git repo. It prepares a
# pipeline-wide working directory downloading/updating external repos.
# TODO: updating tests and tpls is not necessary anymore since pipelines are
# now using unique directories so repo are never shared with another pipeline.
# This is not memory efficient (we keep a lot of data), hence this reminder.
# Setup
setup:
tags:
- shell
- quartz
stage: setup
variables:
GIT_STRATEGY: none
script:
- mkdir -p ${BUILD_ROOT} && cd ${BUILD_ROOT}
- if [ ! -d data ]; then git clone ${MFEM_DATA_REPO}; fi
needs: []
# The setup_baseline job in setup stage_baseline doesn't rely on MFEM git repo.
# It prepares a pipeline-wide working directory downloading/updating external
# repos. TODO: updating tests and tpls is not necessary anymore since pipelines
# are now using unique directories so repo are never shared with another
# pipeline. This is not memory efficient (we keep a lot of data), hence this
# reminder.
setup_baseline:
tags:
- shell
- quartz
stage: setup_baseline
variables:
GIT_STRATEGY: none
script:
+1 -1
View File
@@ -25,7 +25,7 @@
.build_and_test_on_lassen:
extends: [.build_blueos_3_ppc64le_ib_script, .on_lassen]
stage: l_build_and_test
needs: [setup]
needs: []
opt_mpi_cuda_xl_16_1_1_8:
variables:
+1 -2
View File
@@ -94,7 +94,6 @@ q_report_failure:
.build_and_test_on_quartz:
extends: [.build_toss_3_x86_64_ib_script, .on_quartz]
stage: q_build_and_test
needs: [setup]
# Build MFEM
debug_ser_gcc_4_9_3:
@@ -140,7 +139,7 @@ opt_par_gcc_6_1_0_pumi:
# Baseline
baselinecheck_mfem_intel_quartz:
extends: [.baselinecheck_mfem, .on_quartz]
needs: [setup_baseline]
needs: [setup]
update_autotest:
extends: [.on_quartz]
+469
View File
@@ -0,0 +1,469 @@
# Copyright (c) 2010-2021, 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.
language: cpp
os: linux
dist: bionic
stages:
- checks
- tests
- optional
env:
global:
- HYPRE_ARCHIVE=v2.19.0.tar.gz
HYPRE_URL=https://github.com/hypre-space/hypre/archive/$HYPRE_ARCHIVE
HYPRE_TOP_DIR=hypre-2.19.0
jobs:
include:
# ========================
# Checks
# ========================
# - code-style
# - documentation
# - gitignore
- stage: checks
os: linux
dist: xenial
name: "code-style"
addons:
apt:
packages:
- astyle=2.05.1-0ubuntu1
script:
- cd ${TRAVIS_BUILD_DIR}
- cd tests/scripts
- ./runtest code-style
- stage: checks
os: linux
name: "documentation"
addons:
apt:
packages:
- doxygen
- graphviz
script:
- cd ${TRAVIS_BUILD_DIR}
- cd tests/scripts
- ./runtest documentation
- stage: checks
os: linux
name: "gitignore"
addons:
apt:
packages:
- mpich
- libmpich-dev
env: MPI=YES
before_script:
- cd ${TRAVIS_BUILD_DIR}
- mpicxx -v
- make config MFEM_USE_MPI=YES MFEM_MPI_NP=2
- make all -j3
- make test-noclean
script:
- cd tests/scripts
- ./runtest gitignore
cache:
ccache: true
directories:
- $TRAVIS_BUILD_DIR/../$HYPRE_TOP_DIR/src/hypre
- $TRAVIS_BUILD_DIR/../metis-4.0
before_cache:
- cd $TRAVIS_BUILD_DIR/../metis-4.0;
mv libmetis.a Lib ..; rm -rf * ; mv ../libmetis.a ../Lib .;
rm -f Lib/*.{c,o}
# ========================
# Optional Checks/Tests
# ========================
# - branch-history
- stage: optional
name: "branch-history"
if: branch != next
# need full git history for the binary/big files check
git:
depth: false
script:
- cd ${TRAVIS_BUILD_DIR}
# update master
- git fetch origin master:master
# checkout a branch (otherwise Travis works in detached head)
- git checkout -b travis_tests
- cd tests/scripts
- ./runtest branch-history
# ========================
# Linux tests
# ========================
# - serial + debug
# - serial
# - parallel + debug
# - parallel
- stage: tests
os: linux
compiler: gcc
name: "Linux: Serial + Debug"
env: DEBUG=YES
MPI=NO
CODECOV=NO
MFEM_TEST_TARGET=check
cache:
ccache: true
- os: linux
compiler: gcc
name: "Linux: Serial"
env: DEBUG=NO
MPI=NO
CODECOV=NO
MFEM_TEST_TARGET=test
cache:
ccache: true
- os: linux
compiler: gcc
name: "Linux: Parallel + Debug"
addons:
apt:
# sources:
# - ubuntu-toolchain-r-test
packages:
# GCC 4.9
# - g++-4.9
# MPICH
- mpich
- libmpich-dev
# OpenMPI
# - openmpi-bin
# - libopenmpi-dev
env: DEBUG=YES
MPI=YES
CODECOV=NO
MFEM_TEST_TARGET=check
NPROCS=2
cache:
ccache: true
directories:
- $TRAVIS_BUILD_DIR/../$HYPRE_TOP_DIR/src/hypre
- $TRAVIS_BUILD_DIR/../metis-4.0
before_cache:
- cd $TRAVIS_BUILD_DIR/../metis-4.0;
mv libmetis.a Lib ..; rm -rf * ; mv ../libmetis.a ../Lib .;
rm -f Lib/*.{c,o}
- os: linux
compiler: gcc
name: "Linux: Parallel"
addons:
apt:
# sources:
# - ubuntu-toolchain-r-test
packages:
# GCC 4.9
# - g++-4.9
# MPICH
- mpich
- libmpich-dev
# OpenMPI
# - openmpi-bin
# - libopenmpi-dev
env: DEBUG=NO
MPI=YES
CODECOV=YES
MFEM_TEST_TARGET=test
NPROCS=2
cache:
ccache: true
directories:
- $TRAVIS_BUILD_DIR/../$HYPRE_TOP_DIR/src/hypre
- $TRAVIS_BUILD_DIR/../metis-4.0
before_cache:
- cd $TRAVIS_BUILD_DIR/../metis-4.0;
mv libmetis.a Lib ..; rm -rf * ; mv ../libmetis.a ../Lib .;
rm -f Lib/*.{c,o}
- os: linux
compiler: gcc
name: "Linux: Parallel (cmake)"
addons:
apt:
packages:
- mpich
- libmpich-dev
env: MPI=YES
NPROCS=2
script:
- cd ${TRAVIS_BUILD_DIR}
- mkdir ${TRAVIS_BUILD_DIR}/build
- cd ${TRAVIS_BUILD_DIR}/build
- cmake ..
-DMFEM_USE_MPI=ON
-DHYPRE_DIR=${TRAVIS_BUILD_DIR}/../$HYPRE_TOP_DIR/src/hypre
-DMFEM_MPI_NP=$NPROCS
- make -j3 mfem examples
- cd ${TRAVIS_BUILD_DIR}/build/tests/unit
- make -j3
- ctest --output-on-failure
cache:
ccache: true
directories:
- $TRAVIS_BUILD_DIR/../$HYPRE_TOP_DIR/src/hypre
- $TRAVIS_BUILD_DIR/../metis-4.0
before_cache:
- cd $TRAVIS_BUILD_DIR/../metis-4.0;
mv libmetis.a Lib ..; rm -rf * ; mv ../libmetis.a ../Lib .;
rm -f Lib/*.{c,o}
# ========================
# Mac OS X tests
# ========================
# - serial + debug
# - serial
# - parallel + debug
# - parallel
- os: osx
osx_image: xcode11.2
compiler: clang
name: "Mac: Serial + Debug"
addons:
homebrew:
packages:
- ccache
env: DEBUG=YES
MPI=NO
CODECOV=NO
MFEM_TEST_TARGET=check
cache:
ccache: true
- os: osx
osx_image: xcode11.2
compiler: clang
name: "Mac: Serial"
addons:
homebrew:
packages:
- ccache
env: DEBUG=NO
MPI=NO
CODECOV=NO
MFEM_TEST_TARGET=test
cache:
ccache: true
- os: osx
osx_image: xcode11.2
compiler: clang
name: "Mac: Parallel + Debug"
addons:
homebrew:
packages:
- ccache
env: DEBUG=YES
MPI=YES
CODECOV=NO
MFEM_TEST_TARGET=check
NPROCS=4
TMPDIR=/tmp
cache:
ccache: true
directories:
- $TRAVIS_BUILD_DIR/../$HYPRE_TOP_DIR/src/hypre
- $TRAVIS_BUILD_DIR/../metis-4.0
- $HOME/local-cached
before_cache:
- cd $TRAVIS_BUILD_DIR/../metis-4.0;
mv libmetis.a Lib ..; rm -rf * ; mv ../libmetis.a ../Lib .;
rm -f Lib/*.{c,o}
- os: osx
osx_image: xcode11.2
compiler: clang
name: "Mac: Parallel"
addons:
homebrew:
packages:
- ccache
env: DEBUG=NO
MPI=YES
CODECOV=YES
MFEM_TEST_TARGET=test
NPROCS=4
TMPDIR=/tmp
cache:
ccache: true
directories:
- $TRAVIS_BUILD_DIR/../$HYPRE_TOP_DIR/src/hypre
- $TRAVIS_BUILD_DIR/../metis-4.0
- $HOME/local-cached
before_cache:
- cd $TRAVIS_BUILD_DIR/../metis-4.0;
mv libmetis.a Lib ..; rm -rf * ; mv ../libmetis.a ../Lib .;
rm -f Lib/*.{c,o}
before_install:
# No addon for brew yet, have to install OSX packages this way.
# - if [ $TRAVIS_OS_NAME == "osx" ] && [ $MPI == "YES" ]; then
# brew install open-mpi;
# fi
# Disable ccache while building dependencies that are cached:
- echo "before \$PATH = $PATH";
export PATH=${PATH//\/usr\/lib\/ccache:/};
echo "after \$PATH = $PATH"
# On Mac OS X, build and cache OpenMPI 2.1.6:
- if [ $TRAVIS_OS_NAME == "osx" ] && [ $MPI == "YES" ]; then
if [ ! -e $HOME/local-cached/bin/mpicc ]; then
mkdir -p $HOME/builds && cd $HOME/builds &&
wget https://download.open-mpi.org/release/open-mpi/v2.1/openmpi-2.1.6.tar.bz2 &&
tar jxf openmpi-2.1.6.tar.bz2 &&
mkdir openmpi-build && cd openmpi-build &&
../openmpi-2.1.6/configure --prefix=$HOME/local-cached &&
make -j3 all && make install;
fi;
PATH=$HOME/local-cached/bin:$PATH;
cd $TRAVIS_BUILD_DIR;
fi
# Update environment to find g++ 4.9 installation first.
# - if [ $TRAVIS_OS_NAME == "linux" ]; then
# mkdir -p latest-gcc-symlinks;
# ln -s /usr/bin/g++-4.9 latest-gcc-symlinks/g++;
# ln -s /usr/bin/gcc-4.9 latest-gcc-symlinks/gcc;
# ln -s /usr/bin/gcov-4.9 latest-gcc-symlinks/gcov;
# export PATH=$PWD/latest-gcc-symlinks:$PATH;
# fi
# Install tool to upload code coverage reports to coveralls.io
- if [ "$CODECOV" == "YES" ]; then
export PYTHONUSERBASE=$HOME/local;
pip install --user cpp-coveralls;
pip install --user pyyaml;
PATH=$HOME/local/bin:$PATH;
fi
install:
# Set MPI compilers, print compiler version
- if [ $MPI == "YES" ]; then
if [ "$TRAVIS_OS_NAME" == "linux" ]; then
export MPICH_CC="$CC";
export MPICH_CXX="$CXX";
else
export OMPI_CC="$CC";
export OMPI_CXX="$CXX";
mpic++ --showme:version;
fi;
mpic++ -v;
else
$CXX -v;
fi
# Back out of the mfem directory to install the libraries
- cd ..
# hypre
- if [ $MPI == "YES" ]; then
if [ ! -e $HYPRE_TOP_DIR/src/hypre/lib/libHYPRE.a ]; then
wget $HYPRE_URL;
rm -rf $HYPRE_TOP_DIR;
tar xvzf $HYPRE_ARCHIVE;
cd $HYPRE_TOP_DIR/src;
./configure --disable-fortran CC=mpicc CXX=mpic++;
make -j3;
cd ../..;
else
echo "Reusing cached $HYPRE_TOP_DIR/";
fi;
ln -s $HYPRE_TOP_DIR hypre;
else
echo "Serial build, not using hypre";
fi
# METIS, use a mirror because the original source server is not always up.
# Original url:
# http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/OLD/metis-4.0.3.tar.gz
- if [ $MPI == "YES" ]; then
if [ ! -e metis-4.0/libmetis.a ]; then
wget https://mfem.github.io/tpls/metis-4.0.3.tar.gz;
tar xvzf metis-4.0.3.tar.gz;
make -j3 -C metis-4.0.3/Lib CC="$CC" OPTFLAGS="-O2";
rm -rf metis-4.0;
mv metis-4.0.3 metis-4.0;
else
echo "Reusing cached metis-4.0/";
fi;
fi
# Re-enable ccache on linux; enable ccache on mac os:
- if [ $TRAVIS_OS_NAME == "linux" ]; then
export PATH="/usr/lib/ccache:$PATH";
else
if [ $TRAVIS_OS_NAME == "osx" ]; then
export PATH="/usr/local/opt/ccache/libexec:$PATH";
fi;
fi
- printf "which \$CC = "; which $CC;
printf "which \$CXX = "; which $CXX
script:
# Compiler
- if [ $MPI == "YES" ]; then
export MYCXX=mpic++;
export MAKE_CXX_FLAG=MPICXX=$MYCXX;
else
export MYCXX="$CXX";
export MAKE_CXX_FLAG=CXX=$MYCXX;
fi
# Print the compiler version
- $MYCXX -v
# Set some variables
- cd $TRAVIS_BUILD_DIR;
CPPFLAGS="";
SKIP_TEST_DIRS="";
if [ "$CODECOV" == "YES" ]; then
CPPFLAGS="--coverage -g";
fi;
if [ "$TRAVIS_OS_NAME" != "linux" ] || [ "$DEBUG" == "YES" ]; then
CPPFLAGS+=" -pedantic -Wall -Werror";
fi
# Configure the library
- make config MFEM_USE_MPI=$MPI MFEM_DEBUG=$DEBUG $MAKE_CXX_FLAG
MFEM_MPI_NP=$NPROCS CPPFLAGS="$CPPFLAGS"
# Show the configuration
- make info
# Build the library
- make -j3
# Build the examples and the miniapps
- make -j3 all
# Run tests
- make $MFEM_TEST_TARGET SKIP_TEST_DIRS="$SKIP_TEST_DIRS"
after_success:
- if [ "$CODECOV" == "YES" ]; then
coveralls --include fem --include general --include linalg --include
mesh --exclude /usr --gcov-options '\-lp' --root $TRAVIS_BUILD_DIR;
fi
+1 -14
View File
@@ -173,11 +173,7 @@ Version 4.2.1 (development)
mixed meshes. The LOR Transfer miniapp (miniapps/tools/lor-transfer.cpp) now
supports meshes with any element geometry.
- Testing improvements:
* Transitioned from Travis to GitHub Action for testing/CI on GitHub.
* Effectively remove Travis from CI.
* Use Spack (and Uberenv) to automate TPL building in LLNL GitLab tests.
* Added a set of suggested git hooks for developers in config/githooks.
- Gitlab CI: use Spack (and Uberenv) to automate the build of TPLs.
- Added new miniapps demonstrating: 1) the use of GSLIB for overlapping grids,
see gslib/schwarz_ex1, and 2) coupling different physics in different domains,
@@ -235,15 +231,6 @@ Version 4.2.1 (development)
- Added makefile rule to generate TAGS table for vi or Emacs users.
API changes
-----------
- Added an abstract interface `mfem::FaceRestriction` for `H1FaceRestriction`
and `L2FaceRestriction`.
In order to conform with the semantic of `MultTranspose` in `mfem::Operator`,
`mfem::FaceRestriction::MultTranspose` now sets instead of adding values, and
`mfem::FaceRestriction::AddMultTranspose` should replace previous calls to
`mfem::FaceRestriction::MultTranspose`.
libCEED integration improvements
--------------------------------
- Refactor the libCEED integration
+16 -33
View File
@@ -4,9 +4,7 @@
<p align="center">
<a href="https://github.com/mfem/mfem/blob/master/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-brightgreen.svg"></a>
<a href="https://github.com/mfem/mfem/actions?query=workflow%3Arepo-check+branch%3Amaster"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=master"></a>
<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuild-analysis+branch%3Amaster"><img alt="Build Analysis" src="https://github.com/mfem/mfem/actions/workflows/mfem-analysis.yml/badge.svg?branch=master"></a>
<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuilds-and-tests+branch%3Amaster"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=master"></a>
<a href="https://travis-ci.org/mfem/mfem"><img alt="Build Status" src="https://travis-ci.org/mfem/mfem.svg?branch=master"></a>
<a href="https://ci.appveyor.com/project/mfem/mfem"><img alt="Build Status" src="https://ci.appveyor.com/api/projects/status/19non9sqm6msi2wy?svg=true"></a>
<a href="https://mfem.github.io/doxygen/html/index.html"><img alt="Doxygen" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></a>
</p>
@@ -65,8 +63,6 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
development branches off `mfem:master`.
- Please follow the [developer guidelines](#developer-guidelines), in particular
with regards to documentation and code styling.
- Please do not commit large/binary files to the central repository (use a fork
instead).
- Pull requests should be issued toward `mfem:master`. Make sure
to check the items off the [Pull Request Checklist](#pull-request-checklist).
- When your contribution is fully working and ready to be reviewed, add
@@ -75,7 +71,6 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
reviewers to evaluate the changes.
- The reviewers have 3 weeks to evaluate the PR and work with the author to
fix issues and implement improvements.
- During review there should be no force pushes/rewriting history in the branch.
- After approval, MFEM developers merge the PR manually in the [mfem:next branch](#masternext-workflow).
- After a week of testing in `mfem:next`, the original PR is merged in `mfem:master`.
- We use [milestones](https://github.com/mfem/mfem/milestones) to coordinate the
@@ -96,9 +91,8 @@ The MFEM source code has the following structure:
```
.
├── config
── cmake
└── ...
│ └── githooks
── cmake
│ └── ...
├── data
├── doc
├── examples
@@ -369,10 +363,6 @@ Before you can start, you need a GitHub account, here are a few suggestions:
two reviewers to evaluate the changes. The reviewers have 3 weeks to evaluate
the PR and work with the author to implement improvements and fix issues.
- Once the `ready-for-review` label has been applied and reviewers have been
assigned, the PR is considered under review. To help with the review process
there should be no force pushes/rewriting history in the branch.
- After approval, the PR is [tested](#masternext-workflow) for a week with
other approved PRs in the `mfem:next` branch.
@@ -380,20 +370,16 @@ Before you can start, you need a GitHub account, here are a few suggestions:
`mfem:next`, see the [README](tests/scripts/README) file in that directory
for more details.
- Track the GitHub Actions and Appveyor [continuous integration](#automated-testing)
- Track the Travis CI, Github Actions and Appveyor [continuous integration](#automated-testing)
builds at the end of the PR. These should generally run clean, so address any
errors as soon as possible. Please ask if you are unsure how to do that.
- Note that some tests, such as the `branch-history` check in GitHub Actions
are safeguards that are allowed to fail in certain cases.
- Note that some tests, such as the `branch-history` check in Travis and Github
Actions are safeguards that are allowed to fail in certain cases.
- Other tests, such as the `code-style`, `documentation` and `gitignore`
checks in GitHub Actions enforce MFEM-specific rules which are explained in
the error messages and the `tests/scripts` directory.
- Also note that the tests `branch-history` and `repos-checks` found in GitHub
Actions can be triggered automatically before each push using git hooks. See
the [git hooks README](config/githooks/README.md) for a detailed explanation.
checks in Travis and Github Actions enforce MFEM-specific rules which are
explained in the error messages and the `tests/scripts` directory.
- If triggered, track the status of the LLNL GitLab tests. If failing, ask
one of the _LLNL developers_ for details.
@@ -413,7 +399,7 @@ Before a PR can be merged, it should satisfy the following:
- [ ] Does `make` or `cmake` have a new target?
- [ ] Did the requirements or the installation process change? *(rare)*
- [ ] Update continuous integration server configurations if necessary (e.g. with new version requirements for each of MFEM's dependencies)
- [ ] `.github`
- [ ] `.travis.yml`
- [ ] `.appveyor.yml`
- [ ] Update `.gitignore`:
- [ ] Check if `make distclean; git status` shows any files that were generated from the source by the project (not an IDE) but we don't want to track in the repository.
@@ -530,7 +516,7 @@ MFEM uses a `master`/`next`-branch workflow as described below:
- [ ] `doc/CodeDocumentation.conf.in`
- [ ] Check that version requirements for each of MFEM's dependencies are documented in `INSTALL` and up-to-date
- [ ] Check that continuous integration server configurations reflect the dependency version requirements of the new release
- [ ] `.github`
- [ ] `.travis.yml`
- [ ] `.appveyor.yml`
- [ ] Update the `CHANGELOG` to organize all release contributions
- [ ] Review the whole source code once over
@@ -592,17 +578,14 @@ MFEM uses a `master`/`next`-branch workflow as described below:
MFEM has several levels of automated testing running on GitHub, as well as on
local Mac and Linux workstations, and Livermore Computing clusters at LLNL.
In addition, developers can set local git hooks to run some quick checks on
commit or push, see the [README](config/githooks/README.md) in the `config/githooks`
directory.
### Linux and Mac smoke tests
We use GitHub Actions to drive the default tests on the `master` and `next`
branches. See the `.github/workflows` files and the logs at
[https://github.com/mfem/mfem/actions](https://github.com/mfem/mfem/actions).
We use Travis CI and Github Actions to drive the default tests on the `master`
and `next` branches. See the `.travis` file and the logs at
[https://travis-ci.org/mfem/mfem](https://travis-ci.org/mfem/mfem).
Testing using GitHub Actions should be kept lightweight, as there is a time
constraint on jobs. Two virtual machines are configured - Mac (OS X) and Linux.
Testing using Travis CI and Github Actions should be kept lightweight, as there
is a time constraint on jobs. Two virtual machines are configured - Mac (OS X)
and Linux.
- Tests on the `master` branch are triggered whenever a PR is issued on this branch.
- Tests on the `next` branch are currently scheduled to run each night.
-41
View File
@@ -1,41 +0,0 @@
Finite Element Discretization Library
__
_ __ ___ / _| ___ _ __ ___
| '_ ` _ \ | |_ / _ \| '_ ` _ \
| | | | | || _|| __/| | | | | |
|_| |_| |_||_| \___||_| |_| |_|
https://mfem.org
This directory contains recommended git hooks, which are scripts that can be
used to improve your development experience with MFEM:
### The hooks
* `pre-commit` is a hook that will be applied before each commit and run
`astyle` on the code. This will ensure that your changes comply with the MFEM
code styling guidelines.
* `pre-push` is a hook that will be applied before each push to run a quick set
of tests that verify that your files headers are in compliance, and that you did
not add any large files to the repo.
### Setup
To setup the git hooks, run `make hooks`, which creates symlinks to the hooks in
the `.git/hooks` directory. Individual hooks can be enabled by manually creating
symlinks.
(You may also copy the scripts directly and customize them further, but this way
you may miss additional updates in the future.)
### Failures
The `branch-history` check can fail in some cases when the history is OK. For
example, when a large number of files were modified for a legitimate reason, or
when a picture was added for documentation.
If that is the case, make sure the failure is indeed justified, and rerun the
push command with the `--no-verify` option. This will skip the hooks, allowing
you to push those changes.
-4
View File
@@ -1,4 +0,0 @@
#!/bin/sh
# Apply automated code formatting
make -C $(git rev-parse --show-toplevel) style
-107
View File
@@ -1,107 +0,0 @@
#!/bin/bash
# Copyright (c) 2010-2021, 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.
option=${1:-""}
if [[ "${option}" == "--help" ]]; then
echo "This script runs checks on the repository."
echo "It has 2 modes: with and without an option."
echo ""
echo "Options are used in GitHub Actions and can be:"
echo " --copyright"
echo " --license"
echo " --release"
echo " --style"
echo " --history"
echo ""
echo "As a githook, the script is used without options."
echo "In that case, it will run all the checks except style."
echo ""
echo "Use --help to print this help message."
fi
cd $(git rev-parse --show-toplevel)
# copyright check
copyright=true
if [[ "${option}" == "--copyright" || "${option}" == "" ]]; then
if git grep -l "^\(#\|//\).*\(\-2020\|\ 2010,\)" > matches.txt; then
echo "Please update the following files to Copyright (c) 2010-2021:"
cat matches.txt
copyright=false
fi
fi
# license check
license=true
if [[ "${option}" == "--license" || "${option}" == "" ]]; then
if git grep -li "^\(#\|//\).*GNU\ Lesser\ General\ Public\ License" > matches.txt; then
echo "Please update the following files to the BSD-3 license:"
cat matches.txt
license=false
fi
fi
# release check
release=true
if [[ "${option}" == "--release" || "${option}" == "" ]]; then
if git grep -l "^\(#\|//\).*LLNL\-CODE\-443211" > matches.txt
then
echo "Please update the following files to LLNL-CODE-806117:"
cat matches.txt
release=false
fi
fi
# wrap-up
code=0
if ! $copyright ; then
echo "copyright check failed, unroll log for details"
code=1
fi
if ! $license ; then
echo "license check failed, unroll log for details"
code=1
fi
if ! $release ; then
echo "release check failed, unroll log for details"
code=1
fi
# `code-style` is not just a check, it will actually reformat the code if
# necessary. This means that if one pushes while the repo is in dirty state
# (changes not staged), those changes may be mixed with format changes.
# To activate this, you will need to hard-copy this hook script in the hook
# directory and uncomment only then. (See README.md)
#
## style check
#if [[ "${option}" == "--style" || "${option}" == "" ]]; then
if [[ "${option}" == "--style" ]]; then
if which astyle && [[ "$(astyle --version)" == "Artistic Style Version 2.05.1" ]]; then
cd tests/scripts
if ! ./runtest code-style; then code=1; fi
cd -
else
echo "Warning: astyle not found or version is not 2.05.1"
fi
fi
# branch-history
if [[ "${option}" == "--history" || "${option}" == "" ]]; then
git fetch origin master:master
cd tests/scripts
if ! ./runtest branch-history; then code=1; fi
cd -
fi
exit $code
+61
View File
@@ -0,0 +1,61 @@
# Copyright (c) 2010-2021, 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.
# Use the MFEM build directory
MFEM_DIR ?= ../..
MFEM_BUILD_DIR ?= ../..
SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/MeshPart/,)
CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
MFEM_LIB_FILE = mfem_is_not_built
-include $(CONFIG_MK)
SEQ_EXAMPLES = test_mesh_partition
PAR_EXAMPLES =
ifeq ($(MFEM_USE_MPI),NO)
EXAMPLES = $(SEQ_EXAMPLES)
else
EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
endif
.SUFFIXES:
.SUFFIXES: .o .cpp .mk
.PHONY: all clean
.PRECIOUS: %.o
COMMON_O= mesh_partition.o
# Remove built-in rules
%: %.cpp
%.o: %.cpp
all: $(EXAMPLES)
# Rules for building the EXAMPLES
%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
# Rules for compiling miniapp dependencies
$(COMMON_O) $($(EXAMPLES)): \
%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
# Generate an error message if the MFEM library is not built and exit
$(MFEM_LIB_FILE):
$(error The MFEM library is not built)
clean:
rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
rm -rf *.dSYM *.TVD.*breakpoints
rm output/*
+301
View File
@@ -0,0 +1,301 @@
#include "mesh_partition.hpp"
Subdomain::Subdomain(const Mesh & mesh0_)
: mesh0(&mesh0_), dim(mesh0->Dimension()), sdim(mesh0->SpaceDimension())
{
// MFEM_VERIFY(dim == 3, "Only 3D domains for now are supported");
if(mesh0->NURBSext)
{
MFEM_ABORT("Nurbs meshes are not supported yet");
}
}
void Subdomain::BuildSubMesh(const Array<int> & elems, const entity_type & etype)
{
// if mesh nodes are defined we use them for the vertices
// otherwise we use the vertices them selfs
cout << "entity type " << etype << endl;
int nv = mesh0->GetNV();
int subelems = elems.Size();
Array<int> vmarker(nv); vmarker = 0;
int numvertices = 0;
for (int ie = 0; ie<subelems; ie++)
{
int el = elems[ie];
Array<int> vertices;
switch (etype)
{
case 0: mesh0->GetElementVertices(el,vertices); break;
case 1: mesh0->GetFaceVertices(el,vertices); break;
default:
MFEM_ABORT("Wrong entity type choice");
break;
}
for (int iv=0; iv<vertices.Size(); iv++)
{
int v = vertices[iv];
if (vmarker[v]) continue;
vmarker[v] = 1;
numvertices++;
}
}
cout << "Num of new vertices: " << numvertices << endl;
// Construct new mesh
Mesh * meshptr = nullptr;
switch (etype)
{
case 0:
mesh = new Mesh(dim,numvertices, subelems);
element_map = elems;
meshptr = mesh;
break;
case 1:
surface_mesh = new Mesh(dim-1,numvertices, subelems,0,sdim);
surface_element_map = elems;
meshptr = surface_mesh;
break;
default:
MFEM_ABORT("Wrong entity type choice");
break;
}
Vector values;
const GridFunction * nodes0 = mesh0->GetNodes();
int vk = 0;
// if (nodes0) // this is NOT NEEDED here
// {
// vcoords.SetSize(sdim, mesh0->GetNV());
// for (int i = 0; i< sdim; i++)
// {
// nodes0->GetNodalValues(values,i+1);
// vcoords.SetRow(i,values);
// cout << "values size = " << values.Size() << endl;
// }
// for (int iv = 0; iv<mesh0->GetNV(); ++iv)
// {
// if (!vmarker[iv]) continue;
// meshptr->AddVertex(vcoords.GetColumn(iv));
// vmarker[iv] = ++vk;
// }
// }
// else
{
for (int iv = 0; iv<mesh0->GetNV(); ++iv)
{
if (!vmarker[iv]) continue;
meshptr->AddVertex(mesh0->GetVertex(iv));
vmarker[iv] = ++vk;
}
}
// Add elements
for (int ie = 0; ie<subelems; ie++)
{
const Element * el = nullptr;
switch (etype)
{
case 0: el = mesh0->GetElement(elems[ie]); break;
case 1: el = mesh0->GetFace(elems[ie]); break;
default: MFEM_ABORT("Wrong entity type choice"); break;
}
Element * nel = meshptr->NewElement(el->GetGeometryType());
int nv0 = el->GetNVertices();
const int * v0 = el->GetVertices();
Array<int> v1(nv0);
for (int i=0; i<nv0; i++)
{
v1[i] = vmarker[v0[i]]-1;
}
nel->SetVertices(v1.GetData());
meshptr->AddElement(nel);
}
meshptr->FinalizeTopology();
if (nodes0)
{
cout << "nodes not null" << endl;
// Extract Nodes GridFunction and determine its type
const FiniteElementSpace * fes0 = nodes0->FESpace();
Ordering::Type ordering = fes0->GetOrdering();
int order = fes0->FEColl()->GetOrder();
bool discont = fes0->IsDGSpace();
cout << "discont = " << discont << endl;
// Set curvature of the same type as original mesh
meshptr->SetCurvature(order, discont, sdim, ordering);
const FiniteElementSpace * fes1 = meshptr->GetNodalFESpace();
GridFunction * nodes = meshptr->GetNodes();
Array<int> vdofs0;
Array<int> vdofs;
Vector loc_vec;
// Copy nodes to submesh
for (int e = 0; e < elems.Size(); e++)
{
fes1->GetElementVDofs(e, vdofs);
switch (etype)
{
case 0:
fes0->GetElementVDofs(elems[e], vdofs0);
nodes0->GetSubVector(vdofs0, loc_vec);
break;
case 1:
if (!discont)
{
fes0->GetFaceVDofs(elems[e], vdofs0);
nodes0->GetSubVector(vdofs0, loc_vec);
}
else
{
const FiniteElement * el = fes1->GetFE(e);
const IntegrationRule & ir = el->GetNodes();
int np = ir.GetNPoints();
FaceElementTransformations * Tr =
const_cast<Mesh *>(mesh0)->GetFaceElementTransformations(elems[e]);
int el1 = Tr->Elem1No;
loc_vec.SetSize(vdofs.Size());
for (int i = 0; i<np; i++)
{
Tr->SetAllIntPoints(&ir[i]);
const IntegrationPoint & ip = Tr->GetElement1IntPoint();
Vector val;
nodes0->GetVectorValue(el1,ip,val);
for (int j = 0; j<val.Size(); j++)
{
loc_vec[i+j*np] = val[j];
}
}
}
break;
default:
MFEM_ABORT("Wrong entity type choice");
break;
}
nodes->SetSubVector(vdofs, loc_vec);
}
}
meshptr->Finalize();
}
void Subdomain::BuildDofMap(const entity_type & etype)
{
Array<int> elems;
FiniteElementSpace * fesptr = nullptr;
const FiniteElementCollection *fec = fes0->FEColl();
switch(etype)
{
case 0:
fesptr = new FiniteElementSpace(mesh,fec);
elems = element_map;
break;
case 1:
fesptr = new FiniteElementSpace(surface_mesh,fec);
elems = surface_element_map;
break;
default:
MFEM_ABORT("Wrong entity type choice");
break;
}
Array<int> dofs(fesptr->GetVSize());
for (int iel = 0; iel<elems.Size(); ++iel)
{
// index in the global mesh
int iel_idx = elems[iel];
// get the dofs of this element
Array<int> ldofs;
Array<int> gdofs;
switch(etype)
{
case 0: fes0->GetElementVDofs(iel_idx,gdofs); break;
case 1: fes0->GetFaceVDofs(iel_idx,gdofs); break;
default: MFEM_ABORT("Wrong entity type"); break;
}
fesptr->GetElementDofs(iel,ldofs);
// the sizes have to match
MFEM_VERIFY(gdofs.Size() == ldofs.Size(),
"Size inconsistency");
// loop through the dofs and take into account the signs;
int ndof = ldofs.Size();
for (int i = 0; i<ndof; ++i)
{
int ldof_ = ldofs[i];
int gdof_ = gdofs[i];
int ldof = (ldof_ >= 0) ? ldof_ : abs(ldof_) - 1;
int gdof = (gdof_ >= 0) ? gdof_ : abs(gdof_) - 1;
dofs[ldof] = gdof;
}
}
switch(etype)
{
case 0:
dof_map = dofs;
fes = fesptr;
break;
case 1:
surface_dof_map = dofs;
surface_fes = fesptr;
break;
default:
MFEM_ABORT("Wrong entity type"); break;
}
}
void Subdomain::BuildProlongationMatrix(const entity_type & etype)
{
Array<int> dofs;
SparseMatrix * Ptr = nullptr;
switch (etype)
{
case 0:
if (!dof_map.Size()) BuildDofMap(etype);
dofs = dof_map;
Ptr = P;
break;
case 1:
if (!surface_dof_map.Size()) BuildDofMap(etype);
dofs = surface_dof_map;
Ptr = Pf;
break;
default:
MFEM_ABORT("Wrong entity type");
break;
}
int height = fes0->GetVSize();
int width = dofs.Size();
Ptr = new SparseMatrix(height,width);
for (int i = 0; i< dofs.Size(); i++)
{
int j = dofs[i];
Ptr->Set(j,i,1.);
}
Ptr->Finalize();
switch (etype)
{
case 0: P = Ptr; break;
case 1: Pf = Ptr; break;
default: MFEM_ABORT("Wrong entity type"); break;
}
}
Subdomain::~Subdomain()
{
delete mesh;
delete surface_mesh;
delete fes;
delete surface_fes;
delete P;
delete Pf;
}
+106
View File
@@ -0,0 +1,106 @@
#pragma once
#include "mfem.hpp"
#include <fstream>
#include <iostream>
using namespace std;
using namespace mfem;
class Subdomain
{
public:
enum entity_type
{
volume,
surface
};
private:
const Mesh *mesh0=nullptr;
int dim, sdim;
const FiniteElementSpace *fes0=nullptr;
DenseMatrix vcoords;
Mesh *mesh=nullptr; // Submesh
Mesh *surface_mesh=nullptr; //Surface mesh
FiniteElementSpace *fes=nullptr; // FE Space on the submesh
FiniteElementSpace *surface_fes=nullptr; // FE Space on the submesh
Array<int> element_map, surface_element_map;
Array<int> dof_map, surface_dof_map;
SparseMatrix * P=nullptr;
SparseMatrix * Pf=nullptr;
void BuildDofMap(const entity_type & etype);
void BuildProlongationMatrix(const entity_type & etype);
void BuildSubMesh(const Array<int> & elems, const entity_type & etype);
public:
Subdomain(const Mesh & mesh_);
Mesh * GetSubMesh(const Array<int> & elems)
{
if(!mesh) BuildSubMesh(elems, entity_type::volume);
return mesh;
}
Mesh * GetSurfaceMesh(const Array<int> & surface_elems)
{
if (!surface_mesh) BuildSubMesh(surface_elems, entity_type::surface);
return surface_mesh;
}
void SetFESpace(const FiniteElementSpace & fes0_)
{
fes0 = &fes0_;
}
void GetElementMap(Array<int> & element_map_)
{
element_map_ = element_map;
}
void GetFaceElementMap(Array<int> & surface_element_map_)
{
surface_element_map_ = surface_element_map;
}
void GetDofMap(Array<int> & dof_map_)
{
if (!dof_map.Size()) BuildDofMap(entity_type::volume);
dof_map_ = dof_map;
}
void GetSurfaceDofMap(Array<int> & surface_dof_map_)
{
if (!surface_dof_map.Size()) BuildDofMap(entity_type::surface);
surface_dof_map_ = surface_dof_map;
}
SparseMatrix * GetProlonationMatrix()
{
if (!P) BuildProlongationMatrix(entity_type::volume);
return P;
}
SparseMatrix * GetSurfaceProlonationMatrix()
{
if (!Pf) BuildProlongationMatrix(entity_type::surface);
return Pf;
}
FiniteElementSpace * GetSubFESpace(const entity_type & etype)
{
switch (etype)
{
case 0:
if (!fes)
{
MFEM_VERIFY(mesh, "Volume mesh not built");
BuildDofMap(etype);
}
return fes;
break;
case 1:
if (!surface_fes)
{
MFEM_VERIFY(surface_mesh, "Surface mesh not built");
BuildDofMap(etype);
}
return surface_fes;
break;
default:
MFEM_ABORT("Wrong entity type");
return nullptr;
break;
}
}
~Subdomain();
};
+96
View File
@@ -0,0 +1,96 @@
#include "mfem.hpp"
#include <fstream>
#include <iostream>
using namespace std;
using namespace mfem;
int main(int argc, char *argv[])
{
// 1. Parse command line options
const char *mesh_file = "../../data/periodic-annulus-sector.msh";
int order = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
args.ParseCheck();
// 2. Read the mesh from the given mesh file, and refine once uniformly.
Mesh orig_mesh(mesh_file);
orig_mesh.CheckElementOrientation(true);
orig_mesh.CheckBdrElementOrientation(true);
// mesh.EnsureNodes();
// mesh.UniformRefinement();
// Array<int> elems({1,3,21,10,20,2,0});
Array<int> elems({0,1,2,3});
// int nel = orig_mesh.GetNE();
int nel = elems.Size();
// Array<int> elems(nel);
// for (int i = 0; i<nel; i++)
// {
// elems[i] = i;
// }
// elems.Print();
Mesh new_mesh = Mesh::ExtractMesh(orig_mesh,elems);
new_mesh.CheckElementOrientation(true);
new_mesh.CheckBdrElementOrientation(true);
{
char vishost[] = "localhost";
int visport = 19916;
socketstream mesh0_sock(vishost, visport);
mesh0_sock.precision(8);
mesh0_sock << "mesh\n" << orig_mesh << "keys n \n" << flush;
socketstream mesh1_sock(vishost, visport);
mesh1_sock.precision(8);
mesh1_sock << "mesh\n" << new_mesh << "keys n \n" << flush;
}
Array<int> faces;
for (int i = 0; i<orig_mesh.GetNBE(); i++)
{
if (orig_mesh.GetBdrAttribute(i) >= 1)
faces.Append(orig_mesh.GetBdrFace(i));
}
// faces.Append(orig_mesh.GetBdrFace(1));
// faces.Append(orig_mesh.GetBdrFace(2));
Mesh surface_mesh = Mesh::ExtractSurfaceMesh(orig_mesh,faces);
surface_mesh.CheckElementOrientation(true);
surface_mesh.CheckBdrElementOrientation(true);
surface_mesh.Print();
if (surface_mesh.Dimension() > 1)
{
char vishost[] = "localhost";
int visport = 19916;
socketstream mesh2_sock(vishost, visport);
mesh2_sock.precision(8);
mesh2_sock << "mesh\n" << surface_mesh << "keys n \n" << flush;
}
else
{
ParaViewDataCollection paraview_dc("surf_mesh", &surface_mesh);
paraview_dc.SetPrefixPath("ParaView");
paraview_dc.SetLevelsOfDetail(3);
paraview_dc.SetCycle(0);
paraview_dc.SetDataFormat(VTKFormat::BINARY);
paraview_dc.SetHighOrderOutput(true);
paraview_dc.SetTime(0.0); // set the time
H1_FECollection fec(order,surface_mesh.Dimension());
FiniteElementSpace fespace(&surface_mesh,&fec);
GridFunction gf(&fespace);
gf.Randomize();
paraview_dc.RegisterField("solution",&gf);
paraview_dc.Save();
}
return 0;
}
+180
View File
@@ -0,0 +1,180 @@
#include "mfem.hpp"
#include <fstream>
#include <iostream>
#include "mesh_partition.hpp"
using namespace std;
using namespace mfem;
double sin_func(const Vector & x);
int main(int argc, char *argv[])
{
// 1. Parse command line options
const char *mesh_file = "../../data/periodic-annulus-sector.msh";
int order = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
args.ParseCheck();
// 2. Read the mesh from the given mesh file, and refine once uniformly.
Mesh mesh(mesh_file);
// mesh.EnsureNodes();
// mesh.UniformRefinement();
// Array<int> elems0({0,1,2,3});
int nel = mesh.GetNE();
// int nel = 5;
Array<int> elems0(nel/2);
for (int i = 0; i<nel/2; i++)
{
elems0[i] = i;
}
elems0.Print();
// elems0.Append(24);
// elems0.Append(23);
// elems0.Append(26);
Subdomain subdomain0(mesh);
Mesh * submesh = subdomain0.GetSubMesh(elems0);
// cout << "number of boundary elements = " << mesh.GetNBE() << endl;
Array<int> faces(mesh.GetNBE()/2);
for (int i = 0; i<mesh.GetNBE()/2; i++)
{
faces[i] = mesh.GetBdrFace(i);
}
Mesh * surfmesh = subdomain0.GetSurfaceMesh(faces);
H1_FECollection fec(order, mesh.Dimension());
FiniteElementSpace fespace(&mesh, &fec);
FunctionCoefficient coeff(sin_func);
GridFunction gf(&fespace);
gf.ProjectCoefficient(coeff);
{
char vishost[] = "localhost";
int visport = 19916;
socketstream mesh_sock(vishost, visport);
mesh_sock.precision(8);
// mesh_sock << "mesh\n" << mesh << "keys n \n" << flush;
mesh_sock << "solution\n" << mesh << gf << "keys jnmR \n"
<< "valuerange 0 1.0 \n" << flush;
// << flush;
}
subdomain0.SetFESpace(fespace);
SparseMatrix * P = subdomain0.GetProlonationMatrix();
FiniteElementSpace * elem_fes =
subdomain0.GetSubFESpace(Subdomain::entity_type::volume);
GridFunction gf_e(elem_fes);
cout << "Size P = " << P->Height() << " x " << P->Width() << endl;
cout << "gf_e.Size = " << gf_e.Size() << endl;
cout << "gf.Size = " << gf.Size() << endl;
P->MultTranspose(gf,gf_e);
SparseMatrix * Pb = subdomain0.GetSurfaceProlonationMatrix();
FiniteElementSpace * bdr_elem_fes =
subdomain0.GetSubFESpace(Subdomain::entity_type::surface);
GridFunction gf_b(bdr_elem_fes);
Pb->MultTranspose(gf,gf_b);
{
char vishost[] = "localhost";
int visport = 19916;
if (submesh)
{
socketstream mesh0_sock(vishost, visport);
mesh0_sock.precision(8);
// mesh0_sock << "mesh\n" << *submesh << "keys n \n" << flush;
mesh0_sock << "solution\n" << *submesh << gf_e << "keys nmR \n"
<< "valuerange 0 1.0 \n" << flush;
// << flush;
}
if (surfmesh && mesh.Dimension() == 3)
{
socketstream mesh1_sock(vishost, visport);
mesh1_sock.precision(8);
// mesh1_sock << "mesh\n" << *bdrmesh0 << "keys n \n" << flush;
mesh1_sock << "solution\n" << *surfmesh << gf_b
<< "valuerange 0 1.0 \n" << flush;
// << flush;
}
}
// ParaViewDataCollection paraview_dc("mesh_partition", surfmesh);
// paraview_dc.SetPrefixPath("ParaView");
// const FiniteElementSpace * fes_ = surfmesh->GetNodalFESpace();
// int ord = (fes_) ? fes_->GetOrder(0) : order;
// paraview_dc.SetLevelsOfDetail(ord);
// paraview_dc.SetCycle(0);
// paraview_dc.SetDataFormat(VTKFormat::BINARY);
// paraview_dc.SetHighOrderOutput(true);
// paraview_dc.SetTime(0.0); // set the time
// paraview_dc.RegisterField("solution",&gf_b);
// paraview_dc.Save();
// // ---------------------------------------------------------
// FiniteElementCollection *fec1 = new H1_FECollection(order, submesh->Dimension());
// FiniteElementSpace fespace1(submesh, fec1);
// Array<int> ess_tdof_list;
// if (submesh->bdr_attributes.Size())
// {
// Array<int> ess_bdr(mesh.bdr_attributes.Max());
// ess_bdr = 1;
// fespace1.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
// }
// LinearForm b(&fespace1);
// ConstantCoefficient one(1.0);
// b.AddDomainIntegrator(new DomainLFIntegrator(one));
// b.Assemble();
// GridFunction x(&fespace1);
// x = 0.0;
// // 9. Set up the bilinear form a(.,.) on the finite element space
// // corresponding to the Laplacian operator -Delta, by adding the Diffusion
// // domain integrator.
// BilinearForm a(&fespace1);
// a.AddDomainIntegrator(new DiffusionIntegrator(one));
// a.Assemble();
// OperatorPtr A;
// Vector B, X;
// a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
// cout << "Size of linear system: " << A->Height() << endl;
// // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
// GSSmoother M((SparseMatrix&)(*A));
// PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
// // 12. Recover the solution as a finite element grid function.
// a.RecoverFEMSolution(X, b, x);
// {
// char vishost[] = "localhost";
// int visport = 19916;
// socketstream sol_sock2(vishost, visport);
// sol_sock2.precision(8);
// sol_sock2 << "solution\n" << *submesh << x << flush;
// }
return 0;
}
double sin_func(const Vector & x)
{
Vector c(x.Size());
c.Randomize();
// double dotp = c*x;
// return (sin(10.0*M_PI*dotp));
// return sin(2.*M_PI*x[0]);
// return 1.-x[1]*x[1]/4.0;
// return (0.5-x[1])*(0.5-x[1]);
return x[1];
// double r = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
// return r;
}
+178
View File
@@ -0,0 +1,178 @@
#include "mfem.hpp"
#include <fstream>
#include <iostream>
#include "mesh_partition.hpp"
using namespace std;
using namespace mfem;
double sin_func(const Vector & x);
int main(int argc, char *argv[])
{
// 1. Parse command line options
const char *mesh_file = "../data/star.mesh";
int order = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
args.ParseCheck();
// 2. Read the mesh from the given mesh file, and refine once uniformly.
Mesh mesh(mesh_file);
// mesh.Print(cout);
mesh.UniformRefinement();
// mesh.EnsureNodes();
// Array<int> elems0({0,4,8,12,16});
// Array<int> elems0({0,4,8,12,16});
// Array<int> elems0({6,7,8});
Array<int> elems0({0,1,2,3,4});
// Array<int> elems0({7,6,17,20,21,22});
// Array<int> elems0({0,1,2,3,4});
// Array<int> elems0({104,103,86,109});
// elems0.Print(cout, elems0.Size());
Subdomain subdomain0(mesh);
Mesh * submesh0 = subdomain0.GetSubMesh(elems0);
// Array<int> bdrelems0({8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23});
cout << "number of boundary elements = " << mesh.GetNBE() << endl;
Array<int> bdrelems0({0,1});
Array<int> faces(bdrelems0.Size());
for (int i = 0; i<bdrelems0.Size(); i++)
{
faces[i] = mesh.GetBdrFace(bdrelems0[i]);
}
Mesh * surfmesh0 = subdomain0.GetSurfaceMesh(faces);
H1_FECollection fec(order, mesh.Dimension());
FiniteElementSpace fespace(&mesh, &fec);
FunctionCoefficient coeff(sin_func);
GridFunction gf(&fespace);
gf.ProjectCoefficient(coeff);
{
char vishost[] = "localhost";
int visport = 19916;
socketstream mesh_sock(vishost, visport);
mesh_sock.precision(8);
mesh_sock << "solution\n" << mesh << gf
<< "valuerange -1.0 1.0 \n" << flush;
}
subdomain0.SetFESpace(fespace);
SparseMatrix * P = subdomain0.GetProlonationMatrix();
FiniteElementSpace * elem_fes =
subdomain0.GetSubFESpace(Subdomain::entity_type::volume);
GridFunction gf_e(elem_fes);
P->MultTranspose(gf,gf_e);
SparseMatrix * Pf = subdomain0.GetSurfaceProlonationMatrix();
FiniteElementSpace * face_elem_fes =
subdomain0.GetSubFESpace(Subdomain::entity_type::surface);
GridFunction gf_f(face_elem_fes);
Pf->MultTranspose(gf,gf_f);
{
char vishost[] = "localhost";
int visport = 19916;
if (submesh0)
{
socketstream mesh0_sock(vishost, visport);
mesh0_sock.precision(8);
// mesh0_sock << "mesh\n" << *submesh0 << "keys n \n" << flush;
mesh0_sock << "solution\n" << *submesh0 << gf_e
<< "valuerange -1.0 1.0 \n" << flush;
}
if (surfmesh0 && mesh.Dimension()==3)
{
socketstream mesh1_sock(vishost, visport);
mesh1_sock.precision(8);
// mesh1_sock << "mesh\n" << *bdrmesh0 << "keys n \n" << flush;
mesh1_sock << "solution\n" << *surfmesh0 << gf_f
<< "valuerange -1.0 1.0 \n" << flush;
}
}
// Array<int> bdr_faces;
// for (int i =0; i<mesh.GetNBE(); i++)
// {
// int attr = mesh.GetBdrAttribute(i);
// if (attr == 4)
// {
// bdr_faces.Append(mesh.GetBdrFace(i));
// }
// }
// H1_FECollection fec(order, mesh.Dimension());
// FiniteElementSpace fespace(&mesh, &fec);
// FunctionCoefficient coeff(sin_func);
// GridFunction gf(&fespace);
// gf.ProjectCoefficient(coeff);
// {
// char vishost[] = "localhost";
// int visport = 19916;
// socketstream mesh_sock(vishost, visport);
// mesh_sock.precision(8);
// // mesh_sock << "mesh\n" << mesh << "keys n \n" << flush;
// mesh_sock << "solution\n" << mesh << gf << flush;
// // << "valuerange -5000.0 5000.0 \n" << flush;
// // << "valuerange -1.0 1.0 \n" << flush;
// }
// Subdomain subdomain1(mesh);
// subdomain1.SetFESpace(fespace);
// Mesh * bdrmesh0 = subdomain1.GetBdrSurfaceMesh(bdr_faces);
// SparseMatrix * Pb = subdomain1.GetBdrProlonationMatrix();
// FiniteElementSpace * bdr_elem_fes =
// subdomain1.GetSubFESpace(Subdomain::entity_type::bdr);
// GridFunction gf_f(bdr_elem_fes);
// Pb->MultTranspose(gf,gf_f);
// // gf.Print();
// // gf_f.Print();
// // bdrmesh0->Print(cout);
// if (bdrmesh0)
// {
// char vishost[] = "localhost";
// int visport = 19916;
// socketstream mesh1_sock(vishost, visport);
// mesh1_sock.precision(8);
// // mesh1_sock << "mesh\n" << *bdrmesh0 << "keys n \n" << flush;
// mesh1_sock << "solution\n" << *bdrmesh0 << gf_f << flush;
// // << "valuerange -5000.0 5000.0 \n" << flush;
// }
ParaViewDataCollection paraview_dc("mesh_partition", surfmesh0);
paraview_dc.SetPrefixPath("ParaView");
const FiniteElementSpace * fes_ = surfmesh0->GetNodalFESpace();
int ord = (fes_) ? fes_->GetOrder(0) : order;
paraview_dc.SetLevelsOfDetail(5);
paraview_dc.SetCycle(0);
paraview_dc.SetDataFormat(VTKFormat::BINARY);
paraview_dc.SetHighOrderOutput(true);
paraview_dc.SetTime(0.0); // set the time
paraview_dc.RegisterField("solution",&gf_f);
paraview_dc.Save();
return 0;
}
double sin_func(const Vector & x)
{
Vector c(x.Size());
c.Randomize();
// double dotp = c*x;
double dotp = x.Sum();
// return (sin(10.0*M_PI*dotp));
// return sin(2.*M_PI*x[0]);
// return 1.-x[1]*x[1]/4.0;
return (0.5-x[1])*(0.5-x[1]);
// double r = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
// return r;
}
+1 -1
View File
@@ -105,7 +105,7 @@ private:
Vector diag(fespace.GetTrueVSize());
bfs.Last()->AssembleDiagonal(diag);
Solver* smoother = new OperatorChebyshevSmoother(*opr, diag,
Solver* smoother = new OperatorChebyshevSmoother(opr.Ptr(), diag,
*essentialTrueDofs.Last(), 2);
AddLevel(opr.Ptr(), smoother, true, true);
}
+1 -1
View File
@@ -115,7 +115,7 @@ private:
Vector diag(fespace.GetTrueVSize());
bfs.Last()->AssembleDiagonal(diag);
Solver* smoother = new OperatorChebyshevSmoother(*opr, diag,
Solver* smoother = new OperatorChebyshevSmoother(opr.Ptr(), diag,
*essentialTrueDofs.Last(), 2, fespace.GetParMesh()->GetComm());
AddLevel(opr.Ptr(), smoother, true, true);
+148 -171
View File
@@ -31,7 +31,7 @@ void BilinearForm::AllocMat()
const Table &elem_dof = fes->GetElementToDofTable();
Table dof_dof;
if (interior_face_integs.Size() > 0)
if (fbfi.Size() > 0)
{
// the sparsity pattern is defined from the map: face->element->dof
Table face_dof, dof_face;
@@ -99,15 +99,15 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
ext = NULL;
// Copy the pointers to the integrators
domain_integs = bf->domain_integs;
dbfi = bf->dbfi;
boundary_integs = bf->boundary_integs;
boundary_integs_marker = bf->boundary_integs_marker;
bbfi = bf->bbfi;
bbfi_marker = bf->bbfi_marker;
interior_face_integs = bf->interior_face_integs;
fbfi = bf->fbfi;
boundary_face_integs = bf->boundary_face_integs;
boundary_face_integs_marker = bf->boundary_face_integs_marker;
bfbfi = bf->bfbfi;
bfbfi_marker = bf->bfbfi_marker;
AllocMat();
}
@@ -234,47 +234,46 @@ void BilinearForm::Finalize (int skip_zeros)
void BilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi)
{
domain_integs.Append(bfi);
domain_integs_marker.Append(NULL); // NULL marker means apply everywhere
dbfi.Append(bfi);
dbfi_marker.Append(NULL); // NULL marker means apply everywhere
}
void BilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi,
Array<int> &elem_marker)
{
domain_integs.Append(bfi);
domain_integs_marker.Append(&elem_marker);
dbfi.Append(bfi);
dbfi_marker.Append(&elem_marker);
}
void BilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi)
{
boundary_integs.Append (bfi);
boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
bbfi.Append (bfi);
bbfi_marker.Append(NULL); // NULL marker means apply everywhere
}
void BilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
Array<int> &bdr_marker)
{
boundary_integs.Append (bfi);
boundary_integs_marker.Append(&bdr_marker);
bbfi.Append (bfi);
bbfi_marker.Append(&bdr_marker);
}
void BilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi)
{
interior_face_integs.Append (bfi);
fbfi.Append (bfi);
}
void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
{
boundary_face_integs.Append(bfi);
// NULL marker means apply everywhere
boundary_face_integs_marker.Append(NULL);
bfbfi.Append(bfi);
bfbfi_marker.Append(NULL); // NULL marker means apply everywhere
}
void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
Array<int> &bdr_marker)
{
boundary_face_integs.Append(bfi);
boundary_face_integs_marker.Append(&bdr_marker);
bfbfi.Append(bfi);
bfbfi_marker.Append(&bdr_marker);
}
void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
@@ -286,14 +285,14 @@ void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
return;
}
if (domain_integs.Size())
if (dbfi.Size())
{
const FiniteElement &fe = *fes->GetFE(i);
ElementTransformation *eltrans = fes->GetElementTransformation(i);
domain_integs[0]->AssembleElementMatrix(fe, *eltrans, elmat);
for (int k = 1; k < domain_integs.Size(); k++)
dbfi[0]->AssembleElementMatrix(fe, *eltrans, elmat);
for (int k = 1; k < dbfi.Size(); k++)
{
domain_integs[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
elmat += elemmat;
}
}
@@ -307,14 +306,14 @@ void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
void BilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
{
if (boundary_integs.Size())
if (bbfi.Size())
{
const FiniteElement &be = *fes->GetBE(i);
ElementTransformation *eltrans = fes->GetBdrElementTransformation(i);
boundary_integs[0]->AssembleElementMatrix(be, *eltrans, elmat);
for (int k = 1; k < boundary_integs.Size(); k++)
bbfi[0]->AssembleElementMatrix(be, *eltrans, elmat);
for (int k = 1; k < bbfi.Size(); k++)
{
boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
elmat += elemmat;
}
}
@@ -408,14 +407,13 @@ void BilinearForm::Assemble(int skip_zeros)
}
#endif
if (domain_integs.Size())
if (dbfi.Size())
{
for (int k = 0; k < domain_integs.Size(); k++)
for (int k = 0; k < dbfi.Size(); k++)
{
if (domain_integs_marker[k] != NULL)
if (dbfi_marker[k] != NULL)
{
MFEM_VERIFY(mesh->attributes.Size() ==
domain_integs_marker[k]->Size(),
MFEM_VERIFY(mesh->attributes.Size() == dbfi_marker[k]->Size(),
"invalid element marker for domain integrator #"
<< k << ", counting from zero");
}
@@ -432,14 +430,14 @@ void BilinearForm::Assemble(int skip_zeros)
else
{
elmat.SetSize(0);
for (int k = 0; k < domain_integs.Size(); k++)
for (int k = 0; k < dbfi.Size(); k++)
{
if ( domain_integs_marker[k] == NULL ||
(*(domain_integs_marker[k]))[elem_attr-1] == 1)
if ( dbfi_marker[k] == NULL ||
(*(dbfi_marker[k]))[elem_attr-1] == 1)
{
const FiniteElement &fe = *fes->GetFE(i);
eltrans = fes->GetElementTransformation(i);
domain_integs[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
dbfi[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
if (elmat.Size() == 0)
{
elmat = elemmat;
@@ -474,20 +472,20 @@ void BilinearForm::Assemble(int skip_zeros)
}
}
if (boundary_integs.Size())
if (bbfi.Size())
{
// Which boundary attributes need to be processed?
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < boundary_integs.Size(); k++)
for (int k = 0; k < bbfi.Size(); k++)
{
if (boundary_integs_marker[k] == NULL)
if (bbfi_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *boundary_integs_marker[k];
Array<int> &bdr_marker = *bbfi_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary integrator #"
<< k << ", counting from zero");
@@ -506,21 +504,21 @@ void BilinearForm::Assemble(int skip_zeros)
fes -> GetBdrElementVDofs (i, vdofs);
eltrans = fes -> GetBdrElementTransformation (i);
int k = 0;
for (; k < boundary_integs.Size(); k++)
for (; k < bbfi.Size(); k++)
{
if (boundary_integs_marker[k] &&
(*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
if (bbfi_marker[k] &&
(*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elmat);
bbfi[k]->AssembleElementMatrix(be, *eltrans, elmat);
k++;
break;
}
for (; k < boundary_integs.Size(); k++)
for (; k < bbfi.Size(); k++)
{
if (boundary_integs_marker[k] &&
(*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
if (bbfi_marker[k] &&
(*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
boundary_integs[k]->AssembleElementMatrix(be, *eltrans, elemmat);
bbfi[k]->AssembleElementMatrix(be, *eltrans, elemmat);
elmat += elemmat;
}
if (!static_cond)
@@ -538,7 +536,7 @@ void BilinearForm::Assemble(int skip_zeros)
}
}
if (interior_face_integs.Size())
if (fbfi.Size())
{
FaceElementTransformations *tr;
Array<int> vdofs2;
@@ -552,19 +550,18 @@ void BilinearForm::Assemble(int skip_zeros)
fes -> GetElementVDofs (tr -> Elem1No, vdofs);
fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
vdofs.Append (vdofs2);
for (int k = 0; k < interior_face_integs.Size(); k++)
for (int k = 0; k < fbfi.Size(); k++)
{
interior_face_integs[k]->
AssembleFaceMatrix(*fes->GetFE(tr->Elem1No),
*fes->GetFE(tr->Elem2No),
*tr, elemmat);
fbfi[k] -> AssembleFaceMatrix (*fes -> GetFE (tr -> Elem1No),
*fes -> GetFE (tr -> Elem2No),
*tr, elemmat);
mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
}
}
}
}
if (boundary_face_integs.Size())
if (bfbfi.Size())
{
FaceElementTransformations *tr;
const FiniteElement *fe1, *fe2;
@@ -573,14 +570,14 @@ void BilinearForm::Assemble(int skip_zeros)
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < boundary_face_integs.Size(); k++)
for (int k = 0; k < bfbfi.Size(); k++)
{
if (boundary_face_integs_marker[k] == NULL)
if (bfbfi_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *boundary_face_integs_marker[k];
Array<int> &bdr_marker = *bfbfi_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary face integrator #"
<< k << ", counting from zero");
@@ -604,14 +601,12 @@ void BilinearForm::Assemble(int skip_zeros)
// but we can't dereference a NULL pointer, and we don't want to
// actually make a fake element.
fe2 = fe1;
for (int k = 0; k < boundary_face_integs.Size(); k++)
for (int k = 0; k < bfbfi.Size(); k++)
{
if (boundary_face_integs_marker[k] &&
(*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
{ continue; }
if (bfbfi_marker[k] &&
(*bfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
elemmat);
bfbfi[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr, elemmat);
mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
}
}
@@ -862,7 +857,7 @@ void BilinearForm::RecoverFEMSolution(const Vector &X,
void BilinearForm::ComputeElementMatrices()
{
if (element_matrices || domain_integs.Size() == 0 || fes->GetNE() == 0)
if (element_matrices || dbfi.Size() == 0 || fes->GetNE() == 0)
{
return;
}
@@ -891,11 +886,11 @@ void BilinearForm::ComputeElementMatrices()
#endif
fes->GetElementTransformation(i, &eltrans);
domain_integs[0]->AssembleElementMatrix(fe, eltrans, elmat);
for (int k = 1; k < domain_integs.Size(); k++)
dbfi[0]->AssembleElementMatrix(fe, eltrans, elmat);
for (int k = 1; k < dbfi.Size(); k++)
{
// note: some integrators may not be thread-safe
domain_integs[k]->AssembleElementMatrix(fe, eltrans, tmp);
dbfi[k]->AssembleElementMatrix(fe, eltrans, tmp);
elmat += tmp;
}
elmat.ClearExternalData();
@@ -1110,12 +1105,10 @@ BilinearForm::~BilinearForm()
if (!extern_bfs)
{
int k;
for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
for (k=0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
for (k=0; k < interior_face_integs.Size(); k++)
{ delete interior_face_integs[k]; }
for (k=0; k < boundary_face_integs.Size(); k++)
{ delete boundary_face_integs[k]; }
for (k=0; k < dbfi.Size(); k++) { delete dbfi[k]; }
for (k=0; k < bbfi.Size(); k++) { delete bbfi[k]; }
for (k=0; k < fbfi.Size(); k++) { delete fbfi[k]; }
for (k=0; k < bfbfi.Size(); k++) { delete bfbfi[k]; }
}
delete ext;
@@ -1148,13 +1141,13 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
ext = NULL;
// Copy the pointers to the integrators
domain_integs = mbf->domain_integs;
boundary_integs = mbf->boundary_integs;
trace_face_integs = mbf->trace_face_integs;
boundary_trace_face_integs = mbf->boundary_trace_face_integs;
dbfi = mbf->dbfi;
bbfi = mbf->bbfi;
tfbfi = mbf->tfbfi;
btfbfi = mbf->btfbfi;
boundary_integs_marker = mbf->boundary_integs_marker;
boundary_trace_face_integs_marker = mbf->boundary_trace_face_integs_marker;
bbfi_marker = mbf->bbfi_marker;
btfbfi_marker = mbf->btfbfi_marker;
assembly = AssemblyLevel::LEGACY;
ext = NULL;
@@ -1243,8 +1236,7 @@ MatrixInverse * MixedBilinearForm::Inverse() const
{
if (assembly != AssemblyLevel::LEGACY)
{
MFEM_WARNING("MixedBilinearForm::Inverse not possible with this "
"assembly level!");
MFEM_WARNING("MixedBilinearForm::Inverse not possible with this assembly level!");
return NULL;
}
else
@@ -1275,39 +1267,38 @@ void MixedBilinearForm::GetBlocks(Array2D<SparseMatrix *> &blocks) const
void MixedBilinearForm::AddDomainIntegrator (BilinearFormIntegrator * bfi)
{
domain_integs.Append (bfi);
dbfi.Append (bfi);
}
void MixedBilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi)
{
boundary_integs.Append (bfi);
boundary_integs_marker.Append(NULL); // NULL marker means apply everywhere
bbfi.Append (bfi);
bbfi_marker.Append(NULL); // NULL marker means apply everywhere
}
void MixedBilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
Array<int> &bdr_marker)
{
boundary_integs.Append (bfi);
boundary_integs_marker.Append(&bdr_marker);
bbfi.Append (bfi);
bbfi_marker.Append(&bdr_marker);
}
void MixedBilinearForm::AddTraceFaceIntegrator (BilinearFormIntegrator * bfi)
{
trace_face_integs.Append (bfi);
tfbfi.Append (bfi);
}
void MixedBilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi)
{
boundary_trace_face_integs.Append(bfi);
// NULL marker means apply everywhere
boundary_trace_face_integs_marker.Append(NULL);
btfbfi.Append(bfi);
btfbfi_marker.Append(NULL); // NULL marker means apply everywhere
}
void MixedBilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi,
Array<int> &bdr_marker)
{
boundary_trace_face_integs.Append(bfi);
boundary_trace_face_integs_marker.Append(&bdr_marker);
btfbfi.Append(bfi);
btfbfi_marker.Append(&bdr_marker);
}
void MixedBilinearForm::Assemble (int skip_zeros)
@@ -1329,37 +1320,37 @@ void MixedBilinearForm::Assemble (int skip_zeros)
mat = new SparseMatrix(height, width);
}
if (domain_integs.Size())
if (dbfi.Size())
{
for (int i = 0; i < test_fes -> GetNE(); i++)
{
trial_fes -> GetElementVDofs (i, tr_vdofs);
test_fes -> GetElementVDofs (i, te_vdofs);
eltrans = test_fes -> GetElementTransformation (i);
for (int k = 0; k < domain_integs.Size(); k++)
for (int k = 0; k < dbfi.Size(); k++)
{
domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
*test_fes -> GetFE(i),
*eltrans, elemmat);
dbfi[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
*test_fes -> GetFE(i),
*eltrans, elemmat);
mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
}
}
}
if (boundary_integs.Size())
if (bbfi.Size())
{
// Which boundary attributes need to be processed?
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < boundary_integs.Size(); k++)
for (int k = 0; k < bbfi.Size(); k++)
{
if (boundary_integs_marker[k] == NULL)
if (bbfi_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *boundary_integs_marker[k];
Array<int> &bdr_marker = *bbfi_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary integrator #"
<< k << ", counting from zero");
@@ -1377,20 +1368,20 @@ void MixedBilinearForm::Assemble (int skip_zeros)
trial_fes -> GetBdrElementVDofs (i, tr_vdofs);
test_fes -> GetBdrElementVDofs (i, te_vdofs);
eltrans = test_fes -> GetBdrElementTransformation (i);
for (int k = 0; k < boundary_integs.Size(); k++)
for (int k = 0; k < bbfi.Size(); k++)
{
if (boundary_integs_marker[k] &&
(*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
if (bbfi_marker[k] &&
(*bbfi_marker[k])[bdr_attr-1] == 0) { continue; }
boundary_integs[k]->AssembleElementMatrix2 (*trial_fes -> GetBE(i),
*test_fes -> GetBE(i),
*eltrans, elemmat);
bbfi[k] -> AssembleElementMatrix2 (*trial_fes -> GetBE(i),
*test_fes -> GetBE(i),
*eltrans, elemmat);
mat -> AddSubMatrix (te_vdofs, tr_vdofs, elemmat, skip_zeros);
}
}
}
if (trace_face_integs.Size())
if (tfbfi.Size())
{
FaceElementTransformations *ftr;
Array<int> te_vdofs2;
@@ -1417,16 +1408,16 @@ void MixedBilinearForm::Assemble (int skip_zeros)
// want to actually make a fake element.
test_fe2 = test_fe1;
}
for (int k = 0; k < trace_face_integs.Size(); k++)
for (int k = 0; k < tfbfi.Size(); k++)
{
trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1,
*test_fe2, *ftr, elemmat);
tfbfi[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
*ftr, elemmat);
mat->AddSubMatrix(te_vdofs, tr_vdofs, elemmat, skip_zeros);
}
}
}
if (boundary_trace_face_integs.Size())
if (btfbfi.Size())
{
FaceElementTransformations *ftr;
Array<int> te_vdofs2;
@@ -1436,17 +1427,17 @@ void MixedBilinearForm::Assemble (int skip_zeros)
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
for (int k = 0; k < btfbfi.Size(); k++)
{
if (boundary_trace_face_integs_marker[k] == NULL)
if (btfbfi_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *boundary_trace_face_integs_marker[k];
Array<int> &bdr_marker = *btfbfi_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary trace face"
"integrator #" << k << ", counting from zero");
"invalid boundary marker for boundary trace face integrator #"
<< k << ", counting from zero");
for (int i = 0; i < bdr_attr_marker.Size(); i++)
{
bdr_attr_marker[i] |= bdr_marker[i];
@@ -1469,16 +1460,13 @@ void MixedBilinearForm::Assemble (int skip_zeros)
// boundaries, but we can't dereference a NULL pointer, and we don't
// want to actually make a fake element.
test_fe2 = test_fe1;
for (int k = 0; k < boundary_trace_face_integs.Size(); k++)
for (int k = 0; k < btfbfi.Size(); k++)
{
if (boundary_trace_face_integs_marker[k] &&
(*boundary_trace_face_integs_marker[k])[bdr_attr-1] == 0)
{ continue; }
if (btfbfi_marker[k] &&
(*btfbfi_marker[k])[bdr_attr-1] == 0) { continue; }
boundary_trace_face_integs[k]->AssembleFaceMatrix(*trial_face_fe,
*test_fe1,
*test_fe2,
*ftr, elemmat);
btfbfi[k]->AssembleFaceMatrix(*trial_face_fe, *test_fe1, *test_fe2,
*ftr, elemmat);
mat->AddSubMatrix(te_vdofs, tr_vdofs, elemmat, skip_zeros);
}
}
@@ -1569,17 +1557,15 @@ void MixedBilinearForm::ConformingAssemble()
void MixedBilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
{
if (domain_integs.Size())
if (dbfi.Size())
{
const FiniteElement &trial_fe = *trial_fes->GetFE(i);
const FiniteElement &test_fe = *test_fes->GetFE(i);
ElementTransformation *eltrans = test_fes->GetElementTransformation(i);
domain_integs[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
elmat);
for (int k = 1; k < domain_integs.Size(); k++)
dbfi[0]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans, elmat);
for (int k = 1; k < dbfi.Size(); k++)
{
domain_integs[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans,
elemmat);
dbfi[k]->AssembleElementMatrix2(trial_fe, test_fe, *eltrans, elemmat);
elmat += elemmat;
}
}
@@ -1594,17 +1580,15 @@ void MixedBilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
void MixedBilinearForm::ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
{
if (boundary_integs.Size())
if (bbfi.Size())
{
const FiniteElement &trial_be = *trial_fes->GetBE(i);
const FiniteElement &test_be = *test_fes->GetBE(i);
ElementTransformation *eltrans = test_fes->GetBdrElementTransformation(i);
boundary_integs[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
elmat);
for (int k = 1; k < boundary_integs.Size(); k++)
bbfi[0]->AssembleElementMatrix2(trial_be, test_be, *eltrans, elmat);
for (int k = 1; k < bbfi.Size(); k++)
{
boundary_integs[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans,
elemmat);
bbfi[k]->AssembleElementMatrix2(trial_be, test_be, *eltrans, elemmat);
elmat += elemmat;
}
}
@@ -1704,10 +1688,10 @@ void MixedBilinearForm::EliminateTestDofs (const Array<int> &bdr_attr_is_ess)
}
}
void MixedBilinearForm::FormRectangularSystemMatrix(
const Array<int> &trial_tdof_list,
const Array<int> &test_tdof_list,
OperatorHandle &A)
void MixedBilinearForm::FormRectangularSystemMatrix(const Array<int>
&trial_tdof_list,
const Array<int> &test_tdof_list,
OperatorHandle &A)
{
if (ext)
@@ -1745,17 +1729,17 @@ void MixedBilinearForm::FormRectangularSystemMatrix(
A.Reset(mat, false);
}
void MixedBilinearForm::FormRectangularLinearSystem(
const Array<int> &trial_tdof_list,
const Array<int> &test_tdof_list,
Vector &x, Vector &b,
OperatorHandle &A,
Vector &X, Vector &B)
void MixedBilinearForm::FormRectangularLinearSystem(const Array<int>
&trial_tdof_list,
const Array<int> &test_tdof_list,
Vector &x, Vector &b,
OperatorHandle &A,
Vector &X, Vector &B)
{
if (ext)
{
ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list,
x, b, A, X, B);
ext->FormRectangularLinearSystem(trial_tdof_list, test_tdof_list, x, b, A, X,
B);
return;
}
@@ -1793,13 +1777,10 @@ MixedBilinearForm::~MixedBilinearForm()
if (!extern_bfs)
{
int i;
for (i = 0; i < domain_integs.Size(); i++) { delete domain_integs[i]; }
for (i = 0; i < boundary_integs.Size(); i++)
{ delete boundary_integs[i]; }
for (i = 0; i < trace_face_integs.Size(); i++)
{ delete trace_face_integs[i]; }
for (i = 0; i < boundary_trace_face_integs.Size(); i++)
{ delete boundary_trace_face_integs[i]; }
for (i = 0; i < dbfi.Size(); i++) { delete dbfi[i]; }
for (i = 0; i < bbfi.Size(); i++) { delete bbfi[i]; }
for (i = 0; i < tfbfi.Size(); i++) { delete tfbfi[i]; }
for (i = 0; i < btfbfi.Size(); i++) { delete btfbfi[i]; }
}
delete ext;
}
@@ -1849,7 +1830,7 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
mat = new SparseMatrix(height, width);
}
if (domain_integs.Size() > 0)
if (dbfi.Size() > 0)
{
for (int i = 0; i < test_fes->GetNE(); i++)
{
@@ -1859,19 +1840,17 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
dom_fe = trial_fes->GetFE(i);
ran_fe = test_fes->GetFE(i);
domain_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
totelmat);
for (int j = 1; j < domain_integs.Size(); j++)
dbfi[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
for (int j = 1; j < dbfi.Size(); j++)
{
domain_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
elmat);
dbfi[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
totelmat += elmat;
}
mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
}
}
if (trace_face_integs.Size())
if (tfbfi.Size())
{
const int nfaces = test_fes->GetMesh()->GetNumFaces();
for (int i = 0; i < nfaces; i++)
@@ -1882,12 +1861,10 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
dom_fe = trial_fes->GetFaceElement(i);
ran_fe = test_fes->GetFaceElement(i);
trace_face_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
totelmat);
for (int j = 1; j < trace_face_integs.Size(); j++)
tfbfi[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, totelmat);
for (int j = 1; j < tfbfi.Size(); j++)
{
trace_face_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
elmat);
tfbfi[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T, elmat);
totelmat += elmat;
}
mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
+30 -36
View File
@@ -84,29 +84,28 @@ protected:
the BilinearForm. */
long sequence;
/** @brief Indicates the BilinearFormIntegrator%s stored in #domain_integs,
#boundary_integs, #interior_face_integs, and #boundary_face_integs are
owned by another BilinearForm. */
/** @brief Indicates the BilinearFormIntegrator%s stored in #dbfi, #bbfi,
#fbfi, and #bfbfi are owned by another BilinearForm. */
int extern_bfs;
/// Set of Domain Integrators to be applied.
Array<BilinearFormIntegrator*> domain_integs;
Array<BilinearFormIntegrator*> dbfi;
/// Element attribute marker (should be of length mesh->attributes)
/// Includes all by default.
/// 0 - ignore attribute
/// 1 - include attribute
Array<Array<int>*> domain_integs_marker;
Array<Array<int>*> dbfi_marker;
/// Set of Boundary Integrators to be applied.
Array<BilinearFormIntegrator*> boundary_integs;
Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned.
Array<BilinearFormIntegrator*> bbfi;
Array<Array<int>*> bbfi_marker; ///< Entries are not owned.
/// Set of interior face Integrators to be applied.
Array<BilinearFormIntegrator*> interior_face_integs;
Array<BilinearFormIntegrator*> fbfi;
/// Set of boundary face Integrators to be applied.
Array<BilinearFormIntegrator*> boundary_face_integs;
Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned.
Array<BilinearFormIntegrator*> bfbfi;
Array<Array<int>*> bfbfi_marker; ///< Entries are not owned.
DenseMatrix elemmat;
Array<int> vdofs;
@@ -232,25 +231,24 @@ public:
void AllocateMatrix() { if (mat == NULL) { AllocMat(); } }
/// Access all the integrators added with AddDomainIntegrator().
Array<BilinearFormIntegrator*> *GetDBFI() { return &domain_integs; }
Array<BilinearFormIntegrator*> *GetDBFI() { return &dbfi; }
/// Access all the integrators added with AddBoundaryIntegrator().
Array<BilinearFormIntegrator*> *GetBBFI() { return &boundary_integs; }
Array<BilinearFormIntegrator*> *GetBBFI() { return &bbfi; }
/** @brief Access all boundary markers added with AddBoundaryIntegrator().
If no marker was specified when the integrator was added, the
corresponding pointer (to Array<int>) will be NULL. */
Array<Array<int>*> *GetBBFI_Marker() { return &boundary_integs_marker; }
Array<Array<int>*> *GetBBFI_Marker() { return &bbfi_marker; }
/// Access all integrators added with AddInteriorFaceIntegrator().
Array<BilinearFormIntegrator*> *GetFBFI() { return &interior_face_integs; }
Array<BilinearFormIntegrator*> *GetFBFI() { return &fbfi; }
/// Access all integrators added with AddBdrFaceIntegrator().
Array<BilinearFormIntegrator*> *GetBFBFI() { return &boundary_face_integs; }
Array<BilinearFormIntegrator*> *GetBFBFI() { return &bfbfi; }
/** @brief Access all boundary markers added with AddBdrFaceIntegrator().
If no marker was specified when the integrator was added, the
corresponding pointer (to Array<int>) will be NULL. */
Array<Array<int>*> *GetBFBFI_Marker()
{ return &boundary_face_integs_marker; }
Array<Array<int>*> *GetBFBFI_Marker() { return &bfbfi_marker; }
/// Returns a reference to: \f$ M_{ij} \f$
const double &operator()(int i, int j) { return (*mat)(i,j); }
@@ -654,25 +652,23 @@ protected:
Partial Assembly (PA), or Matrix Free assembly (MF). */
MixedBilinearFormExtension *ext;
/** @brief Indicates the BilinearFormIntegrator%s stored in #domain_integs,
#boundary_integs, #trace_face_integs and #boundary_trace_face_integs
are owned by another MixedBilinearForm. */
/** @brief Indicates the BilinearFormIntegrator%s stored in #dbfi, #bbfi,
#tfbfi and #btfbfi are owned by another MixedBilinearForm. */
int extern_bfs;
/// Domain integrators.
Array<BilinearFormIntegrator*> domain_integs;
Array<BilinearFormIntegrator*> dbfi;
/// Boundary integrators.
Array<BilinearFormIntegrator*> boundary_integs;
Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned.
Array<BilinearFormIntegrator*> bbfi;
Array<Array<int>*> bbfi_marker;///< Entries are not owned.
/// Trace face (skeleton) integrators.
Array<BilinearFormIntegrator*> trace_face_integs;
Array<BilinearFormIntegrator*> tfbfi;
/// Boundary trace face (skeleton) integrators.
Array<BilinearFormIntegrator*> boundary_trace_face_integs;
/// Entries are not owned.
Array<Array<int>*> boundary_trace_face_integs_marker;
Array<BilinearFormIntegrator*> btfbfi;
Array<Array<int>*> btfbfi_marker;///< Entries are not owned.
DenseMatrix elemmat;
Array<int> trial_vdofs, test_vdofs;
@@ -766,26 +762,24 @@ public:
Array<int> &bdr_marker);
/// Access all integrators added with AddDomainIntegrator().
Array<BilinearFormIntegrator*> *GetDBFI() { return &domain_integs; }
Array<BilinearFormIntegrator*> *GetDBFI() { return &dbfi; }
/// Access all integrators added with AddBoundaryIntegrator().
Array<BilinearFormIntegrator*> *GetBBFI() { return &boundary_integs; }
Array<BilinearFormIntegrator*> *GetBBFI() { return &bbfi; }
/** @brief Access all boundary markers added with AddBoundaryIntegrator().
If no marker was specified when the integrator was added, the
corresponding pointer (to Array<int>) will be NULL. */
Array<Array<int>*> *GetBBFI_Marker() { return &boundary_integs_marker; }
Array<Array<int>*> *GetBBFI_Marker() { return &bbfi_marker; }
/// Access all integrators added with AddTraceFaceIntegrator().
Array<BilinearFormIntegrator*> *GetTFBFI() { return &trace_face_integs; }
Array<BilinearFormIntegrator*> *GetTFBFI() { return &tfbfi; }
/// Access all integrators added with AddBdrTraceFaceIntegrator().
Array<BilinearFormIntegrator*> *GetBTFBFI()
{ return &boundary_trace_face_integs; }
Array<BilinearFormIntegrator*> *GetBTFBFI() { return &btfbfi; }
/** @brief Access all boundary markers added with AddBdrTraceFaceIntegrator().
If no marker was specified when the integrator was added, the
corresponding pointer (to Array<int>) will be NULL. */
Array<Array<int>*> *GetBTFBFI_Marker()
{ return &boundary_trace_face_integs_marker; }
Array<Array<int>*> *GetBTFBFI_Marker() { return &btfbfi_marker; }
/// Sets all sparse values of \f$ M \f$ to @a a.
void operator=(const double a) { *mat = a; }
@@ -1010,7 +1004,7 @@ public:
{ AddTraceFaceIntegrator(di); }
/// Access all interpolators added with AddDomainInterpolator().
Array<BilinearFormIntegrator*> *GetDI() { return &domain_integs; }
Array<BilinearFormIntegrator*> *GetDI() { return &dbfi; }
/// Set the desired assembly level. The default is AssemblyLevel::FULL.
/** This method must be called before assembly. */
+12 -12
View File
@@ -160,7 +160,7 @@ void MFBilinearFormExtension::Mult(const Vector &x, Vector &y) const
{
intFaceIntegrators[i]->AddMultMF(faceIntX, faceIntY);
}
int_face_restrict_lex->AddMultTranspose(faceIntY, y);
int_face_restrict_lex->MultTranspose(faceIntY, y);
}
}
@@ -176,7 +176,7 @@ void MFBilinearFormExtension::Mult(const Vector &x, Vector &y) const
{
bdrFaceIntegrators[i]->AddMultMF(faceBdrX, faceBdrY);
}
bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
bdr_face_restrict_lex->MultTranspose(faceBdrY, y);
}
}
}
@@ -217,7 +217,7 @@ void MFBilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
{
intFaceIntegrators[i]->AddMultTransposeMF(faceIntX, faceIntY);
}
int_face_restrict_lex->AddMultTranspose(faceIntY, y);
int_face_restrict_lex->MultTranspose(faceIntY, y);
}
}
@@ -233,7 +233,7 @@ void MFBilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
{
bdrFaceIntegrators[i]->AddMultTransposeMF(faceBdrX, faceBdrY);
}
bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
bdr_face_restrict_lex->MultTranspose(faceBdrY, y);
}
}
}
@@ -417,7 +417,7 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
{
intFaceIntegrators[i]->AddMultPA(faceIntX, faceIntY);
}
int_face_restrict_lex->AddMultTranspose(faceIntY, y);
int_face_restrict_lex->MultTranspose(faceIntY, y);
}
}
@@ -433,7 +433,7 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
{
bdrFaceIntegrators[i]->AddMultPA(faceBdrX, faceBdrY);
}
bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
bdr_face_restrict_lex->MultTranspose(faceBdrY, y);
}
}
}
@@ -474,7 +474,7 @@ void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
{
intFaceIntegrators[i]->AddMultTransposePA(faceIntX, faceIntY);
}
int_face_restrict_lex->AddMultTranspose(faceIntY, y);
int_face_restrict_lex->MultTranspose(faceIntY, y);
}
}
@@ -490,7 +490,7 @@ void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
{
bdrFaceIntegrators[i]->AddMultTransposePA(faceBdrX, faceBdrY);
}
bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
bdr_face_restrict_lex->MultTranspose(faceBdrY, y);
}
}
}
@@ -657,7 +657,7 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
Y(j, 0, f) += res;
});
// Apply the Interior Face Restriction transposed
int_face_restrict_lex->AddMultTranspose(faceIntY, y);
int_face_restrict_lex->MultTranspose(faceIntY, y);
}
}
@@ -688,7 +688,7 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
Y(j, f) += res;
});
// Apply the Boundary Face Restriction transposed
bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
bdr_face_restrict_lex->MultTranspose(faceBdrY, y);
}
}
}
@@ -783,7 +783,7 @@ void EABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
Y(j, 0, f) += res;
});
// Apply the Interior Face Restriction transposed
int_face_restrict_lex->AddMultTranspose(faceIntY, y);
int_face_restrict_lex->MultTranspose(faceIntY, y);
}
}
@@ -814,7 +814,7 @@ void EABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
Y(j, f) += res;
});
// Apply the Boundary Face Restriction transposed
bdr_face_restrict_lex->AddMultTranspose(faceBdrY, y);
bdr_face_restrict_lex->MultTranspose(faceBdrY, y);
}
}
}
+4 -4
View File
@@ -72,8 +72,8 @@ protected:
mutable Vector faceIntX, faceIntY;
mutable Vector faceBdrX, faceBdrY;
const Operator *elem_restrict; // Not owned
const FaceRestriction *int_face_restrict_lex; // Not owned
const FaceRestriction *bdr_face_restrict_lex; // Not owned
const Operator *int_face_restrict_lex; // Not owned
const Operator *bdr_face_restrict_lex; // Not owned
public:
PABilinearFormExtension(BilinearForm*);
@@ -143,8 +143,8 @@ protected:
mutable Vector faceIntX, faceIntY;
mutable Vector faceBdrX, faceBdrY;
const Operator *elem_restrict; // Not owned
const FaceRestriction *int_face_restrict_lex; // Not owned
const FaceRestriction *bdr_face_restrict_lex; // Not owned
const Operator *int_face_restrict_lex; // Not owned
const Operator *bdr_face_restrict_lex; // Not owned
public:
MFBilinearFormExtension(BilinearForm *form);
+2 -29
View File
@@ -175,11 +175,6 @@ void BilinearFormIntegrator::AssembleFaceVector(
elmat.Mult(elfun, elvect);
}
void TransposeIntegrator::SetIntRule(const IntegrationRule *ir)
{
IntRule = ir;
bfi->SetIntRule(ir);
}
void TransposeIntegrator::AssembleElementMatrix (
const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
@@ -207,12 +202,6 @@ void TransposeIntegrator::AssembleFaceMatrix (
elmat.Transpose (bfi_elmat);
}
void LumpedIntegrator::SetIntRule(const IntegrationRule *ir)
{
IntRule = ir;
bfi->SetIntRule(ir);
}
void LumpedIntegrator::AssembleElementMatrix (
const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
{
@@ -220,12 +209,6 @@ void LumpedIntegrator::AssembleElementMatrix (
elmat.Lump();
}
void InverseIntegrator::SetIntRule(const IntegrationRule *ir)
{
IntRule = ir;
integrator->SetIntRule(ir);
}
void InverseIntegrator::AssembleElementMatrix(
const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
{
@@ -233,15 +216,6 @@ void InverseIntegrator::AssembleElementMatrix(
elmat.Invert();
}
void SumIntegrator::SetIntRule(const IntegrationRule *ir)
{
IntRule = ir;
for (int i = 0; i < integrators.Size(); i++)
{
integrators[i]->SetIntRule(ir);
}
}
void SumIntegrator::AssembleElementMatrix(
const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
{
@@ -1777,16 +1751,15 @@ void DerivativeIntegrator::AssembleElementMatrix2 (
int dim = trial_fe.GetDim();
int trial_nd = trial_fe.GetDof();
int test_nd = test_fe.GetDof();
int spaceDim = Trans.GetSpaceDim();
int i, l;
double det;
elmat.SetSize (test_nd,trial_nd);
dshape.SetSize (trial_nd,dim);
dshapedxt.SetSize(trial_nd, spaceDim);
dshapedxt.SetSize(trial_nd,dim);
dshapedxi.SetSize(trial_nd);
invdfdx.SetSize(dim, spaceDim);
invdfdx.SetSize(dim);
shape.SetSize (test_nd);
const IntegrationRule *ir = IntRule;
-8
View File
@@ -261,8 +261,6 @@ public:
TransposeIntegrator (BilinearFormIntegrator *bfi_, int own_bfi_ = 1)
{ bfi = bfi_; own_bfi = own_bfi_; }
virtual void SetIntRule(const IntegrationRule *ir);
virtual void AssembleElementMatrix(const FiniteElement &el,
ElementTransformation &Trans,
DenseMatrix &elmat);
@@ -330,8 +328,6 @@ public:
LumpedIntegrator (BilinearFormIntegrator *bfi_, int own_bfi_ = 1)
{ bfi = bfi_; own_bfi = own_bfi_; }
virtual void SetIntRule(const IntegrationRule *ir);
virtual void AssembleElementMatrix(const FiniteElement &el,
ElementTransformation &Trans,
DenseMatrix &elmat);
@@ -350,8 +346,6 @@ public:
InverseIntegrator(BilinearFormIntegrator *integ, int own_integ = 1)
{ integrator = integ; own_integrator = own_integ; }
virtual void SetIntRule(const IntegrationRule *ir);
virtual void AssembleElementMatrix(const FiniteElement &el,
ElementTransformation &Trans,
DenseMatrix &elmat);
@@ -370,8 +364,6 @@ private:
public:
SumIntegrator(int own_integs = 1) { own_integrators = own_integs; }
virtual void SetIntRule(const IntegrationRule *ir);
void AddIntegrator(BilinearFormIntegrator *integ)
{ integrators.Append(integ); }
+1 -1
View File
@@ -143,7 +143,7 @@ Solver *BuildSmootherFromCeed(ConstrainedOperator &op, bool chebyshev)
if (chebyshev)
{
const int cheb_order = 3;
out = new OperatorChebyshevSmoother(op, t_diag, ess_tdofs, cheb_order);
out = new OperatorChebyshevSmoother(&op, t_diag, ess_tdofs, cheb_order);
}
else
{
+8 -58
View File
@@ -17,7 +17,6 @@
#include <cerrno> // errno
#include <sstream>
#include <regex>
#ifndef _WIN32
#include <sys/stat.h> // mkdir
@@ -765,8 +764,7 @@ ParaViewDataCollection::ParaViewDataCollection(const std::string&
: DataCollection(collection_name, mesh_),
levels_of_detail(1),
pv_data_format(VTKFormat::BINARY),
high_order_output(false),
restart_mode(false)
high_order_output(false)
{
#ifdef MFEM_USE_ZLIB
compression = -1; // default zlib compression level, equivalent to 6
@@ -844,60 +842,17 @@ void ParaViewDataCollection::Save()
}
// the directory is created
// create pvd file if needed. If we are not in restart mode, a new pvd file
// is always created. In restart mode, we keep any previously defined
// timestep values as long as they are less than the currently defined time.
// create pvd file if needed
if (myid == 0 && !pvd_stream.is_open())
{
std::string dpath=GenerateCollectionPath();
std::string pvdname=dpath+"/"+GeneratePVDFileName();
std::ifstream pvd_in;
if (restart_mode && (pvd_in.open(pvdname,std::ios::binary),pvd_in.good()))
{
// PVD file exists and restart mode enabled: preserve existing time
// steps less than the current time.
std::fstream::pos_type pos_begin = pvd_in.tellg();
std::fstream::pos_type pos_end = pos_begin;
std::regex regexp("timestep=\"([^[:space:]]+)\".*file=\"Cycle(\\d+)");
std::smatch match;
std::string line;
while (getline(pvd_in,line))
{
if (regex_search(line,match,regexp))
{
MFEM_ASSERT(match.size() == 3, "Unable to parse DataSet");
double tvalue = std::stod(match[1]);
if (tvalue >= GetTime()) { break; }
int cvalue = std::stoi(match[2]);
MFEM_VERIFY(cvalue < GetCycle(), "Cycle " << GetCycle() <<
" is too small for restart mode: trying to overwrite"
" existing data.");
pos_end = pvd_in.tellg();
}
}
size_t count = pos_end - pos_begin;
std::vector<char> buf(count);
pvd_in.clear();
pvd_in.seekg(pos_begin);
pvd_in.read(buf.data(), count);
pvd_in.close();
pvd_stream.open(pvdname.c_str(),std::ios::out);
pvd_stream.write(buf.data(), count);
}
else
{
// initialize new pvd file
pvd_stream.open(pvdname.c_str(),std::ios::out);
// initialize the file
pvd_stream << "<?xml version=\"1.0\"?>\n";
pvd_stream << "<VTKFile type=\"Collection\" version=\"0.1\"";
pvd_stream << " byte_order=\"" << VTKByteOrder() << "\">\n";
pvd_stream << "<Collection>" << std::endl;
}
pvd_stream.open(pvdname.c_str(),std::ios::out);
// initialize the file
pvd_stream << "<?xml version=\"1.0\"?>\n";
pvd_stream << "<VTKFile type=\"Collection\" version=\"0.1\"";
pvd_stream << " byte_order=\"" << VTKByteOrder() << "\">\n";
pvd_stream << "<Collection>" << std::endl;
}
// define the vtu file
@@ -1136,11 +1091,6 @@ void ParaViewDataCollection::SetCompression(bool compression_)
}
}
void ParaViewDataCollection::UseRestartMode(bool restart_mode_)
{
restart_mode = restart_mode_;
}
const char *ParaViewDataCollection::GetDataFormatString() const
{
if (pv_data_format == VTKFormat::ASCII)
-6
View File
@@ -488,7 +488,6 @@ private:
std::fstream pvd_stream;
VTKFormat pv_data_format;
bool high_order_output;
bool restart_mode;
protected:
void SaveDataVTU(std::ostream &out, int ref);
@@ -546,11 +545,6 @@ public:
/// by default). Reading high-order data requires ParaView 5.5 or later.
void SetHighOrderOutput(bool high_order_output_);
/// Enable or disable restart mode. If restart is enabled, new writes will
/// preserve timestep metadata for any solutions prior to the currently
/// defined time.
void UseRestartMode(bool restart_mode_);
/// Load the collection - not implemented in the ParaView writer
virtual void Load(int cycle_ = 0) override;
};
+4 -23
View File
@@ -7948,27 +7948,7 @@ VectorTensorFiniteElement::VectorTensorFiniteElement(const int dims,
p, M, FunctionSpace::Qk),
TensorBasisElement(dims, p, VerifyNodal(cbtype), dmtype),
cbasis1d(poly1d.GetBasis(p, VerifyClosed(cbtype))),
obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(obtype)))
{
MFEM_VERIFY(dims > 1, "Constructor for VectorTensorFiniteElement with both "
"open and closed bases is not valid for 1D elements.");
}
VectorTensorFiniteElement::VectorTensorFiniteElement(const int dims,
const int d,
const int p,
const int obtype,
const int M,
const DofMapType dmtype)
: VectorFiniteElement(dims, GetTensorProductGeometry(dims), d,
p, M, FunctionSpace::Pk),
TensorBasisElement(dims, p, obtype, dmtype),
cbasis1d(poly1d.GetBasis(p, VerifyOpen(obtype))),
obasis1d(poly1d.GetBasis(p, VerifyOpen(obtype)))
{
MFEM_VERIFY(dims == 1, "Constructor for VectorTensorFiniteElement without "
"closed basis is only valid for 1D elements.");
}
obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(obtype))) { }
H1_SegmentElement::H1_SegmentElement(const int p, const int btype)
: NodalTensorFiniteElement(1, p, VerifyClosed(btype), H1_DOF_MAP)
@@ -13075,8 +13055,9 @@ void ND_TriangleElement::CalcCurlShape(const IntegrationPoint &ip,
const double ND_SegmentElement::tk[1] = { 1. };
ND_SegmentElement::ND_SegmentElement(const int p, const int ob_type)
: VectorTensorFiniteElement(1, p, p - 1, ob_type, H_CURL,
DofMapType::L2_DOF_MAP),
: VectorFiniteElement(1, Geometry::SEGMENT, p, p - 1,
H_CURL, FunctionSpace::Pk),
obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type))),
dof2tk(dof)
{
if (obasis1d.IsIntegratedType()) { is_nodal = false; }
+3 -6
View File
@@ -2239,11 +2239,6 @@ public:
const int cbtype, const int obtype,
const int M, const DofMapType dmtype);
// For 1D elements: there is only an "open basis", no "closed basis"
VectorTensorFiniteElement(const int dims, const int d, const int p,
const int obtype, const int M,
const DofMapType dmtype);
const DofToQuad &GetDofToQuad(const IntegrationRule &ir,
DofToQuad::Mode mode) const;
@@ -3316,9 +3311,11 @@ public:
/// Arbitrary order Nedelec elements in 1D on a segment
class ND_SegmentElement : public VectorTensorFiniteElement
class ND_SegmentElement : public VectorFiniteElement
{
static const double tk[1];
Poly_1D::Basis &obasis1d;
Array<int> dof2tk;
public:
+2 -2
View File
@@ -1225,7 +1225,7 @@ const Operator *FiniteElementSpace::GetElementRestriction(
return L2E_nat.Ptr();
}
const FaceRestriction *FiniteElementSpace::GetFaceRestriction(
const Operator *FiniteElementSpace::GetFaceRestriction(
ElementDofOrdering e_ordering, FaceType type, L2FaceValues mul) const
{
const bool is_dg_space = IsDGSpace();
@@ -1239,7 +1239,7 @@ const FaceRestriction *FiniteElementSpace::GetFaceRestriction(
}
else
{
FaceRestriction *res;
Operator* res;
if (is_dg_space)
{
res = new L2FaceRestriction(*this, e_ordering, type, m);
+2 -2
View File
@@ -164,7 +164,7 @@ protected:
+ 8 * (int)std::get<3>(k);
}
};
using map_L2F = std::unordered_map<const key_face,FaceRestriction*,key_hash>;
using map_L2F = std::unordered_map<const key_face,Operator*,key_hash>;
mutable map_L2F L2F;
mutable Array<QuadratureInterpolator*> E2Q_array;
@@ -488,7 +488,7 @@ public:
const Operator *GetElementRestriction(ElementDofOrdering e_ordering) const;
/// Return an Operator that converts L-vectors to E-vectors on each face.
virtual const FaceRestriction *GetFaceRestriction(
virtual const Operator *GetFaceRestriction(
ElementDofOrdering e_ordering, FaceType,
L2FaceValues mul = L2FaceValues::DoubleValued) const;
+4 -9
View File
@@ -34,7 +34,6 @@ IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry)
nx = irx.GetNPoints();
ny = iry.GetNPoints();
SetSize(nx * ny);
SetPointIndices();
for (j = 0; j < ny; j++)
{
@@ -49,6 +48,8 @@ IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry)
ip.weight = ipx.weight * ipy.weight;
}
}
SetPointIndices();
}
IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry,
@@ -58,7 +59,6 @@ IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry,
const int ny = iry.GetNPoints();
const int nz = irz.GetNPoints();
SetSize(nx*ny*nz);
SetPointIndices();
for (int iz = 0; iz < nz; ++iz)
{
@@ -78,6 +78,8 @@ IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry,
}
}
}
SetPointIndices();
}
const Array<double> &IntegrationRule::GetWeights() const
@@ -123,7 +125,6 @@ void IntegrationRule::GrundmannMollerSimplexRule(int s, int n)
}
np /= f;
SetSize(np);
SetPointIndices();
int pt = 0;
for (int i = 0; i <= s; i++)
@@ -374,7 +375,6 @@ public:
void QuadratureFunctions1D::GaussLegendre(const int np, IntegrationRule* ir)
{
ir->SetSize(np);
ir->SetPointIndices();
switch (np)
{
@@ -477,7 +477,6 @@ void QuadratureFunctions1D::GaussLobatto(const int np, IntegrationRule* ir)
*/
ir->SetSize(np);
ir->SetPointIndices();
if ( np == 1 )
{
ir->IntPoint(0).Set1w(0.5, 1.0);
@@ -577,7 +576,6 @@ void QuadratureFunctions1D::GaussLobatto(const int np, IntegrationRule* ir)
void QuadratureFunctions1D::OpenUniform(const int np, IntegrationRule* ir)
{
ir->SetSize(np);
ir->SetPointIndices();
// The Newton-Cotes quadrature is based on weights that integrate exactly the
// interpolatory polynomial through the equally spaced quadrature points.
@@ -593,7 +591,6 @@ void QuadratureFunctions1D::ClosedUniform(const int np,
IntegrationRule* ir)
{
ir->SetSize(np);
ir->SetPointIndices();
if ( np == 1 ) // allow this case as "closed"
{
ir->IntPoint(0).Set1w(0.5, 1.0);
@@ -611,7 +608,6 @@ void QuadratureFunctions1D::ClosedUniform(const int np,
void QuadratureFunctions1D::OpenHalfUniform(const int np, IntegrationRule* ir)
{
ir->SetSize(np);
ir->SetPointIndices();
// Open half points: the centers of np uniform intervals
for (int i = 0; i < np ; ++i)
@@ -625,7 +621,6 @@ void QuadratureFunctions1D::OpenHalfUniform(const int np, IntegrationRule* ir)
void QuadratureFunctions1D::ClosedGL(const int np, IntegrationRule* ir)
{
ir->SetSize(np);
ir->SetPointIndices();
ir->IntPoint(0).x = 0.0;
ir->IntPoint(np-1).x = 1.0;
+3 -5
View File
@@ -96,6 +96,9 @@ private:
by request with the method GetWeights(). */
mutable Array<double> weights;
/// Sets the indices of each quadrature point on initialization.
void SetPointIndices();
/// Define n-simplex rule (triangle/tetrahedron for n=2/3) of order (2s+1)
void GrundmannMollerSimplexRule(int s, int n = 3);
@@ -224,11 +227,6 @@ public:
}
}
/// Sets the indices of each quadrature point on initialization.
/** Note that most calls to IntegrationRule::SetSize should be paired with a
call to SetPointIndices in order for the indices to be set correctly. */
void SetPointIndices();
/// Tensor product of two 1D integration rules
IntegrationRule(IntegrationRule &irx, IntegrationRule &iry);
+64 -74
View File
@@ -26,14 +26,14 @@ LinearForm::LinearForm(FiniteElementSpace *f, LinearForm *lf)
extern_lfs = 1;
// Copy the pointers to the integrators
domain_integs = lf->domain_integs;
dlfi = lf->dlfi;
domain_delta_integs = lf->domain_delta_integs;
dlfi_delta = lf->dlfi_delta;
boundary_integs = lf->boundary_integs;
blfi = lf->blfi;
boundary_face_integs = lf->boundary_face_integs;
boundary_face_integs_marker = lf->boundary_face_integs_marker;
flfi = lf->flfi;
flfi_marker = lf->flfi_marker;
}
void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi)
@@ -42,13 +42,13 @@ void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi)
dynamic_cast<DeltaLFIntegrator *>(lfi);
if (!maybe_delta || !maybe_delta->IsDelta())
{
domain_integs.Append(lfi);
dlfi.Append(lfi);
}
else
{
domain_delta_integs.Append(maybe_delta);
dlfi_delta.Append(maybe_delta);
}
domain_integs_marker.Append(NULL);
dlfi_marker.Append(NULL);
}
void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi,
@@ -58,45 +58,44 @@ void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi,
dynamic_cast<DeltaLFIntegrator *>(lfi);
if (!maybe_delta || !maybe_delta->IsDelta())
{
domain_integs.Append(lfi);
dlfi.Append(lfi);
}
else
{
domain_delta_integs.Append(maybe_delta);
dlfi_delta.Append(maybe_delta);
}
domain_integs_marker.Append(&elem_marker);
dlfi_marker.Append(&elem_marker);
}
void LinearForm::AddBoundaryIntegrator (LinearFormIntegrator * lfi)
{
boundary_integs.Append (lfi);
boundary_integs_marker.Append(NULL); // NULL -> all attributes are active
blfi.Append (lfi);
blfi_marker.Append(NULL); // NULL -> all attributes are active
}
void LinearForm::AddBoundaryIntegrator (LinearFormIntegrator * lfi,
Array<int> &bdr_attr_marker)
{
boundary_integs.Append (lfi);
boundary_integs_marker.Append(&bdr_attr_marker);
blfi.Append (lfi);
blfi_marker.Append(&bdr_attr_marker);
}
void LinearForm::AddBdrFaceIntegrator (LinearFormIntegrator * lfi)
{
boundary_face_integs.Append(lfi);
// NULL -> all attributes are active
boundary_face_integs_marker.Append(NULL);
flfi.Append(lfi);
flfi_marker.Append(NULL); // NULL -> all attributes are active
}
void LinearForm::AddBdrFaceIntegrator(LinearFormIntegrator *lfi,
Array<int> &bdr_attr_marker)
{
boundary_face_integs.Append(lfi);
boundary_face_integs_marker.Append(&bdr_attr_marker);
flfi.Append(lfi);
flfi_marker.Append(&bdr_attr_marker);
}
void LinearForm::AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
{
interior_face_integs.Append(lfi);
iflfi.Append(lfi);
}
void LinearForm::Assemble()
@@ -113,14 +112,14 @@ void LinearForm::Assemble()
// The first use of AddElementVector() below will move it back to host
// because both 'vdofs' and 'elemvect' are on host.
if (domain_integs.Size())
if (dlfi.Size())
{
for (int k = 0; k < domain_integs.Size(); k++)
for (int k = 0; k < dlfi.Size(); k++)
{
if (domain_integs_marker[k] != NULL)
if (dlfi_marker[k] != NULL)
{
MFEM_VERIFY(fes->GetMesh()->attributes.Size() ==
domain_integs_marker[k]->Size(),
dlfi_marker[k]->Size(),
"invalid element marker for domain linear form "
"integrator #" << k << ", counting from zero");
}
@@ -129,15 +128,14 @@ void LinearForm::Assemble()
for (i = 0; i < fes -> GetNE(); i++)
{
int elem_attr = fes->GetMesh()->GetAttribute(i);
for (int k = 0; k < domain_integs.Size(); k++)
for (int k = 0; k < dlfi.Size(); k++)
{
if ( domain_integs_marker[k] == NULL ||
(*(domain_integs_marker[k]))[elem_attr-1] == 1 )
if ( dlfi_marker[k] == NULL ||
(*(dlfi_marker[k]))[elem_attr-1] == 1 )
{
fes -> GetElementVDofs (i, vdofs);
eltrans = fes -> GetElementTransformation (i);
domain_integs[k]->AssembleRHSElementVect(*fes->GetFE(i),
*eltrans, elemvect);
dlfi[k]->AssembleRHSElementVect(*fes->GetFE(i), *eltrans, elemvect);
AddElementVector (vdofs, elemvect);
}
}
@@ -145,7 +143,7 @@ void LinearForm::Assemble()
}
AssembleDelta();
if (boundary_integs.Size())
if (blfi.Size())
{
Mesh *mesh = fes->GetMesh();
@@ -153,14 +151,14 @@ void LinearForm::Assemble()
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < boundary_integs.Size(); k++)
for (int k = 0; k < blfi.Size(); k++)
{
if (boundary_integs_marker[k] == NULL)
if (blfi_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *boundary_integs_marker[k];
Array<int> &bdr_marker = *blfi_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary integrator #"
<< k << ", counting from zero");
@@ -176,19 +174,18 @@ void LinearForm::Assemble()
if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
fes -> GetBdrElementVDofs (i, vdofs);
eltrans = fes -> GetBdrElementTransformation (i);
for (int k=0; k < boundary_integs.Size(); k++)
for (int k=0; k < blfi.Size(); k++)
{
if (boundary_integs_marker[k] &&
(*boundary_integs_marker[k])[bdr_attr-1] == 0) { continue; }
if (blfi_marker[k] &&
(*blfi_marker[k])[bdr_attr-1] == 0) { continue; }
boundary_integs[k]->AssembleRHSElementVect(*fes->GetBE(i),
*eltrans, elemvect);
blfi[k]->AssembleRHSElementVect(*fes->GetBE(i), *eltrans, elemvect);
AddElementVector (vdofs, elemvect);
}
}
}
if (boundary_face_integs.Size())
if (flfi.Size())
{
FaceElementTransformations *tr;
Mesh *mesh = fes->GetMesh();
@@ -197,14 +194,14 @@ void LinearForm::Assemble()
Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
mesh->bdr_attributes.Max() : 0);
bdr_attr_marker = 0;
for (int k = 0; k < boundary_face_integs.Size(); k++)
for (int k = 0; k < flfi.Size(); k++)
{
if (boundary_face_integs_marker[k] == NULL)
if (flfi_marker[k] == NULL)
{
bdr_attr_marker = 1;
break;
}
Array<int> &bdr_marker = *boundary_face_integs_marker[k];
Array<int> &bdr_marker = *flfi_marker[k];
MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
"invalid boundary marker for boundary face integrator #"
<< k << ", counting from zero");
@@ -223,26 +220,24 @@ void LinearForm::Assemble()
if (tr != NULL)
{
fes -> GetElementVDofs (tr -> Elem1No, vdofs);
for (int k = 0; k < boundary_face_integs.Size(); k++)
for (int k = 0; k < flfi.Size(); k++)
{
if (boundary_face_integs_marker[k] &&
(*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
{ continue; }
if (flfi_marker[k] &&
(*flfi_marker[k])[bdr_attr-1] == 0) { continue; }
boundary_face_integs[k]->
AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
*tr, elemvect);
flfi[k] -> AssembleRHSElementVect (*fes->GetFE(tr -> Elem1No),
*tr, elemvect);
AddElementVector (vdofs, elemvect);
}
}
}
}
if (interior_face_integs.Size())
if (iflfi.Size())
{
Mesh *mesh = fes->GetMesh();
for (int k = 0; k < interior_face_integs.Size(); k++)
for (int k = 0; k < iflfi.Size(); k++)
{
for (i = 0; i < mesh->GetNumFaces(); i++)
{
@@ -254,10 +249,9 @@ void LinearForm::Assemble()
Array<int> vdofs2;
fes -> GetElementVDofs (tr -> Elem2No, vdofs2);
vdofs.Append(vdofs2);
interior_face_integs[k]->
AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
*fes->GetFE(tr->Elem2No),
*tr, elemvect);
iflfi[k] -> AssembleRHSElementVect (*fes->GetFE(tr -> Elem1No),
*fes->GetFE(tr -> Elem2No),
*tr, elemvect);
AddElementVector (vdofs, elemvect);
}
}
@@ -283,41 +277,40 @@ void LinearForm::MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
void LinearForm::AssembleDelta()
{
if (domain_delta_integs.Size() == 0) { return; }
if (dlfi_delta.Size() == 0) { return; }
if (!HaveDeltaLocations())
{
int sdim = fes->GetMesh()->SpaceDimension();
Vector center;
DenseMatrix centers(sdim, domain_delta_integs.Size());
DenseMatrix centers(sdim, dlfi_delta.Size());
for (int i = 0; i < centers.Width(); i++)
{
centers.GetColumnReference(i, center);
domain_delta_integs[i]->GetDeltaCenter(center);
dlfi_delta[i]->GetDeltaCenter(center);
MFEM_VERIFY(center.Size() == sdim,
"Point dim " << center.Size() <<
" does not match space dim " << sdim);
}
fes->GetMesh()->FindPoints(centers, domain_delta_integs_elem_id,
domain_delta_integs_ip);
fes->GetMesh()->FindPoints(centers, dlfi_delta_elem_id, dlfi_delta_ip);
}
Array<int> vdofs;
Vector elemvect;
for (int i = 0; i < domain_delta_integs.Size(); i++)
for (int i = 0; i < dlfi_delta.Size(); i++)
{
int elem_id = domain_delta_integs_elem_id[i];
int elem_id = dlfi_delta_elem_id[i];
// The delta center may be outside of this sub-domain, or
// (Par)Mesh::FindPoints() failed to find this point:
if (elem_id < 0) { continue; }
const IntegrationPoint &ip = domain_delta_integs_ip[i];
const IntegrationPoint &ip = dlfi_delta_ip[i];
ElementTransformation &Trans = *fes->GetElementTransformation(elem_id);
Trans.SetIntPoint(&ip);
fes->GetElementVDofs(elem_id, vdofs);
domain_delta_integs[i]->AssembleDeltaElementVect(*fes->GetFE(elem_id),
Trans, elemvect);
dlfi_delta[i]->AssembleDeltaElementVect(*fes->GetFE(elem_id), Trans,
elemvect);
AddElementVector(vdofs, elemvect);
}
}
@@ -340,14 +333,11 @@ LinearForm::~LinearForm()
if (!extern_lfs)
{
int k;
for (k=0; k < domain_delta_integs.Size(); k++)
{ delete domain_delta_integs[k]; }
for (k=0; k < domain_integs.Size(); k++) { delete domain_integs[k]; }
for (k=0; k < boundary_integs.Size(); k++) { delete boundary_integs[k]; }
for (k=0; k < boundary_face_integs.Size(); k++)
{ delete boundary_face_integs[k]; }
for (k=0; k < interior_face_integs.Size(); k++)
{ delete interior_face_integs[k]; }
for (k=0; k < dlfi_delta.Size(); k++) { delete dlfi_delta[k]; }
for (k=0; k < dlfi.Size(); k++) { delete dlfi[k]; }
for (k=0; k < blfi.Size(); k++) { delete blfi[k]; }
for (k=0; k < flfi.Size(); k++) { delete flfi[k]; }
for (k=0; k < iflfi.Size(); k++) { delete iflfi[k]; }
}
}
+19 -25
View File
@@ -26,46 +26,43 @@ protected:
/// FE space on which the LinearForm lives. Not owned.
FiniteElementSpace *fes;
/** @brief Indicates the LinearFormIntegrator%s stored in #domain_integs,
#domain_delta_integs, #boundary_integs, and #boundary_face_integs are
owned by another LinearForm. */
/** @brief Indicates the LinearFormIntegrator%s stored in #dlfi, #dlfi_delta,
#blfi, and #flfi are owned by another LinearForm. */
int extern_lfs;
/// Set of Domain Integrators to be applied.
Array<LinearFormIntegrator*> domain_integs;
Array<LinearFormIntegrator*> dlfi;
/// Element attribute marker (should be of length mesh->attributes)
/// Includes all by default.
/// 0 - ignore attribute
/// 1 - include attribute
Array<Array<int>*> domain_integs_marker;
Array<Array<int>*> dlfi_marker;
/// Separate array for integrators with delta function coefficients.
Array<DeltaLFIntegrator*> domain_delta_integs;
Array<DeltaLFIntegrator*> dlfi_delta;
/// Set of Boundary Integrators to be applied.
Array<LinearFormIntegrator*> boundary_integs;
/// Entries are not owned.
Array<Array<int>*> boundary_integs_marker;
Array<LinearFormIntegrator*> blfi;
Array<Array<int>*> blfi_marker; ///< Entries are not owned.
/// Set of Boundary Face Integrators to be applied.
Array<LinearFormIntegrator*> boundary_face_integs;
Array<Array<int>*> boundary_face_integs_marker; ///< Entries not owned.
Array<LinearFormIntegrator*> flfi;
Array<Array<int>*> flfi_marker; ///< Entries are not owned.
/// Set of Internal Face Integrators to be applied.
Array<LinearFormIntegrator*> interior_face_integs;
Array<LinearFormIntegrator*> iflfi;
/// The element ids where the centers of the delta functions lie
Array<int> domain_delta_integs_elem_id;
Array<int> dlfi_delta_elem_id;
/// The reference coordinates where the centers of the delta functions lie
Array<IntegrationPoint> domain_delta_integs_ip;
Array<IntegrationPoint> dlfi_delta_ip;
/// If true, the delta locations are not (re)computed during assembly.
bool HaveDeltaLocations()
{ return (domain_delta_integs_elem_id.Size() != 0); }
bool HaveDeltaLocations() { return (dlfi_delta_elem_id.Size() != 0); }
/// Force (re)computation of delta locations.
void ResetDeltaLocations() { domain_delta_integs_elem_id.SetSize(0); }
void ResetDeltaLocations() { dlfi_delta_elem_id.SetSize(0); }
private:
/// Copy construction is not supported; body is undefined.
@@ -153,25 +150,22 @@ public:
/** @brief Access all integrators added with AddDomainIntegrator() which are
not DeltaLFIntegrator%s or they are DeltaLFIntegrator%s with non-delta
coefficients. */
Array<LinearFormIntegrator*> *GetDLFI() { return &domain_integs; }
Array<LinearFormIntegrator*> *GetDLFI() { return &dlfi; }
/** @brief Access all integrators added with AddDomainIntegrator() which are
DeltaLFIntegrator%s with delta coefficients. */
Array<DeltaLFIntegrator*> *GetDLFI_Delta() { return &domain_delta_integs; }
Array<DeltaLFIntegrator*> *GetDLFI_Delta() { return &dlfi_delta; }
/// Access all integrators added with AddBoundaryIntegrator().
Array<LinearFormIntegrator*> *GetBLFI() { return &boundary_integs; }
Array<LinearFormIntegrator*> *GetBLFI() { return &blfi; }
/// Access all integrators added with AddBdrFaceIntegrator().
Array<LinearFormIntegrator*> *GetFLFI() { return &boundary_face_integs; }
/// Access all integrators added with AddInteriorFaceIntegrator().
Array<LinearFormIntegrator*> *GetIFLFI() { return &interior_face_integs; }
Array<LinearFormIntegrator*> *GetFLFI() { return &flfi; }
/** @brief Access all boundary markers added with AddBdrFaceIntegrator().
If no marker was specified when the integrator was added, the
corresponding pointer (to Array<int>) will be NULL. */
Array<Array<int>*> *GetFLFI_Marker() { return &boundary_face_integs_marker; }
Array<Array<int>*> *GetFLFI_Marker() { return &flfi_marker; }
/// Assembles the linear form i.e. sums over all domain/bdr integrators.
void Assemble();
-2
View File
@@ -979,8 +979,6 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
for (int k = 0; k < bfnfi.Size(); ++k)
{
if (bfnfi_marker[k] &&
(*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }
bfnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
for (int l=0; l<fes.Size(); ++l)
{
+2 -2
View File
@@ -40,10 +40,10 @@ protected:
public:
/** @brief Prescribe a fixed IntegrationRule to use (when @a ir != NULL) or
let the integrator choose (when @a ir == NULL). */
virtual void SetIntRule(const IntegrationRule *ir) { IntRule = ir; }
void SetIntRule(const IntegrationRule *ir) { IntRule = ir; }
/// Prescribe a fixed IntegrationRule to use.
void SetIntegrationRule(const IntegrationRule &ir) { SetIntRule(&ir); }
void SetIntegrationRule(const IntegrationRule &irule) { IntRule = &irule; }
/// Set the memory type used for GeometricFactors and other large allocations
/// in PA extensions.
+8 -10
View File
@@ -130,7 +130,7 @@ void ParBilinearForm::ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local)
OperatorHandle dA(A.Type()), Ph(A.Type()), hdA;
if (interior_face_integs.Size() == 0)
if (fbfi.Size() == 0)
{
// construct a parallel block-diagonal matrix 'A' based on 'a'
dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
@@ -214,12 +214,11 @@ void ParBilinearForm::AssembleSharedFaces(int skip_zeros)
}
}
vdofs_all.Append(vdofs2);
for (int k = 0; k < interior_face_integs.Size(); k++)
for (int k = 0; k < fbfi.Size(); k++)
{
interior_face_integs[k]->
AssembleFaceMatrix(*pfes->GetFE(T->Elem1No),
*pfes->GetFaceNbrFE(Elem2NbrNo),
*T, elemmat);
fbfi[k]->AssembleFaceMatrix(*pfes->GetFE(T->Elem1No),
*pfes->GetFaceNbrFE(Elem2NbrNo),
*T, elemmat);
if (keep_nbr_block)
{
mat->AddSubMatrix(vdofs_all, vdofs_all, elemmat, skip_zeros);
@@ -234,7 +233,7 @@ void ParBilinearForm::AssembleSharedFaces(int skip_zeros)
void ParBilinearForm::Assemble(int skip_zeros)
{
if (interior_face_integs.Size())
if (fbfi.Size())
{
pfes->ExchangeFaceNbrData();
if (!ext && mat == NULL)
@@ -245,7 +244,7 @@ void ParBilinearForm::Assemble(int skip_zeros)
BilinearForm::Assemble(skip_zeros);
if (!ext && interior_face_integs.Size() > 0)
if (!ext && fbfi.Size() > 0)
{
AssembleSharedFaces(skip_zeros);
}
@@ -317,8 +316,7 @@ ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const double a)
const
{
MFEM_VERIFY(interior_face_integs.Size() == 0,
"the case of interior face integrators is not"
MFEM_VERIFY(fbfi.Size() == 0, "the case of interior face integrators is not"
" implemented");
if (X.ParFESpace() != pfes)
+2 -2
View File
@@ -515,7 +515,7 @@ const FiniteElement *ParFiniteElementSpace::GetFE(int i) const
else { return FiniteElementSpace::GetFE(i); }
}
const FaceRestriction *ParFiniteElementSpace::GetFaceRestriction(
const Operator *ParFiniteElementSpace::GetFaceRestriction(
ElementDofOrdering e_ordering, FaceType type, L2FaceValues mul) const
{
const bool is_dg_space = IsDGSpace();
@@ -529,7 +529,7 @@ const FaceRestriction *ParFiniteElementSpace::GetFaceRestriction(
}
else
{
FaceRestriction *res;
Operator* res;
if (is_dg_space)
{
res = new ParL2FaceRestriction(*this, e_ordering, type, m);
+1 -1
View File
@@ -299,7 +299,7 @@ public:
presence of shared faces. Shared faces are treated as interior faces,
the returned operator handles the communication needed to get the
shared face values from other MPI ranks */
virtual const FaceRestriction *GetFaceRestriction(
virtual const Operator *GetFaceRestriction(
ElementDofOrdering e_ordering, FaceType type,
L2FaceValues mul = L2FaceValues::DoubleValued) const;
+6 -7
View File
@@ -47,7 +47,7 @@ void ParLinearForm::Assemble()
{
LinearForm::Assemble();
if (interior_face_integs.Size())
if (iflfi.Size())
{
pfes->ExchangeFaceNbrData();
AssembleSharedFaces();
@@ -59,10 +59,10 @@ void ParLinearForm::AssembleSharedFaces()
Array<int> vdofs;
Vector elemvect;
if (interior_face_integs.Size())
if (iflfi.Size())
{
ParMesh *pmesh = pfes->GetParMesh();
for (int k = 0; k < interior_face_integs.Size(); k++)
for (int k = 0; k < iflfi.Size(); k++)
{
for (int i = 0; i < pmesh->GetNSharedFaces(); i++)
{
@@ -73,10 +73,9 @@ void ParLinearForm::AssembleSharedFaces()
{
int Elem2Nbr = tr->Elem2No - pmesh->GetNE();
fes -> GetElementVDofs (tr -> Elem1No, vdofs);
interior_face_integs[k]->
AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
*pfes->GetFaceNbrFE(Elem2Nbr),
*tr, elemvect);
iflfi[0] -> AssembleRHSElementVect (*fes->GetFE(tr -> Elem1No),
*pfes->GetFaceNbrFE(Elem2Nbr),
*tr, elemvect);
AddElementVector (vdofs, elemvect);
}
}
+5 -5
View File
@@ -847,7 +847,7 @@ void H1FaceRestriction::Mult(const Vector& x, Vector& y) const
});
}
void H1FaceRestriction::AddMultTranspose(const Vector& x, Vector& y) const
void H1FaceRestriction::MultTranspose(const Vector& x, Vector& y) const
{
// Assumes all elements have the same number of dofs
const int nd = dof;
@@ -856,7 +856,7 @@ void H1FaceRestriction::AddMultTranspose(const Vector& x, Vector& y) const
auto d_offsets = offsets.Read();
auto d_indices = gather_indices.Read();
auto d_x = Reshape(x.Read(), nd, vd, nf);
auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
MFEM_FORALL(i, ndofs,
{
const int offset = d_offsets[i];
@@ -1267,7 +1267,7 @@ void L2FaceRestriction::Mult(const Vector& x, Vector& y) const
}
}
void L2FaceRestriction::AddMultTranspose(const Vector& x, Vector& y) const
void L2FaceRestriction::MultTranspose(const Vector& x, Vector& y) const
{
// Assumes all elements have the same number of dofs
const int nd = dof;
@@ -1280,7 +1280,7 @@ void L2FaceRestriction::AddMultTranspose(const Vector& x, Vector& y) const
if (m == L2FaceValues::DoubleValued)
{
auto d_x = Reshape(x.Read(), nd, vd, 2, nf);
auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
MFEM_FORALL(i, ndofs,
{
const int offset = d_offsets[i];
@@ -1304,7 +1304,7 @@ void L2FaceRestriction::AddMultTranspose(const Vector& x, Vector& y) const
else
{
auto d_x = Reshape(x.Read(), nd, vd, nf);
auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
MFEM_FORALL(i, ndofs,
{
const int offset = d_offsets[i];
+15 -121
View File
@@ -21,6 +21,10 @@ namespace mfem
class FiniteElementSpace;
enum class ElementDofOrdering;
/** An enum type to specify if only e1 value is requested (SingleValued) or both
e1 and e2 (DoubleValued). */
enum class L2FaceValues : bool {SingleValued, DoubleValued};
/// Operator that converts FiniteElementSpace L-vectors to E-vectors.
/** Objects of this type are typically created and owned by FiniteElementSpace
objects, see FiniteElementSpace::GetElementRestriction(). */
@@ -100,75 +104,10 @@ public:
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const;
};
/** An enum type to specify if only e1 value is requested (SingleValued) or both
e1 and e2 (DoubleValued). */
enum class L2FaceValues : bool {SingleValued, DoubleValued};
/** @brief Base class for operators that extracts Face degrees of freedom.
In order to compute quantities on the faces of a mesh, it is often useful to
extract the degrees of freedom on the faces of the elements. This class
provides an interface for such operations.
If the FiniteElementSpace is ordered by Ordering::byVDIM, then the expected
format for the L-vector is (vdim x ndofs), otherwise if Ordering::byNODES
the expected format is (ndofs x vdim), where ndofs is the total number of
degrees of freedom.
Since FiniteElementSpace can either be continuous or discontinuous, the
degrees of freedom on a face can either be single valued or double valued,
this is what we refer to as the multiplicity and is represented by the
L2FaceValues enum type.
The format of the output face E-vector of degrees of freedom is
(face_dofs x vdim x multiplicity x nfaces), where face_dofs is the number of
degrees of freedom on each face, and nfaces the number of faces of the
requested FaceType (see FiniteElementSpace::GetNFbyType).
@note Objects of this type are typically created and owned by
FiniteElementSpace objects, see FiniteElementSpace::GetFaceRestriction(). */
class FaceRestriction : public Operator
{
public:
FaceRestriction(): Operator() { }
FaceRestriction(int h, int w): Operator(h, w) { }
virtual ~FaceRestriction() { }
/** @brief Extract the face degrees of freedom from @a x into @a y.
@param[in] x The L-vector of degrees of freedom.
@param[out] y The degrees of freedom on the face, corresponding to a face
E-vector.
*/
void Mult(const Vector &x, Vector &y) const override = 0;
/** @brief Add the face degrees of freedom @a x to the element degrees of
freedom @a y.
@param[in] x The face degrees of freedom on the face.
@param[in,out] y The L-vector of degrees of freedom to which we add the
face degrees of freedom.
*/
virtual void AddMultTranspose(const Vector &x, Vector &y) const = 0;
/** @brief Set the face degrees of freedom in the element degrees of freedom
@a y to the values given in @a x.
@param[in] x The face degrees of freedom on the face.
@param[in,out] y The L-vector of degrees of freedom to which we add the
face degrees of freedom.
*/
void MultTranspose(const Vector &x, Vector &y) const override
{
y = 0.0;
AddMultTranspose(x, y);
}
};
/// Operator that extracts Face degrees of freedom for H1 FiniteElementSpaces.
/// Operator that extracts Face degrees of freedom.
/** Objects of this type are typically created and owned by FiniteElementSpace
objects, see FiniteElementSpace::GetFaceRestriction(). */
class H1FaceRestriction : public FaceRestriction
class H1FaceRestriction : public Operator
{
protected:
const FiniteElementSpace &fes;
@@ -183,42 +122,16 @@ protected:
Array<int> gather_indices;
public:
/** @brief Constructor for a H1FaceRestriction.
@param[in] fes The FiniteElementSpace on which this H1FaceRestriction
operates.
@param[in] ordering The requested output ordering of the
H1FaceRestriction, either Native or Lexicographic.
@param[in] type The requested type of faces on which this operator
extracts the degrees of freedom, either Interior or
Boundary.
*/
H1FaceRestriction(const FiniteElementSpace& fes,
const ElementDofOrdering ordering,
const FaceType type);
/** @brief Extract the face degrees of freedom from @a x into @a y.
@param[in] x The L-vector of degrees of freedom.
@param[out] y The degrees of freedom on the face, corresponding to a face
E-vector.
*/
void Mult(const Vector &x, Vector &y) const override;
/** @brief Add the face degrees of freedom @a x to the element degrees of
freedom @a y.
@param[in] x The face degrees of freedom on the face.
@param[in,out] y The L-vector of degrees of freedom to which we add the
face degrees of freedom.
*/
void AddMultTranspose(const Vector &x, Vector &y) const override;
H1FaceRestriction(const FiniteElementSpace&, const ElementDofOrdering,
const FaceType);
void Mult(const Vector &x, Vector &y) const;
void MultTranspose(const Vector &x, Vector &y) const;
};
/// Operator that extracts Face degrees of freedom on L2 FiniteElementSpaces.
/// Operator that extracts Face degrees of freedom.
/** Objects of this type are typically created and owned by FiniteElementSpace
objects, see FiniteElementSpace::GetFaceRestriction(). */
class L2FaceRestriction : public FaceRestriction
class L2FaceRestriction : public Operator
{
protected:
const FiniteElementSpace &fes;
@@ -241,38 +154,19 @@ protected:
const L2FaceValues m = L2FaceValues::DoubleValued);
public:
L2FaceRestriction(const FiniteElementSpace&,
const ElementDofOrdering,
L2FaceRestriction(const FiniteElementSpace&, const ElementDofOrdering,
const FaceType,
const L2FaceValues m = L2FaceValues::DoubleValued);
/** @brief Extract the face degrees of freedom from @a x into @a y.
@param[in] x The L-vector of degrees of freedom.
@param[out] y The degrees of freedom on the face, corresponding to a face
E-vector.
*/
void Mult(const Vector &x, Vector &y) const override;
/** @brief Add the face degrees of freedom @a x to the element degrees of
freedom @a y.
@param[in] x The face degrees of freedom on the face.
@param[in,out] y The L-vector of degrees of freedom to which we add the
face degrees of freedom.
*/
void AddMultTranspose(const Vector &x, Vector &y) const override;
virtual void Mult(const Vector &x, Vector &y) const;
void MultTranspose(const Vector &x, Vector &y) const;
/** Fill the I array of SparseMatrix corresponding to the sparsity pattern
given by this L2FaceRestriction. */
virtual void FillI(SparseMatrix &mat, const bool keep_nbr_block = false) const;
/** Fill the J and Data arrays of SparseMatrix corresponding to the sparsity
pattern given by this L2FaceRestriction, and the values of ea_data. */
virtual void FillJAndData(const Vector &ea_data,
SparseMatrix &mat,
const bool keep_nbr_block = false) const;
/// This methods adds the DG face matrices to the element matrices.
void AddFaceMatricesToElementMatrices(Vector &fea_data,
Vector &ea_data) const;
-2
View File
@@ -34,7 +34,6 @@ void TMOP_Combo_QualityMetric::EvalP(const DenseMatrix &Jpt,
DenseMatrix &P) const
{
DenseMatrix Pt(P.Size());
P = 0.0;
for (int i = 0; i < tmop_q_arr.Size(); i++)
{
tmop_q_arr[i]->EvalP(Jpt, Pt);
@@ -51,7 +50,6 @@ void TMOP_Combo_QualityMetric::AssembleH(const DenseMatrix &Jpt,
DenseMatrix At(A.Size());
for (int i = 0; i < tmop_q_arr.Size(); i++)
{
At = 0.0;
tmop_q_arr[i]->AssembleH(Jpt, DS, weight, At);
At *= wt_arr[i];
A += At;
-48
View File
@@ -371,8 +371,6 @@ public:
AddQualityMetric(sh_metric, 1.-gamma_);
AddQualityMetric(sz_metric, gamma_);
}
virtual int Id() const { return 80; }
double GetGamma() const { return gamma; }
virtual ~TMOP_Metric_080() { delete sh_metric; delete sz_metric; }
};
@@ -592,52 +590,6 @@ public:
virtual int Id() const { return 321; }
};
/// 3D barrier Shape+Size (VS) metric (polyconvex).
class TMOP_Metric_332 : public TMOP_Combo_QualityMetric
{
protected:
double gamma;
TMOP_QualityMetric *sh_metric, *sz_metric;
public:
TMOP_Metric_332(double gamma_) : gamma(gamma_),
sh_metric(new TMOP_Metric_302),
sz_metric(new TMOP_Metric_315)
{
// (1-gamma) mu_302 + gamma mu_315
AddQualityMetric(sh_metric, 1.-gamma_);
AddQualityMetric(sz_metric, gamma_);
}
virtual int Id() const { return 332; }
double GetGamma() const { return gamma; }
virtual ~TMOP_Metric_332() { delete sh_metric; delete sz_metric; }
};
/// 3D barrier Shape+Size (VS) metric (polyconvex).
class TMOP_Metric_333 : public TMOP_Combo_QualityMetric
{
protected:
double gamma;
TMOP_QualityMetric *sh_metric, *sz_metric;
public:
TMOP_Metric_333(double gamma_) : gamma(gamma_),
sh_metric(new TMOP_Metric_302),
sz_metric(new TMOP_Metric_316)
{
// (1-gamma) mu_302 + gamma mu_316
AddQualityMetric(sh_metric, 1.-gamma_);
AddQualityMetric(sz_metric, gamma_);
}
virtual int Id() const { return 333; }
double GetGamma() const { return gamma; }
virtual ~TMOP_Metric_333() { delete sh_metric; delete sz_metric; }
};
/// Shifted barrier form of 3D metric 16 (volume, ideal barrier metric), 3D
class TMOP_Metric_352 : public TMOP_QualityMetric
{
+2 -46
View File
@@ -150,49 +150,9 @@ void EvalH_077(const int e, const int qx, const int qy,
}
}
static MFEM_HOST_DEVICE inline
void EvalH_080(const int e, const int qx, const int qy,
const double weight, const double gamma, const double *Jpt,
DeviceTensor<7,double> H)
{
// h_80 = (1-gamma) h_2 + gamma h_77.
constexpr int DIM = 2;
double ddI1[4], ddI1b[4], dI2[4], dI2b[4], ddI2[4];
kernels::InvariantsEvaluator2D ie(Args()
.J(Jpt)
.dI2(dI2)
.ddI1(ddI1)
.ddI1b(ddI1b)
.dI2b(dI2b)
.ddI2(ddI2));
const double I2 = ie.Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
ConstDeviceMatrix di2(ie.Get_dI2(),DIM,DIM);
for (int i = 0; i < DIM; i++)
{
for (int j = 0; j < DIM; j++)
{
ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
ConstDeviceMatrix ddi2(ie.Get_ddI2(i,j),DIM,DIM);
for (int r = 0; r < DIM; r++)
{
for (int c = 0; c < DIM; c++)
{
H(r,c,i,j,qx,qy,e) =
(1.0 - gamma) * 0.5 * weight * ddi1b(r,c) +
gamma * ( weight * 0.5 * (1.0 - I2inv_sq) * ddi2(r,c) +
weight * (I2inv_sq / I2) * di2(r,c) * di2(i,j) );
}
}
}
}
}
MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_2D,
const Vector &x_,
const double metric_normal,
const double metric_param,
const int mid,
const int NE,
const Array<double> &w_,
@@ -203,7 +163,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_2D,
const int d1d,
const int q1d)
{
MFEM_VERIFY(mid == 1 || mid == 2 || mid == 7 || mid == 77 || mid == 80,
MFEM_VERIFY(mid == 1 || mid == 2 || mid == 7 || mid == 77,
"Metric not yet implemented!");
constexpr int DIM = 2;
@@ -262,7 +222,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_2D,
if (mid == 2) { EvalH_002(e,qx,qy,weight,Jpt,H); }
if (mid == 7) { EvalH_007(e,qx,qy,weight,Jpt,H); }
if (mid == 77) { EvalH_077(e,qx,qy,weight,Jpt,H); }
if (mid == 80) { EvalH_080(e,qx,qy,weight,metric_param,Jpt,H); }
} // qx
} // qy
});
@@ -282,10 +241,7 @@ void TMOP_Integrator::AssembleGradPA_2D(const Vector &X) const
const Array<double> &G = PA.maps->G;
Vector &H = PA.H;
double mp = 0.0;
if (auto m = dynamic_cast<TMOP_Metric_080 *>(metric)) { mp = m->GetGamma(); }
MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_2D,id,X,mn,mp,M,N,W,B,G,J,H);
MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_2D,id,X,mn,M,N,W,B,G,J,H);
}
} // namespace mfem
+3 -57
View File
@@ -181,58 +181,8 @@ void EvalH_321(const int e, const int qx, const int qy, const int qz,
}
}
// H_332 = (1-gamma) H_302 + gamma H_315
static MFEM_HOST_DEVICE inline
void EvalH_332(const int e, const int qx, const int qy, const int qz,
const double weight, const double gamma,
const double *J, DeviceTensor<8,double> dP)
{
double B[9];
double dI1b[9], ddI1b[9];
double dI2[9], dI2b[9], ddI2[9], ddI2b[9];
double dI3b[9], ddI3b[9];
constexpr int DIM = 3;
kernels::InvariantsEvaluator3D ie(Args()
.J(J).B(B)
.dI1b(dI1b).ddI1b(ddI1b)
.dI2(dI2).dI2b(dI2b).ddI2(ddI2).ddI2b(ddI2b)
.dI3b(dI3b).ddI3b(ddI3b));
double sign_detJ;
const double c1 = weight/9.;
const double I1b = ie.Get_I1b();
const double I2b = ie.Get_I2b();
const double I3b = ie.Get_I3b(sign_detJ);
ConstDeviceMatrix di1b(ie.Get_dI1b(),DIM,DIM);
ConstDeviceMatrix di2b(ie.Get_dI2b(),DIM,DIM);
ConstDeviceMatrix di3b(ie.Get_dI3b(sign_detJ),DIM,DIM);
for (int i = 0; i < DIM; i++)
{
for (int j = 0; j < DIM; j++)
{
ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
ConstDeviceMatrix ddi2b(ie.Get_ddI2b(i,j),DIM,DIM);
ConstDeviceMatrix ddi3b(ie.Get_ddI3b(i,j),DIM,DIM);
for (int r = 0; r < DIM; r++)
{
for (int c = 0; c < DIM; c++)
{
const double dp_302 =
(di2b(r,c)*di1b(i,j) + di1b(r,c)*di2b(i,j))
+ ddi2b(r,c)*I1b
+ ddi1b(r,c)*I2b;
const double dp_315 = 2.0 * weight * (I3b - 1.0) * ddi3b(r,c) +
2.0 * weight * di3b(r,c) * di3b(i,j);
dP(r,c,i,j,qx,qy,qz,e) = (1.0 - gamma) * c1 * dp_302 +
gamma * dp_315;
}
}
}
}
}
MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_3D,
const double metric_normal,
const double metric_param,
const int mid,
const Vector &x_,
const int NE,
@@ -244,8 +194,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_3D,
const int d1d,
const int q1d)
{
MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 ||
mid == 321 || mid == 332, "3D metric not yet implemented!");
MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 || mid == 321 ,
"3D metric not yet implemented!");
constexpr int DIM = 3;
const int D1D = T_D1D ? T_D1D : d1d;
@@ -305,7 +255,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_3D,
if (mid == 303) { EvalH_303(e,qx,qy,qz,weight,Jpt,H); }
if (mid == 315) { EvalH_315(e,qx,qy,qz,weight,Jpt,H); }
if (mid == 321) { EvalH_321(e,qx,qy,qz,weight,Jpt,H); }
if (mid == 332) { EvalH_332(e,qx,qy,qz,weight,metric_param,Jpt,H); }
} // qx
} // qy
} // qz
@@ -326,10 +275,7 @@ void TMOP_Integrator::AssembleGradPA_3D(const Vector &X) const
const Array<double> &G = PA.maps->G;
Vector &H = PA.H;
double mp = 0.0;
if (auto m = dynamic_cast<TMOP_Metric_332 *>(metric)) { mp = m->GetGamma(); }
MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_3D,id,mn,mp,M,X,N,W,B,G,J,H);
MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_3D,id,mn,M,X,N,W,B,G,J,H);
}
} // namespace mfem
+2 -22
View File
@@ -58,24 +58,8 @@ void EvalP_077(const double *Jpt, double *P)
kernels::Set(2,2, 0.5 * (1.0 - 1.0 / (I2 * I2)), ie.Get_dI2(), P);
}
static MFEM_HOST_DEVICE inline
void EvalP_080(const double *Jpt, double gamma, double *P)
{
// p_80 = (1-gamma) p_2 + gamma p_77.
double dI1b[4], dI2[4], dI2b[4];
kernels::InvariantsEvaluator2D ie(Args().J(Jpt).
dI1b(dI1b).dI2(dI2).dI2b(dI2b));
kernels::Set(2,2, (1.0 - gamma) * 1./2., ie.Get_dI1b(), P);
const double I2 = ie.Get_I2();
kernels::Add(2,2, gamma * 0.5 * (1.0 - 1.0 / (I2 * I2)), ie.Get_dI2(), P);
}
MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_2D,
const double metric_normal,
const double metric_param,
const int mid,
const int NE,
const DenseTensor &j_,
@@ -87,7 +71,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_2D,
const int d1d,
const int q1d)
{
MFEM_VERIFY(mid == 1 || mid == 2 || mid == 7 || mid == 77 || mid == 80,
MFEM_VERIFY(mid == 1 || mid == 2 || mid == 7 || mid == 77,
"Metric not yet implemented!");
constexpr int DIM = 2;
@@ -148,7 +132,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_2D,
if (mid == 2) { EvalP_002(Jpt, P); }
if (mid == 7) { EvalP_007(Jpt, P); }
if (mid == 77) { EvalP_077(Jpt, P); }
if (mid == 80) { EvalP_080(Jpt, metric_param, P); }
for (int i = 0; i < 4; i++) { P[i] *= weight; }
// PMatO += DS . P^t += DSh . (Jrt . P^t)
@@ -177,10 +160,7 @@ void TMOP_Integrator::AddMultPA_2D(const Vector &X, Vector &Y) const
const Array<double> &G = PA.maps->G;
const double mn = metric_normal;
double mp = 0.0;
if (auto m = dynamic_cast<TMOP_Metric_080 *>(metric)) { mp = m->GetGamma(); }
MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_2D,id,mn,mp,M,N,J,W,B,G,X,Y);
MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_2D,id,mn,M,N,J,W,B,G,X,Y);
}
} // namespace mfem
+7 -32
View File
@@ -75,29 +75,8 @@ void EvalP_321(const double *J, double *P)
kernels::Add(3,3, ie.Get_dI1(), P);
}
// P_332 = (1-gamma) P_302 + gamma P_315.
static MFEM_HOST_DEVICE inline
void EvalP_332(const double *J, double gamma, double *P)
{
double B[9];
double dI1b[9], dI2[9], dI2b[9], dI3b[9];
kernels::InvariantsEvaluator3D ie(Args()
.J(J).B(B)
.dI1b(dI1b)
.dI2(dI2).dI2b(dI2b)
.dI3b(dI3b));
const double alpha = (1.0 - gamma) * ie.Get_I1b()/9.;
const double beta = (1.0 - gamma) * ie.Get_I2b()/9.;
kernels::Add(3,3, alpha, ie.Get_dI2b(), beta, ie.Get_dI1b(), P);
double sign_detJ;
const double I3b = ie.Get_I3b(sign_detJ);
kernels::Add(3,3, gamma * 2.0 * (I3b - 1.0), ie.Get_dI3b(sign_detJ), P);
}
MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_3D,
const double metric_normal,
double metric_param,
const int mid,
const int NE,
const DenseTensor &j_,
@@ -109,8 +88,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_3D,
const int d1d,
const int q1d)
{
MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 ||
mid == 321 || mid == 332, "3D metric not yet implemented!");
MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 || mid == 321 ,
"3D metric not yet implemented!");
constexpr int DIM = 3;
const int D1D = T_D1D ? T_D1D : d1d;
@@ -167,11 +146,10 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_3D,
// metric->EvalP(Jpt, P);
double P[9];
if (mid == 302) { EvalP_302(Jpt, P); }
if (mid == 303) { EvalP_303(Jpt, P); }
if (mid == 315) { EvalP_315(Jpt, P); }
if (mid == 321) { EvalP_321(Jpt, P); }
if (mid == 332) { EvalP_332(Jpt, metric_param, P); }
if (mid == 302) { EvalP_302(Jpt,P); }
if (mid == 303) { EvalP_303(Jpt,P); }
if (mid == 315) { EvalP_315(Jpt,P); }
if (mid == 321) { EvalP_321(Jpt,P); }
for (int i = 0; i < 9; i++) { P[i] *= weight; }
// Y += DS . P^t += DSh . (Jrt . P^t)
@@ -202,10 +180,7 @@ void TMOP_Integrator::AddMultPA_3D(const Vector &X, Vector &Y) const
const Array<double> &G = PA.maps->G;
const double mn = metric_normal;
double mp = 0.0;
if (auto m = dynamic_cast<TMOP_Metric_332 *>(metric)) { mp = m->GetGamma(); }
MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_3D,id,mn,mp,M,N,J,W,B,G,X,Y);
MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_3D,id,mn,M,N,J,W,B,G,X,Y);
}
} // namespace mfem
+4 -15
View File
@@ -50,15 +50,8 @@ double EvalW_077(const double *Jpt)
return 0.5*(I2b*I2b + 1./(I2b*I2b) - 2.);
}
static MFEM_HOST_DEVICE inline
double EvalW_080(const double *Jpt, double gamma)
{
return (1.0 - gamma) * EvalW_002(Jpt) + gamma * EvalW_077(Jpt);
}
MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_2D,
const double metric_normal,
const double metric_param,
const int mid,
const int NE,
const DenseTensor &j_,
@@ -71,7 +64,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_2D,
const int d1d,
const int q1d)
{
MFEM_VERIFY(mid == 1 || mid == 2 || mid == 7 || mid == 77 || mid == 80,
MFEM_VERIFY(mid == 1 || mid == 2 || mid == 7 || mid == 77,
"2D metric not yet implemented!");
constexpr int DIM = 2;
@@ -132,8 +125,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_2D,
mid == 1 ? EvalW_001(Jpt) :
mid == 2 ? EvalW_002(Jpt) :
mid == 7 ? EvalW_007(Jpt) :
mid == 77 ? EvalW_077(Jpt) :
mid == 80 ? EvalW_080(Jpt, metric_param) : 0.0;
mid == 77 ? EvalW_077(Jpt) : 0.0;
E(qx,qy,e) = weight * EvalW;
}
@@ -149,7 +141,7 @@ double TMOP_Integrator::GetLocalStateEnergyPA_2D(const Vector &X) const
const int D1D = PA.maps->ndof;
const int Q1D = PA.maps->nqpt;
const int id = (D1D << 4 ) | Q1D;
const double mn = metric_normal;
const double m = metric_normal;
const DenseTensor &J = PA.Jtr;
const Array<double> &W = PA.ir->GetWeights();
const Array<double> &B = PA.maps->B;
@@ -157,10 +149,7 @@ double TMOP_Integrator::GetLocalStateEnergyPA_2D(const Vector &X) const
const Vector &O = PA.O;
Vector &E = PA.E;
double mp = 0.0;
if (auto m = dynamic_cast<TMOP_Metric_080 *>(metric)) { mp = m->GetGamma(); }
MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_2D,id,mn,mp,M,N,J,W,B,G,X,O,E);
MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_2D,id,m,M,N,J,W,B,G,X,O,E);
}
} // namespace mfem
+4 -15
View File
@@ -58,15 +58,8 @@ double EvalW_321(const double *J)
return ie.Get_I1() + ie.Get_I2()/ie.Get_I3() - 6.0;
}
static MFEM_HOST_DEVICE inline
double EvalW_332(const double *J, double gamma)
{
return (1.0 - gamma) * EvalW_302(J) + gamma * EvalW_315(J);
}
MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_3D,
const double metric_normal,
const double metric_param,
const int mid,
const int NE,
const DenseTensor &j_,
@@ -79,8 +72,8 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_3D,
const int d1d,
const int q1d)
{
MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 ||
mid == 321 || mid == 332, "3D metric not yet implemented!");
MFEM_VERIFY(mid == 302 || mid == 303 || mid == 315 || mid == 321 ,
"3D metric not yet implemented!");
constexpr int DIM = 3;
const int D1D = T_D1D ? T_D1D : d1d;
@@ -141,8 +134,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_3D,
mid == 302 ? EvalW_302(Jpt) :
mid == 303 ? EvalW_303(Jpt) :
mid == 315 ? EvalW_315(Jpt) :
mid == 321 ? EvalW_321(Jpt) :
mid == 332 ? EvalW_332(Jpt, metric_param) : 0.0;
mid == 321 ? EvalW_321(Jpt) : 0.0;
E(qx,qy,qz,e) = weight * EvalW;
}
@@ -167,10 +159,7 @@ double TMOP_Integrator::GetLocalStateEnergyPA_3D(const Vector &X) const
const Vector &O = PA.O;
Vector &E = PA.E;
double mp = 0.0;
if (auto m = dynamic_cast<TMOP_Metric_332 *>(metric)) { mp = m->GetGamma(); }
MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_3D,id,mn,mp,M,N,J,W,B,G,O,X,E);
MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_3D,id,mn,M,N,J,W,B,G,O,X,E);
}
} // namespace mfem
+3 -5
View File
@@ -71,7 +71,7 @@ static const unsigned char b64table[] =
255,255,255,255,255,255,255,255,255,255,255,255,255,255,255,255
};
void DecodeBase64(const char *src, size_t len, std::vector<char> &buf)
void DecodeBase64(const char *src, size_t len, std::vector<unsigned char> &buf)
{
const unsigned char *in = (const unsigned char *)src;
buf.clear();
@@ -79,7 +79,7 @@ void DecodeBase64(const char *src, size_t len, std::vector<char> &buf)
for (size_t i=0; i<len; ++i) { if (b64table[in[i]] != 255) { ++count; } }
if (count % 4 != 0) { return; }
buf.resize(3*len/4);
unsigned char *out = (unsigned char *)buf.data();
unsigned char *out = buf.data();
count = 0;
int pad = 0;
unsigned char c[4];
@@ -97,10 +97,8 @@ void DecodeBase64(const char *src, size_t len, std::vector<char> &buf)
count = pad = 0;
}
}
buf.resize(out - (unsigned char *)buf.data());
buf.resize(out - buf.data());
}
size_t NumBase64Chars(size_t nbytes) { return ((4*nbytes/3) + 3) & ~3; }
} // namespace mfem::bin_io
} // namespace mfem
+2 -11
View File
@@ -50,7 +50,6 @@ inline T read(const char *buf)
return value;
}
/// Append the binary representation of @a val to the byte buffer @a vec.
template <typename T>
void AppendBytes(std::vector<char> &vec, const T &val)
{
@@ -58,17 +57,9 @@ void AppendBytes(std::vector<char> &vec, const T &val)
vec.insert(vec.end(), ptr, ptr + sizeof(T));
}
/// Given a buffer @a buf of length @a nbytes, encode the data in base-64
/// format, and write the encoded data to the output stream @a out.
void WriteBase64(std::ostream &out, const void *bytes, size_t nbytes);
void WriteBase64(std::ostream &out, const void *bytes, size_t length);
/// Decode @a len base-64 encoded characters in the buffer @a src, and store the
/// resulting decoded data in @a buf. @a buf will be resized as needed.
void DecodeBase64(const char *src, size_t len, std::vector<char> &buf);
/// Return the number of characters needed to encode @a nbytes in base-64. This
/// is equal to 4*nbytes/3, rounded up to the nearest multiple of 4.
size_t NumBase64Chars(size_t nbytes);
void DecodeBase64(const char *src, size_t len, std::vector<unsigned char> &buf);
} // namespace mfem::bin_io
-2
View File
@@ -17,7 +17,6 @@ list(APPEND SRCS
complex_operator.cpp
constraints.cpp
densemat.cpp
fdsolver.cpp
symmat.cpp
handle.cpp
matrix.cpp
@@ -40,7 +39,6 @@ list(APPEND HDRS
dinvariants.hpp
symmat.hpp
dtensor.hpp
fdsolver.hpp
handle.hpp
invariants.hpp
kernels.hpp
+8 -10
View File
@@ -208,14 +208,11 @@ void AmgXSolver::DefaultParameters(const AMGX_MODE amgxMode_,
" \"config_version\": 2, \n"
" \"solver\": { \n"
" \"solver\": \"AMG\", \n"
" \"scope\": \"main\", \n"
" \"smoother\": \"JACOBI_L1\", \n"
" \"presweeps\": 1, \n"
" \"interpolator\": \"D2\", \n"
" \"max_row_sum\" : 0.9, \n"
" \"strength_threshold\" : 0.25, \n"
" \"postsweeps\": 1, \n"
" \"max_iters\": 1, \n"
" \"interpolator\": \"D2\", \n"
" \"max_iters\": 2, \n"
" \"convergence\": \"ABSOLUTE\", \n"
" \"cycle\": \"V\"";
if (verbose)
{
@@ -242,21 +239,22 @@ void AmgXSolver::DefaultParameters(const AMGX_MODE amgxMode_,
" \"solver\": \"AMG\", \n"
" \"smoother\": { \n"
" \"scope\": \"jacobi\", \n"
" \"solver\": \"JACOBI_L1\" \n"
" \"solver\": \"BLOCK_JACOBI\", \n"
" \"relaxation_factor\": 0.7 \n"
" }, \n"
" \"presweeps\": 1, \n"
" \"interpolator\": \"D2\", \n"
" \"max_row_sum\" : 0.9, \n"
" \"strength_threshold\" : 0.25, \n"
" \"max_iters\": 1, \n"
" \"max_iters\": 2, \n"
" \"scope\": \"amg\", \n"
" \"max_levels\": 100, \n"
" \"cycle\": \"V\", \n"
" \"postsweeps\": 1 \n"
" }, \n"
" \"solver\": \"PCG\", \n"
" \"max_iters\": 150, \n"
" \"convergence\": \"RELATIVE_INI_CORE\", \n"
" \"max_iters\": 100, \n"
" \"convergence\": \"RELATIVE_MAX\", \n"
" \"scope\": \"main\", \n"
" \"tolerance\": 1e-12, \n"
" \"monitor_residual\": 1, \n"
+17 -276
View File
@@ -50,10 +50,6 @@ dsyevr_(char *JOBZ, char *RANGE, char *UPLO, int *N, double *A, int *LDA,
double *W, double *Z, int *LDZ, int *ISUPPZ, double *WORK, int *LWORK,
int *IWORK, int *LIWORK, int *INFO);
extern "C" void
dgeev_(const char * jobvl, const char * jobvr, int *n, double * A, int * lda,
double * wr, double * wl, double * vl, int * ldvl, double * vr, int * ldvr,
double * work, int * lwork, int * info);
extern "C" void
dsyev_(char *JOBZ, char *UPLO, int *N, double *A, int *LDA, double *W,
double *WORK, int *LWORK, int *INFO);
extern "C" void
@@ -178,14 +174,7 @@ const double &DenseMatrix::Elem(int i, int j) const
void DenseMatrix::Mult(const double *x, double *y) const
{
const double *data = Read();
const int h = height;
const int w = width;
MFEM_FORALL(i, 1,
{
kernels::Mult(h, w, data, x, y);
});
kernels::Mult(height, width, Data(), x, y);
}
void DenseMatrix::Mult(const Vector &x, Vector &y) const
@@ -193,9 +182,7 @@ void DenseMatrix::Mult(const Vector &x, Vector &y) const
MFEM_ASSERT(height == y.Size() && width == x.Size(),
"incompatible dimensions");
const double *dx = x.Read();
double *dy = y.ReadWrite();
Mult(dx, dy);
Mult((const double *)x, (double *)y);
}
double DenseMatrix::operator *(const DenseMatrix &m) const
@@ -2016,7 +2003,7 @@ void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
MFEM_ASSERT(a.Height() == b.Height() && a.Width() == c.Width() &&
b.Width() == c.Height(), "incompatible dimensions");
#if defined(MFEM_USE_LAPACK) && !defined(MFEM_USE_CUDA) && !defined(MFEM_USE_HIP)
#ifdef MFEM_USE_LAPACK
static char transa = 'N', transb = 'N';
static double alpha = 1.0, beta = 0.0;
int m = b.Height(), n = c.Width(), k = b.Width();
@@ -2027,13 +2014,10 @@ void Mult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
const int ah = a.Height();
const int aw = a.Width();
const int bw = b.Width();
double *ad = a.ReadWrite();
const double *bd = b.Read();
const double *cd = c.Read();
MFEM_FORALL(i, 1,
{
kernels::Mult(ah, aw, bw, bd, cd, ad);
});
double *ad = a.Data();
const double *bd = b.Data();
const double *cd = c.Data();
kernels::Mult(ah,aw,bw,bd,cd,ad);
#endif
}
@@ -2875,155 +2859,6 @@ void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
}
}
void KronProd(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C)
{
const int ah = A.Height();
const int aw = A.Width();
const int bh = B.Height();
const int bw = B.Width();
C.SetSize(ah*bh,aw*bw);
const double * ad = A.Read();
const double * bd = B.Read();
double * cd = C.ReadWrite();
MFEM_FORALL(i, 1,
{
for (int ja = 0; ja<aw; ++ja)
for (int jb = 0; jb<bw; ++jb)
for (int ia = 0; ia<ah; ++ia)
for (int ib = 0; ib<bh; ++ib)
cd[bh*ia + ib + ah*bh*(bw*ja + jb)]
= ad[ia + ja * ah] * bd[ib + jb*bh];
});
}
#if 0 // this is a finer level parallel KronProd
void KronProd2(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C)
{
const int ah = A.Height();
const int aw = A.Width();
const int bh = B.Height();
const int bw = B.Width();
const int ch = ah * bh;
const int cw = aw * bw;
C.SetSize(ch, cw);
const double *ad = A.Read();
const double *bd = B.Read();
double *cd = C.ReadWrite();
MFEM_FORALL(i, ch * cw,
{
const int jc = i / ch;
const int ic = i - jc * ch;
const int ja = jc / bw;
const int jb = jc - ja * bw;
const int ia = ic / bh;
const int ib = ic - ia * bh;
cd[jc * ch + ic] = ad[ja * ah + ia] * bd[jb * bh + ib];
});
}
#endif
void KronMult(const DenseMatrix &A, const DenseMatrix &B, const Vector &r,
Vector &z)
{
const int nA = A.Height();
const int mA = A.Width();
const int nB = B.Height();
const int mB = B.Width();
const int nr = r.Size();
MFEM_VERIFY(nr == mA*mB, "Wrong size of Vector r");
z.SetSize(nA*nB);
#if !defined(MFEM_USE_CUDA)
DenseMatrix R(r.GetData(),mB,mA);
DenseMatrix X(nB,mA);
DenseMatrix Y(z.GetData(),nB,nA);
Mult(B,R,X);
MultABt(X,A,Y);
#else
const double *ad = A.Read();
const double *bd = B.Read();
const double *rd = r.Read();
double *zd = z.Write();
MFEM_FORALL(i, 1,
{
kernels::KronMult(nA, mA, ad, nB, mB, bd, rd, zd);
});
#endif
}
void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &R,
DenseMatrix & Z)
{
const int nA = A.Height();
const int nB = B.Height();
const int nR = R.Height();
const int mR = R.Width();
Z.SetSize(nA*nB,mR);
Vector r,z;
double * dataR = R.Data();
for (int i = 0; i<mR; i++)
{
r.SetDataAndSize(&dataR[i*nR],nR);
KronMult(A,B,r,z);
Z.SetCol(i,z);
}
}
void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &C,
const Vector &r, Vector &z)
{
const int nA = A.Height();
const int mA = A.Width();
const int nB = B.Height();
const int mB = B.Width();
const int nC = C.Height();
const int mC = C.Width();
const int nr = r.Size();
MFEM_VERIFY(nr == mA*mB*mC, "Wrong size of Vector r");
z.SetSize(nA*nB*nC);
#if !defined(MFEM_USE_CUDA)
double * dataR = r.GetData();
DenseMatrix R(dataR,mC,mA*mB);
DenseMatrix X(nC,mA*mB);
Mult(C,R,X);
X.Transpose();
DenseMatrix Z(z.GetData(),mA*mB,nC);
KronMult(A,B,X,Z);
Z.Transpose();
#else
const double *ad = A.Read();
const double *bd = B.Read();
const double *cd = C.Read();
const double *rd = r.Read();
double *zd = z.Write();
MFEM_FORALL(i, 1,
{
kernels::KronMult(nA, mA, ad, nB, mB, bd, nC, mC, cd, rd, zd);
});
#endif
}
void KronMult(const Array<DenseMatrix *> & A, const Vector & r, Vector & z)
{
int dim = A.Size();
if (dim == 2)
{
KronMult(*A[0],*A[1],r,z);
}
else if (dim == 3)
{
KronMult(*A[0],*A[1],*A[2], r,z);
}
else
{
MFEM_ABORT("KronMult::Wrong dimension");
}
}
bool LUFactors::Factor(int m, double TOL)
{
@@ -3475,105 +3310,23 @@ DenseMatrixInverse::~DenseMatrixInverse()
delete [] lu.ipiv;
}
void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
const Vector &r, Vector & z)
{
// A and B are square matrices
z.SetSize(r.Size());
int nA = A.Height();
int nB = B.Height();
DenseMatrix R(r.GetData(),nB,nA);
DenseMatrix X(nB,nA);
B.Mult(R,X);
X.Transpose();
DenseMatrix Y(z.GetData(),nA,nB);
A.Mult(X,Y);
Y.Transpose();
}
void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
const DenseMatrix &R, DenseMatrix & Z)
{
// A and B are square matrices
int nR = R.Height();
int mR = R.Width();
Z.SetSize(nR,mR);
Vector r(nR);
Vector z(nR);
double * dataR = R.GetData();
double * dataZ = Z.GetData();
for (int i = 0; i<mR; i++)
{
r.SetData(&dataR[i*nR]);
z.SetData(&dataZ[i*nR]);
KronMult(A,B,r,z);
}
}
void KronMult(const DenseMatrixInverse &A, const DenseMatrixInverse &B,
const DenseMatrixInverse &C, const Vector &r, Vector & z)
{
// A, B and C are square matrices
int n = r.Size();
z.SetSize(n);
int nA = A.Height();
int nB = B.Height();
int nC = C.Height();
double * dataR = r.GetData();
DenseMatrix R(dataR,nC,nA*nB);
DenseMatrix X(nC,nA*nB);
C.Mult(R,X);
X.Transpose();
DenseMatrix Z(z.GetData(),nC,nA*nB);
KronMult(A,B,X,Z);
Z.Transpose();
}
void KronMult(const Array<DenseMatrixInverse *> & A, const Vector & r,
Vector & z)
{
int dim = A.Size();
if (dim == 2)
{
KronMult(*A[0],*A[1],r,z);
}
else if (dim == 3)
{
KronMult(*A[0],*A[1],*A[2], r,z);
}
else
{
MFEM_ABORT("KronMult::Wrong dimension");
}
}
DenseMatrixEigensystem::DenseMatrixEigensystem(DenseMatrix &m, bool sym_)
DenseMatrixEigensystem::DenseMatrixEigensystem(DenseMatrix &m)
: mat(m)
{
n = mat.Width();
EVal.SetSize(n);
EVali.SetSize(n);
EVect.SetSize(n);
ev.SetDataAndSize(NULL, n);
#ifdef MFEM_USE_LAPACK
sym = sym_;
jobz = 'V';
uplo = 'U';
lwork = -1;
double qwork;
if (sym)
{
uplo = 'U';
dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
&qwork, &lwork, &info);
}
else
{
char jobvl = 'N';
int ldvl = 1;
dgeev_(&jobvl,&jobz,&n, mat.GetData(), &n, EVal.GetData(), EVali.GetData(),
nullptr, &ldvl, EVect.GetData(), &n, &qwork, &lwork, &info);
}
dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
&qwork, &lwork, &info);
lwork = (int) qwork;
work = new double[lwork];
#endif
@@ -3585,7 +3338,6 @@ DenseMatrixEigensystem::DenseMatrixEigensystem(
n(other.n)
{
#ifdef MFEM_USE_LAPACK
sym = other.sym;
jobz = other.jobz;
uplo = other.uplo;
lwork = other.lwork;
@@ -3604,24 +3356,13 @@ void DenseMatrixEigensystem::Eval()
#endif
#ifdef MFEM_USE_LAPACK
if (sym)
{
EVect = mat;
dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
work, &lwork, &info);
}
else
{
char jobvl = 'N';
int ldvl = 1;
DenseMatrix T = mat; // mat is overwritten by dgeev
dgeev_(&jobvl,&jobz,&n, T.GetData(), &n, EVal.GetData(), EVali.GetData(),
nullptr, &ldvl, EVect.GetData(), &n, work, &lwork, &info);
}
EVect = mat;
dsyev_(&jobz, &uplo, &n, EVect.Data(), &n, EVal.GetData(),
work, &lwork, &info);
if (info != 0)
{
string lpck = (sym) ? "DSYEV" : "DGEEV";
mfem::err << "DenseMatrixEigensystem::Eval(): " << lpck << "error code: "
mfem::err << "DenseMatrixEigensystem::Eval(): DSYEV error code: "
<< info << endl;
mfem_error();
}
+2 -36
View File
@@ -523,23 +523,6 @@ void AddMult_a_VWt(const double a, const Vector &v, const Vector &w,
/// VVt += a * v v^t
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt);
/// C = A ⊗ B
void KronProd(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C);
void KronProd2(const DenseMatrix & A, const DenseMatrix & B, DenseMatrix & C);
/// z = (A ⊗ B) r = vec(B R A^T), where R := vec^-1 (r)
void KronMult(const DenseMatrix &A, const DenseMatrix &B, const Vector &r,
Vector & z);
/// z = (A ⊗ B) R
void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &R,
DenseMatrix & Z);
/// z = ( A ⊗ B ⊗ C ) r
void KronMult(const DenseMatrix &A, const DenseMatrix &B, const DenseMatrix &C,
const Vector &r, Vector & z);
void KronMult(const Array<DenseMatrix *> & A, const Vector & r, Vector & z);
/** Class that can compute LU factorization of external data and perform various
operations with the factored data. */
@@ -702,33 +685,16 @@ public:
virtual ~DenseMatrixInverse();
};
/// z = (A^-1 ⊗ B^-1) r = vec(B^-1 R A^-T), where R := vec^-1 (r)
void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
const Vector &r, Vector & z);
/// z = (A^-1 ⊗ B^-1) R
void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
const DenseMatrix &R, DenseMatrix & Z);
/// z = ( A^-1 ⊗ B^-1 ⊗ C^-1 ) r
void KronMult(const DenseMatrixInverse &Ainv, const DenseMatrixInverse &Binv,
const DenseMatrixInverse &Cinv, const Vector &r, Vector & z);
void KronMult(const Array<DenseMatrixInverse *> & A, const Vector & r,
Vector & z);
class DenseMatrixEigensystem
{
DenseMatrix &mat;
Vector EVal;
// Possible non zero imaginary part of Eigenvalues
Vector EVali;
DenseMatrix EVect;
Vector ev;
int n;
#ifdef MFEM_USE_LAPACK
bool sym;
double *work;
char jobz, uplo;
int lwork, info;
@@ -736,10 +702,10 @@ class DenseMatrixEigensystem
public:
DenseMatrixEigensystem(DenseMatrix &m, bool sym_ = false);
DenseMatrixEigensystem(DenseMatrix &m);
DenseMatrixEigensystem(const DenseMatrixEigensystem &other);
void Eval();
Vector &Eigenvalues(bool imag = false) { return imag ? EVali : EVal; }
Vector &Eigenvalues() { return EVal; }
DenseMatrix &Eigenvectors() { return EVect; }
double Eigenvalue(int i) { return EVal(i); }
const Vector &Eigenvector(int i)
-143
View File
@@ -1,143 +0,0 @@
// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "linalg.hpp"
namespace mfem
{
void KronProdInvDiag(const Vector & a, const Vector & b, Vector & dinv)
{
int n = a.Size(), m = b.Size();
dinv.SetSize(n*m);
for (int j = 0; j<m; j++)
for (int i = 0; i<n; i++)
{
dinv(i*m+j) = 1./(a(i) + b(j));
}
}
void KronProdInvDiag(const Vector & a, const Vector & b,
const Vector & c, Vector & dinv)
{
int n = a.Size(), m = b.Size(), l = c.Size();
dinv.SetSize(n*m*l);
for (int k = 0; k<l; k++)
for (int j = 0; j<m; j++)
for (int i = 0; i<n; i++)
{
dinv(i*m*l+j*l+k) = 1./(a(i) + b(j) + c(k));
}
}
void KronProdInvDiag(const Array<Vector *> & X, Vector & dinv)
{
int dim = X.Size();
if (dim == 1)
{
int n = X[0]->Size();
dinv.SetSize(n);
for (int i = 0; i<n; i++) { dinv(i) = 1./(*X[0])(i); }
}
else if (dim == 2)
{
KronProdInvDiag(*X[0], *X[1], dinv);
}
else if (dim == 3)
{
KronProdInvDiag(*X[0], *X[1], *X[2], dinv);
}
else
{
MFEM_ABORT("KronProdInvDiag::Wrong dimension");
}
}
#ifdef MFEM_USE_LAPACK
FDSolver::FDSolver(const Array<DenseMatrix *> & A,
const Array<DenseMatrix *> & B)
{
MFEM_ASSERT(A.Size() == B.Size(), "DenseFDSolver: Incompatible Dimensions");
dim = A.Size();
int solver_size = 1;
for (int i = 0; i<dim; i++)
{
MFEM_ASSERT(A[i]->Height() == A[i]->Width(),
"DenseFDSolver: Matrix is not square");
MFEM_ASSERT(B[i]->Height() == B[i]->Width(),
"DenseFDSolver: Matrix is not square");
MFEM_ASSERT(A[i]->Height() == B[i]->Height(),
"DenseFDSolver: Matrices A and B have incompatible size");
solver_size *= A[i]->Height();
}
this->height = solver_size;
this->width = solver_size;
if (solver_size) { Setup(A,B); }
}
void FDSolver::Setup(const Array<DenseMatrix *> & A,
const Array<DenseMatrix *> & B)
{
EigSystem.SetSize(dim);
eigv.SetSize(dim);
Array<Vector *> evalues(dim);
SQ.SetSize(dim);
DenseMatrix D;
for (int i = 0; i<dim; i++)
{
DenseMatrixInverse Minv(*B[i]);
Minv.Mult(*A[i],D);
EigSystem[i] = new DenseMatrixEigensystem(D);
EigSystem[i]->Eval();
evalues[i] = &EigSystem[i]->Eigenvalues();
eigv[i] = &EigSystem[i]->Eigenvectors();
DenseMatrixInverse Qinv(*eigv[i]);
DenseMatrix Sdinv;
Minv.GetInverseMatrix(Sdinv);
SQ[i] = new DenseMatrix;
Qinv.Mult(Sdinv,*SQ[i]);
}
KronProdInvDiag(evalues,dinv);
}
void FDSolver::Mult(const Vector & r,Vector & z) const
{
MFEM_ASSERT(height == r.Size(),
"DenseFDSolver::Mult: Inconsistent vector size");
if (r.Size() == 0) { return; }
Vector rtemp;
KronMult(SQ,r,rtemp);
// 2. Diagonal solve;
rtemp *= dinv;
// 3. Modify RHS; z <-- (Q1 x Q2) rtemp
KronMult(eigv,rtemp,z);
}
FDSolver::~FDSolver()
{
if (height)
{
for (int i=0; i<dim; i++)
{
delete SQ[i];
delete EigSystem[i];
}
}
}
#endif // MFEM_USE_LAPACK
} // namespace mfem
-60
View File
@@ -1,60 +0,0 @@
// Copyright (c) 2010-2021, 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_FDSOLVER
#define MFEM_FDSOLVER
#include "../config/config.hpp"
#include "densemat.hpp"
namespace mfem
{
/// Computes the inverse diagonal dinv = (a⊗I + I⊗b)^-1
/// where a, b are diagonal matrices and I is the identity of the
/// appropriate size
void KronProdInvDiag(const Vector & a, const Vector & b, Vector & dinv);
/// Computes the inverse diagonal dinv = (a⊗I⊗I + I⊗b⊗I + I⊗I⊗c)^-1
/// where a, b, c are diagonal matrices and I is the identity of the
/// appropriate size
void KronProdInvDiag(const Vector & a, const Vector & b,
const Vector & c, Vector & dinv);
void KronProdInvDiag(const Array<Vector *> & X, Vector & dinv);
#ifdef MFEM_USE_LAPACK
/// In 2D it solves the system (A_0 ⊗ B_1 + B_0 ⊗ A_1) z = r
/// In 3D it solves the system
/// (A_0 ⊗ B_1 ⊗ B_2 + B_0 ⊗ A_1 ⊗ B_2 + B_0 ⊗ B_1 ⊗ A_2) z = r
class FDSolver: public Solver
{
private:
int dim = 2;
Array<DenseMatrixEigensystem *> EigSystem;
Array<DenseMatrix *> eigv; // eigenvectors
Array<DenseMatrix *> SQ;
mutable Vector dinv;
void Setup(const Array<DenseMatrix *> & A, const Array<DenseMatrix *> & B);
public:
FDSolver(const Array<DenseMatrix *> & A, const Array<DenseMatrix *> & B);
virtual void SetOperator(const Operator &op) {}
virtual void Mult(const Vector &r, Vector &z) const;
virtual ~FDSolver();
};
#endif // MFEM_USE_LAPACK
} // mfem name space
#endif // MFEM_FDSOLVER
+1 -8
View File
@@ -516,14 +516,7 @@ public:
/// Initialize all entries with value.
HypreParMatrix &operator=(double value)
{
#if MFEM_HYPRE_VERSION < 22200
internal::hypre_ParCSRMatrixSetConstantValues(A, value);
#else
hypre_ParCSRMatrixSetConstantValues(A, value);
#endif
return *this;
}
{ internal::hypre_ParCSRMatrixSetConstantValues(A, value); return *this; }
/** Perform the operation `*this += B`, assuming that both matrices use the
same row and column partitions and the same col_map_offd arrays, or B has
+2 -2
View File
@@ -1942,8 +1942,8 @@ HYPRE_Int
hypre_ParCSRMatrixSetConstantValues(hypre_ParCSRMatrix *A,
HYPRE_Complex value)
{
internal::hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixDiag(A), value);
internal::hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixOffd(A), value);
hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixDiag(A), value);
hypre_CSRMatrixSetConstantValues(hypre_ParCSRMatrixOffd(A), value);
return 0;
}
-10
View File
@@ -198,16 +198,6 @@ hypre_CSRMatrixSum(hypre_CSRMatrix *A,
HYPRE_Complex beta,
hypre_CSRMatrix *B);
#if MFEM_HYPRE_VERSION >= 22200
/** Provide an overloaded function for code consistency between HYPRE API
versions. */
inline hypre_CSRMatrix *hypre_CSRMatrixAdd(hypre_CSRMatrix *A,
hypre_CSRMatrix *B)
{
return ::hypre_CSRMatrixAdd(1.0, A, 1.0, B);
}
#endif
/** Return a new matrix containing the sum of A and B, assuming that both
matrices use the same row and column partitions. The col_map_offd do not
need to be the same, but a more efficient algorithm is used if that's the
+3 -63
View File
@@ -160,7 +160,7 @@ double Norml2(const int size, const T *data)
data of the input and output vectors. */
template<typename TA, typename TX, typename TY>
MFEM_HOST_DEVICE inline
void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
void Mult(const int height, const int width, TA *data, const TX *x, TY *y)
{
if (width == 0)
{
@@ -170,8 +170,7 @@ void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
}
return;
}
TA *d_col = (TA *) data;
TA *d_col = data;
TX x_col = x[0];
for (int row = 0; row < height; row++)
{
@@ -189,52 +188,6 @@ void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
}
}
template<typename TA, typename TB, typename TR, typename TZ>
MFEM_HOST_DEVICE inline
void KronMult(const int ah, const int aw, const TA *ad, const int bh, const int bw, const TB *bd, TR *r, TZ *z)
{
for (int i = 0; i < bh; i++)
{
for (int l = 0; l < ah; l++)
{
TZ t1 = 0.0;
for (int j = 0; j < bw; j++)
{
const TB t2 = bd[i + j * bh];
for (int k = 0; k < aw; k++)
{
t1 += t2 * r[j + k * bw] * ad[l + k * ah];
}
}
z[i + l * bh] = t1;
}
}
}
template<typename TA, typename TB, typename TC, typename TR, typename TZ>
MFEM_HOST_DEVICE inline
void KronMult(const int ah, const int aw, const TA *ad, const int bh, const int bw, const TB *bd, const int ch, const int cw, const TC *cd, TR *r, TZ *z)
{
for (int i = 0; i < ch; i++)
{
for (int l = 0; l < ah * bh; l++)
{
TZ t1 = 0.0;
for (int j = 0; j < cw; j++)
{
const TB t2 = cd[i + j * ch];
for (int k = 0; k < aw * bw; k++)
{
const TA ta = ad[(l / bh) + (k / bw) * ah];
const TB tb = bd[(l % bh) + (k % bw) * bh];
t1 += t2 * r[j + k * cw] * ta * tb;
}
}
z[i + l * ch] = t1;
}
}
}
/// Symmetrize a square matrix with given @a size and @a data: A -> (A+A^T)/2.
template<typename T>
MFEM_HOST_DEVICE inline
@@ -323,20 +276,6 @@ void Add(const int height, const int width, const TA *Adata, TB *Bdata)
}
}
/** @brief Compute B +=alpha*A, where the matrices A and B are of size
@a height x @a width with data @a Adata and @a Bdata. */
template<typename TA, typename TB>
MFEM_HOST_DEVICE inline
void Add(const int height, const int width,
const double alpha, const TA *Adata, TB *Bdata)
{
const int m = height * width;
for (int i = 0; i < m; i++)
{
Bdata[i] += alpha * Adata[i];
}
}
/** @brief Compute B = alpha*A, where the matrices A and B are of size
@a height x @a width with data @a Adata and @a Bdata. */
template<typename TA, typename TB>
@@ -351,6 +290,7 @@ void Set(const int height, const int width,
}
}
/** @brief Matrix-matrix multiplication: A = B * C, where the matrices A, B and
C are of sizes @a Aheight x @a Awidth, @a Aheight x @a Bwidth and @a Bwidth
x @a Awidth, respectively. */
-1
View File
@@ -31,7 +31,6 @@
#include "invariants.hpp"
#include "constraints.hpp"
#include "auxiliary.hpp"
#include "fdsolver.hpp"
#ifdef MFEM_USE_AMGX
#include "amgxsolver.hpp"
+6 -28
View File
@@ -232,7 +232,7 @@ void OperatorJacobiSmoother::Mult(const Vector &x, Vector &y) const
MFEM_FORALL(i, height, Y[i] += DI[i] * R[i]; );
}
OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
OperatorChebyshevSmoother::OperatorChebyshevSmoother(Operator* oper_,
const Vector &d,
const Array<int>& ess_tdofs,
int order_, double max_eig_estimate_)
@@ -246,15 +246,15 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
coeffs(order),
ess_tdof_list(ess_tdofs),
residual(N),
oper(&oper_) { Setup(); }
oper(oper_) { Setup(); }
#ifdef MFEM_USE_MPI
OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
OperatorChebyshevSmoother::OperatorChebyshevSmoother(Operator* oper_,
const Vector &d,
const Array<int>& ess_tdofs,
int order_, MPI_Comm comm, int power_iterations, double power_tolerance)
#else
OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
OperatorChebyshevSmoother::OperatorChebyshevSmoother(Operator* oper_,
const Vector &d,
const Array<int>& ess_tdofs,
int order_, int power_iterations, double power_tolerance)
@@ -267,7 +267,7 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
coeffs(order),
ess_tdof_list(ess_tdofs),
residual(N),
oper(&oper_)
oper(oper_)
{
OperatorJacobiSmoother invDiagOperator(diag, ess_tdofs, 1.0);
ProductOperator diagPrecond(&invDiagOperator, oper, false, false);
@@ -284,28 +284,6 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
Setup();
}
OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator* oper_,
const Vector &d,
const Array<int>& ess_tdofs,
int order_, double max_eig_estimate_)
: OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, max_eig_estimate_) { }
#ifdef MFEM_USE_MPI
OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator* oper_,
const Vector &d,
const Array<int>& ess_tdofs,
int order_, MPI_Comm comm, int power_iterations, double power_tolerance)
: OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, comm,
power_iterations, power_tolerance) { }
#else
OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator* oper_,
const Vector &d,
const Array<int>& ess_tdofs,
int order_, int power_iterations, double power_tolerance)
: OperatorChebyshevSmoother(*oper_, d, ess_tdofs, order_, power_iterations,
power_tolerance) { }
#endif
void OperatorChebyshevSmoother::Setup()
{
// Invert diagonal
@@ -2555,7 +2533,7 @@ BlockILU::BlockILU(int block_size_,
reordering(reordering_)
{ }
BlockILU::BlockILU(const Operator &op,
BlockILU::BlockILU(Operator &op,
int block_size_,
Reordering reordering_,
int k_fill_)
+4 -25
View File
@@ -208,13 +208,7 @@ public:
the matrix-free setting. The estimated largest eigenvalue of the
diagonally preconditoned operator must be provided via
max_eig_estimate. */
OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
const Array<int>& ess_tdof_list,
int order, double max_eig_estimate);
/// Deprecated: see pass-by-reference version above
MFEM_DEPRECATED
OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
OperatorChebyshevSmoother(Operator* oper_, const Vector &d,
const Array<int>& ess_tdof_list,
int order, double max_eig_estimate);
@@ -226,28 +220,13 @@ public:
accuracy of the estimated eigenvalue may be controlled via
power_iterations and power_tolerance. */
#ifdef MFEM_USE_MPI
OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
const Array<int>& ess_tdof_list,
int order, MPI_Comm comm = MPI_COMM_NULL,
int power_iterations = 10,
double power_tolerance = 1e-8);
/// Deprecated: see pass-by-reference version above
MFEM_DEPRECATED
OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
OperatorChebyshevSmoother(Operator* oper_, const Vector &d,
const Array<int>& ess_tdof_list,
int order, MPI_Comm comm = MPI_COMM_NULL,
int power_iterations = 10,
double power_tolerance = 1e-8);
#else
OperatorChebyshevSmoother(const Operator &oper_, const Vector &d,
const Array<int>& ess_tdof_list,
int order, int power_iterations = 10,
double power_tolerance = 1e-8);
/// Deprecated: see pass-by-reference version above
MFEM_DEPRECATED
OperatorChebyshevSmoother(const Operator* oper_, const Vector &d,
OperatorChebyshevSmoother(Operator* oper_, const Vector &d,
const Array<int>& ess_tdof_list,
int order, int power_iterations = 10,
double power_tolerance = 1e-8);
@@ -752,7 +731,7 @@ public:
* case that @a op is a HypreParMatrix, the ILU factorization is performed
* on the diagonal blocks of the parallel decomposition.
*/
BlockILU(const Operator &op, int block_size_ = 1,
BlockILU(Operator &op, int block_size_ = 1,
Reordering reordering_ = Reordering::MINIMUM_DISCARDED_FILL,
int k_fill_ = 0);
-12
View File
@@ -37,7 +37,6 @@ MFEM makefile targets:
make distclean
make style
make tags
make hooks
Examples:
@@ -97,8 +96,6 @@ make style
make tags
Generate a vi or Emacs compatible TAGS file in ${MFEM_DIR}/TAGS. Requires
functional "etags" and "egrep" in the user ${PATH}.
make hooks
Creates symlinks to the hooks in the `.git/hooks` directory.
endef
# Save the MAKEOVERRIDES for cases where we explicitly want to pass the command
@@ -762,15 +759,6 @@ endif
@cd $(MFEM_REAL_DIR) && $(ETAGS_BIN) --class-qualify \
--declarations -o $(MFEM_REAL_DIR)/TAGS $(MFEM_TRACKED_SOURCE)
# Creates symlinks to the hooks in the `.git/hooks` directory. Individual
# hooks can be enabled by manually creating symlinks. Hooks can be customized
# using hard copies (trading off with automated updates).
.PHONY: hooks
hooks:
@cd $(MFEM_DIR)/.git/hooks && \
ln -s ../../config/githooks/pre-commit pre-commit; \
ln -s ../../config/githooks/pre-push pre-push;
# Print the contents of a makefile variable, e.g.: 'make print-MFEM_LIBS'.
print-%:
$(info [ variable name]: $*)
+309 -4
View File
@@ -4440,6 +4440,311 @@ void Mesh::MakeSimplicial_(const Mesh &orig_mesh, int *vglobal)
MFEM_ASSERT(CheckBdrElementOrientation(false) == 0, "");
}
Mesh Mesh::ExtractMesh(const Mesh &orig_mesh, const Array<int> & elems)
{
Mesh mesh;
mesh.ExtractMesh_(orig_mesh, elems);
return mesh;
}
void Mesh::ExtractMesh_(const Mesh &orig_mesh, const Array<int> & elems)
{
int dim = orig_mesh.Dimension();
int sdim = orig_mesh.SpaceDimension();
int nv = orig_mesh.GetNV();
int nf = orig_mesh.GetNumFaces();
// vertex marker
Array<int> vmarker(nv); vmarker = 0;
int new_nv = 0;
int new_ne = elems.Size();
// Count and mark the vertices to be added to the new mesh
Array<int> vertices;
for (int iel=0; iel<new_ne; ++iel)
{
int el = elems[iel];
orig_mesh.GetElementVertices(el,vertices);
for (int iv=0; iv<vertices.Size(); ++iv)
{
int v = vertices[iv];
if (vmarker[v]) { continue; }
vmarker[v] = 1;
new_nv++;
}
}
// Count the bdry elements to be added to the new mesh
Array<int> BoundaryMarker(nf); BoundaryMarker = 0;
int new_nbe = 0;
for (int i = 0; i<new_ne; i++)
{
int el = elems[i];
Array<int> faces, ori;
switch (dim)
{
case 1: orig_mesh.GetElementVertices(el,faces); break;
case 2: orig_mesh.GetElementEdges(el,faces,ori); break;
default: orig_mesh.GetElementFaces(el,faces,ori); break;
}
for (int f=0; f<faces.Size(); ++f) { BoundaryMarker[faces[f]]++; }
}
for (int f=0; f<nf; ++f) { if (BoundaryMarker[f] == 1) { new_nbe++; } }
InitMesh(dim,sdim,new_nv,new_ne,new_nbe);
// Add vertices
int vk = 0;
for (int iv = 0; iv<nv; ++iv)
{
if (!vmarker[iv]) { continue; }
AddVertex(orig_mesh.GetVertex(iv));
// save the vertex unique id to the marker
vmarker[iv] = ++vk;
}
// add BdrElements
int attr = 0;
for (int f=0; f<nf; ++f)
{
if (BoundaryMarker[f] != 1) { continue; }
const Element * el = orig_mesh.GetFace(f);
Element * new_bel = nullptr;
if (dim == 1)
{
int pt = vmarker[f]-1;
new_bel = (Element*)new Point(&pt);
new_bel->SetAttribute(++attr);
}
else
{
new_bel = NewElement(el->GetGeometryType());
int nv0 = el->GetNVertices();
const int * v0 = el->GetVertices();
Array<int> v1(nv0);
for (int i=0; i<nv0; i++)
{
v1[i] = vmarker[v0[i]]-1;
}
new_bel->SetVertices(v1.GetData());
new_bel->SetAttribute(++attr);
}
AddBdrElement(new_bel);
}
// add elements
for (int e=0; e<new_ne; ++e)
{
const Element * el = orig_mesh.GetElement(elems[e]);
Element * new_el = NewElement(el->GetGeometryType());
int nv0 = el->GetNVertices();
const int * v0 = el->GetVertices();
Array<int> v1(nv0);
for (int i=0; i<nv0; i++) { v1[i] = vmarker[v0[i]]-1; }
new_el->SetVertices(v1.GetData());
new_el->SetAttribute(el->GetAttribute());
AddElement(new_el);
}
FinalizeTopology();
const GridFunction * orig_nodes = orig_mesh.GetNodes();
if (orig_nodes)
{
const FiniteElementSpace * orig_fes = orig_nodes->FESpace();
Ordering::Type ordering = orig_fes->GetOrdering();
int order = orig_fes->FEColl()->GetOrder();
bool discont = orig_fes->IsDGSpace();
SetCurvature(order, discont, sdim, ordering);
const FiniteElementSpace * new_fes = GetNodalFESpace();
GridFunction * new_nodes = GetNodes();
Array<int> orig_vdofs;
Array<int> new_vdofs;
Vector vec;
// Copy nodes to submesh
for (int e = 0; e < new_ne; e++)
{
new_fes->GetElementVDofs(e, new_vdofs);
orig_fes->GetElementVDofs(elems[e], orig_vdofs);
orig_nodes->GetSubVector(orig_vdofs, vec);
new_nodes->SetSubVector(new_vdofs, vec);
}
}
Finalize();
}
Mesh Mesh::ExtractSurfaceMesh(const Mesh &orig_mesh, const Array<int> & faces)
{
Mesh mesh;
mesh.ExtractSurfaceMesh_(orig_mesh, faces);
return mesh;
}
void Mesh::ExtractSurfaceMesh_(const Mesh &orig_mesh, const Array<int> & faces)
{
int dim = orig_mesh.Dimension()-1;
MFEM_VERIFY(dim > 0, "Only dim > 1 is supported");
int sdim = orig_mesh.SpaceDimension();
int nv = orig_mesh.GetNV();
int nf = (dim == 2) ? orig_mesh.GetNEdges() : nv;
// vertex marker
Array<int> vmarker(nv); vmarker = 0;
int new_nv = 0;
int new_ne = faces.Size();
// Count and mark the vertices to be added to the new mesh
Array<int> vertices;
for (int f=0; f<new_ne; ++f)
{
int el = faces[f];
orig_mesh.GetFaceVertices(el,vertices);
for (int iv=0; iv<vertices.Size(); ++iv)
{
int v = vertices[iv];
if (vmarker[v]) { continue; }
vmarker[v] = 1;
new_nv++;
}
}
// Count the bdry elements to be added to the new mesh
Array<int> BoundaryMarker(nf); BoundaryMarker = 0;
int new_nbe = 0;
for (int i = 0; i<new_ne; i++)
{
int el = faces[i];
Array<int> edges, ori;
switch (dim)
{
case 1: orig_mesh.GetEdgeVertices(el,edges); break; // these are vertices
case 2: orig_mesh.GetFaceEdges(el,edges,ori); break;
default: MFEM_ABORT("Unreachable"); break;
}
for (int f=0; f<edges.Size(); ++f) { BoundaryMarker[edges[f]]++; }
}
for (int f=0; f<nf; ++f) { if (BoundaryMarker[f] == 1) { new_nbe++; } }
InitMesh(dim,sdim,new_nv,new_ne,new_nbe);
// Add vertices
int vk = 0;
for (int iv = 0; iv<nv; ++iv)
{
if (!vmarker[iv]) { continue; }
AddVertex(orig_mesh.GetVertex(iv));
// save the vertex unique id to the marker
vmarker[iv] = ++vk;
}
// add BdrElements
int attr = 0;
for (int f=0; f<nf; ++f)
{
if (BoundaryMarker[f] != 1) { continue; }
Element * new_bel = nullptr;
if (dim == 1)
{
int pt = vmarker[f]-1;
new_bel = (Element*)new Point(&pt);
new_bel->SetAttribute(++attr);
}
else
{
new_bel = NewElement(mfem::Geometry::SEGMENT);
int nv0 = 2;
Array<int> vert;
orig_mesh.GetEdgeVertices(f,vert);
Array<int> v1(nv0);
for (int i=0; i<nv0; i++)
{
v1[i] = vmarker[vert[i]]-1;
}
new_bel->SetVertices(v1.GetData());
new_bel->SetAttribute(++attr);
}
AddBdrElement(new_bel);
}
// add elements
for (int e=0; e<new_ne; ++e)
{
const Element * el = orig_mesh.GetFace(faces[e]);
Element * new_el = NewElement(el->GetGeometryType());
int nv0 = el->GetNVertices();
const int * v0 = el->GetVertices();
Array<int> v1(nv0);
for (int i=0; i<nv0; i++) { v1[i] = vmarker[v0[i]]-1; }
new_el->SetVertices(v1.GetData());
new_el->SetAttribute(el->GetAttribute());
AddElement(new_el);
}
FinalizeTopology();
const GridFunction * orig_nodes = orig_mesh.GetNodes();
if (orig_nodes)
{
const FiniteElementSpace * orig_fes = orig_nodes->FESpace();
Ordering::Type ordering = orig_fes->GetOrdering();
int order = orig_fes->FEColl()->GetOrder();
bool discont = orig_fes->IsDGSpace();
SetCurvature(order, discont, sdim, ordering);
const FiniteElementSpace * new_fes = GetNodalFESpace();
GridFunction * new_nodes = GetNodes();
Array<int> orig_vdofs;
Array<int> new_vdofs;
Vector vec;
// Copy nodes to submesh
for (int e = 0; e < new_ne; e++)
{
new_fes->GetElementVDofs(e, new_vdofs);
if (!discont)
{
orig_fes->GetFaceVDofs(faces[e], orig_vdofs);
orig_nodes->GetSubVector(orig_vdofs, vec);
}
else
{
const FiniteElement * el = new_fes->GetFE(e);
const IntegrationRule & ir = el->GetNodes();
int np = ir.GetNPoints();
FaceElementTransformations * Tr =
const_cast<Mesh *>(&orig_mesh)->GetFaceElementTransformations(faces[e]);
int el1 = Tr->Elem1No;
vec.SetSize(new_vdofs.Size());
for (int i = 0; i<np; i++)
{
Tr->SetAllIntPoints(&ir[i]);
const IntegrationPoint & ip = Tr->GetElement1IntPoint();
Vector val;
orig_nodes->GetVectorValue(el1,ip,val);
for (int j = 0; j<val.Size(); j++)
{
vec[i+j*np] = val[j];
}
}
}
new_nodes->SetSubVector(new_vdofs, vec);
}
}
Finalize();
}
Mesh Mesh::MakePeriodic(const Mesh &orig_mesh, const std::vector<int> &v2v)
{
Mesh periodic_mesh(orig_mesh, true); // Make a copy of the original mesh
@@ -11478,10 +11783,10 @@ FaceGeometricFactors::FaceGeometricFactors(const Mesh *mesh,
const int NF = fespace->GetNFbyType(type);
const int NQ = ir.GetNPoints();
const FaceRestriction *face_restr = fespace->GetFaceRestriction(
ElementDofOrdering::LEXICOGRAPHIC,
type,
L2FaceValues::SingleValued );
const Operator *face_restr = fespace->GetFaceRestriction(
ElementDofOrdering::LEXICOGRAPHIC,
type,
L2FaceValues::SingleValued );
Vector Fnodes(face_restr->Height());
face_restr->Mult(*nodes, Fnodes);
+16
View File
@@ -486,6 +486,12 @@ protected:
// Internal helper used in MakeSimplicial (and ParMesh::MakeSimplicial).
void MakeSimplicial_(const Mesh &orig_mesh, int *vglobal);
/// Internal helper used in ExtractMesh
void ExtractMesh_(const Mesh &orig_mesh, const Array<int> & elems);
/// Internal helper used in ExtractMesh
void ExtractSurfaceMesh_(const Mesh &orig_mesh, const Array<int> & faces);
public:
Mesh() { SetEmpty(); }
@@ -579,6 +585,16 @@ public:
new mesh. Periodic meshes are not supported. */
static Mesh MakeSimplicial(const Mesh &orig_mesh);
/** Create a mesh by extracting from @a orig_mesh the given list of
@a elems. Periodic and high-order meshes are also supported.
@warning NonConforming meshes are not supported. */
static Mesh ExtractMesh(const Mesh &orig_mesh, const Array<int> & elems);
/** Create a surface mesh by extracting from @a orig_mesh the given list of
@a faces. Periodic and high-order meshes are also supported.
@warning NonConforming meshes are not supported. */
static Mesh ExtractSurfaceMesh(const Mesh &orig_mesh, const Array<int> & faces);
/// Create a periodic mesh by identifying vertices of @a orig_mesh.
/** Each vertex @a i will be mapped to vertex @a v2v[i], such that all
vertices that are coincident under the periodic mapping get mapped to
+34 -63
View File
@@ -691,31 +691,16 @@ struct BufferReader : BufferReaderBase
BufferReader(bool compressed_, HeaderType header_type_)
: compressed(compressed_), header_type(header_type_) { }
/// Return the number of bytes of each header entry.
size_t HeaderEntrySize() const
{
return header_type == UINT64_HEADER ? sizeof(uint64_t) : sizeof(uint32_t);
}
/// Return the value of the header entry pointer to by @a header_buf. The
/// value is stored as either uint32_t or uint64_t, according to the @a
/// header_type, and is returned as uint64_t.
uint64_t ReadHeaderEntry(const char *header_buf) const
{
return (header_type == UINT64_HEADER) ? bin_io::read<uint64_t>(header_buf)
: bin_io::read<uint32_t>(header_buf);
}
/// Return the number of bytes in the header. The header consists of one
/// integer if the data is uncompressed, and @a N + 3 integers if the data is
/// compressed, where @a N is the number of blocks. The integers are either
/// 32 or 64 bytes depending on the value of @a header_type. The number of
/// blocks is determined by reading the first integer (of type @a
/// header_type) pointed to by @a header_buf.
int NumHeaderBytes(const char *header_buf) const
/// integer if the data is uncompressed, and four integers if the data is
/// compressed. The integers are either 32 or 64 bytes depending on the value
/// of @a header_type.
int NumHeaderBytes() const
{
if (!compressed) { return HeaderEntrySize(); }
return (3 + ReadHeaderEntry(header_buf))*HeaderEntrySize();
int num_entries = compressed ? 4 : 1;
int entry_size = (header_type == UINT64_HEADER)
? sizeof(uint64_t) : sizeof(uint32_t);
return num_entries*entry_size;
}
/// Read @a n elements of type @a F from the source buffer @a buf into the
@@ -726,42 +711,32 @@ struct BufferReader : BufferReaderBase
void ReadBinaryWithHeader(const char *header_buf, const char *buf,
void *dest_void, int n) const
{
std::vector<char> uncompressed_data;
std::vector<unsigned char> uncompressed_data;
T *dest = static_cast<T*>(dest_void);
if (compressed)
{
#ifdef MFEM_USE_ZLIB
// The header has format (where header_t is uint32_t or uint64_t):
// header_t number_of_blocks;
// header_t uncompressed_block_size;
// header_t uncompressed_last_block_size;
// header_t compressed_size[number_of_blocks];
int header_entry_size = HeaderEntrySize();
int nblocks = ReadHeaderEntry(header_buf);
header_buf += header_entry_size;
std::vector<int> header(nblocks + 2);
for (int i=0; i<nblocks+2; ++i)
uint64_t header[4];
if (header_type == UINT32_HEADER)
{
header[i] = ReadHeaderEntry(header_buf);
header_buf += header_entry_size;
uint32_t *header_32 = (uint32_t *)header_buf;
for (int i=0; i<4; ++i) { header[i] = header_32[i]; }
}
uncompressed_data.resize((nblocks-1)*header[0] + header[1]);
Bytef *dest_ptr = (Bytef *)uncompressed_data.data();
Bytef *dest_start = dest_ptr;
const Bytef *source_ptr = (const Bytef *)buf;
for (int i=0; i<nblocks; ++i)
else
{
uLongf source_len = header[i+2];
uLong dest_len = (i == nblocks-1) ? header[1] : header[0];
int res = uncompress(dest_ptr, &dest_len, source_ptr, source_len);
MFEM_VERIFY(res == Z_OK, "Error uncompressing");
dest_ptr += dest_len;
source_ptr += source_len;
uint64_t *header_64 = (uint64_t *)header_buf;
for (int i=0; i<4; ++i) { header[i] = header_64[i]; }
}
MFEM_VERIFY(int(sizeof(F)*n) == (dest_ptr - dest_start),
"AppendedData: wrong data size");
buf = uncompressed_data.data();
MFEM_VERIFY(header[0] == 1, "Multiple compressed blocks not supported");
uLongf dest_len = header[1];
uncompressed_data.resize(dest_len);
int res = uncompress(uncompressed_data.data(), &dest_len,
(const Bytef *)buf, header[3]);
MFEM_VERIFY(res == Z_OK, "Error uncompressing");
MFEM_VERIFY(sizeof(F)*n == dest_len, "AppendedData: wrong data size");
buf = (const char *)uncompressed_data.data();
#else
MFEM_ABORT("MFEM must be compiled with zlib enabled to uncompress.")
#endif
@@ -804,7 +779,7 @@ struct BufferReader : BufferReaderBase
/// buffer @a dest. The input buffer contains both the header and the data.
void ReadBinary(const char *buf, void *dest, int n) const override
{
ReadBinaryWithHeader(buf, buf + NumHeaderBytes(buf), dest, n);
ReadBinaryWithHeader(buf, buf + NumHeaderBytes(), dest, n);
}
/// Read @a n elements of type @a F from base-64 encoded source buffer into
@@ -821,25 +796,21 @@ struct BufferReader : BufferReaderBase
}
if (compressed)
{
// Decode the first entry of the header, which we need to determine
// how long the rest of the header is.
std::vector<char> nblocks_buf;
int nblocks_b64 = bin_io::NumBase64Chars(HeaderEntrySize());
bin_io::DecodeBase64(txt, nblocks_b64, nblocks_buf);
std::vector<char> data, header;
std::vector<unsigned char> data, header;
// Compute number of characters needed to encode header in base 64,
// then round to nearest multiple of 4 to take padding into account.
int header_b64 = bin_io::NumBase64Chars(NumHeaderBytes(nblocks_buf.data()));
int b64_header = ((4*NumHeaderBytes()/3) + 3) & ~3;
// If data is compressed, header is encoded separately
bin_io::DecodeBase64(txt, header_b64, header);
bin_io::DecodeBase64(txt + header_b64, strlen(txt)-header_b64, data);
ReadBinaryWithHeader(header.data(), data.data(), dest, n);
bin_io::DecodeBase64(txt, b64_header, header);
bin_io::DecodeBase64(txt + b64_header, strlen(txt)-b64_header, data);
ReadBinaryWithHeader((const char *)header.data(),
(const char *)data.data(), dest, n);
}
else
{
std::vector<char> data;
std::vector<unsigned char> data;
bin_io::DecodeBase64(txt, strlen(txt), data);
ReadBinary(data.data(), dest, n);
ReadBinary((const char *)data.data(), dest, n);
}
}
};
+2 -20
View File
@@ -4321,7 +4321,7 @@ const CoarseFineTransformations& NCMesh::GetRefinementTransforms()
if (!transforms.embeddings.Size())
{
transforms.Clear();
transforms.embeddings.SetSize(NElements);
transforms.embeddings.SetSize(leaf_elements.Size());
std::string ref_path;
ref_path.reserve(100);
@@ -4455,8 +4455,7 @@ struct RefType
void CoarseFineTransformations::GetCoarseToFineMap(
const mfem::Mesh &fine_mesh, Table &coarse_to_fine,
Array<int> &coarse_to_ref_type, Table &ref_type_to_matrix,
Array<mfem::Geometry::Type> &ref_type_to_geom,
bool get_coarse_to_fine_only) const
Array<mfem::Geometry::Type> &ref_type_to_geom) const
{
const int fine_ne = embeddings.Size();
int coarse_ne = -1;
@@ -4496,11 +4495,6 @@ void CoarseFineTransformations::GetCoarseToFineMap(
coarse_to_fine.GetJ()[i] = cf_j[i].two;
}
if (get_coarse_to_fine_only) { return; }
MFEM_VERIFY(fine_mesh.GetLastOperation() != Mesh::Operation::DEREFINE,
"GetCoarseToFineMap is not fully supported for derefined meshes."
" Set 'get_coarse_to_fine_only=true'.")
using internal::RefType;
using std::map;
using std::pair;
@@ -4542,18 +4536,6 @@ void CoarseFineTransformations::GetCoarseToFineMap(
ref_type_to_matrix.ShiftUpI();
}
void CoarseFineTransformations::GetCoarseToFineMap(const Mesh &fine_mesh,
Table &coarse_to_fine) const
{
Array<int> coarse_to_ref_type;
Table ref_type_to_matrix;
Array<mfem::Geometry::Type> ref_type_to_geom;
bool get_coarse_to_fine_only = true;
GetCoarseToFineMap(fine_mesh, coarse_to_fine, coarse_to_ref_type,
ref_type_to_matrix, ref_type_to_geom,
get_coarse_to_fine_only);
}
void NCMesh::ClearTransforms()
{
coarse_elements.DeleteAll();
+1 -5
View File
@@ -68,11 +68,7 @@ struct CoarseFineTransformations
Table &coarse_to_fine,
Array<int> &coarse_to_ref_type,
Table &ref_type_to_matrix,
Array<Geometry::Type> &ref_type_to_geom,
bool get_coarse_to_fine_only = false) const;
void GetCoarseToFineMap(const Mesh &fine_mesh,
Table &coarse_to_fine) const;
Array<Geometry::Type> &ref_type_to_geom) const;
void Clear();
bool IsInitialized() const;
+3
View File
@@ -297,6 +297,9 @@ protected: // implementation
virtual void Update();
virtual int GetNumGhostElements() const { return NGhostElements; }
virtual int GetNumGhostVertices() const { return NGhostVertices; }
/// Return the processor number for a global element number.
int Partition(long index, long total_elements) const
{ return index * NRanks / total_elements; }
+3 -7
View File
@@ -109,19 +109,15 @@ int main(int argc, char *argv[])
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
args.PrintUsage(cout);
return 1;
}
if (myid == 0) { args.PrintOptions(cout); }
args.PrintOptions(cout);
// Enable hardware devices such as GPUs, and programming models such as CUDA,
// OCCA, RAJA and OpenMP based on command line options.
Device device("cpu");
if (myid == 0) { device.Print(); }
device.Print();
// Refine the mesh.
Mesh mesh(mesh_file, 1, 1);
+2 -2
View File
@@ -96,7 +96,7 @@ SparseMatrix* AggToInteriorDof(const Array<int>& bdr_truedofs,
const SparseMatrix& agg_elem,
const SparseMatrix& elem_dof,
const HypreParMatrix& dof_truedof,
Array<HYPRE_Int>& agg_starts)
Array<int>& agg_starts)
{
OperatorPtr agg_dof(Mult(agg_elem, elem_dof));
SparseMatrix& agg_dof_ref = *agg_dof.As<SparseMatrix>();
@@ -131,7 +131,7 @@ SparseMatrix* AggToInteriorDof(const Array<int>& bdr_truedofs,
void DFSSpaces::MakeDofRelationTables(int level)
{
Array<HYPRE_Int> agg_starts(Array<HYPRE_Int>(l2_0_fes_->GetDofOffsets(), 2));
Array<int> agg_starts(Array<int>(l2_0_fes_->GetDofOffsets(), 2));
auto& elem_agg = (const SparseMatrix&)*l2_0_fes_->GetUpdateOperator();
OperatorPtr agg_elem(Transpose(elem_agg));
SparseMatrix& agg_el = *agg_elem.As<SparseMatrix>();
-1
View File
@@ -110,7 +110,6 @@ then
echo "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"
cp ${hostconfig_path} ${project_dir}/config/
ln -sf ${build_root}/data ${project_dir}/../
make all -j ${threads}
fi
-1
View File
@@ -28,7 +28,6 @@ set(UNIT_TESTS_SRCS
linalg/test_complex_operator.cpp
linalg/test_constrainedsolver.cpp
linalg/test_direct_solvers.cpp
linalg/test_fdsolver.cpp
linalg/test_hypre_ilu.cpp
linalg/test_ilu.cpp
linalg/test_matrix_block.cpp
+38 -4
View File
@@ -11,14 +11,48 @@
#define CATCH_CONFIG_RUNNER
#include "mfem.hpp"
#include "run_unit_tests.hpp"
#include "unit_tests.hpp"
bool launch_all_non_regression_tests = false;
std::string mfem_data_dir;
int main(int argc, char *argv[])
{
mfem::Device device("cuda");
// Include only tests labeled with CUDA. Exclude parallel tests.
return RunCatchSession(argc, argv, {"[CUDA]", "~[Parallel]"});
// There must be exactly one instance.
Catch::Session session;
// Build a new command line parser on top of Catch's
using namespace Catch::clara;
auto cli = session.cli() |
Opt(launch_all_non_regression_tests) ["--all"] ("all tests");
session.cli(cli);
// For floating point comparisons, print 8 digits for single precision
// values, and 16 digits for double precision values.
Catch::StringMaker<float>::precision = 8;
Catch::StringMaker<double>::precision = 16;
// Apply provided command line arguments.
int r = session.applyCommandLine(argc, argv);
if (r != 0) { return r; }
auto cfg = session.configData();
cfg.testsOrTags.push_back("[CUDA]");
#ifdef MFEM_USE_MPI
// Exclude tests marked as Parallel in a serial run, even when compiled with
// MPI. This is done because there is no MPI session initialized.
cfg.testsOrTags.push_back("~[Parallel]");
#endif
std::cout << "INFO: Test filter: [CUDA] ~[Parallel]" << std::endl;
device.Print();
session.useConfigData(cfg);
int result = session.run();
return result;
}
-84
View File
@@ -11,7 +11,6 @@
#include "mfem.hpp"
#include "unit_tests.hpp"
#include "general/tinyxml2.h"
#include <stdio.h>
#ifndef _WIN32
@@ -236,86 +235,3 @@ TEST_CASE("Save and load from collections", "[DataCollection]")
}
}
void SaveDataCollection(DataCollection &dc, int cycle, double t)
{
dc.SetCycle(cycle);
dc.SetTime(t);
dc.Save();
}
TEST_CASE("ParaView restart mode", "[ParaView]")
{
Mesh mesh = Mesh::MakeCartesian2D(2, 3, Element::QUADRILATERAL);
H1_FECollection fec(1, mesh.Dimension());
FiniteElementSpace fes(&mesh, &fec);
GridFunction u(&fes);
u = 0.0;
// Write initial dataset with three timesteps: 0, 1, 2.
{
ParaViewDataCollection dc("ParaView", &mesh);
dc.RegisterField("u", &u);
SaveDataCollection(dc, 0, 0);
SaveDataCollection(dc, 1, 1);
SaveDataCollection(dc, 2, 2);
}
// Using restart mode, append to the existing dataset, overwriting timesteps
// 1 and 2 with 1 and 1.5.
{
ParaViewDataCollection dc("ParaView", &mesh);
dc.UseRestartMode(true);
dc.RegisterField("u", &u);
SaveDataCollection(dc, 1, 1.0);
SaveDataCollection(dc, 2, 1.5);
}
// Parse the resulting PVD file, and verify that the structure is correct,
// and that it contains three timesteps: 0, 1, and 1.5.
using namespace tinyxml2;
auto StringCompare = [](const char *s1, const char *s2)
{
if (s1 == NULL || s2 == NULL) { return false; }
return strcmp(s1, s2) == 0;
};
auto VerifyDataset = [StringCompare](const XMLElement *ds, double t_ref)
{
REQUIRE(ds);
REQUIRE(StringCompare(ds->Name(), "DataSet"));
const char *timestep = ds->Attribute("timestep");
REQUIRE(timestep);
double t = std::stod(timestep);
REQUIRE(t == MFEM_Approx(t_ref));
};
XMLDocument xml;
xml.LoadFile("ParaView/ParaView.pvd");
REQUIRE(xml.ErrorID() == XML_SUCCESS);
const XMLElement *vtkfile = xml.FirstChildElement();
REQUIRE(vtkfile);
REQUIRE(StringCompare(vtkfile->Name(), "VTKFile"));
const XMLElement *collection = vtkfile->FirstChildElement();
REQUIRE(collection);
REQUIRE(StringCompare(collection->Name(), "Collection"));
const XMLElement *dataset = collection->FirstChildElement();
VerifyDataset(dataset, 0.0);
dataset = dataset->NextSiblingElement();
VerifyDataset(dataset, 1.0);
dataset = dataset->NextSiblingElement();
VerifyDataset(dataset, 1.5);
REQUIRE(dataset->NextSiblingElement() == NULL);
// Clean up
for (int c=0; c<=2; ++c)
{
std::string prefix = "ParaView/Cycle00000" + std::to_string(c);
REQUIRE(remove((prefix + "/data.pvtu").c_str()) == 0);
REQUIRE(remove((prefix + "/proc000000.vtu").c_str()) == 0);
REQUIRE(rmdir(prefix.c_str()) == 0);
}
REQUIRE(remove("ParaView/ParaView.pvd") == 0);
REQUIRE(rmdir("ParaView") == 0);
}
+7 -7
View File
@@ -167,9 +167,9 @@ TEST_CASE("DG SumIntegrator", "[SumIntegrator][PartialAssembly]")
integ2.AssemblePAInteriorFaces(fes);
integ_sum.AssemblePAInteriorFaces(fes);
const FaceRestriction *R_int = fes.GetFaceRestriction(
ElementDofOrdering::LEXICOGRAPHIC,
FaceType::Interior);
const Operator *R_int = fes.GetFaceRestriction(
ElementDofOrdering::LEXICOGRAPHIC,
FaceType::Interior);
int n_int = R_int->Height();
Vector x(n_int), y1(n_int), y2(n_int);
@@ -198,10 +198,10 @@ TEST_CASE("DG SumIntegrator", "[SumIntegrator][PartialAssembly]")
integ2.AssemblePABoundaryFaces(fes);
integ_sum.AssemblePABoundaryFaces(fes);
const FaceRestriction *R_bdr = fes.GetFaceRestriction(
ElementDofOrdering::LEXICOGRAPHIC,
FaceType::Boundary,
L2FaceValues::DoubleValued);
const Operator *R_bdr = fes.GetFaceRestriction(
ElementDofOrdering::LEXICOGRAPHIC,
FaceType::Boundary,
L2FaceValues::DoubleValued);
int n_bdr = R_bdr->Height();
x.SetSize(n_bdr);
+1 -1
View File
@@ -38,7 +38,7 @@ TEST_CASE("OperatorChebyshevSmoother", "[Chebyshev symmetry]")
Vector diag(fespace.GetTrueVSize());
aform.AssembleDiagonal(diag);
Solver* smoother = new OperatorChebyshevSmoother(*opr, diag, ess_tdof_list,
Solver* smoother = new OperatorChebyshevSmoother(opr.Ptr(), diag, ess_tdof_list,
cheb_order);
int n = smoother->Width();
-126
View File
@@ -1,126 +0,0 @@
// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
// LICENSE and NOTICE for details. LLNL-CODE-806117.
//
// This file is part of the MFEM library. For more information and source code
// availability visit https://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the BSD-3 license. We welcome feedback and contributions, see file
// CONTRIBUTING.md for details.
#include "mfem.hpp"
#include "unit_tests.hpp"
using namespace mfem;
#ifdef MFEM_USE_LAPACK
TEST_CASE("FDSolver",
"[FDSolver]")
{
double tol = 1e-10;
// SPD matrices
DenseMatrix A0(
{
{
1.29919, 0.61256, 0.82545
},
{0.61256, 0.57891, 0.39662},
{0.82545, 0.39662, 0.57541}
});
DenseMatrix A1(
{
{0.748236, 0.701663, 0.607517, 0.236740},
{0.701663, 0.809316, 0.713186, 0.256070},
{0.607517, 0.713186, 0.794221, 0.233943},
{0.236740, 0.256070, 0.233943, 0.083129}
});
DenseMatrix B0(
{
{0.13483, 0.51389, 0.43052},
{0.51389, 2.26750, 1.86331},
{0.43052, 1.86331, 1.59869}
});
DenseMatrix B1(
{
{0.94177, 1.02400, 1.14743, 0.35723},
{1.02400, 1.79087, 1.78708, 0.78304},
{1.14743, 1.78708, 2.06259, 0.80837},
{0.35723, 0.78304, 0.80837, 1.01798}
});
SECTION("2D")
{
Array<DenseMatrix *> A(2), B(2);
A[0] = &A0; A[1] = &A1;
B[0] = &B0; B[1] = &B1;
Vector y(12); y.Randomize(1);
Vector x(12), diff(12);
FDSolver S(A,B);
S.Mult(y,x);
DenseMatrix C1, C;
KronProd(A0, B1, C1);
KronProd(B0, A1, C);
C.Add(1., C1);
DenseMatrixInverse Cinv(C);
Cinv.Mult(y,diff);
diff-=x;
REQUIRE(diff.Norml2() < tol);
}
SECTION("3D")
{
DenseMatrix A2(
{
{1.14593, 0.76119},
{0.76119, 0.78993}
});
DenseMatrix B2(
{
{0.88088, 0.37899},
{0.37899, 0.45096}
});
Array<DenseMatrix *> A(3), B(3);
A[0] = &A0; A[1] = &A1; A[2] = &A2;
B[0] = &B0; B[1] = &B1; B[2] = &B2;
Vector y(24); y.Randomize(1);
Vector x(24), diff(24);
FDSolver S(A,B);
S.Mult(y,x);
DenseMatrix Temp, C0, C1, C;
KronProd(A0, B1, Temp);
KronProd(Temp, B2, C0);
KronProd(B0, A1, Temp);
KronProd(Temp, B2, C1);
KronProd(B0, B1, Temp);
KronProd(Temp, A2, C);
C.Add(1.,C0);
C.Add(1.,C1);
DenseMatrixInverse Cinv(C);
Cinv.Mult(y,diff);
diff-=x;
REQUIRE(diff.Norml2() < tol);
}
}
#endif // if MFEM_USE_LAPACK
-290
View File
@@ -12,8 +12,6 @@
#include "mfem.hpp"
#include "unit_tests.hpp"
#include "linalg/dtensor.hpp"
#include "nvToolsExt.h"
#include "nvToolsExtCudaRt.h"
using namespace mfem;
@@ -241,212 +239,6 @@ TEST_CASE("DenseMatrix A*B^T methods",
}
}
TEST_CASE("KronMult methods",
"[DenseMatrix], [CUDA]")
{
double tol = 1e-12;
int nA = 3, mA = 4;
int nB = 5, mB = 6;
DenseMatrix A(nA,mA);
DenseMatrix B(nB,mB);
for (int i = 0; i<nA; i++)
for (int j = 0; j<mA; j++)
{
A(i,j) = ((double)rand()/(double)RAND_MAX);
}
for (int i = 0; i<nB; i++)
for (int j = 0; j<mB; j++)
{
B(i,j) = ((double)rand()/(double)RAND_MAX);
}
DenseMatrix AB;
KronProd(A,B,AB);
AB.HostRead();
// (A ⊗ B) r
SECTION("KronMultABr")
{
nvtxRangePush("KronMultABr");
Vector r(mA*mB);
MFEM_VERIFY(r.Size() == AB.Width(), "Check r size");
r.HostReadWrite();
r.Randomize();
Vector z0(AB.Height());
AB.Mult(r,z0);
//z0.HostRead();
Vector z1;
KronMult(A,B,r,z1);
MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
z0-=z1;
REQUIRE(z0.Norml2() < tol);
nvtxRangePop();
}
// (A ⊗ B) R
SECTION("KronMultABR")
{
nvtxRangePush("KronMultABR");
int nR = mA*mB;
int mR = 7;
DenseMatrix R(nR, mR);
R.HostReadWrite();
for (int i = 0; i<nR; i++)
for (int j = 0; j<mR; j++)
{
R(i,j) = ((double)rand()/(double)RAND_MAX);
}
DenseMatrix Z0(nA*nB,mR);
Mult(AB,R,Z0);
DenseMatrix Z1;
KronMult(A,B,R,Z1);
MFEM_VERIFY(Z0.Height() == Z1.Height() &&
Z0.Width() == Z1.Width(), "Check z1 size");
Z0-=Z1;
REQUIRE(Z0.MaxMaxNorm() < tol);
nvtxRangePop();
}
// (A ⊗ B ⊗ C) r
SECTION("KronMultABCr")
{
nvtxRangePush("KronMultABCr");
int nC = 7, mC = 2;
DenseMatrix C(nC, mC);
C.HostReadWrite();
for (int i = 0; i<nC; i++)
for (int j = 0; j<mC; j++)
{
C(i,j) = ((double)rand()/(double)RAND_MAX);
}
DenseMatrix ABC;
KronProd(AB,C,ABC);
Vector r(mA*mB*mC);
r.HostReadWrite();
r.Randomize();
MFEM_VERIFY(r.Size() == ABC.Width(), "Check r size");
Vector z0(nA*nB*nC);
ABC.Mult(r,z0);
Vector z1;
KronMult(A,B,C,r,z1);
MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
z0-=z1;
REQUIRE(z0.Norml2() < tol);
nvtxRangePop();
}
}
TEST_CASE("KronMultInv methods",
"[DenseMatrixInverse]")
{
double tol = 1e-12;
int nA = 3;
int nB = 2;
DenseMatrix A(
{
{ 1.0, 0.2, 3.4},
{-2.0, -1.0, 3.1},
{ 0.7, 1.4,-0.9}
});
DenseMatrix B(
{
{-10.1, 5.7},
{-3.0, 4.2}
});
DenseMatrixInverse Ainv(A);
DenseMatrixInverse Binv(B);
DenseMatrix AB;
KronProd(A,B,AB);
// (A^-1 ⊗ B^-1) r
SECTION("KronMultInvABr")
{
Vector r(nA*nB); r.Randomize();
MFEM_VERIFY(r.Size() == AB.Width(), "Check r size");
Vector z0(AB.Height());
DenseMatrixInverse ABinv(AB);
ABinv.Mult(r,z0);
Vector z1;
KronMult(Ainv,Binv,r,z1);
MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
z0-=z1;
REQUIRE(z0.Norml2() < tol);
}
// (A^-1 ⊗ B^-1) R
SECTION("KronMultInvABR")
{
int nR = nA*nB;
int mR = 7;
DenseMatrix R(nR, mR);
for (int i = 0; i<nR; i++)
for (int j = 0; j<mR; j++)
{
R(i,j) = ((double)rand()/(double)RAND_MAX);
}
DenseMatrixInverse ABinv(AB);
DenseMatrix Z0(nA*nB,mR);
ABinv.Mult(R,Z0);
DenseMatrix Z1;
KronMult(Ainv,Binv,R,Z1);
MFEM_VERIFY(Z0.Height() == Z1.Height() &&
Z0.Width() == Z1.Width(), "Check z1 size");
Z0-=Z1;
REQUIRE(Z0.MaxMaxNorm() < tol);
}
// (A^-1 ⊗ B^-1 ⊗ C^-1) r
SECTION("KronMultInvABCr")
{
int nC = 4;
DenseMatrix C(
{
{-2.1, 1.6, -3.4, 17.5},
{-7.1, 1.3, -7.5, -12.5},
{ 0.5, 5.7, -6.0, -0.5},
{ 9.2, 0.3, -1.4, -14.9}
});
DenseMatrix ABC;
KronProd(AB,C,ABC);
DenseMatrixInverse ABCInv(ABC);
Vector r(nA*nB*nC); r.Randomize();
MFEM_VERIFY(r.Size() == ABC.Width(), "Check r size");
Vector z0(nA*nB*nC);
ABCInv.Mult(r,z0);
DenseMatrixInverse Cinv(C);
Vector z1;
KronMult(Ainv,Binv,Cinv,r,z1);
MFEM_VERIFY(z0.Size() == z1.Size(), "Check z1 size");
z0-=z1;
REQUIRE(z0.Norml2() < tol);
}
}
TEST_CASE("LUFactors RightSolve", "[DenseMatrix]")
{
@@ -522,88 +314,6 @@ TEST_CASE("DenseTensor LinearSolve methods",
}
}
#ifdef MFEM_USE_LAPACK
TEST_CASE("EigenSystem methods",
"[DenseMatrix]")
{
double tol = 1e-12;
SECTION("SPD Matrix")
{
DenseMatrix A({{0.56806, 0.29211, 0.48315, 0.70024},
{0.29211, 0.85147, 0.68123, 0.70689},
{0.48315, 0.68123, 1.07229, 1.02681},
{0.70024, 0.70689, 1.02681, 1.15468}
});
DenseMatrix V, AV(4);
Vector Lambda;
for (bool sym: { false, true })
{
DenseMatrixEigensystem EigA(A,sym);
EigA.Eval();
V = EigA.Eigenvectors();
Lambda = EigA.Eigenvalues();
Mult(A,V,AV);
V.RightScaling(Lambda);
AV -= V;
REQUIRE(AV.MaxMaxNorm() < tol);
}
}
SECTION("Indefinite Matrix")
{
DenseMatrix A({{0.486278, 0.041135, 0.480727, 0.616026},
{0.523599, 0.119827, 0.087808, 0.415241},
{0.214454, 0.661631, 0.909626, 0.744259},
{0.107007, 0.630604, 0.077862, 0.221006}
});
DenseMatrixEigensystem EigA(A);
EigA.Eval();
Vector Lambda_r, Lambda_i;
// Real part of eigenvalues
Lambda_r = EigA.Eigenvalues();
// Imag part of eigenvalues
Lambda_i = EigA.Eigenvalues(true);
DenseMatrix V;
V = EigA.Eigenvectors();
// Real part of eigenvectors
DenseMatrix Vr(4), Vi(4);
Vr.SetCol(0,V.GetColumn(0));
Vr.SetCol(1,V.GetColumn(1));
Vr.SetCol(2,V.GetColumn(1));
Vr.SetCol(3,V.GetColumn(3));
// Imag part of eigenvectors
Vector vi(4); V.GetColumn(2,vi);
Vi.SetCol(0,0.);
Vi.SetCol(1,vi); vi *= -1.;
Vi.SetCol(2,vi);
Vi.SetCol(3,0.);
// Check that A*V = V * Lambda
// or A * (V_r + i V_i ) = (V_r + i V_i)*(Lamda_r + i Lambda_i)
// or A * V_r = V_r * Lambda_r - V_i * Lambda_i
// and A * V_i = V_r ( Lambda_i + V_i * Lambda_r
DenseMatrix AVr(4), AVi(4);
Mult(A,Vr, AVr);
Mult(A,Vi, AVi);
DenseMatrix Vrlr = Vr; Vrlr.RightScaling(Lambda_r);
DenseMatrix Vrli = Vr; Vrli.RightScaling(Lambda_i);
DenseMatrix Vilr = Vi; Vilr.RightScaling(Lambda_r);
DenseMatrix Vili = Vi; Vili.RightScaling(Lambda_i);
AVr -= Vrlr; AVr+= Vili;
AVi -= Vrli; AVi-= Vilr;
REQUIRE(AVr.MaxMaxNorm() < tol);
REQUIRE(AVi.MaxMaxNorm() < tol);
}
}
#endif // if MFEM_USE_LAPACK
TEST_CASE("DenseTensor copy", "[DenseMatrix][DenseTensor]")
{
DenseTensor t1(2,3,4);
+8 -13
View File
@@ -214,27 +214,21 @@ $(DATA_DIR):
MFEM_TESTS = UNIT_TESTS
include $(MFEM_TEST_MK)
ifeq (,$(wildcard $(MFEM_DIR)/../data))
MFEM_DATA_FLAG =
else
MFEM_DATA_FLAG = --data $(MFEM_DIR)/../data
endif
%-test-seq: %
@$(call mfem-test,$<,, Unit tests,$(MFEM_DATA_FLAG),SKIP-NO-VIS)
@$(call mfem-test,$<,, Unit tests,,SKIP-NO-VIS)
ceed_tests-test-seq: ceed_tests
@$(call mfem-test,$<,, CEED Unit tests (cpu),--device ceed-cpu $(MFEM_DATA_FLAG),SKIP-NO-VIS)
@$(call mfem-test,$<,, CEED Unit tests (cpu),--device ceed-cpu,SKIP-NO-VIS)
ifeq ($(MFEM_USE_CUDA),YES)
@$(call mfem-test,$<,, CEED Unit tests (cuda-ref),--device ceed-cuda:/gpu/cuda/ref $(MFEM_DATA_FLAG),SKIP-NO-VIS)
@$(call mfem-test,$<,, CEED Unit tests (cuda-shared),--device ceed-cuda:/gpu/cuda/shared $(MFEM_DATA_FLAG),SKIP-NO-VIS)
@$(call mfem-test,$<,, CEED Unit tests (cuda-gen),--device ceed-cuda:/gpu/cuda/gen $(MFEM_DATA_FLAG),SKIP-NO-VIS)
@$(call mfem-test,$<,, CEED Unit tests (cuda-ref),--device ceed-cuda:/gpu/cuda/ref,SKIP-NO-VIS)
@$(call mfem-test,$<,, CEED Unit tests (cuda-shared),--device ceed-cuda:/gpu/cuda/shared,SKIP-NO-VIS)
@$(call mfem-test,$<,, CEED Unit tests (cuda-gen),--device ceed-cuda:/gpu/cuda/gen,SKIP-NO-VIS)
endif
RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP)
%-test-par: %
@$(call mfem-test,$<, $(RUN_MPI) 1, Parallel unit tests,$(MFEM_DATA_FLAG),SKIP-NO-VIS)
@$(call mfem-test,$<, $(RUN_MPI) $(MFEM_MPI_NP), Parallel unit tests,$(MFEM_DATA_FLAG),SKIP-NO-VIS)
@$(call mfem-test,$<, $(RUN_MPI) 1, Parallel unit tests,,SKIP-NO-VIS)
@$(call mfem-test,$<, $(RUN_MPI) $(MFEM_MPI_NP), Parallel unit tests,,SKIP-NO-VIS)
# Generate an error message if the MFEM library is not built and exit
$(MFEM_LIB_FILE):
@@ -243,3 +237,4 @@ $(MFEM_LIB_FILE):
clean:
rm -f $(SEQ_UNIT_TESTS) $(PAR_UNIT_TESTS) *.o */*.o */*~ *~
rm -rf *.dSYM output_meshes
-18
View File
@@ -41,21 +41,3 @@ TEST_CASE("VTU XML Reader", "[Mesh][VTU][XML]")
REQUIRE(mesh.GetNumGeometries(2) == 1);
}
}
TEST_CASE("VTU XML Compressed Blocks", "[VTU][XML][MFEMData]")
{
auto filename = GENERATE(
"bracket_appended_compressed.vtu",
"bracket_appended_encoded_compressed.vtu",
"bracket_inline_compressed.vtu"
);
std::string mesh_path = mfem_data_dir + "/vtk/" + filename;
Mesh mesh = Mesh::LoadFromFile(mesh_path.c_str());
REQUIRE(mesh.Dimension() == 3);
REQUIRE(mesh.GetNE() == 206208);
REQUIRE(mesh.GetNV() == 50000);
REQUIRE(mesh.HasGeometry(Geometry::TETRAHEDRON));
REQUIRE(mesh.GetNumGeometries(3) == 1);
}
-15
View File
@@ -186,12 +186,10 @@ int tmop(int id, Req &res, int argc, char *argv[])
case 2: metric = new TMOP_Metric_002; break;
case 7: metric = new TMOP_Metric_007; break;
case 77: metric = new TMOP_Metric_077; break;
case 80: metric = new TMOP_Metric_080(0.5); break;
case 302: metric = new TMOP_Metric_302; break;
case 303: metric = new TMOP_Metric_303; break;
case 315: metric = new TMOP_Metric_315; break;
case 321: metric = new TMOP_Metric_321; break;
case 332: metric = new TMOP_Metric_332(0.5); break;
default:
{
if (id == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
@@ -742,13 +740,6 @@ static void tmop_tests(int id = 0, bool all = false)
POR({1,2}).QOR({2,4}).
TID({4}).MID({1,2})).Run(id,all);
Launch(Launch::Args("Square01 + Adapted discrete size").
MESH("../../miniapps/meshing/square01.mesh").REFINE(1).
NORMALIZATION(true).
POR({1,2}).QOR({4,6}).
LINEAR_ITERATIONS(150).
TID({5}).MID({80}).LS({3})).Run(id,all);
Launch(Launch::Args("Blade").
MESH("../../miniapps/meshing/blade.mesh").
POR({1,2}).QOR({2,4}).
@@ -777,12 +768,6 @@ static void tmop_tests(int id = 0, bool all = false)
POR({1,2}).QOR({4,2}).
TID({7}).MID({302,321})).Run(id,all);
Launch(Launch::Args("Cube + Discrete size + normalization").
MESH("../../miniapps/meshing/cube.mesh").
NORMALIZATION(true).
POR({1,2}).QOR({4,2}).
TID({5}).MID({332})).Run(id,all);
// Note: order 1 has no interior nodes, so all residuals are zero and the
// Newton iteration exits immediately.
Launch(Launch::Args("Toroid-Hex").
+36 -7
View File
@@ -11,10 +11,9 @@
#define CATCH_CONFIG_RUNNER
#include "mfem.hpp"
#include "run_unit_tests.hpp"
#include "unit_tests.hpp"
bool launch_all_non_regression_tests = false;
std::string mfem_data_dir;
#ifdef MFEM_USE_MPI
mfem::MPI_Session *GlobalMPISession;
@@ -25,14 +24,44 @@ mfem::MPI_Session *GlobalMPISession;
int main(int argc, char *argv[])
{
mfem::Device device("cuda");
// There must be exactly one instance.
Catch::Session session;
// Build a new command line parser on top of Catch's
using namespace Catch::clara;
auto cli = session.cli() |
Opt(launch_all_non_regression_tests) ["--all"] ("all tests");
session.cli(cli);
// For floating point comparisons, print 8 digits for single precision
// values, and 16 digits for double precision values.
Catch::StringMaker<float>::precision = 8;
Catch::StringMaker<double>::precision = 16;
// Apply provided command line arguments.
int r = session.applyCommandLine(argc, argv);
if (r != 0) { return r; }
#ifdef MFEM_USE_MPI
mfem::MPI_Session mpi;
GlobalMPISession = &mpi;
bool root = mpi.Root();
#else
bool root = true;
// Exclude all tests that are not labeled with Parallel and CUDA.
auto cfg = session.configData();
cfg.testsOrTags.push_back("[Parallel]");
cfg.testsOrTags.push_back("[CUDA]");
session.useConfigData(cfg);
if (mpi.Root())
{
std::cout << "INFO: Test filter: [Parallel] [CUDA]" << std::endl;
device.Print();
}
#endif
// Include only tests that are labeled with both CUDA and Parallel.
return RunCatchSession(argc, argv, {"[CUDA]","[Parallel]"}, root);
int result = session.run();
return result;
}
+35 -7
View File
@@ -11,10 +11,9 @@
#define CATCH_CONFIG_RUNNER
#include "mfem.hpp"
#include "run_unit_tests.hpp"
#include "unit_tests.hpp"
bool launch_all_non_regression_tests = false;
std::string mfem_data_dir;
#ifdef MFEM_USE_MPI
mfem::MPI_Session *GlobalMPISession;
@@ -24,14 +23,43 @@ mfem::MPI_Session *GlobalMPISession;
int main(int argc, char *argv[])
{
// There must be exactly one instance.
Catch::Session session;
// Build a new command line parser on top of Catch's
using namespace Catch::clara;
auto cli = session.cli() |
Opt(launch_all_non_regression_tests) ["--all"] ("all tests");
session.cli(cli);
// For floating point comparisons, print 8 digits for single precision
// values, and 16 digits for double precision values.
Catch::StringMaker<float>::precision = 8;
Catch::StringMaker<double>::precision = 16;
// Apply provided command line arguments.
int r = session.applyCommandLine(argc, argv);
if (r != 0) { return r; }
#ifdef MFEM_USE_MPI
mfem::MPI_Session mpi;
GlobalMPISession = &mpi;
bool root = mpi.Root();
#else
bool root = true;
// Exclude all tests that are not labeled with Parallel.
auto cfg = session.configData();
cfg.testsOrTags.push_back("[Parallel]");
session.useConfigData(cfg);
// NOTE: tests marked with "[CUDA]" (in addition to "[Parallel]") are still
// run with the default device.
if (mpi.Root())
{
std::cout << "INFO: Test filter: [Parallel]" << std::endl;
}
#endif
// Only run tests that are labeled with Parallel.
return RunCatchSession(argc, argv, {"[Parallel]"}, root);
int result = session.run();
return result;
}
-60
View File
@@ -1,60 +0,0 @@
// Copyright (c) 2010-2021, 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_RUN_UNIT_TEST
#define MFEM_RUN_UNIT_TEST
#include "unit_tests.hpp"
static int RunCatchSession(int argc, char *argv[],
const std::vector<std::string> &testsOrTags,
bool root=true)
{
// There must be exactly one instance.
Catch::Session session;
// Build a new command line parser on top of Catch's
using namespace Catch::clara;
auto cli = session.cli()
| Opt(launch_all_non_regression_tests) ["--all"] ("all tests")
| Opt(mfem_data_dir, "") ["--data"] ("mfem/data repository");
session.cli(cli);
// For floating point comparisons, print 8 digits for single precision
// values, and 16 digits for double precision values.
Catch::StringMaker<float>::precision = 8;
Catch::StringMaker<double>::precision = 16;
// Apply provided command line arguments.
int r = session.applyCommandLine(argc, argv);
if (r != 0) { return r; }
auto cfg = session.configData();
cfg.testsOrTags.insert(cfg.testsOrTags.end(), testsOrTags.begin(), testsOrTags.end());
if (mfem_data_dir == "") { cfg.testsOrTags.push_back("~[MFEMData]"); }
session.useConfigData(cfg);
if (root)
{
std::cout << "INFO: Test filter: ";
for (std::string &filter : cfg.testsOrTags)
{
std::cout << filter << " ";
}
std::cout << std::endl;
}
int result = session.run();
return result;
}
#endif
+34 -4
View File
@@ -11,13 +11,43 @@
#define CATCH_CONFIG_RUNNER
#include "mfem.hpp"
#include "run_unit_tests.hpp"
#include "unit_tests.hpp"
bool launch_all_non_regression_tests = false;
std::string mfem_data_dir;
int main(int argc, char *argv[])
{
// Exclude parallel tests.
return RunCatchSession(argc, argv, {"~[Parallel]"});
// There must be exactly one instance.
Catch::Session session;
// Build a new command line parser on top of Catch's
using namespace Catch::clara;
auto cli = session.cli() |
Opt(launch_all_non_regression_tests) ["--all"] ("all tests");
session.cli(cli);
// For floating point comparisons, print 8 digits for single precision
// values, and 16 digits for double precision values.
Catch::StringMaker<float>::precision = 8;
Catch::StringMaker<double>::precision = 16;
// Apply provided command line arguments.
int r = session.applyCommandLine(argc, argv);
if (r != 0) { return r; }
#ifdef MFEM_USE_MPI
// Exclude tests marked as Parallel in a serial run, even when compiled with
// MPI. This is done because there is no MPI session initialized.
auto cfg = session.configData();
cfg.testsOrTags.push_back("~[Parallel]");
session.useConfigData(cfg);
#endif
// NOTE: tests marked with "[CUDA]" are still run using the default device.
std::cout << "INFO: Test filter: ~[Parallel]" << std::endl;
int result = session.run();
return result;
}
-5
View File
@@ -17,11 +17,6 @@
/// Command line '--all' option to launch all non-regression tests.
extern bool launch_all_non_regression_tests;
/// Command line '--data' argument for path to mfem/data repo.
/** If no --data path is provided, then mfem_data_dir will be the empty string,
and tests tagged with [MFEMData] will be skipped. */
extern std::string mfem_data_dir;
/** @brief MFEM_Approx can be used to compare floating point values within an
absolute tolerance of @a abs_tol (default value 1e-12) and relative
tolerance of @a rel_tol (default value 1e-12). */