diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json index 79e76a9d24..46d12694bd 100644 --- a/.vscode/c_cpp_properties.json +++ b/.vscode/c_cpp_properties.json @@ -23,7 +23,7 @@ ], "cStandard": "c11", "cppStandard": "c++14", - "intelliSenseMode": "macos-clang-x64" + "intelliSenseMode": "macos-clang-arm64" }, { "name": "linux-gcc", diff --git a/CHANGES.md b/CHANGES.md index e2f725c72f..7514e18f1d 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -7,16 +7,19 @@ The version numbers are in the form `MAJOR.MINOR.PATCH`, where major releases in ## Version 4.2.0 +* **Changed** + * Updated `examples/strikeslip-2d` Steps 4-7 to use a more realistic slip distribution and mesh refinement in output. + * Updated to PETSc 3.22.0 + * Switch CI from Azure Pipelines to GitHub Actions. * **Added** * Default filenames for progress monitor and parameters file are set from the simulation name like the other output files. * Consistency check for material properties and scales used in nondimensionalization (currently just the shear modulus). * Add section to User Guide on troubleshooting solver issues. -* **Changed** - * Switch CI from Azure Pipelines to GitHub Actions. + * Allow output on a finer resolution mesh than used in the simulation to facilitate accurate representation of fields with a basis order of 2 or greater. * **Fixed** * Fix inconsistency in normal direction on fault surfaces. Orientation was correct but direction was flipped at some locations. This affected local slip direction and the resulting deformation close to the fault. This bug fix was not in version 4.1.3. - * Update autoconf macros for numpy for compatibility with location of include files in numpy version 2.x. + * Update autoconf numpy macros for compatibility with location of include files in numpy version 2.x. ## Version 4.1.3 (2024/07/31) diff --git a/docs/Makefile.am b/docs/Makefile.am index 02c7cc9ec0..e52421aba0 100644 --- a/docs/Makefile.am +++ b/docs/Makefile.am @@ -473,13 +473,15 @@ dist_noinst_DATA = \ user/examples/strikeslip-2d/figs/step04-solution.jpg \ user/examples/strikeslip-2d/figs/step05-diagram.pdf \ user/examples/strikeslip-2d/figs/step05-diagram.svg \ + user/examples/strikeslip-2d/figs/step05_greensfns-impulses.pdf \ + user/examples/strikeslip-2d/figs/step05_greensfns-impulses.svg \ user/examples/strikeslip-2d/figs/step05-solution.jpg \ - user/examples/strikeslip-2d/figs/step06-solution.pdf \ - user/examples/strikeslip-2d/figs/step06-solution.svg \ - user/examples/strikeslip-2d/figs/step07a-solution.pdf \ - user/examples/strikeslip-2d/figs/step07a-solution.svg \ - user/examples/strikeslip-2d/figs/step07b-solution.pdf \ - user/examples/strikeslip-2d/figs/step07b-solution.svg \ + user/examples/strikeslip-2d/figs/step06_inversion-results.pdf \ + user/examples/strikeslip-2d/figs/step06_inversion-results.svg \ + user/examples/strikeslip-2d/figs/step07a_catmip-results.pdf \ + user/examples/strikeslip-2d/figs/step07a_catmip-results.svg \ + user/examples/strikeslip-2d/figs/step07b_cfcatmip-results.pdf \ + user/examples/strikeslip-2d/figs/step07b_cfcatmip-results.svg \ user/examples/subduction-2d/index.md \ user/examples/subduction-2d/common-information.md \ user/examples/subduction-2d/exercises.md \ diff --git a/docs/user/components/meshio/OutputObserver.md b/docs/user/components/meshio/OutputObserver.md index 3010d8583d..5909c20867 100644 --- a/docs/user/components/meshio/OutputObserver.md +++ b/docs/user/components/meshio/OutputObserver.md @@ -21,4 +21,8 @@ Abstract base class for managing output of solution information. - **default value**: 1 - **current value**: 1, from {default} - **validator**: (in [0, 1]) +* `refine_levels`=\: Number of mesh refinement levels for output. + - **default value**: 0 + - **current value**: 0, from {default} + - **validator**: (greater than or equal to 0) diff --git a/docs/user/components/meshio/OutputPhysics.md b/docs/user/components/meshio/OutputPhysics.md index dceacafa6e..55073b734e 100644 --- a/docs/user/components/meshio/OutputPhysics.md +++ b/docs/user/components/meshio/OutputPhysics.md @@ -33,6 +33,10 @@ Implements `OutputObserver`. - **default value**: 1 - **current value**: 1, from {default} - **validator**: (in [0, 1]) +* `refine_levels`=\: Number of mesh refinement levels for output. + - **default value**: 0 + - **current value**: 0, from {default} + - **validator**: (greater than or equal to 0) ## Example diff --git a/docs/user/components/meshio/OutputSoln.md b/docs/user/components/meshio/OutputSoln.md index 7527d9e769..0021d71aa7 100644 --- a/docs/user/components/meshio/OutputSoln.md +++ b/docs/user/components/meshio/OutputSoln.md @@ -26,4 +26,8 @@ Implements `OutputObserver`. - **default value**: 1 - **current value**: 1, from {default} - **validator**: (in [0, 1]) +* `refine_levels`=\: Number of mesh refinement levels for output. + - **default value**: 0 + - **current value**: 0, from {default} + - **validator**: (greater than or equal to 0) diff --git a/docs/user/components/meshio/OutputSolnBoundary.md b/docs/user/components/meshio/OutputSolnBoundary.md index b707c66a2a..f47f7ee5c0 100644 --- a/docs/user/components/meshio/OutputSolnBoundary.md +++ b/docs/user/components/meshio/OutputSolnBoundary.md @@ -29,7 +29,7 @@ Implements `OutputSoln`. * `label`=\: Name of label identifier for external boundary. - **default value**: '' - **current value**: '', from {default} - - **validator**: + - **validator**: * `label_value`=\: Value of label identifier for external boundary (tag of physical group in Gmsh files). - **default value**: 1 - **current value**: 1, from {default} @@ -37,6 +37,10 @@ Implements `OutputSoln`. - **default value**: 1 - **current value**: 1, from {default} - **validator**: (in [0, 1]) +* `refine_levels`=\: Number of mesh refinement levels for output. + - **default value**: 0 + - **current value**: 0, from {default} + - **validator**: (greater than or equal to 0) ## Example diff --git a/docs/user/components/meshio/OutputSolnDomain.md b/docs/user/components/meshio/OutputSolnDomain.md index 7cbbe2ad11..d10e3828c6 100644 --- a/docs/user/components/meshio/OutputSolnDomain.md +++ b/docs/user/components/meshio/OutputSolnDomain.md @@ -30,6 +30,10 @@ Implements `OutputSoln`. - **default value**: 1 - **current value**: 1, from {default} - **validator**: (in [0, 1]) +* `refine_levels`=\: Number of mesh refinement levels for output. + - **default value**: 0 + - **current value**: 0, from {default} + - **validator**: (greater than or equal to 0) ## Example diff --git a/docs/user/components/meshio/OutputSolnPoints.md b/docs/user/components/meshio/OutputSolnPoints.md index 3beca50750..9588ace460 100644 --- a/docs/user/components/meshio/OutputSolnPoints.md +++ b/docs/user/components/meshio/OutputSolnPoints.md @@ -36,6 +36,10 @@ Implements `OutputSoln`. - **default value**: 1 - **current value**: 1, from {default} - **validator**: (in [0, 1]) +* `refine_levels`=\: Number of mesh refinement levels for output. + - **default value**: 0 + - **current value**: 0, from {default} + - **validator**: (greater than or equal to 0) ## Example diff --git a/docs/user/examples/strikeslip-2d/figs/step04-slip.pdf b/docs/user/examples/strikeslip-2d/figs/step04-slip.pdf index 6609b620b9..36c3d5deec 100644 Binary files a/docs/user/examples/strikeslip-2d/figs/step04-slip.pdf and b/docs/user/examples/strikeslip-2d/figs/step04-slip.pdf differ diff --git a/docs/user/examples/strikeslip-2d/figs/step04-slip.svg b/docs/user/examples/strikeslip-2d/figs/step04-slip.svg index 950d98c193..b81ec524f3 100644 --- a/docs/user/examples/strikeslip-2d/figs/step04-slip.svg +++ b/docs/user/examples/strikeslip-2d/figs/step04-slip.svg @@ -1,236 +1,224 @@ - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/user/examples/strikeslip-2d/figs/step07a-solution.pdf b/docs/user/examples/strikeslip-2d/figs/step05_greensfns-impulses.pdf similarity index 57% rename from docs/user/examples/strikeslip-2d/figs/step07a-solution.pdf rename to docs/user/examples/strikeslip-2d/figs/step05_greensfns-impulses.pdf index 1223f609fa..6723ec1469 100644 Binary files a/docs/user/examples/strikeslip-2d/figs/step07a-solution.pdf and b/docs/user/examples/strikeslip-2d/figs/step05_greensfns-impulses.pdf differ diff --git a/docs/user/examples/strikeslip-2d/figs/step05_greensfns-impulses.svg b/docs/user/examples/strikeslip-2d/figs/step05_greensfns-impulses.svg new file mode 100644 index 0000000000..1f33602be2 --- /dev/null +++ b/docs/user/examples/strikeslip-2d/figs/step05_greensfns-impulses.svg @@ -0,0 +1,286 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/user/examples/strikeslip-2d/figs/step06-solution.svg b/docs/user/examples/strikeslip-2d/figs/step06-solution.svg deleted file mode 100644 index 906e7383cb..0000000000 --- a/docs/user/examples/strikeslip-2d/figs/step06-solution.svg +++ /dev/null @@ -1,1285 +0,0 @@ - - - - - - - - 2024-06-06T05:30:41.253651 - image/svg+xml - - - Matplotlib v3.8.2, https://matplotlib.org/ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/user/examples/strikeslip-2d/figs/step06-solution.pdf b/docs/user/examples/strikeslip-2d/figs/step06_inversion-results.pdf similarity index 55% rename from docs/user/examples/strikeslip-2d/figs/step06-solution.pdf rename to docs/user/examples/strikeslip-2d/figs/step06_inversion-results.pdf index bf40381433..3c901af6d8 100644 Binary files a/docs/user/examples/strikeslip-2d/figs/step06-solution.pdf and b/docs/user/examples/strikeslip-2d/figs/step06_inversion-results.pdf differ diff --git a/docs/user/examples/strikeslip-2d/figs/step06_inversion-results.svg b/docs/user/examples/strikeslip-2d/figs/step06_inversion-results.svg new file mode 100644 index 0000000000..e95411c9c3 --- /dev/null +++ b/docs/user/examples/strikeslip-2d/figs/step06_inversion-results.svg @@ -0,0 +1,336 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/user/examples/strikeslip-2d/figs/step07a-solution.svg b/docs/user/examples/strikeslip-2d/figs/step07a-solution.svg deleted file mode 100644 index 71104fae95..0000000000 --- a/docs/user/examples/strikeslip-2d/figs/step07a-solution.svg +++ /dev/null @@ -1,1293 +0,0 @@ - - - - - - - - 2023-06-09T21:21:51.387141 - image/svg+xml - - - Matplotlib v3.7.1, https://matplotlib.org/ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/user/examples/strikeslip-2d/figs/step07b-solution.pdf b/docs/user/examples/strikeslip-2d/figs/step07a_catmip-results.pdf similarity index 54% rename from docs/user/examples/strikeslip-2d/figs/step07b-solution.pdf rename to docs/user/examples/strikeslip-2d/figs/step07a_catmip-results.pdf index 85af0efc20..f465d6754b 100644 Binary files a/docs/user/examples/strikeslip-2d/figs/step07b-solution.pdf and b/docs/user/examples/strikeslip-2d/figs/step07a_catmip-results.pdf differ diff --git a/docs/user/examples/strikeslip-2d/figs/step07a_catmip-results.svg b/docs/user/examples/strikeslip-2d/figs/step07a_catmip-results.svg new file mode 100644 index 0000000000..2b1e9ccef6 --- /dev/null +++ b/docs/user/examples/strikeslip-2d/figs/step07a_catmip-results.svg @@ -0,0 +1,322 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/user/examples/strikeslip-2d/figs/step07b-solution.svg b/docs/user/examples/strikeslip-2d/figs/step07b-solution.svg deleted file mode 100644 index df18bc2978..0000000000 --- a/docs/user/examples/strikeslip-2d/figs/step07b-solution.svg +++ /dev/null @@ -1,1293 +0,0 @@ - - - - - - - - 2023-06-09T21:22:00.438095 - image/svg+xml - - - Matplotlib v3.7.1, https://matplotlib.org/ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/user/examples/strikeslip-2d/figs/step07b_cfcatmip-results.pdf b/docs/user/examples/strikeslip-2d/figs/step07b_cfcatmip-results.pdf new file mode 100644 index 0000000000..95ed95d2dc Binary files /dev/null and b/docs/user/examples/strikeslip-2d/figs/step07b_cfcatmip-results.pdf differ diff --git a/docs/user/examples/strikeslip-2d/figs/step07b_cfcatmip-results.svg b/docs/user/examples/strikeslip-2d/figs/step07b_cfcatmip-results.svg new file mode 100644 index 0000000000..b393040559 --- /dev/null +++ b/docs/user/examples/strikeslip-2d/figs/step07b_cfcatmip-results.svg @@ -0,0 +1,322 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/user/examples/strikeslip-2d/step04-varslip.md b/docs/user/examples/strikeslip-2d/step04-varslip.md index abb340a14d..1119f0cb34 100644 --- a/docs/user/examples/strikeslip-2d/step04-varslip.md +++ b/docs/user/examples/strikeslip-2d/step04-varslip.md @@ -24,7 +24,7 @@ Boundary conditions for static coseismic slip. We set the x and y displacement to zero on the +x and -x boundaries and prescribe left-lateral slip that varies along strike. ::: -For greater accuracy in modeling the spatial variation in slip, we use a basis order of 2 for the solution subfields. +For greater accuracy in modeling the spatial variation in slip, we refine the mesh by a factor of 2 and use a basis order of 2 for the solution subfields. ```{code-block} cfg --- @@ -93,12 +93,15 @@ Prescribed left-lateral slip that varies along the strike of the fault. A strike of 0 corresponds to y=0. ::: -We use a `SimpleDB` to define the spatial variation in slip. +We use a `SimpleGridDB` to define the spatial variation in slip. ```{code-block} cfg --- -caption: Prescribed slip parameters for Step 4. +caption: Prescribed slip parameters for Step 4. We refine the fault mesh by a factor of 8 (3 levels of refinement by a factor of 2) so that the output, which uses a basis order of 1, better captures the discretization of slip, which uses a basis order of 2. --- +[pylithapp.problem.interfaces.fault] +observers.observer.refine_levels = 3 + [pylithapp.problem.interfaces.fault.eq_ruptures.rupture] db_auxiliary_field = spatialdata.spatialdb.SimpleDB db_auxiliary_field.description = Fault rupture auxiliary field spatial database diff --git a/docs/user/examples/strikeslip-2d/step05-greensfns.md b/docs/user/examples/strikeslip-2d/step05-greensfns.md index 239ff0a66e..a050fccf03 100644 --- a/docs/user/examples/strikeslip-2d/step05-greensfns.md +++ b/docs/user/examples/strikeslip-2d/step05-greensfns.md @@ -25,14 +25,12 @@ In the fault interfaces section we set the fault type to `FaultCohesiveImpulses` We also use a spatial database to limit the section of the fault where we impose the fault slip impulses to -25 km $\le$ y $\le$ +25 km. :::{important} -**Currently, a basis order of 1 (default) for the slip auxiliary subfield is the only choice that gives accurate results in a slip inversion due to the factors described here.** +**Currently, we recommend a basis order of 1 (default) for the slip auxiliary subfield.** The basis order for the slip auxiliary subfield controls the representation of the slip field for the impulses. For a given impulse, a basis order of 1 will impose unit slip at a vertex with zero slip at all other vertices. -Likewise, a basis order of 0 will attempt to impose unit slip over a cell with zero slip in all other cells; however, this creates a jump in slip at the cell boundaries that cannot be accurately represented by the finite-element solution. -As a result, you should not use a basis order of 0 for the slip auxiliary field. -A basis order of 2 will impose slip at vertices as well as edge degrees of freedom in the cell. -Because PyLith output decimates the basis order to 0 or 1, you should avoid this choice of basis order as well until we provide better ways to output fields discretized with higher order basis functions. +A basis order of 0 will impose unit slip over a cell with zero slip in all other cells; however, this creates jumps in the slip at the cell boundaries and reduces the accuracy of the slip inversion for a given mesh resolution. +A basis order of 2 will impose slip at vertices as well as edge degrees of freedom in the cell; regularization becomes trickier with the higher order discretization, and a basis order of 2 does not always give more accurate results for a given number of impulses. ::: ```{code-block} cfg @@ -112,16 +110,24 @@ We get warnings about unused PETSc options because we do not use time stepping. ## Visualizing the results +In {numref}`fig:example:strikeslip:2d:step05:impulses` we use the `viz/plot_slip_impulses.py` Pyhon script to visualize the slip impulses. In {numref}`fig:example:strikeslip:2d:step05:solution` we use the `pylith_viz` utility to visualize the y displacement field. You can move the slider or use the `p` and `n` keys to change the increment or decrement slip impulses (shown as different time stamps). ```{code-block} console --- -caption: Visualize PyLith output using `pylith_viz`. +caption: Plot the slip impulses and visualize PyLith output using `pylith_viz`. --- +viz/plot_slip_impulses.py pylith_viz --filename=output/step05_greensfns-domain.h5 warp_grid --component=y ``` +:::{figure-md} fig:example:strikeslip:2d:step05:impulses +Slip impulses for the Green's functions in Step 5. + +Slip impulses for the Green's functions in Step 5. +::: + :::{figure-md} fig:example:strikeslip:2d:step05:solution Solution for Step 5. The colors indicate the y displacement, and the deformation is exaggerated by a factor of 1000. diff --git a/docs/user/examples/strikeslip-2d/step06-inversion-leastsquares.md b/docs/user/examples/strikeslip-2d/step06-inversion-leastsquares.md index 8fdd519ae3..5aabd81916 100644 --- a/docs/user/examples/strikeslip-2d/step06-inversion-leastsquares.md +++ b/docs/user/examples/strikeslip-2d/step06-inversion-leastsquares.md @@ -16,7 +16,7 @@ $ ./invert_slip.py $ ./invert_slip.py --help ``` -By default, the inversion code will write the results of the inversion to `output/step06_greensfns-inversion_results.txt`. +By default, the inversion code will write the results of the inversion to `output/step06_inversion-results.txt`. ## Plotting the results @@ -26,18 +26,18 @@ If you are using the PyLith binary, which includes the `matplotlib` Python modul --- caption: Plot the inversion results using the `matplotlib` Python module. --- -$ ./viz/plot_inversion_results.py +$ viz/plot_inversion_results.py ``` :::{figure-md} fig:example:strikeslip:2d:step06:solution -Results of slip inversion in Step 6. +Results of slip inversion in Step 6. Results of slip inversion in Step 6. -The thick black line shows the prescribed slip in Step 4. +The thick black line shows the prescribed slip in Step 4, and the thick gray line shows the raw (exact) slip. The thin colored lines show the slip from the inversion with different penalty factors. ::: ::{tip} -You can pass`--no-gui` as a command line argument to the plotting script turn off displaying the plot window. +You can pass`--no-gui` as a command line argument to the plotting script to turn off displaying the plot window. This is useful if you do not have a matplotlib GUI backend. ::: diff --git a/docs/user/examples/strikeslip-2d/step07-inversion-catmip.md b/docs/user/examples/strikeslip-2d/step07-inversion-catmip.md index c242e813ed..9a5edf6ee6 100644 --- a/docs/user/examples/strikeslip-2d/step07-inversion-catmip.md +++ b/docs/user/examples/strikeslip-2d/step07-inversion-catmip.md @@ -74,10 +74,10 @@ viz/plot_catmip_results.py --catmip-theta=output/step07a_catmip-theta20.bin ``` :::{figure-md} fig:example:strikeslip:2d:step07a:solution -Results of slip inversion in Step 7a. +Results of slip inversion in Step 7a. Results of slip inversion in Step 7a. -The thick black line shows the prescribed slip in Step 4. +The thick black line shows the prescribed slip in Step 4, and the thick gray line shows the raw (exact) slip. The solid orange line shows the median slip, the shaded orange regions shows the median plus and minus one standard deviation, and the dashed lines show the minimum and maximum values. The median values are almost identical to the results with minimal smoothing in Step 6. ::: @@ -157,10 +157,10 @@ viz/plot_catmip_results.py --catmip-theta=output/step07b_cfcatmip-theta1.bin ``` :::{figure-md} fig:example:strikeslip:2d:step07b:solution -Results of slip inversion in Step 7b. +Results of slip inversion in Step 7b. Results of slip inversion in Step 7b. -The thick black line shows the prescribed slip in Step 4. +The thick black line shows the prescribed slip in Step 4, and the thick gray line shows the raw (exact) slip. The solid orange line shows the median slip, the shaded orange regions shows the median plus and minus one standard deviation, and the dashed lines show the minimum and maximum values. The results are almost identical to those in Step 7a. ::: diff --git a/docs/user/problems/output.md b/docs/user/problems/output.md index 15bb50402c..06517920ef 100644 --- a/docs/user/problems/output.md +++ b/docs/user/problems/output.md @@ -79,6 +79,13 @@ a_z &= a_n n_z + a_h h_z + a_v v_z This transformation is useful for plotting fault tractions and slip in 3D visualization tools such as ParaView that require vectors in the model coordinate system. +*New in v4.2.0* + +Output observers can output information on a finer resolution mesh than the one used in the simulaiton. +This is useful when using a basis order of 2 or greater for the discretization of the solution and writing output at a basis order of 1. +In this case, output using a mesh that is 4-8 times finer and fields discretized using a basis order of 1 (usually required by visualization tools) will do a good job representing the fields with a basis order of 2 on the mesh used in the simulation. +This feature is used in [`examples/strikeslip-2d/step04-varslip`](../examples/strikeslip-2d/step04-varslip.md). + (sec-user-solution-observers)= ### Solution Observers diff --git a/examples/strikeslip-2d/Makefile.am b/examples/strikeslip-2d/Makefile.am index 12c567f400..b62d9e970f 100644 --- a/examples/strikeslip-2d/Makefile.am +++ b/examples/strikeslip-2d/Makefile.am @@ -12,7 +12,12 @@ dist_noinst_PYTHON = \ generate_gmsh.py \ generate_cubit.py \ generate_gnssstations.py \ - invert_slip.py + generate_slip.py \ + invert_slip.py \ + viz/plot_slip.py \ + viz/plot_slip_impulses.py \ + viz/plot_inversion_results.py \ + viz/plot_catmip_results.py dist_noinst_DATA = \ README.md \ @@ -31,9 +36,7 @@ dist_noinst_DATA = \ disprate_bc_xpos.spatialdb \ gnss_stations.txt \ slip_impulses.spatialdb \ - slip_variable.spatialdb \ - viz/plot_inversion_results.py \ - viz/plot_catmip_results.py + slip_variable.spatialdb # End of file diff --git a/examples/strikeslip-2d/cfcatmip_parameters.txt b/examples/strikeslip-2d/cfcatmip_parameters.txt index 43a5bd85e0..26ff163eaa 100644 --- a/examples/strikeslip-2d/cfcatmip_parameters.txt +++ b/examples/strikeslip-2d/cfcatmip_parameters.txt @@ -1,8 +1,8 @@ -# Sample parameter file for pylith_catmip model. +# Sample parameter file for pylith_cfcatmip model. # Data files -filename_observations = output/step04_varslip-gps_stations.h5 -filename_greens_fns = output/step05_greensfns-gps_stations.h5 +filename_observations = output/step04_varslip-gnss_stations.h5 +filename_greens_fns = output/step05_greensfns-gnss_stations.h5 # Model parameters rake_parallel_prior = logistic diff --git a/examples/strikeslip-2d/generate_slip.py b/examples/strikeslip-2d/generate_slip.py new file mode 100755 index 0000000000..91ba7e5569 --- /dev/null +++ b/examples/strikeslip-2d/generate_slip.py @@ -0,0 +1,362 @@ +#!/usr/bin/env nemesis +"""Python script for generating slip profile. + +REQUIRES: scipy (not included in PyLith binary distribution) +""" + +import math +import abc + +import scipy +import numpy + +from pythia.pyre.units.length import m, km + + +def cli(): + """Command line interface. + """ + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument("--eq-magnitude", action="store", type=float, dest="eq_magnitude", default=6.5, help="Earthquake magnitude.") + parser.add_argument("--total-length", action="store", type=float, dest="total_length", default=40.0, help="Rupture length (km).") + parser.add_argument("--taper-length", action="store", type=float, dest="taper_length", default=0.5, help="Taper length (km).") + parser.add_argument("--resolution", action="store", type=float, dest="resolution", default=100.0, help="Resolution (m).") + parser.add_argument("--db-filename", action="store", dest="filename_db", default="slip_variable.spatialdb", help="Filename for slip spatial database.") + parser.add_argument("--raw-filename", action="store", dest="filename_raw", default="output/slip_variable.txt", help="Filename for raw slip data.") + args = parser.parse_args() + + generator = SlipGenerator( + eq_magnitude=args.eq_magnitude, + total_length=args.total_length*km, + taper_length=args.taper_length*km, + dx=args.resolution*m, + ) + (x, slip) = generator.generate() + SpatialDBWriter.write(args.filename_db, x, slip) + RawWriter.write(args.filename_raw, x, slip) + + +class SlipGenerator: + """Object for generating 1D slip profile. + """ + + def __init__(self, eq_magnitude: float, total_length: float, taper_length: float, dx: float): + """Generate 1D slip profile. + + Args: + eq_magnitude: Earthquake magnitude. + total_length: Rupture length (m) + taper_length: Length of taper (m) at ends of rupture. + dx: Discretization size for slip profile. + """ + self.eq_magnitude = eq_magnitude + self.total_length = total_length + self.taper_length = taper_length + self.dx = dx + + def generate(self) -> tuple: + """Generate slip profile. + + Returns: + (numpy.ndarray, numpy.ndarray): (distance along fault (m), slip along fault (m)) + """ + MAX_ITERATIONS = 100 + END_TOLERANCE = 0.02 + ZERO_TOLERANCE = 0.05 + SPLINE_SPACING = 4 # Ratio of discretization size in raw slip profile to discretization size in spline + + dx_points = SPLINE_SPACING*self.dx # Discretization size of raw slip profile + randomizer = RandomizerLA( + eq_magnitude=self.eq_magnitude, + total_length=self.total_length, + dx=dx_points, + ) + total_length_m = self.total_length.value + dx_points_m = dx_points.value + x_pts = numpy.arange(-0.5*total_length_m, +0.5*(total_length_m+dx_points_m), dx_points_m) + for i in range(MAX_ITERATIONS): + # Iterate until we find slip profile that is nearly zero at the ends + slip_pts = randomizer.create(x_pts) + if slip_pts[0] < END_TOLERANCE and numpy.sum(slip_pts == 0.0) < (1.0-ZERO_TOLERANCE)*len(slip_pts): + break + else: + print("Failed to generated slip distribution.") + + # Apply taper to force slip profile to zero at the ends. + taper = Taper( + total_length=self.total_length, + taper_length=self.taper_length, + ) + taper.apply(x_pts, slip_pts) + + # Fit spline to slip profile + dx_spline_m = self.dx.value + x = numpy.arange(-0.5*total_length_m, +0.5*(total_length_m+dx_spline_m), dx_spline_m) + spline = scipy.interpolate.CubicSpline(x_pts, slip_pts) + slip_spline = spline(x) + slip_spline[0] = 0.0 + slip_spline[-1] = 0.0 + slip_spline[slip_spline < 0.0] = 0.0 + return (x, slip_spline) + + +class SpatialDBWriter: + """Object for writing rupture information to spatial database file. + """ + + @staticmethod + def write(filename: str, y: numpy.ndarray, slip: numpy.ndarray): + """Write rupture information to SimpleGridDB ASCII file. + + Args: + filename: Name of spatial database file. + y: y coordinate (m) + slip: Slip (m) along fault + """ + from spatialdata.spatialdb.SimpleGridAscii import SimpleGridAscii + from spatialdata.geocoords.CSCart import CSCart + + DOMAIN_Y = 150.0*km + + writer = SimpleGridAscii() + writer.filename = filename + + # Add points to encompass entire fault (regions beyond rupture) + y_ext = numpy.concatenate(([-0.5*DOMAIN_Y.value], y, [+0.5*DOMAIN_Y.value])) + slip_ext = numpy.concatenate(([0.0], slip, [0.0])) + + cs = CSCart() + cs.spaceDim = 2 + cs._configure() + points = numpy.zeros((len(y_ext), 2)) + points[:,1] = y_ext + + data = { + "points": points, + 'x': numpy.array([0]), + 'y': y_ext, + 'z': numpy.array([0]), + "coordsys": cs, + "data_dim": 1, + "values": [{ + "name": "final_slip_left_lateral", + "units": "m", + "data": slip_ext, + },{ + "name": "final_slip_opening", + "units": "m", + "data": 0*slip_ext, + },{ + "name": "initiation_time", + "units": "s", + "data": 0*slip_ext, + }] + } + writer.write(data) + + +class RawWriter: + """Write slip profile to ASCII file for easy plotting. + """ + + @staticmethod + def write(filename: str, x: numpy.ndarray, slip: numpy.ndarray): + """Write rupture information to SimpleGridDB ASCII file. + + Args: + filename: Name of spatial database file. + x: Along-fault coordinate (m) + slip: Slip (m) along fault + """ + data = numpy.stack((x, slip)).T + numpy.savetxt(filename, data, fmt=("%8.1f", "%5.3f")) + + +class Randomizer(abc.ABC): + """Abstract base class for generating random distribution. + """ + SHEAR_MODULUS = 3.3e+10 # Pa + SEED = 234 + + def __init__(self, eq_magnitude: float, total_length: float, dx: float): + """Constructor. + + Args: + eq_magnitude: Earthquake magnitude. + total_length: Total length of rupture. + dx: Discretization size for slip profile along fault. + """ + self.eq_magnitude = eq_magnitude + self.total_length = total_length + self.dx = dx + numpy.random.seed(self.SEED) + + @staticmethod + def _inverse_fft(amplitude: numpy.ndarray, phase: numpy.ndarray) -> numpy.ndarray: + """Compute inverse FFT to get slip profile from amplitude and phase. + + Args: + amplitude: Slip amplitude. + phase: Slip phase. + + Returns: + Slip profile. + """ + amplitudeFFT = amplitude * (numpy.cos(phase) + numpy.sin(phase)*complex(0,1)) + return numpy.real(numpy.fft.ifft(amplitudeFFT)) + + @staticmethod + def _avg_slip(eq_magnitude: float) -> float: + """Compute average slip given earthquake magnitude. + + Args: + eq_magnitude: Earthquake magnitude. + + Returns: + Average slip (m). + """ + area = 10**(0.8*(eq_magnitude - 3.30)) + area *= 1.0e+6 # km**2 -> m**2 + seismic_moment = 10**(1.5*(eq_magnitude+10.7)-7) + avg_slip = seismic_moment / (area * Randomizer.SHEAR_MODULUS) + return avg_slip + + @staticmethod + def _correlation_coef(eq_magnitude: float) -> float: + """Compute correlation coefficient. + + Args: + eq_magnitude: Earthquake magnitude. + + Returns: + Correlation coefficient (m). + """ + return 1*km * 10**(0.5*eq_magnitude - 2.5) + + +class RandomizerLA(Randomizer): + """Random slip profile based on Lavrentiadis and Abrahamson + (2019, doi: 10.1785/0120180252). + + We use a slightly larger value for phi_0. + """ + + Np = 1.24 + #PHI0 = 0.265 # Original value used by Lavrentiadis and Abrahamson + PHI0 = 0.35 # Modified value + + def create(self, x: numpy.ndarray) -> numpy.ndarray: + """Create 1D slip profile. + + Args: + x: Along-fault coordinate. + + Returns: + Slip profile. + """ + mag = self._create_amplitude(x) + phase = self._create_phase(x) + + slip = numpy.abs(self._inverse_fft(mag, phase)) + i_min = numpy.argmin(slip) + slip = numpy.concatenate((slip[i_min:], slip[:i_min])) + return slip + + def _create_amplitude(self, x: numpy.ndarray) -> numpy.ndarray: + """Generate slip amplitude. + + Args: + x: Along-fault coordinate. + + Returns: + Slip amplitude. + """ + avg_slip_m = self._avg_slip(self.eq_magnitude) + + kx = numpy.fft.fftfreq(len(x), self.dx / km) + k = numpy.abs(kx) + kc = self._kc() + scale0 = avg_slip_m * len(x) + scale = scale0 / numpy.sqrt(1.0+(k/kc)**(2*self.Np)) + + amplitude = scale * 10**(numpy.random.normal(loc=0.0, scale=self.PHI0, size=len(scale))) + amplitude[k==0] = scale0 + return amplitude + + def _create_phase(self, x: numpy.ndarray) -> numpy.ndarray: + """Generate slip phase. + + Args: + x: Along-fault coordinate. + + Returns: + Slip phase. + """ + total_length_km = self.total_length / km + dx_km = self.dx / km + + mu = -total_length_km * math.pi + phase_der = numpy.random.logistic(loc=mu, scale=self._s(), size=len(x)) + + dk = 1.0 / (dx_km * len(x)) + phase = scipy.integrate.cumulative_simpson(phase_der, dx=dk, initial=0.0) + phase = phase % (2.0*math.pi) + return phase + + def _kc(self) -> float: + """Compute kc parameter. + + Returns: + kc parameter. + """ + c1 = -2.03 + c2 = -1.0 + c3 = 1.6 + total_length_km = self.total_length / km + return 10**(c1 + c2*(math.log10(total_length_km) - c3)) + + def _s(self) -> float: + """Compute asa parameter. + + Returns: + `s` parameter. + """ + d1 = 1.49 + d2 = 1.0 + d3 = 1.6 + total_length_km = self.total_length / km + return 10**(d1 + d2*(math.log10(total_length_km) - d3)) + + +class Taper: + """Object for applying linear taper at ends of slip profile. + """ + + def __init__(self, total_length: float, taper_length: float): + """Constructor. + + Args: + total_length: Total rupture length. + taper_length: Length of linear taper at each end. + """ + self.total_length = total_length + self.taper_length = taper_length + + def apply(self, x: numpy.ndarray, slip: numpy.ndarray): + """Apply linear taper at ends of slip profile. + + Args: + x: Along-fault coordinate. + slip: Slip profile to apply taper to. + """ + l_taper = self.taper_length.value + maskL = x <= x[0] + l_taper + maskR = x >= x[-1] - l_taper + taper = numpy.logical_and(~maskL, ~maskR) + maskL * (x-x[0]) / l_taper + maskR * (x[-1] - x) / l_taper + slip *= taper + + +if __name__ == "__main__": + cli() diff --git a/examples/strikeslip-2d/invert_slip.py b/examples/strikeslip-2d/invert_slip.py index b1a3e2a9be..f5b3aa3d4a 100755 --- a/examples/strikeslip-2d/invert_slip.py +++ b/examples/strikeslip-2d/invert_slip.py @@ -11,40 +11,64 @@ consistent ordering. """ -# Import argparse Python module (standard Python) -import argparse +import pathlib # Import numpy and h5py modules (included in PyLith installation) import numpy import h5py +OUTPUT_DIR = pathlib.Path("output") + +def cli(): + """Command line interface. + """ + import argparse + parser = argparse.ArgumentParser() + + parser.add_argument("--impulse-data", action="store", type=str, dest="filename_fault", default="step05_greensfns-fault.h5", help="Fault HDF5 data file from Green's function simulation.") + parser.add_argument("--impulse-responses", action="store", type=str, dest="filename_responses", default="step05_greensfns-gnss_stations.h5", help="Station HDF5 data file from Green's function simulation.") + parser.add_argument("--observed-data", action="store", type=str, dest="filename_observed", default="step04_varslip-gnss_stations.h5", help="Station HDF5 data file from variable slip simulation.") + parser.add_argument("--penalties", action="store", type=str, dest="penalties", default="0.02,0.1,0.4", help="Comma separated list of penalties.") + parser.add_argument("--output", action="store", type=str, dest="filename_output", default="step06_inversion-results.txt", help="Name of output file with inversion results.") + + args = parser.parse_args() + app = InvertSlipApp(output_dir=OUTPUT_DIR) + app.run( + filename_fault=args.filename_fault, + filename_responses=args.filename_responses, + filename_observed=args.filename_observed, + filename_output=args.filename_output, + penalties=[float(p) for p in args.penalties.split(",")], + ) + + class InvertSlipApp: """Application to invert for fault slip using PyLith-generated Green's functions. """ - def __init__(self): - self.filename_fault = "output/step05_greensfns-fault.h5" - self.filename_responses = "output/step05_greensfns-gnss_stations.h5" - self.filename_observed = "output/step04_varslip-gnss_stations.h5" - self.filename_output = "output/step06_inversion-results.txt" - self.penalties = [0.01, 0.1, 1.0] - - def main(self): - """Entry point for running application. + def __init__(self, output_dir: str): + """Constructor + + Args: + output_dir: Directory with simulation output. """ - args = self._parse_command_line() + self.output_dir = output_dir - self.get_fault_impulses(args.filename_fault) - self.get_station_responses(args.filename_responses) - self.get_station_observed(args.filename_observed) + def run(self, filename_fault: str, filename_responses: str, filename_observed: str, filename_output: str, penalties: list): + """Entry point for running application. + """ + self._get_fault_impulses(self.output_dir / filename_fault) + self._get_station_responses(self.output_dir / filename_responses) + self._get_station_observed(self.output_dir / filename_observed) - self.penalties = list(map(float, args.penalties.split(","))) - results = self.invert() - self.write_results(self.filename_output, results) + self.penalties = penalties + results = self._invert() + self.write_results(self.output_dir / filename_output, results) - def get_fault_impulses(self, filename): + def _get_fault_impulses(self, filename: str): """Get coordinates, amplitude, and order of impulses. Fault points are sorted by y coordinate. - :param filename: Name of HDF5 file with fault data from PyLith Green's function simulation (Step 5). + Args: + filename: Name of HDF5 file with fault data from PyLith Green's function simulation (Step 5). """ h5 = h5py.File(filename, "r") y = h5['geometry/vertices'][:,1] @@ -57,11 +81,12 @@ def get_fault_impulses(self, filename): self.impulse_y = y[reorder] self.impulse_slip = slip[:,reorder] - def get_station_responses(self, filename): + def _get_station_responses(self, filename): """ Get coordinates and displacements at stations for Green's function responses. - :param filename: Name of HDF5 file with point data from PyLith Green's function simulation (Step 5_). + Args: + filename: Name of HDF5 file with point data from PyLith Green's function simulation (Step 5_). """ h5 = h5py.File(filename, "r") xy = h5['geometry/vertices'][:,0:2] @@ -72,11 +97,12 @@ def get_station_responses(self, filename): reorder = numpy.argsort(xy[:,1]) self.station_responses = displacement[:,reorder,:] - def get_station_observed(self, filename): + def _get_station_observed(self, filename): """ Get coordinates and displacements at stations for fake observations. - :param filename: Name of HDF5 file with point data from PyLith forward simulation (Step 4). + Args: + filename: Name of HDF5 file with point data from PyLith forward simulation (Step 4). """ h5 = h5py.File(filename, "r") xy = h5['geometry/vertices'][:,0:2] @@ -87,7 +113,7 @@ def get_station_observed(self, filename): reorder = numpy.argsort(xy[:,1]) self.station_observed = displacement[:,reorder,:] - def invert(self): + def _invert(self) -> numpy.ndarray: """Invert observations for amplitude of impulses and fault slip. """ # Determine matrix sizes and set up A-matrix. @@ -129,11 +155,12 @@ def invert(self): print(f"Residual norm: {residual_norm}") return results - def write_results(self, filename, results): + def write_results(self, filename: str, results: numpy.ndarray): """Write inversion results to file. - :param filename: Name of output file. - :param results: Inversion results as Numpy array. + Args: + filename: Name of output file. + results: Inversion results as Numpy array. """ header = "# y" for penalty in self.penalties: @@ -144,22 +171,9 @@ def write_results(self, filename, results): fout.write(header) numpy.savetxt(fout, results, fmt="%14.6e") - def _parse_command_line(self): - """Parse command line arguments. - - :returns: Command line arugment information as argparse.Namespace. - """ - parser = argparse.ArgumentParser() - - parser.add_argument("--impulse-data", action="store", type=str, dest="filename_fault", default=self.filename_fault, help="Fault HDF5 data file from Green's function simulation.") - parser.add_argument("--impulse-responses", action="store", type=str, dest="filename_responses", default=self.filename_responses, help="Station HDF5 data file from Green's function simulation.") - parser.add_argument("--observed-data", action="store", type=str, dest="filename_observed", default=self.filename_observed, help="Station HDF5 data file from variable slip simulation.") - parser.add_argument("--penalties", action="store", type=str, dest="penalties", default="0.01,0.1,1.0", help="Comma separated list of penalties.") - parser.add_argument("--output", action="store", type=str, dest="filename_output", default=self.filename_output, help="Name of output file with inversion results.") - return parser.parse_args() if __name__ == "__main__": - InvertSlipApp().main() + cli() # End of file diff --git a/examples/strikeslip-2d/slip_variable.spatialdb b/examples/strikeslip-2d/slip_variable.spatialdb index 1180c4d334..ef864f28ae 100644 --- a/examples/strikeslip-2d/slip_variable.spatialdb +++ b/examples/strikeslip-2d/slip_variable.spatialdb @@ -1,33 +1,424 @@ -// This spatial database specifies the distribution of slip on the -// fault surface for Step 4 (forward simulation to generate fake observations). -// The slip extends from -20 <= y <= 20 km. -#SPATIAL.ascii 1 -SimpleDB { - num-values = 3 - value-names = final_slip_left_lateral final_slip_opening initiation_time - value-units = cm cm year - num-locs = 11 - data-dim = 1 // Data is specified along a line. +#SPATIAL_GRID.ascii 1 +SimpleGridDB { + num-x = 1 + num-y = 403 + num-z = 0 space-dim = 2 + num-values = 3 + value-names = final_slip_left_lateral final_slip_opening initiation_time + value-units = m m s cs-data = cartesian { - to-meters = 1.0e+3 // Specify coordinates in km for convenience. - space-dim = 2 - } // cs-data -} // SimpleDB -// Columns are -// (1) x coordinate (km) -// (2) y coordinate (km) -// (3) left-lateral slip (cm) -// (4) fault opening (cm) -// (5) slip time relative to origin time (year) -0.0 99.0 0.0 0.0 0.0 -0.0 20.0 0.0 0.0 0.0 -0.0 15.0 30.0 0.0 0.0 -0.0 10.0 40.0 0.0 0.0 -0.0 5.0 20.0 0.0 0.0 -0.0 0.0 60.0 0.0 0.0 -0.0 -5.0 80.0 0.0 0.0 -0.0 -10.0 30.0 0.0 0.0 -0.0 -15.0 10.0 0.0 0.0 -0.0 -20.0 0.0 0.0 0.0 -0.0 -99.0 0.0 0.0 0.0 + to-meters = 1 + space-dim = 2 +} +} +// x-coordinates + 0.000000e+00 +// y-coordinates + -7.500000e+04 -2.000000e+04 -1.990000e+04 -1.980000e+04 -1.970000e+04 -1.960000e+04 -1.950000e+04 -1.940000e+04 -1.930000e+04 -1.920000e+04 -1.910000e+04 -1.900000e+04 -1.890000e+04 -1.880000e+04 -1.870000e+04 -1.860000e+04 -1.850000e+04 -1.840000e+04 -1.830000e+04 -1.820000e+04 -1.810000e+04 -1.800000e+04 -1.790000e+04 -1.780000e+04 -1.770000e+04 -1.760000e+04 -1.750000e+04 -1.740000e+04 -1.730000e+04 -1.720000e+04 -1.710000e+04 -1.700000e+04 -1.690000e+04 -1.680000e+04 -1.670000e+04 -1.660000e+04 -1.650000e+04 -1.640000e+04 -1.630000e+04 -1.620000e+04 -1.610000e+04 -1.600000e+04 -1.590000e+04 -1.580000e+04 -1.570000e+04 -1.560000e+04 -1.550000e+04 -1.540000e+04 -1.530000e+04 -1.520000e+04 -1.510000e+04 -1.500000e+04 -1.490000e+04 -1.480000e+04 -1.470000e+04 -1.460000e+04 -1.450000e+04 -1.440000e+04 -1.430000e+04 -1.420000e+04 -1.410000e+04 -1.400000e+04 -1.390000e+04 -1.380000e+04 -1.370000e+04 -1.360000e+04 -1.350000e+04 -1.340000e+04 -1.330000e+04 -1.320000e+04 -1.310000e+04 -1.300000e+04 -1.290000e+04 -1.280000e+04 -1.270000e+04 -1.260000e+04 -1.250000e+04 -1.240000e+04 -1.230000e+04 -1.220000e+04 -1.210000e+04 -1.200000e+04 -1.190000e+04 -1.180000e+04 -1.170000e+04 -1.160000e+04 -1.150000e+04 -1.140000e+04 -1.130000e+04 -1.120000e+04 -1.110000e+04 -1.100000e+04 -1.090000e+04 -1.080000e+04 -1.070000e+04 -1.060000e+04 -1.050000e+04 -1.040000e+04 -1.030000e+04 -1.020000e+04 -1.010000e+04 -1.000000e+04 -9.900000e+03 -9.800000e+03 -9.700000e+03 -9.600000e+03 -9.500000e+03 -9.400000e+03 -9.300000e+03 -9.200000e+03 -9.100000e+03 -9.000000e+03 -8.900000e+03 -8.800000e+03 -8.700000e+03 -8.600000e+03 -8.500000e+03 -8.400000e+03 -8.300000e+03 -8.200000e+03 -8.100000e+03 -8.000000e+03 -7.900000e+03 -7.800000e+03 -7.700000e+03 -7.600000e+03 -7.500000e+03 -7.400000e+03 -7.300000e+03 -7.200000e+03 -7.100000e+03 -7.000000e+03 -6.900000e+03 -6.800000e+03 -6.700000e+03 -6.600000e+03 -6.500000e+03 -6.400000e+03 -6.300000e+03 -6.200000e+03 -6.100000e+03 -6.000000e+03 -5.900000e+03 -5.800000e+03 -5.700000e+03 -5.600000e+03 -5.500000e+03 -5.400000e+03 -5.300000e+03 -5.200000e+03 -5.100000e+03 -5.000000e+03 -4.900000e+03 -4.800000e+03 -4.700000e+03 -4.600000e+03 -4.500000e+03 -4.400000e+03 -4.300000e+03 -4.200000e+03 -4.100000e+03 -4.000000e+03 -3.900000e+03 -3.800000e+03 -3.700000e+03 -3.600000e+03 -3.500000e+03 -3.400000e+03 -3.300000e+03 -3.200000e+03 -3.100000e+03 -3.000000e+03 -2.900000e+03 -2.800000e+03 -2.700000e+03 -2.600000e+03 -2.500000e+03 -2.400000e+03 -2.300000e+03 -2.200000e+03 -2.100000e+03 -2.000000e+03 -1.900000e+03 -1.800000e+03 -1.700000e+03 -1.600000e+03 -1.500000e+03 -1.400000e+03 -1.300000e+03 -1.200000e+03 -1.100000e+03 -1.000000e+03 -9.000000e+02 -8.000000e+02 -7.000000e+02 -6.000000e+02 -5.000000e+02 -4.000000e+02 -3.000000e+02 -2.000000e+02 -1.000000e+02 0.000000e+00 1.000000e+02 2.000000e+02 3.000000e+02 4.000000e+02 5.000000e+02 6.000000e+02 7.000000e+02 8.000000e+02 9.000000e+02 1.000000e+03 1.100000e+03 1.200000e+03 1.300000e+03 1.400000e+03 1.500000e+03 1.600000e+03 1.700000e+03 1.800000e+03 1.900000e+03 2.000000e+03 2.100000e+03 2.200000e+03 2.300000e+03 2.400000e+03 2.500000e+03 2.600000e+03 2.700000e+03 2.800000e+03 2.900000e+03 3.000000e+03 3.100000e+03 3.200000e+03 3.300000e+03 3.400000e+03 3.500000e+03 3.600000e+03 3.700000e+03 3.800000e+03 3.900000e+03 4.000000e+03 4.100000e+03 4.200000e+03 4.300000e+03 4.400000e+03 4.500000e+03 4.600000e+03 4.700000e+03 4.800000e+03 4.900000e+03 5.000000e+03 5.100000e+03 5.200000e+03 5.300000e+03 5.400000e+03 5.500000e+03 5.600000e+03 5.700000e+03 5.800000e+03 5.900000e+03 6.000000e+03 6.100000e+03 6.200000e+03 6.300000e+03 6.400000e+03 6.500000e+03 6.600000e+03 6.700000e+03 6.800000e+03 6.900000e+03 7.000000e+03 7.100000e+03 7.200000e+03 7.300000e+03 7.400000e+03 7.500000e+03 7.600000e+03 7.700000e+03 7.800000e+03 7.900000e+03 8.000000e+03 8.100000e+03 8.200000e+03 8.300000e+03 8.400000e+03 8.500000e+03 8.600000e+03 8.700000e+03 8.800000e+03 8.900000e+03 9.000000e+03 9.100000e+03 9.200000e+03 9.300000e+03 9.400000e+03 9.500000e+03 9.600000e+03 9.700000e+03 9.800000e+03 9.900000e+03 1.000000e+04 1.010000e+04 1.020000e+04 1.030000e+04 1.040000e+04 1.050000e+04 1.060000e+04 1.070000e+04 1.080000e+04 1.090000e+04 1.100000e+04 1.110000e+04 1.120000e+04 1.130000e+04 1.140000e+04 1.150000e+04 1.160000e+04 1.170000e+04 1.180000e+04 1.190000e+04 1.200000e+04 1.210000e+04 1.220000e+04 1.230000e+04 1.240000e+04 1.250000e+04 1.260000e+04 1.270000e+04 1.280000e+04 1.290000e+04 1.300000e+04 1.310000e+04 1.320000e+04 1.330000e+04 1.340000e+04 1.350000e+04 1.360000e+04 1.370000e+04 1.380000e+04 1.390000e+04 1.400000e+04 1.410000e+04 1.420000e+04 1.430000e+04 1.440000e+04 1.450000e+04 1.460000e+04 1.470000e+04 1.480000e+04 1.490000e+04 1.500000e+04 1.510000e+04 1.520000e+04 1.530000e+04 1.540000e+04 1.550000e+04 1.560000e+04 1.570000e+04 1.580000e+04 1.590000e+04 1.600000e+04 1.610000e+04 1.620000e+04 1.630000e+04 1.640000e+04 1.650000e+04 1.660000e+04 1.670000e+04 1.680000e+04 1.690000e+04 1.700000e+04 1.710000e+04 1.720000e+04 1.730000e+04 1.740000e+04 1.750000e+04 1.760000e+04 1.770000e+04 1.780000e+04 1.790000e+04 1.800000e+04 1.810000e+04 1.820000e+04 1.830000e+04 1.840000e+04 1.850000e+04 1.860000e+04 1.870000e+04 1.880000e+04 1.890000e+04 1.900000e+04 1.910000e+04 1.920000e+04 1.930000e+04 1.940000e+04 1.950000e+04 1.960000e+04 1.970000e+04 1.980000e+04 1.990000e+04 2.000000e+04 7.500000e+04 +// z-coordinates + +// data + 0.000000e+00 -7.500000e+04 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.000000e+04 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.990000e+04 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.980000e+04 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.970000e+04 2.412684e-03 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.960000e+04 2.256749e-02 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.950000e+04 4.867393e-02 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.940000e+04 7.849781e-02 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.930000e+04 1.098049e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.920000e+04 1.403611e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.910000e+04 1.682456e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.900000e+04 1.927916e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.890000e+04 2.136458e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.880000e+04 2.304548e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.870000e+04 2.433468e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.860000e+04 2.543759e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.850000e+04 2.660779e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.840000e+04 2.809884e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.830000e+04 3.005525e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.820000e+04 3.218530e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.810000e+04 3.408820e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.800000e+04 3.536317e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.790000e+04 3.577180e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.780000e+04 3.572513e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.770000e+04 3.579658e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.760000e+04 3.655955e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.750000e+04 3.843875e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.740000e+04 4.126400e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.730000e+04 4.471640e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.720000e+04 4.847705e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.710000e+04 5.221622e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.700000e+04 5.556079e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.690000e+04 5.812681e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.680000e+04 5.953032e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.670000e+04 5.960158e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.660000e+04 5.902772e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.650000e+04 5.871010e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.640000e+04 5.955007e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.630000e+04 6.215960e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.620000e+04 6.599316e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.610000e+04 7.021583e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.600000e+04 7.399272e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.590000e+04 7.665415e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.580000e+04 7.819146e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.570000e+04 7.876126e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.560000e+04 7.852013e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.550000e+04 7.765686e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.540000e+04 7.648898e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.530000e+04 7.536622e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.520000e+04 7.463831e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.510000e+04 7.459639e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.500000e+04 7.529715e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.490000e+04 7.673874e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.480000e+04 7.891925e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.470000e+04 8.171729e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.460000e+04 8.453338e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.450000e+04 8.664853e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.440000e+04 8.734374e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.430000e+04 8.615923e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.420000e+04 8.367204e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.410000e+04 8.071842e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.400000e+04 7.813462e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.390000e+04 7.655544e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.380000e+04 7.580990e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.370000e+04 7.552559e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.360000e+04 7.533008e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.350000e+04 7.495864e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.340000e+04 7.457725e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.330000e+04 7.445958e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.320000e+04 7.487930e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.310000e+04 7.599690e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.300000e+04 7.752009e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.290000e+04 7.904341e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.280000e+04 8.016140e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.270000e+04 8.054088e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.260000e+04 8.013786e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.250000e+04 7.898066e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.240000e+04 7.709756e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.230000e+04 7.458569e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.220000e+04 7.181742e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.210000e+04 6.923394e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.200000e+04 6.727643e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.190000e+04 6.625286e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.180000e+04 6.593831e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.170000e+04 6.597463e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.160000e+04 6.600367e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.150000e+04 6.575607e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.140000e+04 6.531759e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.130000e+04 6.486277e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.120000e+04 6.456617e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.110000e+04 6.451082e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.100000e+04 6.441376e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.090000e+04 6.390052e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.080000e+04 6.259663e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.070000e+04 6.027037e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.060000e+04 5.726102e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.050000e+04 5.405060e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.040000e+04 5.112111e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.030000e+04 4.881081e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.020000e+04 4.688279e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.010000e+04 4.495638e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.000000e+04 4.265090e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.900000e+03 3.973335e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.800000e+03 3.656145e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.700000e+03 3.364061e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.600000e+03 3.147624e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.500000e+03 3.041734e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.400000e+03 3.018743e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.300000e+03 3.035361e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.200000e+03 3.048298e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.100000e+03 3.026867e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.000000e+03 2.990775e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.900000e+03 2.972333e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.800000e+03 3.003848e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.700000e+03 3.105374e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.600000e+03 3.247942e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.500000e+03 3.390328e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.400000e+03 3.491309e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.300000e+03 3.520572e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.200000e+03 3.491455e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.100000e+03 3.428207e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.000000e+03 3.355079e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.900000e+03 3.293136e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.800000e+03 3.250702e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.700000e+03 3.232917e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.600000e+03 3.244921e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.500000e+03 3.288110e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.400000e+03 3.348907e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.300000e+03 3.409993e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.200000e+03 3.454046e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.100000e+03 3.468976e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.000000e+03 3.463594e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.900000e+03 3.451944e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.800000e+03 3.448065e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.700000e+03 3.462225e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.600000e+03 3.489599e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.500000e+03 3.521587e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.400000e+03 3.549590e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.300000e+03 3.567982e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.200000e+03 3.583036e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.100000e+03 3.603998e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.000000e+03 3.640112e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.900000e+03 3.696177e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.800000e+03 3.759192e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.700000e+03 3.811710e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.600000e+03 3.836282e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.500000e+03 3.822439e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.400000e+03 3.787627e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.300000e+03 3.756270e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.200000e+03 3.752794e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.100000e+03 3.791992e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.000000e+03 3.850131e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.900000e+03 3.893848e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.800000e+03 3.889777e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.700000e+03 3.817088e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.600000e+03 3.705073e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.500000e+03 3.595561e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.400000e+03 3.530375e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.300000e+03 3.541422e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.200000e+03 3.620918e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.100000e+03 3.751156e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.000000e+03 3.914430e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.900000e+03 4.092948e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.800000e+03 4.268563e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.700000e+03 4.423039e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.600000e+03 4.538143e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.500000e+03 4.602813e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.400000e+03 4.634682e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.300000e+03 4.658553e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.200000e+03 4.699234e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.100000e+03 4.775196e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.000000e+03 4.879574e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.900000e+03 4.999171e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.800000e+03 5.120790e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.700000e+03 5.233784e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.600000e+03 5.337717e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.500000e+03 5.434702e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.400000e+03 5.526852e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.300000e+03 5.616372e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.200000e+03 5.705822e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.100000e+03 5.797851e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.000000e+03 5.895108e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.900000e+03 5.998409e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.800000e+03 6.101233e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.700000e+03 6.195224e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.600000e+03 6.272026e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.500000e+03 6.325783e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.400000e+03 6.360630e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.300000e+03 6.383200e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.200000e+03 6.400127e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.100000e+03 6.416080e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.000000e+03 6.427876e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.000000e+02 6.430370e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.000000e+02 6.418415e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.000000e+02 6.387166e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.000000e+02 6.332977e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.000000e+02 6.252501e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.000000e+02 6.142392e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.000000e+02 6.003878e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.000000e+02 5.856475e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.000000e+02 5.724272e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 0.000000e+00 5.631357e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.000000e+02 5.595985e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.000000e+02 5.613062e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.000000e+02 5.671661e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.000000e+02 5.760856e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.000000e+02 5.869534e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.000000e+02 5.985844e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.000000e+02 6.097753e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.000000e+02 6.193226e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.000000e+02 6.264477e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.000000e+03 6.320722e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.100000e+03 6.375423e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.200000e+03 6.442043e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.300000e+03 6.529204e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.400000e+03 6.626162e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.500000e+03 6.717328e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.600000e+03 6.787116e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.700000e+03 6.825109e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.800000e+03 6.841578e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.900000e+03 6.851964e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.000000e+03 6.871706e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.100000e+03 6.913161e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.200000e+03 6.976341e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.300000e+03 7.058170e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.400000e+03 7.155576e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.500000e+03 7.261780e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.600000e+03 7.355185e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.700000e+03 7.410489e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.800000e+03 7.402391e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 2.900000e+03 7.316842e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.000000e+03 7.184801e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.100000e+03 7.048480e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.200000e+03 6.950088e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.300000e+03 6.919010e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.400000e+03 6.933313e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.500000e+03 6.958235e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.600000e+03 6.959015e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.700000e+03 6.910968e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.800000e+03 6.829714e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 3.900000e+03 6.740951e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.000000e+03 6.670375e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.100000e+03 6.635449e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.200000e+03 6.620703e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.300000e+03 6.602436e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.400000e+03 6.556943e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.500000e+03 6.468512e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.600000e+03 6.353387e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.700000e+03 6.235803e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.800000e+03 6.139996e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 4.900000e+03 6.082016e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.000000e+03 6.045185e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.100000e+03 6.004639e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.200000e+03 5.935514e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.300000e+03 5.823580e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.400000e+03 5.697130e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.500000e+03 5.595090e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.600000e+03 5.556385e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.700000e+03 5.606441e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.800000e+03 5.716692e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 5.900000e+03 5.845071e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.000000e+03 5.949512e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.100000e+03 5.998187e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.200000e+03 6.000212e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.300000e+03 5.974942e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.400000e+03 5.941732e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.500000e+03 5.917137e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.600000e+03 5.906526e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.700000e+03 5.912469e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.800000e+03 5.937534e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 6.900000e+03 5.983378e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.000000e+03 6.047993e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.100000e+03 6.128457e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.200000e+03 6.221850e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.300000e+03 6.325833e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.400000e+03 6.440407e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.500000e+03 6.566156e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.600000e+03 6.703664e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.700000e+03 6.851582e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.800000e+03 7.000816e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 7.900000e+03 7.140342e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.000000e+03 7.259133e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.100000e+03 7.349952e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.200000e+03 7.420721e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.300000e+03 7.483151e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.400000e+03 7.548955e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.500000e+03 7.622325e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.600000e+03 7.677374e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.700000e+03 7.680695e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.800000e+03 7.598884e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 8.900000e+03 7.411567e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.000000e+03 7.150500e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.100000e+03 6.860474e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.200000e+03 6.586278e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.300000e+03 6.364743e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.400000e+03 6.200865e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.500000e+03 6.091678e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.600000e+03 6.034221e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.700000e+03 6.025935e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.800000e+03 6.065898e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 9.900000e+03 6.153590e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.000000e+04 6.288495e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.010000e+04 6.464937e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.020000e+04 6.656608e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.030000e+04 6.832041e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.040000e+04 6.959767e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.050000e+04 7.020821e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.060000e+04 7.046232e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.070000e+04 7.079531e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.080000e+04 7.164250e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.090000e+04 7.331307e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.100000e+04 7.561178e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.110000e+04 7.821728e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.120000e+04 8.080822e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.130000e+04 8.310616e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.140000e+04 8.500442e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.150000e+04 8.643920e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.160000e+04 8.734674e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.170000e+04 8.770231e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.180000e+04 8.763741e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.190000e+04 8.732262e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.200000e+04 8.692847e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.210000e+04 8.657004e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.220000e+04 8.614041e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.230000e+04 8.547717e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.240000e+04 8.441791e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.250000e+04 8.283698e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.260000e+04 8.075575e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.270000e+04 7.823238e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.280000e+04 7.532500e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.290000e+04 7.212912e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.300000e+04 6.888964e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.310000e+04 6.588886e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.320000e+04 6.340903e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.330000e+04 6.164106e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.340000e+04 6.041035e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.350000e+04 5.945089e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.360000e+04 5.849670e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.370000e+04 5.732284e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.380000e+04 5.586845e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.390000e+04 5.411375e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.400000e+04 5.203892e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.410000e+04 4.966009e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.420000e+04 4.713706e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.430000e+04 4.466558e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.440000e+04 4.244137e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.450000e+04 4.063868e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.460000e+04 3.934576e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.470000e+04 3.862939e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.480000e+04 3.855634e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.490000e+04 3.909772e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.500000e+04 3.984205e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.510000e+04 4.028217e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.520000e+04 3.991094e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.530000e+04 3.838987e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.540000e+04 3.605505e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.550000e+04 3.341123e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.560000e+04 3.096315e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.570000e+04 2.909893e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.580000e+04 2.774015e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.590000e+04 2.669174e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.600000e+04 2.575866e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.610000e+04 2.481058e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.620000e+04 2.397616e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.630000e+04 2.344880e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.640000e+04 2.342190e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.650000e+04 2.400258e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.660000e+04 2.495287e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.670000e+04 2.594855e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.680000e+04 2.666536e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.690000e+04 2.688999e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.700000e+04 2.685271e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.710000e+04 2.689472e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.720000e+04 2.735718e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.730000e+04 2.848830e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.740000e+04 3.016425e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.750000e+04 3.216823e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.760000e+04 3.428342e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.770000e+04 3.630299e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.780000e+04 3.806006e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.790000e+04 3.939771e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.800000e+04 4.015906e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.810000e+04 4.021254e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.820000e+04 3.952800e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.830000e+04 3.810066e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.840000e+04 3.592572e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.850000e+04 3.305635e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.860000e+04 2.977753e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.870000e+04 2.643223e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.880000e+04 2.336340e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.890000e+04 2.084238e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.900000e+04 1.885409e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.910000e+04 1.731183e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.920000e+04 1.612889e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.930000e+04 1.520061e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.940000e+04 1.435048e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.950000e+04 1.338403e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.960000e+04 1.210676e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.970000e+04 1.032420e-01 0.000000e+00 0.000000e+00 + 0.000000e+00 1.980000e+04 7.841873e-02 0.000000e+00 0.000000e+00 + 0.000000e+00 1.990000e+04 4.465299e-02 0.000000e+00 0.000000e+00 + 0.000000e+00 2.000000e+04 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 7.500000e+04 0.000000e+00 0.000000e+00 0.000000e+00 diff --git a/examples/strikeslip-2d/step04_varslip.cfg b/examples/strikeslip-2d/step04_varslip.cfg index 81e7cbc236..78762a868f 100644 --- a/examples/strikeslip-2d/step04_varslip.cfg +++ b/examples/strikeslip-2d/step04_varslip.cfg @@ -40,6 +40,12 @@ features = [ # output filenames. The default directory for output is 'output'. defaults.name = step04_varslip +# ---------------------------------------------------------------------- +# mesh_generator +# ---------------------------------------------------------------------- +[pylithapp.mesh_generator] +refiner = pylith.topology.RefineUniform + # ---------------------------------------------------------------------- # problem # ---------------------------------------------------------------------- @@ -94,11 +100,15 @@ derived_subfields.cauchy_stress.basis_order = 1 # ---------------------------------------------------------------------- # fault # ---------------------------------------------------------------------- +# Refine output by a factor of 2**3=8 to capture higher order discretization. +[pylithapp.problem.interfaces.fault] +observers.observer.refine_levels = 3 + # Specify slip on the fault using a `SimpleDB` spatial database. [pylithapp.problem.interfaces.fault.eq_ruptures.rupture] -db_auxiliary_field = spatialdata.spatialdb.SimpleDB +db_auxiliary_field = spatialdata.spatialdb.SimpleGridDB db_auxiliary_field.description = Fault rupture auxiliary field spatial database -db_auxiliary_field.iohandler.filename = slip_variable.spatialdb +db_auxiliary_field.filename = slip_variable.spatialdb db_auxiliary_field.query_type = linear diff --git a/examples/strikeslip-2d/step05_greensfns.cfg b/examples/strikeslip-2d/step05_greensfns.cfg index 44e495cd28..1d607436c3 100644 --- a/examples/strikeslip-2d/step05_greensfns.cfg +++ b/examples/strikeslip-2d/step05_greensfns.cfg @@ -52,6 +52,7 @@ greensfns = 1 # ---------------------------------------------------------------------- [pylithapp.mesh_generator] refiner = pylith.topology.RefineUniform +refiner.levels = 1 # ---------------------------------------------------------------------- # problem diff --git a/examples/strikeslip-2d/viz/plot_catmip_results.py b/examples/strikeslip-2d/viz/plot_catmip_results.py index 684bd65040..67c275ed77 100755 --- a/examples/strikeslip-2d/viz/plot_catmip_results.py +++ b/examples/strikeslip-2d/viz/plot_catmip_results.py @@ -4,38 +4,64 @@ the true solution for Step 7 (CATMIP inversions). """ -# Import argparse and re (regular expressions) Python modules(standard Python) -import argparse +# Standard Python modules import pathlib -# Import numpy, h5py, and matplotlib modules (included in PyLith binary installation) +# Python modules included in PyLith binary installation import numpy import h5py import matplotlib.pyplot as pyplot +from pythia.pyre.units.length import km -class PlotApp: +OUTPUT_DIR = pathlib.Path("output") - FILENAME_IMPULSES = "output/step05_greensfns-fault.h5" - FILENAME_OBSERVED = "output/step04_varslip-fault.h5" - - COLOR_OBSERVED = "black" - COLOR_MODEL = "orange" - XLIM = (-25.0, 25.0) +FILENAME_IMPULSES = OUTPUT_DIR / "step05_greensfns-fault.h5" +FILENAME_PRESCRIBED = OUTPUT_DIR / "step04_varslip-fault.h5" +FILENAME_RAW = OUTPUT_DIR / "slip_variable.txt" - def __init__(self): - self.observed_coords = None - self.observed_slip = None - self.impulse_coords = None - self.impulse_slip = None +def cli(): + """Command line interface. + """ + import argparse - self.inversion_coefs = None - - def main(self, filename_theta: str, filename_observed: str=None, filename_impulses: str=None, show_plot: bool=False): - filename_observed = filename_observed or self.FILENAME_OBSERVED - filename_impulses = filename_impulses or self.FILENAME_IMPULSES - - self._load_observed(filename_observed) + parser = argparse.ArgumentParser() + parser.add_argument("--catmip-theta", action="store", dest="filename_theta", type=str, required=True, help="Filename of output from CATMIP inversion.") + parser.add_argument("--pylith-impulses", action="store", dest="filename_impulses", type=str, default=FILENAME_IMPULSES, help="Name of HDF5 file with fault slip inpulses from Green's functions simulation.") + parser.add_argument("--prescribed-slip", action="store", dest="filename_prescribed", type=str, default=FILENAME_PRESCRIBED, help="Name of HDF5 file for fault from prescribed slip simulation.") + parser.add_argument("--raw", action="store", dest="filename_raw", type=str, default=FILENAME_RAW, help="Name of ASCII file with raw slip profile.") + parser.add_argument("--no-gui", action="store_false", dest="show_plot", default=True, help="Do not display plot.") + + args = parser.parse_args() + PlotApp().run( + filename_prescribed=args.filename_prescribed, + filename_impulses=args.filename_impulses, + filename_theta=args.filename_theta, + filename_raw=args.filename_raw, + show_plot=args.show_plot, + ) + + +class PlotApp: + """Application to plot CATMIP inversion results. + """ + COLOR_RAW = "gray" + COLOR_OBSERVED = "black" + COLOR_MODEL = "orange" + XLIM = (-25.0, 35.0) + + def run(self, filename_theta: str, filename_prescribed: str, filename_impulses: str, filename_raw: str, show_plot: bool): + """Run plotting application. + + Args: + filename_theta: Filename of output from inversion. + filename_prescribed: Name of HDF5 file with station output with fake observations. + filename_raw: Name of ASCII file with raw slip values. + prescribed_only: Show only prescribed slip. + show_plot: Show plot window. + """ + self._load_observed(filename_prescribed) + self._load_raw(filename_raw) self._load_impulses(filename_impulses) self._load_inversion_results(filename_theta) @@ -45,9 +71,9 @@ def main(self, filename_theta: str, filename_observed: str=None, filename_impuls slip_median = self._calc_median() slip_stddev = self._calc_stddev() slip_minmax = self._calc_minmax() - self.plot(slip_median, slip_stddev, slip_minmax, filename_plot, show_plot) + self._plot(slip_median, slip_stddev, slip_minmax, filename_plot, show_plot) - def _load_observed(self, filename): + def _load_observed(self, filename: str): h5 = h5py.File(filename, "r") observed_coords = h5['geometry/vertices'][:] observed_slip = h5['vertex_fields/slip'][:,:,:].squeeze() @@ -58,7 +84,12 @@ def _load_observed(self, filename): self.observed_coords = observed_coords[reorder,:] self.observed_slip = observed_slip[reorder,:] - def _load_impulses(self, filename): + def _load_raw(self, filename: str): + data = numpy.loadtxt(filename) + self.raw_y = data[:,0] + self.raw_slip = data[:,1] + + def _load_impulses(self, filename: str): h5 = h5py.File(filename, "r") impulse_coords = h5['geometry/vertices'][:] impulse_slip = h5['vertex_fields/slip'][:,:,:].squeeze() @@ -69,41 +100,37 @@ def _load_impulses(self, filename): self.impulse_coords = impulse_coords[reorder,:] self.impulse_slip = impulse_slip[:,reorder,:] - def _load_inversion_results(self, filename): + def _load_inversion_results(self, filename: str): with open(filename, "rb") as fin: coefs = numpy.frombuffer(fin.read(), dtype=numpy.float64) nimpulses = self.impulse_slip.shape[0] self.inversion_coefs = coefs.reshape((nimpulses,-1)) - def _calc_median(self): - median_coefs = numpy.median(self.inversion_coefs, axis=1) - return _calc_slip(median_coefs) - - def _calc_median(self): + def _calc_median(self) -> numpy.ndarray: median_coefs = numpy.median(self.inversion_coefs, axis=1) return self._calc_slip(median_coefs) - def _calc_stddev(self): + def _calc_stddev(self) -> numpy.ndarray: median_coefs = numpy.median(self.inversion_coefs, axis=1) stddev_coefs = numpy.std(self.inversion_coefs, axis=1) low_slip = self._calc_slip(median_coefs - stddev_coefs) high_slip = self._calc_slip(median_coefs + stddev_coefs) return (low_slip, high_slip) - def _calc_minmax(self): + def _calc_minmax(self) -> tuple: min_coefs = numpy.min(self.inversion_coefs, axis=1) max_coefs = numpy.max(self.inversion_coefs, axis=1) min_slip = self._calc_slip(min_coefs) max_slip = self._calc_slip(max_coefs) return (min_slip, max_slip) - def _calc_slip(self, slip_coefs): + def _calc_slip(self, slip_coefs) -> numpy.ndarray: nimpulses, nlocs = self.impulse_slip.shape[0:2] slip = numpy.dot(slip_coefs.reshape((1,-1)), self.impulse_slip.reshape((nimpulses,-1))) return slip.reshape((nlocs,-1)) - def plot(self, median, stddev, minmax, filename, show_plot): - figure = pyplot.figure(figsize=(7.0, 4.0), dpi=150, layout="tight") + def _plot(self, median: numpy.ndarray, stddev: numpy.ndarray, minmax: tuple, filename: str, show_plot: bool): + figure = pyplot.figure(figsize=(6.5, 4.0), layout="tight") axes = self._setup_axes(figure) self._plot_observed(axes) self._plot_median(median, axes) @@ -118,20 +145,22 @@ def plot(self, median, stddev, minmax, filename, show_plot): def _setup_axes(self, figure): axes = figure.add_subplot() - axes.set_xlabel("Distance along Strike (km)") - axes.set_ylabel("Left-lateral Slip (m)") + axes.spines[["top", "right"]].set_visible(False) + axes.set_xlabel("Distance along strike, km") + axes.set_ylabel("Left-lateral slip, m") axes.set_xlim(self.XLIM) return axes def _plot_observed(self, axes): - axes.plot(self.observed_coords[:,1] / 1.0e+3, self.observed_slip[:,1], linewidth=2, color=self.COLOR_OBSERVED, label="Observed") + axes.plot(self.raw_y / km.value, self.raw_slip, linewidth=2, color=self.COLOR_RAW, label="Exact slip") + axes.plot(self.observed_coords[:,1] / 1.0e+3, self.observed_slip[:,1], linewidth=2, color=self.COLOR_OBSERVED, label="Prescribed slip") def _plot_median(self, median, axes): axes.plot(self.impulse_coords[:,1] / 1.0e+3, median[:,1], linewidth=2, color=self.COLOR_MODEL, label="Model median") def _plot_stddev(self, stddev, axes): low_slip, high_slip = stddev - axes.fill_between(self.impulse_coords[:,1] / 1.0e+3, low_slip[:,1], high_slip[:,1], color=self.COLOR_MODEL, lw=0, alpha=0.5, label="median$\pm$stddev") + axes.fill_between(self.impulse_coords[:,1] / 1.0e+3, low_slip[:,1], high_slip[:,1], color=self.COLOR_MODEL, lw=0, alpha=0.5, label=r"median$\pm$stddev") def _plot_minmax(self, minmax, axes): min_slip, max_slip = minmax @@ -139,24 +168,6 @@ def _plot_minmax(self, minmax, axes): axes.plot(self.impulse_coords[:,1] / 1.0e+3, max_slip[:,1], color=self.COLOR_MODEL, lw=1, ls="--", label=None) - -def cli(): - parser = argparse.ArgumentParser() - parser.add_argument("--catmip-theta", action="store", dest="filename_theta", type=str, required=True, help="Filename of output from CATMIP inversion.") - parser.add_argument("--pylith-impulses", action="store", dest="filename_impulses", type=str, default=PlotApp.FILENAME_IMPULSES, help="Name of HDF5 file with fault slip inpulses from Green's functions simulation.") - parser.add_argument("--observed", action="store", dest="filename_observed", type=str, default=PlotApp.FILENAME_OBSERVED, help="Name of HDF5 file with station output with fake observations.") - parser.add_argument("--no-gui", action="store_false", dest="show_plot", default=True, help="Do not display plot.") - - args = parser.parse_args() - kwargs = { - "filename_observed": args.filename_observed, - "filename_impulses": args.filename_impulses, - "filename_theta": args.filename_theta, - "show_plot": args.show_plot, - } - PlotApp().main(**kwargs) - - if __name__ == "__main__": cli() diff --git a/examples/strikeslip-2d/viz/plot_inversion_results.py b/examples/strikeslip-2d/viz/plot_inversion_results.py index fa9b09de9e..2f73b88762 100755 --- a/examples/strikeslip-2d/viz/plot_inversion_results.py +++ b/examples/strikeslip-2d/viz/plot_inversion_results.py @@ -4,34 +4,69 @@ the true solution for Step 6. """ -# Import argparse and re (regular expressions) Python modules(standard Python) -import argparse +# Standard Python modules +import pathlib import re -# Import numpy, h5py, and matplotlib modules (included in PyLith binary installation) +# Python modules included in PyLith binary installation import numpy import h5py import matplotlib.pyplot as pyplot +from pythia.pyre.units.length import km -class PlotApp: +OUTPUT_DIR = pathlib.Path("output") - FILENAME_MODELS = "output/step06_inversion-results.txt" - FILENAME_OBSERVED = "output/step04_varslip-fault.h5" +FILENAME_MODELS = OUTPUT_DIR / "step06_inversion-results.txt" +FILENAME_PRESCRIBED = OUTPUT_DIR / "step04_varslip-fault.h5" +FILENAME_RAW = OUTPUT_DIR / "slip_variable.txt" - - COLOR_OBSERVED = "black" - COLORS_MODEL = ["orange", "purple", "red", "cyan"] - XLIM = (-25.0, 25.0) - def main(self, filename_models: str=None, filename_observed: str=None, observed_only: bool=False, show_plot: bool=False): - filename_models = filename_models or self.FILENAME_MODELS - filename_observed = filename_observed or self.FILENAME_OBSERVED +def cli(): + """Command line interface. + """ + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument("--models", action="store", dest="filename_models", type=str, default=FILENAME_MODELS, help="Filename of output from inversion.") + parser.add_argument("--prescribed-slip", action="store", dest="filename_prescribed", type=str, default=FILENAME_PRESCRIBED, help="Name of HDF5 file for fault from prescribed slip simulation.") + parser.add_argument("--raw", action="store", dest="filename_raw", type=str, default=FILENAME_RAW, help="Name of ASCII file with raw slip profile.") + parser.add_argument("--prescribed-only", action="store_true", dest="prescribed_only", default=False, help="Show only observed slip.") + parser.add_argument("--no-gui", action="store_false", dest="show_plot", default=True, help="Do not display plot.") + + args = parser.parse_args() + PlotApp().run( + filename_models=args.filename_models, + filename_prescribed=args.filename_prescribed, + filename_raw=args.filename_raw, + prescribed_only=args.prescribed_only, + show_plot=args.show_plot, + ) - self.load_observed(filename_observed) - self.load_models(filename_models) - self.plot(observed_only, show_plot) - def load_observed(self, filename): +class PlotApp: + """Application to plot simple inversion results. + """ + COLOR_RAW = "gray" + COLOR_OBSERVED = "black" + COLORS_MODEL = ["orange", "purple", "red", "cyan"] + XLIM = (-25.0, 35.0) + + def run(self, filename_models: str, filename_prescribed: str, filename_raw: str, prescribed_only: bool, show_plot: bool): + """Run plotting application. + + Args: + filename_models: Filename of output from inversion. + filename_prescribed: Name of HDF5 file with station output with fake observations. + filename_raw: Name of ASCII file with raw slip values. + prescribed_only: Show only prescribed slip. + show_plot: Show plot window. + """ + self._load_observed(filename_prescribed) + self._load_raw(filename_raw) + self._load_models(filename_models) + self._plot(prescribed_only, show_plot) + + def _load_observed(self, filename: str): h5 = h5py.File(filename, "r") observed_coords = h5['geometry/vertices'][:] observed_slip = h5['vertex_fields/slip'][:,:,:].squeeze() @@ -42,7 +77,12 @@ def load_observed(self, filename): self.observed_coords = observed_coords[reorder,:] self.observed_slip = observed_slip[reorder,:] - def load_models(self, filename): + def _load_raw(self, filename: str): + data = numpy.loadtxt(filename) + self.raw_y = data[:,0] + self.raw_slip = data[:,1] + + def _load_models(self, filename: str): model = numpy.loadtxt(filename) self.models_coords = model[:,0] self.models_slip = model[:, 1:] @@ -51,53 +91,38 @@ def load_models(self, filename): header = fin.readline() self.models_penalty = list(map(float, re.findall(r"penalty=(\d+\.\d+)", header))) - def plot(self, observed_only=False, show_plot=False): - figure = pyplot.figure(figsize=(7.0, 3.5), dpi=150, layout="tight") + def _plot(self, prescribed_only=False, show_plot=False): + figure = pyplot.figure(figsize=(6.5, 4.0), layout="tight") axes = self._setup_axes(figure) self._plot_observed(axes) - if not observed_only: + if not prescribed_only: self._plot_models(axes) axes.legend(loc="upper right") if show_plot: pyplot.show() - filename = "step04-slip.pdf" if observed_only else "step06_inversion-results.pdf" + filename = "step04-slip.pdf" if prescribed_only else "step06_inversion-results.pdf" figure.savefig(filename) def _setup_axes(self, figure): axes = figure.add_subplot() - axes.set_xlabel("Distance along Strike (km)") - axes.set_ylabel("Left-lateral Slip (m)") + axes.spines[["top", "right"]].set_visible(False) + axes.set_xlabel("Distance along strike, km") + axes.set_ylabel("Left-lateral slip, m") axes.set_xlim(self.XLIM) return axes def _plot_observed(self, axes): - axes.plot(self.observed_coords[:,1] / 1.0e+3, self.observed_slip[:,1], linewidth=2, color=self.COLOR_OBSERVED, label="Observed") + axes.plot(self.raw_y / km.value, self.raw_slip, linewidth=2, color=self.COLOR_RAW, label="Exact slip") + axes.plot(self.observed_coords[:,1] / km.value, self.observed_slip[:,1], linewidth=2, color=self.COLOR_OBSERVED, label="Prescribed slip") def _plot_models(self, axes): nmodels = self.models_slip.shape[1] - coords = self.models_coords / 1.0e+3 + coords = self.models_coords / km.value for imodel in range(nmodels): axes.plot(coords, self.models_slip[:,imodel], "--", color=self.COLORS_MODEL[imodel], label=f"Model penalty={self.models_penalty[imodel]}") -def cli(): - parser = argparse.ArgumentParser() - parser.add_argument("--models", action="store", dest="filename_models", type=str, default=PlotApp.FILENAME_MODELS, help="Filename of output from inversion.") - parser.add_argument("--observed", action="store", dest="filename_observed", type=str, default=PlotApp.FILENAME_OBSERVED, help="Name of HDF5 file with station output with fake observations.") - parser.add_argument("--observed-only", action="store_true", dest="observed_only", default=False, help="Show only observed slip.") - parser.add_argument("--no-gui", action="store_false", dest="show_plot", default=True, help="Do not display plot.") - - args = parser.parse_args() - kwargs = { - "filename_models": args.filename_models, - "filename_observed": args.filename_observed, - "observed_only": args.observed_only, - "show_plot": args.show_plot, - } - PlotApp().main(**kwargs) - - if __name__ == "__main__": cli() diff --git a/examples/strikeslip-2d/viz/plot_slip.py b/examples/strikeslip-2d/viz/plot_slip.py new file mode 100755 index 0000000000..9b1744f14d --- /dev/null +++ b/examples/strikeslip-2d/viz/plot_slip.py @@ -0,0 +1,25 @@ +#!/usr/bin/env nemesis + +from matplotlib import pyplot +import numpy + +from pythia.pyre.units.length import km + +FILENAME = "output/slip_variable.txt" + +data = numpy.loadtxt(FILENAME) +x = data[:,0] / km.value +slip = data[:,1] + +figure = pyplot.figure(figsize=(6.5, 4.0)) +ax = figure.add_axes(rect=(0.09, 0.12, 0.88, 0.86)) +ax.spines[["top", "right"]].set_visible(False) +ax.plot(x, slip) +ax.set_xlabel("Distance along strike, km") +ax.set_ylabel("Left-lateral slip, m") + +ax.set_xlim(numpy.min(x), numpy.max(x)) +ax.set_ylim(0, 1.1*numpy.max(slip)) + +pyplot.savefig("step04-slip.pdf") +print(f"Mean slip: {numpy.mean(data[:,1]):.2f}") diff --git a/examples/strikeslip-2d/viz/plot_slip_impulses.py b/examples/strikeslip-2d/viz/plot_slip_impulses.py new file mode 100755 index 0000000000..2a7e4be564 --- /dev/null +++ b/examples/strikeslip-2d/viz/plot_slip_impulses.py @@ -0,0 +1,35 @@ +#!/usr/bin/env nemesis + +import pathlib + +from matplotlib import pyplot + +import numpy +import h5py + + +OUTPUT_DIR = pathlib.Path("output") +SIM_NAME = "step05_greensfns" + +filename = f"{SIM_NAME}-fault.h5" + +h5 = h5py.File(OUTPUT_DIR / filename) +y = h5["geometry/vertices"][:,1] +slip = h5["vertex_fields/slip"][:,:,1] +n_impulses = slip.shape[0] + +i_sort = numpy.argsort(y) +y = y[i_sort] / 1.0e+3 +slip = slip[:,i_sort] + +figure = pyplot.figure(figsize=(6.5, 4.0)) +ax = figure.add_axes(rect=(0.09, 0.13, 0.9, 0.83)) +ax.spines[["top", "right"]].set_visible(False) +for i_impulse in range(n_impulses): + ax.plot(y, slip[i_impulse,:]) +ax.set_xlabel("Distance along strike, km") +ax.set_ylabel("Impulse amplitude, m") +ax.autoscale(axis="both", enable=True, tight="True") + +pyplot.savefig(f"{SIM_NAME}-impulses.pdf") +pyplot.close(figure) diff --git a/libsrc/pylith/Makefile.am b/libsrc/pylith/Makefile.am index d56d79d51f..b790dd9349 100644 --- a/libsrc/pylith/Makefile.am +++ b/libsrc/pylith/Makefile.am @@ -136,6 +136,7 @@ libpylith_la_SOURCES = \ topology/Distributor.cc \ topology/ReverseCuthillMcKee.cc \ topology/RefineUniform.cc \ + topology/RefineInterpolator.cc \ utils/EventLogger.cc \ utils/PyreComponent.cc \ utils/GenericComponent.cc \ diff --git a/libsrc/pylith/faults/TopologyOps.cc b/libsrc/pylith/faults/TopologyOps.cc index 11560d6e49..13402d8bed 100644 --- a/libsrc/pylith/faults/TopologyOps.cc +++ b/libsrc/pylith/faults/TopologyOps.cc @@ -222,7 +222,7 @@ pylith::faults::TopologyOps::create(pylith::topology::Mesh* mesh, err = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);PYLITH_CHECK_ERROR(err); err = DMPlexSetScale(sdm, PETSC_UNIT_LENGTH, lengthScale);PYLITH_CHECK_ERROR(err); err = DMViewFromOptions(sdm, NULL, "-pylith_cohesive_dm_view");PYLITH_CHECK_ERROR(err); - mesh->setDM(sdm); + mesh->setDM(sdm, "domain"); } // create diff --git a/libsrc/pylith/meshio/DataWriterHDF5.cc b/libsrc/pylith/meshio/DataWriterHDF5.cc index fe3830072c..b331620fd7 100644 --- a/libsrc/pylith/meshio/DataWriterHDF5.cc +++ b/libsrc/pylith/meshio/DataWriterHDF5.cc @@ -178,7 +178,7 @@ pylith::meshio::DataWriterHDF5::writeVertexField(const PylithScalar t, const int istep = _timesteps[name]; // Add time stamp to "/time" if necessary. MPI_Comm comm; - err = PetscObjectGetComm((PetscObject)subfield.getDM(), &comm);PYLITH_CHECK_ERROR(err); + err = PetscObjectGetComm((PetscObject)subfield.getOutputDM(), &comm);PYLITH_CHECK_ERROR(err); PetscMPIInt commRank; err = MPI_Comm_rank(comm, &commRank);PYLITH_CHECK_ERROR(err); if (_tstampIndex == istep) { @@ -189,7 +189,7 @@ pylith::meshio::DataWriterHDF5::writeVertexField(const PylithScalar t, err = PetscViewerHDF5PushTimestepping(_viewer);PYLITH_CHECK_ERROR(err); err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err); - PetscVec vector = subfield.getVector();assert(vector); + PetscVec vector = subfield.getOutputVector();assert(vector); DataWriter::_writeVec(vector, _viewer);PYLITH_CHECK_ERROR(err); err = PetscViewerHDF5PopTimestepping(_viewer);PYLITH_CHECK_ERROR(err); err = PetscViewerHDF5PopGroup(_viewer);PYLITH_CHECK_ERROR(err); @@ -240,8 +240,7 @@ pylith::meshio::DataWriterHDF5::writeCellField(const PylithScalar t, } const int istep = _timesteps[name]; // Add time stamp to "/time" if necessary. - MPI_Comm comm; - err = PetscObjectGetComm((PetscObject)subfield.getDM(), &comm);PYLITH_CHECK_ERROR(err); + MPI_Comm comm = PetscObjectComm((PetscObject)subfield.getOutputDM()); PetscMPIInt commRank; err = MPI_Comm_rank(comm, &commRank);PYLITH_CHECK_ERROR(err); if (_tstampIndex == istep) { @@ -252,7 +251,7 @@ pylith::meshio::DataWriterHDF5::writeCellField(const PylithScalar t, err = PetscViewerHDF5PushTimestepping(_viewer);PYLITH_CHECK_ERROR(err); err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err); - PetscVec vector = subfield.getVector();assert(vector); + PetscVec vector = subfield.getOutputVector();assert(vector); DataWriter::_writeVec(vector, _viewer); err = PetscViewerHDF5PopTimestepping(_viewer);PYLITH_CHECK_ERROR(err); err = PetscViewerHDF5PopGroup(_viewer);PYLITH_CHECK_ERROR(err); diff --git a/libsrc/pylith/meshio/DataWriterHDF5Ext.cc b/libsrc/pylith/meshio/DataWriterHDF5Ext.cc index f9dbe3854f..e8d1863cf9 100644 --- a/libsrc/pylith/meshio/DataWriterHDF5Ext.cc +++ b/libsrc/pylith/meshio/DataWriterHDF5Ext.cc @@ -148,7 +148,7 @@ pylith::meshio::DataWriterHDF5Ext::writeVertexField(const PylithScalar t, const char* name = subfield.getDescription().label.c_str(); try { - PetscDM dmMesh = subfield.getDM();assert(dmMesh); + PetscDM dmMesh = subfield.getOutputDM();assert(dmMesh); PetscErrorCode err; MPI_Comm comm; @@ -176,7 +176,7 @@ pylith::meshio::DataWriterHDF5Ext::writeVertexField(const PylithScalar t, } // else assert(binaryViewer); - PetscVec vector = subfield.getVector();assert(vector); + PetscVec vector = subfield.getOutputVector();assert(vector); DataWriter::_writeVec(vector, binaryViewer); ExternalDataset& datasetInfo = _datasets[name]; @@ -288,7 +288,7 @@ pylith::meshio::DataWriterHDF5Ext::writeCellField(const PylithScalar t, try { PetscErrorCode err; - PetscDM dmMesh = subfield.getDM();assert(dmMesh); + PetscDM dmMesh = subfield.getOutputDM();assert(dmMesh); MPI_Comm comm; PetscMPIInt commRank; err = PetscObjectGetComm((PetscObject) dmMesh, &comm);PYLITH_CHECK_ERROR(err); @@ -314,7 +314,7 @@ pylith::meshio::DataWriterHDF5Ext::writeCellField(const PylithScalar t, } // else assert(binaryViewer); - PetscVec vector = subfield.getVector();assert(vector); + PetscVec vector = subfield.getOutputVector();assert(vector); DataWriter::_writeVec(vector, binaryViewer); ExternalDataset& datasetInfo = _datasets[name]; @@ -333,7 +333,7 @@ pylith::meshio::DataWriterHDF5Ext::writeCellField(const PylithScalar t, if (createdExternalDataset) { // Get cell information PetscSection section = NULL; - err = DMGetLocalSection(subfield.getDM(), §ion);assert(section); + err = DMGetLocalSection(subfield.getOutputDM(), §ion);assert(section); PetscInt dof = 0, numLocalCells = 0, numCells, cellHeight, cStart, cEnd; PetscIS globalCellNumbers; diff --git a/libsrc/pylith/meshio/DataWriterVTK.cc b/libsrc/pylith/meshio/DataWriterVTK.cc index 651749c220..0eee256e75 100644 --- a/libsrc/pylith/meshio/DataWriterVTK.cc +++ b/libsrc/pylith/meshio/DataWriterVTK.cc @@ -212,8 +212,8 @@ pylith::meshio::DataWriterVTK::writeVertexField(const PylithScalar t, PetscViewerVTKFieldType ft = subfield.getDescription().vectorFieldType == pylith::topology::FieldBase::VECTOR ? PETSC_VTK_POINT_VECTOR_FIELD : PETSC_VTK_POINT_FIELD; err = PetscViewerVTKAddField(_viewer, (PetscObject)_dm, DMPlexVTKWriteAll, PETSC_DEFAULT, ft, - PETSC_TRUE, (PetscObject)subfield.getVector());PYLITH_CHECK_ERROR(err); - err = PetscObjectReference((PetscObject) subfield.getVector());PYLITH_CHECK_ERROR(err); // Viewer destroys Vec + PETSC_TRUE, (PetscObject)subfield.getOutputVector());PYLITH_CHECK_ERROR(err); + err = PetscObjectReference((PetscObject) subfield.getOutputVector());PYLITH_CHECK_ERROR(err); // Viewer destroys Vec _wroteVertexHeader = true; @@ -234,8 +234,8 @@ pylith::meshio::DataWriterVTK::writeCellField(const PylithScalar t, PetscViewerVTKFieldType ft = subfield.getDescription().vectorFieldType == pylith::topology::FieldBase::VECTOR ? PETSC_VTK_CELL_VECTOR_FIELD : PETSC_VTK_CELL_FIELD; err = PetscViewerVTKAddField(_viewer, (PetscObject)_dm, DMPlexVTKWriteAll, PETSC_DEFAULT, ft, - PETSC_TRUE, (PetscObject)subfield.getVector());PYLITH_CHECK_ERROR(err); - err = PetscObjectReference((PetscObject) subfield.getVector());PYLITH_CHECK_ERROR(err); // Viewer destroys Vec + PETSC_TRUE, (PetscObject)subfield.getOutputVector());PYLITH_CHECK_ERROR(err); + err = PetscObjectReference((PetscObject) subfield.getOutputVector());PYLITH_CHECK_ERROR(err); // Viewer destroys Vec _wroteCellHeader = true; diff --git a/libsrc/pylith/meshio/MeshBuilder.cc b/libsrc/pylith/meshio/MeshBuilder.cc index dcd26039b9..404e9dd901 100644 --- a/libsrc/pylith/meshio/MeshBuilder.cc +++ b/libsrc/pylith/meshio/MeshBuilder.cc @@ -127,7 +127,7 @@ pylith::meshio::MeshBuilder::buildMesh(topology::Mesh* mesh, err = DMPlexInvertCell(ct, (int *) &cells[coff]);PYLITH_CHECK_ERROR(err); } err = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, &cells[0], spaceDim, &(*coordinates)[0], &dmMesh);PYLITH_CHECK_ERROR(err); - mesh->setDM(dmMesh); + mesh->setDM(dmMesh, "domain"); _MeshBuilder::Events::logger.eventEnd(_MeshBuilder::Events::buildMesh); PYLITH_METHOD_END; diff --git a/libsrc/pylith/meshio/MeshIOPetsc.cc b/libsrc/pylith/meshio/MeshIOPetsc.cc index d6c38058ad..f8268ba9cd 100644 --- a/libsrc/pylith/meshio/MeshIOPetsc.cc +++ b/libsrc/pylith/meshio/MeshIOPetsc.cc @@ -148,7 +148,7 @@ pylith::meshio::MeshIOPetsc::_read(void) { _MeshIOPetsc::fixMaterialLabel(&dmMesh); _MeshIOPetsc::fixBoundaryLabels(&dmMesh); - _mesh->setDM(dmMesh); + _mesh->setDM(dmMesh, "domain"); } catch (...) { DMDestroy(&dmMesh); throw; diff --git a/libsrc/pylith/meshio/OutputObserver.cc b/libsrc/pylith/meshio/OutputObserver.cc index 6aeff61cff..25ae96a226 100644 --- a/libsrc/pylith/meshio/OutputObserver.cc +++ b/libsrc/pylith/meshio/OutputObserver.cc @@ -67,6 +67,7 @@ pylith::meshio::_OutputObserver::Events::init(void) { // Constructor pylith::meshio::OutputObserver::OutputObserver(void) : _timeScale(1.0), + _outputMesh(NULL), _writer(NULL), _trigger(NULL), _outputBasisOrder(1) { @@ -95,6 +96,7 @@ pylith::meshio::OutputObserver::deallocate(void) { delete iter->second;iter->second = NULL; } // for _subfields.clear(); + delete _outputMesh;_outputMesh = NULL; _writer = NULL; // :TODO: Use shared pointer _trigger = NULL; // :TODO: Use shared pointer @@ -152,6 +154,26 @@ pylith::meshio::OutputObserver::setOutputBasisOrder(const int value) { } // setOutputBasisOrder +// ------------------------------------------------------------------------------------------------ + +// Set number of refinement levels for output. +void +pylith::meshio::OutputObserver::setRefineLevels(const int value) { + PYLITH_METHOD_BEGIN; + PYLITH_COMPONENT_DEBUG("OutputObserver::setRefineLevels(value="<setDM(dmOutput); + } // if + + PYLITH_METHOD_RETURN(_outputMesh); +} // _getOutputMesh + + // ------------------------------------------------------------------------------------------------ // Get output subfield, creating if necessary. pylith::meshio::OutputSubfield* @@ -176,7 +215,7 @@ pylith::meshio::OutputObserver::_getSubfield(const pylith::topology::Field& fiel _OutputObserver::Events::logger.eventBegin(_OutputObserver::Events::getSubfield); if (0 == _subfields.count(name) ) { - _subfields[name] = OutputSubfield::create(field, submesh, name, _outputBasisOrder); + _subfields[name] = OutputSubfield::create(field, submesh, name, _outputBasisOrder, _refineLevels); } // if _OutputObserver::Events::logger.eventEnd(_OutputObserver::Events::getSubfield); diff --git a/libsrc/pylith/meshio/OutputObserver.hh b/libsrc/pylith/meshio/OutputObserver.hh index c7279d3009..0a29fabf80 100644 --- a/libsrc/pylith/meshio/OutputObserver.hh +++ b/libsrc/pylith/meshio/OutputObserver.hh @@ -57,6 +57,12 @@ public: */ void setOutputBasisOrder(const int value); + /** Set number of mesh refinement levels for output. + * + * @param[in] value Number of mesh refinement levels for output. + */ + void setRefineLevels(const int value); + /** Set time scale. * * @param[in] value Time scale for dimensionalizing time. @@ -73,6 +79,13 @@ protected: */ void _setContext(const pylith::topology::Mesh & mesh); + /** Get mesh associated with subfield output. + * + * @param[in] subfield Subfield for output. + * @returns Mesh associated with output. + */ + pylith::topology::Mesh* _getOutputMesh(const pylith::meshio::OutputSubfield& subfield); + /** Get output subfield, creating if necessary. * * @param[in] field Field containing subfields. @@ -97,9 +110,11 @@ protected: PylithReal _timeScale; ///< Time scale for dimentionalizing time. std::map _subfields; ///< Subfields extracted for output. + pylith::topology::Mesh* _outputMesh; ///< Mesh associated with output.ß DataWriter* _writer; ///< Writer for data. OutputTrigger* _trigger; ///< Trigger for deciding how often to write output. int _outputBasisOrder; ///< Basis order for output. + int _refineLevels; ///< Number of mesh refinement levels for output. // NOT IMPLEMENTED //////////////////////////////////////////////////////////////////////////// private: diff --git a/libsrc/pylith/meshio/OutputPhysics.cc b/libsrc/pylith/meshio/OutputPhysics.cc index 597f5a1c44..e28eb7a7cb 100644 --- a/libsrc/pylith/meshio/OutputPhysics.cc +++ b/libsrc/pylith/meshio/OutputPhysics.cc @@ -268,13 +268,11 @@ pylith::meshio::OutputPhysics::_writeInfo(void) { assert(_physics); const pylith::topology::Field* auxiliaryField = _physics->getAuxiliaryField(); const pylith::topology::Field* diagnosticField = _physics->getDiagnosticField(); + const pylith::topology::Mesh& domainMesh = _physics->getPhysicsDomainMesh(); const pylith::string_vector& infoNames = _expandInfoFieldNames(auxiliaryField, diagnosticField); const bool isInfo = true; - const pylith::topology::Mesh& domainMesh = _physics->getPhysicsDomainMesh(); - _open(domainMesh, isInfo); - _openDataStep(0.0, domainMesh); if (auxiliaryField) { auxiliaryField->scatterLocalToOutput(); } PetscVec auxiliaryVector = (auxiliaryField) ? auxiliaryField->getOutputVector() : NULL; @@ -284,20 +282,28 @@ pylith::meshio::OutputPhysics::_writeInfo(void) { const size_t numInfoFields = infoNames.size(); for (size_t i = 0; i < numInfoFields; i++) { + OutputSubfield* subfield = NULL; if (auxiliaryField->hasSubfield(infoNames[i].c_str())) { - OutputSubfield* subfield = _getSubfield(*auxiliaryField, domainMesh, infoNames[i].c_str()); + subfield = _getSubfield(*auxiliaryField, domainMesh, infoNames[i].c_str()); subfield->project(auxiliaryVector); - _appendField(0.0, *subfield); } else if (diagnosticField->hasSubfield(infoNames[i].c_str())) { - OutputSubfield* subfield = _getSubfield(*diagnosticField, domainMesh, infoNames[i].c_str()); + subfield = _getSubfield(*diagnosticField, domainMesh, infoNames[i].c_str()); subfield->project(diagnosticVector); - _appendField(0.0, *subfield); } else { std::ostringstream msg; msg << "Internal Error: Could not find subfield '" << infoNames[i] << "' for info output."; PYLITH_COMPONENT_ERROR(msg.str()); throw std::runtime_error(msg.str()); } // if/else + + if (0 == i) { + // Need output mesh from subfield (which may be refined). + assert(subfield); + pylith::topology::Mesh* outputMesh = _getOutputMesh(*subfield); + _open(*outputMesh, isInfo); + _openDataStep(0.0, *outputMesh); + } // if + OutputObserver::_appendField(0.0, *subfield); } // for _closeDataStep(); @@ -395,8 +401,6 @@ pylith::meshio::OutputPhysics::_writeDataStep(const PylithReal t, const pylith::string_vector& dataNames = _expandDataFieldNames(solution, auxiliaryField, derivedField); - _openDataStep(t, domainMesh); - if (auxiliaryField) { auxiliaryField->scatterLocalToOutput(); } PetscVec auxiliaryVector = (auxiliaryField) ? auxiliaryField->getOutputVector() : NULL; @@ -429,6 +433,12 @@ pylith::meshio::OutputPhysics::_writeDataStep(const PylithReal t, throw std::runtime_error(msg.str()); } // if/else + if (0 == i) { + // Need output mesh from subfield (which may be refined). + assert(subfield); + pylith::topology::Mesh* outputMesh = _getOutputMesh(*subfield); + _openDataStep(t, *outputMesh); + } // if OutputObserver::_appendField(t, *subfield); } // for _closeDataStep(); diff --git a/libsrc/pylith/meshio/OutputSolnBoundary.cc b/libsrc/pylith/meshio/OutputSolnBoundary.cc index 54677772a1..262bfd602b 100644 --- a/libsrc/pylith/meshio/OutputSolnBoundary.cc +++ b/libsrc/pylith/meshio/OutputSolnBoundary.cc @@ -116,7 +116,6 @@ pylith::meshio::OutputSolnBoundary::_writeSolnStep(const PylithReal t, const pylith::string_vector& subfieldNames = pylith::topology::FieldOps::getSubfieldNamesDomain(solution); PetscVec solutionVector = solution.getOutputVector();assert(solutionVector); - _openSolnStep(t, *_boundaryMesh); const size_t numSubfieldNames = subfieldNames.size(); for (size_t iField = 0; iField < numSubfieldNames; iField++) { assert(solution.hasSubfield(subfieldNames[iField].c_str())); @@ -125,6 +124,12 @@ pylith::meshio::OutputSolnBoundary::_writeSolnStep(const PylithReal t, subfield = OutputObserver::_getSubfield(solution, *_boundaryMesh, subfieldNames[iField].c_str());assert(subfield); subfield->project(solutionVector); + if (0 == iField) { + // Need output mesh from subfield (which may be refined). + assert(subfield); + pylith::topology::Mesh* outputMesh = _getOutputMesh(*subfield); + _openSolnStep(t, *outputMesh); + } // if OutputObserver::_appendField(t, *subfield); } // for _closeSolnStep(); diff --git a/libsrc/pylith/meshio/OutputSolnDomain.cc b/libsrc/pylith/meshio/OutputSolnDomain.cc index c6a59d2bad..309b5f2ceb 100644 --- a/libsrc/pylith/meshio/OutputSolnDomain.cc +++ b/libsrc/pylith/meshio/OutputSolnDomain.cc @@ -44,7 +44,6 @@ pylith::meshio::OutputSolnDomain::_writeSolnStep(const PylithReal t, const pylith::string_vector& subfieldNames = pylith::topology::FieldOps::getSubfieldNamesDomain(solution); PetscVec solutionVector = solution.getOutputVector();assert(solutionVector); - _openSolnStep(t, solution.getMesh()); const size_t numSubfieldNames = subfieldNames.size(); for (size_t iField = 0; iField < numSubfieldNames; iField++) { assert(solution.hasSubfield(subfieldNames[iField].c_str())); @@ -53,6 +52,12 @@ pylith::meshio::OutputSolnDomain::_writeSolnStep(const PylithReal t, subfield = OutputObserver::_getSubfield(solution, solution.getMesh(), subfieldNames[iField].c_str());assert(subfield); subfield->project(solutionVector); + if (0 == iField) { + // Need output mesh from subfield (which may be refined). + assert(subfield); + pylith::topology::Mesh* outputMesh = _getOutputMesh(*subfield); + _openSolnStep(t, *outputMesh); + } // if OutputObserver::_appendField(t, *subfield); } // for _closeSolnStep(); diff --git a/libsrc/pylith/meshio/OutputSubfield.cc b/libsrc/pylith/meshio/OutputSubfield.cc index bfe71e6185..a65015edb5 100644 --- a/libsrc/pylith/meshio/OutputSubfield.cc +++ b/libsrc/pylith/meshio/OutputSubfield.cc @@ -16,6 +16,7 @@ #include "pylith/topology/Mesh.hh" // USES Mesh #include "pylith/topology/FieldOps.hh" // USES FieldOps #include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh +#include "pylith/topology/RefineInterpolator.hh" // USES RefineInterpolator #include "pylith/fekernels/Solution.hh" // USES Solution::passThruSubfield #include "pylith/utils/error.hh" // USES PYLITH_CHECK_ERROR @@ -75,9 +76,14 @@ pylith::meshio::_OutputSubfield::Events::init(void) { // ------------------------------------------------------------------------------------------------ // Constructor pylith::meshio::OutputSubfield::OutputSubfield(void) : - _dm(NULL), - _vector(NULL), + _subfieldIndex(-1), + _projectDM(PETSC_NULLPTR), + _projectVector(PETSC_NULLPTR), + _projectVectorInterp(PETSC_NULLPTR), _fn(pylith::fekernels::Solution::passThruSubfield), + _outputDM(PETSC_NULLPTR), + _outputVector(PETSC_NULLPTR), + _interpolator(NULL), _label(NULL), _labelValue(0) { _OutputSubfield::Events::init(); @@ -95,9 +101,14 @@ pylith::meshio::OutputSubfield::~OutputSubfield(void) { // Deallocate PETSc and local data structures. void pylith::meshio::OutputSubfield::deallocate(void) { + delete _interpolator;_interpolator = NULL; + PetscErrorCode err; - err = DMDestroy(&_dm);PYLITH_CHECK_ERROR(err); - err = VecDestroy(&_vector);PYLITH_CHECK_ERROR(err); + err = VecDestroy(&_projectVector);PYLITH_CHECK_ERROR(err); + err = VecDestroy(&_projectVectorInterp);PYLITH_CHECK_ERROR(err); + err = DMDestroy(&_projectDM);PYLITH_CHECK_ERROR(err); + err = VecDestroy(&_outputVector);PYLITH_CHECK_ERROR(err); + err = DMDestroy(&_outputDM);PYLITH_CHECK_ERROR(err); _label = NULL; // Destroyed by DMDestroy() } // deallocate @@ -109,7 +120,8 @@ pylith::meshio::OutputSubfield* pylith::meshio::OutputSubfield::create(const pylith::topology::Field& field, const pylith::topology::Mesh& mesh, const char* name, - const int basisOrder) { + const int basisOrder, + const int refineLevels) { PYLITH_METHOD_BEGIN; // _OutputSubfield::Events::logger.eventBegin(_OutputSubfield::Events::create); @@ -118,27 +130,56 @@ pylith::meshio::OutputSubfield::create(const pylith::topology::Field& field, const pylith::topology::Field::SubfieldInfo& info = field.getSubfieldInfo(name); subfield->_subfieldIndex = info.index; subfield->_description = info.description; - subfield->_discretization = info.fe; - subfield->_discretization.dimension = mesh.getDimension(); - // Basis order of output should be less than or equai to the basis order of the computed field. - subfield->_discretization.basisOrder = std::min(basisOrder, info.fe.basisOrder); + const int outputBasisOrder = std::min(basisOrder, info.fe.basisOrder); + + // Discretization for projection + pylith::topology::FieldBase::Discretization projectDiscretization = info.fe; + projectDiscretization.dimension = mesh.getDimension(); + projectDiscretization.basisOrder = refineLevels ? info.fe.basisOrder : outputBasisOrder; + + // Discretization for output + subfield->_discretization = projectDiscretization; + subfield->_discretization.basisOrder = outputBasisOrder; PetscErrorCode err = PETSC_SUCCESS; - err = DMClone(mesh.getDM(), &subfield->_dm);PYLITH_CHECK_ERROR(err); - err = DMReorderSectionSetDefault(subfield->_dm, DM_REORDER_DEFAULT_FALSE);PYLITH_CHECK_ERROR(err); - err = DMReorderSectionSetType(subfield->_dm, NULL);PYLITH_CHECK_ERROR(err); - err = PetscObjectSetName((PetscObject)subfield->_dm, name);PYLITH_CHECK_ERROR(err); - - PetscFE fe = pylith::topology::FieldOps::createFE(subfield->_discretization, subfield->_dm, - info.description.numComponents);assert(fe); - err = PetscFESetName(fe, info.description.label.c_str());PYLITH_CHECK_ERROR(err); - err = DMSetField(subfield->_dm, 0, NULL, (PetscObject)fe);PYLITH_CHECK_ERROR(err); - err = DMSetFieldAvoidTensor(subfield->_dm, 0, PETSC_TRUE);PYLITH_CHECK_ERROR(err); - err = PetscFEDestroy(&fe);PYLITH_CHECK_ERROR(err); - err = DMCreateDS(subfield->_dm);PYLITH_CHECK_ERROR(err); - - err = DMCreateGlobalVector(subfield->_dm, &subfield->_vector);PYLITH_CHECK_ERROR(err); - err = PetscObjectSetName((PetscObject)subfield->_vector, name);PYLITH_CHECK_ERROR(err); + + // Setup PETSc DM for projection + err = DMClone(mesh.getDM(), &subfield->_projectDM);PYLITH_CHECK_ERROR(err); + err = PetscObjectSetName((PetscObject)subfield->_projectDM, name);PYLITH_CHECK_ERROR(err); + err = DMReorderSectionSetDefault(subfield->_projectDM, DM_REORDER_DEFAULT_FALSE);PYLITH_CHECK_ERROR(err); + err = DMReorderSectionSetType(subfield->_projectDM, NULL);PYLITH_CHECK_ERROR(err); + err = DMPlexReorderSetDefault(subfield->_projectDM, DM_REORDER_DEFAULT_FALSE); + + // Setup PETSc FE (discretization) for projection + PetscFE projectFE = pylith::topology::FieldOps::createFE(projectDiscretization, subfield->_projectDM, + info.description.numComponents);assert(projectFE); + err = PetscFESetName(projectFE, info.description.label.c_str());PYLITH_CHECK_ERROR(err); + err = DMSetField(subfield->_projectDM, 0, NULL, (PetscObject)projectFE);PYLITH_CHECK_ERROR(err); + err = DMSetFieldAvoidTensor(subfield->_projectDM, 0, PETSC_TRUE);PYLITH_CHECK_ERROR(err); + err = PetscFEDestroy(&projectFE);PYLITH_CHECK_ERROR(err); + err = DMCreateDS(subfield->_projectDM);PYLITH_CHECK_ERROR(err); + + if (!refineLevels) { + subfield->_outputDM = subfield->_projectDM; + err = PetscObjectReference((PetscObject)subfield->_outputDM);PYLITH_CHECK_ERROR(err); + } else { + delete subfield->_interpolator;subfield->_interpolator = new pylith::topology::RefineInterpolator(); + assert(subfield->_interpolator); + subfield->_interpolator->initialize(subfield->_projectDM, refineLevels, outputBasisOrder, info.description, subfield->_discretization); + subfield->_outputDM = subfield->_interpolator->getOutputDM(); + err = PetscObjectReference((PetscObject)subfield->_outputDM);PYLITH_CHECK_ERROR(err); + } // if/else + + err = DMCreateGlobalVector(subfield->_projectDM, &subfield->_projectVector);PYLITH_CHECK_ERROR(err); + err = PetscObjectSetName((PetscObject)subfield->_projectVector, name);PYLITH_CHECK_ERROR(err); + if (refineLevels) { + err = DMCreateGlobalVector(subfield->_outputDM, &subfield->_outputVector);PYLITH_CHECK_ERROR(err); + err = PetscObjectSetName((PetscObject)subfield->_outputVector, name);PYLITH_CHECK_ERROR(err); + err = VecDuplicate(subfield->_projectVector, &subfield->_projectVectorInterp);PYLITH_CHECK_ERROR(err); + } else { + subfield->_outputVector = subfield->_projectVector; + err = PetscObjectReference((PetscObject)subfield->_outputVector);PYLITH_CHECK_ERROR(err); + } // _OutputSubfield::Events::logger.eventEnd(_OutputSubfield::Events::create); PYLITH_METHOD_RETURN(subfield); @@ -160,11 +201,11 @@ pylith::meshio::OutputSubfield::create(const pylith::topology::Field& field, subfield->_subfieldIndex = info.index; subfield->_description = info.description; - PetscErrorCode err; - err = DMClone(mesh.getDM(), &subfield->_dm);PYLITH_CHECK_ERROR(err); - err = DMReorderSectionSetDefault(subfield->_dm, DM_REORDER_DEFAULT_FALSE);PYLITH_CHECK_ERROR(err); - err = DMReorderSectionSetType(subfield->_dm, NULL);PYLITH_CHECK_ERROR(err); - err = PetscObjectSetName((PetscObject)subfield->_dm, name);PYLITH_CHECK_ERROR(err); + PetscErrorCode err = PETSC_SUCCESS; + err = DMClone(mesh.getDM(), &subfield->_projectDM);PYLITH_CHECK_ERROR(err);assert(subfield->_projectDM); + err = DMReorderSectionSetDefault(subfield->_projectDM, DM_REORDER_DEFAULT_FALSE);PYLITH_CHECK_ERROR(err); + err = DMReorderSectionSetType(subfield->_projectDM, NULL);PYLITH_CHECK_ERROR(err); + err = PetscObjectSetName((PetscObject)subfield->_projectDM, name);PYLITH_CHECK_ERROR(err); pylith::topology::VecVisitorMesh fieldVisitor(field, name); @@ -178,10 +219,15 @@ pylith::meshio::OutputSubfield::create(const pylith::topology::Field& field, err = PetscSectionSetDof(subfieldSection, point, numDof);PYLITH_CHECK_ERROR(err); offset += numDof; } // for - err = DMSetLocalSection(subfield->_dm, subfieldSection);PYLITH_CHECK_ERROR(err); + err = DMSetLocalSection(subfield->_projectDM, subfieldSection);PYLITH_CHECK_ERROR(err); err = PetscSectionDestroy(&subfieldSection);PYLITH_CHECK_ERROR(err); - err = DMCreateGlobalVector(subfield->_dm, &subfield->_vector);PYLITH_CHECK_ERROR(err); - err = PetscObjectSetName((PetscObject)subfield->_vector, name);PYLITH_CHECK_ERROR(err); + err = DMCreateGlobalVector(subfield->_projectDM, &subfield->_projectVector);PYLITH_CHECK_ERROR(err); + err = PetscObjectSetName((PetscObject)subfield->_projectVector, name);PYLITH_CHECK_ERROR(err); + + subfield->_outputDM = subfield->_projectDM; + err = PetscObjectReference((PetscObject) subfield->_outputDM);PYLITH_CHECK_ERROR(err); + subfield->_outputVector = subfield->_projectVector; + err = PetscObjectReference((PetscObject) subfield->_outputVector);PYLITH_CHECK_ERROR(err); _OutputSubfield::Events::logger.eventEnd(_OutputSubfield::Events::createBasisOrder); PYLITH_METHOD_RETURN(subfield); @@ -200,8 +246,8 @@ pylith::meshio::OutputSubfield::setLabel(const char* name, PYLITH_METHOD_END; } // if PetscErrorCode err; - err = DMGetLabel(_dm, name, &_label);PYLITH_CHECK_ERROR(err); - err = DMPlexLabelComplete(_dm, _label); + err = DMGetLabel(_projectDM, name, &_label);PYLITH_CHECK_ERROR(err); + err = DMPlexLabelComplete(_projectDM, _label); _labelValue = value; @@ -229,16 +275,16 @@ pylith::meshio::OutputSubfield::getBasisOrder(void) const { // ------------------------------------------------------------------------------------------------ // Get filtered PETSc global vector. PetscVec -pylith::meshio::OutputSubfield::getVector(void) const { - return _vector; +pylith::meshio::OutputSubfield::getOutputVector(void) const { + return _outputVector; } // ------------------------------------------------------------------------------------------------ // Get PETSc DM for filtered vector. PetscDM -pylith::meshio::OutputSubfield::getDM(void) const { - return _dm; +pylith::meshio::OutputSubfield::getOutputDM(void) const { + return _outputDM; } @@ -249,13 +295,18 @@ pylith::meshio::OutputSubfield::project(const PetscVec& fieldVector) { PYLITH_METHOD_BEGIN; _OutputSubfield::Events::logger.eventBegin(_OutputSubfield::Events::project); assert(fieldVector); - assert(_vector); + assert(_projectVector); + assert(_outputVector); PetscErrorCode err; - const PetscReal t = PetscReal(_subfieldIndex) + 0.01; // :KLUDGE: Easiest way to get subfield to extract into fn. + const PetscReal t = PetscReal(_subfieldIndex) + 0.01; // :KLUDGE: Easiest way to get subfield to extract into fn - err = DMProjectField(_dm, t, fieldVector, &_fn, INSERT_VALUES, _vector);PYLITH_CHECK_ERROR(err); - err = VecScale(_vector, _description.scale);PYLITH_CHECK_ERROR(err); + err = DMProjectField(_projectDM, t, fieldVector, &_fn, INSERT_VALUES, _projectVector);PYLITH_CHECK_ERROR(err); + if (_interpolator) { + pylith::topology::FieldOps::transformVector(&_projectVectorInterp, _interpolator->getInputDM(), _projectVector, _projectDM); + _interpolator->interpolate(&_outputVector, _projectVectorInterp); + } // if + err = VecScale(_outputVector, _description.scale);PYLITH_CHECK_ERROR(err); _OutputSubfield::Events::logger.eventEnd(_OutputSubfield::Events::project); PYLITH_METHOD_END; @@ -269,14 +320,19 @@ pylith::meshio::OutputSubfield::projectWithLabel(const PetscVec& fieldVector) { PYLITH_METHOD_BEGIN; _OutputSubfield::Events::logger.eventBegin(_OutputSubfield::Events::projectWithLabel); assert(fieldVector); - assert(_vector); + assert(_projectVector); + assert(_outputVector); assert(_label); PetscErrorCode err; - const PetscReal t = PetscReal(_subfieldIndex) + 0.01; // :KLUDGE: Easiest way to get subfield to extract into fn. + const PetscReal t = PetscReal(_subfieldIndex) + 0.01; // :KLUDGE: Easiest way to get subfield to extract into fn - err = DMProjectFieldLabel(_dm, t, _label, 1, &_labelValue, PETSC_DETERMINE, NULL, fieldVector, &_fn, INSERT_VALUES, _vector);PYLITH_CHECK_ERROR(err); - err = VecScale(_vector, _description.scale);PYLITH_CHECK_ERROR(err); + err = DMProjectFieldLabel(_projectDM, t, _label, 1, &_labelValue, PETSC_DETERMINE, NULL, fieldVector, &_fn, INSERT_VALUES, _projectVector);PYLITH_CHECK_ERROR(err); + if (_interpolator) { + pylith::topology::FieldOps::transformVector(&_projectVectorInterp, _interpolator->getInputDM(), _projectVector, _projectDM); + _interpolator->interpolate(&_outputVector, _projectVectorInterp); + } // if + err = VecScale(_outputVector, _description.scale);PYLITH_CHECK_ERROR(err); _OutputSubfield::Events::logger.eventEnd(_OutputSubfield::Events::projectWithLabel); PYLITH_METHOD_END; @@ -297,7 +353,7 @@ pylith::meshio::OutputSubfield::extractSubfield(const pylith::topology::Field& f err = PetscSectionGetField(field.getLocalSection(), subfieldIndex, &subfieldSection);PYLITH_CHECK_ERROR(err); err = PetscSectionGetStorageSize(subfieldSection, &storageSize);PYLITH_CHECK_ERROR(err); - PetscVec subfieldVector = this->getVector(); + PetscVec subfieldVector = this->getOutputVector(); PetscInt subfieldSize = 0; err = VecGetLocalSize(subfieldVector, &subfieldSize);PYLITH_CHECK_ERROR(err); assert(subfieldSize == storageSize); diff --git a/libsrc/pylith/meshio/OutputSubfield.hh b/libsrc/pylith/meshio/OutputSubfield.hh index fe6c012f7f..9501a87832 100644 --- a/libsrc/pylith/meshio/OutputSubfield.hh +++ b/libsrc/pylith/meshio/OutputSubfield.hh @@ -15,7 +15,7 @@ #include "pylith/topology/FieldBase.hh" // HASA Description, Discretization -#include "pylith/topology/topologyfwd.hh" // USES Field +#include "pylith/topology/topologyfwd.hh" // USES Field, RefineInterpolator #include "pylith/utils/petscfwd.h" // HASA PetscVec class pylith::meshio::OutputSubfield : public pylith::utils::GenericComponent { @@ -30,12 +30,14 @@ public: * @param[in] mesh Mesh for subfield. * @param[in] name Name of subfield that will be extracted. * @param[in] basisOrder Basis order for subfield. + * @param[in] refineLevels Number of levels of mesh refinement. */ static OutputSubfield* create(const pylith::topology::Field& field, const pylith::topology::Mesh& mesh, const char* name, - const int basisOrder); + const int basisOrder, + const int refineLevels); /** Create OutputSubfield from Field. * @@ -76,17 +78,17 @@ public: */ int getBasisOrder(void) const; - /** Get PETSc global vector for projected subfield. + /** Get PETSc global vector for output subfield. * * @returns PETSc global vector. */ - PetscVec getVector(void) const; + PetscVec getOutputVector(void) const; - /** Get PETSc DM for projected subfield. + /** Get PETSc DM for output subfield. * * @returns PETSc DM. */ - PetscDM getDM(void) const; + PetscDM getOutputDM(void) const; /** Project PETSc vector to subfield. * @@ -121,10 +123,14 @@ protected: pylith::topology::FieldBase::Description _description; ///< Description of subfield. pylith::topology::FieldBase::Discretization _discretization; ///< Discretization of subfield. - PetscDM _dm; ///< PETSc DM for subfield. - PetscVec _vector; ///< PETSc global vector for subfield. - PetscPointFunc _fn; ///< PETSc point function for projection. PetscInt _subfieldIndex; ///< Index of subfield in fields. + PetscDM _projectDM; ///< PETSc global vector for subfield projection. + PetscVec _projectVector; ///< PETSc global vector for subfield projection. + PetscVec _projectVectorInterp; ///< PETSc global vector for subfield projection transformed for interpolation. + PetscPointFunc _fn; ///< PETSc point function for projection. + PetscDM _outputDM; ///< PETSc DM for subfield output. + PetscVec _outputVector; ///< PETSc global vector for subfield output. + pylith::topology::RefineInterpolator* _interpolator; ///< Interpolator for refined output. PetscDMLabel _label; ///< PETSc label associated with subfield. PetscInt _labelValue; ///< Value of PETSc label associated with subfield. diff --git a/libsrc/pylith/topology/Distributor.cc b/libsrc/pylith/topology/Distributor.cc index bf3f6e933d..1a67437848 100644 --- a/libsrc/pylith/topology/Distributor.cc +++ b/libsrc/pylith/topology/Distributor.cc @@ -109,7 +109,7 @@ pylith::topology::Distributor::distribute(pylith::topology::Mesh* const newMesh, err = DMPlexDistributeSetDefault(dmNew, PETSC_FALSE);PYLITH_CHECK_ERROR(err); err = DMPlexReorderCohesiveSupports(dmNew);PYLITH_CHECK_ERROR(err); err = DMViewFromOptions(dmNew, NULL, "-pylith_dist_dm_view");PYLITH_CHECK_ERROR(err); - newMesh->setDM(dmNew); + newMesh->setDM(dmNew, "domain"); PYLITH_METHOD_END; } // distribute @@ -171,8 +171,9 @@ pylith::topology::Distributor::write(meshio::DataWriter* const writer, partitionField.scatterLocalToOutput(); const int basisOrder = 0; + const int refineLevels = 0; pylith::meshio::OutputSubfield* outputField = - pylith::meshio::OutputSubfield::create(partitionField, mesh, "partition", basisOrder); + pylith::meshio::OutputSubfield::create(partitionField, mesh, "partition", basisOrder, refineLevels); outputField->project(partitionField.getOutputVector()); const PylithScalar t = 0.0; diff --git a/libsrc/pylith/topology/FieldOps.cc b/libsrc/pylith/topology/FieldOps.cc index 1ce3269fab..7f13738272 100644 --- a/libsrc/pylith/topology/FieldOps.cc +++ b/libsrc/pylith/topology/FieldOps.cc @@ -318,4 +318,100 @@ pylith::topology::FieldOps::createOutputLabel(const pylith::topology::Field* fie } +// ------------------------------------------------------------------------------------------------ +void +pylith::topology::FieldOps::transformVector(PetscVec* outputVector, + const PetscDM& outputDM, + const PetscVec& inputVector, + const PetscDM& inputDM) { + PYLITH_METHOD_BEGIN; + + assert(outputVector); + assert(outputDM); + assert(inputVector); + assert(inputDM); + + PetscErrorCode err = PETSC_SUCCESS; + PetscSection inputSection = PETSC_NULLPTR, outputSection = PETSC_NULLPTR; + err = DMGetGlobalSection(inputDM, &inputSection);PYLITH_CHECK_ERROR(err); + err = DMGetGlobalSection(outputDM, &outputSection);PYLITH_CHECK_ERROR(err); + + // Verify sizes + PetscInt outputNumFields = 0, inputNumFields = 0; + err = PetscSectionGetNumFields(inputSection, &inputNumFields);PYLITH_CHECK_ERROR(err); + err = PetscSectionGetNumFields(outputSection, &outputNumFields);PYLITH_CHECK_ERROR(err); + assert(inputNumFields == outputNumFields); + + PetscInt inputSize = 0, outputSize = 0; + err = PetscSectionGetStorageSize(inputSection, &inputSize);PYLITH_CHECK_ERROR(err); + err = PetscSectionGetStorageSize(outputSection, &outputSize);PYLITH_CHECK_ERROR(err); + assert(inputSize == outputSize); + + // Copy values from input vector to output vector + const PetscScalar* inputArray = PETSC_NULLPTR; + PetscScalar* outputArray = PETSC_NULLPTR; + err = VecGetArrayRead(inputVector, &inputArray);PYLITH_CHECK_ERROR(err); + err = VecGetArray(*outputVector, &outputArray);PYLITH_CHECK_ERROR(err); + + // Check whether inputDM and outputDM have the same names for points. + // We will assume that if inputDM and outputDM have the same range of points, + // then the names are the same and we ignore the subpoint map. + PetscDMLabel subpointMap = PETSC_NULLPTR; // Mapping of points in output DM back to input DM + PetscInt pStartIn, pEndIn, pStartOut, pEndOut; + PetscBool renamePoints = PETSC_FALSE; + err = DMPlexGetSubpointMap(outputDM, &subpointMap);PYLITH_CHECK_ERROR(err); + err = DMPlexGetChart(inputDM, &pStartIn, &pEndIn); + err = DMPlexGetChart(outputDM, &pStartOut, &pEndOut); + renamePoints = subpointMap && ((pStartIn != pStartOut) || (pEndIn != pEndOut)) ? PETSC_TRUE : PETSC_FALSE; + if (renamePoints) { + PetscIS subpointIS = PETSC_NULLPTR; // Mapping of points in output DM back to input DM + PetscInt subpointISSize = 0; + const PetscInt* subpointISPoints = PETSC_NULLPTR; + err = DMPlexGetSubpointIS(outputDM, &subpointIS);PYLITH_CHECK_ERROR(err); + err = ISGetSize(subpointIS, &subpointISSize);PYLITH_CHECK_ERROR(err); + err = ISGetIndices(subpointIS, &subpointISPoints);PYLITH_CHECK_ERROR(err); + + PetscInt pStart = 0, pEnd = 0; + err = PetscSectionGetChart(outputSection, &pStart, &pEnd); + for (PetscInt iPoint = 0; iPoint < subpointISSize; ++iPoint) { + PetscInt inputDof = 0, inputOffset = 0; + PetscInt outputDof = 0, outputOffset = 0; + const PetscInt inputPoint = subpointISPoints[iPoint]; + const PetscInt outputPoint = pStart + iPoint; + err = PetscSectionGetDof(inputSection, inputPoint, &inputDof);PYLITH_CHECK_ERROR(err); + err = PetscSectionGetOffset(inputSection, inputPoint, &inputOffset);PYLITH_CHECK_ERROR(err); + err = PetscSectionGetDof(outputSection, outputPoint, &outputDof);PYLITH_CHECK_ERROR(err); + err = PetscSectionGetOffset(outputSection, outputPoint, &outputOffset);PYLITH_CHECK_ERROR(err); + + assert(inputDof == outputDof); + for (PetscInt iDof = 0; iDof < inputDof; ++iDof) { + outputArray[outputOffset + iDof] = inputArray[inputOffset + iDof]; + } // for + } // for + } else { + PetscInt pStart = 0, pEnd = 0; + err = PetscSectionGetChart(outputSection, &pStart, &pEnd); + for (PetscInt point = pStart; point < pEnd; ++point) { + PetscInt inputDof = 0, inputOffset = 0; + PetscInt outputDof = 0, outputOffset = 0; + const PetscInt inputPoint = point; + const PetscInt outputPoint = point; + err = PetscSectionGetDof(inputSection, inputPoint, &inputDof);PYLITH_CHECK_ERROR(err); + err = PetscSectionGetOffset(inputSection, inputPoint, &inputOffset);PYLITH_CHECK_ERROR(err); + err = PetscSectionGetDof(outputSection, outputPoint, &outputDof);PYLITH_CHECK_ERROR(err); + err = PetscSectionGetOffset(outputSection, outputPoint, &outputOffset);PYLITH_CHECK_ERROR(err); + + assert(inputDof == outputDof); + for (PetscInt iDof = 0; iDof < inputDof; ++iDof) { + outputArray[outputOffset + iDof] = inputArray[inputOffset + iDof]; + } // for + } // for + } // if/else + err = VecRestoreArrayRead(inputVector, &inputArray);PYLITH_CHECK_ERROR(err); + err = VecRestoreArray(*outputVector, &outputArray);PYLITH_CHECK_ERROR(err); + + PYLITH_METHOD_END; +} + + // End of file diff --git a/libsrc/pylith/topology/FieldOps.hh b/libsrc/pylith/topology/FieldOps.hh index b434d08ba3..32307d088f 100644 --- a/libsrc/pylith/topology/FieldOps.hh +++ b/libsrc/pylith/topology/FieldOps.hh @@ -120,6 +120,21 @@ public: static void createOutputLabel(const pylith::topology::Field* field); + /** Transform vector using mapping from one PETSc DM to another. + * + * @note PETSc DMs must be compatible, including discretization. + * + * @param[out] outputVector Output vector. + * @param[in] outputDM PETSc DM associated with output vector. + * @param[out] outputVector Output vector. + * @param[in] outputDM PETSc DM associated with output vector. + */ + static + void transformVector(PetscVec* outputVector, + const PetscDM& outputDM, + const PetscVec& inputVector, + const PetscDM& inputDM); + /** Free saved PetscFE objects. */ static diff --git a/libsrc/pylith/topology/Makefile.am b/libsrc/pylith/topology/Makefile.am index 0143890fdb..566e69c864 100644 --- a/libsrc/pylith/topology/Makefile.am +++ b/libsrc/pylith/topology/Makefile.am @@ -31,6 +31,7 @@ subpkginclude_HEADERS = \ VisitorSubmesh.hh \ VisitorSubmesh.icc \ RefineUniform.hh \ + RefineInterpolator.hh \ topologyfwd.hh diff --git a/libsrc/pylith/topology/Mesh.cc b/libsrc/pylith/topology/Mesh.cc index cd64a6e38b..d417be8a3a 100644 --- a/libsrc/pylith/topology/Mesh.cc +++ b/libsrc/pylith/topology/Mesh.cc @@ -108,10 +108,12 @@ pylith::topology::Mesh::setDM(PetscDM dm, const char* label) { PYLITH_METHOD_BEGIN; - PetscErrorCode err; + PetscErrorCode err = PETSC_SUCCESS; err = DMDestroy(&_dm);PYLITH_CHECK_ERROR(err); _dm = dm; - err = PetscObjectSetName((PetscObject) _dm, label);PYLITH_CHECK_ERROR(err); + if (label) { + err = PetscObjectSetName((PetscObject) _dm, label);PYLITH_CHECK_ERROR(err); + } // if PYLITH_METHOD_END; } diff --git a/libsrc/pylith/topology/Mesh.hh b/libsrc/pylith/topology/Mesh.hh index e08cb52b37..006925d309 100644 --- a/libsrc/pylith/topology/Mesh.hh +++ b/libsrc/pylith/topology/Mesh.hh @@ -66,7 +66,7 @@ public: * @param label Label for mesh. */ void setDM(PetscDM dm, - const char* label="domain"); + const char* label=NULL); /** Set coordinate system. * diff --git a/libsrc/pylith/topology/MeshOps.cc b/libsrc/pylith/topology/MeshOps.cc index 3da68afd20..2ad9fff789 100644 --- a/libsrc/pylith/topology/MeshOps.cc +++ b/libsrc/pylith/topology/MeshOps.cc @@ -113,7 +113,7 @@ pylith::topology::MeshOps::createSubdomainMesh(const pylith::topology::Mesh& mes throw std::runtime_error(msg.str()); } // if - /* TODO: Add creation of pointSF for submesh */ + /* :TODO: Add creation of pointSF for submesh */ PetscDMLabel dmLabel = NULL; err = DMGetLabel(dmDomain, labelName, &dmLabel);PYLITH_CHECK_ERROR(err);assert(dmLabel); PetscBool hasLabelValue = PETSC_FALSE; @@ -146,9 +146,6 @@ pylith::topology::MeshOps::createSubdomainMesh(const pylith::topology::Mesh& mes throw std::runtime_error(msg.str()); } // if - // Set name - err = PetscObjectSetName((PetscObject) dmSubdomain, descriptiveLabel);PYLITH_CHECK_ERROR(err); - // Set lengthscale PylithScalar lengthScale; err = DMPlexGetScale(dmDomain, PETSC_UNIT_LENGTH, &lengthScale);PYLITH_CHECK_ERROR(err); @@ -156,7 +153,7 @@ pylith::topology::MeshOps::createSubdomainMesh(const pylith::topology::Mesh& mes pylith::topology::Mesh* submesh = new pylith::topology::Mesh(true);assert(submesh); submesh->setCoordSys(mesh.getCoordSys()); - submesh->setDM(dmSubdomain); + submesh->setDM(dmSubdomain, descriptiveLabel); _MeshOps::Events::logger.eventEnd(_MeshOps::Events::createSubdomainMesh); PYLITH_METHOD_RETURN(submesh); @@ -223,18 +220,15 @@ pylith::topology::MeshOps::createLowerDimMesh(const pylith::topology::Mesh& mesh throw std::runtime_error(msg.str()); } // if - // Set name - std::string meshLabel = "subdomain_" + std::string(labelName); - err = PetscObjectSetName((PetscObject) dmSubmesh, meshLabel.c_str());PYLITH_CHECK_ERROR(err); - // Set lengthscale PylithScalar lengthScale; err = DMPlexGetScale(dmDomain, PETSC_UNIT_LENGTH, &lengthScale);PYLITH_CHECK_ERROR(err); err = DMPlexSetScale(dmSubmesh, PETSC_UNIT_LENGTH, lengthScale);PYLITH_CHECK_ERROR(err); - pylith::topology::Mesh* submesh = new pylith::topology::Mesh(true);assert(submesh); submesh->setCoordSys(mesh.getCoordSys()); - submesh->setDM(dmSubmesh); + + std::string meshLabel = "subdomain_" + std::string(labelName); + submesh->setDM(dmSubmesh, meshLabel.c_str()); // Check topology MeshOps::checkTopology(*submesh); @@ -349,6 +343,55 @@ pylith::topology::MeshOps::nondimensionalize(Mesh* const mesh, } // nondimensionalize +// --------------------------------------------------------------------------------------------------------------------- +// Strip out "ghost" cells hanging off mesh +PetscDM +pylith::topology::MeshOps::removeHangingCells(const PetscDM& dmMesh) { + PYLITH_METHOD_BEGIN; + + PetscErrorCode err = PETSC_SUCCESS; + PetscDM dmClean = PETSC_NULLPTR; + + MPI_Comm comm = PetscObjectComm((PetscObject) dmMesh); + pylith::topology::Stratum cells(dmMesh, pylith::topology::Stratum::HEIGHT, 0); + DMPolytopeType cellType; + err = DMPlexGetCellType(dmMesh, cells.begin(), &cellType);PYLITH_CHECK_ERROR(err); + if (DMPolytopeTypeGetDim(cellType) < 0) { + // Hanging cells have dim == -1 + + // Create label over cells 1 dimension lower + PetscDMLabel labelInclude = PETSC_NULLPTR; + const PetscInt labelValue = 1; + err = DMLabelCreate(comm, "no_hanging_cells", &labelInclude);PYLITH_CHECK_ERROR(err); + pylith::topology::Stratum faces(dmMesh, pylith::topology::Stratum::HEIGHT, 1); + for (PetscInt face = faces.begin(); face < faces.end(); ++face) { + err = DMLabelSetValue(labelInclude, face, labelValue);PYLITH_CHECK_ERROR(err); + } // for + + err = DMPlexFilter(dmMesh, labelInclude, labelValue, PETSC_FALSE, PETSC_FALSE, PETSC_NULLPTR, &dmClean);PYLITH_CHECK_ERROR(err); + err = DMLabelDestroy(&labelInclude);PYLITH_CHECK_ERROR(err); + + // Create section using subpoint map to ensure sections are consistent. + PetscIS subpointIS = PETSC_NULLPTR; + PetscSection sectionOld = PETSC_NULLPTR, sectionNew = PETSC_NULLPTR; + err = DMPlexGetSubpointIS(dmClean, &subpointIS);PYLITH_CHECK_ERROR(err); + err = DMGetLocalSection(dmMesh, §ionOld);PYLITH_CHECK_ERROR(err); + err = PetscSectionCreateSubmeshSection(sectionOld, subpointIS, §ionNew);PYLITH_CHECK_ERROR(err); + err = DMSetLocalSection(dmClean, sectionNew);PYLITH_CHECK_ERROR(err); + err = PetscSectionDestroy(§ionNew);PYLITH_CHECK_ERROR(err); + + PetscReal lengthScale = 0.0; + err = DMPlexGetScale(dmMesh, PETSC_UNIT_LENGTH, &lengthScale);PYLITH_CHECK_ERROR(err); + err = DMPlexSetScale(dmClean, PETSC_UNIT_LENGTH, lengthScale);PYLITH_CHECK_ERROR(err); + } else { + dmClean = dmMesh; + err = PetscObjectReference((PetscObject) dmClean); + } // if/else + + PYLITH_METHOD_RETURN(dmClean); +} + + // --------------------------------------------------------------------------------------------------------------------- // Check topology of mesh. void diff --git a/libsrc/pylith/topology/MeshOps.hh b/libsrc/pylith/topology/MeshOps.hh index dbe17a507c..b04e5f736c 100644 --- a/libsrc/pylith/topology/MeshOps.hh +++ b/libsrc/pylith/topology/MeshOps.hh @@ -66,6 +66,15 @@ public: const PylithReal lengthScale, MPI_Comm comm); + /** Remove cells hanging off mesh. + * + * @param[in] dmMesh PETSc DM to clean. + * + * @returns PETSc DM without hanging cells. + */ + static + PetscDM removeHangingCells(const PetscDM& dmMesh); + /** Nondimensionalize the finite-element mesh. * * @param[in] mesh Finite-element mesh. diff --git a/libsrc/pylith/topology/RefineInterpolator.cc b/libsrc/pylith/topology/RefineInterpolator.cc new file mode 100644 index 0000000000..1ff604f27f --- /dev/null +++ b/libsrc/pylith/topology/RefineInterpolator.cc @@ -0,0 +1,191 @@ +// ================================================================================================= +// This code is part of PyLith, developed through the Computational Infrastructure +// for Geodynamics (https://github.com/geodynamics/pylith). +// +// Copyright (c) 2010-2024, University of California, Davis and the PyLith Development Team. +// All rights reserved. +// +// See https://mit-license.org/ and LICENSE.md and for license information. +// ================================================================================================= + +#include + +#include "pylith/topology/RefineInterpolator.hh" // implementation of class methods + +#include "pylith/topology/Mesh.hh" // USES Mesh +#include "pylith/topology/MeshOps.hh" // USES MeshOps +#include "pylith/topology/Field.hh" // USES Field +#include "pylith/topology/FieldOps.hh" // USES FieldOps + +#include "pylith/utils/error.hh" // USES PYLITH_METHOD_* +#include "pylith/utils/EventLogger.hh" // USES EventLogger + +#include // USES assert() + +// ------------------------------------------------------------------------------------------------ +namespace pylith { + namespace topology { + class _RefineInterpolator { +public: + + class Events { +public: + + static + void init(void); + + static pylith::utils::EventLogger logger; + static PylithInt initialize; + static PylithInt interpolate; + }; + + }; // _RefineInterpolator + } // topology +} // pylith + +pylith::utils::EventLogger pylith::topology::_RefineInterpolator::Events::logger; +PylithInt pylith::topology::_RefineInterpolator::Events::initialize; +PylithInt pylith::topology::_RefineInterpolator::Events::interpolate; + +// ------------------------------------------------------------------------------------------------ +void +pylith::topology::_RefineInterpolator::Events::init(void) { + logger.setClassName("RefineInterpolator"); + logger.initialize(); + initialize = logger.registerEvent("PL:RefineInterpolator:initialize"); + interpolate = logger.registerEvent("PL:RefineInterpolator:interpolate"); +} + + +// ------------------------------------------------------------------------------------------------ +// Constructor +pylith::topology::RefineInterpolator::RefineInterpolator(void) { + _RefineInterpolator::Events::init(); +} + + +// ------------------------------------------------------------------------------------------------ +// Destructor +pylith::topology::RefineInterpolator::~RefineInterpolator(void) { + deallocate(); +} + + +// ------------------------------------------------------------------------------------------------ +// Deallocate data structures. +void +pylith::topology::RefineInterpolator::deallocate(void) { + for (auto level : _levels) { + DMDestroy(&level.dm); + MatDestroy(&level.interpolateMatrix); + VecDestroy(&level.vector); + } // for + _levels.resize(0); +} // deallocate + + +// ------------------------------------------------------------------------------------------------ +// Get PETSc DM for input (coarsest level) +PetscDM +pylith::topology::RefineInterpolator::getInputDM(void) { + PetscDM dmStart = PETSC_NULLPTR; + if (_levels.size() > 0) { + PetscErrorCode err = PETSC_SUCCESS; + err = DMGetCoarseDM(_levels[0].dm, &dmStart);PYLITH_CHECK_ERROR(err); + } // if + return dmStart; +} + + +// ------------------------------------------------------------------------------------------------ +// Get PETSc DM for output (finest level) +PetscDM +pylith::topology::RefineInterpolator::getOutputDM(void) { + return (_levels.size() > 0) ? _levels[_levels.size()-1].dm : PETSC_NULLPTR; +} + + +// ------------------------------------------------------------------------------------------------ +// Initialize interpolation to refined mesh. +void +pylith::topology::RefineInterpolator::initialize(const PetscDM& dmMesh, + const int refineLevels, + const int outputBasisOrder, + const pylith::topology::FieldBase::Description& description, + const pylith::topology::FieldBase::Discretization& discretization) { + PYLITH_METHOD_BEGIN; + _RefineInterpolator::Events::logger.eventBegin(_RefineInterpolator::Events::initialize); + + _levels.resize(refineLevels); + PetscErrorCode err = PETSC_SUCCESS; + + PetscDM dmStart = pylith::topology::MeshOps::removeHangingCells(dmMesh);assert(dmStart); + err = DMCopyDisc(dmMesh, dmStart);PYLITH_CHECK_ERROR(err); + + PetscDM dmPrev = dmStart; + PetscReal lengthScale = 1.0; + MPI_Comm comm = PetscObjectComm((PetscObject) dmMesh); + for (size_t iLevel = 0; iLevel < _levels.size(); ++iLevel) { + _levels[iLevel].dm = PETSC_NULLPTR; + _levels[iLevel].interpolateMatrix = PETSC_NULLPTR; + _levels[iLevel].vector = PETSC_NULLPTR; + + err = DMPlexSetRefinementUniform(dmPrev, PETSC_TRUE);PYLITH_CHECK_ERROR(err); + err = DMRefine(dmPrev, comm, &_levels[iLevel].dm);PYLITH_CHECK_ERROR(err); + err = DMSetCoarseDM(_levels[iLevel].dm, dmPrev);PYLITH_CHECK_ERROR(err); + err = DMPlexGetScale(dmPrev, PETSC_UNIT_LENGTH, &lengthScale);PYLITH_CHECK_ERROR(err); + err = DMPlexSetScale(_levels[iLevel].dm, PETSC_UNIT_LENGTH, lengthScale);PYLITH_CHECK_ERROR(err); + err = DMPlexReorderSetDefault(_levels[iLevel].dm, DM_REORDER_DEFAULT_FALSE); + +#if 0 // needed for higher order coordinates (not needed for affine coordinates) + PetscCall(DMPlexCreateCoordinateSpace(rdm, rd, PETSC_FALSE, NULL)); + PetscCall(PetscObjectSetName((PetscObject)rdm, "Refined Mesh with Linear Coordinates")); + PetscCall(DMGetCoordinateDM(odm, &cdm)); + PetscCall(DMGetCoordinateDM(rdm, &rcdm)); + PetscCall(DMGetCoordinatesLocal(odm, &cl)); + PetscCall(DMGetCoordinatesLocal(rdm, &rcl)); +#endif + + if (iLevel < _levels.size()-1) { + err = DMCopyDisc(dmPrev, _levels[iLevel].dm);PYLITH_CHECK_ERROR(err); + } else { + const PetscInt minBasisOrder = PETSC_DETERMINE; + const PetscInt maxBasisOrder = outputBasisOrder; + err = DMCopyFields(dmPrev, minBasisOrder, maxBasisOrder, _levels[iLevel].dm);PYLITH_CHECK_ERROR(err); + err = DMCopyDS(dmPrev, minBasisOrder, maxBasisOrder, _levels[iLevel].dm);PYLITH_CHECK_ERROR(err); + } // else + err = DMCreateGlobalVector(_levels[iLevel].dm, &_levels[iLevel].vector);PYLITH_CHECK_ERROR(err); + err = DMCreateInterpolation(dmPrev, _levels[iLevel].dm, &_levels[iLevel].interpolateMatrix, NULL);PYLITH_CHECK_ERROR(err); + + dmPrev = _levels[iLevel].dm; + } // for + err = DMDestroy(&dmStart);PYLITH_CHECK_ERROR(err); + + _RefineInterpolator::Events::logger.eventEnd(_RefineInterpolator::Events::initialize); + PYLITH_METHOD_END; +} // initialize + + +// ------------------------------------------------------------------------------------------------ +// Interpolate field to fine mesh level. +void +pylith::topology::RefineInterpolator::interpolate(const PetscVec* vectorOut, + const PetscVec& vectorIn) { + PYLITH_METHOD_BEGIN; + _RefineInterpolator::Events::logger.eventBegin(_RefineInterpolator::Events::interpolate); + assert(vectorOut); + + PetscVec vectorPrev = vectorIn; + PetscErrorCode err = PETSC_SUCCESS; + for (auto level : _levels) { + err = MatMult(level.interpolateMatrix, vectorPrev, level.vector);PYLITH_CHECK_ERROR(err); + vectorPrev = level.vector; + } // for + err = VecCopy(vectorPrev, *vectorOut); + + _RefineInterpolator::Events::logger.eventEnd(_RefineInterpolator::Events::interpolate); + PYLITH_METHOD_END; +} // interpolate + + +// End of file diff --git a/libsrc/pylith/topology/RefineInterpolator.hh b/libsrc/pylith/topology/RefineInterpolator.hh new file mode 100644 index 0000000000..7614560a89 --- /dev/null +++ b/libsrc/pylith/topology/RefineInterpolator.hh @@ -0,0 +1,84 @@ +// ================================================================================================= +// This code is part of PyLith, developed through the Computational Infrastructure +// for Geodynamics (https://github.com/geodynamics/pylith). +// +// Copyright (c) 2010-2024, University of California, Davis and the PyLith Development Team. +// All rights reserved. +// +// See https://mit-license.org/ and LICENSE.md and for license information. +// ================================================================================================= +#pragma once + +#include "pylith/topology/topologyfwd.hh" // forward declarations +#include "pylith/utils/petscfwd.h" // USES PetscDM + +#include "pylith/topology/FieldBase.hh" // USES FieldBase + +#include // HASA std::vector + +/// @brief Interpolate fields to uniformly refined mesh. +class pylith::topology::RefineInterpolator { + friend class TestRefineUniform; // unit testing + + // PUBLIC METHODS ///////////////////////////////////////////////////////////////////////////// +public: + + /// Constructor + RefineInterpolator(void); + + /// Destructor + ~RefineInterpolator(void); + + /// Deallocate data structures. + void deallocate(void); + + /** Get PETSc DM for input (coarsest level) + * + */ + PetscDM getInputDM(void); + + /** Get PETSc DM for output (finest level) + * + */ + PetscDM getOutputDM(void); + + /** Initialize interpolation to refined mesh. + * + * @param[in] dmMesh PETSc DM for starting point of refinement. + * @param[in] refineLevels Number of levels of mesh refinement. + * @param[in] outputBasisOrder Basis order for output. + */ + void initialize(const PetscDM& dmMesh, + const int refineLevels, + const int outputBasisOrder, + const pylith::topology::FieldBase::Description& description, + const pylith::topology::FieldBase::Discretization& discretization); + + /** Interpolate field to fine mesh level. + * + * @param[out] vectorOut PETSc Vec with interpolated field. + * @param[in] vectorIn PETSc Vec with field to interpolate. + */ + void interpolate(const PetscVec* vectorOut, + const PetscVec& vectorIn); + + // PRIvATE MEMBERS //////////////////////////////////////////////////////////////////////////// +private: + + struct Level { + PetscDM dm; + PetscMat interpolateMatrix; + PetscVec vector; + }; + + std::vector _levels; ///< Information at each refinement level. + + // NOT IMPLEMENTED //////////////////////////////////////////////////////////////////////////// +private: + + RefineInterpolator(const RefineInterpolator&); ///< Not implemented + const RefineInterpolator& operator=(const RefineInterpolator&); ///< Not implemented + +}; // RefineInterpolator + +// End of file diff --git a/libsrc/pylith/topology/RefineUniform.cc b/libsrc/pylith/topology/RefineUniform.cc index c4830712d5..d0250b7b9c 100644 --- a/libsrc/pylith/topology/RefineUniform.cc +++ b/libsrc/pylith/topology/RefineUniform.cc @@ -123,7 +123,7 @@ pylith::topology::RefineUniform::refine(Mesh* const newMesh, } // for err = DMPlexReorderCohesiveSupports(dmNew);PYLITH_CHECK_ERROR(err); - newMesh->setDM(dmNew); + newMesh->setDM(dmNew, "domain"); _RefineUniform::Events::logger.eventBegin(_RefineUniform::Events::refineFixCellLabel); // Remove all non-cells from cells label diff --git a/libsrc/pylith/topology/ReverseCuthillMcKee.cc b/libsrc/pylith/topology/ReverseCuthillMcKee.cc index d4e8e9fadf..efadce8fb6 100644 --- a/libsrc/pylith/topology/ReverseCuthillMcKee.cc +++ b/libsrc/pylith/topology/ReverseCuthillMcKee.cc @@ -32,7 +32,7 @@ pylith::topology::ReverseCuthillMcKee::reorder(topology::Mesh* mesh) { err = DMPlexGetOrdering(dmOrig, MATORDERINGRCM, dmLabel, &permutation);PYLITH_CHECK_ERROR(err); err = DMPlexPermute(dmOrig, permutation, &dmNew);PYLITH_CHECK_ERROR(err); err = ISDestroy(&permutation);PYLITH_CHECK_ERROR(err); - mesh->setDM(dmNew); + mesh->setDM(dmNew, "domain"); // Verify that all material points (cells) are consecutive. PetscIS valuesIS = NULL; diff --git a/libsrc/pylith/topology/Stratum.icc b/libsrc/pylith/topology/Stratum.icc index 7cf76a5672..321a347b9d 100644 --- a/libsrc/pylith/topology/Stratum.icc +++ b/libsrc/pylith/topology/Stratum.icc @@ -16,7 +16,7 @@ inline pylith::topology::Stratum::Stratum(const PetscDM dmMesh, const StratumEnum stype, - const int level) { // constructor + const int level) { assert(dmMesh); PetscErrorCode err = 0; switch (stype) { @@ -37,7 +37,7 @@ pylith::topology::Stratum::Stratum(const PetscDM dmMesh, // ---------------------------------------------------------------------- // Default destructor inline -pylith::topology::Stratum::~Stratum(void) { // destructor +pylith::topology::Stratum::~Stratum(void) { _begin = 0; _end = 0; } // destructor @@ -47,7 +47,7 @@ pylith::topology::Stratum::~Stratum(void) { // destructor // Get starting point. inline PetscInt -pylith::topology::Stratum::begin(void) const { // begin +pylith::topology::Stratum::begin(void) const { return _begin; } // begin diff --git a/libsrc/pylith/topology/topologyfwd.hh b/libsrc/pylith/topology/topologyfwd.hh index d0eac0df64..6d6b2e18c0 100644 --- a/libsrc/pylith/topology/topologyfwd.hh +++ b/libsrc/pylith/topology/topologyfwd.hh @@ -29,6 +29,7 @@ namespace pylith { class FieldFactory; class FieldOps; class FieldQuery; + class RefineInterpolator; class MatVisitorMesh; class MatVisitorSubmesh; diff --git a/modulesrc/meshio/OutputObserver.i b/modulesrc/meshio/OutputObserver.i index 0d4483ef62..c3b65855da 100644 --- a/modulesrc/meshio/OutputObserver.i +++ b/modulesrc/meshio/OutputObserver.i @@ -54,6 +54,12 @@ public: */ void setOutputBasisOrder(const int value); + /** Set number of mesh refinement levels for output. + * + * @param[in] value Number of mesh refinement levels for output. + */ + void setRefineLevels(const int value); + /** Set time scale. * * @param[in] value Time scale for dimensionalizing time. diff --git a/modulesrc/topology/Mesh.i b/modulesrc/topology/Mesh.i index 01828c443d..640982036c 100644 --- a/modulesrc/topology/Mesh.i +++ b/modulesrc/topology/Mesh.i @@ -45,20 +45,6 @@ public: */ Mesh* clone(void) const; - /** Get DMPlex mesh. - * - * @returns DMPlex mesh. - */ - PetscDM getDM(void) const; - - /** Set DMPlex mesh. - * - * @param DMPlex mesh. - * @param label Label for mesh. - */ - void setDM(PetscDM dm, - const char* label="domain"); - /** Set coordinate system. * * @param cs Coordinate system. diff --git a/pylith/meshio/OutputObserver.py b/pylith/meshio/OutputObserver.py index 1e9995569f..29769acec7 100644 --- a/pylith/meshio/OutputObserver.py +++ b/pylith/meshio/OutputObserver.py @@ -30,6 +30,9 @@ class OutputObserver(PetscComponent, ModuleOutputObserver): outputBasisOrder = pythia.pyre.inventory.int("output_basis_order", default=1, validator=pythia.pyre.inventory.choice([0,1])) outputBasisOrder.meta['tip'] = "Basis order for output." + refineLevels = pythia.pyre.inventory.int("refine_levels", default=0, validator=pythia.pyre.inventory.greaterEqual(0)) + refineLevels.meta['tip'] = "Number of mesh refinement levels for output." + def __init__(self, name="outputobserver"): """Constructor. """ @@ -50,6 +53,7 @@ def preinitialize(self, problem): else: outputBasisOrder = self.outputBasisOrder ModuleOutputObserver.setOutputBasisOrder(self, outputBasisOrder) + ModuleOutputObserver.setRefineLevels(self, self.refineLevels) self.writer.preinitialize() ModuleOutputObserver.setWriter(self, self.writer) diff --git a/pylith/meshio/OutputSolnDomain.py b/pylith/meshio/OutputSolnDomain.py index 184bf8e265..8e92c88c9a 100644 --- a/pylith/meshio/OutputSolnDomain.py +++ b/pylith/meshio/OutputSolnDomain.py @@ -35,7 +35,10 @@ class OutputSolnDomain(OutputSoln, ModuleOutputSolnDomain): writer = pylith.meshio.DataWriterHDF5 writer.filename = domain.h5 + # Output with a basis order of 1 and refine mesh 3x (cells are 1/8 size). + # Refining the output mesh is useful with a basis order of 2 or greater in the solution. output_basis_order = 1 + refine_levels = 3 """ } diff --git a/pylith/meshio/gmsh_utils.py b/pylith/meshio/gmsh_utils.py index 98f792f1d5..9ded3b6a06 100644 --- a/pylith/meshio/gmsh_utils.py +++ b/pylith/meshio/gmsh_utils.py @@ -101,7 +101,7 @@ def main(self): """ args = self._parse_command_line() - self.initialize(args.name) + self.initialize(args) if args.geometry: self.create_geometry() if args.mark: @@ -112,13 +112,12 @@ def main(self): self.write(args.filename, args.binary) self.finalize(args.gui) - def initialize(self, name: str): + def initialize(self, args): """Initialize Gmsh. :param name: Name for mesh. """ gmsh.initialize() - gmsh.model.add(name) def finalize(self, gui=False): """Finalize Gmsh. @@ -174,7 +173,7 @@ def _parse_command_line(self): parser.add_argument("--gui", action="store_true", dest="gui", help="Show GUI after running steps.") - GenerateMesh._add_arguments(parser) + self._add_arguments(parser) args = parser.parse_args() if args.write: args.generate = True diff --git a/tests/fullscale/linearelasticity/faults-2d/Makefile.am b/tests/fullscale/linearelasticity/faults-2d/Makefile.am index c111ba7e10..f33c46e043 100644 --- a/tests/fullscale/linearelasticity/faults-2d/Makefile.am +++ b/tests/fullscale/linearelasticity/faults-2d/Makefile.am @@ -43,6 +43,9 @@ dist_noinst_DATA = \ shearnoslip.cfg \ shearnoslip_quad.cfg \ shearnoslip_tri.cfg \ + shearnoslip_refineoutput.cfg \ + shearnoslip_refineoutput_quad.cfg \ + shearnoslip_refineoutput_tri.cfg \ zeroslipfn.timedb diff --git a/tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py b/tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py index c5caa5de66..6cc3641fe4 100644 --- a/tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py +++ b/tests/fullscale/linearelasticity/faults-2d/TestShearNoSlip.py @@ -103,11 +103,37 @@ def setUp(self): return +# ------------------------------------------------------------------------------------------------- +class TestQuadGmshRefineOutput(TestCase): + + def setUp(self): + self.name = "shearnoslip_refineoutput_quad" + self.mesh = meshes.QuadGmshRefineOutput() + super().setUp() + + TestCase.run_pylith(self, self.name, ["shearnoslip.cfg", "shearnoslip_refineoutput.cfg", "shearnoslip_refineoutput_quad.cfg"]) + return + + +# ------------------------------------------------------------------------------------------------- +class TestTriGmshRefineOutput(TestCase): + + def setUp(self): + self.name = "shearnoslip_refineoutput_tri" + self.mesh = meshes.TriGmshRefineOutput() + super().setUp() + + TestCase.run_pylith(self, self.name, ["shearnoslip.cfg", "shearnoslip_refineoutput.cfg", "shearnoslip_refineoutput_tri.cfg"]) + return + + # ------------------------------------------------------------------------------------------------- def test_cases(): return [ TestQuadGmsh, TestTriGmsh, + TestQuadGmshRefineOutput, + TestTriGmshRefineOutput, ] diff --git a/tests/fullscale/linearelasticity/faults-2d/meshes.py b/tests/fullscale/linearelasticity/faults-2d/meshes.py index cd990f37cf..c6dd00d6f8 100644 --- a/tests/fullscale/linearelasticity/faults-2d/meshes.py +++ b/tests/fullscale/linearelasticity/faults-2d/meshes.py @@ -63,6 +63,58 @@ class QuadGmsh(object): } +class TriGmshRefineOutput(object): + """Mesh information for tri mesh using Gmsh with 1 level of refinement for output. + """ + ENTITIES = { + "domain": MeshEntity(ncells=164*4, ncorners=3, nvertices=361+2*(9+8)), + "points": MeshEntity(ncells=3, ncorners=1, nvertices=3), + + # Materials + "mat_xneg": MeshEntity(ncells=38*4, ncorners=3, nvertices=97), + "mat_xmid": MeshEntity(ncells=38, ncorners=3, nvertices=30), + "mat_xposypos": MeshEntity(ncells=44, ncorners=3, nvertices=31), + "mat_xposyneg": MeshEntity(ncells=44, ncorners=3, nvertices=31), + + # Faults + "fault_xmid": MeshEntity(ncells=8*2, ncorners=2, nvertices=9+8), + "fault_xneg": MeshEntity(ncells=8, ncorners=2, nvertices=9), + + # Boundaries + "bc_xneg": MeshEntity(ncells=8, ncorners=2, nvertices=9), + "bc_xpos": MeshEntity(ncells=8*2, ncorners=2, nvertices=9+8), + "bc_yneg": MeshEntity(ncells=8, ncorners=2, nvertices=9+2), + "bc_ypos": MeshEntity(ncells=8*2, ncorners=2, nvertices=9+8+2), + "boundary_ypos": MeshEntity(ncells=16, ncorners=2, nvertices=17+2), + } + + +class QuadGmshRefineOutput(object): + """Mesh information for quad mesh using Gmsh with 1 level of refinement for output. + """ + ENTITIES = { + "domain": MeshEntity(ncells=100*4, ncorners=4, nvertices=(11+10)**2 + 2*(11+10)), + "points": MeshEntity(ncells=3, ncorners=1, nvertices=3), + + # Materials + "mat_xneg": MeshEntity(ncells=30*4, ncorners=3, nvertices=147), + "mat_xmid": MeshEntity(ncells=20, ncorners=3, nvertices=33), + "mat_xposypos": MeshEntity(ncells=25, ncorners=3, nvertices=36), + "mat_xposyneg": MeshEntity(ncells=25, ncorners=3, nvertices=36), + + # Faults + "fault_xmid": MeshEntity(ncells=20*2, ncorners=2, nvertices=11+10), + "fault_xneg": MeshEntity(ncells=10, ncorners=2, nvertices=11), + + # Boundaries + "bc_xneg": MeshEntity(ncells=10, ncorners=2, nvertices=11), + "bc_xpos": MeshEntity(ncells=20*2, ncorners=2, nvertices=11+10), + "bc_yneg": MeshEntity(ncells=11, ncorners=2, nvertices=11+2), + "bc_ypos": MeshEntity(ncells=11*2, ncorners=2, nvertices=11+10+2), + "boundary_ypos": MeshEntity(ncells=11*2, ncorners=2, nvertices=11+10+2), + } + + class TriCubit(object): """Mesh information for tri mesh using Cubit. """ diff --git a/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput.cfg b/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput.cfg new file mode 100644 index 0000000000..329d370471 --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput.cfg @@ -0,0 +1,31 @@ +# ---------------------------------------------------------------------- +# problem +# ---------------------------------------------------------------------- +[pylithapp.problem.solution_observers] +domain.refine_levels = 1 +boundary_ypos.refine_levels = 1 + + +# ---------------------------------------------------------------------- +# materials +# ---------------------------------------------------------------------- +[pylithapp.problem.materials] +mat_xneg.observers.observer.refine_levels = 1 + + +# ---------------------------------------------------------------------- +# faults +# ---------------------------------------------------------------------- +[pylithapp.problem.interfaces] +fault_xmid.observers.observer.refine_levels = 1 + + +# ---------------------------------------------------------------------- +# boundary conditions +# ---------------------------------------------------------------------- +[pylithapp.problem.bc] +bc_xpos.observers.observer.refine_levels = 1 +bc_ypos.observers.observer.refine_levels = 1 + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput_quad.cfg b/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput_quad.cfg new file mode 100644 index 0000000000..79aaa60f3d --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput_quad.cfg @@ -0,0 +1,16 @@ +[pylithapp.metadata] +base = [pylithapp.cfg, shearnoslip.cfg, shearnoslip_refineoutput.cfg] +keywords = [quadrilateral cells] +arguments = [shearnoslip_refineoutput.cfg, shearnoslip_refineoutput_quad.cfg] + +[pylithapp.problem] +defaults.name = shearnoslip_refineoutput_quad + +# ---------------------------------------------------------------------- +# mesh_generator +# ---------------------------------------------------------------------- +[pylithapp.mesh_generator.reader] +filename = mesh_quad.msh + + +# End of file diff --git a/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput_tri.cfg b/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput_tri.cfg new file mode 100644 index 0000000000..9652448ad5 --- /dev/null +++ b/tests/fullscale/linearelasticity/faults-2d/shearnoslip_refineoutput_tri.cfg @@ -0,0 +1,16 @@ +[pylithapp.metadata] +base = [pylithapp.cfg, shearnoslip.cfg, shearnoslip_refineoutput.cfg] +keywords = [triangular cells] +arguments = [shearnoslip_refineoutput.cfg, shearnoslip_refineoutput_tri.cfg] + +[pylithapp.problem] +defaults.name = shearnoslip_refineoutput_tri + +# ---------------------------------------------------------------------- +# mesh_generator +# ---------------------------------------------------------------------- +[pylithapp.mesh_generator.reader] +filename = mesh_tri.msh + + +# End of file diff --git a/tests/libtests/meshio/TestDataWriterHDF5ExtMaterial.cc b/tests/libtests/meshio/TestDataWriterHDF5ExtMaterial.cc index db0d2dc91f..48d265542b 100644 --- a/tests/libtests/meshio/TestDataWriterHDF5ExtMaterial.cc +++ b/tests/libtests/meshio/TestDataWriterHDF5ExtMaterial.cc @@ -86,7 +86,7 @@ pylith::meshio::TestDataWriterHDF5ExtMaterial::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_materialMesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_materialMesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -127,7 +127,7 @@ pylith::meshio::TestDataWriterHDF5ExtMaterial::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_materialMesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_materialMesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestDataWriterHDF5ExtMesh.cc b/tests/libtests/meshio/TestDataWriterHDF5ExtMesh.cc index c8bf8e15de..f4fb71b22f 100644 --- a/tests/libtests/meshio/TestDataWriterHDF5ExtMesh.cc +++ b/tests/libtests/meshio/TestDataWriterHDF5ExtMesh.cc @@ -156,7 +156,7 @@ pylith::meshio::TestDataWriterHDF5ExtMesh::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_mesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_mesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -197,7 +197,7 @@ pylith::meshio::TestDataWriterHDF5ExtMesh::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_mesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_mesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestDataWriterHDF5ExtSubmesh.cc b/tests/libtests/meshio/TestDataWriterHDF5ExtSubmesh.cc index ff03340964..c237cbec9a 100644 --- a/tests/libtests/meshio/TestDataWriterHDF5ExtSubmesh.cc +++ b/tests/libtests/meshio/TestDataWriterHDF5ExtSubmesh.cc @@ -87,7 +87,7 @@ pylith::meshio::TestDataWriterHDF5ExtSubmesh::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_submesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_submesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -128,7 +128,7 @@ pylith::meshio::TestDataWriterHDF5ExtSubmesh::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_submesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_submesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestDataWriterHDF5Material.cc b/tests/libtests/meshio/TestDataWriterHDF5Material.cc index 403f5ef896..7dd65469cc 100644 --- a/tests/libtests/meshio/TestDataWriterHDF5Material.cc +++ b/tests/libtests/meshio/TestDataWriterHDF5Material.cc @@ -86,7 +86,7 @@ pylith::meshio::TestDataWriterHDF5Material::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_materialMesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_materialMesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -127,7 +127,7 @@ pylith::meshio::TestDataWriterHDF5Material::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_materialMesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_materialMesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestDataWriterHDF5Mesh.cc b/tests/libtests/meshio/TestDataWriterHDF5Mesh.cc index f0be68e6cd..f3aa359122 100644 --- a/tests/libtests/meshio/TestDataWriterHDF5Mesh.cc +++ b/tests/libtests/meshio/TestDataWriterHDF5Mesh.cc @@ -129,7 +129,7 @@ pylith::meshio::TestDataWriterHDF5Mesh::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_mesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_mesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -170,7 +170,7 @@ pylith::meshio::TestDataWriterHDF5Mesh::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_mesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_mesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestDataWriterHDF5Submesh.cc b/tests/libtests/meshio/TestDataWriterHDF5Submesh.cc index 2ca0d280d6..d185943a96 100644 --- a/tests/libtests/meshio/TestDataWriterHDF5Submesh.cc +++ b/tests/libtests/meshio/TestDataWriterHDF5Submesh.cc @@ -87,7 +87,7 @@ pylith::meshio::TestDataWriterHDF5Submesh::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_submesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_submesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -128,7 +128,7 @@ pylith::meshio::TestDataWriterHDF5Submesh::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_submesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_submesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestDataWriterMaterial.cc b/tests/libtests/meshio/TestDataWriterMaterial.cc index 7b8d7bd06e..2bf614659b 100644 --- a/tests/libtests/meshio/TestDataWriterMaterial.cc +++ b/tests/libtests/meshio/TestDataWriterMaterial.cc @@ -270,7 +270,7 @@ pylith::meshio::TestDataWriterMaterial::_initialize(void) { } // if delete _materialMesh; - _materialMesh = pylith::topology::MeshOps::createSubdomainMesh(*_domainMesh, pylith::topology::Mesh::cells_label_name, data->materialId, ":UNKNOWN:"); + _materialMesh = pylith::topology::MeshOps::createSubdomainMesh(*_domainMesh, pylith::topology::Mesh::cells_label_name, data->materialId, "material"); assert(_materialMesh); PYLITH_METHOD_END; diff --git a/tests/libtests/meshio/TestDataWriterVTKMaterial.cc b/tests/libtests/meshio/TestDataWriterVTKMaterial.cc index 096a6a71bc..0addb28da3 100644 --- a/tests/libtests/meshio/TestDataWriterVTKMaterial.cc +++ b/tests/libtests/meshio/TestDataWriterVTKMaterial.cc @@ -104,7 +104,7 @@ pylith::meshio::TestDataWriterVTKMaterial::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_materialMesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_materialMesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -148,7 +148,7 @@ pylith::meshio::TestDataWriterVTKMaterial::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_materialMesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_materialMesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestDataWriterVTKMesh.cc b/tests/libtests/meshio/TestDataWriterVTKMesh.cc index 4604f55da0..ff1fe16fc8 100644 --- a/tests/libtests/meshio/TestDataWriterVTKMesh.cc +++ b/tests/libtests/meshio/TestDataWriterVTKMesh.cc @@ -163,7 +163,7 @@ pylith::meshio::TestDataWriterVTKMesh::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_mesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_mesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -206,7 +206,7 @@ pylith::meshio::TestDataWriterVTKMesh::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_mesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_mesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestDataWriterVTKSubmesh.cc b/tests/libtests/meshio/TestDataWriterVTKSubmesh.cc index 6e6e4601f4..5b3d8b7c98 100644 --- a/tests/libtests/meshio/TestDataWriterVTKSubmesh.cc +++ b/tests/libtests/meshio/TestDataWriterVTKSubmesh.cc @@ -101,7 +101,7 @@ pylith::meshio::TestDataWriterVTKSubmesh::testWriteVertexField(void) { const pylith::string_vector& subfieldNames = vertexField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(vertexField, *_submesh, subfieldNames[i].c_str(), 1); + OutputSubfield* subfield = OutputSubfield::create(vertexField, *_submesh, subfieldNames[i].c_str(), 1, 0); assert(subfield); subfield->project(vertexField.getOutputVector()); writer.writeVertexField(t, *subfield); @@ -146,7 +146,7 @@ pylith::meshio::TestDataWriterVTKSubmesh::testWriteCellField(void) { const pylith::string_vector& subfieldNames = cellField.getSubfieldNames(); const size_t numFields = subfieldNames.size(); for (size_t i = 0; i < numFields; ++i) { - OutputSubfield* subfield = OutputSubfield::create(cellField, *_submesh, subfieldNames[i].c_str(), 0); + OutputSubfield* subfield = OutputSubfield::create(cellField, *_submesh, subfieldNames[i].c_str(), 0, 0); assert(subfield); subfield->project(cellField.getOutputVector()); writer.writeCellField(t, *subfield); diff --git a/tests/libtests/meshio/TestMeshIO.cc b/tests/libtests/meshio/TestMeshIO.cc index 3a9dac285b..c9a44ca317 100644 --- a/tests/libtests/meshio/TestMeshIO.cc +++ b/tests/libtests/meshio/TestMeshIO.cc @@ -114,7 +114,7 @@ pylith::meshio::TestMeshIO::_createMesh(void) { } // for err = DMPlexCreateFromCellListPetsc(_mesh->getComm(), _data->cellDim, _data->numCells, _data->numVertices, _data->numCorners, interpolateMesh, cells, _data->spaceDim, _data->vertices, &dmMesh);PYLITH_CHECK_ERROR(err); delete [] cells; - _mesh->setDM(dmMesh); + _mesh->setDM(dmMesh, "domain"); // Material ids PylithInt cStart, cEnd; diff --git a/tests/libtests/meshio/data/hex8_mat_cell_t10.vtk b/tests/libtests/meshio/data/hex8_mat_cell_t10.vtk index 13db520ad4..fc526eb108 100644 --- a/tests/libtests/meshio/data/hex8_mat_cell_t10.vtk +++ b/tests/libtests/meshio/data/hex8_mat_cell_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +material ASCII DATASET UNSTRUCTURED_GRID POINTS 8 double diff --git a/tests/libtests/meshio/data/hex8_mat_vertex_t10.vtk b/tests/libtests/meshio/data/hex8_mat_vertex_t10.vtk index bba4d463ba..5874d5b6d3 100644 --- a/tests/libtests/meshio/data/hex8_mat_vertex_t10.vtk +++ b/tests/libtests/meshio/data/hex8_mat_vertex_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +material ASCII DATASET UNSTRUCTURED_GRID POINTS 8 double diff --git a/tests/libtests/meshio/data/hex8_surf_cell_t10.vtk b/tests/libtests/meshio/data/hex8_surf_cell_t10.vtk index 5141dd8d63..d0bd7cdd40 100644 --- a/tests/libtests/meshio/data/hex8_surf_cell_t10.vtk +++ b/tests/libtests/meshio/data/hex8_surf_cell_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +subdomain_top ASCII DATASET UNSTRUCTURED_GRID POINTS 6 double diff --git a/tests/libtests/meshio/data/hex8_surf_vertex_t10.vtk b/tests/libtests/meshio/data/hex8_surf_vertex_t10.vtk index ac0ddd8d61..56bcc135f2 100644 --- a/tests/libtests/meshio/data/hex8_surf_vertex_t10.vtk +++ b/tests/libtests/meshio/data/hex8_surf_vertex_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +subdomain_top ASCII DATASET UNSTRUCTURED_GRID POINTS 6 double diff --git a/tests/libtests/meshio/data/quad4_mat_cell_t10.vtk b/tests/libtests/meshio/data/quad4_mat_cell_t10.vtk index c94c13a22f..c37c7dff7b 100644 --- a/tests/libtests/meshio/data/quad4_mat_cell_t10.vtk +++ b/tests/libtests/meshio/data/quad4_mat_cell_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +material ASCII DATASET UNSTRUCTURED_GRID POINTS 4 double diff --git a/tests/libtests/meshio/data/quad4_mat_vertex_t10.vtk b/tests/libtests/meshio/data/quad4_mat_vertex_t10.vtk index e0324b73ca..c440a56d1d 100644 --- a/tests/libtests/meshio/data/quad4_mat_vertex_t10.vtk +++ b/tests/libtests/meshio/data/quad4_mat_vertex_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +material ASCII DATASET UNSTRUCTURED_GRID POINTS 4 double diff --git a/tests/libtests/meshio/data/quad4_surf_cell_t10.vtk b/tests/libtests/meshio/data/quad4_surf_cell_t10.vtk index 3da7d562ee..fda62a7282 100644 --- a/tests/libtests/meshio/data/quad4_surf_cell_t10.vtk +++ b/tests/libtests/meshio/data/quad4_surf_cell_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +subdomain_bc3 ASCII DATASET UNSTRUCTURED_GRID POINTS 3 double diff --git a/tests/libtests/meshio/data/quad4_surf_vertex_t10.vtk b/tests/libtests/meshio/data/quad4_surf_vertex_t10.vtk index a8c9e25451..0f2d7b5c76 100644 --- a/tests/libtests/meshio/data/quad4_surf_vertex_t10.vtk +++ b/tests/libtests/meshio/data/quad4_surf_vertex_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +subdomain_bc3 ASCII DATASET UNSTRUCTURED_GRID POINTS 3 double diff --git a/tests/libtests/meshio/data/tet4_mat_cell_t10.vtk b/tests/libtests/meshio/data/tet4_mat_cell_t10.vtk index 48ca2e35a1..0cf84701b4 100644 --- a/tests/libtests/meshio/data/tet4_mat_cell_t10.vtk +++ b/tests/libtests/meshio/data/tet4_mat_cell_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +material ASCII DATASET UNSTRUCTURED_GRID POINTS 8 double diff --git a/tests/libtests/meshio/data/tet4_mat_vertex_t10.vtk b/tests/libtests/meshio/data/tet4_mat_vertex_t10.vtk index 68a5601186..a0a936cbf1 100644 --- a/tests/libtests/meshio/data/tet4_mat_vertex_t10.vtk +++ b/tests/libtests/meshio/data/tet4_mat_vertex_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +material ASCII DATASET UNSTRUCTURED_GRID POINTS 8 double diff --git a/tests/libtests/meshio/data/tet4_surf_cell_t10.vtk b/tests/libtests/meshio/data/tet4_surf_cell_t10.vtk index 0e394235a6..82dbfd2ba3 100644 --- a/tests/libtests/meshio/data/tet4_surf_cell_t10.vtk +++ b/tests/libtests/meshio/data/tet4_surf_cell_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +subdomain_boundary ASCII DATASET UNSTRUCTURED_GRID POINTS 4 double diff --git a/tests/libtests/meshio/data/tet4_surf_vertex_t10.vtk b/tests/libtests/meshio/data/tet4_surf_vertex_t10.vtk index 4226f0d8db..5eb51ac2c8 100644 --- a/tests/libtests/meshio/data/tet4_surf_vertex_t10.vtk +++ b/tests/libtests/meshio/data/tet4_surf_vertex_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +subdomain_boundary ASCII DATASET UNSTRUCTURED_GRID POINTS 4 double diff --git a/tests/libtests/meshio/data/tri3_mat_cell_t10.vtk b/tests/libtests/meshio/data/tri3_mat_cell_t10.vtk index 1fdbc33632..fb31945e3e 100644 --- a/tests/libtests/meshio/data/tri3_mat_cell_t10.vtk +++ b/tests/libtests/meshio/data/tri3_mat_cell_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +material ASCII DATASET UNSTRUCTURED_GRID POINTS 3 double diff --git a/tests/libtests/meshio/data/tri3_mat_vertex_t10.vtk b/tests/libtests/meshio/data/tri3_mat_vertex_t10.vtk index 6cf3589c45..7cffc794df 100644 --- a/tests/libtests/meshio/data/tri3_mat_vertex_t10.vtk +++ b/tests/libtests/meshio/data/tri3_mat_vertex_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +material ASCII DATASET UNSTRUCTURED_GRID POINTS 3 double diff --git a/tests/libtests/meshio/data/tri3_surf_cell_t10.vtk b/tests/libtests/meshio/data/tri3_surf_cell_t10.vtk index 95b7f6a899..177cd6b5b7 100644 --- a/tests/libtests/meshio/data/tri3_surf_cell_t10.vtk +++ b/tests/libtests/meshio/data/tri3_surf_cell_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +subdomain_bc ASCII DATASET UNSTRUCTURED_GRID POINTS 2 double diff --git a/tests/libtests/meshio/data/tri3_surf_vertex_t10.vtk b/tests/libtests/meshio/data/tri3_surf_vertex_t10.vtk index 35d4129bf0..6eebf45eeb 100644 --- a/tests/libtests/meshio/data/tri3_surf_vertex_t10.vtk +++ b/tests/libtests/meshio/data/tri3_surf_vertex_t10.vtk @@ -1,5 +1,5 @@ # vtk DataFile Version 2.0 -domain +subdomain_bc ASCII DATASET UNSTRUCTURED_GRID POINTS 2 double diff --git a/tests/libtests/topology/TestReverseCuthillMcKee.cc b/tests/libtests/topology/TestReverseCuthillMcKee.cc index 51f4db0f01..69dc221a77 100644 --- a/tests/libtests/topology/TestReverseCuthillMcKee.cc +++ b/tests/libtests/topology/TestReverseCuthillMcKee.cc @@ -65,7 +65,7 @@ pylith::topology::TestReverseCuthillMcKee::testReorder(void) { const PetscDM dmOrig = _mesh->getDM(); PetscObjectReference((PetscObject) dmOrig); Mesh meshOrig; - meshOrig.setDM(dmOrig); + meshOrig.setDM(dmOrig, "domain"); ReverseCuthillMcKee::reorder(_mesh);