Compare commits

..
Author SHA1 Message Date
Stowell, Mark L 44fbd87bff Extending a unit test to include pyramids 2022-02-14 12:54:39 -08:00
Stowell, Mark L e456727573 Clean up 2022-02-01 20:14:38 -08:00
Stowell, Mark L ec8f2baeec Fixing the higher order integration rules for pyramids 2022-02-01 20:14:24 -08:00
Stowell, Mark L 47dbbf0869 Experimenting with arbitrary order H1 on pyramids 2022-01-31 17:22:19 -08:00
Tzanio Kolev bba2c08025 Merge pull request #2740 from mfem/vsx
vsx128 GNU fix [vsx]
2022-01-20 08:08:21 -08:00
Tzanio Kolev b23a3921bf Merge pull request #2743 from mfem/shadow-examples-dev
Fixing shadow variables in examples [shadow-examples-dev]
2022-01-19 09:56:40 -08:00
Will Pazner e4bcd7c58c Merge pull request #2746 from mfem/shadow-tests-dev
Fixing shadow variables in unit tests [shadow-tests-dev]
2022-01-18 15:08:30 -08:00
Veselin Dobrev 30a3c76325 Merge pull request #2771 from mfem/stopwatch-fix-comment
Fix erroneous comment in StopWatch
2022-01-18 12:12:35 -08:00
Stowell, Mark L 99e1a076ba Switching to ostream named os 2022-01-18 11:15:11 -08:00
Will Pazner 1ae61b0ada Merge remote-tracking branch 'origin/master' into shadow-tests-dev 2022-01-18 10:29:05 -08:00
Tzanio Kolev f0398ff391 Merge pull request #2732 from mfem/shadow-tmop-dev
Fixing shadowed variables in tmop files [shadow-tmop-dev]
2022-01-18 10:26:08 -08:00
Tzanio Kolev df1deca20e Merge pull request #2736 from mfem/shadow-gen-dev
Fixing shadow variables in general classes [shadow-gen-dev]
2022-01-17 17:53:37 -08:00
Tzanio Kolev d1b79a7fc3 Merge pull request #2744 from mfem/shadow-miniapps-dev
Fixing shadow variables in miniapps [shadow-miniapps-dev]
2022-01-17 17:53:00 -08:00
Tzanio Kolev 426a77be6c Merge pull request #2754 from mfem/shifted-bug-fix
Shifted bug fix
2022-01-17 17:52:17 -08:00
Will Pazner 20dced51fd Fix erroneous comment in StopWatch 2022-01-14 17:59:13 -08:00
Vladimir Z Tomov 00c4f59c76 Minor. 2022-01-13 18:47:38 -08:00
Ketan Mittal 5a6e0dad2e minor 2022-01-11 15:56:56 -08:00
Ketan Mittal 774794595b move temp_elvect->elvect outside the loop 2022-01-11 13:56:42 -08:00
camierjs cfa82dccc0 In device.?pp, use device_id for shadowed variable and switch out to os 2022-01-11 11:09:12 -08:00
camierjs 30803b8f00 Rename Memory<T>::New input arguments 2022-01-11 11:05:13 -08:00
Nabil Atallah 92112fccaa changed '=' sign to '+=' in line 644 otherwise at every integration point the values get overwritten 2022-01-05 16:24:51 -08:00
Nabil Atallah 70bd112ea4 added the sbm_aux.hpp file to the CMakeLists.txt and also added it to line 56 since the object files depend on it 2022-01-05 13:59:31 -08:00
Nabil Atallah 6ba7e82f91 Merge branch 'master' of https://github.com/mfem/mfem
Fixed a bug in the sbm_solver.cpp file of the shifted miniapp.
2022-01-05 13:38:44 -08:00
Nabil Atallah ba83609892 reset the temp_elvect vector to zero at the beginning of the loop over each integration point in sbm_solver.cpp belonging to the shifted miniapp, also fixed a typo in the comments related to the sign of the SBM penalty term 2022-01-05 11:28:49 -08:00
Stowell, Mark L 74fc2de695 Adjusting names in function declarations 2022-01-04 17:34:47 -08:00
Stowell, Mark L df6e73654a Fixing shadow variables in autodiff miniapps 2021-12-30 11:33:21 -08:00
Stowell, Mark L a07fc07206 Changes suggested by @camierjs 2021-12-30 11:04:44 -08:00
Stowell, Mark L 5e672f26bc Removing unneeded function after consulting @camierjs 2021-12-29 10:24:18 -08:00
Vladimir Z Tomov f3f37e8e5a Minor. 2021-12-28 15:19:24 -08:00
Stowell, Mark L 475d3d8102 Fixing shadow variables in unit tests 2021-12-28 11:42:51 -08:00
Veselin Dobrev 0f239c22e0 Merge pull request #2730 from mfem/ScaledOperatorFix
Fix size of ScaledOperator and add MultTranspose
2021-12-28 10:28:17 -08:00
Tzanio Kolev 783f0e0304 Merge pull request #2709 from mfem/new-compilers-warnings-fix
Fix some compiler warnings from new clang and gcc versions
2021-12-28 09:58:40 -08:00
Stowell, Mark L dd03eb6d79 Adding missing shadow vdofs 2021-12-27 14:38:20 -08:00
Stowell, Mark L 0f4bd5925a Fixing shadow variables in miniapps 2021-12-27 12:46:22 -08:00
Stowell, Mark L ca3e2b3098 make style 2021-12-27 10:42:01 -08:00
Stowell, Mark L b14bd1f237 Fixing shadow variables in examples 2021-12-27 10:41:09 -08:00
Stowell, Mark L fa25a70da6 Fixing shadow variables in general classes 2021-12-23 16:21:31 -08:00
camierjs e9b9152f59 vsx128 fix 2021-12-23 16:20:05 -08:00
Stowell, Mark L bd7527e922 Swapping leading underscores for trailing 2021-12-23 09:49:01 -08:00
Stowell, Mark L 60a1d3d1d4 Fixing shadowed variables in tmop files 2021-12-22 15:46:33 -08:00
Syun'ichi Shiraiwa 07ba5bd336 make style 2021-12-21 21:55:09 -05:00
Syun'ichi Shiraiwa 4027c86797 fix size of ScaledOperator and add MultTranspose 2021-12-21 21:52:54 -05:00
Tzanio Kolev 8300a563ae Merge pull request #2722 from mfem/pncmesh-getessvdof-fix
Fix `ParFiniteElementSpace::GetEssentialVDofs()` for non-conforming meshes
2021-12-21 08:36:14 -08:00
Tzanio Kolev 753f351934 Merge pull request #2519 from mfem/bicgstab-print-lvl-fix
Modernize print level semantics
2021-12-21 07:47:49 -08:00
Veselin Dobrev 9310e05a15 Override appropriately the SetPrintLevel methods in class KINSolver 2021-12-20 16:20:27 -08:00
termi-official 985fff881c Treat 'no convergence' as a warning. 2021-12-20 21:33:54 +01:00
termi-official 5b841fb239 Remove deprecation attribute from IterativeSolver::print_level together with all of its support code. 2021-12-20 20:53:59 +01:00
Tzanio Kolev 6f7b29ac14 Merge branch 'master' into bicgstab-print-lvl-fix 2021-12-19 17:58:19 -08:00
Tzanio Kolev 589d2775cb Merge pull request #2629 from mfem/simd-unary-plus
Add unary + to AutoSIMD
2021-12-19 12:14:36 -08:00
Tzanio Kolev 8a565cad67 Merge pull request #2677 from mfem/jacobi-abs-diag
Positive diagonal in Jacobi smoothers
2021-12-16 17:58:31 -08:00
Tzanio Kolev 5f9ee51d40 Merge pull request #2684 from mfem/tmop-renameandrefactor
TMOP clean up
2021-12-16 17:57:45 -08:00
Tzanio Kolev acb85daaff Merge pull request #2656 from mfem/lor-gridfunction-coeff
GridFunction coefficients that work with LOR preconditioning
2021-12-16 17:57:04 -08:00
Tzanio Kolev 36ec075849 Merge pull request #2702 from mfem/add-code-of-conduct
MFEM Code of Conduct
2021-12-16 14:12:31 -08:00
Veselin Dobrev 8b0262bd29 Fix an issue where some boundary dofs are not properly marked
as essential in ParFiniteElementSpace::GetEssentialVDofs() when
using a parallel non-conforming mesh.

This should address GitHub issue #2559.
2021-12-16 02:53:09 -08:00
Tzanio d9a18088a2 Updated to Contributor Covenant Code 2.1 2021-12-15 08:02:04 -08:00
Tzanio KolevandVeselin Dobrev 830ea90cf8 Update CODE_OF_CONDUCT.md
Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>
2021-12-15 07:28:56 -08:00
Tzanio KolevandVeselin Dobrev 9377a43d28 Update CODE_OF_CONDUCT.md
Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>
2021-12-15 07:28:51 -08:00
Tzanio KolevandVeselin Dobrev 87cd94a1fe Update CODE_OF_CONDUCT.md
Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>
2021-12-15 07:28:44 -08:00
Tzanio KolevandVeselin Dobrev 1aa1e0633b Update CODE_OF_CONDUCT.md
Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>
2021-12-15 07:28:36 -08:00
Tzanio KolevandVeselin Dobrev 57f3462a00 Update CODE_OF_CONDUCT.md
Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>
2021-12-15 07:28:31 -08:00
Tzanio KolevandVeselin Dobrev 186f65ccc1 Update CODE_OF_CONDUCT.md
Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>
2021-12-15 07:28:23 -08:00
Tzanio KolevandVeselin Dobrev 8e11743052 Update CODE_OF_CONDUCT.md
Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>
2021-12-15 07:28:16 -08:00
Tzanio KolevandVeselin Dobrev 649163a36f Update CODE_OF_CONDUCT.md
Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>
2021-12-15 07:28:09 -08:00
Veselin Dobrev a58567695e Fix an error when building with HIP which may also be causing
problems with CUDA.

In OperatorJacobiSmoother, use absolute value of diagonal, if
requested, during setup instead of during action.

A few tweaks in SparseMatrix::Jacobi and SparseMatrix::DiagScale.

For the GitLab CI on Lassen, disable ATS.
2021-12-14 22:37:53 -08:00
Tzanio Kolev 53c04ef171 Update CONTRIBUTING.md 2021-12-14 21:52:00 -08:00
Veselin Dobrev 2c09747196 Fix some compiler warnings from new clang and gcc versions 2021-12-14 21:49:31 -08:00
Tzanio e9feadbbfc Small updates 2021-12-14 08:27:27 -08:00
Tzanio 81dceb94b6 Adressing comments 2021-12-13 13:21:20 -08:00
Tzanio Kolev 39022bce0f Merge pull request #2698 from mfem/nvwarnings
Remove nvcc warnings about partially overridden virtual functions [nvwarnings]
2021-12-13 13:13:38 -08:00
Tzanio 15242ffcf6 Editorial adjustments in CHANGELOG. 2021-12-13 13:08:53 -08:00
Will Pazner 9a5cafafc6 Merge remote-tracking branch 'origin/master' into simd-unary-plus 2021-12-13 11:37:29 -08:00
Will Pazner 11275857a3 Update CHANGELOG 2021-12-13 11:37:18 -08:00
Will Pazner d06a528958 Merge remote-tracking branch 'origin/master' into lor-gridfunction-coeff 2021-12-13 10:49:50 -08:00
Will Pazner f54b8bd6d8 Update CHANGELOG to mention LOR GridFunction coefficient support 2021-12-13 10:49:39 -08:00
Will Pazner 239c672988 Merge pull request #2655 from mfem/pmesh-mem-leak-fix
Fix Memory leaks in PMesh
2021-12-13 10:40:51 -08:00
Vladimir Z Tomov aabf4ce84d Fixed a wrong coefficient in the adaptive limiting. 2021-12-12 16:27:40 -08:00
Tzanio Kolev f95c2e156b Create CODE_OF_CONDUCT.md 2021-12-10 12:11:10 -08:00
camierjs 80fa63cb64 Remove nvcc warnings about partially overridden virtual functions 2021-12-08 10:46:29 -08:00
Veselin Dobrev 97796f0beb Merge branch 'master' into bicgstab-print-lvl-fix
Resolved conflict:
   CHANGELOG
2021-12-07 12:44:31 -08:00
Will Pazner c1562552f1 Make performance ex1 and ex1p templated on dimension 2021-12-02 10:17:06 -08:00
termi-official a44a8640d1 Change output behavior of warnings to mfem::out for consistency and update PrintLevel docs with the associated stream. 2021-12-01 17:24:26 +01:00
Ketan Mittal 001f1a8b79 add missing PC 3D metrics 2021-11-30 13:59:26 -08:00
Ketan Mittal 92e42d4332 fix spacing etc 2021-11-30 12:45:19 -08:00
Ketan Mittal 9c7150c93e Merge branch 'master' of https://github.com/mfem/mfem into tmop-renameandrefactor 2021-11-30 10:24:14 -08:00
Ketan Mittal 19a2ad26e1 minor 2021-11-30 10:24:12 -08:00
Veselin Dobrev 568562a7d5 Restore a comment that was deleted unintentionallly. 2021-11-26 22:14:42 -08:00
Veselin Dobrev aa8be1b59c In class IterativeSolver, silence the depracation warnings for
`print_level` in a way that works for both Clang and GCC.
2021-11-26 21:34:33 -08:00
Veselin Dobrev 85a79cbea5 Try to suppress deprecation warning with GCC from the
`IterativeSolver` constructors.
2021-11-26 20:17:19 -08:00
Veselin Dobrev 3bf7e07457 Tweak the Doxygen documentation for class IterativeSolver.
Suppress deprecation warnings only when modifying the `print_level`
data member of class `IterativeSolver`.
2021-11-26 20:05:44 -08:00
Vladimir Z Tomov c94cd73dd1 empty line. 2021-11-26 16:02:27 -08:00
Vladimir Z Tomov ff427a04a6 Minor. 2021-11-26 16:01:41 -08:00
Vladimir Z Tomov 2d10dd0abe Options to use abs values of the diagonal in OperatorJacobi, DSmoother. 2021-11-26 15:49:40 -08:00
termi-official a53353776d Update convergence criterion docs. 2021-11-26 17:08:39 +01:00
termi-official 2f1c03796a Adapt tmop print policy. 2021-11-26 16:39:06 +01:00
termi-official 199d19c195 Add documentation about warning suppression in solver. 2021-11-26 16:28:32 +01:00
Will Pazner eb6ef4a537 Minor 2021-11-24 17:47:27 -08:00
Will Pazner af7c58aa0e Minor solver output consistency fixes 2021-11-24 17:20:31 -08:00
Will Pazner d4f404cf74 Don't repeat last iteration in GMRES output 2021-11-24 10:53:53 -08:00
Will Pazner 2080fde2a1 Small changes to reduce diffs with old solver output 2021-11-24 10:44:38 -08:00
Will Pazner 7ff5874904 Typo 2021-11-21 13:30:50 -08:00
Will Pazner 757ee1a24a Merge remote-tracking branch 'origin/master' into lor-gridfunction-coeff 2021-11-20 13:13:17 -08:00
Will Pazner 861c7c4a04 Add unit tests for GridFunction coefficients on refined meshes 2021-11-20 13:13:16 -08:00
Will Pazner 60b5110031 Allow the same GridFunctionCoefficient (and related classes) to work on both coarse and refined meshes
Useful for LOR preconditioning with GridFunction coefficients. The same
Coefficient can be used on both the coarse and refined mesh.
2021-11-09 16:53:31 -08:00
Will Pazner fe08c6bd36 Add mesh data member to ElementTransformation 2021-11-09 15:55:36 -08:00
Tom Stitt 977e978ad7 fix one more leak 2021-11-09 15:50:22 -08:00
Tom Stitt db43873b7d delete face_nbr_el_to_face in ParMesh::Destroy 2021-11-09 13:43:48 -08:00
Will Pazner 40f672c230 Add unary + SIMD specializations 2021-11-02 13:33:38 -07:00
termi-official e49a744288 Apply style... 2021-11-02 18:50:49 +01:00
termi-official 388350029b Propagate first_and_last option into nonlinear iterative solvers. 2021-11-02 18:48:43 +01:00
termi-official 3c3d220460 Last MINRES iteration was sometimes missing. 2021-11-02 18:37:22 +01:00
termi-official 0279dd9787 Fix FGMRES one-by-off. 2021-11-02 18:34:29 +01:00
termi-official 46d800d05e Report correct norm in BiCGStab. 2021-11-02 18:29:20 +01:00
termi-official 4e1420f1e2 Fix GMRES final iteration report. 2021-11-02 18:29:02 +01:00
termi-official a9e93adb3f Deduplicate MINRES output. 2021-11-02 17:48:18 +01:00
termi-official 3bd47c26fd Make FGMRES output more consistent. 2021-11-02 17:21:01 +01:00
termi-official 3348ba2dbd Make BiCGStab output more consistent. 2021-11-02 17:07:46 +01:00
Ketan Mittal d8df06ab70 minor - put back accidentally removed lines 2021-11-01 10:33:07 -07:00
Ketan Mittal 5dcd85c7c7 make style 2021-11-01 09:18:20 -07:00
Ketan Mittal 047ca50acb Merge branch 'master' of https://github.com/mfem/mfem into tmop-renameandrefactor 2021-11-01 09:17:34 -07:00
Tzanio Kolev 4d1a0ac51d Update CHANGELOG 2021-10-31 14:59:55 -07:00
Tzanio c887336fa4 Updated CHANGELOG 2021-10-31 14:59:06 -07:00
Tzanio Kolev cde981b141 Merge branch 'master' into bicgstab-print-lvl-fix 2021-10-31 14:50:34 -07:00
Will Pazner 8d1780147d Add 2D compile-time option to performance/ex1[p]
Adds sample runs testing the 2D templated code to avoid regressions
2021-10-29 17:40:27 -07:00
Will Pazner f452882d29 Add unary + operator to AutoSIMD 2021-10-29 16:22:13 -07:00
Ketan Mittal 0d45eab46d initial commit with some refactoring and renaming 2021-10-20 10:38:36 -07:00
Will Pazner 15fbf80e57 Remove unneeded comment 2021-10-12 21:27:21 -07:00
Will Pazner e3450da53e Minor changes to MINRESSolver print level 2021-10-12 20:49:35 -07:00
Will Pazner 1c4b657f63 Minor changes to SLISolver printing 2021-10-12 20:37:12 -07:00
Will Pazner a19e625fc9 Delete empty comments 2021-10-12 18:10:40 -07:00
Will Pazner 37bd0ec75d Edits to IterativeSolver Doxygen comments 2021-10-12 18:03:59 -07:00
Will Pazner 020d2173c3 Minor changes to "FirstAndLast" print options 2021-10-12 11:39:54 -07:00
termi-official b889952837 Fix new macros in MSVC. 2021-09-29 16:54:26 +02:00
termi-official 8d96e40489 Make style... 2021-09-29 01:40:03 +02:00
termi-official d493efa295 Tweak docs. 2021-09-29 01:32:34 +02:00
termi-official a43de945bd Forgot a file. 2021-09-29 01:16:36 +02:00
termi-official 32a7939e4e Incorporate Will's changes. 2021-09-29 01:15:02 +02:00
termi-official 1f17768e33 Print warnings to mfem::err 2021-09-29 00:44:45 +02:00
termi-official 91f59a269c Revamp 'first and last' print option from legacy print level 3. 2021-09-29 00:41:59 +02:00
termi-official f8d501562f Exchange PrintLevel enum with struct+builder pattern. 2021-09-29 00:16:25 +02:00
termi-official b6530be3f9 Make style... 2021-09-28 21:45:51 +02:00
termi-official 3cbfee625a Add portable way to skip warnings generated by internal usage of deprecated variables. 2021-09-28 21:42:18 +02:00
Dennis Ogiermann 0ec9789849 Merge pull request #2572 from mfem/deprecate-print-level
Fix deprecation warnings for IterativeSolver::print_level
2021-09-28 18:25:50 +02:00
Will Pazner 7a20d17f90 Add pragma to ignore deprecated warnings in IterativeSolver::IterativeSolver 2021-09-27 15:47:04 -07:00
Will Pazner 317032fdbd Add bitwise or overload for PrintLevel 2021-09-27 15:34:20 -07:00
Will Pazner d6a0a5a13d Fix deprecation warnings for IterativeSolver::print_level 2021-09-27 15:21:07 -07:00
termi-official c1c39a51b5 Make style... 2021-09-27 23:40:41 +02:00
termi-official c5a7b179d8 Refactor setter/getter logic for print levels. 2021-09-27 23:36:47 +02:00
termi-official e88b99c2e2 Incorporate Will's suggestion: struct PrintOptions -> enum PrintLevel + deprecate old mechanism 2021-09-27 23:24:39 +02:00
Tzanio 793222a7e1 Merge branch 'master' into bicgstab-print-lvl-fix
Conflicts:
	CHANGELOG
