[cig-commits] r20036 - in short/3D/PyLith/branches/v1.7-trunk: examples/3d/hex8 examples/3d/hex8/spatialdb libsrc/pylith/faults modulesrc/faults pylith/faults unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Thu May 3 15:26:48 PDT 2012
Author: brad
Date: 2012-05-03 15:26:48 -0700 (Thu, 03 May 2012)
New Revision: 20036
Added:
short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/spatialdb/tractions_opening.spatialdb
short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/step20.cfg
Modified:
short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/README
short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/step19.cfg
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
short/3D/PyLith/branches/v1.7-trunk/modulesrc/faults/FaultCohesiveDyn.i
short/3D/PyLith/branches/v1.7-trunk/pylith/faults/FaultCohesiveDyn.py
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh
short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
Log:
Added flag to allow initial tractions even when fault opens. Added examples/3d/hex8/step20 to illustrate its use for a dike intrusion.
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/Makefile.am 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/Makefile.am 2012-05-03 22:26:48 UTC (rev 20036)
@@ -58,7 +58,8 @@
spatialdb/powerlaw/powerlaw_gendb.cfg \
spatialdb/powerlaw/powerlaw_params.spatialdb \
spatialdb/powerlaw/powerlaw_points.txt \
- spatialdb/powerlaw/temperature.spatialdb
+ spatialdb/powerlaw/temperature.spatialdb \
+ spatialdb/tractions_opening.spatialdb
SUBDIRS = \
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/README
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/README 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/README 2012-05-03 22:26:48 UTC (rev 20036)
@@ -1,4 +1,4 @@
-The examples in this directory form a step-by-step sequence of 19 problems,
+The examples in this directory form a step-by-step sequence of 20 problems,
each building on the one before (for the most part). All of the examples
use the same mesh, which was created by Cubit. The mesh is 6 km x 6 km x 4
km with linear hexahedral cells that have edges 1.0 km long.
@@ -68,6 +68,7 @@
rheology (quasi-static)
step18: Axial traction surface load on top surface (static)
step19: Time dependent axial traction surface load on top surface (quasi-static)
+step20: Dike opening via initial tractions and fault constitutive model
----------------------------------------
mesh directory
@@ -160,3 +161,7 @@
desired. The powerlaw_gendb.py utility code is used from this
directory to create spatialdb/mat_powerlaw.spatialdb.
+tractions_opening.spatialdb
+
+ Spatial database defining initial tractions for example step20.
+
Added: short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/spatialdb/tractions_opening.spatialdb
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/spatialdb/tractions_opening.spatialdb (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/spatialdb/tractions_opening.spatialdb 2012-05-03 22:26:48 UTC (rev 20036)
@@ -0,0 +1,44 @@
+// -*- C++ -*- (tell Emacs to use C++ mode for syntax highlighting)
+//
+// This spatial database specifies the distribution of initial
+// tractions on the fault surface. We want to mimic a dike intrusion,
+// so we impose tensile (positive) normal tractions in the central
+// portion of the fault surface. We want to keep the rest of the fault
+// surface in contact with zero slip so we impose large compressive
+// (negative) normal tractions.
+//
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 3
+ value-names = traction-shear-leftlateral traction-shear-updip traction-normal
+ value-units = MPa MPa MPa
+ num-locs = 12
+ data-dim = 2 // Locations of data points are on a surface.
+ space-dim = 3
+ cs-data = cartesian {
+ to-meters = 1.0e+3 // Specify coordinates in km for convenience.
+ space-dim = 3
+ } // cs-data
+} // SimpleDB
+// Columns are
+// (1) x coordinate (km)
+// (2) y coordinate (km)
+// (3) z coordinate (km)
+// (4) left-lateral-slip (m) (right-lateral is negative)
+// (5) reverse-slip (m)
+// (6) fault-opening (m)
+0.0 -1.1 0.0 0.00 0.00 -100.00
+0.0 -1.1 -2.0 0.00 0.00 -100.00
+0.0 -1.1 -2.1 0.00 0.00 -100.00
+
+0.0 -1.0 0.0 0.00 0.00 10.00
+0.0 -1.0 -2.0 0.00 0.00 10.00
+0.0 -1.0 -2.1 0.00 0.00 -100.00
+
+0.0 +1.0 0.0 0.00 0.00 10.00
+0.0 +1.0 -2.0 0.00 0.00 10.00
+0.0 +1.0 -2.1 0.00 0.00 -100.00
+
+0.0 +1.1 0.0 0.00 0.00 -100.00
+0.0 +1.1 -2.0 0.00 0.00 -100.00
+0.0 +1.1 -2.1 0.00 0.00 -100.00
Modified: short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/step19.cfg
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/step19.cfg 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/step19.cfg 2012-05-03 22:26:48 UTC (rev 20036)
@@ -1,20 +1,20 @@
[pylithapp]
-
+#
# ----------------------------------------------------------------------
# PROBLEM DESCRIPTION
# ----------------------------------------------------------------------
-
#
+#
# This is a purely elastic quasi-static problem using time-dependent
# Neumann (traction) boundary conditions. We apply normal tractions to
# the top surface that increase and then decrease while applying
# roller (fixed normal but free lateral motion) boundary conditions on
# the lateral sides and bottom surfaces.
-
+#
# ----------------------------------------------------------------------
# RUNNING THE SIMULATON
# ----------------------------------------------------------------------
-
+#
# This is not a self-contained simulation configuration file. This
# file only specifies parameters specific to tutorial step19.
# The general parameters are specificed in the pylithapp.cfg
@@ -24,7 +24,7 @@
# pylith step19.cfg
#
# Output will be directed to directory output.
-
+#
# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
Added: short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/step20.cfg
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/step20.cfg (rev 0)
+++ short/3D/PyLith/branches/v1.7-trunk/examples/3d/hex8/step20.cfg 2012-05-03 22:26:48 UTC (rev 20036)
@@ -0,0 +1,165 @@
+[pylithapp]
+#
+# ----------------------------------------------------------------------
+# PROBLEM DESCRIPTION
+# ----------------------------------------------------------------------
+#
+# This is a pure elastic static problem that demonstrates how to use
+# initial tractions with the FaultCohesiveDyn object to create fault
+# opening consistent with a fluid intrusion, such as a dike.
+#
+# ----------------------------------------------------------------------
+# RUNNING THE SIMULATON
+# ----------------------------------------------------------------------
+#
+# This is not a self-contained simulation configuration file. This
+# file only specifies parameters specific to tutorial step20.
+# The general parameters are specificed in the pylithapp.cfg
+# file which PyLith reads by default.
+#
+# To run the simulation:
+# pylith step20.cfg
+#
+# Output will be directed to directory output.
+#
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[pylithapp.timedependent.formulation.time_step]
+# Define the total time for the simulation and the default time step size.
+total_time = 0.0*s ; total time of simulation
+# Define an appropriat time step for simulations. Important for
+# nondimensionalization of velocities and slip rates.
+dt = 1.0*year
+
+[pylithapp.timedependent]
+# Set bc to an array of 3 boundary conditions: 'x_pos','x_neg', and 'z_neg'.
+bc = [x_pos,x_neg,z_neg]
+
+# Set interfaces to an array of 1 fault: 'fault'.
+interfaces = [fault]
+
+[pylithapp.timedependent.implicit]
+# Set the output to an array of 2 output managers.
+# We will output the solution over the domain and the ground surface.
+output = [domain,subdomain]
+
+# Set subdomain component to OutputSolnSubset (boundary of the domain).
+output.subdomain = pylith.meshio.OutputSolnSubset
+
+# Fault friction is a nonlinear problem so we need to use the nonlinear
+# solver.
+solver = pylith.problems.SolverNonlinear
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+# Set the parameters for Dirichlet boundary conditions applied on the
+# -x, +x, and -z faces of the box.
+#
+# We fix the x and y degrees of freedom on the -x and +x faces, and
+# fix the z degree of freedom on the bottom (-z) face.
+#
+# For all boundaries, we retain the default ZeroDispDB, which specifies
+# a zero value.
+#
+
+# The label corresponds to the name of the nodeset in CUBIT.
+
+# +x face
+[pylithapp.timedependent.bc.x_pos]
+bc_dof = [0,1]
+label = face_xpos
+db_initial.label = Dirichlet BC on +x
+
+# -x face
+[pylithapp.timedependent.bc.x_neg]
+bc_dof = [0,1]
+label = face_xneg
+db_initial.label = Dirichlet BC on -x
+
+# -z face
+[pylithapp.timedependent.bc.z_neg]
+bc_dof = [2]
+label = face_zneg_nofault
+db_initial.label = Dirichlet BC on -z
+
+# ----------------------------------------------------------------------
+# faults
+# ----------------------------------------------------------------------
+[pylithapp.timedependent.interfaces]
+# Change fault to dynamic fault interface.
+fault = pylith.faults.FaultCohesiveDyn
+
+[pylithapp.timedependent.interfaces.fault]
+# The label corresponds to the name of the nodeset in CUBIT.
+label = fault
+
+# Use the static friction model.
+friction = pylith.friction.StaticFriction
+friction.label = Static friction
+
+# We must define the quadrature information for fault cells.
+# The fault cells are 2D (surface).
+quadrature.cell = pylith.feassemble.FIATLagrange
+quadrature.cell.dimension = 2
+
+# Set static friction model parameters using a uniform DB. Set the
+# static coefficient of friction to 0.6 and cohesion to 0.0 Pa.
+friction.db_properties = spatialdata.spatialdb.UniformDB
+friction.db_properties.label = Static friction
+friction.db_properties.values = [friction-coefficient,cohesion]
+friction.db_properties.data = [0.1,0.0*Pa]
+
+db_initial_tractions = spatialdata.spatialdb.SimpleDB
+db_initial_tractions.label = Initial fault tractions
+db_initial_tractions.iohandler.filename = spatialdb/tractions_opening.spatialdb
+db_initial_tractions.query_type = nearest
+
+# ----------------------------------------------------------------------
+# PETSc settings
+# ----------------------------------------------------------------------
+# NOTE: There are additional settings specific to fault friction.
+[pylithapp.petsc]
+
+# 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 = 25
+friction_ksp_gmres_restart = 30
+# Uncomment to view details of friction sensitivity solve.
+#friction_ksp_monitor = true
+#friction_ksp_view = true
+friction_ksp_converged_reason = true
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+# Give basename for VTK domain output of solution over domain.
+[pylithapp.problem.formulation.output.domain.writer]
+filename = output/step20.vtk
+
+# Give basename for VTK domain output of solution over ground surface.
+[pylithapp.problem.formulation.output.subdomain]
+# Name of nodeset for ground surface.
+label = face_zpos
+writer.filename = output/step20-groundsurf.vtk
+
+# Give basename for VTK fault output.
+[pylithapp.problem.interfaces.fault.output]
+vertex_info_fields = [initial_traction]
+writer.filename = output/step20-fault.vtk
+
+# Give basename for VTK output of upper_crust state variables.
+[pylithapp.timedependent.materials.upper_crust.output]
+# Average values over quadrature points.
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = output/step20-upper_crust.vtk
+
+# Give basename for VTK output of lower_crust state variables.
+[pylithapp.timedependent.materials.lower_crust.output]
+# Average values over quadrature points.
+cell_filter = pylith.meshio.CellFilterAvgMesh
+writer.filename = output/step20-lower_crust.vtk
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-05-03 22:26:48 UTC (rev 20036)
@@ -63,6 +63,7 @@
// Default constructor.
pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(void) :
_zeroTolerance(1.0e-10),
+ _openFreeSurf(true),
_dbInitialTract(0),
_friction(0),
_jacobian(0),
@@ -87,7 +88,7 @@
_friction = 0; // :TODO: Use shared pointer
delete _jacobian; _jacobian = 0;
- if (0 != _ksp) {
+ if (_ksp) {
PetscErrorCode err = KSPDestroy(&_ksp); _ksp = 0;
CHECK_PETSC_ERROR(err);
} // if
@@ -125,14 +126,23 @@
} // zeroTolerance
// ----------------------------------------------------------------------
+// Set flag used to determine when fault is traction free when it
+// opens or it still imposes any initial tractions.
+void
+pylith::faults::FaultCohesiveDyn::openFreeSurf(const bool value)
+{ // openFreeSurf
+ _openFreeSurf = value;
+} // openFreeSurf
+
+// ----------------------------------------------------------------------
// Initialize fault. Determine orientation and setup boundary
void
pylith::faults::FaultCohesiveDyn::initialize(const topology::Mesh& mesh,
const PylithScalar upDir[3])
{ // initialize
- assert(0 != upDir);
- assert(0 != _quadrature);
- assert(0 != _normalizer);
+ assert(upDir);
+ assert(_quadrature);
+ assert(_normalizer);
FaultCohesiveLagrange::initialize(mesh, upDir);
@@ -140,14 +150,14 @@
_setupInitialTractions();
// Setup fault constitutive model.
- assert(0 != _friction);
- assert(0 != _faultMesh);
- assert(0 != _fields);
+ assert(_friction);
+ assert(_faultMesh);
+ assert(_fields);
_friction->normalizer(*_normalizer);
_friction->initialize(*_faultMesh, _quadrature);
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- assert(0 != cs);
+ assert(cs);
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("FaultFields");
@@ -172,9 +182,9 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateResidual
- assert(0 != fields);
- assert(0 != _fields);
- assert(0 != _logger);
+ assert(fields);
+ assert(_fields);
+ assert(_logger);
// Cohesive cells with conventional vertices N and P, and constraint
// vertex L make contributions to the assembled residual:
@@ -332,7 +342,10 @@
residualVertexN = 0.0;
residualVertexL = 0.0;
- if (slipNormal < _zeroTolerance) { // if no opening
+ if (slipNormal < _zeroTolerance || !_openFreeSurf) {
+ // if no opening or flag indicates to still impose initial
+ // tractions when fault is open.
+ //
// Initial (external) tractions oppose (internal) tractions
// associated with Lagrange multiplier.
residualVertexN = areaVertex * (dispTpdtVertexL - initialTractionsVertex);
@@ -381,8 +394,8 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // updateStateVars
- assert(0 != fields);
- assert(0 != _fields);
+ assert(fields);
+ assert(_fields);
_updateRelMotion(*fields);
@@ -1127,8 +1140,8 @@
const scalar_array&,
const bool);
- assert(0 != fields);
- assert(0 != _quadrature);
+ assert(fields);
+ assert(_quadrature);
// Cohesive cells with conventional vertices i and j, and constraint
// vertex k require three adjustments to the solution:
@@ -1462,11 +1475,11 @@
pylith::faults::FaultCohesiveDyn::vertexField(const char* name,
const topology::SolutionFields* fields)
{ // vertexField
- assert(0 != _faultMesh);
- assert(0 != _quadrature);
- assert(0 != _normalizer);
- assert(0 != _fields);
- assert(0 != _friction);
+ assert(_faultMesh);
+ assert(_quadrature);
+ assert(_normalizer);
+ assert(_fields);
+ assert(_friction);
const int cohesiveDim = _faultMesh->dimension();
const int spaceDim = _quadrature->spaceDim();
@@ -1541,7 +1554,7 @@
return buffer;
} else if (0 == strcasecmp("initial_traction", name)) {
- assert(0 != _dbInitialTract);
+ assert(_dbInitialTract);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer =
_fields->get("buffer (vector)");
@@ -1552,7 +1565,7 @@
return buffer;
} else if (0 == strcasecmp("traction", name)) {
- assert(0 != fields);
+ assert(fields);
const topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer =
@@ -1574,7 +1587,7 @@
throw std::logic_error("Unknown field in FaultCohesiveDyn::vertexField().");
// Satisfy return values
- assert(0 != _fields);
+ assert(_fields);
const topology::Field<topology::SubMesh>& buffer = _fields->get(
"buffer (vector)");
@@ -1585,13 +1598,13 @@
void
pylith::faults::FaultCohesiveDyn::_setupInitialTractions(void)
{ // _setupInitialTractions
- assert(0 != _normalizer);
+ assert(_normalizer);
// If no initial tractions specified, leave method
if (0 == _dbInitialTract)
return;
- assert(0 != _normalizer);
+ assert(_normalizer);
const PylithScalar pressureScale = _normalizer->pressureScale();
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -1616,7 +1629,7 @@
assert(!orientationSection.isNull());
const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
- assert(0 != cs);
+ assert(cs);
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
@@ -1627,7 +1640,7 @@
assert(!coordsSection.isNull());
- assert(0 != _dbInitialTract);
+ assert(_dbInitialTract);
_dbInitialTract->open();
switch (spaceDim) { // switch
case 1: {
@@ -1713,10 +1726,10 @@
topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& dispT)
{ // _calcTractions
- assert(0 != tractions);
- assert(0 != _faultMesh);
- assert(0 != _fields);
- assert(0 != _normalizer);
+ assert(tractions);
+ assert(_faultMesh);
+ assert(_fields);
+ assert(_normalizer);
// Fiber dimension of tractions matches spatial dimension.
const int spaceDim = _quadrature->spaceDim();
@@ -1791,7 +1804,7 @@
void
pylith::faults::FaultCohesiveDyn::_updateRelMotion(const topology::SolutionFields& fields)
{ // _updateRelMotion
- assert(0 != _fields);
+ assert(_fields);
const int spaceDim = _quadrature->spaceDim();
@@ -1885,8 +1898,8 @@
void
pylith::faults::FaultCohesiveDyn::_sensitivitySetup(const topology::Jacobian& jacobian)
{ // _sensitivitySetup
- assert(0 != _fields);
- assert(0 != _quadrature);
+ assert(_fields);
+ assert(_quadrature);
const int spaceDim = _quadrature->spaceDim();
@@ -1934,7 +1947,7 @@
// Setup Jacobian sparse matrix for sensitivity solve.
if (0 == _jacobian)
_jacobian = new topology::Jacobian(solution, jacobian.matrixType());
- assert(0 != _jacobian);
+ assert(_jacobian);
_jacobian->zero();
// Setup PETSc KSP linear solver.
@@ -1970,8 +1983,8 @@
const topology::Jacobian& jacobian,
const topology::SolutionFields& fields)
{ // _sensitivityUpdateJacobian
- assert(0 != _quadrature);
- assert(0 != _fields);
+ assert(_quadrature);
+ assert(_fields);
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
@@ -1997,7 +2010,7 @@
// Visitor for Jacobian matrix associated with domain.
scalar_array jacobianSubCell(submatrixSize);
const PetscMat jacobianDomainMatrix = jacobian.matrix();
- assert(0 != jacobianDomainMatrix);
+ assert(jacobianDomainMatrix);
const ALE::Obj<SieveMesh::order_type>& globalOrderDomain =
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solutionDomainSection);
assert(!globalOrderDomain.isNull());
@@ -2020,9 +2033,9 @@
assert(!solutionFaultSection.isNull());
// Visitor for Jacobian matrix associated with fault.
- assert(0 != _jacobian);
+ assert(_jacobian);
const PetscMat jacobianFaultMatrix = _jacobian->matrix();
- assert(0 != jacobianFaultMatrix);
+ assert(jacobianFaultMatrix);
const ALE::Obj<SieveSubMesh::order_type>& globalOrderFault =
faultSieveMesh->getFactory()->getGlobalOrder(faultSieveMesh, "default", solutionFaultSection);
assert(!globalOrderFault.isNull());
@@ -2041,7 +2054,7 @@
const int coneSize = ncV.getSize();
assert(coneSize == 3*numBasis);
const SieveMesh::point_type *cohesiveCone = ncV.getPoints();
- assert(0 != cohesiveCone);
+ assert(cohesiveCone);
const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
jacobianSubCell = 0.0;
@@ -2196,9 +2209,9 @@
void
pylith::faults::FaultCohesiveDyn::_sensitivitySolve(void)
{ // _sensitivitySolve
- assert(0 != _fields);
- assert(0 != _jacobian);
- assert(0 != _ksp);
+ assert(_fields);
+ assert(_jacobian);
+ assert(_ksp);
const topology::Field<topology::SubMesh>& residual =
_fields->get("sensitivity residual");
@@ -2232,8 +2245,8 @@
void
pylith::faults::FaultCohesiveDyn::_sensitivityUpdateSoln(const bool negativeSide)
{ // _sensitivityUpdateSoln
- assert(0 != _fields);
- assert(0 != _quadrature);
+ assert(_fields);
+ assert(_quadrature);
const int spaceDim = _quadrature->spaceDim();
@@ -2530,7 +2543,7 @@
const scalar_array& tractionTpdt,
const bool iterating)
{ // _constrainSolnSpace1D
- assert(0 != dTractionTpdt);
+ assert(dTractionTpdt);
if (fabs(slip[0]) < _zeroTolerance) {
// if compression, then no changes to solution
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-05-03 22:26:48 UTC (rev 20036)
@@ -73,6 +73,16 @@
*/
void zeroTolerance(const PylithScalar value);
+ /** Set flag used to determine when fault is traction free when it
+ * opens or it still imposes any initial tractions.
+ *
+ * If true, acts as a frictional contact. If false, one can simulate
+ * a dike opening.
+ *
+ * @param value Nondimensional tolerance
+ */
+ void openFreeSurf(const bool value);
+
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
*
@@ -265,6 +275,11 @@
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
+ /// Flag to control whether to continue to impose initial tractions
+ /// on the fault surface when it opens. If it is a frictional
+ /// contact, then it should be a free surface.
+ bool _openFreeSurf;
+
/// Minimum resolvable value accounting for roundoff errors
PylithScalar _zeroTolerance;
Modified: short/3D/PyLith/branches/v1.7-trunk/modulesrc/faults/FaultCohesiveDyn.i
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/modulesrc/faults/FaultCohesiveDyn.i 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/modulesrc/faults/FaultCohesiveDyn.i 2012-05-03 22:26:48 UTC (rev 20036)
@@ -59,6 +59,16 @@
*/
void zeroTolerance(const PylithScalar value);
+ /** Set flag used to determine when fault is traction free when it
+ * opens or it still imposes any initial tractions.
+ *
+ * If true, acts as a frictional contact. If false, one can simulate
+ * a dike opening.
+ *
+ * @param value Nondimensional tolerance
+ */
+ void openFreeSurf(const bool value);
+
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
*
Modified: short/3D/PyLith/branches/v1.7-trunk/pylith/faults/FaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/pylith/faults/FaultCohesiveDyn.py 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/pylith/faults/FaultCohesiveDyn.py 2012-05-03 22:26:48 UTC (rev 20036)
@@ -43,6 +43,9 @@
\b Properties
@li \b zero_tolerance Tolerance for detecting zero values.
+ @li \b open_free_surface If True, enforce traction free surface when
+ the fault opens, otherwise use initial tractions even when the
+ fault opens.
\b Facilities
@li \b db_initial_tractions Spatial database for initial tractions.
@@ -60,6 +63,11 @@
validator=pyre.inventory.greaterEqual(0.0))
zeroTolerance.meta['tip'] = "Tolerance for detecting zero values."
+ openFreeSurf = pyre.inventory.bool("open_free_surface", default=False)
+ openFreeSurf.meta['tip'] = "If True, enforce traction free surface when " \
+ "the fault opens, otherwise use initial tractions even when the " \
+ "fault opens."
+
db = pyre.inventory.facility("db_initial_tractions", family="spatial_database",
factory=NullComponent)
db.meta['tip'] = "Spatial database for initial tractions."
@@ -205,6 +213,7 @@
ModuleFaultCohesiveDyn.dbInitialTract(self, self.inventory.db)
ModuleFaultCohesiveDyn.frictionModel(self, self.inventory.friction)
ModuleFaultCohesiveDyn.zeroTolerance(self, self.inventory.zeroTolerance)
+ ModuleFaultCohesiveDyn.openFreeSurf(self, self.inventory.openFreeSurf)
self.output = self.inventory.output
return
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-05-03 22:26:48 UTC (rev 20036)
@@ -110,6 +110,20 @@
} // zeroTolerance
// ----------------------------------------------------------------------
+// Test openFreeSurf().
+void
+pylith::faults::TestFaultCohesiveDyn::testOpenFreeSurf(void)
+{ // testOpenFreeSurf
+ FaultCohesiveDyn fault;
+
+ CPPUNIT_ASSERT_EQUAL(true, fault._openFreeSurf); // default
+
+ const bool value = false;
+ fault.openFreeSurf(value);
+ CPPUNIT_ASSERT_EQUAL(value, fault._openFreeSurf);
+ } // testOpenFreeSurf
+
+// ----------------------------------------------------------------------
// Test initialize().
void
pylith::faults::TestFaultCohesiveDyn::testInitialize(void)
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2012-05-03 22:26:48 UTC (rev 20036)
@@ -54,6 +54,7 @@
CPPUNIT_TEST( testConstructor );
CPPUNIT_TEST( testDBInitialTract );
CPPUNIT_TEST( testZeroTolerance );
+ CPPUNIT_TEST( testOpenFreeSurf );
// Tests in derived classes:
// testInitialize()
@@ -93,6 +94,9 @@
/// Test zeroTolerance().
void testZeroTolerance(void);
+ /// Test openFreeSurf().
+ void testOpenFreeSurf(void);
+
/// Test initialize().
void testInitialize(void);
Modified: short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc 2012-05-03 21:40:02 UTC (rev 20035)
+++ short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc 2012-05-03 22:26:48 UTC (rev 20036)
@@ -263,7 +263,7 @@
{ // Integrate residual with disp increment.
fault.integrateResidual(residual, t, &fields);
- residual.view("RESIDUAL"); // DEBUGGING
+ //residual.view("RESIDUAL"); // DEBUGGING
// Check values
const PylithScalar* valsE = _data->residualIncr;
More information about the CIG-COMMITS
mailing list