[cig-commits] r19381 - in short/3D/PyLith/branches/pylith-scecdynrup: . doc/releasenotes examples/3d/hex8 libsrc/pylith/faults libsrc/pylith/meshio libsrc/pylith/problems libsrc/pylith/topology pylith/tests tests/2d/frictionslide tests/3d/cyclicfriction tests_auto/2d/quad4 unittests/libtests/faults unittests/libtests/faults/data unittests/libtests/friction unittests/libtests/meshio unittests/pytests/meshio
brad at geodynamics.org
brad at geodynamics.org
Wed Jan 18 13:51:13 PST 2012
Author: brad
Date: 2012-01-18 13:51:13 -0800 (Wed, 18 Jan 2012)
New Revision: 19381
Added:
short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/announce_v1.6.2.txt
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/fieldsplit.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/output/
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/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.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/libsrc/pylith/faults/CohesiveTopology.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/DataWriterVTK.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/Formulation.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/MeshRefiner.cc
short/3D/PyLith/branches/pylith-scecdynrup/pylith/tests/Fault.py
short/3D/PyLith/branches/pylith-scecdynrup/setup.py
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/plot_friction.py
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/pylithapp.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/ratestate.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/geometry.jou
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/hex8_1000m.exo
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/hex8_1000m.jou
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/initial_tractions.spatialdb
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/pylithapp.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/friction_opening.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/slipweakening_opening.cfg
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynData.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynData.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/friction/TestFrictionModel.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestOutputSolnPoints.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/Makefile.am
Log:
Merge from trunk.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/README
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/README 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/README 2012-01-18 21:51:13 UTC (rev 19381)
@@ -83,7 +83,56 @@
- Fixed bug in writing tensor data for Xdmf files. Switched Tensor
to Tensor6 to account for symmetry.
+ - Fixed bug in writing HDF5 files in parallel when one processor
+ does not write any information (e.g., faults and boundary
+ conditions).
+ - Added dimensioning of time dataset in HDF5 files. The units are
+ now seconds rather than nondimensional time.
+
+ - Fixed memory allocation error (std::bad_alloc) when a processor
+ did not contain cells for a boundary condition or output. This bug
+ did not show up on all architectures.
+
+ - Increased robustness of spontaneous rupture (fault friction)
+ implementation to broaden the range of conditions it can
+ handle. The implementation now properly handles cases with fault
+ opening and cases with zero shear or normal tractions.
+
+
+* Internal changes
+
+ - Fault implementation
+
+ Several changes have been made to the fault implementation, but
+ none of these affect the user interface. The runtime performance
+ is nearly identical with improved accuracy for spontaneous rupture
+ (fault friction) simulations. These changes involved switching to
+ using tractions (non-integrated quantities) for the Lagrange
+ multipliers in the global coordinate system rather than integrated
+ quantities in the fault coordinate system. Additionally, initial
+ fault tractions are associated with the fault vertices and their
+ interpolation uses the finite-element basis functions.
+
+ - Distribution of mesh among processors
+
+ The data structures used to distribute the mesh among processors
+ have been improved. This reduces memory use and runtime for this
+ stage of the simulations.
+
+
+KNOWN ISSUES
+
+ The custom line search used with the PETSc nonlinear solver (SNES)
+ has difficulty handling some loading cases. In cases where the
+ direction of the line search tends to be nearly orthogonal to the
+ residual, the rate of convergence in the SNES iterations is
+ extremely slow. In other cases the nonlinear solver gets stuck in a
+ local minimum. We plan to improve the line search algorithm in a
+ future release in order to resolve this issue and improve the rate
+ of convergence in spontaneous rupture simulations.
+
+
----------------------------------------------------------------------
Version 1.6.1
----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-01-18 21:51:13 UTC (rev 19381)
@@ -10,15 +10,10 @@
* Module building
- + We want the _wrap.cxx files to work independent of PETSC_ARCH
- specific details. This was try until we included petscconf.h in
- pylithtypes.i to resolve single/double precision issues.
+ + We want the _wrap.cxx files to work independent of PyLith and
+ PETSc configure options (for example, with/without HDF5 and
+ single/double precision).
-* GPU utilization
-
- + Implicit elasticity finite-element integration
- + Explicit elasticity finite-element integration
-
* Green's functions
* 2-D materials
@@ -36,6 +31,15 @@
Need to redo Maxwell time calculation. Use ratio.
Create benchmark for test (compare against fk)
+* GPU utilization
+
+ Modeled on snes/examples/tutorials/ex52.
+
+ Refactor integration routine so that it uses batches rather than individual cells.
+
+ + Implicit elasticity finite-element integration
+ + Explicit elasticity finite-element integration
+
* Reimplement parameters to use PackedFields. [BRAD]
Fields
Modified: short/3D/PyLith/branches/pylith-scecdynrup/configure.ac
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/configure.ac 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/configure.ac 2012-01-18 21:51:13 UTC (rev 19381)
@@ -17,7 +17,7 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ(2.59)
-AC_INIT([PyLith], [1.6.2], [cig-short at geodynamics.org])
+AC_INIT([PyLith], [1.7.0], [cig-short at geodynamics.org])
AC_CONFIG_AUX_DIR([./aux-config])
AC_CONFIG_HEADER([portinfo])
AC_CONFIG_MACRO_DIR([m4])
@@ -325,6 +325,8 @@
tests/3d/Makefile
tests/3d/matprops/Makefile
tests/3d/slipdir/Makefile
+ tests/3d/cyclicfriction/Makefile
+ tests/3d/cyclicfriction/output/Makefile
tests/3d/plasticity/Makefile
tests/3d/plasticity/threehex27/Makefile
tests/3d/plasticity/threehex27/config/Makefile
Modified: short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/Makefile.am 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/Makefile.am 2012-01-18 21:51:13 UTC (rev 19381)
@@ -28,7 +28,9 @@
announce_v1.5.0.txt \
announce_v1.5.1.txt \
announce_v1.5.2.txt \
- announce_v1.6.0.txt
+ announce_v1.6.0.txt \
+ announce_v1.6.1.txt \
+ announce_v1.6.2.txt
# End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/announce_v1.6.2.txt (from rev 19380, short/3D/PyLith/trunk/doc/releasenotes/announce_v1.6.2.txt)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/announce_v1.6.2.txt (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/doc/releasenotes/announce_v1.6.2.txt 2012-01-18 21:51:13 UTC (rev 19381)
@@ -0,0 +1,74 @@
+Greetings,
+
+I am pleased to announce the release of PyLith 1.6.2, a finite-element
+code designed to solve dynamic elastic problems and quasi-static
+viscoelastic problems in tectonic deformation.
+
+This release fixes several bugs in PyLith v1.6.1. We strongly
+recommend all users of previous PyLith releases switch to this latest
+release, especially anyone running in parallel or using fault friction.
+
+You can download the source code and binaries from
+
+ http://geodynamics.org/cig/software/packages/short/pylith
+
+Detailed installation instructions for the binary packages are in the
+User Manual with detailed building instructions for a few platforms
+in the INSTALL file bundled with the PyLith Installer utility.
+
+
+RELEASE NOTES
+
+* Bug fixes
+
+ - Fixed bug in writing tensor data for Xdmf files. Switched Tensor
+ to Tensor6 to account for symmetry.
+
+ - Fixed bug in writing HDF5 files in parallel when one processor
+ does not write any information (e.g., faults and boundary
+ conditions).
+
+ - Added dimensioning of time dataset in HDF5 files. The units are
+ now seconds rather than nondimensional time.
+
+ - Fixed memory allocation error (std::bad_alloc) when a processor
+ did not contain cells for a boundary condition or output. This bug
+ did not show up on all architectures.
+
+ - Increased robustness of spontaneous rupture (fault friction)
+ implementation to broaden the range of conditions it can
+ handle. The implementation now properly handles cases with fault
+ opening and cases with zero shear or normal tractions.
+
+
+* Internal changes
+
+ - Fault implementation
+
+ Several changes have been made to the fault implementation, but
+ none of these affect the user interface. The runtime performance
+ is nearly identical with improved accuracy for spontaneous rupture
+ (fault friction) simulations. These changes involved switching to
+ using tractions (non-integrated quantities) for the Lagrange
+ multipliers in the global coordinate system rather than integrated
+ quantities in the fault coordinate system. Additionally, initial
+ fault tractions are associated with the fault vertices and their
+ interpolation uses the finite-element basis functions.
+
+ - Distribution of mesh among processors
+
+ The data structures used to distribute the mesh among processors
+ have been improved. This reduces memory use and runtime for this
+ stage of the simulations.
+
+
+KNOWN ISSUES
+
+ The custom line search used with the PETSc nonlinear solver (SNES)
+ has difficulty handling some loading cases. In cases where the
+ direction of the line search tends to be nearly orthogonal to the
+ residual, the rate of convergence in the SNES iterations is
+ extremely slow. In other cases the nonlinear solver gets stuck in a
+ local minimum. We plan to improve the line search algorithm in a
+ future release in order to resolve this issue and improve the rate
+ of convergence in spontaneous rupture simulations.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/pylithapp.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -30,6 +30,10 @@
# Change the default mesh reader to the CUBIT reader.
reader = pylith.meshio.MeshIOCubit
+# Optimize ordering of mesh cells and vertices using reverse
+# Cuthill-KcKee algorithm.
+reorder_mesh = True
+
[pylithapp.mesh_generator.reader]
# Set filename of mesh to import.
filename = mesh/box_hex8_1000m.exo
@@ -78,7 +82,7 @@
sub_pc_factor_shift_type = nonzero
# Convergence parameters.
-ksp_rtol = 1.0e-8
+ksp_rtol = 1.0e-10
ksp_atol = 1.0e-20
ksp_max_it = 100
ksp_gmres_restart = 50
@@ -89,10 +93,11 @@
ksp_converged_reason = true
# Nonlinear solver monitoring options.
-snes_rtol = 1.0e-8
-snes_atol = 1.0e-18
+snes_rtol = 1.0e-10
+snes_atol = 1.0e-9
snes_max_it = 100
snes_monitor = true
+snes_ls_monitor = true
snes_view = true
snes_converged_reason = true
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step12.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step12.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step12.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -88,7 +88,7 @@
db_initial = spatialdata.spatialdb.UniformDB
db_initial.label = Dirichlet BC on +x
db_initial.values = [displacement-x,displacement-y]
-db_initial.data = [-1.0*m,0.0*m]
+db_initial.data = [-0.5*m,0.0*m]
db_rate = spatialdata.spatialdb.UniformDB
db_rate.label = Dirichlet rate BC on +x
@@ -156,10 +156,6 @@
#friction_ksp_view = true
friction_ksp_converged_reason = true
-# Reduce convergence tolerances.
-ksp_rtol = 1.0e-12
-ksp_atol = 1.0e-15
-
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step13.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step13.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step13.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -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.5,0.2*m,0.0*Pa]
+friction.db_properties.data = [0.6,0.4,0.2*m,0.0*Pa]
# ----------------------------------------------------------------------
# PETSc settings
@@ -156,10 +156,6 @@
#friction_ksp_view = true
friction_ksp_converged_reason = true
-# Reduce convergence tolerances.
-ksp_rtol = 1.0e-12
-ksp_atol = 1.0e-15
-
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step14.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step14.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/3d/hex8/step14.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -118,6 +118,7 @@
[pylithapp.timedependent.interfaces.fault]
# The label corresponds to the name of the nodeset in CUBIT.
label = fault
+zero_tolerance = 1.0e-12
# Use the rate-and-state aging friction model.
friction = pylith.friction.RateStateAgeing
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/CohesiveTopology.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/CohesiveTopology.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -661,11 +661,16 @@
for (SieveSubMesh::label_sequence::iterator v_iter = fVertices2Begin;
v_iter != fVertices2EndNew;
++v_iter) {
+ assert(coordinates->getFiberDimension(*v_iter) ==
+ coordinates->getFiberDimension(vertexRenumber[*v_iter]));
coordinates->updatePoint(vertexRenumber[*v_iter],
coordinates->restrictPoint(*v_iter));
- if (constraintCell)
+ if (constraintCell) {
+ assert(coordinates->getFiberDimension(*v_iter) ==
+ coordinates->getFiberDimension(vertexLagrangeRenumber[*v_iter]));
coordinates->updatePoint(vertexLagrangeRenumber[*v_iter],
coordinates->restrictPoint(*v_iter));
+ } // if
} // for
if (debug)
coordinates->view("Coordinates with shadow vertices");
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -175,15 +175,15 @@
assert(0 != _fields);
assert(0 != _logger);
- // Initial fault tractions have been assembled, so they do not need
- // assembling across processors.
+ // Cohesive cells with conventional vertices N and P, and constraint
+ // vertex L make contributions to the assembled residual:
+ //
+ // DOF P: \int_{S_f^+} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p dS
+ // DOF N: -\int_{S_f^+} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p dS
+ // DOF L: \int_S_f \tensor{N}_p^T ( \tensor{R} \cdot \vec{d}
+ // -\tensor{N}_{n^+} \cdot \vec{u}_{n^+}
+ // +\tensor{N}_{n^-} \cdot \vec{u}_{n^-} dS
- FaultCohesiveLagrange::integrateResidual(residual, t, fields);
-
- // No contribution if no initial tractions are specified.
- if (0 == _dbInitialTract)
- return;
-
const int setupEvent = _logger->eventId("FaIR setup");
const int geometryEvent = _logger->eventId("FaIR geometry");
const int computeEvent = _logger->eventId("FaIR compute");
@@ -198,20 +198,31 @@
// Get sections associated with cohesive cells
scalar_array residualVertexN(spaceDim);
scalar_array residualVertexP(spaceDim);
+ scalar_array residualVertexL(spaceDim);
const ALE::Obj<RealSection>& residualSection = residual.section();
assert(!residualSection.isNull());
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
+
+ const ALE::Obj<RealSection>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+
+ scalar_array dispTpdtVertexN(spaceDim);
+ scalar_array dispTpdtVertexP(spaceDim);
+ scalar_array dispTpdtVertexL(spaceDim);
+
+ scalar_array initialTractionsVertex(spaceDim);
+ ALE::Obj<RealSection> initialTractionsSection;
+ if (_dbInitialTract) {
+ initialTractionsSection = _fields->get("initial traction").section();
+ assert(!initialTractionsSection.isNull());
+ } // if
+
const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
assert(!areaSection.isNull());
- const ALE::Obj<RealSection>& initialTractionsSection =
- _fields->get("initial traction").section();
- assert(!initialTractionsSection.isNull());
-
- const ALE::Obj<RealSection>& dispRelSection =
- _fields->get("relative disp").section();
- assert(!dispRelSection.isNull());
-
const ALE::Obj<RealSection>& orientationSection =
_fields->get("orientation").section();
assert(!orientationSection.isNull());
@@ -246,53 +257,89 @@
#endif
// Get initial tractions at fault vertex.
- assert(spaceDim == initialTractionsSection->getFiberDimension(v_fault));
- const PylithScalar* initialTractionsVertex =
- initialTractionsSection->restrictPoint(v_fault);
- assert(initialTractionsVertex);
+ if (_dbInitialTract) {
+ initialTractionsSection->restrictPoint(v_fault,
+ &initialTractionsVertex[0],
+ initialTractionsVertex.size());
+ } else {
+ initialTractionsVertex = 0.0;
+ } // if/else
- // Get relative dislplacement at fault vertex.
- assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
- const PylithScalar* dispRelVertex = dispRelSection->restrictPoint(v_fault);
- assert(dispRelVertex);
+ // Get orientation associated with fault vertex.
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+ const PylithScalar* orientationVertex =
+ orientationSection->restrictPoint(v_fault);
+ assert(orientationVertex);
// Get area associated with fault vertex.
assert(1 == areaSection->getFiberDimension(v_fault));
assert(areaSection->restrictPoint(v_fault));
const PylithScalar areaVertex = *areaSection->restrictPoint(v_fault);
- // Get orientation associated with fault vertex.
- assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
- const PylithScalar* orientationVertex =
- orientationSection->restrictPoint(v_fault);
- assert(orientationVertex);
+ // Get disp(t) at conventional vertices and Lagrange vertex.
+ assert(spaceDim == dispTSection->getFiberDimension(v_negative));
+ const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
+ assert(dispTVertexN);
+ assert(spaceDim == dispTSection->getFiberDimension(v_positive));
+ const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
+ assert(dispTVertexP);
+
+ assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
+ const PylithScalar* dispTVertexL = dispTSection->restrictPoint(v_lagrange);
+ assert(dispTVertexL);
+
+ // Get dispIncr(t->t+dt) at conventional vertices and Lagrange vertex.
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_negative));
+ const PylithScalar* dispTIncrVertexN =
+ dispTIncrSection->restrictPoint(v_negative);
+ assert(dispTIncrVertexN);
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_positive));
+ const PylithScalar* dispTIncrVertexP =
+ dispTIncrSection->restrictPoint(v_positive);
+ assert(dispTIncrVertexP);
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
+ const PylithScalar* dispTIncrVertexL =
+ dispTIncrSection->restrictPoint(v_lagrange);
+ assert(dispTIncrVertexL);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
#endif
- // Initial (external) tractions oppose (internal) tractions
- // associated with Lagrange multiplier, so these terms have the
- // opposite sign as the integration of the Lagrange multipliers in
- // FaultCohesiveLagrange.
+ // Compute current estimate of displacement at time t+dt using
+ // solution increment.
for (int iDim=0; iDim < spaceDim; ++iDim) {
- residualVertexP[iDim] = areaVertex * initialTractionsVertex[iDim];
+ dispTpdtVertexN[iDim] = dispTVertexN[iDim] + dispTIncrVertexN[iDim];
+ dispTpdtVertexP[iDim] = dispTVertexP[iDim] + dispTIncrVertexP[iDim];
+ dispTpdtVertexL[iDim] = dispTVertexL[iDim] + dispTIncrVertexL[iDim];
} // for
- residualVertexN = -residualVertexP;
-
- // Only apply initial tractions if there is no opening.
- // If there is opening, zero out initial tractions
+
+ // Compute slip (in fault coordinates system) from displacements.
PylithScalar slipNormal = 0.0;
+ PylithScalar tractionNormal = 0.0;
const int indexN = spaceDim - 1;
for (int jDim=0; jDim < spaceDim; ++jDim) {
- slipNormal += orientationVertex[indexN*spaceDim+jDim]*dispRelVertex[jDim];
+ slipNormal += orientationVertex[indexN*spaceDim+jDim] *
+ (dispTpdtVertexP[jDim] - dispTpdtVertexN[jDim]);
+ tractionNormal +=
+ orientationVertex[indexN*spaceDim+jDim] * dispTpdtVertexL[jDim];
} // for
+
+ residualVertexN = 0.0;
+ residualVertexL = 0.0;
+ if (slipNormal < _zeroTolerance) { // if no opening
+ // Initial (external) tractions oppose (internal) tractions
+ // associated with Lagrange multiplier.
+ residualVertexN = areaVertex * (dispTpdtVertexL - initialTractionsVertex);
- if (slipNormal > 0.0) {
- residualVertexN = 0.0;
- residualVertexP = 0.0;
- } // if
+ } else { // opening, normal traction should be zero
+ assert(fabs(tractionNormal) < _zeroTolerance);
+ } // if/else
+ residualVertexP = -residualVertexN;
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
@@ -308,11 +355,15 @@
residualSection->getFiberDimension(v_positive));
residualSection->updateAddPoint(v_positive, &residualVertexP[0]);
+ assert(residualVertexL.size() ==
+ residualSection->getFiberDimension(v_lagrange));
+ residualSection->updateAddPoint(v_lagrange, &residualVertexL[0]);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
- PetscLogFlops(numVertices*spaceDim*(2+spaceDim*2));
+ PetscLogFlops(numVertices*spaceDim*8);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
@@ -329,7 +380,7 @@
assert(0 != fields);
assert(0 != _fields);
- _updateVelRel(*fields);
+ _updateRelMotion(*fields);
const int spaceDim = _quadrature->spaceDim();
@@ -462,16 +513,19 @@
assert(0 != _fields);
assert(0 != _friction);
- _updateVelRel(*fields);
_sensitivitySetup(jacobian);
// Update time step in friction (can vary).
_friction->timeStep(_dt);
+ const PylithScalar dt = _dt;
const int spaceDim = _quadrature->spaceDim();
+ const int indexN = spaceDim - 1;
// Allocate arrays for vertex values
scalar_array tractionTpdtVertex(spaceDim);
+ scalar_array dTractionTpdtVertex(spaceDim);
+ scalar_array dDispRelVertex(spaceDim);
// Get sections
scalar_array slipVertex(spaceDim);
@@ -493,9 +547,9 @@
scalar_array dDispTIncrVertexN(spaceDim);
scalar_array dDispTIncrVertexP(spaceDim);
- const ALE::Obj<RealSection>& dispTIncrSection =
+ const ALE::Obj<RealSection>& dispIncrSection =
fields->get("dispIncr(t->t+dt)").section();
- assert(!dispTIncrSection.isNull());
+ assert(!dispIncrSection.isNull());
scalar_array dLagrangeTpdtVertex(spaceDim);
scalar_array dLagrangeTpdtVertexGlobal(spaceDim);
@@ -526,24 +580,36 @@
#if 0 // DEBUGGING
dispRelSection->view("BEFORE RELATIVE DISPLACEMENT");
- dispTIncrSection->view("BEFORE DISP INCR (t->t+dt)");
+ dispIncrSection->view("BEFORE DISP INCR (t->t+dt)");
#endif
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
+ const int v_negative = _cohesiveVertices[iVertex].negative;
+ const int v_positive = _cohesiveVertices[iVertex].positive;
- // Get relative displacement
- assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
- const PylithScalar* dispRelVertex = dispRelSection->restrictPoint(v_fault);
- assert(dispRelVertex);
+ // Get displacement values
+ assert(spaceDim == dispTSection->getFiberDimension(v_negative));
+ const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
+ assert(dispTVertexN);
- // Get relative velocity
- assert(spaceDim == velRelSection->getFiberDimension(v_fault));
- const PylithScalar* velRelVertex = velRelSection->restrictPoint(v_fault);
- assert(velRelVertex);
+ assert(spaceDim == dispTSection->getFiberDimension(v_positive));
+ const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
+ assert(dispTVertexP);
+ // Get displacement increment values.
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
+ const PylithScalar* dispIncrVertexN =
+ dispIncrSection->restrictPoint(v_negative);
+ assert(dispIncrVertexN);
+
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
+ const PylithScalar* dispIncrVertexP =
+ dispIncrSection->restrictPoint(v_positive);
+ assert(dispIncrVertexP);
+
// Get orientation
assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
const PylithScalar* orientationVertex =
@@ -554,11 +620,15 @@
const PylithScalar* lagrangeTVertex = dispTSection->restrictPoint(v_lagrange);
assert(lagrangeTVertex);
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_lagrange));
const PylithScalar* lagrangeTIncrVertex =
- dispTIncrSection->restrictPoint(v_lagrange);
+ dispIncrSection->restrictPoint(v_lagrange);
assert(lagrangeTIncrVertex);
+ // Step 1: Prevent nonphysical trial solutions. The product of the
+ // normal traction and normal slip must be nonnegative (forbid
+ // interpenetration with tension or opening with compression).
+
// Compute slip, slip rate, and Lagrange multiplier at time t+dt
// in fault coordinate system.
slipVertex = 0.0;
@@ -567,14 +637,60 @@
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
slipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
- dispRelVertex[jDim];
+ (dispTVertexP[jDim] + dispIncrVertexP[jDim]
+ - dispTVertexN[jDim] - dispIncrVertexN[jDim]);
slipRateVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
- velRelVertex[jDim];
+ (dispIncrVertexP[jDim] - dispIncrVertexN[jDim]) / dt;
tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
(lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
} // for
+ if (fabs(slipRateVertex[iDim]) < _zeroTolerance) {
+ slipRateVertex[iDim] = 0.0;
+ } // if
} // for
+ if (fabs(slipVertex[indexN]) < _zeroTolerance) {
+ slipVertex[indexN] = 0.0;
+ } // if
+ PylithScalar dSlipVertexNormal = 0.0;
+ PylithScalar dTractionTpdtVertexNormal = 0.0;
+ if (slipVertex[indexN]*tractionTpdtVertex[indexN] < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 1 CORRECTING NONPHYSICAL SLIP/TRACTIONS"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipVertex[indexN]
+ << ", tractionNormal: " << tractionTpdtVertex[indexN]
+ << std::endl;
+#endif
+ // Don't know what behavior is appropriate so set smaller of
+ // traction and slip to zero (should be appropriate if problem
+ // is nondimensionalized correctly).
+ if (fabs(slipVertex[indexN]) > fabs(tractionTpdtVertex[indexN])) {
+ // slip is bigger, so force normal traction back to zero
+ dTractionTpdtVertexNormal = -tractionTpdtVertex[indexN];
+ tractionTpdtVertex[indexN] = 0.0;
+ } else {
+ // traction is bigger, so force slip back to zero
+ dSlipVertexNormal = -slipVertex[indexN];
+ slipVertex[indexN] = 0.0;
+ } // if/else
+ } // if
+ if (slipVertex[indexN] < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 1 CORRECTING INTERPENETRATION"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipVertex[indexN]
+ << ", tractionNormal: " << tractionTpdtVertex[indexN]
+ << std::endl;
+#endif
+ dSlipVertexNormal = -slipVertex[indexN];
+ slipVertex[indexN] = 0.0;
+ } // if
+
+ // Step 2: Apply friction criterion to trial solution to get
+ // change in Lagrange multiplier (dLagrangeTpdtVertex) in fault
+ // coordinate system.
+
// Get friction properties and state variables.
_friction->retrievePropsStateVars(v_fault);
@@ -588,13 +704,18 @@
tractionTpdtVertex,
iterating);
- // Rotate traction back to global coordinate system.
+ // Rotate increment in traction back to global coordinate system.
dLagrangeTpdtVertexGlobal = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
dLagrangeTpdtVertexGlobal[iDim] +=
orientationVertex[jDim*spaceDim+iDim] * dLagrangeTpdtVertex[jDim];
} // for
+
+ // Add in potential contribution from adjusting Lagrange
+ // multiplier for fault normal DOF of trial solution in Step 1.
+ dLagrangeTpdtVertexGlobal[iDim] +=
+ orientationVertex[indexN*spaceDim+iDim] * dTractionTpdtVertexNormal;
} // for
#if 0 // debugging
@@ -622,11 +743,39 @@
std::cout << std::endl;
#endif
+ // Set change in Lagrange multiplier
assert(dLagrangeTpdtVertexGlobal.size() ==
dLagrangeTpdtSection->getFiberDimension(v_fault));
dLagrangeTpdtSection->updatePoint(v_fault, &dLagrangeTpdtVertexGlobal[0]);
+
+ // Update displacement in trial solution (if necessary) so that it
+ // conforms to physical constraints.
+ if (0.0 != dSlipVertexNormal) {
+ // Compute relative displacement from slip.
+ dDispRelVertex = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ dDispRelVertex[iDim] +=
+ orientationVertex[indexN*spaceDim+iDim] * dSlipVertexNormal;
+
+ dDispTIncrVertexN[iDim] = -0.5*dDispRelVertex[iDim];
+ dDispTIncrVertexP[iDim] = +0.5*dDispRelVertex[iDim];
+ } // for
+
+ // Update displacement field
+ assert(dDispTIncrVertexN.size() ==
+ dispIncrSection->getFiberDimension(v_negative));
+ dispIncrSection->updateAddPoint(v_negative, &dDispTIncrVertexN[0]);
+
+ assert(dDispTIncrVertexP.size() ==
+ dispIncrSection->getFiberDimension(v_positive));
+ dispIncrSection->updateAddPoint(v_positive, &dDispTIncrVertexP[0]);
+ } // if
+
} // for
+ // Step 3: Calculate change in displacement field corresponding to
+ // change in Lagrange multipliers imposed by friction criterion.
+
// Solve sensitivity problem for negative side of the fault.
bool negativeSide = true;
_sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
@@ -641,13 +790,14 @@
_sensitivitySolve();
_sensitivityUpdateSoln(negativeSide);
+ // Step 4: Update Lagrange multipliers and displacement fields based
+ // on changes imposed by friction criterion in Step 2 (change in
+ // Lagrange multipliers) and Step 3 (slip associated with change in
+ // Lagrange multipliers).
scalar_array dSlipVertex(spaceDim);
- scalar_array dDispRelVertex(spaceDim);
scalar_array dispRelVertex(spaceDim);
- // Update slip field based on solution of sensitivity problem and
- // increment in Lagrange multipliers.
const ALE::Obj<RealSection>& sensDispRelSection =
_fields->get("sensitivity relative disp").section();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
@@ -656,6 +806,10 @@
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+ // Get change in Lagrange multiplier computed from friction criterion.
+ dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
+ dLagrangeTpdtVertex.size());
+
// Get change in relative displacement from sensitivity solve.
assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
const PylithScalar* sensDispRelVertex =
@@ -672,7 +826,7 @@
orientationSection->restrictPoint(v_fault);
assert(orientationVertex);
- // Get displacements from disp(t) and dispIncr(t).
+ // Get displacement.
assert(spaceDim == dispTSection->getFiberDimension(v_negative));
const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
assert(dispTVertexN);
@@ -681,94 +835,122 @@
const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
assert(dispTVertexP);
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_negative));
- const PylithScalar* dispTIncrVertexN =
- dispTIncrSection->restrictPoint(v_negative);
- assert(dispTIncrVertexN);
+ // Get displacement increment (trial solution).
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
+ const PylithScalar* dispIncrVertexN =
+ dispIncrSection->restrictPoint(v_negative);
+ assert(dispIncrVertexN);
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_positive));
- const PylithScalar* dispTIncrVertexP =
- dispTIncrSection->restrictPoint(v_positive);
- assert(dispTIncrVertexP);
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
+ const PylithScalar* dispIncrVertexP =
+ dispIncrSection->restrictPoint(v_positive);
+ assert(dispIncrVertexP);
// Get Lagrange multiplier at time t
assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
const PylithScalar* lagrangeTVertex = dispTSection->restrictPoint(v_lagrange);
assert(lagrangeTVertex);
- // Get Lagrange multiplier increment at time t
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
+ // Get Lagrange multiplier increment (trial solution)
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_lagrange));
const PylithScalar* lagrangeTIncrVertex =
- dispTIncrSection->restrictPoint(v_lagrange);
+ dispIncrSection->restrictPoint(v_lagrange);
assert(lagrangeTIncrVertex);
- // Get change in Lagrange multiplier.
- dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
- dLagrangeTpdtVertex.size());
+ // Step 4a: Prevent nonphysical trial solutions. The product of the
+ // normal traction and normal slip must be nonnegative (forbid
+ // interpenetration with tension or opening with compression).
- // Only update slip, disp, etc if Lagrange multiplier is
- // changing. If Lagrange multiplier is the same, there is no slip
- // so disp, etc should be unchanged.
- PylithScalar dLagrangeMag = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim)
- dLagrangeMag += dLagrangeTpdtVertex[iDim]*dLagrangeTpdtVertex[iDim];
- if (0.0 == dLagrangeMag) {
- continue; // No change, so continue
- } // if
-
- // Compute slip and change in slip in fault coordinates.
+ // Compute slip, change in slip, and tractions in fault coordinates.
dSlipVertex = 0.0;
slipVertex = 0.0;
+ tractionTpdtVertex = 0.0;
+ dTractionTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
dSlipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
sensDispRelVertex[jDim];
slipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
(dispTVertexP[jDim] - dispTVertexN[jDim] +
- dispTIncrVertexP[jDim] - dispTIncrVertexN[jDim]);
+ dispIncrVertexP[jDim] - dispIncrVertexN[jDim]);
+ tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ (lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
+ dTractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ dLagrangeTpdtVertex[jDim];
} // for
} // for
+ if (fabs(slipVertex[indexN]) < _zeroTolerance) {
+ slipVertex[indexN] = 0.0;
+ } // if
+ if (fabs(dSlipVertex[indexN]) < _zeroTolerance) {
+ dSlipVertex[indexN] = 0.0;
+ } // if
- // Compute normal traction in fault coordinates.
- PylithScalar tractionNormal = 0.0;
- const int indexN = spaceDim - 1;
- for (int jDim=0; jDim < spaceDim; ++jDim) {
- tractionNormal += orientationVertex[indexN*spaceDim+jDim] *
- (lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim] +
- dLagrangeTpdtVertex[jDim]);
- } // for
+ if ((slipVertex[indexN] + dSlipVertex[indexN]) *
+ (tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN])
+ < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 4a CORRECTING NONPHYSICAL SLIP/TRACTIONS"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipVertex[indexN] + dSlipVertex[indexN]
+ << ", tractionNormal: " << tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN]
+ << std::endl;
+#endif
+ // Don't know what behavior is appropriate so set smaller of
+ // traction and slip to zero (should be appropriate if problem
+ // is nondimensionalized correctly).
+ if (fabs(slipVertex[indexN] + dSlipVertex[indexN]) >
+ fabs(tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN])) {
+ // slip is bigger, so force normal traction back to zero
+ dTractionTpdtVertex[indexN] = -tractionTpdtVertex[indexN];
+ } else {
+ // traction is bigger, so force slip back to zero
+ dSlipVertex[indexN] = -slipVertex[indexN];
+ } // if/else
- // Do not allow fault interpenetration and set fault opening to
- // zero if fault is under compression.
- if (tractionNormal < -_zeroTolerance ||
- slipVertex[indexN] + dSlipVertex[indexN] < 0.0) {
+ } // if
+ // Do not allow fault interpenetration.
+ if (slipVertex[indexN] + dSlipVertex[indexN] < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 4a CORRECTING INTERPENETATION"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipVertex[indexN] + dSlipVertex[indexN]
+ << std::endl;
+#endif
dSlipVertex[indexN] = -slipVertex[indexN];
} // if
- // Compute current estimate of slip.
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- slipVertex[iDim] = slipVertex[iDim] + dSlipVertex[iDim];
- } // for
+ // Update current estimate of slip from t to t+dt.
+ slipVertex += dSlipVertex;
- // Update relative displacement from slip.
+ // Compute relative displacement from slip.
dispRelVertex = 0.0;
dDispRelVertex = 0.0;
+ dLagrangeTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
dispRelVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
slipVertex[jDim];
dDispRelVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
dSlipVertex[jDim];
+ dLagrangeTpdtVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
+ dTractionTpdtVertex[jDim];
} // for
dDispTIncrVertexN[iDim] = -0.5*dDispRelVertex[iDim];
dDispTIncrVertexP[iDim] = +0.5*dDispRelVertex[iDim];
-
} // for
#if 0 // debugging
- std::cout << "dLagrangeTpdtVertex: ";
+ std::cout << "v_fault: " << v_fault;
+ std::cout << ", tractionTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
+ std::cout << " " << tractionTpdtVertex[iDim];
+ std::cout << ", dTractionTpdtVertex: ";
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ std::cout << " " << dTractionTpdtVertex[iDim];
+ std::cout << ", dLagrangeTpdtVertex: ";
+ for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << dLagrangeTpdtVertex[iDim];
std::cout << ", slipVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
@@ -786,27 +968,27 @@
assert(dispRelVertex.size() ==
dispRelSection->getFiberDimension(v_fault));
dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
-
+
// Update Lagrange multiplier increment.
assert(dLagrangeTpdtVertex.size() ==
- dispTIncrSection->getFiberDimension(v_lagrange));
- dispTIncrSection->updateAddPoint(v_lagrange, &dLagrangeTpdtVertex[0]);
+ dispIncrSection->getFiberDimension(v_lagrange));
+ dispIncrSection->updateAddPoint(v_lagrange, &dLagrangeTpdtVertex[0]);
// Update displacement field
assert(dDispTIncrVertexN.size() ==
- dispTIncrSection->getFiberDimension(v_negative));
- dispTIncrSection->updateAddPoint(v_negative, &dDispTIncrVertexN[0]);
+ dispIncrSection->getFiberDimension(v_negative));
+ dispIncrSection->updateAddPoint(v_negative, &dDispTIncrVertexN[0]);
assert(dDispTIncrVertexP.size() ==
- dispTIncrSection->getFiberDimension(v_positive));
- dispTIncrSection->updateAddPoint(v_positive, &dDispTIncrVertexP[0]);
+ dispIncrSection->getFiberDimension(v_positive));
+ dispIncrSection->updateAddPoint(v_positive, &dDispTIncrVertexP[0]);
} // for
#if 0 // DEBUGGING
//dLagrangeTpdtSection->view("AFTER dLagrange");
dispRelSection->view("AFTER RELATIVE DISPLACEMENT");
- dispTIncrSection->view("AFTER DISP INCR (t->t+dt)");
+ dispIncrSection->view("AFTER DISP INCR (t->t+dt)");
#endif
} // constrainSolnSpace
@@ -885,16 +1067,16 @@
const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- scalar_array dispTIncrVertexN(spaceDim);
- scalar_array dispTIncrVertexP(spaceDim);
+ scalar_array dispIncrVertexN(spaceDim);
+ scalar_array dispIncrVertexP(spaceDim);
scalar_array lagrangeTIncrVertex(spaceDim);
- const ALE::Obj<RealSection>& dispTIncrSection =
+ const ALE::Obj<RealSection>& dispIncrSection =
fields->get("dispIncr(t->t+dt)").section();
- assert(!dispTIncrSection.isNull());
+ assert(!dispIncrSection.isNull());
- const ALE::Obj<RealSection>& dispTIncrAdjSection = fields->get(
+ const ALE::Obj<RealSection>& dispIncrAdjSection = fields->get(
"dispIncr adjust").section();
- assert(!dispTIncrAdjSection.isNull());
+ assert(!dispIncrAdjSection.isNull());
const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
assert(!jacobianSection.isNull());
@@ -972,11 +1154,11 @@
assert(lagrangeTVertex);
// Get dispIncr(t) at cohesive cell's vertices.
- dispTIncrSection->restrictPoint(v_negative, &dispTIncrVertexN[0],
- dispTIncrVertexN.size());
- dispTIncrSection->restrictPoint(v_positive, &dispTIncrVertexP[0],
- dispTIncrVertexP.size());
- dispTIncrSection->restrictPoint(v_lagrange, &lagrangeTIncrVertex[0],
+ dispIncrSection->restrictPoint(v_negative, &dispIncrVertexN[0],
+ dispIncrVertexN.size());
+ dispIncrSection->restrictPoint(v_positive, &dispIncrVertexP[0],
+ dispIncrVertexP.size());
+ dispIncrSection->restrictPoint(v_lagrange, &lagrangeTIncrVertex[0],
lagrangeTIncrVertex.size());
// Get relative displacement at fault vertex.
@@ -1010,14 +1192,14 @@
assert(S > 0.0);
lagrangeTIncrVertex[iDim] = 1.0/S *
(-residualVertexL[iDim] +
- areaVertex * (dispTIncrVertexP[iDim] - dispTIncrVertexN[iDim]));
+ areaVertex * (dispIncrVertexP[iDim] - dispIncrVertexN[iDim]));
assert(jacobianVertexN[iDim] > 0.0);
- dispTIncrVertexN[iDim] =
+ dispIncrVertexN[iDim] =
+areaVertex / jacobianVertexN[iDim]*lagrangeTIncrVertex[iDim];
assert(jacobianVertexP[iDim] > 0.0);
- dispTIncrVertexP[iDim] =
+ dispIncrVertexP[iDim] =
-areaVertex / jacobianVertexP[iDim]*lagrangeTIncrVertex[iDim];
} // for
@@ -1061,12 +1243,12 @@
} // for
#if 0 // debugging
- std::cout << "dispTIncrP: ";
+ std::cout << "dispIncrP: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dispTIncrVertexP[iDim];
- std::cout << ", dispTIncrN: ";
+ std::cout << " " << dispIncrVertexP[iDim];
+ std::cout << ", dispIncrN: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dispTIncrVertexN[iDim];
+ std::cout << " " << dispIncrVertexN[iDim];
std::cout << ", slipVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << slipVertex[iDim];
@@ -1099,9 +1281,9 @@
assert(jacobianVertexP[iDim] > 0.0);
assert(jacobianVertexN[iDim] > 0.0);
- dispTIncrVertexN[iDim] +=
+ dispIncrVertexN[iDim] +=
areaVertex * dLagrangeTpdtVertexGlobal[iDim] / jacobianVertexN[iDim];
- dispTIncrVertexP[iDim] -=
+ dispIncrVertexP[iDim] -=
areaVertex * dLagrangeTpdtVertexGlobal[iDim] / jacobianVertexP[iDim];
// Set increment in relative displacement.
@@ -1122,13 +1304,13 @@
if (globalOrder->isLocal(v_lagrange)) {
// Adjust displacements to account for Lagrange multiplier values
// (assumed to be zero in preliminary solve).
- assert(dispTIncrVertexN.size() ==
- dispTIncrAdjSection->getFiberDimension(v_negative));
- dispTIncrAdjSection->updateAddPoint(v_negative, &dispTIncrVertexN[0]);
+ assert(dispIncrVertexN.size() ==
+ dispIncrAdjSection->getFiberDimension(v_negative));
+ dispIncrAdjSection->updateAddPoint(v_negative, &dispIncrVertexN[0]);
- assert(dispTIncrVertexP.size() ==
- dispTIncrAdjSection->getFiberDimension(v_positive));
- dispTIncrAdjSection->updateAddPoint(v_positive, &dispTIncrVertexP[0]);
+ assert(dispIncrVertexP.size() ==
+ dispIncrAdjSection->getFiberDimension(v_positive));
+ dispIncrAdjSection->updateAddPoint(v_positive, &dispIncrVertexP[0]);
} // if
// The Lagrange multiplier and relative displacement are NOT
@@ -1137,8 +1319,8 @@
// Set Lagrange multiplier value. Value from preliminary solve is
// bogus due to artificial diagonal entry in Jacobian of 1.0.
assert(lagrangeTIncrVertex.size() ==
- dispTIncrSection->getFiberDimension(v_lagrange));
- dispTIncrSection->updatePoint(v_lagrange, &lagrangeTIncrVertex[0]);
+ dispIncrSection->getFiberDimension(v_lagrange));
+ dispIncrSection->updatePoint(v_lagrange, &lagrangeTIncrVertex[0]);
// Update the relative displacement estimate based on adjustment
// to the Lagrange multiplier values.
@@ -1154,14 +1336,13 @@
9 + // updates
spaceDim*9));
-
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
#endif
#if 0 // DEBUGGING
//dLagrangeTpdtSection->view("AFTER dLagrange");
- //dispTIncrSection->view("AFTER DISP INCR (t->t+dt)");
+ //dispIncrSection->view("AFTER DISP INCR (t->t+dt)");
dispRelSection->view("AFTER RELATIVE DISPLACEMENT");
//velRelSection->view("AFTER RELATIVE VELOCITY");
#endif
@@ -1496,16 +1677,30 @@
} // _calcTractions
// ----------------------------------------------------------------------
-// Update slip rate associated with Lagrange vertex k corresponding
-// to diffential velocity between conventional vertices i and j.
+// Update relative displacement and velocity (slip and slip rate)
+// associated with Lagrange vertex k corresponding to diffential
+// velocity between conventional vertices i and j.
void
-pylith::faults::FaultCohesiveDyn::_updateVelRel(const topology::SolutionFields& fields)
-{ // _updateVelRel
+pylith::faults::FaultCohesiveDyn::_updateRelMotion(const topology::SolutionFields& fields)
+{ // _updateRelMotion
assert(0 != _fields);
const int spaceDim = _quadrature->spaceDim();
// Get section information
+ const ALE::Obj<RealSection>& dispTSection =
+ fields.get("disp(t)").section();
+ assert(!dispTSection.isNull());
+
+ const ALE::Obj<RealSection>& dispIncrSection =
+ fields.get("dispIncr(t->t+dt)").section();
+ assert(!dispIncrSection.isNull());
+
+ scalar_array dispRelVertex(spaceDim);
+ const ALE::Obj<RealSection>& dispRelSection =
+ _fields->get("relative disp").section();
+ assert(!dispRelSection.isNull());
+
const ALE::Obj<RealSection>& velocitySection =
fields.get("velocity(t)").section();
assert(!velocitySection.isNull());
@@ -1513,22 +1708,54 @@
scalar_array velRelVertex(spaceDim);
const ALE::Obj<RealSection>& velRelSection =
_fields->get("relative velocity").section();
+ assert(!velRelSection.isNull());
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
- // Get values
+ // Get displacement values
+ assert(spaceDim == dispTSection->getFiberDimension(v_negative));
+ const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
+ assert(dispTVertexN);
+
+ assert(spaceDim == dispTSection->getFiberDimension(v_positive));
+ const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
+ assert(dispTVertexP);
+
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
+ const PylithScalar* dispIncrVertexN =
+ dispIncrSection->restrictPoint(v_negative);
+ assert(dispIncrVertexN);
+
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
+ const PylithScalar* dispIncrVertexP =
+ dispIncrSection->restrictPoint(v_positive);
+ assert(dispIncrVertexP);
+
+ // Compute relative displacememt
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const PylithScalar value =
+ dispTVertexP[iDim] + dispIncrVertexP[iDim]
+ - dispTVertexN[iDim] - dispIncrVertexN[iDim];
+ dispRelVertex[iDim] = fabs(value) > _zeroTolerance ? value : 0.0;
+ } // for
+
+ // Update relative displacement field.
+ assert(dispRelVertex.size() ==
+ dispRelSection->getFiberDimension(v_fault));
+ dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
+
+ // Get velocity values
+ assert(spaceDim == velocitySection->getFiberDimension(v_negative));
const PylithScalar* velocityVertexN = velocitySection->restrictPoint(v_negative);
- assert(0 != velocityVertexN);
- assert(spaceDim == velocitySection->getFiberDimension(v_negative));
+ assert(velocityVertexN);
+ assert(spaceDim == velocitySection->getFiberDimension(v_positive));
const PylithScalar* velocityVertexP = velocitySection->restrictPoint(v_positive);
- assert(0 != velocityVertexP);
- assert(spaceDim == velocitySection->getFiberDimension(v_positive));
+ assert(velocityVertexP);
// Compute relative velocity
for (int iDim=0; iDim < spaceDim; ++iDim) {
@@ -1536,14 +1763,14 @@
velRelVertex[iDim] = fabs(value) > _zeroTolerance ? value : 0.0;
} // for
- // Update slip rate field.
+ // Update relative velocity field.
assert(velRelVertex.size() ==
velRelSection->getFiberDimension(v_fault));
velRelSection->updatePoint(v_fault, &velRelVertex[0]);
} // for
PetscLogFlops(numVertices*spaceDim*spaceDim*4);
-} // _updateVelRel
+} // _updateRelMotion
// ----------------------------------------------------------------------
// Setup sensitivity problem to compute change in slip given change in Lagrange multipliers.
@@ -1613,7 +1840,7 @@
int maxIters = 0;
err = KSPGetTolerances(_ksp, &rtol, &atol, &dtol, &maxIters);
CHECK_PETSC_ERROR(err);
- rtol = 1.0e-2*_zeroTolerance;
+ rtol = 1.0e-3*_zeroTolerance;
atol = 1.0e-5*_zeroTolerance;
err = KSPSetTolerances(_ksp, rtol, atol, dtol, maxIters);
CHECK_PETSC_ERROR(err);
@@ -1934,16 +2161,16 @@
{ // _constrainSolnSpace1D
assert(0 != dLagrangeTpdt);
- if (tractionTpdt[0] < 0) {
- // if compression, then no changes to solution
- } else {
- // if tension, then traction is zero.
-
- const PylithScalar dlp = -tractionTpdt[0];
- (*dLagrangeTpdt)[0] = dlp;
- } // else
-
- PetscLogFlops(2);
+ if (fabs(slip[0]) < _zeroTolerance) {
+ // if compression, then no changes to solution
+ } else {
+ // if tension, then traction is zero.
+
+ const PylithScalar dlp = -tractionTpdt[0];
+ (*dLagrangeTpdt)[0] = dlp;
+ } // else
+
+ PetscLogFlops(2);
} // _constrainSolnSpace1D
// ----------------------------------------------------------------------
@@ -1964,7 +2191,7 @@
const PylithScalar tractionShearMag = fabs(tractionTpdt[0]);
#if !defined(NO_FAULT_OPENING)
- if (tractionNormal < 0 && 0.0 == slip[1]) {
+ if (fabs(slip[1]) < _zeroTolerance && tractionNormal < -_zeroTolerance) {
#endif
// if in compression and no opening
const PylithScalar frictionStress = _friction->calcFriction(slipMag, slipRateMag,
@@ -1987,7 +2214,9 @@
} else {
// friction exceeds value necessary to stick
// no changes to solution
- assert(0.0 == slipRateMag);
+ if (iterating) {
+ assert(0.0 == slipRateMag);
+ } // if
} // if/else
#if !defined(NO_FAULT_OPENING)
} else {
@@ -2022,7 +2251,7 @@
tractionTpdt[1] * tractionTpdt[1]);
#if !defined(NO_FAULT_OPENING)
- if (tractionNormal < 0.0 && 0.0 == slip[2]) {
+ if (fabs(slip[2]) < _zeroTolerance && tractionNormal < -_zeroTolerance) {
#endif
// if in compression and no opening
const PylithScalar frictionStress =
@@ -2051,7 +2280,9 @@
} else {
// else friction exceeds value necessary, so stick
// no changes to solution
- assert(0.0 == slipRateMag);
+ if (iterating) {
+ assert(0.0 == slipRateMag);
+ } // if
} // if/else
#if !defined (NO_FAULT_OPENING)
} else {
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-01-18 21:51:13 UTC (rev 19381)
@@ -152,13 +152,13 @@
void _calcTractions(topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& solution);
- /** Update relative velocity associated with Lagrange vertex k
- * corresponding to diffential velocity between conventional
- * vertices i and j.
+ /** Update relative displacement and velocity associated with
+ * Lagrange vertex k corresponding to diffential velocity between
+ * conventional vertices i and j.
*
* @param fields Solution fields.
*/
- void _updateVelRel(const topology::SolutionFields& fields);
+ void _updateRelMotion(const topology::SolutionFields& fields);
/** Setup sensitivity problem to compute change in slip given change
* in Lagrange multipliers.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/DataWriterVTK.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/DataWriterVTK.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -198,12 +198,12 @@
assert(!section.isNull());
assert(!sieveMesh->getLabelStratum(labelName, 0).isNull());
- const int localFiberDim =
+ int fiberDimLocal =
(sieveMesh->getLabelStratum(labelName, 0)->size() > 0) ?
section->getFiberDimension(*sieveMesh->getLabelStratum(labelName, 0)->begin()) : 0;
int fiberDim = 0;
- MPI_Allreduce((void *) &localFiberDim, (void *) &fiberDim, 1,
- MPI_INT, MPI_MAX, field.mesh().comm());
+ MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
+ field.mesh().comm());
assert(fiberDim > 0);
const int enforceDim =
(field.vectorFieldType() != topology::FieldBase::VECTOR) ? fiberDim : 3;
@@ -266,12 +266,12 @@
const ALE::Obj<RealSection>& section = field.section();
assert(!section.isNull());
- const int localFiberDim =
+ int fiberDimLocal =
(sieveMesh->getLabelStratum(labelName, depth)->size() > 0) ?
section->getFiberDimension(*sieveMesh->getLabelStratum(labelName, depth)->begin()) : 0;
int fiberDim = 0;
- MPI_Allreduce((void *) &localFiberDim, (void *) &fiberDim, 1,
- MPI_INT, MPI_MAX, field.mesh().comm());
+ MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
+ field.mesh().comm());
assert(fiberDim > 0);
const int enforceDim =
(field.vectorFieldType() != topology::FieldBase::VECTOR) ? fiberDim : 3;
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/meshio/OutputSolnPoints.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -52,15 +52,17 @@
{ // deallocate
OutputManager<topology::Mesh, topology::Field<topology::Mesh> >::deallocate();
-#if 0 // :MATT: Need DM wrapper for mesh
if (_interpolator) {
assert(_mesh);
const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
- PetscErrorCode err =
- DMMeshInterpolationDestroy(sieveMesh, &_interpolator);
- CHECK_PETSC_ERROR(err);
+ DM dm;
+ PetscErrorCode err = 0;
+
+ err = DMMeshCreate(sieveMesh->comm(), &dm);CHECK_PETSC_ERROR(err);
+ err = DMMeshSetMesh(dm, sieveMesh);CHECK_PETSC_ERROR(err);
+ err = DMMeshInterpolationDestroy(dm, &_interpolator);CHECK_PETSC_ERROR(err);
+ err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
} // if
-#endif
_mesh = 0; // :TODO: Use shared pointer
delete _pointsMesh; _pointsMesh = 0;
@@ -92,7 +94,7 @@
// Create mesh without cells for points.
const int meshDim = 0;
delete _pointsMesh; _pointsMesh = new topology::Mesh(meshDim);
- _pointsMesh->createSieveMesh(0);
+ _pointsMesh->createSieveMesh(meshDim);
assert(_pointsMesh);
scalar_array pointsArray(points, numPoints*spaceDim);
@@ -108,26 +110,25 @@
interpolate);
_pointsMesh->coordsys(_mesh->coordsys());
-#if 0 // :MATT: Need DM wrapper for mesh
// Setup interpolator object
+ DM dm;
PetscErrorCode err = 0;
- err = DMMeshInterpolationCreate(_mesh->sieveMesh(), &_interpolator);
- CHECK_PETSC_ERROR(err);
+ assert(!_mesh->sieveMesh().isNull());
+ err = DMMeshCreate(_mesh->sieveMesh()->comm(), &dm);CHECK_PETSC_ERROR(err);
+ err = DMMeshSetMesh(dm, _mesh->sieveMesh());CHECK_PETSC_ERROR(err);
+ err = DMMeshInterpolationCreate(dm, &_interpolator);CHECK_PETSC_ERROR(err);
- const spatialdata::geocoords::CoordSys* cs = _pointsMesh().coordsys();
- assert(0 != cs);
+ const spatialdata::geocoords::CoordSys* cs = _pointsMesh->coordsys();
+ assert(cs);
assert(cs->spaceDim() == spaceDim);
- err = DMMeshInterpolationSetDim(_mesh->sieveMesh(), spaceDim,
- _interpolator); CHECK_PETSC_ERROR(err);
+ err = DMMeshInterpolationSetDim(dm, spaceDim, _interpolator);CHECK_PETSC_ERROR(err);
+ err = DMMeshInterpolationAddPoints(dm, numPoints, (PetscReal*) points, _interpolator);CHECK_PETSC_ERROR(err);
+ err = DMMeshInterpolationSetUp(dm, _interpolator);CHECK_PETSC_ERROR(err);
+ err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
+ CHECK_PETSC_ERROR(err);
- err = DMMeshInterpolationAddPoints(_mesh->sieveMesh(), numPoints, points,
- _interpolator); CHECK_PETSC_ERROR(err);
-
- err = DMMeshInterpolationSetUp(_mesh->sieveMesh(), _interpolator);
- CHECK_PETSC_ERROR(err);
-#endif
} // createPointsMesh
@@ -176,16 +177,22 @@
const ALE::Obj<SieveMesh>& sieveMesh = _mesh->sieveMesh();
assert(!sieveMesh.isNull());
-#if 0 // :MATT: Need DM wrapper for mesh
+ DM dm;
+ SectionReal section;
PetscErrorCode err = 0;
- err = DMMeshInterpolationSetDof(sieveMesh, fiberDim, _interpolator);
- CHECK_PETSC_ERROR(err);
-
- err = DMMeshInterpolationEvaluate(sieveMesh, field.section(), fieldInterpVec,
- _interpolator); CHECK_PETSC_ERROR(err);
-#endif
+ err = DMMeshCreate(sieveMesh->comm(), &dm);CHECK_PETSC_ERROR(err);
+ err = DMMeshSetMesh(dm, sieveMesh);CHECK_PETSC_ERROR(err);
+ err = DMMeshInterpolationSetDof(dm, fiberDim, _interpolator);CHECK_PETSC_ERROR(err);
+ err = SectionRealCreate(sieveMesh->comm(), §ion);CHECK_PETSC_ERROR(err);
+ err = SectionRealSetSection(section, field.section());CHECK_PETSC_ERROR(err);
+ err = SectionRealSetBundle(section, sieveMesh);CHECK_PETSC_ERROR(err);
+
+ err = DMMeshInterpolationEvaluate(dm, section, fieldInterpVec, _interpolator);CHECK_PETSC_ERROR(err);
+ err = SectionRealDestroy(§ion);CHECK_PETSC_ERROR(err);
+ err = DMDestroy(&dm);CHECK_PETSC_ERROR(err);
+
OutputManager<topology::Mesh, topology::Field<topology::Mesh> >::appendVertexField(t, fieldInterp, mesh);
} // appendVertexField
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -425,5 +425,52 @@
solution += adjust;
} // adjustSolnLumped
+#include "pylith/meshio/DataWriterHDF5.hh"
+// ----------------------------------------------------------------------
+void
+pylith::problems::Formulation::printState(PetscVec* solutionVec,
+ PetscVec* residualVec,
+ PetscVec* solution0Vec,
+ PetscVec* searchDirVec)
+{ // printState
+ assert(solutionVec);
+ assert(residualVec);
+ assert(searchDirVec);
+ meshio::DataWriterHDF5<topology::Mesh,topology::Field<topology::Mesh> > writer;
+
+ const topology::Mesh& mesh = _fields->mesh();
+
+ writer.filename("state.h5");
+ const int numTimeSteps = 1;
+ writer.open(mesh, numTimeSteps);
+
+
+ topology::Field<topology::Mesh>& solution = _fields->solution();
+ solution.scatterVectorToSection(*solutionVec);
+ writer.writeVertexField(0.0, solution, mesh);
+ solution.view("DIVERGED_SOLUTION");
+ const char* label = solution.label();
+
+ solution.label("solution_0");
+ solution.scatterVectorToSection(*solution0Vec);
+ writer.writeVertexField(0.0, solution, mesh);
+ solution.view("DIVERGED_SOLUTION0");
+ solution.label(label);
+
+ topology::Field<topology::Mesh>& residual = _fields->get("residual");
+ residual.scatterVectorToSection(*residualVec);
+ writer.writeVertexField(0.0, residual, mesh);
+ residual.view("DIVERGED_RESIDUAL");
+
+ residual.label("search_dir");
+ residual.scatterVectorToSection(*searchDirVec);
+ writer.writeVertexField(0.0, residual, mesh);
+ residual.view("DIVERGED_SEARCHDIR");
+
+ writer.close();
+} // printState
+
+
+
// End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.hh 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Formulation.hh 2012-01-18 21:51:13 UTC (rev 19381)
@@ -186,6 +186,12 @@
virtual
void calcRateFields(void) = 0;
+ /// Write state of system
+ void printState(PetscVec* solutionVec,
+ PetscVec* residualVec,
+ PetscVec* solution0Vec,
+ PetscVec* searchDirVec);
+
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -362,6 +362,10 @@
ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
}
*flag = PETSC_FALSE; // DIVERGED_LINE_SEARCH
+ assert(lsctx);
+ Formulation* formulation = (Formulation*) lsctx;
+ assert(formulation);
+ formulation->printState(&w, &g, &x, &y);
break;
}
t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
@@ -380,11 +384,11 @@
// range. Necessary due to underflow and overflow of a, b, c,
// and d.
#if 0
- lambdatemp = (-b + sqrt(d))/(3.0*a);
+ lambdatemp = (-b + PetscSqrtScalar(d))/(3.0*a);
#else
- if ((-b + sqrt(d) > 0.0 && a > 0.0) ||
- (-b + sqrt(d) < 0.0 && a < 0.0)) {
- lambdatemp = (-b + sqrt(d))/(3.0*a);
+ if ((-b + PetscSqrtScalar(d) > 0.0 && a > 0.0) ||
+ (-b + PetscSqrtScalar(d) < 0.0 && a < 0.0)) {
+ lambdatemp = (-b + PetscSqrtScalar(d))/(3.0*a);
} else {
lambdatemp = 0.05*lambda;
} // else
@@ -454,9 +458,9 @@
// ======================================================================
// Code to constrain solution space.
- assert(0 != lsctx);
+ assert(lsctx);
Formulation* formulation = (Formulation*) lsctx;
- assert(0 != formulation);
+ assert(formulation);
formulation->constrainSolnSpace(&w);
ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/MeshRefiner.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/MeshRefiner.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/topology/MeshRefiner.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -799,6 +799,12 @@
// We have to do flexible assembly since we add the new vertices separately
newSendOverlap->assemble();
newRecvOverlap->assemble();
+ if (mesh->debug()) {
+ sendOverlap->view("OLD SEND OVERLAP");
+ recvOverlap->view("OLD RECV OVERLAP");
+ newSendOverlap->view("NEW SEND OVERLAP");
+ newRecvOverlap->view("NEW RECV OVERLAP");
+ }
} // _calcNewOverlap
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/tests/Fault.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/tests/Fault.py 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/tests/Fault.py 2012-01-18 21:51:13 UTC (rev 19381)
@@ -58,7 +58,7 @@
testcase.assertEqual(dimE, dim)
scale = 1.0
- if name == "traction_change":
+ if name == "traction_change" or name == "traction":
scale *= normalizer.pressureScale().value
for i in xrange(dim):
Modified: short/3D/PyLith/branches/pylith-scecdynrup/setup.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/setup.py 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/setup.py 2012-01-18 21:51:13 UTC (rev 19381)
@@ -22,7 +22,7 @@
setup(
name = 'PyLith',
- version = '1.6.2',
+ version = '1.7.0',
zip_safe = False,
packages = find_packages(),
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-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/plot_friction.py 2012-01-18 21:51:13 UTC (rev 19381)
@@ -1,6 +1,6 @@
#!/usr/bin/env python
-sim = "ratestate_stable"
+sim = "ratestate_weak"
# ======================================================================
import tables
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/pylithapp.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/pylithapp.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -78,7 +78,7 @@
[pylithapp.timedependent.interfaces]
fault = pylith.faults.FaultCohesiveDyn
-fault.zero_tolerance = 1.0e-14
+fault.zero_tolerance = 1.0e-12
[pylithapp.timedependent.interfaces.fault]
id = 100
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/ratestate.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/ratestate.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/ratestate.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -46,6 +46,7 @@
db_initial_tractions.values = [traction-shear,traction-normal]
db_initial_tractions.data = [0.0*MPa, -10.0*MPa]
+zero_tolerance = 1.0e-14
# ----------------------------------------------------------------------
# PETSc
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/Makefile.am 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/Makefile.am 2012-01-18 21:51:13 UTC (rev 19381)
@@ -17,9 +17,16 @@
#
dist_noinst_DATA = \
+ geometry.jou \
+ hex8_1000m.jou \
+ hex8_1000m.exo \
pylithapp.cfg \
- matprops.spatialdb \
+ fieldsplit.cfg \
+ initial_tractions.spatialdb \
cycle.timedb
+SUBDIRS = \
+ output
+
# End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/fieldsplit.cfg (from rev 19380, short/3D/PyLith/trunk/tests/3d/cyclicfriction/fieldsplit.cfg)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/fieldsplit.cfg (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/fieldsplit.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -0,0 +1,18 @@
+[pylithapp.timedependent.formulation]
+split_fields = True
+matrix_type = aij
+use_custom_constraint_pc = True
+
+[pylithapp.petsc]
+ksp_gmres_restart = 500
+fs_pc_type = fieldsplit
+fs_pc_fieldsplit_real_diagonal =
+fs_pc_fieldsplit_type = multiplicative
+fs_fieldsplit_0_pc_type = ml
+fs_fieldsplit_1_pc_type = ml
+fs_fieldsplit_2_pc_type = ml
+fs_fieldsplit_3_pc_type = ml
+fs_fieldsplit_0_ksp_type = preonly
+fs_fieldsplit_1_ksp_type = preonly
+fs_fieldsplit_2_ksp_type = preonly
+fs_fieldsplit_3_ksp_type = preonly
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/geometry.jou
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/geometry.jou 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/geometry.jou 2012-01-18 21:51:13 UTC (rev 19381)
@@ -25,6 +25,7 @@
# Create interface surfaces
# ----------------------------------------------------------------------
create planar surface with plane xplane offset 0
+rotate surface 7 angle -30 about z include_merged
create sheet extended from surface 7 intersecting volume 1
delete body 2
surface 8 name "fault_surface"
@@ -52,3 +53,4 @@
imprint all with volume all
merge all
+
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/hex8_1000m.exo
===================================================================
(Binary files differ)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/hex8_1000m.jou
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/hex8_1000m.jou 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/hex8_1000m.jou 2012-01-18 21:51:13 UTC (rev 19381)
@@ -6,6 +6,7 @@
# ----------------------------------------------------------------------
# Set discretization size
# ----------------------------------------------------------------------
+surface 10 17 16 12 scheme pave
volume all size 1000
# ----------------------------------------------------------------------
@@ -101,3 +102,4 @@
# Export exodus file
# ----------------------------------------------------------------------
export mesh "hex8_1000m.exo" dimension 3 overwrite
+
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/initial_tractions.spatialdb
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/initial_tractions.spatialdb 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/initial_tractions.spatialdb 2012-01-18 21:51:13 UTC (rev 19381)
@@ -12,4 +12,4 @@
}
}
0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 -6.0 0.0 0.0 -60.0
+0.0 0.0 -6.0 0.0 0.0 -25.0
Copied: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/output (from rev 19380, short/3D/PyLith/trunk/tests/3d/cyclicfriction/output)
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/pylithapp.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/pylithapp.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -38,6 +38,7 @@
[pylithapp.timedependent.normalizer]
relaxation_time = 1.0*hour
+length_scale = 1.0*km
[pylithapp.timedependent.implicit]
solver = pylith.problems.SolverNonlinear
@@ -73,7 +74,7 @@
db_change = spatialdata.spatialdb.UniformDB
db_change.label = Amplitude of Dirichlet BC on +x
db_change.values = [displacement-x,displacement-y,displacement-z,change-start-time]
-db_change.data = [-0.2*m,0.5*m,0.0*m,0.0*hour]
+db_change.data = [-0.2*m,-0.4*m,0.0*m,0.0*hour]
th_change = spatialdata.spatialdb.TimeHistory
th_change.label = Time history of tidal load
@@ -88,7 +89,7 @@
db_change = spatialdata.spatialdb.UniformDB
db_change.label = Amplitude of Dirichlet BC on -x
db_change.values = [displacement-x,displacement-y,displacement-z,change-start-time]
-db_change.data = [+0.2*m,-0.5*m,0.0*m,0.0*hour]
+db_change.data = [+0.2*m,+0.4*m,0.0*m,0.0*hour]
th_change = spatialdata.spatialdb.TimeHistory
th_change.label = Time history of tidal load
@@ -113,6 +114,7 @@
[pylithapp.timedependent.interfaces.fault]
label = fault
+zero_tolerance = 1.0e-8
friction = pylith.friction.SlipWeakening
friction.label = Slip weakening
@@ -170,8 +172,8 @@
# Convergence parameters.
ksp_rtol = 1.0e-8
ksp_atol = 1.0e-10
-ksp_max_it = 200
-ksp_gmres_restart = 70
+ksp_max_it = 100
+ksp_gmres_restart = 50
# Linear solver monitoring options.
ksp_monitor = true
@@ -180,13 +182,26 @@
# Nonlinear solver monitoring options.
snes_rtol = 1.0e-8
-snes_atol = 1.0e-8
-snes_max_it = 200
+snes_atol = 1.0e-7
+snes_max_it = 100
snes_monitor = true
snes_ls_monitor = true
#snes_view = true
snes_converged_reason = true
-#snes_monitor_solution_update = true
+snes_monitor_solution_update = true
+snes_monitor_residual = true
#info =
+# Friction sensitivity solve used to compute the increment in slip
+# associated with changes in the Lagrange multiplier imposed by the
+# fault constitutive model.
+friction_pc_type = asm
+friction_sub_pc_factor_shift_type = nonzero
+friction_ksp_max_it = 200
+friction_ksp_gmres_restart = 50
+# Uncomment to view details of friction sensitivity solve.
+#friction_ksp_monitor = true
+#friction_ksp_view = true
+friction_ksp_converged_reason = true
+
log_summary = true
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/friction_opening.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/friction_opening.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/friction_opening.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -124,6 +124,7 @@
sub_pc_factor_shift_type = nonzero
ksp_rtol = 1.0e-8
+ksp_atol = 1.0e-12
ksp_max_it = 100
ksp_gmres_restart = 50
#ksp_monitor = true
@@ -131,8 +132,10 @@
#ksp_converged_reason = true
snes_rtol = 1.0e-8
+snes_atol = 1.0e-9
snes_max_it = 100
#snes_monitor = true
+#snes_ls_monitor = true
#snes_view = true
#snes_converged_reason = true
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/slipweakening_opening.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/slipweakening_opening.cfg 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests_auto/2d/quad4/slipweakening_opening.cfg 2012-01-18 21:51:13 UTC (rev 19381)
@@ -123,8 +123,12 @@
sub_pc_factor_shift_type = nonzero
ksp_rtol = 1.0e-8
+ksp_atol = 1.0e-12
ksp_max_it = 100
ksp_gmres_restart = 50
+
+snes_rtol = 1.0e-8
+snes_atol = 1.0e-9
snes_max_it = 200
#ksp_monitor = true
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -277,23 +277,37 @@
fault.vertexField("slip").section();
CPPUNIT_ASSERT(!slipSection.isNull());
- const PylithScalar valE = 0.0; // slip should be zero
+ //slipSection->view("SLIP"); // DEBUGGING
+
+ // Get expected values
+ CPPUNIT_ASSERT(_data->slipStickE);
+ const PylithScalar* valsE = _data->slipStickE;
int iVertex = 0; // variable to use as index into valsE array
const int fiberDimE = spaceDim; // number of values per point
const PylithScalar tolerance = 1.0e-06;
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
- != verticesEnd;
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
+ v_iter != verticesEnd;
++v_iter, ++iVertex) { // loop over fault vertices
// Check fiber dimension (number of values at point)
const int fiberDim = slipSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ CPPUNIT_ASSERT(vals);
// Check values at point
for (int i = 0; i < fiberDimE; ++i) {
const int index = iVertex * spaceDim + i;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
+ const PylithScalar valE = valsE[index];
+#if 0 // DEBUGGING
+ std::cout << "SLIP valE: " << valE
+ << ", val: " << vals[i]
+ << ", error: " << fabs(1.0-vals[i]/valE)
+ << std::endl;
+#endif // DEBUGGING
+ if (fabs(valE) > tolerance)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
} // for
} // for
} // Check slip values
@@ -321,6 +335,54 @@
fault.timeStep(dt);
fault.constrainSolnSpace(&fields, t, jacobian);
+ { // Check solution values
+ // Lagrange multipliers should be adjusted according to friction
+ // as reflected in the fieldIncrSlipE data member.
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ // Get section containing solution (disp + Lagrange multipliers)
+ const ALE::Obj<RealSection>& dispIncrSection =
+ fields.get("dispIncr(t->t+dt)").section();
+ CPPUNIT_ASSERT(!dispIncrSection.isNull());
+
+ // Get expected values
+ const PylithScalar* valsE = _data->fieldIncrSlipE; // Expected values for dispIncr
+ int iVertex = 0; // variable to use as index into valsE array
+ const int fiberDimE = spaceDim; // number of values per point
+ const PylithScalar tolerance = 1.0e-06;
+ for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter, ++iVertex) { // loop over all vertices in mesh
+ // Check fiber dimension (number of values at point)
+ const int fiberDim = dispIncrSection->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const PylithScalar* vals = dispIncrSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ // Check values at point
+ for (int i = 0; i < fiberDimE; ++i) {
+ const int index = iVertex * spaceDim + i;
+ const PylithScalar valE = valsE[index];
+#if 0 // DEBUGGING
+ std::cout << "SOLUTION valE: " << valE
+ << ", val: " << vals[i]
+ << ", error: " << fabs(1.0-vals[i]/valE)
+ << std::endl;
+#endif // DEBUGGING
+ if (fabs(valE) > tolerance)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ else
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
+ } // for
+ } // for
+ } // Check solution values
+
{ // Check slip values
// Slip values should be adjusted based on the change in the
// Lagrange multipliers as reflected in the slipSlipE data member.
@@ -374,54 +436,6 @@
} // for
} // Check slip values
- { // Check solution values
- // Lagrange multipliers should be adjusted according to friction
- // as reflected in the fieldIncrSlipE data member.
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
-
- // Get section containing solution (disp + Lagrange multipliers)
- const ALE::Obj<RealSection>& dispIncrSection =
- fields.get("dispIncr(t->t+dt)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
-
- // Get expected values
- const PylithScalar* valsE = _data->fieldIncrSlipE; // Expected values for dispIncr
- int iVertex = 0; // variable to use as index into valsE array
- const int fiberDimE = spaceDim; // number of values per point
- const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-05;
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) { // loop over all vertices in mesh
- // Check fiber dimension (number of values at point)
- const int fiberDim = dispIncrSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = dispIncrSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
-
- // Check values at point
- for (int i = 0; i < fiberDimE; ++i) {
- const int index = iVertex * spaceDim + i;
- const PylithScalar valE = valsE[index];
-#if 0 // DEBUGGING
- std::cout << "SOLUTION valE: " << valE
- << ", val: " << vals[i]
- << ", error: " << fabs(1.0-vals[i]/valE)
- << std::endl;
-#endif // DEBUGGING
- if (fabs(valE) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
- } // for
- } // for
- } // Check solution values
-
} // testConstrainSolnSpaceSlip
// ----------------------------------------------------------------------
@@ -737,7 +751,7 @@
mesh->coordsys(&cs);
mesh->nondimensionalize(normalizer);
- const PylithScalar upDir[] = { 0.0, 0.0, 1.0 };
+ const PylithScalar upDir[3] = { 0.0, 0.0, 1.0 };
fault->initialize(*mesh, upDir);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynData.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynData.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynData.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -42,6 +42,7 @@
orientation(0),
initialTractions(0),
area(0),
+ slipStickE(0),
fieldIncrSlipE(0),
slipSlipE(0),
fieldIncrOpenE(0),
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynData.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynData.hh 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynData.hh 2012-01-18 21:51:13 UTC (rev 19381)
@@ -78,6 +78,7 @@
PylithScalar* orientation; ///< Expected values for fault orientation.
PylithScalar* area; ///< Expected values for fault area.
PylithScalar* initialTractions; ///< Expected values for initial tractions.
+ PylithScalar* slipStickE; ///< Expected values for slip for sticking case.
PylithScalar* fieldIncrSlipE; ///< Expected values for solution increment for slipping case.
PylithScalar* slipSlipE; ///< Expected values for slip for slipping case.
PylithScalar* fieldIncrOpenE; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -1128,18 +1128,23 @@
0.0, 2.0, 1.0,
1.1, 3.1, 2.1,
1.2, 3.2, 2.2,
- 1.3, 3.3, 2.3, // 14
- 1.5, 3.5, 2.5, // 15
- 1.7, 3.7, 2.7, // 16
- 1.9, 3.9, 2.9, // 17
- 41.4, 3.4, 2.4, // 18
- 41.6, 3.6, 2.6, // 19
- 41.8, 3.8, 2.8, // 20
- 41.0, 3.0, 2.0, // 21
+ 0.5, 3.3, 2.3, // 14
+ 0.6, 3.5, 2.5, // 15
+ 0.7, 3.7, 2.7, // 16
+ 0.8, 3.9, 2.9, // 17
+ 12.266666666667, 3.6, 4.6, // 18
+ 12.066666666667, 5.4, 2.4, // 19
+ 16.866666666667, 2.2, 8.2, // 20
+ 17.666666666667, 10.0, 2.0, // 21
};
-// No change in fieldIncr
-// Zero slip
+// No slip
+const PylithScalar pylith::faults::CohesiveDynDataHex8::_slipStickE[] = {
+ -0.8, 0.8, 0.0,
+ -0.9, 0.9, 0.0,
+ -1.0, 1.0, 0.0,
+ -1.1, 1.1, 0.0,
+};
// ----------------------------------------------------------------------
// Slip case
@@ -1285,6 +1290,7 @@
// Stick
fieldIncrStick = const_cast<PylithScalar*>(_fieldIncrStick);
+ slipStickE = const_cast<PylithScalar*>(_slipStickE);
// Slip
fieldIncrSlip = const_cast<PylithScalar*>(_fieldIncrSlip);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.hh 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.hh 2012-01-18 21:51:13 UTC (rev 19381)
@@ -68,6 +68,7 @@
static const PylithScalar _orientation[]; ///< Expected values for fault orientation.
static const PylithScalar _area[]; ///< Expected values for fault area.
static const PylithScalar _initialTractions[]; ///< Expected values for initial tractions.
+ static const PylithScalar _slipStickE[]; ///< Expected values for slip for stick case.
static const PylithScalar _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const PylithScalar _slipSlipE[]; ///< Expected values for slip for slip case.
static const PylithScalar _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -302,14 +302,17 @@
1.4, 2.4, // 5
1.5, 2.5,
1.6, 2.6,
- 1.7, 2.7, // 8
- 1.9, 2.9, // 9
- 21.8, 2.8, // 10
- 21.0, 2.0, // 11
+ 1.3, 2.7, // 8
+ 1.4, 2.9, // 9
+ 11.2, 2.20, // 10
+ 16.0, 5.40, // 11
};
-// No change in fieldIncr
-// Zero slip
+// No change in slip
+const PylithScalar pylith::faults::CohesiveDynDataQuad4::_slipStickE[] = {
+ -0.4, 0.0,
+ -0.5, 0.0,
+};
// ----------------------------------------------------------------------
// Slip case
@@ -411,6 +414,7 @@
// Stick
fieldIncrStick = const_cast<PylithScalar*>(_fieldIncrStick);
+ slipStickE = const_cast<PylithScalar*>(_slipStickE);
// Slip
fieldIncrSlip = const_cast<PylithScalar*>(_fieldIncrSlip);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.hh 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.hh 2012-01-18 21:51:13 UTC (rev 19381)
@@ -68,6 +68,7 @@
static const PylithScalar _orientation[]; ///< Expected values for fault orientation.
static const PylithScalar _area[]; ///< Expected values for fault area.
static const PylithScalar _initialTractions[]; ///< Expected values for initial tractions.
+ static const PylithScalar _slipStickE[]; ///< Expected values for slip for stick case.
static const PylithScalar _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const PylithScalar _slipSlipE[]; ///< Expected values for slip for slip case.
static const PylithScalar _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -506,16 +506,20 @@
1.3, 2.3, 3.3, // 4
1.4, 2.4, 3.4, // 5
1.5, 2.5, 3.5,
- 1.6, 2.6, 3.6, // 7
- 1.8, 2.8, 3.8, // 8
- 1.0, 2.0, 3.0, // 9
+ 1.2, 2.2, 3.2, // 7
+ 1.3, 2.3, 3.3, // 8
+ 1.4, 2.4, 3.4, // 9
41.7, 2.7, 3.7, // 10
41.9, 2.9, 3.9, // 11
41.1, 2.1, 3.1, // 12
};
-// No change in fieldIncr
-// Zero slip
+// No slip
+const PylithScalar pylith::faults::CohesiveDynDataTet4::_slipStickE[] = {
+ 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0,
+ 0.0, 0.0, 0.0,
+};
// ----------------------------------------------------------------------
// Slip case
@@ -623,6 +627,7 @@
// Stick
fieldIncrStick = const_cast<PylithScalar*>(_fieldIncrStick);
+ slipStickE = const_cast<PylithScalar*>(_slipStickE);
// Slip
fieldIncrSlip = const_cast<PylithScalar*>(_fieldIncrSlip);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.hh 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.hh 2012-01-18 21:51:13 UTC (rev 19381)
@@ -68,6 +68,7 @@
static const PylithScalar _orientation[]; ///< Expected values for fault orientation.
static const PylithScalar _area[]; ///< Expected values for fault area.
static const PylithScalar _initialTractions[]; ///< Expected values for initial tractions.
+ static const PylithScalar _slipStickE[]; ///< Expected values for slip for stick case.
static const PylithScalar _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const PylithScalar _slipSlipE[]; ///< Expected values for slip for slip case.
static const PylithScalar _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -276,8 +276,11 @@
21.8, 2.8, // 9
};
-// No change in fieldIncr
-// Zero slip
+// No slip
+const PylithScalar pylith::faults::CohesiveDynDataTri3::_slipStickE[] = {
+ 0.0, 0.0,
+ 0.0, 0.0,
+};
// ----------------------------------------------------------------------
// Slip case
@@ -371,6 +374,7 @@
// Stick
fieldIncrStick = const_cast<PylithScalar*>(_fieldIncrStick);
+ slipStickE = const_cast<PylithScalar*>(_slipStickE);
// Slip
fieldIncrSlip = const_cast<PylithScalar*>(_fieldIncrSlip);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.hh 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.hh 2012-01-18 21:51:13 UTC (rev 19381)
@@ -68,6 +68,7 @@
static const PylithScalar _orientation[]; ///< Expected values for fault orientation.
static const PylithScalar _area[]; ///< Expected values for fault area.
static const PylithScalar _initialTractions[]; ///< Expected values for initial tractions.
+ static const PylithScalar _slipStickE[]; ///< Expected values for slip for stick case.
static const PylithScalar _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const PylithScalar _slipSlipE[]; ///< Expected values for slip for slip case.
static const PylithScalar _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -454,16 +454,20 @@
1.4, 2.4,
1.5, 2.5, // 8
1.6, 2.6,
- 1.7, 2.7, // 10
- 1.9, 2.9, // 11
- 2.1, 2.1, // 12
+ 1.2, 2.2, // 10
+ 1.3, 2.3, // 11
+ 1.5, 2.5, // 12
21.8, 22.8, // 13
21.0, 2.0, // 14
2.2, 22.2, // 15
};
-// No change in fieldIncr
-// Zero slip
+// No slip
+const PylithScalar pylith::faults::CohesiveDynDataTri3d::_slipStickE[] = {
+ 0.0, 0.0,
+ 0.0, 0.0,
+ 0.0, 0.0,
+};
// ----------------------------------------------------------------------
// Slip case
@@ -575,6 +579,7 @@
// Stick
fieldIncrStick = const_cast<PylithScalar*>(_fieldIncrStick);
+ slipStickE = const_cast<PylithScalar*>(_slipStickE);
// Slip
fieldIncrSlip = const_cast<PylithScalar*>(_fieldIncrSlip);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.hh 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.hh 2012-01-18 21:51:13 UTC (rev 19381)
@@ -68,6 +68,7 @@
static const PylithScalar _orientation[]; ///< Expected values for fault orientation.
static const PylithScalar _area[]; ///< Expected values for fault area.
static const PylithScalar _initialTractions[]; ///< Expected values for initial tractions.
+ static const PylithScalar _slipStickE[]; ///< Expected values for slip for stick case.
static const PylithScalar _fieldIncrSlipE[]; ///< Expected values for solution increment for slip case.
static const PylithScalar _slipSlipE[]; ///< Expected values for slip for slip case.
static const PylithScalar _fieldIncrOpenE[]; ///< Expected values for solution increment for opening case.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/friction/TestFrictionModel.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/friction/TestFrictionModel.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/friction/TestFrictionModel.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -784,10 +784,10 @@
const int numCorners = 2;
const int numQuadPts = 2;
const int spaceDim = 2;
- const PylithScalar basis[] = { 0.75, 0.25, 0.25, 0.75 };
- const PylithScalar basisDeriv[] = { -0.5, 0.5, -0.5, 0.5 };
- const PylithScalar quadPtsRef[] = { -0.5, 0.5 };
- const PylithScalar quadWts[] = { 1.0, 1.0 };
+ const PylithScalar basis[4] = { 1.0, 0.0, 0.0, 1.0 };
+ const PylithScalar basisDeriv[4] = { -0.5, 0.5, -0.5, 0.5 };
+ const PylithScalar quadPtsRef[2] = { -1.0, 1.0 };
+ const PylithScalar quadWts[2] = { 1.0, 1.0 };
quadrature.initialize(basis, numQuadPts, numCorners,
basisDeriv, numQuadPts, numCorners, cellDim,
quadPtsRef, numQuadPts, cellDim,
@@ -818,8 +818,9 @@
friction->normalizer(normalizer);
fault->frictionModel(friction);
- const PylithScalar upDir[] = { 0.0, 0.0, 1.0 };
+ const PylithScalar upDir[3] = { 0.0, 0.0, 1.0 };
fault->initialize(*mesh, upDir);
+ fault->verifyConfiguration(*mesh);
} // _initialize
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestOutputSolnPoints.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestOutputSolnPoints.cc 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/meshio/TestOutputSolnPoints.cc 2012-01-18 21:51:13 UTC (rev 19381)
@@ -54,11 +54,11 @@
const char* filename = "data/quad4.mesh";
const int numPoints = 5;
const PylithScalar points[10] = {
- 0.0, 0.1,
- 0.3, 0.4,
- 0.6, 0.7,
- 1.0, 1.1,
- 1.3, 1.4,
+ 0.0, 0.1,
+ 0.3, 0.4,
+ -0.6, -0.7,
+ -1.0, 0.9,
+ -0.3, 0.8,
};
const int nvertices = numPoints;
const int verticesE[5] = { 5, 6, 7, 8, 9 };
@@ -72,11 +72,9 @@
cs.setSpaceDim(spaceDim);
cs.initialize();
mesh.coordsys(&cs);
-#if 0
MeshIOAscii iohandler;
iohandler.filename("data/quad4.mesh");
iohandler.read(&mesh);
-#endif
OutputSolnPoints output;
output.setupInterpolator(&mesh, points, numPoints, spaceDim);
@@ -129,11 +127,11 @@
const char* filename = "data/quad4.mesh";
const int numPoints = 5;
const PylithScalar points[15] = {
- 0.0, 0.1, 0.2,
- 0.3, 0.4, 0.5,
- 0.6, 0.7, 0.8,
- 1.0, 1.1, 1.2,
- 1.3, 1.4, 1.5,
+ 0.0, 0.1, 0.2,
+ 0.3, -0.4, 0.5,
+ 0.6, 0.7, -0.8,
+ -1.0, 0.1, -0.3,
+ 0.3, 0.8, 0.5,
};
const int nvertices = numPoints;
const int verticesE[5] = { 5, 6, 7, 8, 9 };
@@ -147,11 +145,9 @@
cs.setSpaceDim(spaceDim);
cs.initialize();
mesh.coordsys(&cs);
-#if 0
MeshIOAscii iohandler;
iohandler.filename("data/hex8.mesh");
iohandler.read(&mesh);
-#endif
OutputSolnPoints output;
output.setupInterpolator(&mesh, points, numPoints, spaceDim);
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/Makefile.am 2012-01-18 01:51:38 UTC (rev 19380)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/meshio/Makefile.am 2012-01-18 21:51:13 UTC (rev 19381)
@@ -37,6 +37,7 @@
TestOutputSolnSubset.py \
TestDataWriterVTK.py \
TestDataWriterHDF5.py \
+ TestDataWriterHDF5Ext.py \
TestSingleOutput.py \
TestXdmf.py
More information about the CIG-COMMITS
mailing list