[cig-commits] r19779 - in short/3D/PyLith/branches/pylith-scecdynrup: . doc/releasenotes examples/3d/hex8 examples/3d/tet4 examples/twocells examples/twocells/twotet4-geoproj libsrc/pylith/faults libsrc/pylith/materials libsrc/pylith/meshio libsrc/pylith/problems libsrc/pylith/topology modulesrc/materials modulesrc/meshio pylith/apps pylith/bc pylith/faults pylith/materials pylith/meshio pylith/problems pylith/topology pylith/utils share tests/2d/frictionslide tests/3d/plasticity tests_auto/2d/quad4 unittests/libtests/materials unittests/pytests/materials unittests/pytests/meshio unittests/pytests/meshio/data
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 14 09:49:00 PDT 2012
Author: brad
Date: 2012-03-14 09:49:00 -0700 (Wed, 14 Mar 2012)
New Revision: 19779
Added:
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/OutputSolnPoints.i
short/3D/PyLith/branches/pylith-scecdynrup/share/debug_malloc.cfg
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/TestOutputSolnPoints.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/data/points.txt
Removed:
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/plasticity/threehex27/
Modified:
short/3D/PyLith/branches/pylith-scecdynrup/README
short/3D/PyLith/branches/pylith-scecdynrup/TODO
short/3D/PyLith/branches/pylith-scecdynrup/configure.ac
short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/announce_v1.6.3.txt
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step01.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step02.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step03.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step04.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step05.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step06.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step07.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step08.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step09.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step10.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step11.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step12.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step13.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step14.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step15.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step16.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step17.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step18.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step19.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step01.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step02.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step03.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step04.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/twocells/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/examples/twocells/twotet4-geoproj/pylithapp.cfg
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/CellFilter.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverLinear.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Field.cc
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPrager3D.i
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPragerPlaneStrain.i
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/meshio.i
short/3D/PyLith/branches/pylith-scecdynrup/pylith/apps/PyLithApp.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/BoundaryCondition.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/Fault.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/DruckerPrager3D.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/Material.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/CellFilter.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/DataWriter.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputManager.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputSolnPoints.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/Solver.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/SolverNonlinear.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/topology/Mesh.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/topology/SubMesh.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/utils/PetscComponent.py
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/plot_friction.py
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/plasticity/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/TestSlipWeakeningShearSliding.py
short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/slipweakening_shear_sliding.cfg
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPrager3D.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPrager3D.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPrager3D.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/data/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/testmeshio.py
Log:
Merge from trunk.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/README
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/README 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/README 2012-03-14 16:49:00 UTC (rev 19779)
@@ -75,6 +75,74 @@
models.
----------------------------------------------------------------------
+Version 1.6.3
+----------------------------------------------------------------------
+
+* Bug fixes
+
+ - Improved error messages for problems encountered during processing
+ of parameters. A backtrace of the object hierarchy is now included
+ to pinpoint in which object the error occurred.
+
+ - Added a line search to the inner friction solve in quasi-static
+ simulations to increase the robustness of the nonlinear
+ solve. Simulations using rate and state friction now converge
+ under a much wider range of circumstances.
+
+ - Fixed bug in updating slip state variable in slip-weakening
+ friction. This caused slight errors in the cumulative slip. We
+ also added a parameter that forces healing to occur in a single
+ time step. This is used to confine slip to a single time step in
+ quasi-static simulations. See examples/3d/hex8/step13.cfg for an
+ example.
+
+ - Tuned parameters in the slip-weakening friction and rate and state
+ friction examples (step13.cfg and step14.cfg, respectively) in
+ examples/3d/hex8 to give stick-slip behavior.
+
+ - Fixed communication issue associated with writing boundary
+ condition information output in parallel.
+
+ - Changed info in Xdmf file for fields that are not scalars,
+ vectors, or tensors so that the each component is extracted,
+ facilitating visualization in ParaView. The corresponding HDF5
+ file remains the same.
+
+ - Added the ability to specify non-derived units (e.g., degree and
+ radian). This is useful in specifying parameters for the
+ Drucker-Prager elastoplastic rheologies. If no units are
+ specified, radians are assumed.
+
+* Internal changes
+
+ - Rate and state friction with ageing law
+
+ The implementation of rate and state friction with ageing law was
+ modified to work better with the iterative solver. We switched to
+ the conventional, unregularized formulation but added a minimum
+ cutoff for the slip rate. Below this cutoff friction has a linear
+ rather than logarithmic dependence on slip rate. As long as this
+ cutoff is close to the SNES solver tolerance, the difference in
+ behavior is negligible while improving the ability of the solver
+ to converge for very small deformations.
+
+KNOWN ISSUES
+
+ The rate and state friction with ageing law has not been tested for
+ dynamic rupture simulations. We plan to run the SCEC Dynamic Rupture
+ benchmarks for rate and state friction as soon as we add a
+ spatial-temporal specification of initial fault tractions, which are
+ required for the benchmark problems.
+
+ Running simulations with more than a million cells and large faults
+ in parallel can result in severe memory imbalances among
+ processors. Some processors around the fault may use 10x more memory
+ than processors away from the fault. We expect this problem to
+ disappear in v1.7 when we switch to new, more efficient Sieve
+ implementation.
+
+
+----------------------------------------------------------------------
Version 1.6.2
----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-03-14 16:49:00 UTC (rev 19779)
@@ -2,15 +2,6 @@
CURRENT ISSUES/PRIORITIES
======================================================================
-* --BUGS--
-
- DruckerPragerPlaneStrain, DruckerPrager3D
- High-frequency oscillation in dynamic simulations
- Extra plastic strain in quasi-static simulations
- Charles' test (which one?) shows same behavior in 2-D and 3-D, so error is
- in both
-
-
* Manual
- Order of tensor components for Xdmf files
@@ -21,6 +12,8 @@
* Approximate for quad4 and hex8 cells
(assume Jacobian is uniform throughout cell =? rhomboid)
(should be much better than nearest)
+ - Examples/Tutorials
+ + Table with concepts/features
* configure
@@ -127,23 +120,6 @@
SECONDARY PRIORITIES
----------------------------------------------------------------------
-* Paper
-
- General paper - focus on fault implementation
- - JGR, GGG
-
- geometry processing (adjusting topology)
- discretization (cohesive cells)
- model (Lagrange multipliers)
- Slip functions
- Fault constitutive models
- solver (saddle point)
- custom preconditioner
- benchmarks
- Savage and Prescott
- Dynamic
- real problem
-
* Fault preconditioner [BRAD]
FaultCohesiveLagrange::calcPreconditioner() [need unit test]
@@ -448,8 +424,6 @@
+ DruckerPragerPlaneStrain (Drucker-Prager plane strain)
+ PowerLawPlaneStrain (power law plane strain )
- 3. Scalable mesh distribution
-
----------------------------------------------------------------------
Release 2.0
----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/configure.ac
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/configure.ac 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/configure.ac 2012-03-14 16:49:00 UTC (rev 19779)
@@ -328,16 +328,12 @@
tests/3d/cyclicfriction/Makefile
tests/3d/cyclicfriction/output/Makefile
tests/3d/plasticity/Makefile
- tests/3d/plasticity/threehex27/Makefile
- tests/3d/plasticity/threehex27/config/Makefile
- tests/3d/plasticity/threehex27/mesh/Makefile
- tests/3d/plasticity/threehex27/results/Makefile
- tests/3d/plasticity/threehex27/spatialdb/Makefile
tests/3d/plasticity/threehex8/Makefile
tests/3d/plasticity/threehex8/config/Makefile
tests/3d/plasticity/threehex8/mesh/Makefile
tests/3d/plasticity/threehex8/results/Makefile
tests/3d/plasticity/threehex8/spatialdb/Makefile
+ tests/3d/plasticity/dynamic/Makefile
tests/topology/Makefile
tests/refinefaulttip/Makefile
doc/Makefile
@@ -372,7 +368,6 @@
examples/greensfns/hex8/gfresponses/Makefile
examples/greensfns/hex8/gfspatialdb/Makefile
examples/twocells/twohex8/Makefile
- examples/twocells/twohex27-cubit/Makefile
examples/twocells/twoquad4/Makefile
examples/twocells/twotet4/Makefile
examples/twocells/twotet4-geoproj/Makefile
Modified: short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/announce_v1.6.3.txt
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/announce_v1.6.3.txt 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/announce_v1.6.3.txt 2012-03-14 16:49:00 UTC (rev 19779)
@@ -21,9 +21,9 @@
* Bug fixes
- - Improved error messages for problems encountered during
- configuration. A backtrace of the object hierarchy is now included
- to pinpoint in which object the rror occurred.
+ - Improved error messages for problems encountered during processing
+ of parameters. A backtrace of the object hierarchy is now included
+ to pinpoint in which object the error occurred.
- Added a line search to the inner friction solve in quasi-static
simulations to increase the robustness of the nonlinear
@@ -35,7 +35,7 @@
also added a parameter that forces healing to occur in a single
time step. This is used to confine slip to a single time step in
quasi-static simulations. See examples/3d/hex8/step13.cfg for an
- exanple.
+ example.
- Tuned parameters in the slip-weakening friction and rate and state
friction examples (step13.cfg and step14.cfg, respectively) in
@@ -44,6 +44,11 @@
- Fixed communication issue associated with writing boundary
condition information output in parallel.
+ - Changed info in Xdmf file for fields that are not scalars,
+ vectors, or tensors so that the each component is extracted,
+ facilitating visualization in ParaView. The corresponding HDF5
+ file remains the same.
+
- Added the ability to specify non-derived units (e.g., degree and
radian). This is useful in specifying parameters for the
Drucker-Prager elastoplastic rheologies. If no units are
Copied: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt (from rev 19778, short/3D/PyLith/trunk/examples/3d/hex8/output_points.txt)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/output_points.txt 2012-03-14 16:49:00 UTC (rev 19779)
@@ -0,0 +1,13 @@
+# This file contains a list of points where we want the solution. The
+# solution will be interpolated to these points, so they need not
+# coincide with vertices of the mesh.
+#
+# We specify the coordinates in the default coordinate system (3-D
+# Cartesian). To specify points in another coordinate system, use the
+# coordsys component of the OutputSolnPoints object. The points will
+# be transformed into the coordinate system of the mesh before
+# interpolation.
+10.0 10.0 -500.0
+10.0 10.0 -1500.0
+10.0 10.0 -2500.0
+10.0 10.0 -3500.0
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step01.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step01.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step01.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step02.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step02.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step02.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step03.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step03.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step03.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step04.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step04.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step04.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step05.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step05.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step05.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step06.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step06.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step06.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step07.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step07.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step07.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step08.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step08.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step08.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,5 +1,7 @@
-# -*- Python -*-
[pylithapp]
+# Permit unknown components so that we can ignore warnings associated
+# with replacing the spatial database for the power-law rheology.
+typos = strict
# ----------------------------------------------------------------------
# PROBLEM DESCRIPTION
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step09.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step09.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step09.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,5 +1,7 @@
-# -*- Python -*-
[pylithapp]
+# Permit unknown components so that we can ignore warnings associated
+# with replacing the spatial database for the power-law rheology.
+typos = strict
# ----------------------------------------------------------------------
# PROBLEM DESCRIPTION
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step10.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step10.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step10.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step11.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step11.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step11.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step12.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step12.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step12.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step13.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step13.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step13.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step14.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step14.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step14.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step15.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step15.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step15.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step16.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step16.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step16.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step17.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step17.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step17.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step18.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step18.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step18.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step19.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step19.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step19.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# ----------------------------------------------------------------------
@@ -36,11 +35,14 @@
[pylithapp.timedependent.implicit]
# Set the output to an array of 2 output managers.
# We will output the solution over the domain and the ground surface.
-output = [domain,subdomain]
+output = [domain,subdomain,points]
# Set subdomain component to OutputSolnSubset (subset of domain).
output.subdomain = pylith.meshio.OutputSolnSubset
+# Set points component to OutputSolnPoints (arbitrary points in domain).
+output.points = pylith.meshio.OutputSolnPoints
+
# Change the total simulation time to 700 years, and set the time
# step to 10 years.
[pylithapp.timedependent.implicit.time_step]
@@ -142,6 +144,17 @@
writer.time_constant = 1.0*year
writer.time_format = %04.0f
+# Give basename for VTK domain output of solution at points.
+[pylithapp.problem.formulation.output.points]
+# We use the default coordinate system for the points. Changing the
+# coordinate system allows the filename with the list of points to use
+# any compatible coordinate system recognized by the spatialdata
+# package.
+reader.filename = output_points.txt
+writer.filename = output/step19-points.vtk
+writer.time_constant = 1.0*year
+writer.time_format = %04.0f
+
# Give basename for VTK output of upper_crust state variables.
[pylithapp.timedependent.materials.upper_crust.output]
cell_filter = pylith.meshio.CellFilterAvgMesh
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step01.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step01.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step01.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# This is not a self-contained simulation configuration file. This
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step02.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step02.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step02.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# This is not a self-contained simulation configuration file. This
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step03.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step03.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step03.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# This is not a self-contained simulation configuration file. This
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step04.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step04.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/tet4/step04.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,4 +1,3 @@
-# -*- Python -*-
[pylithapp]
# This is not a self-contained simulation configuration file. This
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/twocells/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/twocells/Makefile.am 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/twocells/Makefile.am 2012-03-14 16:49:00 UTC (rev 19779)
@@ -18,7 +18,6 @@
SUBDIRS = \
twohex8 \
- twohex27-cubit \
twoquad4 \
twotet4 \
twotet4-geoproj \
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/twocells/twotet4-geoproj/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/twocells/twotet4-geoproj/pylithapp.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/twocells/twotet4-geoproj/pylithapp.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -52,7 +52,7 @@
[pylithapp.mesh_generator.reader.coordsys.projector]
projection = utm
-proj-options = +zone=11
+proj_options = +zone=11
# ----------------------------------------------------------------------
# problem
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -581,11 +581,11 @@
topology::Jacobian* const jacobian,
topology::SolutionFields* const fields)
{ // calcPreconditioner
- assert(0 != precondMatrix);
- assert(0 != jacobian);
- assert(0 != fields);
- assert(0 != _fields);
- assert(0 != _logger);
+ assert(precondMatrix);
+ assert(jacobian);
+ assert(fields);
+ assert(_fields);
+ assert(_logger);
/** We have A = [K L^T]
* [L 0]
@@ -737,7 +737,7 @@
_logger->eventEnd(updateEvent);
#endif
} // for
- err = MatDestroy(&jacobianNP); CHECK_PETSC_ERROR(err);
+ err = MatDestroy(&jacobianNP);CHECK_PETSC_ERROR(err);
PetscLogFlops(numVertices*spaceDim*6);
#if !defined(DETAILED_EVENT_LOGGING)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -154,6 +154,7 @@
_DruckerPrager3D::dbStateVars,
_DruckerPrager3D::numDBStateVars)),
_fitMohrCoulomb(MOHR_COULOMB_INSCRIBED),
+ _allowTensileYield(false),
_calcElasticConstsFn(0),
_calcStressFn(0),
_updateStateVarsFn(0)
@@ -168,6 +169,14 @@
} // destructor
// ----------------------------------------------------------------------
+// Set boolean for whether tensile yield is allowed.
+void
+pylith::materials::DruckerPrager3D::allowTensileYield(const bool flag)
+{ // allowTensileYield
+ _allowTensileYield = flag;
+} // allowTensileYield
+
+// ----------------------------------------------------------------------
// Set fit to Mohr-Coulomb surface.
void
pylith::materials::DruckerPrager3D::fitMohrCoulomb(FitMohrCoulombEnum value)
@@ -246,7 +255,8 @@
switch (_fitMohrCoulomb) { // switch
case MOHR_COULOMB_INSCRIBED: {
const PylithScalar denomFriction = sqrt(3.0) * (3.0 - sin(frictionAngle));
- const PylithScalar denomDilatation = sqrt(3.0) * (3.0 - sin(dilatationAngle));
+ const PylithScalar denomDilatation = sqrt(3.0) *
+ (3.0 - sin(dilatationAngle));
alphaYield = 2.0 * sin(frictionAngle)/denomFriction;
beta = 6.0 * cohesion * cos(frictionAngle)/denomFriction;
alphaFlow = 2.0 * sin(dilatationAngle)/denomDilatation;
@@ -260,7 +270,8 @@
} // MOHR_COULOMB_MIDDLE
case MOHR_COULOMB_CIRCUMSCRIBED: {
const PylithScalar denomFriction = sqrt(3.0) * (3.0 + sin(frictionAngle));
- const PylithScalar denomDilatation = sqrt(3.0) * (3.0 + sin(dilatationAngle));
+ const PylithScalar denomDilatation = sqrt(3.0) *
+ (3.0 + sin(dilatationAngle));
alphaYield = 2.0 * sin(frictionAngle)/denomFriction;
beta = 6.0 * cohesion * cos(frictionAngle)/denomFriction;
alphaFlow = 2.0 * sin(dilatationAngle)/denomDilatation;
@@ -560,12 +571,16 @@
const PylithScalar e22 = totalStrain[1];
const PylithScalar e33 = totalStrain[2];
const PylithScalar meanStrainTpdt = (e11 + e22 + e33)/3.0;
- const PylithScalar meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT - meanStrainInitial;
+ const PylithScalar meanStrainPPTpdt =
+ meanStrainTpdt - meanPlasticStrainT - meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSize] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] -
+ devStrainInitial[2],
totalStrain[3] - devPlasticStrainT[3] - devStrainInitial[3],
totalStrain[4] - devPlasticStrainT[4] - devStrainInitial[4],
totalStrain[5] - devPlasticStrainT[5] - devStrainInitial[5],
@@ -581,9 +596,12 @@
strainPPTpdt[4]/ae + devStressInitial[4],
strainPPTpdt[5]/ae + devStressInitial[5],
};
- const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
- const PylithScalar stressInvar2 = sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
- const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+ const PylithScalar trialMeanStress =
+ meanStrainPPTpdt/am + meanStressInitial;
+ const PylithScalar stressInvar2 =
+ sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
+ const PylithScalar yieldFunction =
+ 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
#if 0 // DEBUGGING
std::cout << "Function _calcStressElastoPlastic: elastic" << std::endl;
std::cout << " alphaYield: " << alphaYield << std::endl;
@@ -603,18 +621,33 @@
const PylithScalar d =
sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
- const PylithScalar plasticMult =
- std::min(sqrt(2.0)*d,
- 2.0 * ae * am * (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
- (6.0 * alphaYield * alphaFlow * ae + am));
+
+ PylithScalar plasticMult = 0.0;
+ if(_allowTensileYield) {
+ plasticMult =
+ std::min(sqrt(2.0)*d,
+ 2.0 * ae * am *
+ (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae)
+ - beta)/ (6.0 * alphaYield * alphaFlow * ae + am));
+ } else {
+ plasticMult =
+ 2.0 * ae * am *
+ (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ } // if/else
+
const PylithScalar meanStressTpdt =
(meanStrainPPTpdt - plasticMult * alphaFlow)/am + meanStressInitial;
PylithScalar deltaDevPlasticStrain = 0.0;
PylithScalar devStressTpdt = 0.0;
- if (d > 0.0) {
+ if (d > 0.0 || !_allowTensileYield) {
for (int iComp=0; iComp < tensorSize; ++iComp) {
- deltaDevPlasticStrain = plasticMult * (strainPPTpdt[iComp] + ae * devStressInitial[iComp]) / (sqrt(2.0) * d);
- devStressTpdt = (strainPPTpdt[iComp] - deltaDevPlasticStrain) / ae + devStressInitial[iComp];
+ deltaDevPlasticStrain =
+ plasticMult * (strainPPTpdt[iComp] + ae * devStressInitial[iComp])/
+ (sqrt(2.0) * d);
+ devStressTpdt =
+ (strainPPTpdt[iComp] - deltaDevPlasticStrain)/ae +
+ devStressInitial[iComp];
stress[iComp] = devStressTpdt + diag[iComp] * meanStressTpdt;
} // for
} else {
@@ -650,12 +683,18 @@
stateVars[s_plasticStrain+5],
};
- const PylithScalar e11 = totalStrain[0] - plasticStrainTpdt[0] - initialStrain[0];
- const PylithScalar e22 = totalStrain[1] - plasticStrainTpdt[1] - initialStrain[1];
- const PylithScalar e33 = totalStrain[2] - plasticStrainTpdt[2] - initialStrain[2];
- const PylithScalar e12 = totalStrain[3] - plasticStrainTpdt[3] - initialStrain[3];
- const PylithScalar e23 = totalStrain[4] - plasticStrainTpdt[4] - initialStrain[4];
- const PylithScalar e13 = totalStrain[5] - plasticStrainTpdt[5] - initialStrain[5];
+ const PylithScalar e11 =
+ totalStrain[0] - plasticStrainTpdt[0] - initialStrain[0];
+ const PylithScalar e22 =
+ totalStrain[1] - plasticStrainTpdt[1] - initialStrain[1];
+ const PylithScalar e33 =
+ totalStrain[2] - plasticStrainTpdt[2] - initialStrain[2];
+ const PylithScalar e12 =
+ totalStrain[3] - plasticStrainTpdt[3] - initialStrain[3];
+ const PylithScalar e23 =
+ totalStrain[4] - plasticStrainTpdt[4] - initialStrain[4];
+ const PylithScalar e13 =
+ totalStrain[5] - plasticStrainTpdt[5] - initialStrain[5];
const PylithScalar traceStrainTpdt = e11 + e22 + e33;
const PylithScalar s123 = lambda * traceStrainTpdt;
@@ -857,9 +896,12 @@
meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSize] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] -
+ devStrainInitial[2],
totalStrain[3] - devPlasticStrainT[3] - devStrainInitial[3],
totalStrain[4] - devPlasticStrainT[4] - devStrainInitial[4],
totalStrain[5] - devPlasticStrainT[5] - devStrainInitial[5],
@@ -904,13 +946,21 @@
(6.0 * alphaYield * alphaFlow * ae + am);
const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
- const PylithScalar testMult = plasticFac *
- (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
- const PylithScalar plasticMult = std::min(sqrt(2.0)*d, testMult);
-
+
+ PylithScalar plasticMult = 0.0;
bool tensileYield = false;
- if (plasticMult == sqrt(2.0) * d)
- tensileYield = true;
+ PylithScalar dFac2 = 0.0;
+ if (_allowTensileYield) {
+ const PylithScalar testMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ tensileYield = (sqrt(2.0)*d < testMult) ? true : false;
+ plasticMult = (tensileYield) ? sqrt(2.0)*d : testMult;
+ dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
+ } else {
+ plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ dFac2 = 1.0/(sqrt(2.0) * d);
+ } // if/else
// Define some constants, vectors, and matrices.
const PylithScalar third = 1.0/3.0;
@@ -929,7 +979,6 @@
strainPPTpdt[4] + ae * devStressInitial[4],
strainPPTpdt[5] + ae * devStressInitial[5]
};
- const PylithScalar dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
PylithScalar dDeltaEdEpsilon = 0.0;
// Compute elasticity matrix.
@@ -1150,9 +1199,12 @@
meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSize] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ totalStrain[2] - meanStrainTpdt - devPlasticStrainT[2] -
+ devStrainInitial[2],
totalStrain[3] - devPlasticStrainT[3] - devStrainInitial[3],
totalStrain[4] - devPlasticStrainT[4] - devStrainInitial[4],
totalStrain[5] - devPlasticStrainT[5] - devStrainInitial[5]
@@ -1169,8 +1221,30 @@
strainPPTpdt[5]/ae + devStressInitial[5]
};
const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
- const PylithScalar stressInvar2 = sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
- const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+ const PylithScalar stressInvar2 =
+ sqrt(0.5 * scalarProduct3D(trialDevStress, trialDevStress));
+ const PylithScalar yieldFunction =
+ 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+
+ if (!_allowTensileYield) {
+ const PylithScalar feasibleFunction =
+ 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+ // If mean stress is too negative to project back to the yield surface,
+ // throw an exception.
+ if (feasibleFunction < 0.0) {
+ std::ostringstream msg;
+ msg << "Infeasible stress state - cannot project back to yield surface.\n"
+ << " alphaYield: " << alphaYield << "\n"
+ << " alphaFlow: " << alphaFlow << "\n"
+ << " beta: " << beta << "\n"
+ << " trialMeanStress: " << trialMeanStress << "\n"
+ << " stressInvar2: " << stressInvar2 << "\n"
+ << " yieldFunction: " << yieldFunction << "\n"
+ << " feasibleFunction: " << feasibleFunction << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+
#if 0 // DEBUGGING
std::cout << "Function _updateStateVarsElastoPlastic:" << std::endl;
std::cout << " alphaYield: " << alphaYield << std::endl;
@@ -1191,13 +1265,25 @@
const PylithScalar d =
sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
scalarProduct3D(devStressInitial, strainPPTpdt) + strainPPTpdtProd);
- const PylithScalar plasticMult =
- std::min(sqrt(2.0)*d,
- 2.0 * ae * am * (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
- (6.0 * alphaYield * alphaFlow * ae + am));
+ const PylithScalar plasticFac = 2.0 * ae * am/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
+ const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
+
+ PylithScalar plasticMult = 0.0;
+ if (_allowTensileYield) {
+ plasticMult = std::min(sqrt(2.0) * d,
+ plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d -
+ beta));
+ } else {
+ plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ } // if/else
+
const PylithScalar deltaMeanPlasticStrain = plasticMult * alphaFlow;
PylithScalar deltaDevPlasticStrain = 0.0;
- if (d > 0.0) {
+ if (d > 0.0 || !_allowTensileYield) {
for (int iComp=0; iComp < tensorSize; ++iComp) {
deltaDevPlasticStrain = plasticMult *(strainPPTpdt[iComp] +
ae * devStressInitial[iComp])/
@@ -1207,7 +1293,8 @@
} // for
} else {
for (int iComp=0; iComp < tensorSize; ++iComp) {
- stateVars[s_plasticStrain+iComp] += diag[iComp] * deltaMeanPlasticStrain;
+ stateVars[s_plasticStrain+iComp] +=
+ diag[iComp] * deltaMeanPlasticStrain;
} // for
} // if/else
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.hh 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPrager3D.hh 2012-03-14 16:49:00 UTC (rev 19779)
@@ -27,7 +27,7 @@
// Include directives ---------------------------------------------------
#include "ElasticMaterial.hh" // ISA ElasticMaterial
-// Powerlaw3D -----------------------------------------------------------
+// DruckerPrager3D ------------------------------------------------------
/** @brief 3-D, isotropic, Drucker-Prager elastic/perfectly plastic material.
*
* The physical properties are specified using density, shear-wave
@@ -67,6 +67,12 @@
*/
void fitMohrCoulomb(FitMohrCoulombEnum value);
+ /** Set flag for whether to allow tensile yield.
+ *
+ * @param flag True if tensile yield is allowed.
+ */
+ void allowTensileYield(const bool flag);
+
/** Set current time step.
*
* @param dt Current time step.
@@ -499,6 +505,9 @@
/// Fit to Mohr Coulomb surface
FitMohrCoulombEnum _fitMohrCoulomb;
+ /// Whether to allow tensile yield
+ bool _allowTensileYield;
+
static const int p_density;
static const int p_mu;
static const int p_lambda;
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -160,6 +160,7 @@
_DruckerPragerPlaneStrain::dbStateVars,
_DruckerPragerPlaneStrain::numDBStateVars)),
_fitMohrCoulomb(MOHR_COULOMB_INSCRIBED),
+ _allowTensileYield(false),
_calcElasticConstsFn(0),
_calcStressFn(0),
_updateStateVarsFn(0)
@@ -174,6 +175,14 @@
} // destructor
// ----------------------------------------------------------------------
+// Set boolean for whether tensile yield is allowed.
+void
+pylith::materials::DruckerPragerPlaneStrain::allowTensileYield(const bool flag)
+{ // allowTensileYield
+ _allowTensileYield = flag;
+} // allowTensileYield
+
+// ----------------------------------------------------------------------
// Set fit to Mohr-Coulomb surface.
void
pylith::materials::DruckerPragerPlaneStrain::fitMohrCoulomb(
@@ -304,8 +313,9 @@
// ----------------------------------------------------------------------
// Nondimensionalize properties.
void
-pylith::materials::DruckerPragerPlaneStrain::_nondimProperties(PylithScalar* const values,
- const int nvalues) const
+pylith::materials::DruckerPragerPlaneStrain::_nondimProperties(
+ PylithScalar* const values,
+ const int nvalues) const
{ // _nondimProperties
assert(_normalizer);
assert(values);
@@ -329,8 +339,9 @@
// ----------------------------------------------------------------------
// Dimensionalize properties.
void
-pylith::materials::DruckerPragerPlaneStrain::_dimProperties(PylithScalar* const values,
- const int nvalues) const
+pylith::materials::DruckerPragerPlaneStrain::_dimProperties(
+ PylithScalar* const values,
+ const int nvalues) const
{ // _dimProperties
assert(_normalizer);
assert(values);
@@ -578,13 +589,16 @@
// Values for current time step
const PylithScalar meanStrainTpdt = (totalStrain[0] + totalStrain[1])/3.0;
- const PylithScalar meanStrainPPTpdt = meanStrainTpdt - meanPlasticStrainT - meanStrainInitial;
+ const PylithScalar meanStrainPPTpdt =
+ meanStrainTpdt - meanPlasticStrainT - meanStrainInitial;
// devStrainPPTpdt
const PylithScalar strainPPTpdt[tensorSizePS] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
totalStrain[2] - devPlasticStrainT[3] - devStrainInitial[3]
};
@@ -596,11 +610,12 @@
strainPPTpdt[2]/ae + devStressInitial[2],
strainPPTpdt[3]/ae + devStressInitial[3]
};
- const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
+ const PylithScalar trialMeanStress =
+ meanStrainPPTpdt/am + meanStressInitial;
const PylithScalar stressInvar2 =
sqrt(0.5 * scalarProduct2DPS(trialDevStress, trialDevStress));
- const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
- stressInvar2 - beta;
+ const PylithScalar yieldFunction =
+ 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
#if 0 // DEBUGGING
std::cout << "Function _calcStressElastoPlastic: elastic" << std::endl;
std::cout << " alphaYield: " << alphaYield << std::endl;
@@ -621,11 +636,23 @@
sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
scalarProduct2DPS(devStressInitial, strainPPTpdt) +
strainPPTpdtProd);
- const PylithScalar plasticMult =
- std::min(sqrt(2.0)*d,
- 2.0 * ae * am * (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
- (6.0 * alphaYield * alphaFlow * ae + am));
- const PylithScalar meanStressTpdt = (meanStrainPPTpdt - plasticMult * alphaFlow) / am + meanStressInitial;
+
+ PylithScalar plasticMult = 0.0;
+ if(_allowTensileYield) {
+ plasticMult =
+ std::min(sqrt(2.0)*d,
+ 2.0 * ae * am *
+ (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae)
+ - beta)/ (6.0 * alphaYield * alphaFlow * ae + am));
+ } else {
+ plasticMult =
+ 2.0 * ae * am *
+ (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ } // if/else
+
+ const PylithScalar meanStressTpdt =
+ (meanStrainPPTpdt - plasticMult * alphaFlow) / am + meanStressInitial;
PylithScalar deltaDevPlasticStrain = 0.0;
PylithScalar devStressTpdt = 0.0;
PylithScalar totalStress[tensorSizePS];
@@ -653,10 +680,14 @@
<< std::endl;
#endif
- if (d > 0.0) {
+ if (d > 0.0 || !_allowTensileYield) {
for (int iComp=0; iComp < tensorSizePS; ++iComp) {
- deltaDevPlasticStrain = plasticMult * (strainPPTpdt[iComp] + ae * devStressInitial[iComp]) / (sqrt(2.0) * d);
- devStressTpdt = (strainPPTpdt[iComp] - deltaDevPlasticStrain)/ae + devStressInitial[iComp];
+ deltaDevPlasticStrain =
+ plasticMult * (strainPPTpdt[iComp] + ae * devStressInitial[iComp])/
+ (sqrt(2.0) * d);
+ devStressTpdt =
+ (strainPPTpdt[iComp] - deltaDevPlasticStrain)/ae +
+ devStressInitial[iComp];
totalStress[iComp] = devStressTpdt + diag[iComp] * meanStressTpdt;
} // for
} else {
@@ -883,9 +914,12 @@
meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSizePS] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ - meanStrainTpdt - devPlasticStrainT[2] -
+ devStrainInitial[2],
totalStrain[2] - devPlasticStrainT[3] - devStrainInitial[3],
};
@@ -926,13 +960,21 @@
(6.0 * alphaYield * alphaFlow * ae + am);
const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
- const PylithScalar testMult = plasticFac *
- (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
- const PylithScalar plasticMult = std::min(sqrt(2.0)*d, testMult);
+ PylithScalar plasticMult = 0.0;
bool tensileYield = false;
- if (plasticMult == sqrt(2.0) * d)
- tensileYield = true;
+ PylithScalar dFac2 = 0.0;
+ if (_allowTensileYield) {
+ const PylithScalar testMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ tensileYield = (sqrt(2.0) * d < testMult) ? true: false;
+ plasticMult = tensileYield ? sqrt(2.0) * d : testMult;
+ dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
+ } else {
+ plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ dFac2 = 1.0/(sqrt(2.0) * d);
+ } // if/else
// Define some constants, vectors, and matrices.
const PylithScalar third = 1.0/3.0;
@@ -946,7 +988,6 @@
strainPPTpdt[3] + ae * devStressInitial[3],
};
- const PylithScalar dFac2 = (d > 0.0) ? 1.0/(sqrt(2.0) * d) : 0.0;
PylithScalar dDeltaEdEpsilon = 0.0;
// Compute elasticity matrix.
@@ -1132,9 +1173,11 @@
meanStrainInitial;
const PylithScalar strainPPTpdt[tensorSizePS] = {
- totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] - devStrainInitial[0],
- totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] - devStrainInitial[1],
- 0.0 - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
+ totalStrain[0] - meanStrainTpdt - devPlasticStrainT[0] -
+ devStrainInitial[0],
+ totalStrain[1] - meanStrainTpdt - devPlasticStrainT[1] -
+ devStrainInitial[1],
+ - meanStrainTpdt - devPlasticStrainT[2] - devStrainInitial[2],
totalStrain[2] - devPlasticStrainT[3] - devStrainInitial[3]
};
@@ -1148,8 +1191,28 @@
const PylithScalar trialMeanStress = meanStrainPPTpdt/am + meanStressInitial;
const PylithScalar stressInvar2 =
sqrt(0.5 * scalarProduct2DPS(trialDevStress, trialDevStress));
- const PylithScalar yieldFunction = 3.0 * alphaYield * trialMeanStress +
- stressInvar2 - beta;
+ const PylithScalar yieldFunction =
+ 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+
+ if (!_allowTensileYield) {
+ const PylithScalar feasibleFunction =
+ 3.0 * alphaYield * trialMeanStress + stressInvar2 - beta;
+ // If mean stress is too negative to project back to the yield surface,
+ // throw an exception.
+ if (feasibleFunction < 0.0) {
+ std::ostringstream msg;
+ msg << "Infeasible stress state - cannot project back to yield surface.\n"
+ << " alphaYield: " << alphaYield << "\n"
+ << " alphaFlow: " << alphaFlow << "\n"
+ << " beta: " << beta << "\n"
+ << " trialMeanStress: " << trialMeanStress << "\n"
+ << " stressInvar2: " << stressInvar2 << "\n"
+ << " yieldFunction: " << yieldFunction << "\n"
+ << " feasibleFunction: " << feasibleFunction << "\n";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+
#if 0 // DEBUGGING
std::cout << "Function _updateStateVarsElastoPlastic:" << std::endl;
std::cout << " alphaYield: " << alphaYield << std::endl;
@@ -1171,13 +1234,25 @@
sqrt(ae * ae * devStressInitialProd + 2.0 * ae *
scalarProduct2DPS(devStressInitial, strainPPTpdt) +
strainPPTpdtProd);
- const PylithScalar plasticMult =
- std::min(sqrt(2.0)*d,
- 2.0 * ae * am * (3.0 * alphaYield * meanStrainPPTpdt/am + d/(sqrt(2.0) * ae) - beta)/
- (6.0 * alphaYield * alphaFlow * ae + am));
+ const PylithScalar plasticFac = 2.0 * ae * am/
+ (6.0 * alphaYield * alphaFlow * ae + am);
+ const PylithScalar meanStrainFac = 3.0 * alphaYield/am;
+ const PylithScalar dFac = 1.0/(sqrt(2.0) * ae);
+
+ PylithScalar plasticMult = 0.0;
+ if (_allowTensileYield) {
+ plasticMult = std::min(sqrt(2.0) * d,
+ plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d -
+ beta));
+ } else {
+ plasticMult = plasticFac *
+ (meanStrainFac * meanStrainPPTpdt + dFac * d - beta);
+ } // if/else
+
const PylithScalar deltaMeanPlasticStrain = plasticMult * alphaFlow;
PylithScalar deltaDevPlasticStrain = 0.0;
- if (d > 0.0) {
+ if (d > 0.0 || !_allowTensileYield) {
for (int iComp=0; iComp < tensorSizePS; ++iComp) {
deltaDevPlasticStrain = plasticMult *(strainPPTpdt[iComp] +
ae * devStressInitial[iComp])/
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/materials/DruckerPragerPlaneStrain.hh 2012-03-14 16:49:00 UTC (rev 19779)
@@ -67,6 +67,12 @@
*/
void fitMohrCoulomb(FitMohrCoulombEnum value);
+ /** Set flag for whether to allow tensile yield.
+ *
+ * @param flag True if tensile yield is allowed.
+ */
+ void allowTensileYield(const bool flag);
+
/** Set current time step.
*
* @param dt Current time step.
@@ -499,6 +505,9 @@
/// Fit to Mohr Coulomb surface
FitMohrCoulombEnum _fitMohrCoulomb;
+ /// Whether to allow tensile yield
+ bool _allowTensileYield;
+
static const int p_density;
static const int p_mu;
static const int p_lambda;
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/CellFilter.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/CellFilter.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/CellFilter.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -42,7 +42,7 @@
pylith::meshio::CellFilter<mesh_type, field_type>::CellFilter(const CellFilter& f) :
_quadrature(0)
{ // copy constructor
- if (0 != f._quadrature)
+ if (f._quadrature)
_quadrature = new feassemble::Quadrature<mesh_type>(*f._quadrature);
} // copy constructor
@@ -62,7 +62,7 @@
pylith::meshio::CellFilter<mesh_type, field_type>::quadrature(const feassemble::Quadrature<mesh_type>* q)
{ // quadrature
delete _quadrature;
- _quadrature = (0 != q) ? new feassemble::Quadrature<mesh_type>(*q) : 0;
+ _quadrature = (q) ? new feassemble::Quadrature<mesh_type>(*q) : 0;
} // quadrature
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -101,6 +101,16 @@
assert(csMesh);
assert(csMesh->spaceDim() == spaceDim);
+#if 1 // DEBUGGING
+ std::cout << "OUTPUT SOLN POINTS" << std::endl;
+ for (int i=0; i < numPoints; ++i) {
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ std::cout << " " << points[i*spaceDim+iDim];
+ } // for
+ std::cout << "\n";
+ } // for
+#endif
+
scalar_array pointsArray(points, numPoints*spaceDim);
int_array cells(numPoints);
for (int i=0; i < numPoints; ++i)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -61,9 +61,9 @@
#if 0 // :KLUDGE: Assume Solver deallocates matrix.
PetscErrorCode err = 0;
- if (0 != _customConstraintPCMat) {
- err = PetscObjectDereference((PetscObject) _customConstraintPCMat);
- _customConstraintPCMat = 0; CHECK_PETSC_ERROR(err);
+ if (_customConstraintPCMat) {
+ err = PetscObjectDereference((PetscObject) _customConstraintPCMat);CHECK_PETSC_ERROR(err);
+ _customConstraintPCMat = 0;
} // if
#else
_customConstraintPCMat = 0;
@@ -319,7 +319,7 @@
// Assemble jacobian.
_jacobian->assemble("final_assembly");
- if (0 != _customConstraintPCMat) {
+ if (_customConstraintPCMat) {
// Recalculate preconditioner.
numIntegrators = _meshIntegrators.size();
for (int i=0; i < numIntegrators; ++i)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -89,12 +89,14 @@
{ // deallocate
_formulation = 0; // Handle only, do not manage memory.
delete _logger; _logger = 0;
- if (0 != _jacobianPC) {
- PetscErrorCode err = MatDestroy(&_jacobianPC);CHECK_PETSC_ERROR(err);
- } // if
- if (0 != _jacobianPCFault) {
- PetscErrorCode err = MatDestroy(&_jacobianPCFault);CHECK_PETSC_ERROR(err);
- } // if
+
+ PetscErrorCode err = 0;
+ err = MatDestroy(&_jacobianPC);CHECK_PETSC_ERROR(err);
+ err = MatDestroy(&_jacobianPCFault);CHECK_PETSC_ERROR(err);
+
+ _ctx.pc = 0; // KSP PC (managed separately)
+ _ctx.A = 0; // Jacobian (managed separately)
+ _ctx.faultA = 0; // Handle to _jacobianPCFault
} // deallocate
// ----------------------------------------------------------------------
@@ -104,7 +106,7 @@
const topology::Jacobian& jacobian,
Formulation* const formulation)
{ // initialize
- assert(0 != formulation);
+ assert(formulation);
_formulation = formulation;
// Make global preconditioner matrix
@@ -123,6 +125,7 @@
PylithInt M, N, m, n;
PetscErrorCode err = 0;
+ err = MatDestroy(&_jacobianPC);CHECK_PETSC_ERROR(err);
err = MatGetSize(jacobianMat, &M, &N);CHECK_PETSC_ERROR(err);
err = MatGetLocalSize(jacobianMat, &m, &n);CHECK_PETSC_ERROR(err);
err = MatCreateShell(fields.mesh().comm(), m, n, M, N, &_ctx, &_jacobianPC);
@@ -178,9 +181,7 @@
solutionSection, spaceDim);
assert(!lagrangeGlobalOrder.isNull());
- if (_jacobianPCFault) {
- err = MatDestroy(&_jacobianPCFault); CHECK_PETSC_ERROR(err);
- } // if
+ err = MatDestroy(&_jacobianPCFault); CHECK_PETSC_ERROR(err);
PylithInt nrows = lagrangeGlobalOrder->getLocalSize();
PylithInt ncols = nrows;
@@ -190,20 +191,11 @@
err = MatSetType(_jacobianPCFault, MATAIJ);
err = MatSetFromOptions(_jacobianPCFault); CHECK_PETSC_ERROR(err);
-#if 1
// Allocate just the diagonal.
err = MatSeqAIJSetPreallocation(_jacobianPCFault, 1,
PETSC_NULL); CHECK_PETSC_ERROR(err);
err = MatMPIAIJSetPreallocation(_jacobianPCFault, 1, PETSC_NULL,
0, PETSC_NULL); CHECK_PETSC_ERROR(err);
-#else
- // Allocate full matrix (overestimate).
- err = MatSeqAIJSetPreallocation(_jacobianPCFault, ncols,
- PETSC_NULL); CHECK_PETSC_ERROR(err);
- err = MatMPIAIJSetPreallocation(_jacobianPCFault, ncols, PETSC_NULL,
- 0, PETSC_NULL); CHECK_PETSC_ERROR(err);
-#endif
-
// Set preconditioning matrix in formulation
formulation->customPCMatrix(_jacobianPCFault);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverLinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverLinear.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverLinear.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -54,10 +54,7 @@
{ // deallocate
Solver::deallocate();
- if (0 != _ksp) {
- PetscErrorCode err = KSPDestroy(&_ksp); _ksp = 0;
- CHECK_PETSC_ERROR(err);
- } // if
+ PetscErrorCode err = KSPDestroy(&_ksp);CHECK_PETSC_ERROR(err);
} // deallocate
// ----------------------------------------------------------------------
@@ -68,23 +65,20 @@
const topology::Jacobian& jacobian,
Formulation* formulation)
{ // initialize
- assert(0 != formulation);
+ assert(formulation);
_initializeLogger();
Solver::initialize(fields, jacobian, formulation);
PetscErrorCode err = 0;
- if (0 != _ksp) {
- err = KSPDestroy(&_ksp); _ksp = 0;
- CHECK_PETSC_ERROR(err);
- } // if
- err = KSPCreate(fields.mesh().comm(), &_ksp); CHECK_PETSC_ERROR(err);
- err = KSPSetInitialGuessNonzero(_ksp, PETSC_FALSE); CHECK_PETSC_ERROR(err);
- err = KSPSetFromOptions(_ksp); CHECK_PETSC_ERROR(err);
+ err = KSPDestroy(&_ksp);CHECK_PETSC_ERROR(err);
+ err = KSPCreate(fields.mesh().comm(), &_ksp);CHECK_PETSC_ERROR(err);
+ err = KSPSetInitialGuessNonzero(_ksp, PETSC_FALSE);CHECK_PETSC_ERROR(err);
+ err = KSPSetFromOptions(_ksp);CHECK_PETSC_ERROR(err);
if (formulation->splitFields()) {
PetscPC pc = 0;
- err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
+ err = KSPGetPC(_ksp, &pc);CHECK_PETSC_ERROR(err);
_setupFieldSplit(&pc, formulation, jacobian, fields);
} // if
} // initialize
@@ -97,15 +91,17 @@
topology::Jacobian* jacobian,
const topology::Field<topology::Mesh>& residual)
{ // solve
- assert(0 != solution);
- assert(0 != jacobian);
- assert(0 != _formulation);
+ assert(solution);
+ assert(jacobian);
+ assert(_formulation);
const int setupEvent = _logger->eventId("SoLi setup");
const int solveEvent = _logger->eventId("SoLi solve");
const int scatterEvent = _logger->eventId("SoLi scatter");
_logger->eventBegin(scatterEvent);
+ //PetscFunctionBegin; // DEBUGGING
+
// Update PetscVector view of field.
residual.scatterSectionToVector();
@@ -116,10 +112,10 @@
const PetscMat jacobianMat = jacobian->matrix();
if (!jacobian->valuesChanged()) {
err = KSPSetOperators(_ksp, jacobianMat, jacobianMat,
- SAME_PRECONDITIONER); CHECK_PETSC_ERROR(err);
+ SAME_PRECONDITIONER);CHECK_PETSC_ERROR(err);
} else {
err = KSPSetOperators(_ksp, jacobianMat, jacobianMat,
- DIFFERENT_NONZERO_PATTERN); CHECK_PETSC_ERROR(err);
+ DIFFERENT_NONZERO_PATTERN);CHECK_PETSC_ERROR(err);
} // else
jacobian->resetValuesChanged();
@@ -144,12 +140,10 @@
assert(solutionSection->getNumSpaces() == num);
MatStructure flag;
- err = KSPGetOperators(ksps[num-1], &A,
- PETSC_NULL, &flag); CHECK_PETSC_ERROR(err);
- err = PetscObjectReference((PetscObject) A); CHECK_PETSC_ERROR(err);
- err = KSPSetOperators(ksps[num-1], A, _jacobianPCFault,
- flag); CHECK_PETSC_ERROR(err);
- err = PetscFree(ksps); CHECK_PETSC_ERROR(err);
+ err = KSPGetOperators(ksps[num-1], &A, PETSC_NULL, &flag);CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) A);CHECK_PETSC_ERROR(err);
+ err = KSPSetOperators(ksps[num-1], A, _jacobianPCFault, flag);CHECK_PETSC_ERROR(err);
+ err = PetscFree(ksps);CHECK_PETSC_ERROR(err);
} // if
#endif
@@ -172,6 +166,8 @@
// Update rate fields to be consistent with current solution.
_formulation->calcRateFields();
+
+ //PetscFunctionReturnVoid(); // DEBUGGING
} // solve
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -62,10 +62,7 @@
{ // deallocate
Solver::deallocate();
- if (0 != _snes) {
- PetscErrorCode err = SNESDestroy(&_snes); _snes = 0;
- CHECK_PETSC_ERROR(err);
- } // if
+ PetscErrorCode err = SNESDestroy(&_snes);CHECK_PETSC_ERROR(err);
} // deallocate
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Field.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/Field.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -75,17 +75,10 @@
for (typename scatter_map_type::iterator s_iter=_scatters.begin();
s_iter != scattersEnd;
++s_iter) {
- if (s_iter->second.vector) {
- err = VecDestroy(&s_iter->second.vector);CHECK_PETSC_ERROR(err);
- } // if
- if (s_iter->second.scatter) {
- err = VecScatterDestroy(&s_iter->second.scatter);CHECK_PETSC_ERROR(err);
- } // if
-
- if (s_iter->second.scatterVec) {
- err = VecDestroy(&s_iter->second.scatterVec);CHECK_PETSC_ERROR(err);
- } // if
+ err = VecDestroy(&s_iter->second.vector);CHECK_PETSC_ERROR(err);
+ err = VecScatterDestroy(&s_iter->second.scatter);CHECK_PETSC_ERROR(err);
+ err = VecDestroy(&s_iter->second.scatterVec);CHECK_PETSC_ERROR(err);
} // for
_scatters.clear();
} // deallocate
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPrager3D.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPrager3D.i 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPrager3D.i 2012-03-14 16:49:00 UTC (rev 19779)
@@ -50,6 +50,12 @@
* @param value Mohr-Coulomb surface match type.
*/
void fitMohrCoulomb(FitMohrCoulombEnum value);
+
+ /** Set boolean for whether tensile yield is allowed.
+ *
+ * @param flag True if tensile yield is allowed, false otherwise.
+ */
+ void allowTensileYield(const bool flag);
/** Set current time step.
*
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPragerPlaneStrain.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPragerPlaneStrain.i 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/materials/DruckerPragerPlaneStrain.i 2012-03-14 16:49:00 UTC (rev 19779)
@@ -51,6 +51,12 @@
*/
void fitMohrCoulomb(FitMohrCoulombEnum value);
+ /** Set boolean for whether tensile yield is allowed.
+ *
+ * @param flag True if tensile yield is allowed, false otherwise.
+ */
+ void allowTensileYield(const bool flag);
+
/** Set current time step.
*
* @param dt Current time step.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/Makefile.am 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/Makefile.am 2012-03-14 16:49:00 UTC (rev 19779)
@@ -37,7 +37,8 @@
DataWriter.i \
DataWriterVTK.i \
OutputManager.i \
- OutputSolnSubset.i
+ OutputSolnSubset.i \
+ OutputSolnPoints.i
swig_generated = \
meshio_wrap.cxx \
Copied: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/OutputSolnPoints.i (from rev 19778, short/3D/PyLith/trunk/modulesrc/meshio/OutputSolnPoints.i)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/OutputSolnPoints.i (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/OutputSolnPoints.i 2012-03-14 16:49:00 UTC (rev 19779)
@@ -0,0 +1,125 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file modulesrc/meshio/OutputSolnPoints.i
+ *
+ * @brief Python interface to C++ OutputSolnPoints object.
+ */
+
+namespace pylith {
+ namespace meshio {
+
+%template(_MeshOutputManager) pylith::meshio::OutputManager<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >;
+
+
+ class pylith::meshio::OutputSolnPoints : public OutputManager<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >
+ { // OutputSolnPoints
+
+ // PUBLIC METHODS ///////////////////////////////////////////////////
+ public :
+
+ /// Constructor
+ OutputSolnPoints(void);
+
+ /// Destructor
+ ~OutputSolnPoints(void);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(void);
+
+ /** Get mesh associated with points.
+ *
+ * @returns Mesh associated with points.
+ */
+ const pylith::topology::Mesh& pointsMesh(void);
+
+ /** Setup interpolator.
+ *
+ * @param mesh Domain mesh.
+ * @param points Array of coordinates for points [numPoints*spaceDim].
+ * @param numPoints Number of points.
+ * @param spaceDim Spatial dimension for coordinates.
+ */
+ %apply(PylithScalar* IN_ARRAY2, int DIM1, int DIM2) {
+ (const PylithScalar* points,
+ const int numPoints,
+ const int spaceDim)
+ };
+ void setupInterpolator(pylith::topology::Mesh* mesh,
+ const PylithScalar* points,
+ const int numPoints,
+ const int spaceDim);
+ %clear(const PylithScalar* points, const int numPoints, const int spaceDim);
+
+ /** Prepare for output.
+ *
+ * @param mesh Finite-element mesh object.
+ * @param numTimeSteps Expected number of time steps.
+ * @param label Name of label defining cells to include in output
+ * (=0 means use all cells in mesh).
+ * @param labelId Value of label defining which cells to include.
+ */
+ void open(const pylith::topology::Mesh& mesh,
+ const int numTimeSteps,
+ const char* label =0,
+ const int labelId =0);
+
+ /** Setup file for writing fields at time step.
+ *
+ * @param t Time of time step.
+ * @param mesh Finite-element mesh object.
+ * @param label Name of label defining cells to include in output
+ * (=0 means use all cells in mesh).
+ * @param labelId Value of label defining which cells to include.
+ */
+ void openTimeStep(const PylithScalar t,
+ const pylith::topology::Mesh& mesh,
+ const char* label =0,
+ const int labelId =0);
+
+ /** Append finite-element vertex field to file.
+ *
+ * @param t Time associated with field.
+ * @param field Vertex field.
+ * @param mesh Mesh for output.
+ */
+ void appendVertexField(const PylithScalar t,
+ pylith::topology::Field<topology::Mesh>& field,
+ const pylith::topology::Mesh& mesh);
+
+ /** Append finite-element cell field to file.
+ *
+ * @param t Time associated with field.
+ * @param field Cell field.
+ * @param label Name of label defining cells to include in output
+ * (=0 means use all cells in mesh).
+ * @param labelId Value of label defining which cells to include.
+ */
+ void appendCellField(const PylithScalar t,
+ pylith::topology::Field<topology::Mesh>& field,
+ const char* label =0,
+ const int labelId =0);
+
+ }; // OutputSolnPoints
+
+ } // meshio
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/meshio.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/meshio.i 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/meshio/meshio.i 2012-03-14 16:49:00 UTC (rev 19779)
@@ -37,6 +37,7 @@
#include "pylith/meshio/DataWriterVTK.hh"
#include "pylith/meshio/OutputManager.hh"
#include "pylith/meshio/OutputSolnSubset.hh"
+#include "pylith/meshio/OutputSolnPoints.hh"
#if defined(ENABLE_HDF5)
#include "pylith/meshio/DataWriterHDF5.hh"
#include "pylith/meshio/DataWriterHDF5Ext.hh"
@@ -60,6 +61,15 @@
%include "typemaps.i"
%include "../include/scalartypemaps.i"
+// Numpy interface stuff
+%{
+#define SWIG_FILE_WITH_INIT
+%}
+%include "../include/numpy.i"
+%init %{
+import_array();
+%}
+
// Interfaces
%include "MeshIOObj.i"
%include "MeshIOAscii.i"
@@ -77,6 +87,7 @@
%include "DataWriterVTK.i"
%include "OutputManager.i"
%include "OutputSolnSubset.i"
+%include "OutputSolnPoints.i"
#if defined(ENABLE_HDF5)
%include "DataWriterHDF5.i"
%include "DataWriterHDF5Ext.i"
@@ -101,6 +112,7 @@
%template(MeshDataWriterVTK) pylith::meshio::DataWriterVTK<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >;
%template(SubMeshDataWriterVTK) pylith::meshio::DataWriterVTK<pylith::topology::SubMesh, pylith::topology::Field<pylith::topology::Mesh> >;
%template(SubSubMeshDataWriterVTK) pylith::meshio::DataWriterVTK<pylith::topology::SubMesh, pylith::topology::Field<pylith::topology::SubMesh> >;
+%template(PointsDataWriterVTK) pylith::meshio::DataWriterVTK<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >;
#if defined(ENABLE_HDF5)
%template(MeshDataWriterHDF5) pylith::meshio::DataWriterHDF5<pylith::topology::Mesh, pylith::topology::Field<pylith::topology::Mesh> >;
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/apps/PyLithApp.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/apps/PyLithApp.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/apps/PyLithApp.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -144,6 +144,7 @@
self.mesher = self.inventory.mesher
self.problem = self.inventory.problem
self.perfLogger = self.inventory.perfLogger
+ self.typos = self.inventory.typos
import journal
self._debug = journal.debug(self.name)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/BoundaryCondition.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/BoundaryCondition.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/BoundaryCondition.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -155,6 +155,14 @@
return
+ def _cleanup(self):
+ """
+ Deallocate locally managed data structures.
+ """
+ self.deallocate()
+ return
+
+
def _createModuleObj(self):
"""
Call constructor for module object for access to C++ object.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/Fault.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/Fault.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/Fault.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -261,4 +261,12 @@
return
+ def _cleanup(self):
+ """
+ Deallocate locally managed data structures.
+ """
+ self.deallocate()
+ return
+
+
# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/DruckerPrager3D.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/DruckerPrager3D.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/DruckerPrager3D.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -57,6 +57,8 @@
##
## \b Properties
## @li \b fit_mohr_coulomb Fit to Mohr-Coulomb yield surface.
+ ## @li \b allow_tensile_yield If true, allow yield beyond tensile strength;
+ ## otherwise an exception occurs for excessive tensile sttress.
##
## \b Facilities
## @li None
@@ -67,6 +69,10 @@
fitMohrCoulomb = pyre.inventory.str("fit_mohr_coulomb", default="inscribed",
validator=validateFitMohrCoulomb)
fitMohrCoulomb.meta['tip'] = "Fit to Mohr-Coulomb yield surface."
+ allowTensileYield = pyre.inventory.bool("allow_tensile_yield",
+ default=False)
+ allowTensileYield.meta['tip'] = "If true, allow yield beyond tensile " \
+ "strength; otherwise, an exception is thrown."
# PUBLIC METHODS /////////////////////////////////////////////////////
@@ -104,6 +110,8 @@
else:
raise ValueError("Unknown fit to Mohr-Coulomb yield surface.")
ModuleDruckerPrager3D.fitMohrCoulomb(self, fitEnum)
+ ModuleDruckerPrager3D.allowTensileYield(self,
+ self.inventory.allowTensileYield)
return
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/Material.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/Material.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/materials/Material.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -237,4 +237,12 @@
return
+ def _cleanup(self):
+ """
+ Deallocate locally managed data structures.
+ """
+ self.deallocate()
+ return
+
+
# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/CellFilter.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/CellFilter.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/CellFilter.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -96,6 +96,14 @@
return
+ def _cleanup(self):
+ """
+ Deallocate locally managed data structures.
+ """
+ self.deallocate()
+ return
+
+
# FACTORIES ////////////////////////////////////////////////////////////
def output_cell_filter():
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/DataWriter.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/DataWriter.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/DataWriter.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -58,4 +58,14 @@
return
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _cleanup(self):
+ """
+ Deallocate locally managed data structures.
+ """
+ self.deallocate()
+ return
+
+
# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputManager.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputManager.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputManager.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -396,6 +396,14 @@
return
+ def _cleanup(self):
+ """
+ Deallocate PETSc and local data structures.
+ """
+ self.deallocate()
+ return
+
+
def _open(self):
raise NotImplementedError("Implement _open() in derived class.")
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputSolnPoints.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputSolnPoints.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/meshio/OutputSolnPoints.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -98,15 +98,6 @@
return
- def verifyConfiguration(self, mesh):
- """
- Verify compatibility of configuration.
- """
- OutputManager.verifyConfiguration(self, mesh)
- ModuleOutputSolnPoints.verifyConfiguration(self, mesh)
- return
-
-
def initialize(self, mesh, normalizer):
"""
Initialize output manager.
@@ -114,10 +105,19 @@
logEvent = "%sinit" % self._loggingPrefix
self._eventLogger.eventBegin(logEvent)
+ OutputManager.initialize(self, normalizer)
+
+ # Read points
points = self.reader.read()
+
+ # Convert to mesh coordinate system
+ from spatialdata.geocoords.Converter import convert
+ convert(points, mesh.coordsys(), self.coordsys)
+ print "LENGTH SCALE",normalizer.lengthScale()
+ points /= normalizer.lengthScale().value
+
ModuleOutputSolnPoints.setupInterpolator(self, mesh, points)
- self.mesh = ModuleOutputSolnPoints.createPointsMesh(self)
- OutputManager.initialize(self, normalizer)
+ self.mesh = ModuleOutputSolnPoints.pointsMesh(self)
self._eventLogger.eventEnd(logEvent)
return
@@ -153,8 +153,6 @@
"""
try:
OutputManager._configure(self)
- ModuleOutputSolnPoints.label(self, self.label)
- ModuleOutputSolnPoints.coordsys(self, self.inventory.coordsys)
ModuleOutputSolnPoints.writer(self, self.inventory.writer)
from pylith.utils.NullComponent import NullComponent
if not isinstance(self.inventory.vertexFilter, NullComponent):
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/Solver.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/Solver.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/Solver.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -96,6 +96,14 @@
return
+ def _cleanup(self):
+ """
+ Deallocate PETSc and local data structures.
+ """
+ self.deallocate()
+ return
+
+
# FACTORIES ////////////////////////////////////////////////////////////
def solver():
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/SolverNonlinear.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/SolverNonlinear.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/problems/SolverNonlinear.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -69,14 +69,6 @@
return
- def _cleanup(self):
- """
- Deallocate PETSc and local data structures.
- """
- self.deallocate()
- return
-
-
# FACTORIES ////////////////////////////////////////////////////////////
def solver():
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/topology/Mesh.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/topology/Mesh.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/topology/Mesh.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -85,4 +85,14 @@
return groups
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _cleanup(self):
+ """
+ Deallocate locally managed data structures.
+ """
+ self.deallocate()
+ return
+
+
# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/topology/SubMesh.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/topology/SubMesh.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/topology/SubMesh.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -54,4 +54,14 @@
return Communicator(ModuleSubMesh.comm(self))
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _cleanup(self):
+ """
+ Deallocate locally managed data structures.
+ """
+ self.deallocate()
+ return
+
+
# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/utils/PetscComponent.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/utils/PetscComponent.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/utils/PetscComponent.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -70,7 +70,7 @@
# Facility arrays are not PetscComponents but have components().
elif hasattr(component, "components"):
- for subcomponent in self.components():
+ for subcomponent in component.components():
if isinstance(subcomponent, PetscComponent):
subcomponent.cleanup()
Copied: short/3D/PyLith/branches/pylith-scecdynrup/share/debug_malloc.cfg (from rev 19778, short/3D/PyLith/trunk/share/debug_malloc.cfg)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/share/debug_malloc.cfg (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/share/debug_malloc.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -0,0 +1,5 @@
+# Settings used to debug memory allocation/deallocation issues.
+[pylithapp.petsc]
+malloc = True
+malloc_debug = True
+malloc_dump = True
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/plot_friction.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/plot_friction.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/plot_friction.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -1,7 +1,15 @@
#!/usr/bin/env python
-sim = "ratestate_stable"
+import sys
+if len(sys.argv) != 2:
+ raise ValueError("usage: plot_friction.py [weak | stable]")
+
+sim = sys.argv[1]
+
+if not sim in ['weak', 'stable']:
+ raise ValueError("Unknown sim '%s'." % sim)
+
# ======================================================================
import tables
import pylab
@@ -12,10 +20,10 @@
dt = 0.01
t = numpy.arange(0.0, 14.001, dt)
mu0 = 0.6
-if sim == "ratestate_stable":
+if sim == "stable":
a = 0.016
b = 0.012
-elif sim == "ratestate_weak":
+elif sim == "weak":
a = 0.008
b = 0.012
else:
@@ -64,7 +72,7 @@
# ----------------------------------------------------------------------
-h5 = tables.openFile("output/%s-fault.h5" % sim, "r")
+h5 = tables.openFile("output/ratestate_%s-fault.h5" % sim, "r")
time = h5.root.time[:].ravel()
slip = h5.root.vertex_fields.slip[:]
slipRate = h5.root.vertex_fields.slip_rate[:]
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/plasticity/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/plasticity/Makefile.am 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/plasticity/Makefile.am 2012-03-14 16:49:00 UTC (rev 19779)
@@ -18,6 +18,6 @@
SUBDIRS = \
threehex8 \
- threehex27
+ dynamic
# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/TestSlipWeakeningShearSliding.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/TestSlipWeakeningShearSliding.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/TestSlipWeakeningShearSliding.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -134,8 +134,8 @@
strikeDir = (0.0, -1.0)
normalDir = (-1.0, 0.0)
initialTraction = (0.0, -1.0e+6)
- staticCoefficient = 0.6
- dynamicCoefficient = 0.59
+ staticCoefficient = 0.61
+ dynamicCoefficient = 0.60
slipWeakeningParameter = 0.2
uy_l = 1.0
@@ -196,6 +196,7 @@
elif name == "previous_slip":
field = numpy.zeros( (nvertices, 1), dtype=numpy.float64)
+ field[:] = cumulativeSlip
else:
raise ValueError("Unknown fault field '%s'." % name)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/slipweakening_shear_sliding.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/slipweakening_shear_sliding.cfg 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/slipweakening_shear_sliding.cfg 2012-03-14 16:49:00 UTC (rev 19779)
@@ -136,7 +136,7 @@
friction.db_properties = spatialdata.spatialdb.UniformDB
friction.db_properties.label = Slip weakening
friction.db_properties.values = [static-coefficient,dynamic-coefficient,slip-weakening-parameter,cohesion]
-friction.db_properties.data = [0.6,0.59,0.2*m,0.0*Pa]
+friction.db_properties.data = [0.61,0.6,0.2*m,0.0*Pa]
# ----------------------------------------------------------------------
# PETSc
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPrager3D.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPrager3D.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPrager3D.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -36,7 +36,11 @@
pylith::materials::TestDruckerPrager3D::setUp(void)
{ // setUp
_material = new DruckerPrager3D();
- _matElastic = new DruckerPrager3D();
+ DruckerPrager3D* matTmp = new DruckerPrager3D();
+ CPPUNIT_ASSERT(matTmp);
+ matTmp->allowTensileYield(true);
+ _matElastic = matTmp;
+
_data = new DruckerPrager3DElasticData();
_dataElastic = new DruckerPrager3DElasticData();
setupNormalizer();
@@ -92,6 +96,18 @@
} // testUseElasticBehavior
// ----------------------------------------------------------------------
+// Test allowTensileYield()
+void
+pylith::materials::TestDruckerPrager3D::testAllowTensileYield(void)
+{ // testAllowTensileYield
+ DruckerPrager3D material;
+ CPPUNIT_ASSERT(!material._allowTensileYield);
+
+ material.allowTensileYield(true);
+ CPPUNIT_ASSERT(material._allowTensileYield);
+} // testAllowTensileYield
+
+// ----------------------------------------------------------------------
// Test usesHasStateVars()
void
pylith::materials::TestDruckerPrager3D::testHasStateVars(void)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPrager3D.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPrager3D.hh 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPrager3D.hh 2012-03-14 16:49:00 UTC (rev 19779)
@@ -59,6 +59,7 @@
// Need to test Drucker-Prager elastoplastic specific behavior.
CPPUNIT_TEST( testTimeStep );
CPPUNIT_TEST( testUseElasticBehavior );
+ CPPUNIT_TEST( testAllowTensileYield );
CPPUNIT_TEST( testHasStateVars );
CPPUNIT_TEST( test_calcStressElastic );
@@ -85,6 +86,9 @@
/// Test useElasticBehavior()
void testUseElasticBehavior(void);
+ /// Test allowTensileYield()
+ void testAllowTensileYield(void);
+
/// Test hasStateVars()
void testHasStateVars(void);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.cc 2012-03-14 16:49:00 UTC (rev 19779)
@@ -36,7 +36,11 @@
pylith::materials::TestDruckerPragerPlaneStrain::setUp(void)
{ // setUp
_material = new DruckerPragerPlaneStrain();
- _matElastic = new DruckerPragerPlaneStrain();
+ DruckerPragerPlaneStrain* matTmp = new DruckerPragerPlaneStrain();
+ CPPUNIT_ASSERT(matTmp);
+ matTmp->allowTensileYield(true);
+ _matElastic = matTmp;
+
_data = new DruckerPragerPlaneStrainElasticData();
_dataElastic = new DruckerPragerPlaneStrainElasticData();
setupNormalizer();
@@ -98,6 +102,18 @@
} // testUseElasticBehavior
// ----------------------------------------------------------------------
+// Test allowTensileYield()
+void
+pylith::materials::TestDruckerPragerPlaneStrain::testAllowTensileYield(void)
+{ // testAllowTensileYield
+ DruckerPragerPlaneStrain material;
+ CPPUNIT_ASSERT(!material._allowTensileYield);
+
+ material.allowTensileYield(true);
+ CPPUNIT_ASSERT(material._allowTensileYield);
+} // testAllowTensileYield
+
+// ----------------------------------------------------------------------
// Test usesHasStateVars()
void
pylith::materials::TestDruckerPragerPlaneStrain::testHasStateVars(void)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/materials/TestDruckerPragerPlaneStrain.hh 2012-03-14 16:49:00 UTC (rev 19779)
@@ -59,6 +59,7 @@
// Need to test Drucker-Prager elastoplastic specific behavior.
CPPUNIT_TEST( testTimeStep );
CPPUNIT_TEST( testUseElasticBehavior );
+ CPPUNIT_TEST( testAllowTensileYield );
CPPUNIT_TEST( testHasStateVars );
CPPUNIT_TEST( test_calcStressElastic );
@@ -85,6 +86,9 @@
/// Test useElasticBehavior()
void testUseElasticBehavior(void);
+ /// Test allowTensileYield()
+ void testAllowTensileYield(void);
+
/// Test hasStateVars()
void testHasStateVars(void);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPrager3D.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPrager3D.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPrager3D.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -46,6 +46,14 @@
return
+ def test_allowTensileYield(self):
+ """
+ Test allowTensileYield().
+ """
+ self.material.allowTensileYield(True)
+ return
+
+
def test_fitMohrCoulomb(self):
"""
Test useElasticBehavior().
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/materials/TestDruckerPragerPlaneStrain.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -46,6 +46,14 @@
return
+ def test_allowTensileYield(self):
+ """
+ Test allowTensileYield().
+ """
+ self.material.allowTensileYield(True)
+ return
+
+
def test_fitMohrCoulomb(self):
"""
Test useElasticBehavior().
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/Makefile.am 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/Makefile.am 2012-03-14 16:49:00 UTC (rev 19779)
@@ -35,6 +35,7 @@
TestOutputManagerMesh.py \
TestOutputManagerSubMesh.py \
TestOutputSolnSubset.py \
+ TestOutputSolnPoints.py \
TestDataWriterVTK.py \
TestDataWriterHDF5.py \
TestDataWriterHDF5Ext.py \
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/TestOutputSolnPoints.py (from rev 19778, short/3D/PyLith/trunk/unittests/pytests/meshio/TestOutputSolnPoints.py)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/TestOutputSolnPoints.py (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/TestOutputSolnPoints.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -0,0 +1,192 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/meshio/TestOutputSolnPoints.py
+
+## @brief Unit testing of Python OutputSolnPoints object.
+
+import unittest
+
+from pylith.meshio.OutputSolnPoints import OutputSolnPoints
+
+# ----------------------------------------------------------------------
+class TestOutputSolnPoints(unittest.TestCase):
+ """
+ Unit testing of Python OutputSolnPoints object.
+ """
+
+ def setUp(self):
+ from pylith.meshio.MeshIOAscii import MeshIOAscii
+ iohandler = MeshIOAscii()
+ filename = "data/twohex8.txt"
+
+ from spatialdata.units.Nondimensional import Nondimensional
+ normalizer = Nondimensional()
+ normalizer._configure()
+
+ from spatialdata.geocoords.CSCart import CSCart
+ iohandler.inventory.filename = filename
+ iohandler.inventory.coordsys = CSCart()
+ iohandler._configure()
+ mesh = iohandler.read(debug=False, interpolate=False)
+
+ from pylith.topology.SolutionFields import SolutionFields
+ fields = SolutionFields(mesh)
+
+ name = "disp(t)"
+ fields.add(name, "displacement")
+ fields.solutionName(name)
+ field = fields.get(name)
+ field.newSection(field.VERTICES_FIELD, mesh.dimension())
+ field.allocate()
+
+ self.mesh = mesh
+ self.fields = fields
+ self.normalizer = normalizer
+ return
+
+
+ def test_constructor(self):
+ """
+ Test constructor.
+ """
+ output = OutputSolnPoints()
+ output.inventory.writer._configure()
+ output.inventory.reader.inventory.filename = "data/points.txt"
+ output.inventory.reader._configure()
+ output._configure()
+ return
+
+
+ def test_preinitialize(self):
+ """
+ Test preinitialize().
+ """
+ output = OutputSolnPoints()
+ output.inventory.reader.inventory.filename = "data/points.txt"
+ output.inventory.reader._configure()
+ output._configure()
+ output.preinitialize()
+
+ self.assertEqual(output, output.dataProvider)
+ return
+
+
+ def test_verifyConfiguration(self):
+ """
+ Test verifyConfiguration().
+ """
+ output = OutputSolnPoints()
+ output.inventory.reader.inventory.filename = "data/points.txt"
+ output.inventory.reader._configure()
+ output._configure()
+ output.preinitialize()
+
+ output.vertexDataFields = ["displacement"]
+ output.verifyConfiguration(self.mesh)
+ return
+
+
+ def test_initialize(self):
+ """
+ Test initialize().
+ """
+ output = OutputSolnPoints()
+ output.inventory.reader.inventory.filename = "data/points.txt"
+ output.inventory.reader._configure()
+ output.inventory.writer.inventory.filename = "test.vtk"
+ output.inventory.writer._configure()
+ output._configure()
+
+ output.preinitialize()
+ output.initialize(self.mesh, self.normalizer)
+ return
+
+
+ def test_openclose(self):
+ """
+ Test open() and close().
+ """
+ output = OutputSolnPoints()
+ output.inventory.reader.inventory.filename = "data/points.txt"
+ output.inventory.reader._configure()
+ output.inventory.writer.inventory.filename = "test.vtk"
+ output.inventory.writer._configure()
+ output._configure()
+
+ output.preinitialize()
+ output.initialize(self.mesh, self.normalizer)
+
+ from pyre.units.time import s
+ output.open(totalTime=5.0*s, numTimeSteps=2)
+ output.close()
+ return
+
+
+ def test_writeInfo(self):
+ """
+ Test writeInfo().
+ """
+ output = OutputSolnPoints()
+ output.inventory.reader.inventory.filename = "data/points.txt"
+ output.inventory.reader._configure()
+ output.inventory.writer.inventory.filename = "output_sub.vtk"
+ output.inventory.writer._configure()
+ output._configure()
+
+ output.preinitialize()
+ output.initialize(self.mesh, self.normalizer)
+
+ output.open(totalTime=5.0, numTimeSteps=2)
+ output.writeInfo()
+ output.close()
+ return
+
+
+ def test_writeData(self):
+ """
+ Test writeData().
+ """
+ output = OutputSolnPoints()
+ output.inventory.reader.inventory.filename = "data/points.txt"
+ output.inventory.reader._configure()
+ output.inventory.writer.inventory.filename = "outputsub.vtk"
+ output.inventory.writer.inventory.timeFormat = "%3.1f"
+ output.inventory.writer._configure()
+ output.inventory.vertexDataFields = ["displacement"]
+ output._configure()
+
+ output.preinitialize()
+ output.initialize(self.mesh, self.normalizer)
+
+ output.open(totalTime=5.0, numTimeSteps=2)
+ output.writeData(2.0, self.fields)
+ output.close()
+ return
+
+
+ def test_factory(self):
+ """
+ Test factory method.
+ """
+ from pylith.meshio.OutputSolnPoints import output_manager
+ o = output_manager()
+ return
+
+
+# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/data/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/data/Makefile.am 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/data/Makefile.am 2012-03-14 16:49:00 UTC (rev 19779)
@@ -22,7 +22,8 @@
cube2_ascii.pset \
cube2.txt \
twohex8.exo \
- twohex8.txt
+ twohex8.txt \
+ points.txt
noinst_TMP = \
cube2_test.txt \
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/data/points.txt (from rev 19778, short/3D/PyLith/trunk/unittests/pytests/meshio/data/points.txt)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/data/points.txt (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/data/points.txt 2012-03-14 16:49:00 UTC (rev 19779)
@@ -0,0 +1,4 @@
+# Here are some points within the volume of the hex8 mesh.
+ 0.5 0.2 0.0
+-0.5 0.2 0.0
+ 0.5 -0.2 0.0
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/testmeshio.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/testmeshio.py 2012-03-14 16:44:02 UTC (rev 19778)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/testmeshio.py 2012-03-14 16:49:00 UTC (rev 19779)
@@ -98,6 +98,9 @@
from TestOutputSolnSubset import TestOutputSolnSubset
suite.addTest(unittest.makeSuite(TestOutputSolnSubset))
+ from TestOutputSolnPoints import TestOutputSolnPoints
+ suite.addTest(unittest.makeSuite(TestOutputSolnPoints))
+
from TestSingleOutput import TestSingleOutput
suite.addTest(unittest.makeSuite(TestSingleOutput))
More information about the CIG-COMMITS
mailing list