2021-09-27 09:13:41 -07:00
termi-official 05bf01378c Add to changelog. 2021-09-27 18:11:53 +02:00
termi-official 71b7497037 Fix math string. 2021-09-27 17:30:49 +02:00
termi-official 5279744c6e Fix out of sync with upstream. 2021-09-27 17:04:05 +02:00
termi-official 19fae213a5 Fix some doxygen issues where no doc has been generated. 2021-09-27 16:53:42 +02:00
termi-official 0768f1940f Specify convergence docs. 2021-09-25 19:32:29 +02:00
termi-official 0c9296ab4b Propagate print_options through all solvers and unify output. 2021-09-25 19:13:40 +02:00
termi-official 7632309ed6 Fix hidden assumption taken when testing (all ranks print out if MPI is not properly initialized). 2021-09-25 17:27:54 +02:00
termi-official 439a3d773c Give recommendation and reasoning for new printing behavior. 2021-09-25 14:38:00 +02:00
termi-official ec2e958e04 Introduce PrintOptions as a new mechanism to set the reporting policy of iterative solvers + some Docs. 2021-09-25 14:28:25 +02:00
termi-official 53eb418526 Merge branch 'master' into bicgstab-print-lvl-fix 2021-09-08 15:10:26 +02:00
termi-official eea189ba98 Change BICGStabSolver print level semantics. 2021-09-08 14:55:03 +02:00
138 changed files with 2628 additions and 7050 deletions
+1 -1
View File
@@ -45,5 +45,5 @@ variables:
- echo ${MFEM_DATA_DIR}
- echo ${SPEC}
# Next script uses 'THREADS': leaving it empty --> it uses 'make all -j'
- lalloc 1 -W 30 -q pdebug tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
- lalloc 1 -W 30 -q pdebug --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
needs: [setup]
+13 -5
View File
@@ -10,10 +10,16 @@
Version 4.3.1 (development)
===========================
- Added support for automatic differentiation. Users can select between
native implementation and external library implementation at the
configuration phase. A parallel and two serial examples are implemented
in the autodiff miniapp directory.
- More explicit and consistent formating of the output of iterative solvers
with the new IterativeSolver::PrintLevel options. See linalg/solvers.hpp.
- Added support for automatic differentiation. Users can select between native
implementation and external library implementation during configuration. One
parallel and two serial examples are implemented in the miniapps/autodiff/
directory.
- GridFunctionCoefficient (and the related vector, gradient, divergence, and
curl classes) now work properly with LORDiscretization and LORSolver.
- Added support for mesh preprocessing to resolve fine scale problem data
before simulation. This feature uses adaptive mesh refinement to control the
@@ -71,13 +77,15 @@ Version 4.3.1 (development)
- Added initial TMOP-based capabilities for surface fitting and tangential
relaxation in the mesh-optimizer and pmesh-optimizer miniapps.
- Added ParMesh Adjaceny Set (adjset) creation support to the Conduit Mesh
Blueprint MFEM wrapper functions in ConduitDataCollection.
- `HypreParVector` and `Vector` now support move semantics, and the copy
constructor for `HypreParVector` now copies the local vector data.
- The HPC versions of ex1 and ex1p (in miniapps/performance) now support
runtime selection of either 2D or 3D meshes.
Version 4.3, released on July 29, 2021
======================================
+133
View File
@@ -0,0 +1,133 @@
# MFEM Code of Conduct
## Our Pledge
We as members, contributors, and leaders pledge to make participation in our
community a harassment-free experience for everyone, regardless of age, body
size, visible or invisible disability, ethnicity, sex characteristics, gender
identity and expression, level of experience, education, socio-economic status,
nationality, personal appearance, race, caste, color, religion, or sexual
identity and orientation.
We pledge to act and interact in ways that contribute to an open, welcoming,
diverse, inclusive, and healthy community.
## Our Standards
Examples of behavior that contributes to a positive environment for our
community include:
* Demonstrating empathy and kindness toward other people
* Being respectful of differing opinions, viewpoints, and experiences
* Giving and gracefully accepting constructive feedback
* Accepting responsibility and apologizing to those affected by our mistakes,
and learning from the experience
* Focusing on what is best not just for us as individuals, but for the overall
community
Examples of unacceptable behavior include:
* The use of sexualized language or imagery, and sexual attention or advances of
any kind
* Trolling, insulting or derogatory comments, and personal or political attacks
* Public or private harassment
* Publishing others' private information, such as a physical or email address,
without their explicit permission
* Other conduct which could reasonably be considered inappropriate in a
professional setting
## Enforcement Responsibilities
Community leaders are responsible for clarifying and enforcing our standards of
acceptable behavior and will take appropriate and fair corrective action in
response to any behavior that they deem inappropriate, threatening, offensive,
or harmful.
Community leaders have the right and responsibility to remove, edit, or reject
comments, commits, code, wiki edits, issues, and other contributions that are
not aligned to this Code of Conduct, and will communicate reasons for moderation
decisions when appropriate.
## Scope
This Code of Conduct applies within all community spaces, and also applies when
an individual is officially representing the community in public spaces.
Examples of representing our community include using an official e-mail address,
posting via an official social media account, or acting as an appointed
representative at an online or offline event.
## Enforcement
Instances of abusive, harassing, or otherwise unacceptable behavior may be
reported to the community leaders responsible for enforcement at mfem@llnl.gov.
All complaints will be reviewed and investigated promptly and fairly.
All community leaders are obligated to respect the privacy and security of the
reporter of any incident. Anyone involved in the reported behavior will recuse
themselves from the investigation and decision making about the resolution of
the complaint.
## Enforcement Guidelines
Community leaders will follow these Community Impact Guidelines in determining
the consequences for any action they deem in violation of this Code of Conduct:
### 1. Correction
**Community Impact**: Use of inappropriate language or other behavior deemed
unprofessional or unwelcome in the community.
**Consequence**: A private, written warning from community leaders, providing
clarity around the nature of the violation and an explanation of why the
behavior was inappropriate. A public apology may be requested.
### 2. Warning
**Community Impact**: A violation through a single incident or series of
actions.
**Consequence**: A warning with consequences for continued behavior. No
interaction with the people involved, including unsolicited interaction with
those enforcing the Code of Conduct, for a specified period of time. This
includes avoiding interactions in community spaces as well as external channels
like social media. Violating these terms may lead to a temporary or permanent
ban.
### 3. Temporary Ban
**Community Impact**: A serious violation of community standards, including
sustained inappropriate behavior.
**Consequence**: A temporary ban from any sort of interaction or public
communication with the community for a specified period of time. No public or
private interaction with the people involved, including unsolicited interaction
with those enforcing the Code of Conduct, is allowed during this period.
Violating these terms may lead to a permanent ban.
### 4. Permanent Ban
**Community Impact**: Demonstrating a pattern of violation of community
standards, including sustained inappropriate behavior, harassment of an
individual, or aggression toward or disparagement of classes of individuals.
**Consequence**: A permanent ban from any sort of public interaction within the
community.
## Attribution
This Code of Conduct is adapted from the [Contributor Covenant][homepage],
version 2.1, available at
[https://www.contributor-covenant.org/version/2/1/code_of_conduct.html][v2.1].
Community Impact Guidelines were inspired by
[Mozilla's code of conduct enforcement ladder][Mozilla CoC].
For answers to common questions about this code of conduct, see the FAQ at
[https://www.contributor-covenant.org/faq][FAQ]. Translations are available at
[https://www.contributor-covenant.org/translations][translations].
[homepage]: https://www.contributor-covenant.org
[v2.1]: https://www.contributor-covenant.org/version/2/1/code_of_conduct.html
[Mozilla CoC]: https://github.com/mozilla/diversity
[FAQ]: https://www.contributor-covenant.org/faq
[translations]: https://www.contributor-covenant.org/translations
+3
View File
@@ -21,6 +21,9 @@ documentation; new examples and miniapps; HPC performance improvements; etc.
MFEM is distributed under the terms of the BSD-3 license. All new contributions
must be made under this license.
Note also that MFEM has a [Code of Conduct](CODE_OF_CONDUCT.md). By participating
in the MFEM community, you agree to abide by its rules.
If you plan on contributing to MFEM, consider reviewing the
[issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
already exists for your desired feature or the bug you ran into. Use a pull
+3 -3
View File
@@ -12,6 +12,9 @@ to enable high-performance scalable finite element discretization research and
application development on a wide variety of platforms, ranging from laptops to
supercomputers.
We welcome contributions and feedback from the community. Please see the file
CONTRIBUTING.md for additional details about our development process.
* For building instructions, see the file INSTALL, or type "make help".
* Copyright and licensing information can be found in files LICENSE and NOTICE.
@@ -19,9 +22,6 @@ supercomputers.
* The best starting point for new users interested in MFEM's features is to
review the examples and miniapps at https://mfem.org/examples.
* Developers interested in contributing to the library, should read the
instructions and documentation in the CONTRIBUTING.md file.
Conceptually, MFEM can be viewed as a finite element toolbox that provides the
building blocks for developing finite element algorithms in a manner similar to
that of MATLAB for linear algebra methods. In particular, MFEM provides support
+13 -12
View File
@@ -149,7 +149,7 @@ void InitialDeformation(const Vector &x, Vector &y);
void InitialVelocity(const Vector &x, Vector &v);
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes,
GridFunction *field, const char *field_name = NULL,
bool init_vis = false);
@@ -376,10 +376,10 @@ int main(int argc, char *argv[])
}
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes,
GridFunction *field, const char *field_name, bool init_vis)
{
if (!out)
if (!os)
{
return;
}
@@ -389,24 +389,25 @@ void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
mesh->SwapNodes(nodes, owns_nodes);
out << "solution\n" << *mesh << *field;
os << "solution\n" << *mesh << *field;
mesh->SwapNodes(nodes, owns_nodes);
if (init_vis)
{
out << "window_size 800 800\n";
out << "window_title '" << field_name << "'\n";
os << "window_size 800 800\n";
os << "window_title '" << field_name << "'\n";
if (mesh->SpaceDimension() == 2)
{
out << "view 0 0\n"; // view from top
out << "keys jl\n"; // turn off perspective and light
os << "view 0 0\n"; // view from top
os << "keys jl\n"; // turn off perspective and light
}
out << "keys cm\n"; // show colorbar and mesh
out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
out << "pause\n";
os << "keys cm\n"; // show colorbar and mesh
// update value-range; keep mesh-extents fixed
os << "autoscale value\n";
os << "pause\n";
}
out << flush;
os << flush;
}
+17 -13
View File
@@ -154,7 +154,8 @@ void InitialDeformation(const Vector &x, Vector &y);
void InitialVelocity(const Vector &x, Vector &v);
void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
void visualize(ostream &os, ParMesh *mesh,
ParGridFunction *deformed_nodes,
ParGridFunction *field, const char *field_name = NULL,
bool init_vis = false);
@@ -438,10 +439,11 @@ int main(int argc, char *argv[])
return 0;
}
void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
void visualize(ostream &os, ParMesh *mesh,
ParGridFunction *deformed_nodes,
ParGridFunction *field, const char *field_name, bool init_vis)
{
if (!out)
if (!os)
{
return;
}
@@ -451,25 +453,27 @@ void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
mesh->SwapNodes(nodes, owns_nodes);
out << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() << "\n";
out << "solution\n" << *mesh << *field;
os << "parallel " << mesh->GetNRanks()
<< " " << mesh->GetMyRank() << "\n";
os << "solution\n" << *mesh << *field;
mesh->SwapNodes(nodes, owns_nodes);
if (init_vis)
{
out << "window_size 800 800\n";
out << "window_title '" << field_name << "'\n";
os << "window_size 800 800\n";
os << "window_title '" << field_name << "'\n";
if (mesh->SpaceDimension() == 2)
{
out << "view 0 0\n"; // view from top
out << "keys jl\n"; // turn off perspective and light
os << "view 0 0\n"; // view from top
os << "keys jl\n"; // turn off perspective and light
}
out << "keys cm\n"; // show colorbar and mesh
out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
out << "pause\n";
os << "keys cm\n"; // show colorbar and mesh
// update value-range; keep mesh-extents fixed
os << "autoscale value\n";
os << "pause\n";
}
out << flush;
os << flush;
}
+10 -10
View File
@@ -32,7 +32,7 @@ private:
mutable DenseTensor flux;
mutable Vector z;
void GetFlux(const DenseMatrix &state, DenseTensor &flux) const;
void GetFlux(const DenseMatrix &state_, DenseTensor &flux_) const;
public:
FE_Evolution(FiniteElementSpace &vfes_,
@@ -256,26 +256,26 @@ inline double ComputeMaxCharSpeed(const Vector &state, const int dim)
}
// Compute the flux at solution nodes.
void FE_Evolution::GetFlux(const DenseMatrix &x, DenseTensor &flux) const
void FE_Evolution::GetFlux(const DenseMatrix &x_, DenseTensor &flux_) const
{
const int dof = flux.SizeI();
const int dim = flux.SizeJ();
const int flux_dof = flux_.SizeI();
const int flux_dim = flux_.SizeJ();
for (int i = 0; i < dof; i++)
for (int i = 0; i < flux_dof; i++)
{
for (int k = 0; k < num_equation; k++) { state(k) = x(i, k); }
ComputeFlux(state, dim, f);
for (int k = 0; k < num_equation; k++) { state(k) = x_(i, k); }
ComputeFlux(state, flux_dim, f);
for (int d = 0; d < dim; d++)
for (int d = 0; d < flux_dim; d++)
{
for (int k = 0; k < num_equation; k++)
{
flux(i, d, k) = f(k, d);
flux_(i, d, k) = f(k, d);
}
}
// Update max char speed
const double mcs = ComputeMaxCharSpeed(state, dim);
const double mcs = ComputeMaxCharSpeed(state, flux_dim);
if (mcs > max_char_speed) { max_char_speed = mcs; }
}
}
+13 -11
View File
@@ -171,7 +171,7 @@ public:
};
// Visualization driver
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes,
GridFunction *field, const char *field_name = NULL,
bool init_vis = false);
@@ -542,10 +542,10 @@ RubberOperator::~RubberOperator()
// Inline visualization
void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
void visualize(ostream &os, Mesh *mesh, GridFunction *deformed_nodes,
GridFunction *field, const char *field_name, bool init_vis)
{
if (!out)
if (!os)
{
return;
}
@@ -555,23 +555,25 @@ void visualize(ostream &out, Mesh *mesh, GridFunction *deformed_nodes,
mesh->SwapNodes(nodes, owns_nodes);
out << "solution\n" << *mesh << *field;
os << "solution\n" << *mesh << *field;
mesh->SwapNodes(nodes, owns_nodes);
if (init_vis)
{
out << "window_size 800 800\n";
out << "window_title '" << field_name << "'\n";
os << "window_size 800 800\n";
os << "window_title '" << field_name << "'\n";
if (mesh->SpaceDimension() == 2)
{
out << "view 0 0\n"; // view from top
out << "keys jlA\n"; // turn off perspective and light, +anti-aliasing
os << "view 0 0\n"; // view from top
// turn off perspective and light, +anti-aliasing
os << "keys jlA\n";
}
out << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
os << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
// update value-range; keep mesh-extents fixed
os << "autoscale value\n";
}
out << flush;
os << flush;
}
void ReferenceConfiguration(const Vector &x, Vector &y)
+17 -12
View File
@@ -185,7 +185,8 @@ public:
};
// Visualization driver
void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
void visualize(ostream &os, ParMesh *mesh,
ParGridFunction *deformed_nodes,
ParGridFunction *field, const char *field_name = NULL,
bool init_vis = false);
@@ -617,10 +618,11 @@ RubberOperator::~RubberOperator()
// Inline visualization
void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
void visualize(ostream &os, ParMesh *mesh,
ParGridFunction *deformed_nodes,
ParGridFunction *field, const char *field_name, bool init_vis)
{
if (!out)
if (!os)
{
return;
}
@@ -630,24 +632,27 @@ void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
mesh->SwapNodes(nodes, owns_nodes);
out << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() << "\n";
out << "solution\n" << *mesh << *field;
os << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() <<
"\n";
os << "solution\n" << *mesh << *field;
mesh->SwapNodes(nodes, owns_nodes);
if (init_vis)
{
out << "window_size 800 800\n";
out << "window_title '" << field_name << "'\n";
os << "window_size 800 800\n";
os << "window_title '" << field_name << "'\n";
if (mesh->SpaceDimension() == 2)
{
out << "view 0 0\n"; // view from top
out << "keys jlA\n"; // turn off perspective and light, +anti-aliasing
os << "view 0 0\n"; // view from top
// turn off perspective and light, +anti-aliasing
os << "keys jlA\n";
}
out << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
os << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
// update value-range; keep mesh-extents fixed
os << "autoscale value\n";
}
out << flush;
os << flush;
}
void ReferenceConfiguration(const Vector &x, Vector &y)
-266
View File
@@ -1,266 +0,0 @@
// MFEM Example 1
//
// Compile with: make ex1
//
// Sample runs: ex1 -m ../data/square-disc.mesh
// ex1 -m ../data/star.mesh
// ex1 -m ../data/escher.mesh
// ex1 -m ../data/fichera.mesh
// ex1 -m ../data/square-disc-p2.vtk -o 2
// ex1 -m ../data/square-disc-p3.mesh -o 3
// ex1 -m ../data/square-disc-nurbs.mesh -o -1
// ex1 -m ../data/disc-nurbs.mesh -o -1
// ex1 -m ../data/pipe-nurbs.mesh -o -1
// ex1 -m ../data/star-surf.mesh
// ex1 -m ../data/square-disc-surf.mesh
// ex1 -m ../data/inline-segment.mesh
// ex1 -m ../data/amr-quad.mesh
// ex1 -m ../data/amr-hex.mesh
// ex1 -m ../data/fichera-amr.mesh
// ex1 -m ../data/mobius-strip.mesh
// ex1 -m ../data/mobius-strip.mesh -o -1 -sc
//
// Description: This example code demonstrates the use of MFEM to define a
// simple finite element discretization of the Laplace problem
// -Delta u = 1 with homogeneous Dirichlet boundary conditions.
// Specifically, we discretize using a FE space of the specified
// order, or if order < 1 using an isoparametric/isogeometric
// space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
// NURBS mesh, etc.)
//
// The example highlights the use of mesh refinement, finite
// element grid functions, as well as linear and bilinear forms
// corresponding to the left-hand side and right-hand side of the
// discrete linear system. We also cover the explicit elimination
// of essential boundary conditions, static condensation, and the
// optional connection to the GLVis tool for visualization.
#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 = "./star-set.mesh";
int order = 1;
int rs = -1;
int ra = 0;
int bt = EntitySets::INVALID;
const char *bs = "Origin";
bool static_cond = false;
bool visualization = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree) or -1 for"
" isoparametric space.");
args.AddOption(&rs, "-rs", "--refine-serial",
"Number of serial refinement levels");
args.AddOption(&ra, "-ra", "--refine-adaptive",
"Number of adaptive refinement levels");
args.AddOption(&bt, "-bt", "--bc-entity-type",
"");
args.AddOption(&bs, "-bs", "--bc-entity-set-name",
"");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
args.PrintUsage(cout);
return 1;
}
args.PrintOptions(cout);
// 2. Read the mesh from the given mesh file. We can handle triangular,
// quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
// the same code.
Mesh *mesh = new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();
// 3. Refine the mesh to increase the resolution. In this example we do
// 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
// largest number that gives a final mesh with no more than 50,000
// elements.
{
int ref_levels = ( rs >= 0 ) ? rs :
(int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
for (int l = 0; l < ref_levels; l++)
{
mesh->UniformRefinement();
}
}
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
if ( ra > 0 )
{
cout << "calling EnsureNCMesh" << endl;
mesh->EnsureNCMesh();
cout << "back from EnsureNCMesh" << endl;
}
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
cout << "Calling RandomRefinement " << ra << " times." << endl;
for (int l = 0; l < ra; l++)
{
mesh->RandomRefinement(0.2);
}
cout << "Done with refinement" << endl;
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
if ( mesh->ncmesh )
{
mesh->ncmesh->PrintStats(cout);
ofstream ofsV("vp.out");
ofstream ofsE("ce.out");
mesh->ncmesh->PrintVertexParents(ofsV);
mesh->ncmesh->PrintCoarseElements(ofsE);
}
// 4. Define a finite element space on the mesh. Here we use continuous
// Lagrange finite elements of the specified order. If order < 1, we
// instead use an isoparametric/isogeometric space.
FiniteElementCollection *fec;
if (order > 0)
{
fec = new H1_FECollection(order, dim);
}
else if (mesh->GetNodes())
{
fec = mesh->GetNodes()->OwnFEC();
cout << "Using isoparametric FEs: " << fec->Name() << endl;
}
else
{
fec = new H1_FECollection(order = 1, dim);
}
FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
cout << "Number of finite element unknowns: "
<< fespace->GetTrueVSize() << endl;
// 5. Determine the list of true (i.e. conforming) essential boundary dofs.
// In this example, the boundary conditions are defined by marking all
// the boundary attributes from the mesh as essential (Dirichlet) and
// converting them to a list of true dofs.
Array<int> ess_tdof_list;
if ( bt == EntitySets::INVALID )
{
if (mesh->bdr_attributes.Size())
{
Array<int> ess_bdr(mesh->bdr_attributes.Max());
ess_bdr = 1;
fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
}
else
{
fespace->GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
ess_tdof_list);
}
cout << "Number of Dirichlet dofs: " << ess_tdof_list.Size() << endl;
// 6. Set up the linear form b(.) which corresponds to the right-hand side of
// the FEM linear system, which in this case is (1,phi_i) where phi_i are
// the basis functions in the finite element fespace.
LinearForm *b = new LinearForm(fespace);
ConstantCoefficient one(1.0);
b->AddDomainIntegrator(new DomainLFIntegrator(one));
b->Assemble();
// 7. Define the solution vector x as a finite element grid function
// corresponding to fespace. Initialize x with initial guess of zero,
// which satisfies the boundary conditions.
GridFunction x(fespace);
x = 0.0;
// 8. 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 = new BilinearForm(fespace);
a->AddDomainIntegrator(new DiffusionIntegrator(one));
// 9. Assemble the bilinear form and the corresponding linear system,
// applying any necessary transformations such as: eliminating boundary
// conditions, applying conforming constraints for non-conforming AMR,
// static condensation, etc.
if (static_cond) { a->EnableStaticCondensation(); }
a->Assemble();
SparseMatrix A;
Vector B, X;
a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
cout << "Size of linear system: " << A.Height() << endl;
#ifndef MFEM_USE_SUITESPARSE
// 10. Define a simple symmetric Gauss-Seidel preconditioner and use it to
// solve the system A X = B with PCG.
GSSmoother M(A);
PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
#else
// 10. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
UMFPackSolver umf_solver;
umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
umf_solver.SetOperator(A);
umf_solver.Mult(B, X);
#endif
// 11. Recover the solution as a finite element grid function.
a->RecoverFEMSolution(X, *b, x);
// 12. Save the refined mesh and the solution. This output can be viewed later
// using GLVis: "glvis -m refined.mesh -g sol.gf".
ofstream mesh_ofs("refined.mesh");
mesh_ofs.precision(8);
mesh->Print(mesh_ofs);
ofstream sol_ofs("sol.gf");
sol_ofs.precision(8);
x.Save(sol_ofs);
// 13. Send the solution by socket to a GLVis server.
if (visualization)
{
char vishost[] = "localhost";
int visport = 19916;
socketstream sol_sock(vishost, visport);
sol_sock.precision(8);
sol_sock << "solution\n" << *mesh << x << flush;
}
// 14. Free the used memory.
delete a;
delete b;
delete fespace;
if (order > 0) { delete fec; }
delete mesh;
return 0;
}
-388
View File
@@ -1,388 +0,0 @@
// MFEM Example 1 - Parallel Version
//
// Compile with: make ex1p
//
// Sample runs: mpirun -np 4 ex1p -m ../data/square-disc.mesh
// mpirun -np 4 ex1p -m ../data/star.mesh
// mpirun -np 4 ex1p -m ../data/escher.mesh
// mpirun -np 4 ex1p -m ../data/fichera.mesh
// mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
// mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
// mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
// mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
// mpirun -np 4 ex1p -m ../data/star-surf.mesh
// mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
// mpirun -np 4 ex1p -m ../data/inline-segment.mesh
// mpirun -np 4 ex1p -m ../data/amr-quad.mesh
// mpirun -np 4 ex1p -m ../data/amr-hex.mesh
// mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
// mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
//
// The following are examples of using EntitySets to define
// homogeneous Dirichlet boundary condition. These examples
// require a modified mesh file and a specialized version of
// example 1 called "ex1p_es".
// mpirun -np 4 ex1p_es -m ./fichera-set.mesh -bt 0 -bs Origin
// mpirun -np 4 ex1p_es -m ./fichera-set.mesh -bt 1 -bs Axes
// mpirun -np 4 ex1p_es -m ./fichera-set.mesh
// -bt 1 -bs "Negative Axes"
// mpirun -np 4 ex1p_es -m ./fichera-set.mesh
// -bt 2 -bs "Interior Corner"
// mpirun -np 4 ex1p_es -m ./fichera-set.mesh
// -bt 2 -bs "Exterior Corner"
// mpirun -np 4 ex1p_es -m ./fichera-set.mesh
// -bt 3 -bs "Interior Corner"
// mpirun -np 4 ex1p_es -m ./fichera-set.mesh
// -bt 3 -bs "Exterior Corner"
// mpirun -np 4 ex1p_es -m ./fichera-set.mesh -bt 3 -bs "Steps"
//
// Description: This example code demonstrates the use of MFEM to define a
// simple finite element discretization of the Laplace problem
// -Delta u = 1 with homogeneous Dirichlet boundary conditions.
// Specifically, we discretize using a FE space of the specified
// order, or if order < 1 using an isoparametric/isogeometric
// space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
// NURBS mesh, etc.)
//
// The example highlights the use of mesh refinement, finite
// element grid functions, as well as linear and bilinear forms
// corresponding to the left-hand side and right-hand side of the
// discrete linear system. We also cover the explicit elimination
// of essential boundary conditions, static condensation, and the
// optional connection to the GLVis tool for visualization.
#include "mfem.hpp"
#include <fstream>
#include <iostream>
using namespace std;
using namespace mfem;
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
int num_procs, myid;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
// 2. Parse command-line options.
const char *mesh_file = "./star-set.mesh";
int order = 1;
int rs = -1;
int rp = 2;
int ra = 0;
int bt = EntitySets::INVALID;
const char *bs = "Origin";
bool static_cond = false;
bool visualization = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree) or -1 for"
" isoparametric space.");
args.AddOption(&rs, "-rs", "--refine-serial",
"Number of serial refinement levels");
args.AddOption(&rp, "-rp", "--refine-parallel",
"Number of parallel refinement levels");
args.AddOption(&ra, "-ra", "--refine-adaptive",
"Number of adaptive refinement levels");
args.AddOption(&bt, "-bt", "--bc-entity-type",
"");
args.AddOption(&bs, "-bs", "--bc-entity-set-name",
"");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
// 3. Read the (serial) mesh from the given mesh file on all processors. We
// can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
// and volume meshes with the same code.
Mesh *mesh = new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();
// 4. Refine the serial mesh on all processors to increase the resolution. In
// this example we do 'ref_levels' of uniform refinement. We choose
// 'ref_levels' to be the largest number that gives a final mesh with no
// more than 10,000 elements.
{
int ref_levels = ( rs >= 0 ) ? rs :
(int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
for (int l = 0; l < ref_levels; l++)
{
if ( myid == 0 ) { cout << "Uniform refinement in serial..."; }
mesh->UniformRefinement();
}
MPI_Barrier(MPI_COMM_WORLD);
if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
}
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
/*
At this point we have a serial mesh containing an EntitySets
object which stores the current node/edge/face/element indices
for each entity in each set. This data is duplicated on each MPI
rank.
*/
if ( ra > 0 )
{
cout << "calling EnsureNCMesh" << endl;
mesh->EnsureNCMesh();
cout << "back from EnsureNCMesh" << endl;
}
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
/*
We now have an NCEntitySets object which stores the node indices
describing each enity in each node/edge/face set and the element
indices for the elements in each element set. This data is
duplicated on each MPI rank.
*/
// 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
// this mesh further in parallel to increase the resolution. Once the
// parallel mesh is defined, the serial mesh can be deleted.
cout << "creating ParMesh from serial mesh" << endl;
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
cout << "done creating ParMesh from serial mesh" << endl;
delete mesh;
if ( pmesh->pent_sets )
{
cout << "pmesh->pent_sets is non NULL" << endl;
pmesh->pent_sets->PrintSetInfo(cout);
}
else
{
cout << "pmesh->pent_sets is NULL" << endl;
}
/*
We now have a ParEntitySets object which marshals the data stored
in EntitySets objects. The data has now been pruned so that each
rank only contains indices of local entities.
The NCEntitySets object remains unchanged...
If we have an NC mesh a different path is taken and the
EntitySets are ignored.
1) ParNCMesh is created from NCMesh
a) Creates a ParNCEntitySets object from ncmesh (every rank contains
information to find every entity)
2) ParNCMesh is pruned which involves renumbering elements and vertices
3) ParMesh is initialized from ParNCMesh
4) ParNCMesh::OnMeshUpdated is called
5) Mesh::GenerateNCFaceInfo is called
*/
{
int par_ref_levels = rp;
for (int l = 0; l < par_ref_levels; l++)
{
if ( myid == 0 ) { cout << "Uniform refinement in parallel..."; }
pmesh->UniformRefinement();
}
MPI_Barrier(MPI_COMM_WORLD);
if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
}
/*
RandomRefinement will end up calling
ParMesh::NonconformingRefinement which will create a new ParMesh
object using the ParNCMesh object and then call
ParMesh::OnMeshUpdated on this new mesh.
*/
for (int l = 0; l < ra; l++)
{
pmesh->RandomRefinement(0.2);
}
if ( ra > 0 )
{
if ( pmesh->pent_sets )
{
cout << "pmesh->pent_sets is non NULL post random refinement" << endl;
pmesh->pent_sets->PrintSetInfo(cout);
}
else
{
cout << "pmesh->pent_sets is NULL post random refinement" << endl;
}
}
// 6. Define a parallel finite element space on the parallel mesh. Here we
// use continuous Lagrange finite elements of the specified order. If
// order < 1, we instead use an isoparametric/isogeometric space.
FiniteElementCollection *fec;
if (order > 0)
{
fec = new H1_FECollection(order, dim);
}
else if (pmesh->GetNodes())
{
fec = pmesh->GetNodes()->OwnFEC();
if (myid == 0)
{
cout << "Using isoparametric FEs: " << fec->Name() << endl;
}
}
else
{
fec = new H1_FECollection(order = 1, dim);
}
ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
HYPRE_Int size = fespace->GlobalTrueVSize();
if (myid == 0)
{
cout << "Number of finite element unknowns: " << size << endl;
}
// 7. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
if ( bt == EntitySets::INVALID )
{
if (pmesh->bdr_attributes.Size())
{
Array<int> ess_bdr(pmesh->bdr_attributes.Max());
ess_bdr = 1;
fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
}
else
{
fespace->GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
ess_tdof_list);
}
for (int i=0; i<num_procs; i++)
{
if (myid == i)
{
cout << "Number of Dirichlet dofs on proc " << i << ": "
<< ess_tdof_list.Size() << endl;
}
}
// 8. Set up the parallel linear form b(.) which corresponds to the
// right-hand side of the FEM linear system, which in this case is
// (1,phi_i) where phi_i are the basis functions in fespace.
ParLinearForm *b = new ParLinearForm(fespace);
ConstantCoefficient one(1.0);
b->AddDomainIntegrator(new DomainLFIntegrator(one));
b->Assemble();
// 9. Define the solution vector x as a parallel finite element grid function
// corresponding to fespace. Initialize x with initial guess of zero,
// which satisfies the boundary conditions.
ParGridFunction x(fespace);
x = 0.0;
// 10. Set up the parallel bilinear form a(.,.) on the finite element space
// corresponding to the Laplacian operator -Delta, by adding the Diffusion
// domain integrator.
ParBilinearForm *a = new ParBilinearForm(fespace);
a->AddDomainIntegrator(new DiffusionIntegrator(one));
// 11. Assemble the parallel bilinear form and the corresponding linear
// system, applying any necessary transformations such as: parallel
// assembly, eliminating boundary conditions, applying conforming
// constraints for non-conforming AMR, static condensation, etc.
if (static_cond) { a->EnableStaticCondensation(); }
a->Assemble();
HypreParMatrix A;
Vector B, X;
a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
if (myid == 0)
{
cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
}
// 12. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
// preconditioner from hypre.
HypreSolver *amg = new HypreBoomerAMG(A);
HyprePCG *pcg = new HyprePCG(A);
pcg->SetTol(1e-12);
pcg->SetMaxIter(200);
pcg->SetPrintLevel(2);
pcg->SetPreconditioner(*amg);
pcg->Mult(B, X);
// 13. Recover the parallel grid function corresponding to X. This is the
// local finite element solution on each processor.
a->RecoverFEMSolution(X, *b, x);
// 14. Save the refined mesh and the solution in parallel. This output can
// be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
{
ostringstream mesh_name, sol_name;
mesh_name << "mesh." << setfill('0') << setw(6) << myid;
sol_name << "sol." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh->Print(mesh_ofs);
ofstream sol_ofs(sol_name.str().c_str());
sol_ofs.precision(8);
x.Save(sol_ofs);
}
// 15. Send the solution by socket to a GLVis server.
if (visualization)
{
char vishost[] = "localhost";
int visport = 19916;
socketstream sol_sock(vishost, visport);
sol_sock << "parallel " << num_procs << " " << myid << "\n";
sol_sock.precision(8);
sol_sock << "solution\n" << *pmesh << x << flush;
}
// 16. Free the used memory.
delete pcg;
delete amg;
delete a;
delete b;
delete fespace;
if (order > 0) { delete fec; }
delete pmesh;
MPI_Finalize();
return 0;
}
+19 -18
View File
@@ -75,7 +75,7 @@ Mesh * GenerateSerialMesh(int ref);
// alpha*n.Grad(sol) + beta*sol - gamma over the same boundary.
double IntegrateBC(const GridFunction &sol, const Array<int> &bdr_marker,
double alpha, double beta, double gamma,
double &err);
double &error);
int main(int argc, char *argv[])
{
@@ -302,26 +302,26 @@ int main(int argc, char *argv[])
{
// Integrate the solution on the Dirichlet boundary and compare to the
// expected value.
double err, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, err);
double error, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, error);
bool hom_dbc = (dbc_val == 0.0);
err /= hom_dbc ? 1.0 : fabs(dbc_val);
error /= hom_dbc ? 1.0 : fabs(dbc_val);
mfem::out << "Average of solution on Gamma_dbc:\t"
<< avg << ", \t"
<< (hom_dbc ? "absolute" : "relative")
<< " error " << err << endl;
<< " error " << error << endl;
}
{
// Integrate n.Grad(u) on the inhomogeneous Neumann boundary and compare
// to the expected value.
double err, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, err);
double error, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, error);
bool hom_nbc = (nbc_val == 0.0);
err /= hom_nbc ? 1.0 : fabs(nbc_val);
error /= hom_nbc ? 1.0 : fabs(nbc_val);
mfem::out << "Average of n.Grad(u) on Gamma_nbc:\t"
<< avg << ", \t"
<< (hom_nbc ? "absolute" : "relative")
<< " error " << err << endl;
<< " error " << error << endl;
}
{
// Integrate n.Grad(u) on the homogeneous Neumann boundary and compare to
@@ -330,25 +330,26 @@ int main(int argc, char *argv[])
nbc0_bdr = 0;
nbc0_bdr[3] = 1;
double err, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, err);
double error, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, error);
bool hom_nbc = true;
mfem::out << "Average of n.Grad(u) on Gamma_nbc0:\t"
<< avg << ", \t"
<< (hom_nbc ? "absolute" : "relative")
<< " error " << err << endl;
<< " error " << error << endl;
}
{
// Integrate n.Grad(u) + a * u on the Robin boundary and compare to the
// expected value.
double err, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, err);
double error;
double avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, error);
bool hom_rbc = (rbc_b_val == 0.0);
err /= hom_rbc ? 1.0 : fabs(rbc_b_val);
error /= hom_rbc ? 1.0 : fabs(rbc_b_val);
mfem::out << "Average of n.Grad(u)+a*u on Gamma_rbc:\t"
<< avg << ", \t"
<< (hom_rbc ? "absolute" : "relative")
<< " error " << err << endl;
<< " error " << error << endl;
}
// 14. Save the refined mesh and the solution. This output can be viewed
@@ -637,11 +638,11 @@ Mesh * GenerateSerialMesh(int ref)
double IntegrateBC(const GridFunction &x, const Array<int> &bdr,
double alpha, double beta, double gamma,
double &err)
double &error)
{
double nrm = 0.0;
double avg = 0.0;
err = 0.0;
error = 0.0;
const bool a_is_zero = alpha == 0.0;
const bool b_is_zero = beta == 0.0;
@@ -705,20 +706,20 @@ double IntegrateBC(const GridFunction &x, const Array<int> &bdr,
// Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
val -= gamma;
err += (val*val) * ip.weight * face_weight;
error += (val*val) * ip.weight * face_weight;
}
}
// Normalize by the length of the boundary
if (std::abs(nrm) > 0.0)
{
err /= nrm;
error /= nrm;
avg /= nrm;
}
// Compute l2 norm of the error in the boundary condition (negative
// quadrature weights may produce negative 'err')
err = (err >= 0.0) ? sqrt(err) : -sqrt(-err);
// quadrature weights may produce negative 'error')
error = (error >= 0.0) ? sqrt(error) : -sqrt(-error);
// Return the average value of alpha * n.Grad(x) + beta * x
return avg;
+17 -16
View File
@@ -75,7 +75,7 @@ Mesh * GenerateSerialMesh(int ref);
// alpha*n.Grad(sol) + beta*sol - gamma over the same boundary.
double IntegrateBC(const ParGridFunction &sol, const Array<int> &bdr_marker,
double alpha, double beta, double gamma,
double &err);
double &error);
int main(int argc, char *argv[])
{
@@ -321,26 +321,26 @@ int main(int argc, char *argv[])
{
// Integrate the solution on the Dirichlet boundary and compare to the
// expected value.
double err, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, err);
double error, avg = IntegrateBC(u, dbc_bdr, 0.0, 1.0, dbc_val, error);
bool hom_dbc = (dbc_val == 0.0);
err /= hom_dbc ? 1.0 : fabs(dbc_val);
error /= hom_dbc ? 1.0 : fabs(dbc_val);
mfem::out << "Average of solution on Gamma_dbc:\t"
<< avg << ", \t"
<< (hom_dbc ? "absolute" : "relative")
<< " error " << err << endl;
<< " error " << error << endl;
}
{
// Integrate n.Grad(u) on the inhomogeneous Neumann boundary and compare
// to the expected value.
double err, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, err);
double error, avg = IntegrateBC(u, nbc_bdr, 1.0, 0.0, nbc_val, error);
bool hom_nbc = (nbc_val == 0.0);
err /= hom_nbc ? 1.0 : fabs(nbc_val);
error /= hom_nbc ? 1.0 : fabs(nbc_val);
mfem::out << "Average of n.Grad(u) on Gamma_nbc:\t"
<< avg << ", \t"
<< (hom_nbc ? "absolute" : "relative")
<< " error " << err << endl;
<< " error " << error << endl;
}
{
// Integrate n.Grad(u) on the homogeneous Neumann boundary and compare to
@@ -349,25 +349,26 @@ int main(int argc, char *argv[])
nbc0_bdr = 0;
nbc0_bdr[3] = 1;
double err, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, err);
double error, avg = IntegrateBC(u, nbc0_bdr, 1.0, 0.0, 0.0, error);
bool hom_nbc = true;
mfem::out << "Average of n.Grad(u) on Gamma_nbc0:\t"
<< avg << ", \t"
<< (hom_nbc ? "absolute" : "relative")
<< " error " << err << endl;
<< " error " << error << endl;
}
{
// Integrate n.Grad(u) + a * u on the Robin boundary and compare to the
// expected value.
double err, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val, err);
double error, avg = IntegrateBC(u, rbc_bdr, 1.0, rbc_a_val, rbc_b_val,
error);
bool hom_rbc = (rbc_b_val == 0.0);
err /= hom_rbc ? 1.0 : fabs(rbc_b_val);
error /= hom_rbc ? 1.0 : fabs(rbc_b_val);
mfem::out << "Average of n.Grad(u)+a*u on Gamma_rbc:\t"
<< avg << ", \t"
<< (hom_rbc ? "absolute" : "relative")
<< " error " << err << endl;
<< " error " << error << endl;
}
// 15. Save the refined mesh and the solution in parallel. This output can be
@@ -667,11 +668,11 @@ double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
double loc_vals[3];
double &nrm = loc_vals[0];
double &avg = loc_vals[1];
double &err = loc_vals[2];
double &error = loc_vals[2];
nrm = 0.0;
avg = 0.0;
err = 0.0;
error = 0.0;
const bool a_is_zero = alpha == 0.0;
const bool b_is_zero = beta == 0.0;
@@ -735,7 +736,7 @@ double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
// Integrate |alpha * n.Grad(x) + beta * x - gamma|^2
val -= gamma;
err += (val*val) * ip.weight * face_weight;
error += (val*val) * ip.weight * face_weight;
}
}
@@ -754,7 +755,7 @@ double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
}
// Compute l2 norm of the error in the boundary condition (negative
// quadrature weights may produce negative 'err')
// quadrature weights may produce negative 'error')
glb_err = (glb_err >= 0.0) ? sqrt(glb_err) : -sqrt(-glb_err);
// Return the average value of alpha * n.Grad(x) + beta * x
+2 -2
View File
@@ -167,9 +167,9 @@ int main(int argc, char *argv[])
// 13. Compute error in the solution and its flux
FunctionCoefficient uCoef(uExact);
double err = x.ComputeL2Error(uCoef);
double error = x.ComputeL2Error(uCoef);
cout << "|u - u_h|_2 = " << err << endl;
cout << "|u - u_h|_2 = " << error << endl;
FiniteElementSpace flux_fespace(mesh, &fec, 3);
GridFunction flux(&flux_fespace);
+2 -2
View File
@@ -197,9 +197,9 @@ int main(int argc, char *argv[])
// 15. Compute error in the solution and its flux
FunctionCoefficient uCoef(uExact);
double err = x.ComputeL2Error(uCoef);
double error = x.ComputeL2Error(uCoef);
if (myid == 0) { cout << "|u - u_h|_2 = " << err << endl; }
if (myid == 0) { cout << "|u - u_h|_2 = " << error << endl; }
ParFiniteElementSpace flux_fespace(&pmesh, &fec, 3);
ParGridFunction flux(&flux_fespace);
+2 -2
View File
@@ -255,10 +255,10 @@ int main(int argc, char *argv[])
// 15. Compute and print the L^2 norm of the error.
{
double err = x.ComputeL2Error(E);
double error = x.ComputeL2Error(E);
if (myid == 0)
{
cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
cout << "\n|| E_h - E ||_{L^2} = " << error << '\n' << endl;
}
}
-411
View File
@@ -1,411 +0,0 @@
// MFEM Example 3 - Parallel Version
//
// Compile with: make ex3p
//
// Sample runs: mpirun -np 4 ex3p -m ../data/star.mesh
// mpirun -np 4 ex3p -m ../data/square-disc.mesh -o 2
// mpirun -np 4 ex3p -m ../data/beam-tet.mesh
// mpirun -np 4 ex3p -m ../data/beam-hex.mesh
// mpirun -np 4 ex3p -m ../data/escher.mesh
// mpirun -np 4 ex3p -m ../data/fichera.mesh
// mpirun -np 4 ex3p -m ../data/fichera-q2.vtk
// mpirun -np 4 ex3p -m ../data/fichera-q3.mesh
// mpirun -np 4 ex3p -m ../data/square-disc-nurbs.mesh
// mpirun -np 4 ex3p -m ../data/beam-hex-nurbs.mesh
// mpirun -np 4 ex3p -m ../data/amr-quad.mesh -o 2
// mpirun -np 4 ex3p -m ../data/amr-hex.mesh
// mpirun -np 4 ex3p -m ../data/star-surf.mesh -o 2
// mpirun -np 4 ex3p -m ../data/mobius-strip.mesh -o 2 -f 0.1
// mpirun -np 4 ex3p -m ../data/klein-bottle.mesh -o 2 -f 0.1
//
// Description: This example code solves a simple electromagnetic diffusion
// problem corresponding to the second order definite Maxwell
// equation curl curl E + E = f with boundary condition
// E x n = <given tangential field>. Here, we use a given exact
// solution E and compute the corresponding r.h.s. f.
// We discretize with Nedelec finite elements in 2D or 3D.
//
// The example demonstrates the use of H(curl) finite element
// spaces with the curl-curl and the (vector finite element) mass
// bilinear form, as well as the computation of discretization
// error when the exact solution is known. Static condensation is
// also illustrated.
//
// We recommend viewing examples 1-2 before viewing this example.
#include "mfem.hpp"
#include <fstream>
#include <iostream>
using namespace std;
using namespace mfem;
// Exact solution, E, and r.h.s., f. See below for implementation.
//void E_exact(const Vector &, Vector &);
//void f_exact(const Vector &, Vector &);
//double freq = 1.0, kappa;
void f_const(const Vector &, Vector &);
int dim;
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
int num_procs, myid;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
// 2. Parse command-line options.
const char *mesh_file = "../data/beam-tet.mesh";
int order = 1;
int rs = -1;
int rp = 2;
int ra = 0;
int bt = EntitySets::INVALID;
const char *bs = "Origin";
bool static_cond = false;
bool visualization = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree).");
/*
args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
" solution.");
*/
args.AddOption(&rs, "-rs", "--refine-serial",
"Number of serial refinement levels");
args.AddOption(&rp, "-rp", "--refine-parallel",
"Number of parallel refinement levels");
args.AddOption(&ra, "-ra", "--refine-adaptive",
"Number of adaptive refinement levels");
args.AddOption(&bt, "-bt", "--bc-entity-type",
"");
args.AddOption(&bs, "-bs", "--bc-entity-set-name",
"");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
// kappa = freq * M_PI;
// 3. Read the (serial) mesh from the given mesh file on all processors. We
// can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
// and volume meshes with the same code.
Mesh *mesh = new Mesh(mesh_file, 1, 1);
dim = mesh->Dimension();
int sdim = mesh->SpaceDimension();
// 4. Refine the serial mesh on all processors to increase the resolution. In
// this example we do 'ref_levels' of uniform refinement. We choose
// 'ref_levels' to be the largest number that gives a final mesh with no
// more than 1,000 elements.
{
int ref_levels = ( rs >= 0 ) ? rs :
(int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
for (int l = 0; l < ref_levels; l++)
{
if ( myid == 0 ) { cout << "Uniform refinement in serial..."; }
mesh->UniformRefinement();
}
MPI_Barrier(MPI_COMM_WORLD);
if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
}
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
/*
At this point we have a serial mesh containing an EntitySets
object which stores the current node/edge/face/element indices
for each entity in each set. This data is duplicated on each MPI
rank.
*/
if ( ra > 0 )
{
cout << "calling EnsureNCMesh" << endl;
mesh->EnsureNCMesh();
cout << "back from EnsureNCMesh" << endl;
}
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
/*
We now have an NCEntitySets object which stores the node indices
describing each enity in each node/edge/face set and the element
indices for the elements in each element set. This data is
duplicated on each MPI rank.
*/
// 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
// this mesh further in parallel to increase the resolution. Once the
// parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
// meshes need to be reoriented before we can define high-order Nedelec
// spaces on them.
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
if ( pmesh->pent_sets )
{
cout << "pmesh->pent_sets is non NULL" << endl;
pmesh->pent_sets->PrintSetInfo(cout);
}
else
{
cout << "pmesh->pent_sets is NULL" << endl;
}
{
int par_ref_levels = rp;
for (int l = 0; l < par_ref_levels; l++)
{
if ( myid == 0 ) { cout << "Uniform refinement in parallel..."; }
pmesh->UniformRefinement();
}
MPI_Barrier(MPI_COMM_WORLD);
if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
}
pmesh->ReorientTetMesh();
pmesh->ent_sets->PrintSetInfo(cout);
for (int l = 0; l < ra; l++)
{
pmesh->RandomRefinement(0.2);
}
if ( ra > 0 )
{
if ( pmesh->pent_sets )
{
cout << "pmesh->pent_sets is non NULL post random refinement" << endl;
pmesh->pent_sets->PrintSetInfo(cout);
}
else
{
cout << "pmesh->pent_sets is NULL post random refinement" << endl;
}
}
// 6. Define a parallel finite element space on the parallel mesh. Here we
// use the Nedelec finite elements of the specified order.
FiniteElementCollection *fec = new ND_FECollection(order, dim);
ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
HYPRE_Int size = fespace->GlobalTrueVSize();
if (myid == 0)
{
cout << "Number of finite element unknowns: " << size << endl;
}
// 7. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
if ( bt == EntitySets::INVALID )
{
if (pmesh->bdr_attributes.Size())
{
Array<int> ess_bdr(pmesh->bdr_attributes.Max());
ess_bdr = 1;
fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
}
else
{
fespace->GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
ess_tdof_list);
}
if (myid == 0)
{
cout << "Number of Dirichlet dofs: " << ess_tdof_list.Size() << endl;
}
// 8. Set up the parallel linear form b(.) which corresponds to the
// right-hand side of the FEM linear system, which in this case is
// (f,phi_i) where f is given by the function f_exact and phi_i are the
// basis functions in the finite element fespace.
VectorFunctionCoefficient f(sdim, f_const);
ParLinearForm *b = new ParLinearForm(fespace);
b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
b->Assemble();
// 9. Define the solution vector x as a parallel finite element grid function
// corresponding to fespace. Initialize x by projecting the exact
// solution. Note that only values from the boundary edges will be used
// when eliminating the non-homogeneous boundary condition to modify the
// r.h.s. vector b.
ParGridFunction x(fespace);
// VectorFunctionCoefficient E(sdim, E_exact);
// x.ProjectCoefficient(E);
x = 0.0;
// 10. Set up the parallel bilinear form corresponding to the EM diffusion
// operator curl muinv curl + sigma I, by adding the curl-curl and the
// mass domain integrators.
Coefficient *muinv = new ConstantCoefficient(1.0);
Coefficient *sigma = new ConstantCoefficient(1.0);
ParBilinearForm *a = new ParBilinearForm(fespace);
a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
// 11. Assemble the parallel bilinear form and the corresponding linear
// system, applying any necessary transformations such as: parallel
// assembly, eliminating boundary conditions, applying conforming
// constraints for non-conforming AMR, static condensation, etc.
if (static_cond) { a->EnableStaticCondensation(); }
a->Assemble();
HypreParMatrix A;
Vector B, X;
a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
if (myid == 0)
{
cout << "Size of linear system: " << A.GetGlobalNumRows() << endl;
}
// 12. Define and apply a parallel PCG solver for AX=B with the AMS
// preconditioner from hypre.
ParFiniteElementSpace *prec_fespace =
(a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
HypreSolver *ams = new HypreAMS(A, prec_fespace);
HyprePCG *pcg = new HyprePCG(A);
pcg->SetTol(1e-12);
pcg->SetMaxIter(500);
pcg->SetPrintLevel(2);
pcg->SetPreconditioner(*ams);
pcg->Mult(B, X);
// 13. Recover the parallel grid function corresponding to X. This is the
// local finite element solution on each processor.
a->RecoverFEMSolution(X, *b, x);
/*
// 14. Compute and print the L^2 norm of the error.
{
double err = x.ComputeL2Error(E);
if (myid == 0)
{
cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
}
}
*/
// 15. Save the refined mesh and the solution in parallel. This output can
// be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
{
ostringstream mesh_name, sol_name;
mesh_name << "mesh." << setfill('0') << setw(6) << myid;
sol_name << "sol." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh->Print(mesh_ofs);
ofstream sol_ofs(sol_name.str().c_str());
sol_ofs.precision(8);
x.Save(sol_ofs);
}
// 16. Send the solution by socket to a GLVis server.
if (visualization)
{
char vishost[] = "localhost";
int visport = 19916;
socketstream sol_sock(vishost, visport);
sol_sock << "parallel " << num_procs << " " << myid << "\n";
sol_sock.precision(8);
sol_sock << "solution\n" << *pmesh << x << flush;
}
// 17. Free the used memory.
delete pcg;
delete ams;
delete a;
delete sigma;
delete muinv;
delete b;
delete fespace;
delete fec;
delete pmesh;
MPI_Finalize();
return 0;
}
/*
void E_exact(const Vector &x, Vector &E)
{
if (dim == 3)
{
E(0) = sin(kappa * x(1));
E(1) = sin(kappa * x(2));
E(2) = sin(kappa * x(0));
}
else
{
E(0) = sin(kappa * x(1));
E(1) = sin(kappa * x(0));
if (x.Size() == 3) { E(2) = 0.0; }
}
}
void f_exact(const Vector &x, Vector &f)
{
if (dim == 3)
{
f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
}
else
{
f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
if (x.Size() == 3) { f(2) = 0.0; }
}
}
*/
void f_const(const Vector &x, Vector &f)
{
if (dim == 3)
{
f(0) = 1.0;
f(1) = 1.0;
f(2) = 1.0;
}
else
{
f(0) = 1.0;
f(1) = 1.0;
if (x.Size() == 3) { f(2) = 0.0; }
}
}
+2 -2
View File
@@ -256,10 +256,10 @@ int main(int argc, char *argv[])
// 15. Compute and print the L^2 norm of the error.
{
double err = x.ComputeL2Error(F);
double error = x.ComputeL2Error(F);
if (myid == 0)
{
cout << "\n|| F_h - F ||_{L^2} = " << err << '\n' << endl;
cout << "\n|| F_h - F ||_{L^2} = " << error << '\n' << endl;
}
}
-438
View File
@@ -1,438 +0,0 @@
// MFEM Example 4 - Parallel Version
//
// Compile with: make ex4p
//
// Sample runs: mpirun -np 4 ex4p -m ../data/square-disc.mesh
// mpirun -np 4 ex4p -m ../data/star.mesh
// mpirun -np 4 ex4p -m ../data/beam-tet.mesh
// mpirun -np 4 ex4p -m ../data/beam-hex.mesh
// mpirun -np 4 ex4p -m ../data/escher.mesh -o 2 -sc
// mpirun -np 4 ex4p -m ../data/fichera.mesh -o 2 -hb
// mpirun -np 4 ex4p -m ../data/fichera-q2.vtk
// mpirun -np 4 ex4p -m ../data/fichera-q3.mesh -o 2 -sc
// mpirun -np 4 ex4p -m ../data/square-disc-nurbs.mesh -o 3
// mpirun -np 4 ex4p -m ../data/beam-hex-nurbs.mesh -o 3
// mpirun -np 4 ex4p -m ../data/periodic-square.mesh -no-bc
// mpirun -np 4 ex4p -m ../data/periodic-cube.mesh -no-bc
// mpirun -np 4 ex4p -m ../data/amr-quad.mesh
// mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -sc
// mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb
// mpirun -np 4 ex4p -m ../data/star-surf.mesh -o 3 -hb
//
// Description: This example code solves a simple 2D/3D H(div) diffusion
// problem corresponding to the second order definite equation
// -grad(alpha div F) + beta F = f with boundary condition F dot n
// = <given normal field>. Here, we use a given exact solution F
// and compute the corresponding r.h.s. f. We discretize with
// Raviart-Thomas finite elements.
//
// The example demonstrates the use of H(div) finite element
// spaces with the grad-div and H(div) vector finite element mass
// bilinear form, as well as the computation of discretization
// error when the exact solution is known. Bilinear form
// hybridization and static condensation are also illustrated.
//
// We recommend viewing examples 1-3 before viewing this example.
#include "mfem.hpp"
#include <fstream>
#include <iostream>
using namespace std;
using namespace mfem;
// Exact solution, F, and r.h.s., f. See below for implementation.
//void F_exact(const Vector &, Vector &);
//void f_exact(const Vector &, Vector &);
//double freq = 1.0, kappa;
void f_const(const Vector &, Vector &);
int dim;
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
int num_procs, myid;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
// 2. Parse command-line options.
const char *mesh_file = "../data/star.mesh";
int order = 1;
int rs = -1;
int rp = 2;
int ra = 0;
int bt = EntitySets::INVALID;
const char *bs = "Origin";
bool set_bc = true;
bool static_cond = false;
bool hybridization = false;
bool visualization = 1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree).");
args.AddOption(&set_bc, "-bc", "--impose-bc", "-no-bc", "--dont-impose-bc",
"Impose or not essential boundary conditions.");
args.AddOption(&rs, "-rs", "--refine-serial",
"Number of serial refinement levels");
args.AddOption(&rp, "-rp", "--refine-parallel",
"Number of parallel refinement levels");
args.AddOption(&ra, "-ra", "--refine-adaptive",
"Number of adaptive refinement levels");
args.AddOption(&bt, "-bt", "--bc-entity-type",
"");
args.AddOption(&bs, "-bs", "--bc-entity-set-name",
"");
// args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
// " solution.");
args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
"--no-static-condensation", "Enable static condensation.");
args.AddOption(&hybridization, "-hb", "--hybridization", "-no-hb",
"--no-hybridization", "Enable hybridization.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
// kappa = freq * M_PI;
// 3. Read the (serial) mesh from the given mesh file on all processors. We
// can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
// and volume, as well as periodic meshes with the same code.
Mesh *mesh = new Mesh(mesh_file, 1, 1);
dim = mesh->Dimension();
int sdim = mesh->SpaceDimension();
// 4. Refine the serial mesh on all processors to increase the resolution. In
// this example we do 'ref_levels' of uniform refinement. We choose
// 'ref_levels' to be the largest number that gives a final mesh with no
// more than 1,000 elements.
{
int ref_levels = ( rs >= 0 ) ? rs :
(int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
for (int l = 0; l < ref_levels; l++)
{
if ( myid == 0 ) { cout << "Uniform refinement in serial..."; }
mesh->UniformRefinement();
}
MPI_Barrier(MPI_COMM_WORLD);
if ( myid == 0 && rs > 0 ) { cout << "Done" << endl; }
}
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
/*
At this point we have a serial mesh containing an EntitySets
object which stores the current node/edge/face/element indices
for each entity in each set. This data is duplicated on each MPI
rank.
*/
if ( ra > 0 )
{
cout << "calling EnsureNCMesh" << endl;
mesh->EnsureNCMesh();
cout << "back from EnsureNCMesh" << endl;
}
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
/*
We now have an NCEntitySets object which stores the node indices
describing each enity in each node/edge/face set and the element
indices for the elements in each element set. This data is
duplicated on each MPI rank.
*/
// 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
// this mesh further in parallel to increase the resolution. Once the
// parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
// meshes need to be reoriented before we can define high-order Nedelec
// spaces on them (this is needed in the ADS solver below).
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
if ( pmesh->pent_sets )
{
cout << "pmesh->pent_sets is non NULL" << endl;
pmesh->pent_sets->PrintSetInfo(cout);
}
else
{
cout << "pmesh->pent_sets is NULL" << endl;
}
{
int par_ref_levels = rp;
for (int l = 0; l < par_ref_levels; l++)
{
if ( myid == 0 ) { cout << "Uniform refinement in parallel..."; }
pmesh->UniformRefinement();
}
}
pmesh->ReorientTetMesh();
for (int l = 0; l < ra; l++)
{
pmesh->RandomRefinement(0.2);
}
if ( ra > 0 )
{
if ( pmesh->pent_sets )
{
cout << "pmesh->pent_sets is non NULL post random refinement" << endl;
pmesh->pent_sets->PrintSetInfo(cout);
}
else
{
cout << "pmesh->pent_sets is NULL post random refinement" << endl;
}
}
// 6. Define a parallel finite element space on the parallel mesh. Here we
// use the Raviart-Thomas finite elements of the specified order.
FiniteElementCollection *fec = new RT_FECollection(order-1, dim);
ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
HYPRE_Int size = fespace->GlobalTrueVSize();
if (myid == 0)
{
cout << "Number of finite element unknowns: " << size << endl;
}
// 7. Determine the list of true (i.e. parallel conforming) essential
// boundary dofs. In this example, the boundary conditions are defined
// by marking all the boundary attributes from the mesh as essential
// (Dirichlet) and converting them to a list of true dofs.
Array<int> ess_tdof_list;
if ( bt == EntitySets::INVALID )
{
if (pmesh->bdr_attributes.Size())
{
Array<int> ess_bdr(pmesh->bdr_attributes.Max());
ess_bdr = set_bc ? 1 : 0;
fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
}
else
{
fespace->GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
ess_tdof_list);
}
if (myid == 0)
{
cout << "Number of Dirichlet dofs: " << ess_tdof_list.Size() << endl;
}
// 8. Set up the parallel linear form b(.) which corresponds to the
// right-hand side of the FEM linear system, which in this case is
// (f,phi_i) where f is given by the function f_exact and phi_i are the
// basis functions in the finite element fespace.
VectorFunctionCoefficient f(sdim, f_const);
ParLinearForm *b = new ParLinearForm(fespace);
b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
b->Assemble();
// 9. Define the solution vector x as a parallel finite element grid function
// corresponding to fespace. Initialize x by projecting the exact
// solution. Note that only values from the boundary faces will be used
// when eliminating the non-homogeneous boundary condition to modify the
// r.h.s. vector b.
ParGridFunction x(fespace);
// VectorFunctionCoefficient F(sdim, F_exact);
// x.ProjectCoefficient(F);
x = 0.0;
// 10. Set up the parallel bilinear form corresponding to the H(div)
// diffusion operator grad alpha div + beta I, by adding the div-div and
// the mass domain integrators.
Coefficient *alpha = new ConstantCoefficient(1.0);
Coefficient *beta = new ConstantCoefficient(1.0);
ParBilinearForm *a = new ParBilinearForm(fespace);
a->AddDomainIntegrator(new DivDivIntegrator(*alpha));
a->AddDomainIntegrator(new VectorFEMassIntegrator(*beta));
// 11. Assemble the parallel bilinear form and the corresponding linear
// system, applying any necessary transformations such as: parallel
// assembly, eliminating boundary conditions, applying conforming
// constraints for non-conforming AMR, static condensation,
// hybridization, etc.
FiniteElementCollection *hfec = NULL;
ParFiniteElementSpace *hfes = NULL;
if (static_cond)
{
a->EnableStaticCondensation();
}
else if (hybridization)
{
hfec = new DG_Interface_FECollection(order-1, dim);
hfes = new ParFiniteElementSpace(pmesh, hfec);
a->EnableHybridization(hfes, new NormalTraceJumpIntegrator(),
ess_tdof_list);
}
a->Assemble();
HypreParMatrix A;
Vector B, X;
a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
HYPRE_Int glob_size = A.GetGlobalNumRows();
if (myid == 0)
{
cout << "Size of linear system: " << glob_size << endl;
}
// 12. Define and apply a parallel PCG solver for A X = B with the 2D AMS or
// the 3D ADS preconditioners from hypre. If using hybridization, the
// system is preconditioned with hypre's BoomerAMG.
HypreSolver *prec = NULL;
CGSolver *pcg = new CGSolver(A.GetComm());
pcg->SetOperator(A);
pcg->SetRelTol(1e-12);
pcg->SetMaxIter(500);
pcg->SetPrintLevel(1);
if (hybridization) { prec = new HypreBoomerAMG(A); }
else
{
ParFiniteElementSpace *prec_fespace =
(a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace);
if (dim == 2) { prec = new HypreAMS(A, prec_fespace); }
else { prec = new HypreADS(A, prec_fespace); }
}
pcg->SetPreconditioner(*prec);
pcg->Mult(B, X);
// 13. Recover the parallel grid function corresponding to X. This is the
// local finite element solution on each processor.
a->RecoverFEMSolution(X, *b, x);
/*
// 14. Compute and print the L^2 norm of the error.
{
double err = x.ComputeL2Error(F);
if (myid == 0)
{
cout << "\n|| F_h - F ||_{L^2} = " << err << '\n' << endl;
}
}
*/
// 15. Save the refined mesh and the solution in parallel. This output can
// be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
{
ostringstream mesh_name, sol_name;
mesh_name << "mesh." << setfill('0') << setw(6) << myid;
sol_name << "sol." << setfill('0') << setw(6) << myid;
ofstream mesh_ofs(mesh_name.str().c_str());
mesh_ofs.precision(8);
pmesh->Print(mesh_ofs);
ofstream sol_ofs(sol_name.str().c_str());
sol_ofs.precision(8);
x.Save(sol_ofs);
}
// 16. Send the solution by socket to a GLVis server.
if (visualization)
{
char vishost[] = "localhost";
int visport = 19916;
socketstream sol_sock(vishost, visport);
sol_sock << "parallel " << num_procs << " " << myid << "\n";
sol_sock.precision(8);
sol_sock << "solution\n" << *pmesh << x << flush;
}
// 17. Free the used memory.
delete pcg;
delete prec;
delete hfes;
delete hfec;
delete a;
delete alpha;
delete beta;
delete b;
delete fespace;
delete fec;
delete pmesh;
MPI_Finalize();
return 0;
}
/*
// The exact solution (for non-surface meshes)
void F_exact(const Vector &p, Vector &F)
{
int dim = p.Size();
double x = p(0);
double y = p(1);
// double z = (dim == 3) ? p(2) : 0.0;
F(0) = cos(kappa*x)*sin(kappa*y);
F(1) = cos(kappa*y)*sin(kappa*x);
if (dim == 3)
{
F(2) = 0.0;
}
}
// The right hand side
void f_exact(const Vector &p, Vector &f)
{
int dim = p.Size();
double x = p(0);
double y = p(1);
// double z = (dim == 3) ? p(2) : 0.0;
double temp = 1 + 2*kappa*kappa;
f(0) = temp*cos(kappa*x)*sin(kappa*y);
f(1) = temp*cos(kappa*y)*sin(kappa*x);
if (dim == 3)
{
f(2) = 0;
}
}
*/
void f_const(const Vector &x, Vector &f)
{
if (dim == 3)
{
f(0) = 1.0;
f(1) = 1.0;
f(2) = 1.0;
}
else
{
f(0) = 1.0;
f(1) = 1.0;
if (x.Size() == 3) { f(2) = 0.0; }
}
}
-325
View File
@@ -1,325 +0,0 @@
// MFEM Example 6 - Parallel Version
//
// Compile with: make ex6p
//
// Sample runs: mpirun -np 4 ex6p -m ../data/square-disc.mesh -o 1
// mpirun -np 4 ex6p -m ../data/square-disc.mesh -o 2
// mpirun -np 4 ex6p -m ../data/square-disc-nurbs.mesh -o 2
// mpirun -np 4 ex6p -m ../data/star.mesh -o 3
// mpirun -np 4 ex6p -m ../data/escher.mesh -o 2
// mpirun -np 4 ex6p -m ../data/fichera.mesh -o 2
// mpirun -np 4 ex6p -m ../data/disc-nurbs.mesh -o 2
// mpirun -np 4 ex6p -m ../data/ball-nurbs.mesh
// mpirun -np 4 ex6p -m ../data/pipe-nurbs.mesh
// mpirun -np 4 ex6p -m ../data/star-surf.mesh -o 2
// mpirun -np 4 ex6p -m ../data/square-disc-surf.mesh -o 2
// mpirun -np 4 ex6p -m ../data/amr-quad.mesh
//
// Description: This is a version of Example 1 with a simple adaptive mesh
// refinement loop. The problem being solved is again the Laplace
// equation -Delta u = 1 with homogeneous Dirichlet boundary
// conditions. The problem is solved on a sequence of meshes which
// are locally refined in a conforming (triangles, tetrahedrons)
// or non-conforming (quadrilaterals, hexahedra) manner according
// to a simple ZZ error estimator.
//
// The example demonstrates MFEM's capability to work with both
// conforming and nonconforming refinements, in 2D and 3D, on
// linear, curved and surface meshes. Interpolation of functions
// from coarse to fine meshes, as well as persistent GLVis
// visualization are also illustrated.
//
// We recommend viewing Example 1 before viewing this example.
#include "mfem.hpp"
#include <fstream>
#include <iostream>
using namespace std;
using namespace mfem;
static int max_dofs = 100000;
int main(int argc, char *argv[])
{
// 1. Initialize MPI.
int num_procs, myid;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
// 2. Parse command-line options.
const char *mesh_file = "./star-set.mesh";
int order = 1;
int bt = EntitySets::INVALID;
const char *bs = "";
bool visualization = true;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&order, "-o", "--order",
"Finite element order (polynomial degree).");
args.AddOption(&max_dofs, "-md", "--max-dofs",
"Maximum number of degrees of freedom.");
args.AddOption(&bt, "-bt", "--bc-entity-type",
"");
args.AddOption(&bs, "-bs", "--bc-entity-set-name",
"");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.Parse();
if (!args.Good())
{
if (myid == 0)
{
args.PrintUsage(cout);
}
MPI_Finalize();
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
// 3. Read the (serial) mesh from the given mesh file on all processors. We
// can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
// and volume meshes with the same code.
Mesh *mesh = new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();
int sdim = mesh->SpaceDimension();
// 4. Refine the serial mesh on all processors to increase the resolution.
// Also project a NURBS mesh to a piecewise-quadratic curved mesh. Make
// sure that the mesh is non-conforming.
if (mesh->NURBSext)
{
mesh->UniformRefinement();
mesh->SetCurvature(2);
}
mesh->EnsureNCMesh();
if ( mesh->ent_sets )
{
cout << "mesh->ent_sets is non NULL" << endl;
mesh->ent_sets->PrintSetInfo(cout);
}
else
{
cout << "mesh->ent_sets is NULL" << endl;
}
// 5. Define a parallel mesh by partitioning the serial mesh.
// Once the parallel mesh is defined, the serial mesh can be deleted.
ParMesh pmesh(MPI_COMM_WORLD, *mesh);
delete mesh;
if ( pmesh.pent_sets )
{
cout << "pmesh->pent_sets is non NULL" << endl;
pmesh.pent_sets->PrintSetInfo(cout);
}
else
{
cout << "pmesh->pent_sets is NULL" << endl;
}
// 6. Define a finite element space on the mesh. The polynomial order is
// one (linear) by default, but this can be changed on the command line.
H1_FECollection fec(order, dim);
ParFiniteElementSpace fespace(&pmesh, &fec);
Array<int> ess_tdof_list;
if ( bt == EntitySets::INVALID )
{
if (pmesh.bdr_attributes.Size())
{
Array<int> ess_bdr(pmesh.bdr_attributes.Max());
ess_bdr = 1;
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
}
else
{
fespace.GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
ess_tdof_list);
}
// 7. As in Example 1p, we set up bilinear and linear forms corresponding to
// the Laplace problem -\Delta u = 1. We don't assemble the discrete
// problem yet, this will be done in the main loop.
ParBilinearForm a(&fespace);
ParLinearForm b(&fespace);
ConstantCoefficient one(1.0);
BilinearFormIntegrator *integ = new DiffusionIntegrator(one);
a.AddDomainIntegrator(integ);
b.AddDomainIntegrator(new DomainLFIntegrator(one));
// 8. The solution vector x and the associated finite element grid function
// will be maintained over the AMR iterations. We initialize it to zero.
ParGridFunction x(&fespace);
x = 0;
// 9. Connect to GLVis.
char vishost[] = "localhost";
int visport = 19916;
socketstream sout;
if (visualization)
{
sout.open(vishost, visport);
if (!sout)
{
if (myid == 0)
{
cout << "Unable to connect to GLVis server at "
<< vishost << ':' << visport << endl;
cout << "GLVis visualization disabled.\n";
}
visualization = false;
}
sout.precision(8);
}
// 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
// with L2 projection in the smoothing step to better handle hanging
// nodes and parallel partitioning. We need to supply a space for the
// discontinuous flux (L2) and a space for the smoothed flux (H(div) is
// used here).
L2_FECollection flux_fec(order, dim);
ParFiniteElementSpace flux_fes(&pmesh, &flux_fec, sdim);
RT_FECollection smooth_flux_fec(order-1, dim);
ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec);
// Another possible option for the smoothed flux space:
// H1_FECollection smooth_flux_fec(order, dim);
// ParFiniteElementSpace smooth_flux_fes(&pmesh, &smooth_flux_fec, dim);
L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fes, smooth_flux_fes);
// 11. A refiner selects and refines elements based on a refinement strategy.
// The strategy here is to refine elements with errors larger than a
// fraction of the maximum element error. Other strategies are possible.
// The refiner will call the given error estimator.
ThresholdRefiner refiner(estimator);
refiner.SetTotalErrorFraction(0.7);
// 12. The main AMR loop. In each iteration we solve the problem on the
// current mesh, visualize the solution, and refine the mesh.
// const int max_dofs = 100000;
for (int it = 0; ; it++)
{
HYPRE_Int global_dofs = fespace.GlobalTrueVSize();
if (myid == 0)
{
cout << "\nAMR iteration " << it << endl;
cout << "Number of unknowns: " << global_dofs << endl;
}
// 13. Assemble the stiffness matrix and the right-hand side. Note that
// MFEM doesn't care at this point that the mesh is nonconforming
// and parallel. The FE space is considered 'cut' along hanging
// edges/faces, and also across processor boundaries.
a.Assemble();
b.Assemble();
// 14. Create the parallel linear system: eliminate boundary conditions,
// constrain hanging nodes and nodes across processor boundaries.
// The system will be solved for true (unconstrained/unique) DOFs only.
// Array<int> ess_tdof_list;
if ( bt == EntitySets::INVALID )
{
if (pmesh.bdr_attributes.Size())
{
Array<int> ess_bdr(pmesh.bdr_attributes.Max());
ess_bdr = 1;
fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
}
}
else
{
fespace.GetEssentialTrueDofs((EntitySets::EntityType)bt, bs,
ess_tdof_list);
}
HypreParMatrix A;
Vector B, X;
const int copy_interior = 1;
a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
// 15. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
// preconditioner from hypre.
HypreBoomerAMG amg;
amg.SetPrintLevel(0);
CGSolver pcg(A.GetComm());
pcg.SetPreconditioner(amg);
pcg.SetOperator(A);
pcg.SetRelTol(1e-6);
pcg.SetMaxIter(200);
pcg.SetPrintLevel(3); // print the first and the last iterations only
pcg.Mult(B, X);
// 16. Extract the parallel grid function corresponding to the finite element
// approximation X. This is the local solution on each processor.
a.RecoverFEMSolution(X, b, x);
// 17. Send the solution by socket to a GLVis server.
if (visualization)
{
sout << "parallel " << num_procs << " " << myid << "\n";
sout << "solution\n" << pmesh << x << flush;
}
if (global_dofs > max_dofs)
{
if (myid == 0)
{
cout << "Reached the maximum number of dofs. Stop." << endl;
}
break;
}
// 18. Call the refiner to modify the mesh. The refiner calls the error
// estimator to obtain element errors, then it selects elements to be
// refined and finally it modifies the mesh. The Stop() method can be
// used to determine if a stopping criterion was met.
refiner.Apply(pmesh);
if (refiner.Stop())
{
if (myid == 0)
{
cout << "Stopping criterion satisfied. Stop." << endl;
}
break;
}
// 19. Update the finite element space (recalculate the number of DOFs,
// etc.) and create a grid function update matrix. Apply the matrix
// to any GridFunctions over the space. In this case, the update
// matrix is an interpolation matrix so the updated GridFunction will
// still represent the same function as before refinement.
fespace.Update();
x.Update();
// 20. Load balance the mesh, and update the space and solution. Currently
// available only for nonconforming meshes.
if (pmesh.Nonconforming())
{
pmesh.Rebalance();
// Update the space and the GridFunction. This time the update matrix
// redistributes the GridFunction among the processors.
fespace.Update();
x.Update();
}
// 21. Inform also the bilinear and linear forms that the space has
// changed.
a.Update();
b.Update();
}
MPI_Finalize();
return 0;
}
+2 -2
View File
@@ -282,10 +282,10 @@ int main(int argc, char *argv[])
delete b;
// 12. Compute and print the L^2 norm of the error.
double err = x.ComputeL2Error(sol_coef);
double error = x.ComputeL2Error(sol_coef);
if (myid == 0)
{
cout << "\nL2 norm of error: " << err << endl;
cout << "\nL2 norm of error: " << error << endl;
}
// 13. Save the refined mesh and the solution. This output can be viewed
-162
View File
@@ -1,162 +0,0 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
# TRIANGLE = 2
# SQUARE = 3
# TETRAHEDRON = 4
# CUBE = 5
# PRISM = 6
#
dimension
3
elements
14
1 4 13 15 21 25
1 4 12 13 15 21
1 4 13 21 22 25
1 4 15 24 21 25
1 4 13 15 25 16
1 5 0 1 4 3 9 10 13 12
1 5 8 9 12 11 17 18 21 20
1 5 2 3 6 5 11 12 15 14
1 6 3 4 6 12 13 15
1 6 4 7 6 13 16 15
1 6 12 13 21 9 10 18
1 6 13 22 21 10 19 18
1 6 11 14 20 12 15 21
1 6 15 21 24 14 20 23
boundary
30
1 3 5 6 3 2
2 2 3 6 4
2 2 4 6 7
3 3 3 4 1 0
4 3 11 12 9 8
5 3 2 3 12 11
6 3 0 1 10 9
7 2 9 10 18
7 2 10 19 18
8 3 8 9 18 17
9 3 1 4 13 10
10 3 4 7 16 13
11 2 13 16 25
11 2 13 25 22
12 3 10 13 22 19
13 3 7 6 15 16
14 3 6 5 14 15
15 3 15 14 23 24
16 2 16 15 25
16 2 15 24 25
17 3 5 2 11 14
18 3 3 0 9 12
19 3 11 8 17 20
20 2 11 20 14
20 2 14 20 23
21 3 17 18 21 20
22 3 18 19 22 21
23 2 21 22 25
23 2 21 25 24
24 3 20 21 24 23
vertices
26
3
0 -1 -1
1 -1 -1
-1 0 -1
0 0 -1
1 0 -1
-1 1 -1
0 1 -1
1 1 -1
-1 -1 0
0 -1 0
1 -1 0
-1 0 0
0 0 0
1 0 0
-1 1 0
0 1 0
1 1 0
-1 -1 1
0 -1 1
1 -1 1
-1 0 1
0 0 1
1 0 1
-1 1 1
0 1 1
1 1 1
MFEM sets v1.0
vertex_sets
1
Origin
1
12
edge_sets
2
Axes
3
12 13
12 15
12 21
Negative Axes
3
12 9
12 11
12 3
face_sets
2
Interior Corner
3
3 11 12 9 8
3 2 3 12 11
3 3 0 9 12
Exterior Corner
15
2 13 16 25
2 13 25 22
2 16 15 25
2 15 24 25
2 21 22 25
2 21 25 24
3 10 13 22 19
3 4 7 16 13
3 1 4 13 10
3 7 6 15 16
3 6 5 14 15
3 15 14 23 24
3 20 21 24 23
3 18 19 22 21
3 17 18 21 20
element_sets
3
Interior Corner
3
5 6 7
Exterior Corner
5
0 1 2 3 4
Steps
3
6 8 9
-145
View File
@@ -1,145 +0,0 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
# TRIANGLE = 2
# SQUARE = 3
# TETRAHEDRON = 4
# CUBE = 5
#
dimension
3
elements
7
1 5 0 1 4 3 9 10 13 12
1 5 3 4 7 6 12 13 16 15
1 5 2 3 6 5 11 12 15 14
1 5 8 9 12 11 17 18 21 20
1 5 9 10 13 12 18 19 22 21
1 5 12 13 16 15 21 22 25 24
1 5 11 12 15 14 20 21 24 23
boundary
24
1 3 5 6 3 2
2 3 6 7 4 3
3 3 3 4 1 0
4 3 11 12 9 8
5 3 2 3 12 11
6 3 0 1 10 9
7 3 9 10 19 18
8 3 8 9 18 17
9 3 1 4 13 10
10 3 4 7 16 13
11 3 13 16 25 22
12 3 10 13 22 19
13 3 7 6 15 16
14 3 6 5 14 15
15 3 15 14 23 24
16 3 16 15 24 25
17 3 5 2 11 14
18 3 3 0 9 12
19 3 11 8 17 20
20 3 14 11 20 23
21 3 17 18 21 20
22 3 18 19 22 21
23 3 21 22 25 24
24 3 20 21 24 23
vertices
26
3
0 -1 -1
1 -1 -1
-1 0 -1
0 0 -1
1 0 -1
-1 1 -1
0 1 -1
1 1 -1
-1 -1 0
0 -1 0
1 -1 0
-1 0 0
0 0 0
1 0 0
-1 1 0
0 1 0
1 1 0
-1 -1 1
0 -1 1
1 -1 1
-1 0 1
0 0 1
1 0 1
-1 1 1
0 1 1
1 1 1
MFEM sets v1.0
vertex_sets
1
Origin
1
12
edge_sets
2
Axes
3
12 13
12 15
12 21
Negative Axes
3
12 9
12 11
12 3
face_sets
2
Interior Corner
3
3 11 12 9 8
3 2 3 12 11
3 3 0 9 12
Exterior Corner
12
3 13 16 25 22
3 16 15 24 25
3 21 22 25 24
3 10 13 22 19
3 4 7 16 13
3 1 4 13 10
3 7 6 15 16
3 6 5 14 15
3 15 14 23 24
3 20 21 24 23
3 18 19 22 21
3 17 18 21 20
element_sets
3
Interior Corner
3
0 2 3
Exterior Corner
1
5
Steps
2
1 3
-158
View File
@@ -1,158 +0,0 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
# TRIANGLE = 2
# SQUARE = 3
# TETRAHEDRON = 4
# CUBE = 5
# PRISM = 6
#
dimension
2
elements
30
1 3 0 11 26 14
1 3 0 14 27 17
1 3 0 17 28 20
1 3 0 20 29 23
1 3 0 23 30 11
1 2 11 1 26
1 2 1 12 26
1 3 26 12 3 13
1 2 26 13 2
1 2 14 26 2
1 2 14 2 27
1 2 2 15 27
1 3 27 15 5 16
1 2 27 16 4
1 2 17 27 4
1 2 17 4 28
1 2 4 18 28
1 3 28 18 7 19
1 2 28 19 6
1 2 20 28 6
1 2 20 6 29
1 2 6 21 29
1 3 29 21 9 22
1 2 29 22 8
1 2 23 29 8
1 2 23 8 30
1 2 8 24 30
1 3 30 24 10 25
1 2 30 25 1
1 2 11 30 1
boundary
20
1 1 13 2
1 1 12 3
1 1 16 4
1 1 15 5
1 1 19 6
1 1 18 7
1 1 22 8
1 1 21 9
1 1 25 1
1 1 24 10
1 1 3 13
1 1 1 12
1 1 5 16
1 1 2 15
1 1 7 19
1 1 4 18
1 1 9 22
1 1 6 21
1 1 10 25
1 1 8 24
vertices
31
2
0 0
1 0
0.309017 0.951057
1.30902 0.951057
-0.809017 0.587785
-0.5 1.53884
-0.809017 -0.587785
-1.61803 0
0.309017 -0.951057
-0.5 -1.53884
1.30902 -0.951057
0.5 0
1.15451 0.475529
0.809019 0.951057
0.154508 0.475529
-0.0954915 1.24495
-0.654508 1.06331
-0.404508 0.293893
-1.21352 0.293893
-1.21352 -0.293892
-0.404508 -0.293893
-0.654508 -1.06331
-0.0954915 -1.24495
0.154508 -0.475529
0.809019 -0.951057
1.15451 -0.475529
0.654509 0.475529
-0.25 0.769421
-0.809016 0
-0.25 -0.76942
0.654509 -0.475529
MFEM sets v1.0
vertex_sets
3
Origin
1
0
Tent
5
1 2 4 6 8
Gazebo
5
3 5 7 9 10
edge_sets
2
Columbine
5
1 11
2 14
4 17
6 20
8 23
Lily
5
0 11
0 14
0 17
0 20
0 23
element_sets
3
Flying Squirrel
3
7 17 27
Sea Lion
4
12 17 22 27
Pinwheel
5
8 13 18 23 28
-143
View File
@@ -1,143 +0,0 @@
MFEM mesh v1.0
#
# MFEM Geometry Types (see mesh/geom.hpp):
#
# POINT = 0
# SEGMENT = 1
# TRIANGLE = 2
# SQUARE = 3
# TETRAHEDRON = 4
# CUBE = 5
#
dimension
2
elements
20
1 3 0 11 26 14
1 3 0 14 27 17
1 3 0 17 28 20
1 3 0 20 29 23
1 3 0 23 30 11
1 3 11 1 12 26
1 3 26 12 3 13
1 3 14 26 13 2
1 3 14 2 15 27
1 3 27 15 5 16
1 3 17 27 16 4
1 3 17 4 18 28
1 3 28 18 7 19
1 3 20 28 19 6
1 3 20 6 21 29
1 3 29 21 9 22
1 3 23 29 22 8
1 3 23 8 24 30
1 3 30 24 10 25
1 3 11 30 25 1
boundary
20
1 1 13 2
1 1 12 3
1 1 16 4
1 1 15 5
1 1 19 6
1 1 18 7
1 1 22 8
1 1 21 9
1 1 25 1
1 1 24 10
1 1 3 13
1 1 1 12
1 1 5 16
1 1 2 15
1 1 7 19
1 1 4 18
1 1 9 22
1 1 6 21
1 1 10 25
1 1 8 24
vertices
31
2
0 0
1 0
0.309017 0.951057
1.30902 0.951057
-0.809017 0.587785
-0.5 1.53884
-0.809017 -0.587785
-1.61803 0
0.309017 -0.951057
-0.5 -1.53884
1.30902 -0.951057
0.5 0
1.15451 0.475529
0.809019 0.951057
0.154508 0.475529
-0.0954915 1.24495
-0.654508 1.06331
-0.404508 0.293893
-1.21352 0.293893
-1.21352 -0.293892
-0.404508 -0.293893
-0.654508 -1.06331
-0.0954915 -1.24495
0.154508 -0.475529
0.809019 -0.951057
1.15451 -0.475529
0.654509 0.475529
-0.25 0.769421
-0.809016 0
-0.25 -0.76942
0.654509 -0.475529
MFEM sets v1.0
vertex_sets
3
Origin
1
0
Tent
5
1 2 4 6 8
Gazebo
5
3 5 7 9 10
edge_sets
2
Columbine
5
1 11
2 14
4 17
6 20
8 23
Lily
5
0 11
0 14
0 17
0 20
0 23
element_sets
2
Flying Squirrel
3
6 12 18
Sea Lion
4
9 12 15 18
+98 -7
View File
@@ -21,6 +21,33 @@ namespace mfem
using namespace std;
// Given an ElementTransformation and IntegrationPoint in a refined mesh,
// return the ElementTransformation of the parent coarse element, and set
// coarse_ip to the location of the original ip within the coarse element.
ElementTransformation *RefinedToCoarse(
Mesh &coarse_mesh, const ElementTransformation &T,
const IntegrationPoint &ip, IntegrationPoint &coarse_ip)
{
Mesh &fine_mesh = *T.mesh;
// Get the element transformation of the coarse element containing the
// fine element.
int fine_element = T.ElementNo;
const CoarseFineTransformations &cf = fine_mesh.GetRefinementTransforms();
int coarse_element = cf.embeddings[fine_element].parent;
ElementTransformation *coarse_T = coarse_mesh.GetElementTransformation(
coarse_element);
// Transform the integration point from fine element coordinates to coarse
// element coordinates.
Geometry::Type geom = T.GetGeometryType();
IntegrationPointTransformation fine_to_coarse;
IsoparametricTransformation &emb_tr = fine_to_coarse.Transf;
emb_tr.SetIdentityTransformation(geom);
emb_tr.SetPointMat(cf.point_matrices[geom](cf.embeddings[fine_element].matrix));
fine_to_coarse.Transform(ip, coarse_ip);
coarse_T->SetIntPoint(&coarse_ip);
return coarse_T;
}
double PWConstCoefficient::Eval(ElementTransformation & T,
const IntegrationPoint & ip)
{
@@ -95,7 +122,17 @@ double FunctionCoefficient::Eval(ElementTransformation & T,
double GridFunctionCoefficient::Eval (ElementTransformation &T,
const IntegrationPoint &ip)
{
return GridF -> GetValue (T, ip, Component);
Mesh *gf_mesh = GridF->FESpace()->GetMesh();
if (T.mesh == gf_mesh)
{
return GridF->GetValue(T, ip, Component);
}
else
{
IntegrationPoint coarse_ip;
ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
return GridF->GetValue(*coarse_T, coarse_ip, Component);
}
}
void TransformedCoefficient::SetTime(double t)
@@ -305,13 +342,30 @@ void VectorGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
void VectorGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip)
{
GridFunc->GetVectorValue(T, ip, V);
Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
if (T.mesh == gf_mesh)
{
GridFunc->GetVectorValue(T, ip, V);
}
else
{
IntegrationPoint coarse_ip;
ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
GridFunc->GetVectorValue(*coarse_T, coarse_ip, V);
}
}
void VectorGridFunctionCoefficient::Eval(
DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
{
GridFunc->GetVectorValues(T, ir, M);
if (T.mesh == GridFunc->FESpace()->GetMesh())
{
GridFunc->GetVectorValues(T, ir, M);
}
else
{
VectorCoefficient::Eval(M, T, ir);
}
}
GradientGridFunctionCoefficient::GradientGridFunctionCoefficient (
@@ -331,13 +385,30 @@ void GradientGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
void GradientGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip)
{
GridFunc->GetGradient(T, V);
Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
if (T.mesh == gf_mesh)
{
GridFunc->GetGradient(T, V);
}
else
{
IntegrationPoint coarse_ip;
ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
GridFunc->GetGradient(*coarse_T, V);
}
}
void GradientGridFunctionCoefficient::Eval(
DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
{
GridFunc->GetGradients(T, ir, M);
if (T.mesh == GridFunc->FESpace()->GetMesh())
{
GridFunc->GetGradients(T, ir, M);
}
else
{
VectorCoefficient::Eval(M, T, ir);
}
}
CurlGridFunctionCoefficient::CurlGridFunctionCoefficient(
@@ -363,7 +434,17 @@ void CurlGridFunctionCoefficient::SetGridFunction(const GridFunction *gf)
void CurlGridFunctionCoefficient::Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip)
{
GridFunc->GetCurl(T, V);
Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
if (T.mesh == gf_mesh)
{
GridFunc->GetCurl(T, V);
}
else
{
IntegrationPoint coarse_ip;
ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
GridFunc->GetCurl(*coarse_T, V);
}
}
DivergenceGridFunctionCoefficient::DivergenceGridFunctionCoefficient (
@@ -375,7 +456,17 @@ DivergenceGridFunctionCoefficient::DivergenceGridFunctionCoefficient (
double DivergenceGridFunctionCoefficient::Eval(ElementTransformation &T,
const IntegrationPoint &ip)
{
return GridFunc->GetDivergence(T);
Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
if (T.mesh == gf_mesh)
{
return GridFunc->GetDivergence(T);
}
else
{
IntegrationPoint coarse_ip;
ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
return GridFunc->GetDivergence(*coarse_T);
}
}
void VectorDeltaCoefficient::SetTime(double t)
+1
View File
@@ -574,6 +574,7 @@ public:
/// Evaluate the coefficient.
virtual void Eval(Vector &V, ElementTransformation &T,
const IntegrationPoint &ip);
using VectorCoefficient::Eval;
};
/// A general vector function coefficient
+1 -1
View File
@@ -330,7 +330,7 @@ ConduitDataCollection::BlueprintMeshToMesh(const Node &n_mesh,
}
else
{
Node &(n_bndry_conn_conv) =
Node &n_bndry_conn_conv =
n_conv["topologies"][bndry_topo_name]["elements/connectivity"];
n_bndry_conn.to_int_array(n_bndry_conn_conv);
bndry_indices = (n_bndry_conn_conv).value();
+1
View File
@@ -243,6 +243,7 @@ public:
void TransformDual(double *v) const;
void InvTransformDual(double *v) const;
using DofTransformation::InvTransformDual;
};
/// DoF transformation implementation for the Nedelec basis on tetrahedra
+2 -1
View File
@@ -21,7 +21,8 @@ ElementTransformation::ElementTransformation()
EvalState(0),
geom(Geometry::INVALID),
Attribute(-1),
ElementNo(-1)
ElementNo(-1),
mesh(nullptr)
{ }
double ElementTransformation::EvalWeight()
+6
View File
@@ -75,6 +75,12 @@ public:
int Attribute, ElementNo, ElementType;
/// The Mesh object containing the element.
/** If the element transformation belongs to a mesh, this will point to the
containing Mesh object. ElementNo will be the number of the element in
this Mesh. This will be NULL if the element does not belong to a mesh. */
class Mesh *mesh;
ElementTransformation();
/** @brief Force the reevaluation of the Jacobian in the next call. */
+51
View File
@@ -1993,6 +1993,57 @@ void Poly_1D::CalcDBinomTerms(const int p, const double x, const double y,
}
}
void Poly_1D::CalcJacobi(const int p, const double alpha, const double beta,
const double x, double *u)
{
// use the recursive definition for [-1,1]:
// 2(n+1)(n+a+b+1)(2n+a+b)P^{(a,b)}_{n+1}(z) =
// (2n+a+b+1)((2n+a+b+2)(2n+a+b)z+a^2-b^2)P^{(a,b)}_n(z)
// -2(n+a)(n+b)(2n+a+b+2)P^{(a,b)}_{n-1}(z)
u[0] = 1.;
if (p == 0) { return; }
double z = 2.*x - 1.;
u[1] = 0.5 * ((alpha + beta + 2.) * z + alpha - beta);
for (int n = 1; n < p; n++)
{
double c2 = 2.*(1. + n)*(1. + n + alpha + beta)*(2.*n + alpha + beta);
double c1 = (1. + 2.*n + alpha + beta) *
((2. + 2.*n + alpha + beta)*
(2.*n + alpha + beta)*z + pow(alpha,2) - pow(beta,2));
double c0 = 2.*(alpha + n)*(beta + n)*(2. + 2.*n + alpha + beta);
u[n+1] = (c1 * u[n] - c0 * u[n-1]) / c2;
}
}
void Poly_1D::CalcJacobi(const int p, const double alpha, const double beta,
const double x, double *u, double *d)
{
// use the recursive definition for [-1,1]:
// 2(n+1)(n+a+b+1)(2n+a+b)P^{(a,b)}_{n+1}(z) =
// (2n+a+b+1)((2n+a+b+2)(2n+a+b)z+a^2-b^2)P^{(a,b)}_n(z)
// -2(n+a)(n+b)(2n+a+b+2)P^{(a,b)}_{n-1}(z)
u[0] = 1.;
d[0] = 0.;
if (p == 0) { return; }
double z = 2.*x - 1.;
u[1] = 0.5 * ((alpha + beta + 2.) * z + alpha - beta);
d[1] = alpha + beta + 2.;
for (int n = 1; n < p; n++)
{
double c2 = 2.*(1. + n)*(1. + n + alpha + beta)*(2.*n + alpha + beta);
double c1 = (1. + 2.*n + alpha + beta) *
((2. + 2.*n + alpha + beta)*
(2.*n + alpha + beta)*z + pow(alpha,2) - pow(beta,2));
double c0 = 2.*(alpha + n)*(beta + n)*(2. + 2.*n + alpha + beta);
double dc1 = 2.*(1. + 2.*n + alpha + beta) *
(2. + 2.*n + alpha + beta)*(2.*n + alpha + beta);
u[n+1] = (c1 * u[n] - c0 * u[n-1]) / c2;
d[n+1] = (c1 * d[n] + dc1 * u[n] - c0 * d[n-1]) / c2;
}
}
void Poly_1D::CalcLegendre(const int p, const double x, double *u)
{
// use the recursive definition for [-1,1]:
+21
View File
@@ -1112,7 +1112,28 @@ public:
static void CalcBernstein(const int p, const double x, double *u, double *d)
{ CalcBinomTerms(p, x, 1. - x, u, d); }
// Evaluate the values of the Jacobi polynomials, P^{(\alpha,\beta)}_n,
// on [0,1] with degrees ranging from 0 to p and parameters alpha and beta
// at a point x. Internally the point is transformed to the traditional
// domain of [-1,1].
static void CalcJacobi(const int p, const double alpha, const double beta,
const double x, double *u);
// Evaluate the values and derivatives of the Jacobi polynomials,
// P^{(\alpha,\beta)}_n, on [0,1] with degrees ranging from 0 to p and
// parameters alpha and beta at a point x. Internally the point is
// transformed to the traditional domain of [-1,1].
static void CalcJacobi(const int p, const double alpha, const double beta,
const double x, double *u, double *d);
// Evaluate the values of the Legendre polynomials, P_n, on [0,1] with
// degrees ranging from 0 to p at a point x. Internally the point is
// transformed to the traditional domain of [-1,1].
static void CalcLegendre(const int p, const double x, double *u);
// Evaluate the values and derivatives of the Legendre polynomials, P_n,
// on [0,1] with degrees ranging from 0 to p at a point x. Internally the
// point is transformed to the traditional domain of [-1,1].
static void CalcLegendre(const int p, const double x, double *u, double *d);
~Poly_1D();
+311
View File
@@ -1040,4 +1040,315 @@ void H1_WedgeElement::CalcDShape(const IntegrationPoint &ip,
}
}
H1_PyramidElement::H1_PyramidElement(const int p, const int btype)
: NodalFiniteElement(3, Geometry::PYRAMID,
(p + 1) * (p + 2) * (2 * p + 3) / 6,
p, FunctionSpace::Qk)
{
const double *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
const double **mcp = new const double*[p];
mcp[0] = NULL;
for (int k=1; k<p; k++)
{
mcp[k] = poly1d.ClosedPoints(k, VerifyNodal(VerifyClosed(btype)));
}
#ifndef MFEM_THREAD_SAFE
shape_x.SetSize(p + 1);
shape_y.SetSize(p + 1);
shape_z.SetSize(p + 1);
dshape_x.SetSize(p + 1);
dshape_y.SetSize(p + 1);
dshape_z.SetSize(p + 1);
ddshape_x.SetSize(p + 1);
ddshape_y.SetSize(p + 1);
ddshape_z.SetSize(p + 1);
u.SetSize(dof);
du.SetSize(dof, dim);
ddu.SetSize(dof, (dim * (dim + 1)) / 2);
#else
Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1);
#endif
// vertices
Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
Nodes.IntPoint(2).Set3(cp[p], cp[p], cp[0]);
Nodes.IntPoint(3).Set3(cp[0], cp[p], cp[0]);
Nodes.IntPoint(4).Set3(cp[0], cp[0], cp[p]);
// edges
int o = 5;
for (int i = 1; i < p; i++) // (0,1)
{
Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
}
for (int i = 1; i < p; i++) // (1,2)
{
Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[0]);
}
for (int i = 1; i < p; i++) // (3,2)
{
Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[0]);
}
for (int i = 1; i < p; i++) // (0,3)
{
Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
}
for (int i = 1; i < p; i++) // (0,4)
{
Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
}
for (int i = 1; i < p; i++) // (1,4)
{
Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
}
for (int i = 1; i < p; i++) // (2,4)
{
Nodes.IntPoint(o++).Set3(cp[p-i], cp[p-i], cp[i]);
}
for (int i = 1; i < p; i++) // (3,4)
{
Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
}
// quadrilateral face
for (int j = 1; j < p; j++)
{
for (int i = 1; i < p; i++)
{
Nodes.IntPoint(o++).Set3(cp[i], cp[j], cp[0]);
}
}
// triangular faces
for (int j = 1; j < p; j++)
for (int i = 1; i + j < p; i++) // (0,1,4)
{
double w = cp[i] + cp[j] + cp[p-i-j];
Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
}
for (int j = 1; j < p; j++)
for (int i = 1; i + j < p; i++) // (1,2,4)
{
double w = cp[i] + cp[j] + cp[p-i-j];
Nodes.IntPoint(o++).Set3(1.0 - cp[j]/w, cp[i]/w, cp[j]/w);
}
for (int j = 1; j < p; j++)
for (int i = 1; i + j < p; i++) // (3,4,2)
{
double w = cp[i] + cp[j] + cp[p-i-j];
Nodes.IntPoint(o++).Set3(cp[j]/w, 1.0 - cp[i]/w, cp[i]/w);
}
for (int j = 1; j < p; j++)
for (int i = 1; i + j < p; i++) // (0,4,3)
{
double w = cp[i] + cp[j] + cp[p-i-j];
Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
}
// interior
for (int k = 1; k < p - 1; k++)
{
double wk = 1.0 - cp[k];
for (int j = 1; j < p - k; j++)
{
for (int i = 1; i < p - k; i++)
{
Nodes.IntPoint(o++).Set3(mcp[p-k][i] * wk,
mcp[p-k][j] * wk,
cp[k]);
}
}
}
MFEM_ASSERT(o == dof,
"Number of nodes does not match the "
"number of degrees of freedom");
DenseMatrix T(dof);
for (int m = 0; m < dof; m++)
{
const IntegrationPoint &ip = Nodes.IntPoint(m);
double oz = 1.0 - ip.z;
double tol = 1e-6;
double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
double z = ip.z;
poly1d.CalcLegendre(p, x, shape_x);
poly1d.CalcLegendre(p, y, shape_y);
o = 0;
for (int i = 0; i <= p; i++)
{
for (int j = 0; j <= p; j++)
{
int maxij = std::max(i, j);
poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z, shape_z);
for (int k = 0; k <= p - maxij; k++)
{
if (oz <= tol)
{
if (maxij == 0)
{
T(o++, m) = shape_z(k);
}
else
{
T(o++, m) = 0.0;
}
}
else
{
T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
pow(1.0 - ip.z, maxij);
}
}
}
}
}
Ti.Factor(T);
delete [] mcp;
}
void H1_PyramidElement::CalcShape(const IntegrationPoint &ip,
Vector &shape) const
{
const int p = order;
#ifdef MFEM_THREAD_SAFE
Vector shape_x(order+1);
Vector shape_y(order+1);
Vector shape_z(order+1);
Vector u(dof);
#endif
double oz = 1.0 - ip.z;
double tol = 1e-6;
double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
double z = ip.z;
poly1d.CalcLegendre(p, x, shape_x);
poly1d.CalcLegendre(p, y, shape_y);
for (int o = 0, i = 0; i <= p; i++)
{
for (int j = 0; j <= p; j++)
{
int maxij = std::max(i, j);
poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z, shape_z);
for (int k = 0; k <= p - maxij; k++)
{
if (oz <= tol)
{
if (maxij == 0)
{
u(o++) = shape_z(k);
}
else
{
u(o++) = 0.0;
}
}
else
{
u(o++) =
shape_x(i) * shape_y(j) * shape_z(k) * pow(1.0 - ip.z, maxij);
}
}
}
}
Ti.Mult(u, shape);
}
void H1_PyramidElement::CalcDShape(const IntegrationPoint &ip,
DenseMatrix &dshape) const
{
const int p = order;
#ifdef MFEM_THREAD_SAFE
Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1);
Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1);
DenseMatrix du(dof, dim);
#endif
double oz = 1.0 - ip.z;
double tol = 1e-6;
if (oz <= tol)
{
mfem::out << "H1_PyramidElement::CalcDShape near apex" << std::endl;
}
double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
double z = ip.z;
poly1d.CalcLegendre(p, x, shape_x, dshape_x);
poly1d.CalcLegendre(p, y, shape_y, dshape_y);
for (int o = 0, i = 0; i <= p; i++)
{
for (int j = 0; j <= p; j++)
{
int maxij = std::max(i, j);
poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z,
shape_z, dshape_z);
for (int k = 0; k <= p - maxij; k++)
{
if (oz > tol)
{
du(o, 0) = dshape_x(i) * shape_y(j) * shape_z(k) *
pow(1.0 - ip.z, maxij - 1);
du(o, 1) = shape_x(i) * dshape_y(j) * shape_z(k) *
pow(1.0 - ip.z, maxij - 1);
du(o, 2) = (ip.x * dshape_x(i) * shape_y(j) * shape_z(k) +
ip.y * shape_x(i) * dshape_y(j) * shape_z(k) -
shape_x(i) * shape_y(j) * (1.0 - ip.z) *
(shape_z(k) * maxij -
(1.0 - ip.z) * dshape_z(k))
) * pow(1.0 - ip.z, maxij - 2);
}
else
{
// The gradient is multivalued at the apex. The following
// values were produced by setting x and y equal to zero
// and computing the limit as z->1.
du(o, 0) = 0.0;
du(o, 1) = 0.0;
du(o, 2) = 0.0;
if (i == 0 && j == 0)
{
du(o, 2) = k * (k + 1) * (k + 2) * (k + 3) / 6;
}
else if (i == 0 && j == 1)
{
du(o, 1) = -(k + 1) * (k + 2) * (k + 3) * (k + 4) / 24;
}
else if (i == 1 && j == 0)
{
du(o, 0) = -(k + 1) * (k + 2) * (k + 3) * (k + 4) / 24;
}
}
o++;
}
}
}
Ti.Mult(du, dshape);
}
}
+25
View File
@@ -148,6 +148,31 @@ public:
DenseMatrix &dshape) const;
};
class H1_PyramidElement : public NodalFiniteElement
{
private:
#ifndef MFEM_THREAD_SAFE
// mutable Vector shape_0, shape_1, shape_2;
// mutable Vector dshape_0_0, dshape_1_0, dshape_2_0;
// mutable Vector dshape_0_1, dshape_1_1, dshape_2_1;
// mutable Vector u;
// mutable DenseMatrix du;
mutable Vector shape_x, shape_y, shape_z;
mutable Vector dshape_x, dshape_y, dshape_z, u;
mutable Vector ddshape_x, ddshape_y, ddshape_z;
mutable DenseMatrix du, ddu;
#endif
DenseMatrixInverse Ti;
public:
H1_PyramidElement(const int p,
const int btype = BasisType::GaussLobatto);
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
virtual void CalcDShape(const IntegrationPoint &ip,
DenseMatrix &dshape) const;
};
} // namespace mfem
#endif
+3 -12
View File
@@ -1749,7 +1749,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
H1_dof[Geometry::CUBE] = QuadDof*pm1;
H1_dof[Geometry::PRISM] = TriDof*pm1;
H1_dof[Geometry::PYRAMID] = 0;
H1_dof[Geometry::PYRAMID] = pm1*pm2*(2*p-3)/6;
if (b_type == BasisType::Positive)
{
H1_Elements[Geometry::TETRAHEDRON] = new H1Pos_TetrahedronElement(p);
@@ -1763,7 +1763,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p, btype);
H1_Elements[Geometry::PRISM] = new H1_WedgeElement(p, btype);
}
H1_Elements[Geometry::PYRAMID] = new LinearPyramidFiniteElement;
H1_Elements[Geometry::PYRAMID] = new H1_PyramidElement(p, btype);
const int &TetDof = H1_dof[Geometry::TETRAHEDRON];
TetDofOrd[0] = new int[24*TetDof];
@@ -1861,16 +1861,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
const FiniteElement *
H1_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
{
if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
{
return H1_Elements[GeomType];
}
else
{
MFEM_ABORT("H1 Pyramid basis functions are not yet supported "
"for order > 1.");
return NULL;
}
return H1_Elements[GeomType];
}
const int *H1_FECollection::DofOrderForOrientation(Geometry::Type GeomType,
-179
View File
@@ -561,155 +561,6 @@ void FiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
}
}
void FiniteElementSpace::GetEssentialVDofs(EntitySets::EntityType type,
int set_index,
Array<int> &ess_vdofs,
int component) const
{
Array<int> vdofs, dofs;
ess_vdofs.SetSize(GetVSize());
ess_vdofs = 0;
MFEM_VERIFY(mesh->ent_sets != NULL, "Mesh object contains no "
"entity set information");
if (!mesh->ent_sets->SetExists(type, set_index))
{
ostringstream oss; oss << "Entity set of type \""
<< EntitySets::GetTypeName(type)
<< "\" and index " << set_index
<< " was not found.";
MFEM_VERIFY(false, oss.str().c_str());
}
set<int>::iterator it;
for (it=(*mesh->ent_sets)(type, set_index).begin();
it!=(*mesh->ent_sets)(type, set_index).end(); it++)
{
int ent_index = *it;
cout << "collecting vdofs for entity " << ent_index << "->";
if (component < 0)
{
switch (type)
{
case EntitySets::VERTEX:
GetVertexVDofs(ent_index, vdofs);
break;
case EntitySets::EDGE:
GetEdgeVDofs(ent_index, vdofs);
break;
case EntitySets::FACE:
GetFaceVDofs(ent_index, vdofs);
break;
case EntitySets::ELEMENT:
GetElementVDofs(ent_index, vdofs);
break;
default:
mfem_error("GetEssentialVDofs: Invalid entity type");
}
vdofs.Print(cout);
mark_dofs(vdofs, ess_vdofs);
}
else
{
switch (type)
{
case EntitySets::VERTEX:
GetVertexDofs(ent_index, dofs);
break;
case EntitySets::EDGE:
GetEdgeDofs(ent_index, dofs);
break;
case EntitySets::FACE:
GetFaceDofs(ent_index, dofs);
break;
case EntitySets::ELEMENT:
GetElementDofs(ent_index, dofs);
break;
default:
mfem_error("GetEssentialDofs: Invalid entity type");
}
for (int d = 0; d < dofs.Size(); d++)
{ dofs[d] = DofToVDof(dofs[d], component); }
mark_dofs(dofs, ess_vdofs);
}
}
if (mesh->ncmesh)
{
Array<int> es_verts, es_edges, es_faces;
mesh->ncmesh->GetEntitySetClosure(type, set_index,
es_verts, es_edges, es_faces);
cout << "returned from get closure" << endl;
for (int i = 0; i < es_verts.Size(); i++)
{
if (es_verts[i] < GetNV())
{
if (component < 0)
{
GetVertexVDofs(es_verts[i], vdofs);
mark_dofs(vdofs, ess_vdofs);
}
else
{
GetVertexDofs(es_verts[i], dofs);
for (int d = 0; d < dofs.Size(); d++)
{ dofs[d] = DofToVDof(dofs[d], component); }
mark_dofs(dofs, ess_vdofs);
}
}
}
for (int i = 0; i < es_edges.Size(); i++)
{
if (es_edges[i] < GetMesh()->GetNEdges())
{
if (component < 0)
{
GetEdgeVDofs(es_edges[i], vdofs);
mark_dofs(vdofs, ess_vdofs);
}
else
{
GetEdgeDofs(es_edges[i], dofs);
for (int d = 0; d < dofs.Size(); d++)
{ dofs[d] = DofToVDof(dofs[d], component); }
mark_dofs(dofs, ess_vdofs);
}
}
}
for (int i = 0; i < es_faces.Size(); i++)
{
if (es_faces[i] < GetMesh()->GetNFaces())
{
if (component < 0)
{
GetFaceVDofs(es_faces[i], vdofs);
mark_dofs(vdofs, ess_vdofs);
}
else
{
GetFaceDofs(es_faces[i], dofs);
for (int d = 0; d < dofs.Size(); d++)
{ dofs[d] = DofToVDof(dofs[d], component); }
mark_dofs(dofs, ess_vdofs);
}
}
}
}
}
void FiniteElementSpace::GetEssentialVDofs(EntitySets::EntityType type,
const string & set_name,
Array<int> &ess_vdofs,
int component) const
{
MFEM_VERIFY(mesh->ent_sets != NULL, "Mesh object contains no "
"entity set information");
GetEssentialVDofs(type, mesh->ent_sets->GetSetIndex(type, set_name),
ess_vdofs, component);
}
void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
Array<int> &ess_tdof_list,
int component)
@@ -728,36 +579,6 @@ void FiniteElementSpace::GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
MarkerToList(ess_tdofs, ess_tdof_list);
}
void FiniteElementSpace::GetEssentialTrueDofs(EntitySets::EntityType type,
int set_index,
Array<int> &ess_tdof_list,
int component)
{
Array<int> ess_vdofs, ess_tdofs;
GetEssentialVDofs(type, set_index, ess_vdofs, component);
const SparseMatrix *R = GetConformingRestriction();
if (!R)
{
ess_tdofs.MakeRef(ess_vdofs);
}
else
{
R->BooleanMult(ess_vdofs, ess_tdofs);
}
MarkerToList(ess_tdofs, ess_tdof_list);
}
void FiniteElementSpace::GetEssentialTrueDofs(EntitySets::EntityType type,
const string & set_name,
Array<int> &ess_tdof_list,
int component)
{
MFEM_VERIFY(mesh->ent_sets != NULL, "Mesh object contains no "
"entity set information");
GetEssentialTrueDofs(type, mesh->ent_sets->GetSetIndex(type, set_name),
ess_tdof_list, component);
}
void FiniteElementSpace::GetBoundaryTrueDofs(Array<int> &boundary_dofs,
int component)
{
-26
View File
@@ -778,19 +778,6 @@ public:
Array<int> &ess_vdofs,
int component = -1) const;
/** Mark degrees of freedom associated with the entity set with the
specified entity type and set index. */
virtual void GetEssentialVDofs(EntitySets::EntityType type, int set_index,
Array<int> &ess_vdofs,
int component = -1) const;
/** Mark degrees of freedom associated with the entity set with the
specified entity type and set index. */
virtual void GetEssentialVDofs(EntitySets::EntityType type,
const std::string & set_name,
Array<int> &ess_vdofs,
int component = -1) const;
/** @brief Get a list of essential true dofs, ess_tdof_list, corresponding to the
boundary attributes marked in the array bdr_attr_is_ess.
For spaces with 'vdim' > 1, the 'component' parameter can be used
@@ -799,19 +786,6 @@ public:
Array<int> &ess_tdof_list,
int component = -1);
/** Get a list of essential true dofs, ess_tdof_list, corresponding to the
entity set specified by the given entity type and set index. */
virtual void GetEssentialTrueDofs(EntitySets::EntityType type, int set_index,
Array<int> &ess_tdof_list,
int component = -1);
/** Get a list of essential true dofs, ess_tdof_list, corresponding to the
entity set specified by the given entity type and set name. */
virtual void GetEssentialTrueDofs(EntitySets::EntityType type,
const std::string & set_name,
Array<int> &ess_tdof_list,
int component = -1);
/** @brief Get a list of all boundary true dofs, @a boundary_dofs. For spaces
with 'vdim' > 1, the 'component' parameter can be used to restricts the
marked tDOFs to the specified component. Equivalent to
+1 -2
View File
@@ -112,7 +112,7 @@ FmsFieldToGridFunction(FmsMesh fms_mesh, FmsField f, Mesh *mesh,
// NOTE: transplanted from the FmsMeshToMesh function
// We should do this work once and save it.
//--------------------------------------------------
FmsInt dim, n_vert, n_elem, space_dim;
FmsInt dim, n_elem, space_dim;
// Find the first component that has coordinates - that will be the new mfem
// mesh.
@@ -144,7 +144,6 @@ FmsFieldToGridFunction(FmsMesh fms_mesh, FmsField f, Mesh *mesh,
n_ents[et] += num_ents;
}
}
n_vert = n_ents[FMS_VERTEX];
//--------------------------------------------------
// Interrogate the field.
+1 -1
View File
@@ -1675,7 +1675,7 @@ IntegrationRule *IntegrationRules::PyramidIntegrationRule(int Order)
ipp.x = ipc.x * (1.0 - ipc.z);
ipp.y = ipc.y * (1.0 - ipc.z);
ipp.z = ipc.z;
ipp.weight = ipc.weight / 3.0;
ipp.weight = 4.0 * ipc.weight * pow(1.0 - ipc.z, 2);
}
return PyramidIntRules[Order];
}
+5 -54
View File
@@ -23,8 +23,6 @@
#include <limits>
#include <list>
using namespace std;
namespace mfem
{
@@ -1012,36 +1010,9 @@ void ParFiniteElementSpace::GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
{
FiniteElementSpace::GetEssentialVDofs(bdr_attr_is_ess, ess_dofs, component);
if (Conforming())
{
// Make sure that processors without boundary elements mark
// their boundary dofs (if they have any).
Synchronize(ess_dofs);
}
}
void ParFiniteElementSpace::GetEssentialVDofs(EntitySets::EntityType type,
int set_index,
Array<int> &ess_dofs,
int component) const
{
FiniteElementSpace::GetEssentialVDofs(type, set_index, ess_dofs, component);
if (Conforming())
{
// Make sure that processors without boundary elements mark
// their boundary dofs (if they have any).
Synchronize(ess_dofs);
}
}
void ParFiniteElementSpace::GetEssentialVDofs(EntitySets::EntityType type,
const string & set_name,
Array<int> &ess_vdofs,
int component) const
{
GetEssentialVDofs(type, pmesh->ent_sets->GetSetIndex(type, set_name),
ess_vdofs, component);
// Make sure that processors without boundary elements mark
// their boundary dofs (if they have any).
Synchronize(ess_dofs);
}
void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
@@ -1067,33 +1038,13 @@ void ParFiniteElementSpace::GetEssentialTrueDofs(const Array<int>
{
if (bool(ted[i]) != bool(true_ess_dofs2[i])) { counter++; }
}
MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter);
MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter
<< ", rank = " << MyRank);
#endif
MarkerToList(true_ess_dofs, ess_tdof_list);
}
void ParFiniteElementSpace::GetEssentialTrueDofs(EntitySets::EntityType type,
int set_index,
Array<int> &ess_tdof_list,
int component)
{
Array<int> ess_dofs, true_ess_dofs;
GetEssentialVDofs(type, set_index, ess_dofs, component);
GetRestrictionMatrix()->BooleanMult(ess_dofs, true_ess_dofs);
MarkerToList(true_ess_dofs, ess_tdof_list);
}
void ParFiniteElementSpace::GetEssentialTrueDofs(EntitySets::EntityType type,
const string & set_name,
Array<int> &ess_tdof_list,
int component)
{
GetEssentialTrueDofs(type, pmesh->ent_sets->GetSetIndex(type, set_name),
ess_tdof_list, component);
}
int ParFiniteElementSpace::GetLocalTDofNumber(int ldof) const
{
if (Nonconforming())
-26
View File
@@ -355,38 +355,12 @@ public:
Array<int> &ess_dofs,
int component = -1) const;
/** Mark degrees of freedom associated with the entity set with the
specified entity type and set index. */
virtual void GetEssentialVDofs(EntitySets::EntityType type, int set_index,
Array<int> &ess_vdofs,
int component = -1) const;
/** Mark degrees of freedom associated with the entity set with the
specified entity type and set index. */
virtual void GetEssentialVDofs(EntitySets::EntityType type,
const std::string & set_name,
Array<int> &ess_vdofs,
int component = -1) const;
/** Get a list of essential true dofs, ess_tdof_list, corresponding to the
boundary attributes marked in the array bdr_attr_is_ess. */
virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
Array<int> &ess_tdof_list,
int component = -1);
/** Get a list of essential true dofs, ess_tdof_list, corresponding to the
entity set specified by the given entity type and set index. */
virtual void GetEssentialTrueDofs(EntitySets::EntityType type, int set_index,
Array<int> &ess_tdof_list,
int component = -1);
/** Get a list of essential true dofs, ess_tdof_list, corresponding to the
entity set specified by the given entity type and set name. */
virtual void GetEssentialTrueDofs(EntitySets::EntityType type,
const std::string & set_name,
Array<int> &ess_tdof_list,
int component = -1);
/** If the given ldof is owned by the current processor, return its local
tdof number, otherwise return -1 */
int GetLocalTDofNumber(int ldof) const;
+4 -4
View File
@@ -65,13 +65,13 @@ static void Values2D(const int NE,
MFEM_SHARED double sm0[NBZ][MDQ*MDQ];
MFEM_SHARED double sm1[NBZ][MDQ*MDQ];
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
ConstDeviceMatrix B(sB, D1D,Q1D);
DeviceMatrix DD(sm0[tidz], MD1, MD1);
DeviceMatrix DQ(sm1[tidz], MD1, MQ1);
DeviceMatrix QQ(sm0[tidz], MQ1, MQ1);
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
for (int c = 0; c < VDIM; c++)
{
kernels::internal::LoadX(e,D1D,c,x,DD);
@@ -126,14 +126,14 @@ static void Values3D(const int NE,
MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
ConstDeviceMatrix B(sB, D1D,Q1D);
DeviceCube DDD(sm0, MD1,MD1,MD1);
DeviceCube DDQ(sm1, MD1,MD1,MQ1);
DeviceCube DQQ(sm0, MD1,MQ1,MQ1);
DeviceCube QQQ(sm1, MQ1,MQ1,MQ1);
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
for (int c = 0; c < VDIM; c++)
{
kernels::internal::LoadX(e,D1D,c,x,DDD);
+275 -267
View File
File diff suppressed because it is too large Load Diff
+30 -27
View File
@@ -1133,7 +1133,7 @@ protected:
// Evaluation of the discrete target specification on different meshes.
// Owned.
AdaptivityEvaluator *adapt_eval;
AdaptivityEvaluator *adapt_lim_eval;
void SetDiscreteTargetBase(const GridFunction &tspec_);
void SetTspecAtIndex(int idx, const GridFunction &tspec_);
@@ -1156,7 +1156,7 @@ public:
#endif
amr_el(-1), lim_min_size(-0.1),
good_tspec(false), good_tspec_grad(false), good_tspec_hess(false),
adapt_eval(NULL) { }
adapt_lim_eval(NULL) { }
virtual ~DiscreteAdaptTC();
@@ -1232,8 +1232,8 @@ public:
void SetAdaptivityEvaluator(AdaptivityEvaluator *ae)
{
if (adapt_eval) { delete adapt_eval; }
adapt_eval = ae;
if (adapt_lim_eval) { delete adapt_lim_eval; }
adapt_lim_eval = ae;
}
const Vector &GetTspecPert1H() { return tspec_pert1h; }
@@ -1315,15 +1315,15 @@ protected:
int integ_order;
// Weight Coefficient multiplying the quality metric term.
Coefficient *coeff1; // not owned, if NULL -> coeff1 is 1.
Coefficient *metric_coeff; // not owned, if NULL -> metric_coeff is 1.
// Normalization factor for the metric term.
double metric_normal;
// Nodes and weight Coefficient used for "limiting" the TMOP_Integrator.
// These are both NULL when there is no limiting.
// The class doesn't own nodes0 and coeff0.
const GridFunction *nodes0;
Coefficient *coeff0;
// The class doesn't own lim_nodes0 and lim_coeff.
const GridFunction *lim_nodes0;
Coefficient *lim_coeff;
// Limiting reference distance. Not owned.
const GridFunction *lim_dist;
// Limiting function. Owned.
@@ -1332,20 +1332,21 @@ protected:
double lim_normal;
// Adaptive limiting.
const GridFunction *zeta_0; // Not owned.
const GridFunction *adapt_lim_gf0; // Not owned.
#ifdef MFEM_USE_MPI
const ParGridFunction *pzeta_0;
const ParGridFunction *adapt_lim_pgf0;
#endif
GridFunction *zeta; // Owned. Updated by adapt_eval.
Coefficient *coeff_zeta; // Not owned.
AdaptivityEvaluator *adapt_eval; // Not owned.
GridFunction *adapt_lim_gf; // Owned. Updated by adapt_lim_eval.
Coefficient *adapt_lim_coeff; // Not owned.
AdaptivityEvaluator *adapt_lim_eval; // Not owned.
// Surface fitting.
GridFunction *sigma, *sigma_bar; // Owned. Updated by sigma_eval.
const Array<bool> *sigma_marker; // Not owned.
Coefficient *coeff_sigma; // Not owned.
AdaptivityEvaluator *sigma_eval; // Not owned.
double sigma_normal;
GridFunction *surf_fit_gf,
*surf_fit_gf_bar; // Owned, Updated by surf_fit_eval.
const Array<bool> *surf_fit_marker; // Not owned.
Coefficient *surf_fit_coeff; // Not owned.
AdaptivityEvaluator *surf_fit_eval; // Not owned.
double surf_fit_normal;
DiscreteAdaptTC *discr_tc;
@@ -1416,7 +1417,7 @@ protected:
void ComputeNormalizationEnergies(const GridFunction &x,
double &metric_energy, double &lim_energy,
double &sigma_energy);
double &surf_fit_gf_energy);
void AssembleElementVectorExact(const FiniteElement &el,
ElementTransformation &T,
@@ -1471,7 +1472,7 @@ protected:
void DisableLimiting()
{
nodes0 = NULL; coeff0 = NULL; lim_dist = NULL;
lim_nodes0 = NULL; lim_coeff = NULL; lim_dist = NULL;
delete lim_func; lim_func = NULL;
}
@@ -1531,12 +1532,14 @@ public:
TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc,
TMOP_QualityMetric *hm)
: h_metric(hm), metric(m), targetC(tc), IntegRules(NULL),
integ_order(-1), coeff1(NULL), metric_normal(1.0),
nodes0(NULL), coeff0(NULL),
integ_order(-1), metric_coeff(NULL), metric_normal(1.0),
lim_nodes0(NULL), lim_coeff(NULL),
lim_dist(NULL), lim_func(NULL), lim_normal(1.0),
zeta_0(NULL), zeta(NULL), coeff_zeta(NULL), adapt_eval(NULL),
sigma(NULL), sigma_bar(NULL), sigma_marker(NULL), coeff_sigma(NULL),
sigma_eval(NULL), sigma_normal(1.0),
adapt_lim_gf0(NULL), adapt_lim_gf(NULL), adapt_lim_coeff(NULL),
adapt_lim_eval(NULL),
surf_fit_gf(NULL), surf_fit_gf_bar(NULL), surf_fit_marker(NULL),
surf_fit_coeff(NULL),
surf_fit_eval(NULL), surf_fit_normal(1.0),
discr_tc(dynamic_cast<DiscreteAdaptTC *>(tc)),
fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false)
{ PA.enabled = false; }
@@ -1564,7 +1567,7 @@ public:
Note that the Coefficient is evaluated in the physical configuration and
not in the target configuration which may be undefined. */
void SetCoefficient(Coefficient &w1) { coeff1 = &w1; }
void SetCoefficient(Coefficient &w1) { metric_coeff = &w1; }
/** @brief Limiting of the mesh displacements (general version).
@@ -1631,7 +1634,7 @@ public:
void GetSurfaceFittingErrors(double &err_avg, double &err_max);
/// Update the original/reference nodes used for limiting.
void SetLimitingNodes(const GridFunction &n0) { nodes0 = &n0; }
void SetLimitingNodes(const GridFunction &n0) { lim_nodes0 = &n0; }
/** @brief Computes the integral of W(Jacobian(Trt)) over a target zone.
@param[in] el Type of FiniteElement.
+21 -21
View File
@@ -38,13 +38,13 @@ void TMOP_Integrator::AssembleGradPA(const Vector &xe,
if (PA.dim == 2)
{
AssembleGradPA_2D(xe);
if (coeff0) { AssembleGradPA_C0_2D(xe); }
if (lim_coeff) { AssembleGradPA_C0_2D(xe); }
}
if (PA.dim == 3)
{
AssembleGradPA_3D(xe);
if (coeff0) { AssembleGradPA_C0_3D(xe); }
if (lim_coeff) { AssembleGradPA_C0_3D(xe); }
}
}
@@ -53,8 +53,8 @@ void TMOP_Integrator::AssemblePA_Limiting()
const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
Device::GetDeviceMemoryType() : pa_mt;
// Return immediately if limiting is not enabled
if (coeff0 == nullptr) { return; }
MFEM_VERIFY(nodes0, "internal error");
if (lim_coeff == nullptr) { return; }
MFEM_VERIFY(lim_nodes0, "internal error");
MFEM_VERIFY(PA.enabled, "AssemblePA_Limiting but PA is not enabled!");
MFEM_VERIFY(lim_func, "No TMOP_LimiterFunction specification!")
@@ -68,14 +68,14 @@ void TMOP_Integrator::AssemblePA_Limiting()
const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
// H0 for coeff0, (dim x dim) Q-vector
// H0 for lim_coeff, (dim x dim) Q-vector
PA.H0.UseDevice(true);
PA.H0.SetSize(PA.dim * PA.dim * PA.nq * NE, mt);
// coeff0 -> PA.C0 (Q-vector)
// lim_coeff -> PA.C0 (Q-vector)
PA.C0.UseDevice(true);
if (ConstantCoefficient* cQ =
dynamic_cast<ConstantCoefficient*>(coeff0))
dynamic_cast<ConstantCoefficient*>(lim_coeff))
{
PA.C0.SetSize(1, Device::GetMemoryType());
PA.C0.HostWrite();
@@ -90,17 +90,17 @@ void TMOP_Integrator::AssemblePA_Limiting()
ElementTransformation& T = *fes->GetElementTransformation(e);
for (int q = 0; q < ir.GetNPoints(); ++q)
{
C0(q,e) = coeff0->Eval(T, ir.IntPoint(q));
C0(q,e) = lim_coeff->Eval(T, ir.IntPoint(q));
}
}
}
// nodes0 -> PA.X0 (E-vector)
MFEM_VERIFY(nodes0->FESpace() == fes, "");
// lim_nodes0 -> PA.X0 (E-vector)
MFEM_VERIFY(lim_nodes0->FESpace() == fes, "");
const Operator *n0_R = fes->GetElementRestriction(ordering);
PA.X0.SetSize(n0_R->Height(), Device::GetMemoryType());
PA.X0.UseDevice(true);
n0_R->Mult(*nodes0, PA.X0);
n0_R->Mult(*lim_nodes0, PA.X0);
// Limiting distances: lim_dist -> PA.LD (E-vector)
// TODO: remove the hack for the case lim_dist == NULL.
@@ -217,8 +217,8 @@ void TMOP_Integrator::AssemblePA(const FiniteElementSpace &fes)
PA.Jtr_needs_update = true;
PA.Jtr_debug_grad = false;
// Limiting: coeff0 -> PA.C0, nodes0 -> PA.X0, lim_dist -> PA.LD, PA.H0
if (coeff0) { AssemblePA_Limiting(); }
// Limiting: lim_coeff -> PA.C0, lim_nodes0 -> PA.X0, lim_dist -> PA.LD, PA.H0
if (lim_coeff) { AssemblePA_Limiting(); }
}
void TMOP_Integrator::AssembleGradDiagonalPA(Vector &de) const
@@ -236,13 +236,13 @@ void TMOP_Integrator::AssembleGradDiagonalPA(Vector &de) const
if (PA.dim == 2)
{
AssembleDiagonalPA_2D(de);
if (coeff0) { AssembleDiagonalPA_C0_2D(de); }
if (lim_coeff) { AssembleDiagonalPA_C0_2D(de); }
}
if (PA.dim == 3)
{
AssembleDiagonalPA_3D(de);
if (coeff0) { AssembleDiagonalPA_C0_3D(de); }
if (lim_coeff) { AssembleDiagonalPA_C0_3D(de); }
}
}
@@ -258,13 +258,13 @@ void TMOP_Integrator::AddMultPA(const Vector &xe, Vector &ye) const
if (PA.dim == 2)
{
AddMultPA_2D(xe,ye);
if (coeff0) { AddMultPA_C0_2D(xe,ye); }
if (lim_coeff) { AddMultPA_C0_2D(xe,ye); }
}
if (PA.dim == 3)
{
AddMultPA_3D(xe,ye);
if (coeff0) { AddMultPA_C0_3D(xe,ye); }
if (lim_coeff) { AddMultPA_C0_3D(xe,ye); }
}
}
@@ -283,13 +283,13 @@ void TMOP_Integrator::AddMultGradPA(const Vector &re, Vector &ce) const
if (PA.dim == 2)
{
AddMultGradPA_2D(re,ce);
if (coeff0) { AddMultGradPA_C0_2D(re,ce); }
if (lim_coeff) { AddMultGradPA_C0_2D(re,ce); }
}
if (PA.dim == 3)
{
AddMultGradPA_3D(re,ce);
if (coeff0) { AddMultGradPA_C0_3D(re,ce); }
if (lim_coeff) { AddMultGradPA_C0_3D(re,ce); }
}
}
@@ -307,13 +307,13 @@ double TMOP_Integrator::GetLocalStateEnergyPA(const Vector &xe) const
if (PA.dim == 2)
{
energy = GetLocalStateEnergyPA_2D(xe);
if (coeff0) { energy += GetLocalStateEnergyPA_C0_2D(xe); }
if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_2D(xe); }
}
if (PA.dim == 3)
{
energy = GetLocalStateEnergyPA_3D(xe);
if (coeff0) { energy += GetLocalStateEnergyPA_C0_3D(xe); }
if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_3D(xe); }
}
return energy;
+3 -3
View File
@@ -181,9 +181,9 @@ template<int T_D1D, int T_Q1D, int T_MAX> return_t kernel(__VA_ARGS__)
if (K##kernel.Find(id)) { return K##kernel.At(id)(__VA_ARGS__,0,0); }\
else {\
constexpr int T_MAX = 4;\
const int D1D = (id>>4)&0xF, Q1D = id&0xF;\
MFEM_VERIFY(D1D <= MAX_D1D && Q1D <= MAX_Q1D, "Max size error!");\
return kernel<0,0,T_MAX>(__VA_ARGS__,D1D,Q1D); }
const int d1d = (id>>4)&0xF, q1d = id&0xF;\
MFEM_VERIFY(d1d <= MAX_D1D && q1d <= MAX_Q1D, "Max size error!");\
return kernel<0,0,T_MAX>(__VA_ARGS__,d1d,q1d); }
} // namespace kernels
+2 -1
View File
@@ -55,6 +55,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, DatcSize,
MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
ConstDeviceMatrix B(sB, D1D, Q1D);
DeviceCube DDD(sm0, MD1,MD1,MD1);
DeviceCube DDQ(sm1, MD1,MD1,MQ1);
@@ -88,7 +90,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, DatcSize,
}
min = min_size[0];
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
kernels::internal::EvalX(D1D,Q1D,B,DDD,DDQ);
kernels::internal::EvalY(D1D,Q1D,B,DDQ,DQQ);
kernels::internal::EvalZ(D1D,Q1D,B,DQQ,QQQ);
+3 -3
View File
@@ -97,9 +97,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, AssembleDiagonalPA_Kernel_2D,
const double *Jtr = &J(0,0,qx,qy,e);
// Jrt = Jtr^{-1}
double j[4];
ConstDeviceMatrix Jrt(j,2,2);
kernels::CalcInverse<2>(Jtr, j);
double jrt_data[4];
ConstDeviceMatrix Jrt(jrt_data,2,2);
kernels::CalcInverse<2>(Jtr, jrt_data);
const double gg = G(qy,dy) * G(qy,dy);
const double gb = G(qy,dy) * B(qy,dy);
+4 -4
View File
@@ -68,8 +68,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultGradPA_Kernel_C0_2D,
double Xh[2];
kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ,Xh);
double B[4];
DeviceMatrix H(B,2,2);
double H_data[4];
DeviceMatrix H(H_data,2,2);
for (int i = 0; i < DIM; i++)
{
for (int j = 0; j < DIM; j++)
@@ -78,9 +78,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultGradPA_Kernel_C0_2D,
}
}
// p2 = B . Xh
// p2 = H . Xh
double p2[2];
kernels::Mult(2,2,B,Xh,p2);
kernels::Mult(2,2,H_data,Xh,p2);
kernels::internal::PushEval<MQ1,NBZ>(Q1D,qx,qy,p2,QQ);
}
}
+4 -4
View File
@@ -70,8 +70,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultGradPA_Kernel_C0_3D,
double Xh[3];
kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ,Xh);
double B[9];
DeviceMatrix H(B,3,3);
double H_data[9];
DeviceMatrix H(H_data,3,3);
for (int i = 0; i < DIM; i++)
{
for (int j = 0; j < DIM; j++)
@@ -80,9 +80,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultGradPA_Kernel_C0_3D,
}
}
// p2 = B . Xh
// p2 = H . Xh
double p2[3];
kernels::Mult(3,3,B,Xh,p2);
kernels::Mult(3,3,H_data,Xh,p2);
kernels::internal::PushEval<MQ1>(Q1D,qx,qy,qz,p2,QQQ);
}
}
+1 -2
View File
@@ -55,6 +55,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
MFEM_SHARED double sBLD[MQ1*MD1];
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
@@ -66,8 +67,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
kernels::internal::LoadX(e,D1D,LD,DDD);
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);
kernels::internal::EvalZ(D1D,Q1D,BLD,DQQ,QQQ);
+1 -1
View File
@@ -62,6 +62,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
MFEM_SHARED double B[MQ1*MD1];
MFEM_SHARED double sBLD[MQ1*MD1];
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
@@ -86,7 +87,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
kernels::internal::LoadX<MD1>(e,D1D,X1,DDD1);
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);
+1 -1
View File
@@ -63,6 +63,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
MFEM_SHARED double B[MQ1*MD1];
MFEM_SHARED double sBLD[MQ1*MD1];
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
@@ -87,7 +88,6 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
kernels::internal::LoadX<MD1>(e,D1D,X1,DDD1);
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);
+10 -10
View File
@@ -837,26 +837,26 @@ void TMOPHRSolver::ParUpdate()
}
#endif
void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh_, NonlinearForm *nlf_)
{
const FiniteElementSpace &fes = *mesh->GetNodalFESpace();
const FiniteElementSpace &fes = *mesh_->GetNodalFESpace();
// Update Nonlinear form and Set Essential BC
nlf->Update();
nlf_->Update();
const int dim = fes.GetFE(0)->GetDim();
if (move_bnd == false)
{
Array<int> ess_bdr(mesh->bdr_attributes.Max());
Array<int> ess_bdr(mesh_->bdr_attributes.Max());
ess_bdr = 1;
nlf->SetEssentialBC(ess_bdr);
nlf_->SetEssentialBC(ess_bdr);
}
else
{
const int nd = fes.GetBE(0)->GetDof();
int n = 0;
for (int i = 0; i < mesh->GetNBE(); i++)
for (int i = 0; i < mesh_->GetNBE(); i++)
{
const int attr = mesh->GetBdrElement(i)->GetAttribute();
const int attr = mesh_->GetBdrElement(i)->GetAttribute();
MFEM_VERIFY(!(dim == 2 && attr == 3),
"Boundary attribute 3 must be used only for 3D meshes. "
"Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
@@ -866,9 +866,9 @@ void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
}
Array<int> ess_vdofs(n), vdofs;
n = 0;
for (int i = 0; i < mesh->GetNBE(); i++)
for (int i = 0; i < mesh_->GetNBE(); i++)
{
const int attr = mesh->GetBdrElement(i)->GetAttribute();
const int attr = mesh_->GetBdrElement(i)->GetAttribute();
fes.GetBdrElementVDofs(i, vdofs);
if (attr == 1) // Fix x components.
{
@@ -891,7 +891,7 @@ void TMOPHRSolver::UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
{ ess_vdofs[n++] = vdofs[j]; }
}
}
nlf->SetEssentialVDofs(ess_vdofs);
nlf_->SetEssentialVDofs(ess_vdofs);
}
}
+16 -10
View File
@@ -444,15 +444,19 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
if (untangling == false && min_detT_out < 0.0)
{
// No untangling, and detJ got negative -- no good.
if (print_level >= 0)
{ mfem::out << "Scale = " << scale << " Neg det(J) found.\n"; }
if (print_options.iterations)
{
mfem::out << "Scale = " << scale << " Neg det(J) found.\n";
}
scale *= detJ_factor; continue;
}
if (untangling == true && min_detT_out < *min_det_ptr)
{
// Untangling, and detJ got even more negative -- no good.
if (print_level >= 0)
{ mfem::out << "Scale = " << scale << " Neg det(J) decreased.\n"; }
if (print_options.iterations)
{
mfem::out << "Scale = " << scale << " Neg det(J) decreased.\n";
}
scale *= detJ_factor; continue;
}
@@ -477,7 +481,7 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
if (energy_out > energy_in + 0.2*fabs(energy_in) ||
std::isnan(energy_out) != 0)
{
if (print_level >= 0)
if (print_options.iterations)
{
mfem::out << "Scale = " << scale << " Increasing energy: "
<< energy_in << " --> " << energy_out << '\n';
@@ -492,7 +496,7 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
if (norm_out > 1.2*norm_in)
{
if (print_level >= 0)
if (print_options.iterations)
{
mfem::out << "Scale = " << scale << " Norm increased: "
<< norm_in << " --> " << norm_out << '\n';
@@ -508,13 +512,15 @@ double TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
if (min_detT_out > 0.0)
{
*min_det_ptr = 0.0;
if (print_level >= 0)
if (print_options.summary || print_options.iterations ||
print_options.first_and_last)
{ mfem::out << "The mesh has been untangled at the used points!\n"; }
}
else { *min_det_ptr = untangle_factor * min_detT_out; }
}
if (print_level >= 0)
if (print_options.summary || print_options.iterations ||
print_options.first_and_last)
{
if (untangling)
{
@@ -568,9 +574,9 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
if (parallel)
{
#ifdef MFEM_USE_MPI
const ParNonlinearForm *nlf =
const ParNonlinearForm *pnlf =
dynamic_cast<const ParNonlinearForm *>(oper);
const ParFiniteElementSpace *pfesc = nlf->ParFESpace();
const ParFiniteElementSpace *pfesc = pnlf->ParFESpace();
Vector x_loc(pfesc->GetVSize());
pfesc->GetProlongationMatrix()->Mult(x, x_loc);
for (int i = 0; i < integs.Size(); i++)
+14 -14
View File
@@ -177,7 +177,7 @@ Device::~Device()
Get().device_mem_class = MemoryClass::HOST;
}
void Device::Configure(const std::string &device, const int dev)
void Device::Configure(const std::string &device, const int device_id)
{
// If a device was configured via the environment, skip the configuration,
// and avoid the 'singleton_device' to destroy the mm.
@@ -240,7 +240,7 @@ void Device::Configure(const std::string &device, const int dev)
#endif
// Perform setup.
Get().Setup(dev);
Get().Setup(device_id);
// Enable the device
Enable();
@@ -276,35 +276,35 @@ void Device::SetMemoryTypes(MemoryType h_mt, MemoryType d_mt)
// the call mm.Configure(...) in UpdateMemoryTypeAndClass()
}
void Device::Print(std::ostream &out)
void Device::Print(std::ostream &os)
{
out << "Device configuration: ";
os << "Device configuration: ";
bool add_comma = false;
for (int i = 0; i < Backend::NUM_BACKENDS; i++)
{
if (backends & internal::backend_list[i])
{
if (add_comma) { out << ','; }
if (add_comma) { os << ','; }
add_comma = true;
out << internal::backend_name[i];
os << internal::backend_name[i];
}
}
out << '\n';
os << '\n';
#ifdef MFEM_USE_CEED
if (Allows(Backend::CEED_MASK))
{
const char *ceed_backend;
CeedGetResource(internal::ceed, &ceed_backend);
out << "libCEED backend: " << ceed_backend << '\n';
os << "libCEED backend: " << ceed_backend << '\n';
}
#endif
out << "Memory configuration: "
<< MemoryTypeName[static_cast<int>(host_mem_type)];
os << "Memory configuration: "
<< MemoryTypeName[static_cast<int>(host_mem_type)];
if (Device::Allows(Backend::DEVICE_MASK))
{
out << ',' << MemoryTypeName[static_cast<int>(device_mem_type)];
os << ',' << MemoryTypeName[static_cast<int>(device_mem_type)];
}
out << std::endl;
os << std::endl;
}
void Device::UpdateMemoryTypeAndClass()
@@ -502,12 +502,12 @@ static void CeedDeviceSetup(const char* ceed_spec)
#endif
}
void Device::Setup(const int device)
void Device::Setup(const int device_id)
{
MFEM_VERIFY(ngpu == -1, "the mfem::Device is already configured!");
ngpu = 0;
dev = device;
dev = device_id;
#ifndef MFEM_USE_CUDA
MFEM_VERIFY(!Allows(Backend::CUDA_MASK),
"the CUDA backends require MFEM built with MFEM_USE_CUDA=YES");
+1 -1
View File
@@ -150,7 +150,7 @@ private:
static Device& Get() { return device_singleton; }
/// Setup switcher based on configuration settings
void Setup(const int dev = 0);
void Setup(const int device_id = 0);
void MarkBackend(Backend::Id b) { backends |= b; }
+4 -4
View File
@@ -1228,12 +1228,12 @@ Graph::reweight(uint k)
// Linearly order graph.
void
Graph::order(Functional* functional, uint iterations, uint window, uint period,
uint seed, Progress* progress)
Graph::order(Functional* functional_, uint iterations, uint window, uint period,
uint seed, Progress* progress_)
{
// Initialize graph.
this->functional = functional;
progress = this->progress = progress ? progress : new Progress;
this->functional = functional_;
progress_ = this->progress = progress_ ? progress_ : new Progress;
for (level = 0; (1u << level) < nodes(); level++);
place();
Float mincost = cost();
+5 -4
View File
@@ -849,13 +849,14 @@ inline void Memory<T>::New(int size, MemoryType mt)
}
template <typename T>
inline void Memory<T>::New(int size, MemoryType h_mt, MemoryType d_mt)
inline void Memory<T>::New(int size, MemoryType host_mt, MemoryType device_mt)
{
capacity = size;
const size_t bytes = size*sizeof(T);
this->h_mt = h_mt;
T *h_tmp = (h_mt == MemoryType::HOST) ? NewHOST(size) : nullptr;
h_ptr = (T*)MemoryManager::New_(h_tmp, bytes, h_mt, d_mt, VALID_HOST, flags);
this->h_mt = host_mt;
T *h_tmp = (host_mt == MemoryType::HOST) ? NewHOST(size) : nullptr;
h_ptr = (T*)MemoryManager::New_(h_tmp, bytes, host_mt, device_mt,
VALID_HOST, flags);
}
template <typename T>
-1
View File
@@ -69,7 +69,6 @@ void IntegerSet::Recreate(const int n, const int *p)
me.Sort();
// Remove duplicate entries
for (j = 0, i = 1; i < n; i++)
if (me[i] != me[j])
{
+2 -7
View File
@@ -36,7 +36,7 @@ public:
IntegerSet(const int n, const int *p) { Recreate(n, p); }
/// Return the size of the set.
int Size() const { return me.Size(); }
int Size() { return me.Size(); }
/// Return a reference to the sorted array of all the set entries.
operator Array<int>& () { return me; }
@@ -50,8 +50,6 @@ public:
/// Return 1 if the sets are equal and 0 otherwise.
int operator==(IntegerSet &s);
inline const int & operator[](int i) const { return me[i]; }
/** @brief Create an integer set from C-array 'p' of 'n' integers.
Overwrites any existing set data. */
void Recreate(const int n, const int *p);
@@ -66,7 +64,7 @@ private:
public:
/// Return the number of integer sets in the list.
int Size() const { return TheList.Size(); }
int Size() { return TheList.Size(); }
/// Return the value of the first element of the ith set.
int PickElementInSet(int i) { return TheList[i]->PickElement(); }
@@ -86,9 +84,6 @@ public:
/// Write the list of sets into table 't'.
void AsTable(Table &t);
inline const IntegerSet & operator[](int i) const { return *TheList[i]; }
inline IntegerSet & operator[](int i) { return *TheList[i]; }
~ListOfIntegerSets();
};
+7 -9
View File
@@ -61,7 +61,7 @@ inline void Sort3 (int &r, int &c, int &f)
}
}
int STable3D::Push (int r, int c, int f, int t)
int STable3D::Push (int r, int c, int f)
{
STable3DNode *node;
@@ -86,7 +86,6 @@ int STable3D::Push (int r, int c, int f, int t)
#endif
node->Column = c;
node->Floor = f;
node->Tier = t;
node->Number = NElem;
node->Prev = Rows[r];
Rows[r] = node;
@@ -110,9 +109,9 @@ int STable3D::operator() (int r, int c, int f) const
}
}
// MFEM_ABORT("(r,c,f) = (" << r << "," << c << "," << f << ")");
MFEM_ABORT("(r,c,f) = (" << r << "," << c << "," << f << ")");
return -1;
return 0;
}
int STable3D::Index (int r, int c, int f) const
@@ -153,13 +152,13 @@ int STable3D::Push4 (int r, int c, int f, int t)
switch (i)
{
case 0:
return Push (c,f,t,r);
return Push (c,f,t);
case 1:
return Push (r,f,t,c);
return Push (r,f,t);
case 2:
return Push (r,c,t,f);
return Push (r,c,t);
case 3:
return Push (r,c,f,t);
return Push (r,c,f);
}
return -1;
@@ -219,7 +218,6 @@ void STable3D::Print(std::ostream & out) const
out << row
<< ' ' << node_p->Column
<< ' ' << node_p->Floor
<< ' ' << node_p->Tier
<< ' ' << node_p->Number
<< endl;
node_p = node_p->Prev;
+3 -22
View File
@@ -15,8 +15,6 @@
#include "mem_alloc.hpp"
#include "../general/globals.hpp"
#include <iostream>
namespace mfem
{
@@ -24,7 +22,7 @@ class STable3DNode
{
public:
STable3DNode *Prev;
int Column, Floor, Tier, Number;
int Column, Floor, Number;
};
/** @brief Symmetric 3D Table stored as an array of rows each of which has a
@@ -49,7 +47,7 @@ public:
/** @brief Check to see if this entry is in the table and add it to the table
if it is not there. Returns the number assigned to the table entry. */
int Push (int r, int c, int f, int t = -1);
int Push (int r, int c, int f);
/// Return the number assigned to the table entry. Abort if it's not there.
int operator() (int r, int c, int f) const;
@@ -68,30 +66,13 @@ public:
not there. */
int operator() (int r, int c, int f, int t) const;
/// Return the number of rows added to the table.
int NumberOfRows() const { return Size; }
/// Return the number of elements added to the table.
int NumberOfElements() const { return NElem; }
int NumberOfElements() { return NElem; }
/// Print out all of the table elements.
void Print(std::ostream &out = mfem::out) const;
~STable3D ();
class RowIterator
{
private:
STable3DNode *n;
public:
RowIterator (const STable3D &t, int r) { n = t.Rows[r]; }
int operator!() { return (n != NULL); }
void operator++() { n = n->Prev; }
int Column() { return (n->Column); }
int Floor() { return (n->Floor); }
int Tier() { return (n->Tier); }
int Index() { return (n->Number); }
};
};
}
+1 -1
View File
@@ -42,7 +42,7 @@ public:
/// Clear the elapsed time on the stopwatch and restart it if it's running.
void Clear();
/// Clear the elapsed time and start the stopwatch.
/// Start the stopwatch. The elapsed time is @b not cleared.
void Start();
/// Stop the stopwatch.
+3 -3
View File
@@ -327,7 +327,7 @@ void EliminationSolver::Mult(const Vector& rhs, Vector& sol) const
krylov->SetMaxIter(max_iter);
krylov->SetRelTol(rel_tol);
krylov->SetAbsTol(abs_tol);
krylov->SetPrintLevel(print_level);
krylov->SetPrintLevel(print_options);
Vector rtilde(rhs.Size());
if (constraint_rhs.Size() > 0)
@@ -447,7 +447,7 @@ void PenaltyConstrainedSolver::Mult(const Vector& b, Vector& x) const
krylov->SetRelTol(rel_tol);
krylov->SetAbsTol(abs_tol);
krylov->SetMaxIter(max_iter);
krylov->SetPrintLevel(print_level);
krylov->SetPrintLevel(print_options);
krylov->Mult(penalized_rhs, x);
final_iter = krylov->GetNumIterations();
final_norm = krylov->GetFinalNorm();
@@ -574,7 +574,7 @@ void SchurConstrainedSolver::LagrangeSystemMult(const Vector& x,
gmres->SetRelTol(rel_tol);
gmres->SetAbsTol(abs_tol);
gmres->SetMaxIter(max_iter);
gmres->SetPrintLevel(print_level);
gmres->SetPrintLevel(print_options);
gmres->SetPreconditioner(
const_cast<BlockDiagonalPreconditioner&>(*block_pc));
+5 -1
View File
@@ -701,11 +701,15 @@ private:
public:
/// Create an operator which is a scalar multiple of A.
explicit ScaledOperator(const Operator *A, double a)
: Operator(A->Width(), A->Height()), A_(*A), a_(a) { }
: Operator(A->Height(), A->Width()), A_(*A), a_(a) { }
/// Operator application
virtual void Mult(const Vector &x, Vector &y) const
{ A_.Mult(x, y); y *= a_; }
/// Application of the transpose.
virtual void MultTranspose(const Vector &x, Vector &y) const
{ A_.MultTranspose(x, y); y *= a_; }
};
+9
View File
@@ -29,6 +29,10 @@ struct MFEM_AUTOSIMD_ALIGN__ AutoSIMD
scalar_t vec[size];
AutoSIMD() = default;
AutoSIMD(const AutoSIMD &) = default;
inline MFEM_ALWAYS_INLINE scalar_t &operator[](int i)
{
return vec[i];
@@ -117,6 +121,11 @@ struct MFEM_AUTOSIMD_ALIGN__ AutoSIMD
return r;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
{
return *this;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
{
AutoSIMD r;
+9 -1
View File
@@ -38,6 +38,10 @@ template <> struct AutoSIMD<double,2,16>
double vec[size];
};
AutoSIMD() = default;
AutoSIMD(const AutoSIMD &) = default;
inline MFEM_ALWAYS_INLINE double &operator[](int i)
{
return vec[i];
@@ -115,6 +119,11 @@ template <> struct AutoSIMD<double,2,16>
return r;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
{
return *this;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
{
AutoSIMD r;
@@ -251,4 +260,3 @@ AutoSIMD<double,2,16> operator/(const double &e,
#endif // __SSE2__
#endif // MFEM_SIMD_M128_HPP
+9
View File
@@ -38,6 +38,10 @@ template <> struct AutoSIMD<double,4,32>
double vec[size];
};
AutoSIMD() = default;
AutoSIMD(const AutoSIMD &) = default;
inline MFEM_ALWAYS_INLINE double &operator[](int i)
{
return vec[i];
@@ -115,6 +119,11 @@ template <> struct AutoSIMD<double,4,32>
return r;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
{
return *this;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
{
AutoSIMD r;
+9
View File
@@ -39,6 +39,10 @@ template <> struct AutoSIMD<double,8,64>
double vec[size];
};
AutoSIMD() = default;
AutoSIMD(const AutoSIMD &) = default;
inline MFEM_ALWAYS_INLINE double &operator[](int i)
{
return vec[i];
@@ -120,6 +124,11 @@ template <> struct AutoSIMD<double,8,64>
return r;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
{
return *this;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
{
AutoSIMD r;
+9
View File
@@ -34,6 +34,10 @@ template <> struct AutoSIMD<double,4,32>
double vec[size];
};
AutoSIMD() = default;
AutoSIMD(const AutoSIMD &) = default;
inline __ATTRS_ai double &operator[](int i) { return vec[i]; }
inline __ATTRS_ai const double &operator[](int i) const { return vec[i]; }
@@ -105,6 +109,11 @@ template <> struct AutoSIMD<double,4,32>
return r;
}
inline __ATTRS_ai AutoSIMD operator+() const
{
return *this;
}
inline __ATTRS_ai AutoSIMD operator+(const AutoSIMD &v) const
{
AutoSIMD r;
+9
View File
@@ -33,6 +33,10 @@ template <> struct MFEM_AUTOSIMD_ALIGN_SVE AutoSIMD<double,8,64>
double vec[size];
AutoSIMD() = default;
AutoSIMD(const AutoSIMD &) = default;
inline MFEM_ALWAYS_INLINE double &operator[](int i)
{
return vec[i];
@@ -123,6 +127,11 @@ template <> struct MFEM_AUTOSIMD_ALIGN_SVE AutoSIMD<double,8,64>
return r;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
{
return *this;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
{
AutoSIMD r;
+17
View File
@@ -17,6 +17,10 @@
#include "../../config/tconfig.hpp"
#include <altivec.h>
#ifdef __GNUC__
#undef bool
#endif
namespace mfem
{
@@ -34,6 +38,10 @@ template <> struct AutoSIMD<double,2,16>
double vec[size];
};
AutoSIMD() = default;
AutoSIMD(const AutoSIMD &) = default;
inline MFEM_ALWAYS_INLINE double &operator[](int i)
{
return vec[i];
@@ -107,10 +115,19 @@ template <> struct AutoSIMD<double,2,16>
inline MFEM_ALWAYS_INLINE AutoSIMD operator-() const
{
AutoSIMD r;
#ifndef __GNUC__
r.vd = vec_neg(vd);
#else
r.vd = vec_splats(0.0) - vd;
#endif
return r;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+() const
{
return *this;
}
inline MFEM_ALWAYS_INLINE AutoSIMD operator+(const AutoSIMD &v) const
{
AutoSIMD r;
+370 -181
View File
File diff suppressed because it is too large Load Diff
+158 -9
View File
@@ -56,7 +56,7 @@ public:
{
}
/** @brief This method is invoked by ItertiveSolver::SetMonitor, informing
/** @brief This method is invoked by IterativeSolver::SetMonitor, informing
the monitor which IterativeSolver is using it. */
void SetIterativeSolver(const IterativeSolver &solver)
{ iter_solver = &solver; }
@@ -65,10 +65,55 @@ public:
/// Abstract base class for iterative solver
class IterativeSolver : public Solver
{
public:
/** @brief Settings for the output behavior of the IterativeSolver.
By default, all output is suppressed. The construction of the desired
print level can be achieved through a builder pattern, for example
PrintLevel().Errors().Warnings()
constructs the print level with only errors and warnings enabled.
*/
struct PrintLevel
{
/** @brief If a fatal problem has been detected the failure will be
reported to @ref mfem::err. */
bool errors = false;
/** @brief If a non-fatal problem has been detected some context-specific
information will be reported to @ref mfem::out */
bool warnings = false;
/** @brief Detailed information about each iteration will be reported to
@ref mfem::out */
bool iterations = false;
/** @brief A summary of the solver process will be reported after the last
iteration to @ref mfem::out */
bool summary = false;
/** @brief Information about the first and last iteration will be printed
to @ref mfem::out */
bool first_and_last = false;
/// Initializes the print level to suppress
PrintLevel() = default;
/** @name Builder
These methods are utilized to construct PrintLevel objects through a
builder approach by chaining the function calls in this group. */
///@{
PrintLevel &None() { *this = PrintLevel(); return *this; }
PrintLevel &Warnings() { warnings=true; return *this; }
PrintLevel &Errors() { errors=true; return *this; }
PrintLevel &Iterations() { iterations=true; return *this; }
PrintLevel &FirstAndLast() { first_and_last=true; return *this; }
PrintLevel &Summary() { summary=true; return *this; }
PrintLevel &All() { return Warnings().Errors().Iterations().FirstAndLast().Summary(); }
///@}
};
#ifdef MFEM_USE_MPI
private:
int dot_prod_type; // 0 - local, 1 - global over 'comm'
MPI_Comm comm;
MPI_Comm comm = MPI_COMM_NULL;
#endif
protected:
@@ -76,13 +121,52 @@ protected:
Solver *prec;
IterativeSolverMonitor *monitor = nullptr;
int max_iter, print_level;
double rel_tol, abs_tol;
/// @name Reporting (protected attributes and member functions)
///@{
// stats
mutable int final_iter, converged;
/** @brief (DEPRECATED) Legacy print level definition, which is left for
compatibility with custom iterative solvers.
@deprecated #print_options should be used instead. */
int print_level = -1;
/** @brief Output behavior for the iterative solver.
This primarily controls the output behavior of the iterative solvers
provided by this library. This member must be synchronized with
#print_level to ensure compatibility with custom iterative solvers. */
PrintLevel print_options;
/// Convert a legacy print level integer to a PrintLevel object
PrintLevel FromLegacyPrintLevel(int);
/// @brief Use some heuristics to guess a legacy print level corresponding to
/// the given PrintLevel.
static int GuessLegacyPrintLevel(PrintLevel);
///@}
/// @name Convergence (protected attributes)
///@{
/// Limit for the number of iterations the solver is allowed to do
int max_iter;
/// Relative tolerance.
double rel_tol;
/// Absolute tolerance.
double abs_tol;
///@}
/// @name Solver statistics (protected attributes)
///@{
mutable int final_iter;
mutable bool converged;
mutable double final_norm;
///@}
double Dot(const Vector &x, const Vector &y) const;
double Norm(const Vector &x) const { return sqrt(Dot(x, x)); }
void Monitor(int it, double norm, const Vector& r, const Vector& x,
@@ -95,20 +179,80 @@ public:
IterativeSolver(MPI_Comm comm_);
#endif
/** @name Convergence
@brief Termination criteria for the iterative solvers.
@details While the convergence criterion is solver specific, most of the
provided iterative solvers use one of the following criteria
\f$ ||r||_X \leq tol_{rel}||r_0||_X \f$,
\f$ ||r||_X \leq tol_{abs} \f$,
\f$ ||r||_X \leq \max\{ tol_{abs}, tol_{rel} ||r_0||_X \} \f$,
where X denotes the space in which the norm is measured. The choice of
X depends on the specific iterative solver.
*/
///@{
void SetRelTol(double rtol) { rel_tol = rtol; }
void SetAbsTol(double atol) { abs_tol = atol; }
void SetMaxIter(int max_it) { max_iter = max_it; }
void SetPrintLevel(int print_lvl);
///@}
/** @name Reporting
These options control the internal reporting behavior into ::mfem::out
and ::mfem::err of the iterative solvers.
*/
///@{
/// @brief Legacy method to set the level of verbosity of the solver output.
/** This is the old way to control what information will be printed to
::mfem::out and ::mfem::err. The behavior for the print level for all
iterative solvers is:
- -1: Suppress all outputs.
- 0: Print information about all detected issues (e.g. no convergence).
- 1: Same as level 0, but with detailed information about each
iteration.
- 2: Print detected issues and a summary when the solver terminates.
- 3: Same as 2, but print also the first and last iterations.
- >3: Custom print options which are dependent on the specific solver.
In parallel, only rank 0 produces output.
@note It is recommended to use @ref SetPrintLevel(PrintLevel) instead.
@note Some derived classes, like KINSolver, redefine this method and use
their own set of print level constants. */
virtual void SetPrintLevel(int print_lvl);
/// @brief Set the level of verbosity of the solver output.
/** In parallel, only rank 0 produces outputs. Errors are output to
::mfem::err and all other information to ::mfem::out.
@note Not all subclasses of IterativeSolver support all possible options.
@note Some derived classes, like KINSolver, disable this method in favor
of SetPrintLevel(int).
@sa PrintLevel for possible options.
*/
virtual void SetPrintLevel(PrintLevel);
///@}
/// @name Solver statistics
///@{
int GetNumIterations() const { return final_iter; }
int GetConverged() const { return converged; }
bool GetConverged() const { return converged; }
double GetFinalNorm() const { return final_norm; }
///@}
/// This should be called before SetOperator
virtual void SetPreconditioner(Solver &pr);
/// Also calls SetOperator for the preconditioner if there is one
virtual void SetOperator(const Operator &op);
virtual void SetOperator(const Operator &op) override;
/// Set the iterative solver monitor
void SetMonitor(IterativeSolverMonitor &m)
@@ -163,6 +307,9 @@ public:
~OperatorJacobiSmoother() {}
/// Replace diagonal entries with their absolute values.
void SetPositiveDiagonal(bool pos_diag = true) { use_abs_diag = pos_diag; }
void Mult(const Vector &x, Vector &y) const;
void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
@@ -184,6 +331,8 @@ private:
const double damping;
const Array<int> *ess_tdof_list; // not owned; may be NULL
mutable Vector residual;
/// Uses absolute values of the diagonal entries.
bool use_abs_diag = false;
const Operator *oper; // not owned
+8 -5
View File
@@ -2382,7 +2382,7 @@ double SparseMatrix::GetJacobiScaling() const
}
void SparseMatrix::Jacobi(const Vector &b, const Vector &x0, Vector &x1,
double sc) const
double sc, bool use_abs_diag) const
{
MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
@@ -2403,7 +2403,8 @@ void SparseMatrix::Jacobi(const Vector &b, const Vector &x0, Vector &x1,
}
if (d >= 0 && A[d] != 0.0)
{
x1(i) = sc * (sum / A[d]) + (1.0 - sc) * x0(i);
const double diag = (use_abs_diag) ? fabs(A[d]) : A[d];
x1(i) = sc * (sum / diag) + (1.0 - sc) * x0(i);
}
else
{
@@ -2412,7 +2413,8 @@ void SparseMatrix::Jacobi(const Vector &b, const Vector &x0, Vector &x1,
}
}
void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
void SparseMatrix::DiagScale(const Vector &b, Vector &x,
double sc, bool use_abs_diag) const
{
MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
@@ -2438,11 +2440,12 @@ void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
}
if (Jp[j] == i)
{
if (!(std::abs(Ap[j]) > 0.0))
const double diag = (use_abs_diag) ? fabs(Ap[j]) : Ap[j];
if (diag == 0.0)
{
MFEM_ABORT_KERNEL("Zero diagonal in SparseMatrix::DiagScale");
}
xp[i] = sc * bp[i] / Ap[j];
xp[i] = sc * bp[i] / diag;
break;
}
}
+7 -3
View File
@@ -446,10 +446,14 @@ public:
/// Determine appropriate scaling for Jacobi iteration
double GetJacobiScaling() const;
/** One scaled Jacobi iteration for the system A x = b.
x1 = x0 + sc D^{-1} (b - A x0) where D is the diag of A. */
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, double sc) const;
x1 = x0 + sc D^{-1} (b - A x0) where D is the diag of A.
Absolute values of D are used when use_abs_diag = true. */
void Jacobi(const Vector &b, const Vector &x0, Vector &x1,
double sc, bool use_abs_diag = false) const;
void DiagScale(const Vector &b, Vector &x, double sc = 1.0) const;
/// x = sc b / A_ii. When use_abs_diag = true, |A_ii| is used.
void DiagScale(const Vector &b, Vector &x,
double sc = 1.0, bool use_abs_diag = false) const;
/** x1 = x0 + sc D^{-1} (b - A x0) where \f$ D_{ii} = \sum_j |A_{ij}| \f$. */
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
+2 -2
View File
@@ -65,7 +65,7 @@ void DSmoother::Mult(const Vector &x, Vector &y) const
{
if (!iterative_mode && type == 0 && iterations == 1)
{
oper->DiagScale(x, y, scale);
oper->DiagScale(x, y, scale, use_abs_diag);
return;
}
@@ -90,7 +90,7 @@ void DSmoother::Mult(const Vector &x, Vector &y) const
{
if (type == 0)
{
oper->Jacobi(x, *p, *r, scale);
oper->Jacobi(x, *p, *r, scale, use_abs_diag);
}
else if (type == 1)
{
+5
View File
@@ -58,6 +58,8 @@ protected:
int type; // 0, 1, 2 - scaled Jacobi, scaled l1-Jacobi, scaled lumped-Jacobi
double scale;
int iterations;
/// Uses abs values of the diagonal entries. Relevant only when type = 0.
bool use_abs_diag = false;
mutable Vector z;
@@ -69,6 +71,9 @@ public:
/// Create Jacobi smoother.
DSmoother(const SparseMatrix &a, int t = 0, double s = 1., int it = 1);
/// Replace diag entries with their abs values. Relevant only when type = 0.
void SetPositiveDiagonal(bool pos_diag = true) { use_abs_diag = pos_diag; }
/// Matrix vector multiplication with Jacobi smoother.
virtual void Mult(const Vector &x, Vector &y) const;
};
+5
View File
@@ -1976,6 +1976,11 @@ void KINSolver::SetMAA(int m_aa)
}
}
void KINSolver::SetPrintLevel(PrintLevel)
{
MFEM_ABORT("this method is not supported! Use SetPrintLevel(int) instead.");
}
// Compute the scaling vectors and solve nonlinear system
void KINSolver::Mult(const Vector&, Vector &x) const
{
+6
View File
@@ -833,6 +833,12 @@ public:
/** @note Only valid in combination with JFNK */
void SetLSMaxRestarts(int m) { maxlrs = m; }
/// Set the print level for the KINSetPrintLevel function.
virtual void SetPrintLevel(int print_lvl) { print_level = print_lvl; }
/// This method is not supported and will throw an error.
virtual void SetPrintLevel(PrintLevel);
/// Solve the nonlinear system \f$ F(x) = 0 \f$.
/** This method computes the x_scale and fx_scale vectors and calls the
other Mult(Vector&, Vector&, Vector&) const method. The x_scale vector
-1295
View File
File diff suppressed because it is too large Load Diff
-219
View File
@@ -1,219 +0,0 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#ifndef MFEM_ENTITY_SETS
#define MFEM_ENTITY_SETS
#include "../config/config.hpp"
#include "../general/table.hpp"
#include "../general/stable3d.hpp"
#include <limits>
#include <map>
#include <set>
#include <string>
#include <vector>
namespace mfem
{
class Mesh;
class NCMesh;
class NCEntitySets;
class EntitySets
{
friend class Mesh;
friend class NCMesh;
friend class NCEntitySets;
public:
enum EntityType {INVALID = -1, VERTEX = 0, EDGE = 1, FACE = 2, ELEMENT = 3};
static std::map<EntityType,std::string> EntityTypeNames;
EntitySets(Mesh & mesh);
EntitySets(const EntitySets & ent_sets);
EntitySets(Mesh & mesh, NCMesh &ncmesh);
virtual ~EntitySets();
static const std::string & GetTypeName(EntityType t);
bool SetExists(EntityType t, unsigned int s) const;
bool SetExists(EntityType t, const std::string & s) const;
void Load(std::istream &input);
void Print(std::ostream &output) const;
virtual void PrintSetInfo(std::ostream &output) const;
inline Mesh *GetMesh() const { return mesh_; }
unsigned int GetNumSets(EntityType t) const;
const std::string & GetSetName(EntityType t, unsigned int s) const;
unsigned int GetNumEntities(EntityType t, unsigned int s) const;
int GetSetIndex(EntityType t, const std::string & s) const;
unsigned int GetNumEntities(EntityType t, const std::string & s) const;
inline std::set<int> & operator()(EntityType t, unsigned int s)
{ return sets_[t][s]; }
inline const std::set<int> & operator()(EntityType t, unsigned int s) const
{ return sets_[t][s]; }
const Table * GetEdgeVertexTable() const { return edge_vertex_; }
const Table * GetFaceVertexTable() const { return face_vertex_; }
const Table * GetFaceEdgeTable() const { return face_edge_; }
// void Prune(int nelems);
protected:
void SetNumSets(EntityType t, unsigned int n)
{ sets_[t].resize(n); set_names_[t].resize(n); }
void SetSetName(EntityType t, int s, const std::string & name)
{ set_names_[t][s] = name; set_index_by_name_[t][name] = s; }
/// Make local copies of edge_vertex, face_vertex, and face_edge tables.
void CopyMeshTables();
/// Refine quadrilateral mesh.
virtual void QuadUniformRefinement();
/// Refine hexahedral mesh.
virtual void HexUniformRefinement();
/// Refine 2D mesh.
virtual void UniformRefinement2D();
/// Refine 3D mesh.
virtual void UniformRefinement3D();
private:
static void skip_comment_lines(std::istream &is, const char comment_char)
{
while (1)
{
is >> std::ws;
if (is.peek() != comment_char) { break; }
is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
}
}
// Check for, and remove, a trailing '\r'.
static void filter_dos(std::string &line)
{
if (!line.empty() && *line.rbegin() == '\r')
{ line.resize(line.size()-1); }
}
static std::map<EntityType,std::string> init_type_names();
void LoadEntitySets(std::istream &input, EntityType t,
const std::string & header);
void PrintEntitySets(std::ostream &output, EntityType t,
const std::string & header) const;
void PrintEdgeSets(std::ostream &output) const;
void PrintFaceSets(std::ostream &output) const;
void PrintEntitySetInfo(std::ostream & output, EntityType t,
const std::string & ent_name) const;
void CopyEntitySets(const EntitySets & ent_sets, EntityType t);
void BuildEntitySets(NCMesh &ncmesh, EntityType t);
protected:
Mesh * mesh_;
Table * edge_vertex_;
Table * face_vertex_;
Table * face_edge_;
int NumOfVertices_;
int NumOfEdges_;
int NumOfElements_;
/** The node/edge/face/element indices needed by the finite element
space to look up DoFs. */
std::vector<std::vector<std::set<int> > > sets_;
/// Names of each entity set
std::vector<std::vector<std::string> > set_names_;
/// Indices of each entity set indexed by set name
std::vector<std::map<std::string, int> > set_index_by_name_;
};
class NCEntitySets
{
friend class EntitySets;
public:
NCEntitySets(const EntitySets & ent_sets, NCMesh &ncmesh);
NCEntitySets(const NCEntitySets & ncent_sets);
bool SetExists(EntitySets::EntityType t, unsigned int s) const;
bool SetExists(EntitySets::EntityType t, const std::string & s) const;
unsigned int GetNumSets(EntitySets::EntityType t) const;
static int GetEntitySize(EntitySets::EntityType t);
const std::string & GetSetName(EntitySets::EntityType t, int s) const;
unsigned int GetNumEntities(EntitySets::EntityType t, int s) const;
void GetEntityIndex(EntitySets::EntityType t, int s,
int i, Array<int> & inds) const;
int GetSetIndex(EntitySets::EntityType t,
const std::string & s) const;
unsigned int GetNumEntities(EntitySets::EntityType t,
const std::string & s) const;
void GetEntityIndex(EntitySets::EntityType t,
const std::string & s, int i,
Array<int> & inds) const;
inline std::vector<int> & operator()(EntitySets::EntityType t, int s)
{ return sets_[t][s]; }
inline const std::vector<int> & operator()(EntitySets::EntityType t,
int s) const
{ return sets_[t][s]; }
inline int & operator()(EntitySets::EntityType t, int s, int i)
{ return sets_[t][s][i]; }
inline int operator()(EntitySets::EntityType t, int s, int i) const
{ return sets_[t][s][i]; }
private:
void CopyNCEntitySets(const NCEntitySets & ncent_sets,
EntitySets::EntityType t);
protected:
NCMesh * ncmesh_;
/// The nodes defining the node/edge/face/element sets
std::vector<std::vector<std::vector<int> > > sets_;
/// Names of each entity set
std::vector<std::vector<std::string> > set_names_;
/// Indices of each entity set indexed by set name
std::vector<std::map<std::string, int> > set_index_by_name_;
/// Number of indices per entity
static const int entity_size_[4];
};
} // namespace mfem
#endif // MFEM_ENTITY_SETS
+11 -70
View File
@@ -350,6 +350,7 @@ void Mesh::GetElementTransformation(int i, IsoparametricTransformation *ElTr)
ElTr->Attribute = GetAttribute(i);
ElTr->ElementNo = i;
ElTr->ElementType = ElementTransformation::ELEMENT;
ElTr->mesh = this;
ElTr->Reset();
if (Nodes == NULL)
{
@@ -382,6 +383,7 @@ void Mesh::GetElementTransformation(int i, const Vector &nodes,
ElTr->Attribute = GetAttribute(i);
ElTr->ElementNo = i;
ElTr->ElementType = ElementTransformation::ELEMENT;
ElTr->mesh = this;
DenseMatrix &pm = ElTr->GetPointMat();
ElTr->Reset();
nodes.HostRead();
@@ -437,6 +439,7 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
ElTr->Attribute = GetBdrAttribute(i);
ElTr->ElementNo = i; // boundary element number
ElTr->ElementType = ElementTransformation::BDR_ELEMENT;
ElTr->mesh = this;
DenseMatrix &pm = ElTr->GetPointMat();
ElTr->Reset();
if (Nodes == NULL)
@@ -481,6 +484,7 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
"Mesh requires nodal Finite Element.");
IntegrationRule eir(face_el->GetDof());
FaceElemTr.Loc1.Transf.ElementNo = elem_id;
FaceElemTr.Loc1.Transf.mesh = this;
FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
@@ -495,6 +499,7 @@ void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
FTr->Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
FTr->ElementNo = FaceNo;
FTr->ElementType = ElementTransformation::FACE;
FTr->mesh = this;
DenseMatrix &pm = FTr->GetPointMat();
FTr->Reset();
if (Nodes == NULL)
@@ -551,6 +556,7 @@ void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
IntegrationRule eir(face_el->GetDof());
FaceElemTr.Loc1.Transf.ElementNo = face_info.Elem1No;
FaceElemTr.Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
FaceElemTr.Loc1.Transf.mesh = this;
FaceElemTr.Loc1.Transform(face_el->GetNodes(), eir);
Nodes->GetVectorValues(FaceElemTr.Loc1.Transf, eir, pm);
@@ -580,6 +586,7 @@ void Mesh::GetEdgeTransformation(int EdgeNo, IsoparametricTransformation *EdTr)
EdTr->Attribute = 1;
EdTr->ElementNo = EdgeNo;
EdTr->ElementType = ElementTransformation::EDGE;
EdTr->mesh = this;
DenseMatrix &pm = EdTr->GetPointMat();
EdTr->Reset();
if (Nodes == NULL)
@@ -1098,6 +1105,7 @@ FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
tr->Attribute = boundary[BdrElemNo]->GetAttribute();
tr->ElementNo = BdrElemNo;
tr->ElementType = ElementTransformation::BDR_FACE;
tr->mesh = this;
return tr;
}
@@ -1177,15 +1185,13 @@ void Mesh::Init()
own_nodes = 1;
NURBSext = NULL;
ncmesh = NULL;
ent_sets = NULL;
last_operation = Mesh::NONE;
}
void Mesh::InitTables()
{
el_to_edge =
el_to_face = el_to_el = bel_to_edge = face_edge =
face_vertex = edge_vertex = NULL;
el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
}
void Mesh::SetEmpty()
@@ -1207,7 +1213,6 @@ void Mesh::DestroyTables()
}
delete face_edge;
delete face_vertex;
delete edge_vertex;
}
@@ -1215,8 +1220,6 @@ void Mesh::DestroyPointers()
{
if (own_nodes) { delete Nodes; }
delete ent_sets;
delete ncmesh;
delete NURBSext;
@@ -3351,12 +3354,6 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
// Copy the edge-to-vertex Table, edge_vertex
edge_vertex = (mesh.edge_vertex) ? new Table(*mesh.edge_vertex) : NULL;
// Copy the face-to-vertex Table, edge_vertex
face_vertex = (mesh.face_vertex) ? new Table(*mesh.face_vertex) : NULL;
// Do not copy any of the coarse (c_*), fine (f_*) or fine/coarse (fc_*)
// data members.
// Copy the attributes and bdr_attributes
mesh.attributes.Copy(attributes);
mesh.bdr_attributes.Copy(bdr_attributes);
@@ -3407,9 +3404,6 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
Nodes = mesh.Nodes;
own_nodes = 0;
}
// Copy entity sets if present in the input mesh
ent_sets = (mesh.ent_sets) ? new EntitySets(*mesh.ent_sets) : NULL;
}
Mesh::Mesh(Mesh &&mesh) : Mesh()
@@ -5782,38 +5776,6 @@ Table *Mesh::GetEdgeVertexTable() const
return edge_vertex;
}
Table *Mesh::GetFaceVertexTable() const
{
if (face_vertex)
{
return face_vertex;
}
STable3D * faces_tbl = GetFacesTable();
int nfaces = faces_tbl->NumberOfElements();
face_vertex = new Table(nfaces, 4);
for (int i = 0; i < NumOfVertices; i++)
{
for (STable3D::RowIterator it(*faces_tbl, i); !it; ++it)
{
int j = it.Index();
face_vertex->Push(j, i);
face_vertex->Push(j, it.Column());
face_vertex->Push(j, it.Floor());
if ( it.Tier() > 0 )
{
face_vertex->Push(j, it.Tier());
}
}
}
face_vertex->Finalize();
delete faces_tbl;
return face_vertex;
}
Table *Mesh::GetVertexToElementTable()
{
int i, j, nv, *v;
@@ -6448,7 +6410,7 @@ void Mesh::GenerateNCFaceInfo()
}
}
STable3D *Mesh::GetFacesTable() const
STable3D *Mesh::GetFacesTable()
{
STable3D *faces_tbl = new STable3D(NumOfVertices);
for (int i = 0; i < NumOfElements; i++)
@@ -7703,11 +7665,6 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
}
if ( ent_sets )
{
ent_sets->CopyMeshTables();
}
int quad_counter = 0;
for (int i = 0; i < NumOfElements; i++)
{
@@ -7843,11 +7800,6 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
if (update_nodes) { UpdateNodes(); }
if ( ent_sets )
{
ent_sets->UniformRefinement2D();
}
#ifdef MFEM_DEBUG
if (!Nodes || update_nodes)
{
@@ -7878,11 +7830,6 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
GetElementToFaceTable();
}
if ( ent_sets )
{
ent_sets->CopyMeshTables();
}
Array<int> f2qf_loc;
Array<int> &f2qf = f2qf_ptr ? *f2qf_ptr : f2qf_loc;
f2qf.SetSize(0);
@@ -8209,6 +8156,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
}
AverageVertices(vv, 4, oface + f2qf[f[fi]]);
}
for (int ei = 0; ei < 9; ei++)
{
for (int k = 0; k < 2; k++)
@@ -8552,11 +8500,6 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
sequence++;
if (update_nodes) { UpdateNodes(); }
if (ent_sets)
{
ent_sets->UniformRefinement3D();
}
}
void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
@@ -9026,8 +8969,6 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
mfem::Swap(geom_factors, other.geom_factors);
mfem::Swap(ent_sets, other.ent_sets);
#ifdef MFEM_USE_MEMALLOC
TetMemory.Swap(other.TetMemory);
#endif
+2 -11
View File
@@ -20,7 +20,6 @@
#include "vertex.hpp"
#include "vtk.hpp"
#include "ncmesh.hpp"
#include "entsets.hpp"
#include "../fem/eltrans.hpp"
#include "../fem/coefficient.hpp"
#include "../general/zstr.hpp"
@@ -55,11 +54,9 @@ class Mesh
#ifdef MFEM_USE_MPI
friend class ParMesh;
friend class ParNCMesh;
friend class ParEntitySets;
#endif
friend class NCMesh;
friend class NURBSExtension;
friend class EntitySets;
#ifdef MFEM_USE_ADIOS2
friend class adios2stream;
@@ -169,7 +166,6 @@ protected:
Array<int> be_to_face;
mutable Table *face_edge;
mutable Table *edge_vertex;
mutable Table *face_vertex;
IsoparametricTransformation Transformation, Transformation2;
IsoparametricTransformation BdrTransformation;
@@ -220,8 +216,6 @@ public:
Array<FaceGeometricFactors*>
face_geom_factors; ///< Optional face geometric factors.
EntitySets *ent_sets;
// Global parameter that can be used to control the removal of unused
// vertices performed when reading a mesh in MFEM format. The default value
// (true) is set in mesh_readers.cpp.
@@ -293,7 +287,7 @@ protected:
void PrepareNodeReorder(DSTable **old_v_to_v, Table **old_elem_vert);
void DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert);
STable3D *GetFacesTable() const;
STable3D *GetFacesTable();
STable3D *GetElementToFaceTable(int ret_ftbl = 0);
/** Red refinement. Element with index i is refined. The default
@@ -1073,12 +1067,9 @@ public:
/// Returns the face-to-edge Table (3D)
Table *GetFaceEdgeTable() const;
/// Returns the edge-to-vertex Table (2D or 3D)
/// Returns the edge-to-vertex Table (3D)
Table *GetEdgeVertexTable() const;
/// Returns the face-to-vertex Table (2d or 3D)
Table *GetFaceVertexTable() const;
/// Return the indices and the orientations of all faces of element i.
void GetElementFaces(int i, Array<int> &faces, Array<int> &ori) const;
-8
View File
@@ -100,14 +100,6 @@ void Mesh::ReadMFEMMesh(std::istream &input, int version, int &curved)
curved = 1;
}
ent_sets = new EntitySets(*this);
ent_sets->Load(input);
if ( ent_sets->GetNumSets(EntitySets::FACE) > 0 && faces.Size() == 0 )
{
GetElementToFaceTable();
GenerateFaces();
}
// When visualizing solutions on non-conforming grids, PETSc
// may dump additional vertices
if (remove_unused_vertices) { RemoveUnusedVertices(); }
+4 -404
View File
@@ -185,10 +185,6 @@ NCMesh::NCMesh(const Mesh *mesh)
face->attribute = be->GetAttribute();
}
// Store entity set information if present in the Mesh
ncent_sets = (mesh->ent_sets) ?
new NCEntitySets(*mesh->ent_sets, *this) : NULL;
// copy top-level vertex coordinates (leave empty if the mesh is curved)
if (!mesh->Nodes)
{
@@ -220,10 +216,6 @@ NCMesh::NCMesh(const NCMesh &other)
other.free_element_ids.Copy(free_element_ids);
other.root_state.Copy(root_state);
other.coordinates.Copy(coordinates);
// Copy the entity set information
ncent_sets = (other.ncent_sets) ? new NCEntitySets(*other.ncent_sets) : NULL;
Update();
}
@@ -262,11 +254,8 @@ NCMesh::~NCMesh()
DeleteUnusedFaces(elemFaces);
}
}
// NOTE: in release mode, we just throw away all faces and nodes at once
#endif
delete ncent_sets;
}
NCMesh::Node::~Node()
@@ -2515,42 +2504,6 @@ void NCMesh::OnMeshUpdated(Mesh *mesh)
if (face->index < 0) { face->index = NFaces + (nghosts++); }
}
MFEM_ASSERT(nghosts == NGhostFaces, "");
if (ncent_sets)
{
std::cout << "NCMesh::OnMeshUpdated ncent_sets is non NULL" << std::endl;
if (!mesh->ent_sets)
{
std::cout << "NCMesh::OnMeshUpdated creating ent_sets from NCMesh" << std::endl;
mesh->ent_sets = new EntitySets(*mesh, *this);
std::cout << "NCMesh::OnMeshUpdated done creating ent_sets from NCMesh" <<
std::endl;
}
}
std::ostringstream ossN;
ossN << "node_on_mesh_updated.out";
std::ofstream ofsN(ossN.str().c_str());
ofsN << nodes.Size() << std::endl;
for (int i=0; i<nodes.Size(); i++)
{
ofsN << i
// << " " << nodes[i].vert_refc
// << " " << nodes[i].edge_refc
<< " " << nodes[i].HasVertex()
<< " " << nodes[i].HasEdge()
<< " " << nodes[i].vert_index
<< " " << nodes[i].edge_index
<< " " << nodes[i].p1
<< " " << nodes[i].p2
<< " " << nodes[i].next << std::endl;
}
ofsN.close();
NEdges = mesh->GetNEdges();
NFaces = mesh->GetNumFaces();
std::cout << "Leaving NCMesh::OnMeshUpdated" << std::endl;
}
@@ -3353,15 +3306,12 @@ const NCMesh::MeshId& NCMesh::NCList::LookUp(int index, int *type) const
void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)
{
int mid = nodes.FindId(v0, v1);
if (mid >= 0)
if (mid >= 0 && nodes[mid].HasVertex())
{
if (nodes[mid].HasVertex())
{
indices.Append(mid);
indices.Append(mid);
CollectEdgeVertices(v0, mid, indices);
CollectEdgeVertices(mid, v1, indices);
}
CollectEdgeVertices(v0, mid, indices);
CollectEdgeVertices(mid, v1, indices);
}
}
@@ -3423,78 +3373,6 @@ void NCMesh::CollectQuadFaceVertices(int v0, int v1, int v2, int v3,
}
}
void NCMesh::CollectElementVertices(int elem_id, Array<int> &indices)
{
Element &el = elements[elem_id];
if (el.ref_type != 0)
{
// This element has been refined so recurse into its children
for (int i = 0; i < 8; i++)
{
if (el.child[i] >= 0 && el.child[i] < elements.Size())
{
CollectElementVertices(el.child[i], indices);
}
}
}
else
{
// This element has not been refined so add its vertices
for (int i=0; i<8; i++)
{
if (el.node[i] >= 0 && el.node[i] < nodes.Size())
{
indices.Append(el.node[i]);
}
}
}
}
void NCMesh::CollectElementEdges(int elem_id, Array<int> &indices)
{
Element &el = elements[elem_id];
if (el.ref_type != 0)
{
// This element has been refined so recurse into its children
for (int i = 0; i < 8; i++)
{
if (el.child[i] >= 0 && el.child[i] < elements.Size())
{
CollectElementEdges(el.child[i], indices);
}
}
}
else
{
int* node = el.node;
GeomInfo& gi = GI[(int) el.geom];
for (int i = 0; i < gi.nv; i++)
{
if (nodes[node[i]].HasEdge())
{
indices.Append(node[i]);
}
}
for (int i = 0; i < gi.ne; i++)
{
const int* ev = gi.edges[i];
int index = nodes.FindId(node[ev[0]], node[ev[1]]);
if (index >= 0)
{
if (nodes[index].HasEdge())
{
indices.Append(index);
}
}
}
}
}
void NCMesh::BuildElementToVertexTable()
{
int nrows = leaf_elements.Size();
@@ -4996,107 +4874,6 @@ int NCMesh::GetElementDepth(int i) const
return depth;
}
void NCMesh::GetRefinedEdges(int vn0, int vn1, BlockArray<int> & edges)
{
std::cout << "entering NCMesh::GetRefinedEdges "
<<"searching for edge with vertices: " << vn0 << " and " << vn1
<< std::endl;
int mid = nodes.FindId(vn0, vn1);
if (mid < 0) { return; }
Node &nd = nodes[mid];
// if ( nd.edge_index < 0 ) { return; }
// edges.Append(nd.edge_index);
if ( nd.HasEdge() )
{
std::cout << " found node " << mid << std::endl;
edges.Append(mid);
}
GetRefinedEdges(vn0, mid, edges);
GetRefinedEdges(mid, vn1, edges);
}
void NCMesh::GetRefinedFaces(int vn0, int vn1, int vn2, int vn3,
BlockArray<int> & face_ids)
{
// Face* fa = faces.Find(vn0, vn1, vn2, vn3);
int face = faces.FindId(vn0, vn1, vn2, vn3);
/*
if (fa)
{
if ( fa->index >= 0 )
{
face_ids.Append(fa->index);
}
return;
}
*/
if (face>=0)
{
if ( faces[face].index >= 0 )
{
face_ids.Append(face);
}
return;
}
// we need to recurse deeper
int mid[4];
int split = QuadFaceSplitType(vn0, vn1, vn2, vn3, mid);
if (split == 1) // "X" split face
{
GetRefinedFaces(vn0, mid[0], mid[2], vn3, face_ids);
GetRefinedFaces(mid[0], vn1, vn2, mid[2], face_ids);
}
else if (split == 2) // "Y" split face
{
GetRefinedFaces(vn0, vn1, mid[1], mid[3], face_ids);
GetRefinedFaces(mid[3], mid[1], vn2, vn3, face_ids);
}
}
void NCMesh::GetRefinedElements(int elem_id, BlockArray<int> & elem_ids)
{
// std::cout << "entering NCMesh::GetRefinedElements searching for element id: "
// << elem_id << std::endl;
Element &el = elements[elem_id];
/*
if (el.index >= 0 && el.rank >= 0)
{
elem_ids.Append(el.index);
return;
}
for (int i = 0; i < 8; i++)
{
if (el.child[i] >= 0 && el.child[i] < elements.Size() )
{
GetRefinedElements(el.child[i], elem_ids);
}
}
*/
if (el.ref_type != 0)
{
// This element has been refined so recurse into its children
for (int i = 0; i < 8; i++)
{
if (el.child[i] >= 0 && el.child[i] < elements.Size() )
{
GetRefinedElements(el.child[i], elem_ids);
}
}
}
else
{
// This element has not been refined so add it
elem_ids.Append(elem_id);
}
}
int NCMesh::GetElementSizeReduction(int i) const
{
int elem = leaf_elements[i];
@@ -5216,183 +4993,6 @@ void NCMesh::GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
bdr_edges.Unique();
}
void NCMesh::GetEntitySetClosure(EntitySets::EntityType type,
int set_index,
Array<int> &es_vertices,
Array<int> &es_edges,
Array<int> &es_faces)
{
es_vertices.SetSize(0);
es_edges.SetSize(0);
es_faces.SetSize(0);
MFEM_VERIFY(ncent_sets != NULL, "NCMesh object contains no "
"entity set information");
if (!ncent_sets->SetExists(type, set_index))
{
std::ostringstream oss; oss << "Entity set of type \""
<< EntitySets::GetTypeName(type)
<< "\" and index " << set_index
<< " was not found.";
MFEM_VERIFY(false, oss.str().c_str());
}
int ni = ncent_sets->GetNumEntities(type ,set_index);
Array<int> inds;
Array<int> coll_inds;
switch (type)
{
case EntitySets::VERTEX:
{
/// Do nothing because vertices cannot hide
}
break;
case EntitySets::EDGE:
{
for (int i=0; i<ni; i++)
{
ncent_sets->GetEntityIndex(type, set_index, i, inds);
// collect vertices
inds.Copy(coll_inds);
this->CollectEdgeVertices(inds[0], inds[1], coll_inds);
for (int j=0; j<coll_inds.Size(); j++)
{
int index = nodes[coll_inds[j]].vert_index;
if (index >= 0)
{
es_vertices.Append(index);
}
}
}
}
break;
case EntitySets::FACE:
{
for (int i=0; i<ni; i++)
{
ncent_sets->GetEntityIndex(type, set_index, i, inds);
// collect vertices
inds.Copy(coll_inds);
if (inds.Size() == 4)
{
this->CollectQuadFaceVertices(inds[0], inds[1], inds[2], inds[3],
coll_inds);
}
else
{
this->CollectTriFaceVertices(inds[0], inds[1], inds[2],
coll_inds);
}
for (int j=0; j<coll_inds.Size(); j++)
{
int index = nodes[coll_inds[j]].vert_index;
if (index >= 0)
{
es_vertices.Append(index);
}
}
}
}
break;
case EntitySets::ELEMENT:
{
for (int i=0; i<ni; i++)
{
int elem_id = (*ncent_sets)(type, set_index, i);
std::cout << "examining element " << elem_id << std::endl;
// collect vertices
coll_inds.SetSize(0);
this->CollectElementVertices(elem_id, coll_inds);
for (int j=0; j<coll_inds.Size(); j++)
{
int index = nodes[coll_inds[j]].vert_index;
if (index >= 0)
{
es_vertices.Append(index);
}
}
// collect edges
coll_inds.SetSize(0);
this->CollectElementEdges(elem_id, coll_inds);
for (int j=0; j<coll_inds.Size(); j++)
{
int index = nodes[coll_inds[j]].edge_index;
if (index >= 0)
{
es_edges.Append(index);
}
}
}
}
break;
default:
MFEM_ABORT("GetEnitySetClosure - Unknown entity set type: \""
<< EntitySets::GetTypeName(type) << "\"");
}
/*
if (Dim == 3)
{
GetFaceList(); // make sure 'boundary_faces' is up to date
for (int i = 0; i < boundary_faces.Size(); i++)
{
int face = boundary_faces[i];
if (bdr_attr_is_ess[faces[face].attribute - 1])
{
int node[4];
FindFaceNodes(face, node);
for (int j = 0; j < 4; j++)
{
bdr_vertices.Append(nodes[node[j]].vert_index);
int enode = nodes.FindId(node[j], node[(j+1) % 4]);
MFEM_ASSERT(enode >= 0 && nodes[enode].HasEdge(), "Edge not found.");
bdr_edges.Append(nodes[enode].edge_index);
while ((enode = GetEdgeMaster(enode)) >= 0)
{
// append master edges that may not be accessible from any
// boundary element, this happens in 3D in re-entrant corners
bdr_edges.Append(nodes[enode].edge_index);
}
}
}
}
}
else if (Dim == 2)
{
GetEdgeList(); // make sure 'boundary_faces' is up to date
for (int i = 0; i < boundary_faces.Size(); i++)
{
int face = boundary_faces[i];
Face &fc = faces[face];
if (bdr_attr_is_ess[fc.attribute - 1])
{
bdr_vertices.Append(nodes[fc.p1].vert_index);
bdr_vertices.Append(nodes[fc.p3].vert_index);
}
}
}
*/
es_vertices.Sort();
es_vertices.Unique();
es_edges.Sort();
es_edges.Unique();
es_faces.Sort();
es_faces.Unique();
}
static int max4(int a, int b, int c, int d)
{
return std::max(std::max(a, b), std::max(c, d));
+2 -35
View File
@@ -19,7 +19,6 @@
#include "../linalg/densemat.hpp"
#include "element.hpp"
#include "vertex.hpp"
#include "entsets.hpp"
#include "../fem/geom.hpp"
#include <vector>
@@ -118,9 +117,6 @@ struct MatrixMap; // for internal use
*/
class NCMesh
{
friend class EntitySets;
friend class NCEntitySets;
public:
//// Initialize with elements from an existing 'mesh'.
explicit NCMesh(const Mesh *mesh);
@@ -347,16 +343,6 @@ public:
Array<int> &bdr_vertices,
Array<int> &bdr_edges);
/** Get a list of vertices (2D/3D), edges (2D/3D), and faces (3D) that
coincide with members of the specified entity set. In 3D this function
also reveals "hidden" edges or faces. In parallel it helps identifying
vertices/edges/faces affected by non-local entities. */
virtual void GetEntitySetClosure(EntitySets::EntityType t,
int set_index,
Array<int> &es_vertices,
Array<int> &es_edges,
Array<int> &es_faces);
/// Return element geometry type. @a index is the Mesh element number.
Geometry::Type GetElementGeometry(int index) const
{ return elements[leaf_elements[index]].Geom(); }
@@ -371,19 +357,6 @@ public:
/// Return the distance of leaf 'i' from the root.
int GetElementDepth(int i) const;
/** Collect edge indices of all refined edges which are children of
the coarse edge defined by the given vertices. */
void GetRefinedEdges(int vn0, int vn1, BlockArray<int> & edge_ids);
/** Collect face indices of all refined faces which are children of
the coarse face defined by the given vertices. */
void GetRefinedFaces(int vn0, int vn1, int vn2, int vn3,
BlockArray<int> & face_ids);
/** Collect element indices of all refined elements which are children of
the coarse element defined by the given element index. */
void GetRefinedElements(int elem_id, BlockArray<int> & elem_ids);
/** Return the size reduction compared to the root element (ignoring local
stretching and curvature). */
int GetElementSizeReduction(int i) const;
@@ -528,7 +501,6 @@ protected: // implementation
Array<double> coordinates;
// secondary data
/** Apart from the primary data structure, which is the element/node/face
@@ -558,8 +530,6 @@ protected: // implementation
Table element_vertex; ///< leaf-element to vertex table, see FindSetNeighbors
// Node/edge/Face/Element sets defined on the coarse mesh
NCEntitySets * ncent_sets;
void UpdateLeafElements();
void UpdateVertices(); ///< update Vertex::index and vertex_nodeId
@@ -741,10 +711,6 @@ protected: // implementation
void CollectTriFaceVertices(int v0, int v1, int v2, Array<int> &indices);
void CollectQuadFaceVertices(int v0, int v1, int v2, int v3,
Array<int> &indices);
void CollectElementVertices(int elem_id, Array<int> &indices);
void CollectElementEdges(int elem_id, Array<int> &indices);
void BuildElementToVertexTable();
void UpdateElementToVertexTable()
@@ -771,6 +737,8 @@ protected: // implementation
Point() { dim = 0; }
Point(const Point &) = default;
Point(double x)
{ dim = 1; coord[0] = x; }
@@ -960,7 +928,6 @@ public:
#endif
friend class ParNCMesh; // for ParNCMesh::ElementSet
friend class ParNCEntitySets;
friend struct MatrixMap;
friend struct PointMatrixHash;
};
-392
View File
@@ -1,392 +0,0 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#include "../config/config.hpp"
#ifdef MFEM_USE_MPI
#include "pentsets.hpp"
#include "pmesh.hpp"
using namespace std;
namespace mfem
{
ParEntitySets::ParEntitySets(const ParEntitySets & ent_sets)
: EntitySets(ent_sets),
pmesh_(ent_sets.GetParMesh())
{
MPI_Comm_size(pmesh_->GetComm(), &NRanks_);
MPI_Comm_rank(pmesh_->GetComm(), &MyRank_);
cout << MyRank_ << ": Entering ParEntitySets copy c'tor" << endl;
cout << MyRank_ << ": Leaving ParEntitySets copy c'tor" << endl;
}
ParEntitySets::ParEntitySets(ParMesh & pmesh, const EntitySets & ent_sets,
int * partitioning,
const Array<int> & vert_global_local)
: EntitySets(ent_sets),
pmesh_(&pmesh)
{
// The copy constructor for EntitySets will initialize this object's
// data with the correct set names, and numbers of sets. However,
// the set entries themselves will need to be recomputed based on
// local numberings and the paritioning.
//
// The EntitySets object will be a copy of the serial object. This
// constructor will have to prune and renumber the data. Once this
// is done the mesh pointer stored in the EntitySets object can be
// replaced with the local portion of the parallel mesh.
MPI_Comm MyComm = pmesh_->GetComm();
MPI_Comm_size(MyComm, &NRanks_);
MPI_Comm_rank(MyComm, &MyRank_);
cout << MyRank_ << ": Entering ParEntitySets(ParMesh, EntitySets, ...) c'tor" <<
endl;
int nelem = mesh_->GetNE();
DSTable v_to_v(vert_global_local.Size());
pmesh_->GetVertexToVertexTable(v_to_v);
STable3D * faces_tbl = NULL;
const Table * serial_edge_vertex = NULL;
const Table * serial_face_vertex = NULL;
if ( ent_sets.GetNumSets(EDGE) > 0 )
{
serial_edge_vertex = ent_sets.GetEdgeVertexTable();
}
if ( ent_sets.GetNumSets(FACE) > 0 )
{
serial_face_vertex = ent_sets.GetFaceVertexTable();
faces_tbl = pmesh_->GetFacesTable();
}
Array<int> elem_global_local(nelem);
elem_global_local = -1;
int elem_counter = 0;
for (int i=0; i<nelem; i++)
{
if ( partitioning[i] == MyRank_ )
{
elem_global_local[i] = elem_counter;
elem_counter++;
}
}
EntityType t;
unsigned int ns;
t = VERTEX;
ns = ent_sets.GetNumSets(t);
for (unsigned int s=0; s<ns; s++)
{
set<int>::iterator it;
sets_[t][s].clear();
for (it=ent_sets(t,s).begin(); it!=ent_sets(t,s).end(); it++)
{
int v0 = vert_global_local[*it];
if ( v0 >= 0 )
{
sets_[t][s].insert(v0);
}
}
}
if ( pmesh_->Dimension() > 1 )
{
t = EDGE;
ns = ent_sets.GetNumSets(t);
for (unsigned int s=0; s<ns; s++)
{
set<int>::iterator it;
sets_[t][s].clear();
for (it=ent_sets(t,s).begin(); it!=ent_sets(t,s).end(); it++)
{
int old_edge = *it;
const int *v = serial_edge_vertex->GetRow(old_edge);
int v0 = vert_global_local[v[0]];
int v1 = vert_global_local[v[1]];
if ( v0 >= 0 && v1 >= 0 )
{
int new_edge = v_to_v(v0,v1);
if ( new_edge >= 0 )
{
sets_[t][s].insert(new_edge);
}
}
}
}
}
if ( pmesh_->Dimension() > 2 )
{
Array<int> v;
t = FACE;
ns = ent_sets.GetNumSets(t);
for (unsigned int s=0; s<ns; s++)
{
set<int>::iterator it;
sets_[t][s].clear();
for (it=ent_sets(t,s).begin(); it!=ent_sets(t,s).end(); it++)
{
int old_face = *it;
int numv = serial_face_vertex->RowSize(old_face);
const int *v = serial_face_vertex->GetRow(old_face);
if ( vert_global_local[v[0]] >= 0 &&
vert_global_local[v[1]] >= 0 &&
vert_global_local[v[2]] >= 0 )
{
int new_face = -1;
if ( numv == 3 )
{
new_face = (*faces_tbl)(vert_global_local[v[0]],
vert_global_local[v[1]],
vert_global_local[v[2]]);
}
else
{
new_face = (*faces_tbl)(vert_global_local[v[0]],
vert_global_local[v[1]],
vert_global_local[v[2]],
vert_global_local[v[3]]);
}
if ( new_face >= 0 )
{
sets_[t][s].insert(new_face);
}
}
}
}
delete faces_tbl;
}
t = ELEMENT;
ns = ent_sets.GetNumSets(t);
for (unsigned int s=0; s<ns; s++)
{
set<int>::iterator it;
sets_[t][s].clear();
for (it=ent_sets(t,s).begin(); it!=ent_sets(t,s).end(); it++)
{
if ( partitioning[*it] == MyRank_ )
{
sets_[t][s].insert(elem_global_local[*it]);
}
}
}
this->mesh_ = (Mesh*)this->pmesh_;
this->CopyMeshTables();
cout << MyRank_ << ": Leaving ParEntitySets(ParMesh, EntitySets, ...) c'tor" <<
endl;
}
ParEntitySets::ParEntitySets(ParMesh & pmesh, ParNCMesh &pncmesh)
: EntitySets(pmesh),
pmesh_(&pmesh)
{
MPI_Comm MyComm = pmesh_->GetComm();
MPI_Comm_size(MyComm, &NRanks_);
MPI_Comm_rank(MyComm, &MyRank_);
cout << MyRank_ << ": Entering ParEntitySets(ParMesh, ParNCMesh) c'tor" << endl;
this->BuildEntitySets(pncmesh, VERTEX);
this->BuildEntitySets(pncmesh, EDGE);
this->BuildEntitySets(pncmesh, FACE);
this->BuildEntitySets(pncmesh, ELEMENT);
cout << MyRank_ << ": Leaving ParEntitySets(ParMesh, ParNCMesh) c'tor" << endl;
}
ParEntitySets::~ParEntitySets()
{
cout << MyRank_ << ": Entering ParEntitySets d'tor" << endl;
cout << MyRank_ << ": Leaving ParEntitySets d'tor" << endl;
}
void
ParEntitySets::PrintSetInfo(std::ostream & output) const
{
if ( MyRank_ == 0 &&
( GetNumSets(VERTEX) > 0 || GetNumSets(EDGE) > 0 ||
GetNumSets(FACE) > 0 || GetNumSets(ELEMENT) > 0 ) )
{
output << "\nMFEM Parallel Entity Sets:\n";
}
this->PrintEntitySetInfo(output, VERTEX, "Vertex");
this->PrintEntitySetInfo(output, EDGE, "Edge");
this->PrintEntitySetInfo(output, FACE, "Face");
this->PrintEntitySetInfo(output, ELEMENT, "Element");
}
void
ParEntitySets::PrintEntitySetInfo(std::ostream & output, EntityType t,
const string & ent_name) const
{
if ( sets_[t].size() > 0 )
{
if ( MyRank_ == 0 )
{
output << " " << ent_name
<< " Sets (Index, Set Name, Global Size):\n";
}
for (unsigned int s=0; s<sets_[t].size(); s++)
{
int loc_size = sets_[t][s].size();
int glb_size = -1;
MPI_Reduce(&loc_size, &glb_size, 1, MPI_INT, MPI_SUM, 0,
pmesh_->GetComm());
if ( MyRank_ == 0 )
{
output << '\t' << s
<< '\t' << set_names_[t][s]
<< '\t' << glb_size
<< '\n';
}
}
if ( MyRank_ == 0 )
{
output << '\n';
}
}
}
void
ParEntitySets::BuildEntitySets(ParNCMesh &pncmesh, EntityType t)
{
cout << MyRank_ << ": BuildEntitySets for type " << GetTypeName(t) << endl;
int es = pncmesh.pncent_sets->GetEntitySize(t);
unsigned int ns = pncmesh.pncent_sets->GetNumSets(t);
cout << MyRank_ << ": num sets " << ns << endl;
Array<int> inds(es);
sets_[t].resize(ns);
set_names_[t].resize(ns);
for (unsigned int s=0; s<ns; s++)
{
int ni = pncmesh.pncent_sets->GetNumEntities(t, s);
set_names_[t][s] = pncmesh.pncent_sets->GetSetName(t, s);
set_index_by_name_[t][set_names_[t][s]] = s;
switch (t)
{
case VERTEX:
for (int i=0; i<ni; i++)
{
int node = (*pncmesh.pncent_sets)(t, s, i);
int index = pncmesh.nodes[node].vert_index;
if (!pncmesh.IsGhost(0,index))
{
sets_[t][s].insert(index);
}
}
break;
case EDGE:
for (int i=0; i<ni; i++)
{
pncmesh.pncent_sets->GetEntityIndex(t, s, i, inds);
BlockArray<int> ind_coll;
pncmesh.GetRefinedEdges(inds[0], inds[1],
ind_coll);
for (int j=0; j<ind_coll.Size(); j++)
{
int edge = ind_coll[j];
int index = pncmesh.nodes[edge].edge_index;
if (index >= 0 && !pncmesh.IsGhost(1, index))
{
sets_[t][s].insert(index);
}
}
}
break;
case FACE:
for (int i=0; i<ni; i++)
{
pncmesh.pncent_sets->GetEntityIndex(t, s, i, inds);
BlockArray<int> ind_coll;
pncmesh.GetRefinedFaces(inds[0], inds[1], inds[2], inds[3],
ind_coll);
for (int j=0; j<ind_coll.Size(); j++)
{
int face = ind_coll[j];
int index = pncmesh.faces[face].index;
if (index >= 0 && !pncmesh.IsGhost(2, index))
{
sets_[t][s].insert(index);
}
}
}
break;
case ELEMENT:
for (int i=0; i<ni; i++)
{
int elem = (*pncmesh.pncent_sets)(t, s, i);
BlockArray<int> ind_coll;
pncmesh.GetRefinedElements(elem, ind_coll);
for (int j=0; j<ind_coll.Size(); j++)
{
sets_[t][s].insert(pncmesh.elements[ind_coll[j]].index);
}
}
break;
default:
MFEM_ABORT("Unknown entity set type: \"" << GetTypeName(t) << "\"");
}
cout << MyRank_ << ": " << set_names_[t][s] << " " << s << " set size " <<
sets_[t][s].size() << "{";
for (set<int>::iterator it=sets_[t][s].begin(); it!=sets_[t][s].end(); it++)
{
cout << " " << *it;
}
cout << "}" << endl;
}
map<string,int>::iterator it;
cout << MyRank_ << ": set index by name ";
for (it=set_index_by_name_[t].begin(); it != set_index_by_name_[t].end(); it++)
{
cout << " " << it->first << "->" << it->second;
}
cout << endl;
cout << MyRank_ << ": done BuildEntitySets for type " << GetTypeName(t) << endl;
}
ParNCEntitySets::ParNCEntitySets(MPI_Comm comm, const NCMesh &ncmesh)
: NCEntitySets(*ncmesh.ncent_sets)
{
MyComm_ = comm;
MPI_Comm_size(MyComm_, &NRanks_);
MPI_Comm_rank(MyComm_, &MyRank_);
if ( MyRank_ == 0 )
{
cout << "Entering ParNCEntitySets(NCMesh) c'tor" << endl;
}
if ( MyRank_ == 0 )
{
cout << "Leaving ParNCEntitySets(NCMesh) c'tor" << endl;
}
}
} // namespace mfem
#endif // MFEM_USE_MPI
-74
View File
@@ -1,74 +0,0 @@
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
#ifndef MFEM_PAR_ENTITY_SETS
#define MFEM_PAR_ENTITY_SETS
#include "../config/config.hpp"
#ifdef MFEM_USE_MPI
#include "entsets.hpp"
#include "../general/communication.hpp"
namespace mfem
{
class ParMesh;
class ParNCMesh;
class ParEntitySets : public EntitySets
{
friend class ParMesh;
public:
ParEntitySets(const ParEntitySets & ent_sets);
ParEntitySets(ParMesh & _mesh, const EntitySets & ent_sets, int * part,
const Array<int> & vert_global_local);
ParEntitySets(ParMesh & mesh, ParNCMesh &ncmesh);
virtual ~ParEntitySets();
virtual void PrintSetInfo(std::ostream &output) const;
inline ParMesh *GetParMesh() const { return pmesh_; }
private:
void PrintEntitySetInfo(std::ostream & output, EntityType t,
const std::string & ent_name) const;
void BuildEntitySets(ParNCMesh &pncmesh, EntityType t);
ParMesh * pmesh_;
int NRanks_;
int MyRank_;
};
class ParNCEntitySets : public NCEntitySets
{
public:
// ParNCEntitySets(MPI_Comm comm, EntitySets &ent_sets, NCMesh &ncmesh);
ParNCEntitySets(MPI_Comm comm, const NCMesh &ncmesh);
// ParNCEntitySets(const ParMesh & pmesh, const ParNCMesh &pncmesh);
// ParNCEntitySets(const ParNCEntitySets & pncent_sets);
private:
MPI_Comm MyComm_;
int NRanks_;
int MyRank_;
};
} // namespace mfem
#endif // MFEM_USE_MPI
#endif // MFEM_PAR_ENTITY_SETS
+8 -37
View File
@@ -91,10 +91,6 @@ ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
*Nodes = *pmesh.Nodes;
own_nodes = 1;
}
// Copy entity sets if present in the input mesh
ent_sets = pent_sets =
(pmesh.pent_sets) ? new ParEntitySets(*pmesh.pent_sets) : NULL;
}
ParMesh::ParMesh(ParMesh &&mesh) : ParMesh()
@@ -114,7 +110,6 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
, glob_elem_offset(-1)
, glob_offset_sequence(-1)
, gtopo(comm)
, pent_sets(NULL)
{
int *partitioning = NULL;
Array<bool> activeBdrElem;
@@ -123,8 +118,6 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
MPI_Comm_size(MyComm, &NRanks);
MPI_Comm_rank(MyComm, &MyRank);
Array<int> vert_global_local;
if (mesh.Nonconforming())
{
if (partitioning_)
@@ -155,10 +148,6 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
mesh.bdr_attributes.Copy(bdr_attributes);
GenerateNCFaceInfo();
// if (mesh.ent_sets)
// NumOfVertices = BuildLocalVertices(mesh, partitioning,
// vert_global_local);
}
else // mesh.Conforming()
{
@@ -179,6 +168,7 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
// re-enumerate the partitions to better map to actual processor
// interconnect topology !?
Array<int> vert_global_local;
NumOfVertices = BuildLocalVertices(mesh, partitioning, vert_global_local);
NumOfElements = BuildLocalElements(mesh, partitioning, vert_global_local);
@@ -250,12 +240,6 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
SetMeshGen();
meshgen = mesh.meshgen; // copy the global 'meshgen'
ent_sets = pent_sets =
(mesh.ent_sets) ? new ParEntitySets(*this, *mesh.ent_sets,
partitioning,
vert_global_local)
: NULL;
}
if (mesh.NURBSext)
@@ -305,12 +289,7 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
// for compatibility (e.g., Mesh::GetVertex())
SetVerticesFromNodes(Nodes);
}
/*
ent_sets = pent_sets =
(mesh.ent_sets) ? new ParEntitySets(*this, *mesh.ent_sets,
partitioning,
vert_global_local) : NULL;
*/
if (partitioning != partitioning_)
{
delete [] partitioning;
@@ -880,7 +859,6 @@ ParMesh::ParMesh(const ParNCMesh &pncmesh)
, glob_offset_sequence(-1)
, gtopo(MyComm)
, pncmesh(NULL)
, pent_sets(NULL)
{
Mesh::InitFromNCMesh(pncmesh);
ReduceMeshGen();
@@ -947,7 +925,6 @@ ParMesh::ParMesh(MPI_Comm comm, istream &input, bool refine)
, glob_elem_offset(-1)
, glob_offset_sequence(-1)
, gtopo(comm)
, pent_sets(NULL)
{
MyComm = comm;
MPI_Comm_size(MyComm, &NRanks);
@@ -1162,8 +1139,7 @@ void ParMesh::MakeRefined_(ParMesh &orig_mesh, int ref_factor, int ref_type)
gtopo = orig_mesh.gtopo;
have_face_nbr_data = false;
pncmesh = NULL;
pent_sets = NULL;
Array<int> ref_factors(orig_mesh.GetNE());
ref_factors = ref_factor;
Mesh::MakeRefined_(orig_mesh, ref_factors, ref_type);
@@ -1914,6 +1890,7 @@ void ParMesh::GetFaceNbrElementTransformation(
ElTr->Attribute = elem->GetAttribute();
ElTr->ElementNo = NumOfElements + i;
ElTr->ElementType = ElementTransformation::ELEMENT;
ElTr->mesh = this;
ElTr->Reset();
if (Nodes == NULL)
@@ -2735,6 +2712,7 @@ STable3D *ParMesh::GetFaceNbrElementToFaceTable(int ret_ftbl)
}
face_nbr_el_to_face->Finalize();
delete sfaces_tbl;
if (ret_ftbl)
{
return faces_tbl;
@@ -3792,13 +3770,6 @@ void ParMesh::NonconformingRefinement(const Array<Refinement> &refinements,
// and this mesh will be the new fine mesh
Mesh::Swap(*pmesh2, false);
// swap entity set information if present
mfem::Swap(pmesh2->pent_sets, this->pent_sets);
if (this->pent_sets)
{
this->pent_sets->pmesh_ = this;
}
delete pmesh2; // NOTE: old face neighbors destroyed here
pncmesh->GetConformingSharedStructures(*this);
@@ -6202,9 +6173,6 @@ void ParMesh::Destroy()
delete pncmesh;
ncmesh = pncmesh = NULL;
delete pent_sets;
ent_sets = pent_sets = NULL;
DeleteFaceNbrData();
for (int i = 0; i < shared_edges.Size(); i++)
@@ -6212,6 +6180,9 @@ void ParMesh::Destroy()
FreeElement(shared_edges[i]);
}
shared_edges.DeleteAll();
delete face_nbr_el_to_face;
face_nbr_el_to_face = NULL;
}
ParMesh::~ParMesh()
-2
View File
@@ -20,7 +20,6 @@
#include "../general/globals.hpp"
#include "mesh.hpp"
#include "pncmesh.hpp"
#include "pentsets.hpp"
#include <iostream>
namespace mfem
@@ -321,7 +320,6 @@ public:
Table send_face_nbr_vertices;
ParNCMesh* pncmesh;
ParEntitySets* pent_sets;
int GetNGroups() const { return gtopo.NGroups(); }
-508
View File
@@ -20,8 +20,6 @@
#include <map>
#include <climits> // INT_MIN, INT_MAX
#include <fstream> // MLS Debugging
namespace mfem
{
@@ -29,7 +27,6 @@ using namespace bin_io;
ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh, int *part)
: NCMesh(ncmesh)
, pncent_sets(NULL)
{
MyComm = comm;
MPI_Comm_size(MyComm, &NRanks);
@@ -44,40 +41,6 @@ ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh, int *part)
Update();
std::ostringstream oss; oss << "elements_" << MyRank << ".out";
std::ofstream ofs(oss.str().c_str());
for (int i=0; i<elements.Size(); i++)
{
ofs << i
<< '\t' << elements[i].index
<< '\t' << elements[i].rank
<< '\t' << elements[i].attribute
<< '\t' << elements[i].parent;
if ( elements[i].ref_type == 0 )
{
ofs << " nodes {";
for (int j=0; j<8; j++)
{
ofs << " " << elements[i].node[j];
}
ofs << "}";
}
else
{
ofs << " children {";
for (int j=0; j<8; j++)
{
ofs << " " << elements[i].child[j];
}
ofs << "}";
}
ofs << std::endl;
}
ncent_sets = pncent_sets =
(ncmesh.ncent_sets) ? new ParNCEntitySets(comm, ncmesh) : NULL;
// note that at this point all processors still have all the leaf elements;
// we however may now start pruning the refinement tree to get rid of
// branches that only contain someone else's leaves (see Prune())
@@ -122,9 +85,6 @@ ParNCMesh::ParNCMesh(const ParNCMesh &other)
ParNCMesh::~ParNCMesh()
{
ClearAuxPM();
delete pncent_sets;
ncent_sets = pncent_sets = NULL;
}
void ParNCMesh::Update()
@@ -155,386 +115,6 @@ void ParNCMesh::Update()
boundary_layer.SetSize(0);
}
/*
void ParNCMesh::AssignLeafIndices()
{
// This is an override of NCMesh::AssignLeafIndices(). The difference is
// that we shift all elements we own to the beginning of the array
// 'leaf_elements' and assign all ghost elements indices >= NElements.
// Also note that the ordering of ghosts and non-ghosts is preserved here,
// which is important for ParNCMesh::GetFaceNeighbors.
// We store the original leaf ordering in 'leaf_glob_order'. This is later
// used (and deleted) in GetConformingSharedStructures
NCMesh::AssignLeafIndices(); // original numbering, for 'leaf_glob_order'
int nleafs = leaf_elements.Size();
Array<int> ghosts;
ghosts.Reserve(nleafs);
NElements = 0;
for (int i = 0; i < nleafs; i++)
{
int elem = leaf_elements[i];
if (elements[elem].rank == MyRank)
{
leaf_elements[NElements++] = elem;
}
else
{
ghosts.Append(elem);
}
}
NGhostElements = ghosts.Size();
leaf_elements.SetSize(NElements);
leaf_elements.Append(ghosts);
// store original (globally consistent) numbering in 'leaf_glob_order'
leaf_glob_order.SetSize(nleafs);
for (int i = 0; i < nleafs; i++)
{
leaf_glob_order[i] = elements[leaf_elements[i]].index;
}
// new numbering with ghost shifted to the back
NCMesh::AssignLeafIndices();
}
void ParNCMesh::UpdateVertices()
{
// This is an override of NCMesh::UpdateVertices. This version first
// assigns vert_index to vertices of elements of our rank. Only these
// vertices then make it to the Mesh in NCMesh::GetMeshComponents.
// The remaining (ghost) vertices are assigned indices greater or equal to
// Mesh::GetNV().
for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
{
if (node->HasVertex()) { node->vert_index = -1; }
}
NVertices = 0;
for (int i = 0; i < leaf_elements.Size(); i++)
{
Element &el = elements[leaf_elements[i]];
if (el.rank == MyRank)
{
for (int j = 0; j < GI[el.Geom()].nv; j++)
{
int &vindex = nodes[el.node[j]].vert_index;
if (vindex < 0) { vindex = NVertices++; }
}
}
}
vertex_nodeId.SetSize(NVertices);
for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
{
if (node->HasVertex() && node->vert_index >= 0)
{
vertex_nodeId[node->vert_index] = node.index();
}
}
NGhostVertices = 0;
for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
{
if (node->HasVertex() && node->vert_index < 0)
{
node->vert_index = NVertices + (NGhostVertices++);
}
}
}
void ParNCMesh::OnMeshUpdated(Mesh *mesh)
{
std::cout << MyRank << ": Entering ParNCMesh::OnMeshUpdated" << std::endl;
// This is an override (or extension of) NCMesh::OnMeshUpdated().
// In addition to getting edge/face indices from 'mesh', we also
// assign indices to ghost edges/faces that don't exist in the 'mesh'.
// clear edge_index and Face::index
for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
{
if (node->HasEdge()) { node->edge_index = -1; }
}
for (face_iterator face = faces.begin(); face != faces.end(); ++face)
{
face->index = -1;
}
// go assign existing edge/face indices
NCMesh::OnMeshUpdated(mesh);
std::cout << MyRank << ": NVertices = " << NVertices << std::endl;
std::ostringstream ossN;
ossN << "node_on_mesh_updated_" << MyRank << ".out";
std::ofstream ofsN(ossN.str().c_str());
ofsN << nodes.Size() << std::endl;
for (int i=0; i<nodes.Size(); i++)
{
ofsN << i
// << " " << nodes[i].vert_refc
// << " " << nodes[i].edge_refc
<< " " << nodes[i].HasVertex()
<< " " << nodes[i].HasEdge()
<< " " << nodes[i].vert_index
<< " " << nodes[i].edge_index
<< " " << nodes[i].p1
<< " " << nodes[i].p2
<< " " << nodes[i].next << std::endl;
}
ofsN.close();
// count ghost edges and assign their indices
NEdges = mesh->GetNEdges();
NGhostEdges = 0;
for (node_iterator node = nodes.begin(); node != nodes.end(); ++node)
{
if (node->HasEdge() && node->edge_index < 0)
{
node->edge_index = NEdges + (NGhostEdges++);
}
}
// count ghost faces
NFaces = mesh->GetNumFaces();
NGhostFaces = 0;
for (face_iterator face = faces.begin(); face != faces.end(); ++face)
{
if (face->index < 0) { NGhostFaces++; }
}
if (Dim == 2)
{
// in 2D we have fake faces because of DG
MFEM_ASSERT(NFaces == NEdges, "");
MFEM_ASSERT(NGhostFaces == NGhostEdges, "");
}
// resize face_geom (default_geom is for slave faces beyond the ghost layer)
Geometry::Type default_geom = Geometry::SQUARE;
face_geom.SetSize(NFaces + NGhostFaces, default_geom);
// update 'face_geom' for ghost faces, assign ghost face indices
int nghosts = 0;
for (int i = 0; i < NGhostElements; i++)
{
Element &el = elements[leaf_elements[NElements + i]]; // ghost element
GeomInfo &gi = GI[el.Geom()];
for (int j = 0; j < gi.nf; j++)
{
const int *fv = gi.faces[j];
Face* face = faces.Find(el.node[fv[0]], el.node[fv[1]],
el.node[fv[2]], el.node[fv[3]]);
MFEM_ASSERT(face, "face not found!");
if (face->index < 0)
{
face->index = NFaces + (nghosts++);
// store the face geometry
static const Geometry::Type types[5] =
{
Geometry::INVALID, Geometry::INVALID,
Geometry::SEGMENT, Geometry::TRIANGLE, Geometry::SQUARE
};
face_geom[face->index] = types[gi.nfv[j]];
}
}
}
// assign valid indices also to faces beyond the ghost layer
for (face_iterator face = faces.begin(); face != faces.end(); ++face)
{
if (face->index < 0) { face->index = NFaces + (nghosts++); }
}
MFEM_ASSERT(nghosts == NGhostFaces, "");
{
/// Debugging output
std::ostringstream oss; oss << "elements_on_mesh_updated_"
<< MyRank << ".out";
std::ofstream ofs(oss.str().c_str());
for (int i=0; i<elements.Size(); i++)
{
ofs << i
<< '\t' << elements[i].index
<< '\t' << elements[i].rank
<< '\t' << elements[i].attribute
<< '\t' << elements[i].parent;
if ( elements[i].ref_type == 0 )
{
ofs << " nodes {";
for (int j=0; j<8; j++)
{
ofs << " " << elements[i].node[j];
}
ofs << "}";
}
else
{
ofs << " children {";
for (int j=0; j<8; j++)
{
ofs << " " << elements[i].child[j];
}
ofs << "}";
}
ofs << std::endl;
}
if (pncent_sets)
{
std::cout << "ParNCMesh::OnMeshUpdated pncent_sets is non NULL" << std::endl;
}
else
{
std::cout << "ParNCMesh::OnMeshUpdated pncent_sets is NULL" << std::endl;
}
if (ncent_sets)
{
std::cout << "ParNCMesh::OnMeshUpdated ncent_sets is non NULL" << std::endl;
}
else
{
std::cout << "ParNCMesh::OnMeshUpdated ncent_sets is NULL" << std::endl;
}
if (mesh->ent_sets)
{
std::cout << "ParNCMesh::OnMeshUpdated mesh->ent_sets is non NULL" << std::endl;
}
else
{
std::cout << "ParNCMesh::OnMeshUpdated mesh->ent_sets is NULL" << std::endl;
}
ParMesh * pmesh = dynamic_cast<ParMesh*>(mesh);
if (pmesh)
{
std::cout << "dynamic cast succeeded: mesh is a ParMesh" << std::endl;
if (pmesh->pent_sets != NULL)
{
std::cout << "ParNCMesh::OnMeshUpdated deleting ParEntitySets object in ParMesh"
<< std::endl;
delete pmesh->pent_sets;
}
else if (pmesh->ent_sets != NULL)
{
std::cout << "ParNCMesh::OnMeshUpdated deleting EntitySets object in ParMesh" <<
std::endl;
delete pmesh->ent_sets;
}
std::cout << "ParNCMesh::OnMeshUpdated creating ParEntitySets object in ParMesh"
<< std::endl;
pmesh->ent_sets = pmesh->pent_sets =
(pncent_sets) ? new ParEntitySets(*pmesh, *this): NULL;
*/
/*
if (pmesh->ent_sets)
{
std::cout << MyRank << ": ParNCMesh::OnMeshUpdated pmesh->ent_sets is non NULL" << std::endl;
pmesh->ent_sets->PrintSetInfo(std::cout);
std::ostringstream oss; oss << "ent_sets_" << MyRank << ".out";
std::ofstream ofs(oss.str().c_str());
pmesh->ent_sets->Print(ofs);
MPI_Barrier(MyComm);
std::cout << MyRank << ": testing " << NElements << std::endl;
//pmesh->ent_sets->Prune(NElements);
}
else
{
std::cout << "ParNCMesh::OnMeshUpdated pmesh->ent_sets is NULL" << std::endl;
}
*/
/*
if (pmesh->pent_sets)
{
std::cout << "ParNCMesh::OnMeshUpdated pmesh->pent_sets is non NULL" <<
std::endl;
}
else
{
std::cout << "ParNCMesh::OnMeshUpdated pmesh->pent_sets is NULL" << std::endl;
}
}
else
{
std::cout << "dynamic cast failed: mesh is not a ParMesh" << std::endl;
}
*/
/*
if (pncent_sets)
{
if (!pmesh->pent_sets)
{
pmesh->pent_sets = new ParEntitySets(*pmesh, *this);
}
}
*/
/*
// Prune the Entity Sets
if ( entity_sets )
{
EntitySets::EntityType t = EntitySets::INVALID;
unsigned int ns = -1;
std::cout << "Processing node sets" << std::endl;
t = EntitySets::VERTEX;
ns = entity_sets->GetNumSets(t);
for (unsigned int s=0; s<ns; s++)
{
unsigned int ni = entity_sets->GetNumEntities(t, s);
int e = 0;
for (unsigned int i=0; i<ni; i++)
{
if ( (*mesh->ent_sets)(t, s, i) < NVertices )
{
(*mesh->ent_sets)(t, s, e) = (*mesh->ent_sets)(t, s, i);
e++;
}
}
(*mesh->ent_sets)(t, s).resize(e);
}
t = EntitySets::EDGE;
ns = entity_sets->GetNumSets(t);
for (unsigned int s=0; s<ns; s++)
{
unsigned int ni = entity_sets->GetNumEntities(t, s);
BlockArray<int> ids;
for (unsigned int i=0; i<ni; i++)
{
if ( (*mesh->ent_sets)(t, s, i) < NEdges )
{
ids.Append((*mesh->ent_sets)(t, s, i));
}
}
(*mesh->ent_sets)(t, s).resize(ids.Size());
for (int i=0; i<ids.Size(); i++)
{
(*mesh->ent_sets)(t, s, i) = ids[i];
}
}
}
*/
/*
std::cout << MyRank << ": Leaving ParNCMesh::OnMeshUpdated" << std::endl;
}
}
*/
void ParNCMesh::ElementSharesFace(int elem, int local, int face)
{
// Analogous to ElementSharesEdge.
@@ -3152,94 +2732,6 @@ void ParNCMesh::GetDebugMesh(Mesh &debug_mesh) const
debug_mesh.ncmesh = copy;
}
void ParNCMesh::GetRefinedEdges(int vn0, int vn1, BlockArray<int> & edges)
{
std::cout << MyRank
<< ": entering ParNCMesh::GetRefinedEdges "
<<"searching for edge with vertices: " << vn0 << " and " << vn1
<< std::endl;
return this->NCMesh::GetRefinedEdges(vn0, vn1, edges);
int mid = nodes.FindId(vn0, vn1);
if (mid < 0) { return; }
/*
Node &nd = nodes[mid];
if ( nd.edge_index < 0 ) { return; }
edges.Append(nd.edge_index);
GetRefinedEdges(vn0, mid, edges);
GetRefinedEdges(mid, vn1, edges);
*/
edges.Append(mid);
GetRefinedEdges(vn0, mid, edges);
GetRefinedEdges(mid, vn1, edges);
}
void ParNCMesh::GetRefinedFaces(int vn0, int vn1, int vn2, int vn3,
BlockArray<int> & face_ids)
{
return this->NCMesh::GetRefinedFaces(vn0, vn1, vn2, vn3, face_ids);
/*
Face* fa = faces.Find(vn0, vn1, vn2, vn3);
if (fa)
{
if ( fa->index >= 0 )
{
face_ids.Append(fa->index);
}
return;
}
// we need to recurse deeper
int mid[4];
int split = FaceSplitType(vn0, vn1, vn2, vn3, mid);
if (split == 1) // "X" split face
{
GetRefinedFaces(vn0, mid[0], mid[2], vn3, face_ids);
GetRefinedFaces(mid[0], vn1, vn2, mid[2], face_ids);
}
else if (split == 2) // "Y" split face
{
GetRefinedFaces(vn0, vn1, mid[1], mid[3], face_ids);
GetRefinedFaces(mid[3], mid[1], vn2, vn3, face_ids);
}
*/
}
void ParNCMesh::GetRefinedElements(int elem_id, BlockArray<int> & elem_ids)
{
// std::cout << MyRank
// << ": entering ParNCMesh::GetRefinedElements "
// <<"searching for element id: " << elem_id << std::endl;
Element &el = elements[elem_id];
if (el.ref_type != 0)
{
// This element has been refined so recurse into its children
for (int i = 0; i < 8; i++)
{
if (el.child[i] >= 0 && el.child[i] < elements.Size() )
{
GetRefinedElements(el.child[i], elem_ids);
}
}
}
else
{
// This element has not been refined so add it if it's a local element
if (el.rank == MyRank)
{
elem_ids.Append(elem_id);
}
}
}
void ParNCMesh::Trim()
{
NCMesh::Trim();
-23
View File
@@ -20,7 +20,6 @@
#include <set>
#include "ncmesh.hpp"
#include "pentsets.hpp"
#include "../general/communication.hpp"
#include "../general/sort_pairs.hpp"
@@ -249,29 +248,9 @@ public:
The debug mesh will have element attributes set to element rank + 1. */
void GetDebugMesh(Mesh &debug_mesh) const;
/** Collect edge indices of all refined edges which are children of
the coarse edge defined by the given vertices. This method
overrides a method in NCMesh and only returns locally owned
edges. */
void GetRefinedEdges(int vn0, int vn1, BlockArray<int> & edge_ids);
/** Collect face indices of all refined faces which are children of
the coarse face defined by the given vertices. This method
overrides a method in NCMesh and only returns locally owned
faces. */
void GetRefinedFaces(int vn0, int vn1, int vn2, int vn3,
BlockArray<int> & face_ids);
/** Collect element indices of all refined elements which are
children of the coarse element defined by the given element
index. This method overrides a method in NCMesh and only
returns locally owned elements. */
void GetRefinedElements(int elem_id, BlockArray<int> & elem_ids);
protected: // interface for ParMesh
friend class ParMesh;
friend class ParEntitySets;
/** For compatibility with conforming code in ParMesh and ParFESpace.
Initializes shared structures in ParMesh: gtopo, shared_*, group_s*, s*_l*.
@@ -561,8 +540,6 @@ protected: // implementation
Array<DenseMatrix*> aux_pm_store;
void ClearAuxPM();
ParNCEntitySets * pncent_sets;
long GroupsMemoryUsage() const;
friend class NeighborRowMessage;
+2 -2
View File
@@ -804,11 +804,11 @@ public:
mfem::ad::ADFloatType aa;
mfem::ad::ADVectorType lvec(ndof);
for (int i = 0; i < ir.GetNPoints(); i++)
for (int q = 0; q < ir.GetNPoints(); q++)
{
lvec=0.0;
const IntegrationPoint &ip = ir.IntPoint(i);
const IntegrationPoint &ip = ir.IntPoint(q);
trans.SetIntPoint(&ip);
w = trans.Weight();
detJ = (square ? w : w * w);
+2 -2
View File
@@ -485,7 +485,7 @@ int main(int argc, char *argv[])
<< std::endl;
}
// Compute the energy
double energy = nr->GetEnergy(*sv);
energy = nr->GetEnergy(*sv);
if (myrank==0)
{
std::cout << "[pp="<<i<<"] The total energy of the system is E=" << energy
@@ -517,7 +517,7 @@ int main(int argc, char *argv[])
<< std::endl;
}
// Compute the energy
double energy = nr->GetEnergy(*sv);
energy = nr->GetEnergy(*sv);
if (myrank==0)
{
std::cout << "[pp="<<pp<<"] The total energy of the system is E=" << energy
+13 -13
View File
@@ -678,22 +678,22 @@ void MagneticDiffusionEOperator::buildA1(double muInv,
dt_A1 = dt;
}
void MagneticDiffusionEOperator::buildA2(MeshDependentCoefficient &InvTcond,
MeshDependentCoefficient &InvTcap,
double dt)
void MagneticDiffusionEOperator::buildA2(MeshDependentCoefficient &InvTcond_,
MeshDependentCoefficient &InvTcap_,
double dt_)
{
if ( a2 != NULL ) { delete a2; }
InvTcap.SetScaleFactor(dt);
InvTcap_.SetScaleFactor(dt_);
a2 = new ParBilinearForm(&HDivFESpace);
a2->AddDomainIntegrator(new VectorFEMassIntegrator(InvTcond));
a2->AddDomainIntegrator(new DivDivIntegrator(InvTcap));
a2->AddDomainIntegrator(new VectorFEMassIntegrator(InvTcond_));
a2->AddDomainIntegrator(new DivDivIntegrator(InvTcap_));
if (STATIC_COND == 1) { a2->EnableStaticCondensation(); }
a2->Assemble();
// Don't finalize or parallel assemble this is done in FormLinearSystem.
dt_A2 = dt;
dt_A2 = dt_;
}
void MagneticDiffusionEOperator::buildM1(MeshDependentCoefficient &Sigma)
@@ -719,13 +719,13 @@ void MagneticDiffusionEOperator::buildM2(MeshDependentCoefficient &Alpha)
// Don't finalize or parallel assemble this is done in FormLinearSystem.
}
void MagneticDiffusionEOperator::buildM3(MeshDependentCoefficient &Tcapacity)
void MagneticDiffusionEOperator::buildM3(MeshDependentCoefficient &Tcapacity_)
{
if ( m3 != NULL ) { delete m3; }
// ConstantCoefficient Sigma(sigma);
m3 = new ParBilinearForm(&L2FESpace);
m3->AddDomainIntegrator(new MassIntegrator(Tcapacity));
m3->AddDomainIntegrator(new MassIntegrator(Tcapacity_));
m3->Assemble();
m3->Finalize();
M3 = m3->ParallelAssemble();
@@ -741,13 +741,13 @@ void MagneticDiffusionEOperator::buildS1(double muInv)
s1->Assemble();
}
void MagneticDiffusionEOperator::buildS2(MeshDependentCoefficient &InvTcap)
void MagneticDiffusionEOperator::buildS2(MeshDependentCoefficient &InvTcap_)
{
if ( s2 != NULL ) { delete s2; }
// ConstantCoefficient param(a);
s2 = new ParBilinearForm(&HDivFESpace);
s2->AddDomainIntegrator(new DivDivIntegrator(InvTcap));
s2->AddDomainIntegrator(new DivDivIntegrator(InvTcap_));
s2->Assemble();
}
@@ -768,13 +768,13 @@ void MagneticDiffusionEOperator::buildCurl(double muInv)
// no ParallelAssemble since this will be applied to GridFunctions
}
void MagneticDiffusionEOperator::buildDiv(MeshDependentCoefficient &InvTcap)
void MagneticDiffusionEOperator::buildDiv(MeshDependentCoefficient &InvTcap_)
{
if ( weakDiv != NULL ) { delete weakDiv; }
if ( weakDivC != NULL ) { delete weakDivC; }
weakDivC = new ParMixedBilinearForm(&HDivFESpace, &L2FESpace);
weakDivC->AddDomainIntegrator(new VectorFEDivergenceIntegrator(InvTcap));
weakDivC->AddDomainIntegrator(new VectorFEDivergenceIntegrator(InvTcap_));
weakDivC->Assemble();
weakDiv = new ParMixedBilinearForm(&HDivFESpace, &L2FESpace);
+4 -4
View File
@@ -1143,8 +1143,8 @@ int main (int argc, char *argv[])
if (mk == 'S')
{
const char mesh_file[] = "mesh-explorer.mesh";
ofstream omesh(mesh_file);
const char omesh_file[] = "mesh-explorer.mesh";
ofstream omesh(omesh_file);
omesh.precision(14);
mesh->Print(omesh);
cout << "New mesh file: " << mesh_file << endl;
@@ -1152,8 +1152,8 @@ int main (int argc, char *argv[])
if (mk == 'V')
{
const char mesh_file[] = "mesh-explorer.vtk";
ofstream omesh(mesh_file);
const char omesh_file[] = "mesh-explorer.vtk";
ofstream omesh(omesh_file);
omesh.precision(14);
mesh->PrintVTK(omesh);
cout << "New VTK mesh file: " << mesh_file << endl;

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