[cig-commits] r20066 - in short/3D/PyLith/branches/pylith-scecdynrup: . examples/bar_shearwave/quad4 libsrc/pylith/bc libsrc/pylith/faults libsrc/pylith/problems modulesrc/bc modulesrc/faults pylith pylith/bc pylith/faults tests/2d/faultstrip tests/2d/frictionslide tests/3d/cyclicfriction unittests/libtests/faults unittests/libtests/faults/data unittests/pytests/faults unittests/pytests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Thu May 10 12:49:30 PDT 2012
Author: brad
Date: 2012-05-10 12:49:30 -0700 (Thu, 10 May 2012)
New Revision: 20066
Added:
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/TractPerturbation.i
short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/TractPerturbation.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestTractPerturbation.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestTractPerturbation.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/tri3_changetract.spatialdb
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/TestTractPerturbation.py
Modified:
short/3D/PyLith/branches/pylith-scecdynrup/TODO
short/3D/PyLith/branches/pylith-scecdynrup/examples/bar_shearwave/quad4/dynamic.cfg
short/3D/PyLith/branches/pylith-scecdynrup/examples/bar_shearwave/quad4/dynamic_slipweakening.cfg
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/bc/DirichletBoundary.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/bc/Neumann.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/faults/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/TractPerturbation.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/TractPerturbation.hh
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/bc/TimeDependent.i
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/FaultCohesiveDyn.i
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/faults.i
short/3D/PyLith/branches/pylith-scecdynrup/pylith/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/DirichletBoundary.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/Neumann.py
short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/FaultCohesiveDyn.py
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/dynamic_slipweakening.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/dynamic_timeweakening.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/pylithapp.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/ratestate.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/tension.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/fieldsplit.cfg
short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/pylithapp.cfg
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.hh
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/TestFaultCohesiveDyn.py
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/data/Makefile.am
short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/testfaults.py
Log:
Merge from trunk (v1.7-trunk).
Modified: short/3D/PyLith/branches/pylith-scecdynrup/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/TODO 2012-05-10 19:49:30 UTC (rev 20066)
@@ -17,6 +17,12 @@
CURRENT ISSUES/PRIORITIES
======================================================================
+* Add pytables to PyLith binary
+
+* Add pytables to PyLith installer
+
+* Setup fieldsplit for SolverNonlinear [Matt]
+
* 2-D materials
+ PowerLawPlaneStrain (power law plane strain ) [Charles]
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/bar_shearwave/quad4/dynamic.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/bar_shearwave/quad4/dynamic.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/bar_shearwave/quad4/dynamic.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -34,16 +34,19 @@
# Specify the initial tractions on the fault using a uniform DB.
# shear: 6.1 MPa (right-lateral)
# normal 10 MPa (compressive)
-db_initial_tractions = spatialdata.spatialdb.UniformDB
-db_initial_tractions.label = Initial fault tractions
-db_initial_tractions.values = [traction-shear,traction-normal]
-db_initial_tractions.data = [-6.1*MPa, -10.0*MPa]
+traction_perturbation = pylith.faults.TractPerturbation
+[pylithapp.timedependent.interfaces.fault.traction_perturbation]
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Initial fault tractions
+db_initial.values = [traction-shear,traction-normal]
+db_initial.data = [-6.1*MPa, -10.0*MPa]
+
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
[pylithapp.timedependent.interfaces.fault.output]
-vertex_info_fields=[initial_traction]
+vertex_info_fields=[traction_initial_value]
vertex_data_fields=[slip,slip_rate,traction]
skip = 1
writer.time_format = %05.2f
Modified: short/3D/PyLith/branches/pylith-scecdynrup/examples/bar_shearwave/quad4/dynamic_slipweakening.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/examples/bar_shearwave/quad4/dynamic_slipweakening.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/examples/bar_shearwave/quad4/dynamic_slipweakening.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -55,7 +55,7 @@
# Add friction model parameters to fault info file in addition to
# default values.
-vertex_info_fields = [strike_dir,normal_dir,initial_traction,static_coefficient,dynamic_coefficient,slip_weakening_parameter,cohesion]
+vertex_info_fields = [strike_dir,normal_dir,traction_initial_value,static_coefficient,dynamic_coefficient,slip_weakening_parameter,cohesion]
# Add output of state variables cumulative slip and previous slip to default
# values.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/bc/DirichletBoundary.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/bc/DirichletBoundary.cc 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/bc/DirichletBoundary.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -104,15 +104,15 @@
logger.stagePop();
- if (0 == strcasecmp(name, "initial-value"))
+ if (0 == strcasecmp(name, "initial_value"))
return _bufferVector("initial", "initial_displacement", lengthScale);
- else if (0 == strcasecmp(name, "rate-of-change"))
+ else if (0 == strcasecmp(name, "rate_of_change"))
return _bufferVector("rate", "velocity", rateScale);
- else if (0 == strcasecmp(name, "change-in-value"))
+ else if (0 == strcasecmp(name, "change_in_value"))
return _bufferVector("change", "displacement_change", lengthScale);
- else if (0 == strcasecmp(name, "rate-start-time"))
+ else if (0 == strcasecmp(name, "rate_start_time"))
return _bufferScalar("rate time", "velocity_start_time", timeScale);
- else if (0 == strcasecmp(name, "change-start-time"))
+ else if (0 == strcasecmp(name, "change_start_time"))
return _bufferScalar("change time", "change_start_time", timeScale);
else {
std::ostringstream msg;
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/bc/Neumann.cc 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/bc/Neumann.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -203,19 +203,19 @@
assert(_parameters);
assert(name);
- if (0 == strcasecmp(name, "initial-value"))
+ if (0 == strcasecmp(name, "initial_value"))
return _parameters->get("initial");
- else if (0 == strcasecmp(name, "rate-of-change"))
+ else if (0 == strcasecmp(name, "rate_of_change"))
return _parameters->get("rate");
- else if (0 == strcasecmp(name, "change-in-value"))
+ else if (0 == strcasecmp(name, "change_in_value"))
return _parameters->get("change");
- else if (0 == strcasecmp(name, "rate-start-time"))
+ else if (0 == strcasecmp(name, "rate_start_time"))
return _parameters->get("rate time");
- else if (0 == strcasecmp(name, "change-start-time"))
+ else if (0 == strcasecmp(name, "change_start_time"))
return _parameters->get("change time");
else {
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -21,6 +21,7 @@
#include "FaultCohesiveDyn.hh" // implementation of object methods
#include "CohesiveTopology.hh" // USES CohesiveTopology
+#include "TractPerturbation.hh" // HOLDSA TractPerturbation
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
@@ -28,6 +29,7 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
#include "pylith/topology/Jacobian.hh" // USES Jacobian
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/friction/FrictionModel.hh" // USES FrictionModel
@@ -55,6 +57,7 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RealUniformSection RealUniformSection;
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
@@ -66,7 +69,7 @@
pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(void) :
_zeroTolerance(1.0e-10),
_openFreeSurf(true),
- _dbInitialTract(0),
+ _tractPerturbation(0),
_friction(0),
_jacobian(0),
_ksp(0)
@@ -86,7 +89,7 @@
{ // deallocate
FaultCohesiveLagrange::deallocate();
- _dbInitialTract = 0; // :TODO: Use shared pointer
+ _tractPerturbation = 0; // :TODO: Use shared pointer
_friction = 0; // :TODO: Use shared pointer
delete _jacobian; _jacobian = 0;
@@ -99,10 +102,10 @@
// ----------------------------------------------------------------------
// Sets the spatial database for the inital tractions
void
-pylith::faults::FaultCohesiveDyn::dbInitialTract(spatialdata::spatialdb::SpatialDB* db)
-{ // dbInitial
- _dbInitialTract = db;
-} // dbInitial
+pylith::faults::FaultCohesiveDyn::tractPerturbation(TractPerturbation* tract)
+{ // tractPerturbation
+ _tractPerturbation = tract;
+} // tractPerturbation
// ----------------------------------------------------------------------
// Get the friction (constitutive) model.
@@ -149,7 +152,10 @@
FaultCohesiveLagrange::initialize(mesh, upDir);
// Get initial tractions using a spatial database.
- _setupInitialTractions();
+ if (_tractPerturbation) {
+ const topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+ _tractPerturbation->initialize(*_faultMesh, orientation, *_normalizer);
+ } // if
// Setup fault constitutive model.
assert(_friction);
@@ -226,11 +232,22 @@
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());
+ scalar_array tractPerturbVertex(spaceDim);
+ ALE::Obj<RealUniformSection> tractPerturbSection;
+ int tractPerturbIndex = 0;
+ int tractPerturbFiberDim = 0;
+ if (_tractPerturbation) {
+ _tractPerturbation->calculate(t);
+
+ const topology::FieldsNew<topology::SubMesh>* params = _tractPerturbation->parameterFields();
+ assert(params);
+ tractPerturbSection = params->section();
+ assert(!tractPerturbSection.isNull());
+
+ tractPerturbFiberDim = params->fiberDim();
+ tractPerturbIndex = params->sectionIndex("value");
+ const int valueFiberDim = params->sectionFiberDim("value");
+ assert(valueFiberDim == spaceDim);
} // if
const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
@@ -269,13 +286,17 @@
_logger->eventBegin(restrictEvent);
#endif
- // Get initial tractions at fault vertex.
- if (_dbInitialTract) {
- initialTractionsSection->restrictPoint(v_fault,
- &initialTractionsVertex[0],
- initialTractionsVertex.size());
+ // Get prescribed traction perturbation at fault vertex.
+ if (_tractPerturbation) {
+ assert(tractPerturbFiberDim == tractPerturbSection->getFiberDimension(v_fault));
+ const PylithScalar* paramsVertex = tractPerturbSection->restrictPoint(v_fault);
+ assert(paramsVertex);
+
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ tractPerturbVertex[iDim] = paramsVertex[tractPerturbIndex+iDim];
+ } // for
} else {
- initialTractionsVertex = 0.0;
+ tractPerturbVertex = 0.0;
} // if/else
// Get orientation associated with fault vertex.
@@ -350,7 +371,7 @@
//
// Initial (external) tractions oppose (internal) tractions
// associated with Lagrange multiplier.
- residualVertexN = areaVertex * (dispTpdtVertexL - initialTractionsVertex);
+ residualVertexN = areaVertex * (dispTpdtVertexL - tractPerturbVertex);
} else { // opening, normal traction should be zero
if (fabs(tractionNormal) > _zeroTolerance) {
@@ -1475,7 +1496,7 @@
// Get vertex field associated with integrator.
const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::FaultCohesiveDyn::vertexField(const char* name,
- const topology::SolutionFields* fields)
+ const topology::SolutionFields* fields)
{ // vertexField
assert(_faultMesh);
assert(_quadrature);
@@ -1557,17 +1578,6 @@
buffer.copy(dirSection);
return buffer;
- } else if (0 == strcasecmp("initial_traction", name)) {
- assert(_dbInitialTract);
- _allocateBufferVectorField();
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
- topology::Field<topology::SubMesh>& tractions =
- _fields->get("initial traction");
- buffer.copy(tractions);
- FaultCohesiveLagrange::globalToFault(&buffer, orientation);
- return buffer;
-
} else if (0 == strcasecmp("traction", name)) {
assert(fields);
const topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
@@ -1580,6 +1590,19 @@
} else if (_friction->hasPropStateVar(name)) {
return _friction->getField(name);
+ } else if (_tractPerturbation && _tractPerturbation->hasParameter(name)) {
+ const topology::Field<topology::SubMesh>& param = _tractPerturbation->vertexField(name, fields);
+ if (param.vectorFieldType() == topology::FieldBase::VECTOR) {
+ _allocateBufferVectorField();
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ buffer.copy(param);
+ FaultCohesiveLagrange::globalToFault(&buffer, orientation);
+ return buffer;
+ } else {
+ return param;
+ } // if/else
+
} else {
std::ostringstream msg;
msg << "Request for unknown vertex field '" << name << "' for fault '"
@@ -1599,131 +1622,6 @@
} // vertexField
// ----------------------------------------------------------------------
-void
-pylith::faults::FaultCohesiveDyn::_setupInitialTractions(void)
-{ // _setupInitialTractions
- assert(_normalizer);
-
- // If no initial tractions specified, leave method
- if (0 == _dbInitialTract)
- return;
-
- assert(_normalizer);
- const PylithScalar pressureScale = _normalizer->pressureScale();
- const PylithScalar lengthScale = _normalizer->lengthScale();
-
- const int spaceDim = _quadrature->spaceDim();
-
- // Create section to hold initial tractions.
- _fields->add("initial traction", "initial_traction");
- topology::Field<topology::SubMesh>& initialTractions =
- _fields->get("initial traction");
- topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
- initialTractions.cloneSection(dispRel);
- initialTractions.scale(pressureScale);
-
- scalar_array initialTractionsVertex(spaceDim);
- scalar_array initialTractionsVertexGlobal(spaceDim);
- const ALE::Obj<RealSection>& initialTractionsSection =
- initialTractions.section();
- assert(!initialTractionsSection.isNull());
-
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
- assert(!orientationSection.isNull());
-
- const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
- assert(cs);
-
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
-
- scalar_array coordsVertex(spaceDim);
- const ALE::Obj<RealSection>& coordsSection =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordsSection.isNull());
-
-
- assert(_dbInitialTract);
- _dbInitialTract->open();
- switch (spaceDim) { // switch
- case 1: {
- const char* valueNames[] = { "traction-normal" };
- _dbInitialTract->queryVals(valueNames, 1);
- break;
- } // case 1
- case 2: {
- const char* valueNames[] = { "traction-shear", "traction-normal" };
- _dbInitialTract->queryVals(valueNames, 2);
- break;
- } // case 2
- case 3: {
- const char* valueNames[] = { "traction-shear-leftlateral",
- "traction-shear-updip", "traction-normal" };
- _dbInitialTract->queryVals(valueNames, 3);
- break;
- } // case 3
- default:
- std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
- throw std::logic_error("Bad spatial dimension in Neumann.");
- } // switch
-
- const int numVertices = _cohesiveVertices.size();
- for (int iVertex=0; iVertex < numVertices; ++iVertex) {
- const int v_fault = _cohesiveVertices[iVertex].fault;
-
- coordsSection->restrictPoint(v_fault, &coordsVertex[0], coordsVertex.size());
-
- assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
- const PylithScalar* orientationVertex =
- orientationSection->restrictPoint(v_fault);
- assert(orientationVertex);
-
- _normalizer->dimensionalize(&coordsVertex[0], coordsVertex.size(),
- lengthScale);
-
- initialTractionsVertex = 0.0;
- int err = _dbInitialTract->query(&initialTractionsVertex[0],
- initialTractionsVertex.size(),
- &coordsVertex[0], coordsVertex.size(), cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find parameters for physical properties at \n" << "(";
- for (int i = 0; i < spaceDim; ++i)
- msg << " " << coordsVertex[i];
- msg << ") in friction model " << label() << "\n"
- << "using spatial database '" << _dbInitialTract->label() << "'.";
- throw std::runtime_error(msg.str());
- } // if
- _normalizer->nondimensionalize(&initialTractionsVertex[0],
- initialTractionsVertex.size(),
- pressureScale);
-
- // Rotate tractions from fault coordinate system to global
- // coordinate system
- initialTractionsVertexGlobal = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- for (int jDim=0; jDim < spaceDim; ++jDim) {
- initialTractionsVertexGlobal[iDim] +=
- orientationVertex[jDim*spaceDim+iDim] *
- initialTractionsVertex[jDim];
- } // for
- } // for
-
- assert(initialTractionsVertexGlobal.size() ==
- initialTractionsSection->getFiberDimension(v_fault));
- initialTractionsSection->updatePoint(v_fault,
- &initialTractionsVertexGlobal[0]);
- } // for
-
- // Close properties database
- _dbInitialTract->close();
-
- //initialTractions.view("INITIAL TRACTIONS"); // DEBUGGING
-} // _setupInitialTractions
-
-// ----------------------------------------------------------------------
// Compute tractions on fault surface using solution.
void
pylith::faults::FaultCohesiveDyn::_calcTractions(
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-05-10 19:49:30 UTC (rev 20066)
@@ -55,11 +55,11 @@
virtual
void deallocate(void);
- /** Sets the spatial database for the inital tractions.
+ /** Sets the traction perturbation for prescribed tractions.
*
- * @param db spatial database for initial tractions
+ * @param tract Spatial and temporal variation of tractions.
*/
- void dbInitialTract(spatialdata::spatialdb::SpatialDB* db);
+ void tractPerturbation(TractPerturbation* tract);
/** Set the friction (constitutive) model.
*
@@ -152,10 +152,6 @@
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
- /** Get initial tractions using a spatial database.
- */
- void _setupInitialTractions(void);
-
/** Compute tractions on fault surface using solution.
*
* @param tractions Field for tractions.
@@ -283,8 +279,8 @@
/// Minimum resolvable value accounting for roundoff errors
PylithScalar _zeroTolerance;
- /// Database for initial tractions.
- spatialdata::spatialdb::SpatialDB* _dbInitialTract;
+ /// Prescribed traction variation.
+ TractPerturbation* _tractPerturbation;
/// To identify constitutive model
friction::FrictionModel* _friction;
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/Makefile.am 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/Makefile.am 2012-05-10 19:49:30 UTC (rev 20066)
@@ -26,6 +26,14 @@
ConstRateSlipFn.icc \
CohesiveTopology.hh \
EqKinSrc.hh \
+ LiuCosSlipFn.hh \
+ LiuCosSlipFn.icc \
+ SlipTimeFn.hh \
+ StepSlipFn.hh \
+ StepSlipFn.icc \
+ TimeHistorySlipFn.hh \
+ TimeHistorySlipFn.icc \
+ TractPerturbation.hh \
Fault.hh \
Fault.icc \
FaultCohesive.hh \
@@ -35,14 +43,6 @@
FaultCohesiveDyn.hh \
FaultCohesiveKin.hh \
FaultCohesiveImpulses.hh \
- LiuCosSlipFn.hh \
- LiuCosSlipFn.icc \
- SlipTimeFn.hh \
- StepSlipFn.hh \
- StepSlipFn.icc \
- TimeHistorySlipFn.hh \
- TimeHistorySlipFn.icc \
- TractPerturbation.hh \
faultsfwd.hh
noinst_HEADERS = \
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/TractPerturbation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/TractPerturbation.cc 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/TractPerturbation.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -75,6 +75,14 @@
} // label
// ----------------------------------------------------------------------
+// Get parameter fields.
+const pylith::topology::FieldsNew<pylith::topology::SubMesh>*
+pylith::faults::TractPerturbation::parameterFields(void) const
+{ // parameterFields
+ return _parameters;
+} // parameterFields
+
+// ----------------------------------------------------------------------
// Initialize traction perturbation function.
void
pylith::faults::TractPerturbation::initialize(const topology::SubMesh& faultMesh,
@@ -99,16 +107,16 @@
// Create section to hold time dependent values
_parameters->add("value", "traction", spaceDim, topology::FieldBase::VECTOR, pressureScale);
if (_dbInitial)
- _parameters->add("initial", "initial_traction", spaceDim, topology::FieldBase::VECTOR, pressureScale);
+ _parameters->add("initial", "traction_initial", spaceDim, topology::FieldBase::VECTOR, pressureScale);
if (_dbRate) {
_parameters->add("rate", "traction_rate", spaceDim, topology::FieldBase::VECTOR, rateScale);
- _parameters->add("rate time", "traction_rate_time", 1, topology::FieldBase::SCALAR, timeScale);
+ _parameters->add("rate time", "rate_start_time", 1, topology::FieldBase::SCALAR, timeScale);
} // if
if (_dbChange) {
- _parameters->add("change", "change_traction", spaceDim, topology::FieldBase::VECTOR, pressureScale);
- _parameters->add("change time", "change_traction_time", 1, topology::FieldBase::SCALAR, timeScale);
+ _parameters->add("change", "traction_change", spaceDim, topology::FieldBase::VECTOR, pressureScale);
+ _parameters->add("change time", "change_start_time", 1, topology::FieldBase::SCALAR, timeScale);
} // if
- _parameters->allocate(topology::FieldBase::VERTICES_FIELD, 1);
+ _parameters->allocate(topology::FieldBase::VERTICES_FIELD, 0);
const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
assert(!parametersSection.isNull());
@@ -127,8 +135,8 @@
break;
} // case 2
case 3 : {
- const char* valueNames[] = {"traction-shear-horiz",
- "traction-shear-vert",
+ const char* valueNames[] = {"traction-shear-leftlateral",
+ "traction-shear-updip",
"traction-normal"};
_dbInitial->queryVals(valueNames, 3);
break;
@@ -160,8 +168,8 @@
break;
} // case 2
case 3 : {
- const char* valueNames[] = {"traction-rate-shear-horiz",
- "traction-rate-shear-vert",
+ const char* valueNames[] = {"traction-rate-shear-leftlateral",
+ "traction-rate-shear-updip",
"traction-rate-normal"};
_dbRate->queryVals(valueNames, 3);
break;
@@ -196,8 +204,8 @@
break;
} // case 2
case 3 : {
- const char* valueNames[] = {"traction-shear-horiz",
- "traction-shear-vert",
+ const char* valueNames[] = {"traction-shear-leftlateral",
+ "traction-shear-updip",
"traction-normal"};
_dbChange->queryVals(valueNames, 3);
break;
@@ -224,88 +232,149 @@
} // initialize
// ----------------------------------------------------------------------
-// Get traction perturbation on fault surface at time t.
+// Calculate temporal and spatial variation of value over the list of Submesh.
void
-pylith::faults::TractPerturbation::traction(topology::Field<topology::SubMesh>* const tractionField,
- const PylithScalar t)
-{ // traction
- assert(tractionField);
+pylith::faults::TractPerturbation::calculate(const PylithScalar t)
+{ // calculate
assert(_parameters);
- _calculateValue(t);
+ const PylithScalar timeScale = _timeScale;
- const spatialdata::geocoords::CoordSys* cs = tractionField->mesh().coordsys();
- assert(cs);
- const int spaceDim = cs->spaceDim();
-
- // Get vertices in fault mesh
- const ALE::Obj<SieveMesh>& sieveMesh = tractionField->mesh().sieveMesh();
+ // Get vertices.
+ const ALE::Obj<SieveSubMesh>& sieveMesh = _parameters->mesh().sieveMesh();
assert(!sieveMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
assert(!vertices.isNull());
const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
- // Get sections
- scalar_array tractionsVertex(spaceDim);
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
+ const spatialdata::geocoords::CoordSys* cs = _parameters->mesh().coordsys();
+ assert(cs);
+ const int spaceDim = cs->spaceDim();
+
+ const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
assert(!parametersSection.isNull());
const int parametersFiberDim = _parameters->fiberDim();
+ scalar_array parametersVertex(parametersFiberDim);
+
const int valueIndex = _parameters->sectionIndex("value");
const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(valueFiberDim == tractionsVertex.size());
- assert(valueIndex+valueFiberDim <= parametersFiberDim);
+ assert(spaceDim == valueFiberDim);
- const ALE::Obj<RealSection>& tractionSection = tractionField->section();
- assert(!tractionSection.isNull());
+ const int initialIndex = (_dbInitial) ? _parameters->sectionIndex("initial") : -1;
+ const int initialFiberDim = (_dbInitial) ? _parameters->sectionFiberDim("initial") : 0;
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
+ const int rateIndex = (_dbRate) ? _parameters->sectionIndex("rate") : -1;
+ const int rateFiberDim = (_dbRate) ? _parameters->sectionFiberDim("rate") : 0;
+ const int rateTimeIndex = (_dbRate) ? _parameters->sectionIndex("rate time") : -1;
+ const int rateTimeFiberDim = (_dbRate) ? _parameters->sectionFiberDim("rate time") : 0;
+
+ const int changeIndex = (_dbChange) ? _parameters->sectionIndex("change") : -1;
+ const int changeFiberDim = (_dbChange) ? _parameters->sectionFiberDim("change") : 0;
+ const int changeTimeIndex = (_dbChange) ? _parameters->sectionIndex("change time") : -1;
+ const int changeTimeFiberDim = (_dbChange) ? _parameters->sectionFiberDim("change time") : 0;
+
+ for(SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter != verticesEnd; ++v_iter) {
assert(parametersFiberDim == parametersSection->getFiberDimension(*v_iter));
- const PylithScalar* parametersVertex = parametersSection->restrictPoint(*v_iter);
- assert(parametersVertex);
- const PylithScalar* tractionVertex = ¶metersVertex[valueIndex];
- assert(tractionVertex);
+ parametersSection->restrictPoint(*v_iter, ¶metersVertex[0], parametersVertex.size());
+ for (int i=0; i < valueFiberDim; ++i)
+ parametersVertex[valueIndex+i] = 0.0;
- // Update field
- assert(spaceDim == tractionSection->getFiberDimension(*v_iter));
- tractionSection->updateAddPoint(*v_iter, &tractionsVertex[0]);
+ // Contribution from initial value
+ if (_dbInitial) {
+ assert(initialIndex >= 0);
+ assert(initialFiberDim == valueFiberDim);
+ for (int i=0; i < initialFiberDim; ++i)
+ parametersVertex[valueIndex+i] += parametersVertex[initialIndex+i];
+ } // if
+
+ // Contribution from rate of change of value
+ if (_dbRate) {
+ assert(rateIndex >= 0);
+ assert(rateFiberDim == valueFiberDim);
+ assert(rateTimeIndex >= 0);
+ assert(rateTimeFiberDim == 1);
+
+ const PylithScalar tRel = t - parametersVertex[rateTimeIndex];
+ if (tRel > 0.0) // rate of change integrated over time
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ parametersVertex[valueIndex+iDim] += parametersVertex[rateIndex+iDim] * tRel;
+ } // for
+ } // if
+
+ // Contribution from change of value
+ if (_dbChange) {
+ assert(changeIndex >= 0);
+ assert(changeFiberDim == valueFiberDim);
+ assert(changeTimeIndex >= 0);
+ assert(changeTimeFiberDim == 1);
+
+ const PylithScalar tRel = t - parametersVertex[changeTimeIndex];
+ if (tRel >= 0) { // change in value over time
+ PylithScalar scale = 1.0;
+ if (_dbTimeHistory) {
+ PylithScalar tDim = tRel*timeScale;
+ const int err = _dbTimeHistory->query(&scale, tDim);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ parametersVertex[valueIndex+iDim] += parametersVertex[changeIndex+iDim] * scale;
+ } // for
+ } // if
+ } // if
+
+ parametersSection->updatePoint(*v_iter, ¶metersVertex[0]);
} // for
+} // calculate
-} // traction
-
+
// ----------------------------------------------------------------------
-// Get parameter fields.
-const pylith::topology::FieldsNew<pylith::topology::SubMesh>*
-pylith::faults::TractPerturbation::parameterFields(void) const
-{ // parameterFields
- return _parameters;
-} // parameterFields
+// Determine if perturbation has a given parameter.
+bool
+pylith::faults::TractPerturbation::hasParameter(const char* name) const
+{ // hasParameter
+ if (0 == strcasecmp(name, "traction_initial_value"))
+ return (0 != _dbInitial);
+ else if (0 == strcasecmp(name, "traction_rate_of_change"))
+ return (0 != _dbRate);
+ else if (0 == strcasecmp(name, "traction_change_in_value"))
+ return (0 != _dbChange);
+ else if (0 == strcasecmp(name, "traction_rate_start_time"))
+ return (0 != _dbRate);
+ else if (0 == strcasecmp(name, "traction_change_start_time"))
+ return (0 != _dbChange);
+ else
+ return false;
+} // hasParameter
// ----------------------------------------------------------------------
// Get vertex field with traction perturbation information.
const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::TractPerturbation::vertexField(const char* name,
- topology::SolutionFields* const fields)
+ const topology::SolutionFields* const fields)
{ // vertexField
assert(_parameters);
assert(name);
- if (0 == strcasecmp(name, "initial-value"))
+ if (0 == strcasecmp(name, "traction_initial_value"))
return _parameters->get("initial");
- else if (0 == strcasecmp(name, "rate-of-change"))
+ else if (0 == strcasecmp(name, "traction_rate_of_change"))
return _parameters->get("rate");
- else if (0 == strcasecmp(name, "change-in-value"))
+ else if (0 == strcasecmp(name, "traction_change_in_value"))
return _parameters->get("change");
- else if (0 == strcasecmp(name, "rate-start-time"))
+ else if (0 == strcasecmp(name, "traction_rate_start_time"))
return _parameters->get("rate time");
- else if (0 == strcasecmp(name, "change-start-time"))
+ else if (0 == strcasecmp(name, "traction_change_start_time"))
return _parameters->get("change time");
else {
@@ -359,8 +428,8 @@
scalar_array coordsVertexGlobal(spaceDim);
// Get sections.
- const ALE::Obj<RealSection>& coordinates = sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
+ const ALE::Obj<RealSection>& coordsSection = sieveMesh->getRealSection("coordinates");
+ assert(!coordsSection.isNull());
const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
assert(!parametersSection.isNull());
@@ -372,6 +441,8 @@
// Loop over cells in boundary mesh and perform queries.
for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter != verticesEnd; ++v_iter) {
+ assert(spaceDim == coordsSection->getFiberDimension(*v_iter));
+ coordsSection->restrictPoint(*v_iter, &coordsVertexGlobal[0], coordsVertexGlobal.size());
normalizer.dimensionalize(&coordsVertexGlobal[0], coordsVertexGlobal.size(), lengthScale);
valuesVertex = 0.0;
@@ -398,107 +469,4 @@
} // for
} // _queryDB
-// ----------------------------------------------------------------------
-// Calculate temporal and spatial variation of value over the list of Submesh.
-void
-pylith::faults::TractPerturbation::_calculateValue(const PylithScalar t)
-{ // _calculateValue
- assert(_parameters);
-
- const PylithScalar timeScale = _timeScale;
-
- // Get vertices.
- const ALE::Obj<SieveSubMesh>& sieveMesh = _parameters->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
- assert(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
-
- const spatialdata::geocoords::CoordSys* cs = _parameters->mesh().coordsys();
- assert(cs);
- const int spaceDim = cs->spaceDim();
-
- const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- scalar_array parametersVertex(parametersFiberDim);
-
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(spaceDim == valueFiberDim);
-
- const int initialIndex = (_dbInitial) ? _parameters->sectionIndex("initial") : -1;
- const int initialFiberDim = (_dbInitial) ? _parameters->sectionFiberDim("initial") : 0;
-
- const int rateIndex = (_dbRate) ? _parameters->sectionIndex("rate") : -1;
- const int rateFiberDim = (_dbRate) ? _parameters->sectionFiberDim("rate") : 0;
- const int rateTimeIndex = (_dbRate) ? _parameters->sectionIndex("rate time") : -1;
- const int rateTimeFiberDim = (_dbRate) ? _parameters->sectionFiberDim("rate time") : 0;
-
- const int changeIndex = (_dbChange) ? _parameters->sectionIndex("change") : -1;
- const int changeFiberDim = (_dbChange) ? _parameters->sectionFiberDim("change") : 0;
- const int changeTimeIndex = (_dbChange) ? _parameters->sectionIndex("change time") : -1;
- const int changeTimeFiberDim = (_dbChange) ? _parameters->sectionFiberDim("change time") : 0;
-
- for(SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter != verticesEnd; ++v_iter) {
- assert(parametersFiberDim == parametersSection->getFiberDimension(*v_iter));
- parametersSection->restrictPoint(*v_iter, ¶metersVertex[0], parametersVertex.size());
- for (int i=0; i < valueFiberDim; ++i)
- parametersVertex[valueIndex+i] = 0.0;
-
- // Contribution from initial value
- if (_dbInitial) {
- assert(initialIndex >= 0);
- assert(initialFiberDim == valueFiberDim);
- for (int i=0; i < initialFiberDim; ++i)
- parametersVertex[valueIndex+i] += parametersVertex[initialIndex+i];
- } // if
-
- // Contribution from rate of change of value
- if (_dbRate) {
- assert(rateIndex >= 0);
- assert(rateFiberDim == valueFiberDim);
- assert(rateTimeIndex >= 0);
- assert(rateTimeFiberDim == 1);
-
- const PylithScalar tRel = t - parametersVertex[rateTimeIndex];
- if (tRel > 0.0) // rate of change integrated over time
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- parametersVertex[valueIndex+iDim] += parametersVertex[rateIndex+iDim] * tRel;
- } // for
- } // if
-
- // Contribution from change of value
- if (_dbChange) {
- assert(changeIndex >= 0);
- assert(changeFiberDim == valueFiberDim);
- assert(changeTimeIndex >= 0);
- assert(changeTimeFiberDim == 1);
-
- const PylithScalar tRel = t - parametersVertex[changeTimeIndex];
- if (tRel >= 0) { // change in value over time
- PylithScalar scale = 1.0;
- if (_dbTimeHistory) {
- PylithScalar tDim = tRel*timeScale;
- const int err = _dbTimeHistory->query(&scale, tDim);
- if (err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- } // if
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- parametersVertex[valueIndex+iDim] += parametersVertex[changeIndex+iDim] * scale;
- } // for
- } // if
- } // if
-
- parametersSection->updatePoint(*v_iter, ¶metersVertex[0]);
- } // for
-} // _calculateValue
-
-
// End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/TractPerturbation.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/TractPerturbation.hh 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/TractPerturbation.hh 2012-05-10 19:49:30 UTC (rev 20066)
@@ -62,6 +62,12 @@
*/
void label(const char* value);
+ /** Get parameter fields.
+ *
+ * @returns Parameter fields.
+ */
+ const topology::FieldsNew<topology::SubMesh>* parameterFields(void) const;
+
/** Initialize slip time function.
*
* @param faultMesh Finite-element mesh of fault.
@@ -72,20 +78,19 @@
const topology::Field<topology::SubMesh>& faultOrientation,
const spatialdata::units::Nondimensional& normalizer);
- /** Get traction perturbation on fault surface at time t.
+ /** Calculate spatial and temporal variation of value.
*
- * @param tractionField Traction field over fault surface [output].
- * @param t Time t.
+ * @param t Current time.
*/
- void traction(topology::Field<topology::SubMesh>* const tractionField,
- const PylithScalar t);
-
- /** Get parameter fields.
+ void calculate(const PylithScalar t);
+
+ /** Determine if perturbation has a given parameter.
*
- * @returns Parameter fields.
+ * @param name Name of parameter field.
+ * @returns True if perturbation has parameter field, false otherwise.
*/
- const topology::FieldsNew<topology::SubMesh>* parameterFields(void) const;
-
+ bool hasParameter(const char* name) const;
+
/** Get vertex field with traction perturbation information.
*
* @param name Name of field.
@@ -95,7 +100,7 @@
*/
const topology::Field<topology::SubMesh>&
vertexField(const char* name,
- topology::SolutionFields* const fields =0);
+ const topology::SolutionFields* const fields =0);
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -120,12 +125,6 @@
const PylithScalar scale,
const spatialdata::units::Nondimensional& normalizer);
- /** Calculate spatial and temporal variation of value.
- *
- * @param t Current time.
- */
- void _calculateValue(const PylithScalar t);
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/Solver.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -253,7 +253,7 @@
if (formulation->splitFields() &&
formulation->useCustomConstraintPC() &&
- solutionSection->getNumSpaces() > sieveMesh->getDimension()) {
+ numSpaces > spaceDim) {
// We have split fields with a custom constraint preconditioner
// and constraints exist.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/problems/SolverNonlinear.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -248,7 +248,7 @@
SETERRQ(((PetscObject)linesearch)->comm, PETSC_ERR_USER, "SNESLineSearchBT requires a Jacobian matrix");
}
/* precheck */
- ierr = SNESLineSearchPreCheck(linesearch, &changed_y);CHKERRQ(ierr);
+ ierr = SNESLineSearchPreCheck(linesearch,X,Y,&changed_y);CHKERRQ(ierr);
ierr = SNESLineSearchSetSuccess(linesearch, PETSC_TRUE);CHKERRQ(ierr);
ierr = VecNorm(Y, NORM_2, &ynorm);CHKERRQ(ierr);
@@ -460,7 +460,7 @@
}
/* postcheck */
- ierr = SNESLineSearchPostCheck(linesearch, &changed_y, &changed_w);CHKERRQ(ierr);
+ ierr = SNESLineSearchPostCheck(linesearch,X,Y,W,&changed_y,&changed_w);CHKERRQ(ierr);
if (changed_y) {
ierr = VecWAXPY(W,-lambda,Y,X);CHKERRQ(ierr);
if (linesearch->ops->viproject) {
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/bc/TimeDependent.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/bc/TimeDependent.i 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/bc/TimeDependent.i 2012-05-10 19:49:30 UTC (rev 20066)
@@ -24,7 +24,7 @@
namespace pylith {
namespace bc {
- class pylith::bc::TimeDependent
+ class TimeDependent
{ // class TimeDependent
// PUBLIC METHODS /////////////////////////////////////////////////
@@ -81,13 +81,6 @@
virtual
const char* _getLabel(void) const = 0;
- /** Get manager of scales used to nondimensionalize problem.
- *
- * @returns Nondimensionalizer.
- */
- virtual
- const spatialdata::units::Nondimensional& _getNormalizer(void) const = 0;
-
}; // class TimeDependent
} // bc
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/FaultCohesiveDyn.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/FaultCohesiveDyn.i 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/FaultCohesiveDyn.i 2012-05-10 19:49:30 UTC (rev 20066)
@@ -41,11 +41,11 @@
virtual
void deallocate(void);
- /** Sets the spatial database for the inital tractions.
+ /** Sets the traction perturbation for prescribed tractions.
*
- * @param db spatial database for initial tractions
+ * @param tract Spatial and temporal variation of tractions.
*/
- void dbInitialTract(spatialdata::spatialdb::SpatialDB* db);
+ void tractPerturbation(TractPerturbation* tract);
/** Set the friction (constitutive) model.
*
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/Makefile.am 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/Makefile.am 2012-05-10 19:49:30 UTC (rev 20066)
@@ -32,6 +32,7 @@
LiuCosSlipFn.i \
TimeHistorySlipFn.i \
EqKinSrc.i \
+ TractPerturbation.i \
Fault.i \
FaultCohesive.i \
FaultCohesiveLagrange.i \
Copied: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/TractPerturbation.i (from rev 20065, short/3D/PyLith/branches/v1.7-trunk/modulesrc/faults/TractPerturbation.i)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/TractPerturbation.i (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/TractPerturbation.i 2012-05-10 19:49:30 UTC (rev 20066)
@@ -0,0 +1,105 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file modulesrc/faults/TractPerturbation.i
+ *
+ * @brief Python interface to C++ TractPerturbation object.
+ */
+
+namespace pylith {
+ namespace faults {
+
+ class TractPerturbation : public pylith::bc::TimeDependent
+ { // class TractPerturbation
+
+ // PUBLIC METHODS /////////////////////////////////////////////////
+ public :
+
+ /// Default constructor.
+ TractPerturbation(void);
+
+ /// Destructor.
+ virtual
+ ~TractPerturbation(void);
+
+ /// Deallocate PETSc and local data structures.
+ virtual
+ void deallocate(void);
+
+ /** Set label for traction perturbation.
+ *
+ * @param value Label.
+ */
+ void label(const char* value);
+
+ /** Get parameter fields.
+ *
+ * @returns Parameter fields.
+ */
+ const pylith::topology::FieldsNew<pylith::topology::SubMesh>* parameterFields(void) const;
+
+ /** Initialize slip time function.
+ *
+ * @param faultMesh Finite-element mesh of fault.
+ * @param faultOrientation Orientation of fault.
+ * @param normalizer Nondimensionalization of scales.
+ */
+ void initialize(const pylith::topology::SubMesh& faultMesh,
+ const pylith::topology::Field<pylith::topology::SubMesh>& faultOrientation,
+ const spatialdata::units::Nondimensional& normalizer);
+
+ /** Calculate spatial and temporal variation of value.
+ *
+ * @param t Current time.
+ */
+ void calculate(const PylithScalar t);
+
+ /** Determine if perturbation has a given parameter.
+ *
+ * @param name Name of parameter field.
+ * @returns True if perturbation has parameter field, false otherwise.
+ */
+ bool hasParameter(const char* name) const;
+
+ /** Get vertex field with traction perturbation information.
+ *
+ * @param name Name of field.
+ * @param fields Solution fields.
+ *
+ * @returns Traction vector field.
+ */
+ const pylith::topology::Field<pylith::topology::SubMesh>&
+ vertexField(const char* name,
+ pylith::topology::SolutionFields* const fields =0);
+
+ // PROTECTED METHODS //////////////////////////////////////////////
+ protected :
+
+ /** Get label of boundary condition surface.
+ *
+ * @returns Label of surface (from mesh generator).
+ */
+ const char* _getLabel(void) const;
+
+ }; // class TractPerturbation
+
+ } // faults
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/faults.i
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/faults.i 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/modulesrc/faults/faults.i 2012-05-10 19:49:30 UTC (rev 20066)
@@ -28,6 +28,7 @@
#include "pylith/faults/LiuCosSlipFn.hh"
#include "pylith/faults/TimeHistorySlipFn.hh"
#include "pylith/faults/EqKinSrc.hh"
+#include "pylith/faults/TractPerturbation.hh"
#include "pylith/faults/Fault.hh"
#include "pylith/faults/FaultCohesive.hh"
#include "pylith/faults/FaultCohesiveLagrange.hh"
@@ -70,6 +71,7 @@
%include "../topology/SubMesh.i" // ISA Integrator<Quadrature<SubMesh> >
%include "../feassemble/Quadrature.i" // ISA Integrator<Quadrature<SubMesh> >
%include "../feassemble/Integrator.i" // ISA Integrator<Quadrature<SubMesh> >
+%include "../bc/TimeDependent.i" // ISA TimeDependent
// Template instatiation
%template(SubMeshIntegrator) pylith::feassemble::Integrator<pylith::feassemble::Quadrature<pylith::topology::SubMesh > >;
@@ -81,6 +83,7 @@
%include "LiuCosSlipFn.i"
%include "TimeHistorySlipFn.i"
%include "EqKinSrc.i"
+%include "TractPerturbation.i"
%include "Fault.i"
%include "FaultCohesive.i"
%include "FaultCohesiveLagrange.i"
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/Makefile.am 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/Makefile.am 2012-05-10 19:49:30 UTC (rev 20066)
@@ -36,17 +36,18 @@
faults/BruneSlipFn.py \
faults/ConstRateSlipFn.py \
faults/EqKinSrc.py \
+ faults/LiuCosSlipFn.py \
+ faults/SingleRupture.py \
+ faults/SlipTimeFn.py \
+ faults/StepSlipFn.py \
+ faults/TimeHistorySlipFn.py \
+ faults/TractPerturbation.py \
faults/Fault.py \
faults/FaultCohesive.py \
faults/FaultCohesiveKin.py \
faults/FaultCohesiveDyn.py \
faults/FaultCohesiveImpulses.py \
faults/FaultCohesiveTract.py \
- faults/LiuCosSlipFn.py \
- faults/SingleRupture.py \
- faults/SlipTimeFn.py \
- faults/StepSlipFn.py \
- faults/TimeHistorySlipFn.py \
feassemble/__init__.py \
feassemble/Constraint.py \
feassemble/ElasticityExplicit.py \
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/DirichletBoundary.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/DirichletBoundary.py 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/DirichletBoundary.py 2012-05-10 19:49:30 UTC (rev 20066)
@@ -26,6 +26,8 @@
from DirichletBC import DirichletBC
from bc import DirichletBoundary as ModuleDirichletBoundary
+from pylith.utils.NullComponent import NullComponent
+
# DirichletBoundary class
class DirichletBoundary(DirichletBC, ModuleDirichletBoundary):
"""
@@ -69,12 +71,7 @@
self._loggingPrefix = "DiBC "
self.availableFields = \
{'vertex': \
- {'info': ["initial-value",
- "rate-of-change",
- "change-in-value",
- "rate-start-time",
- "change-start-time",
- ],
+ {'info': [],
'data': []},
'cell': \
{'info': [],
@@ -88,6 +85,15 @@
"""
DirichletBC.preinitialize(self, mesh)
self.output.preinitialize(self)
+
+ fields = []
+ if not isinstance(self.inventory.dbInitial, NullComponent):
+ fields += ["initial_value"]
+ if not isinstance(self.inventory.dbRate, NullComponent):
+ fields += ["rate_of_change", "rate_start_time"]
+ if not isinstance(self.inventory.dbChange, NullComponent):
+ fields += ["change_in_value", "change_start_time"]
+ self.availableFields['vertex']['info'] += fields
return
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/Neumann.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/Neumann.py 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/bc/Neumann.py 2012-05-10 19:49:30 UTC (rev 20066)
@@ -27,6 +27,8 @@
from pylith.feassemble.Integrator import Integrator
from bc import Neumann as ModuleNeumann
+from pylith.utils.NullComponent import NullComponent
+
# Neumann class
class Neumann(BoundaryCondition,
TimeDependent,
@@ -68,12 +70,7 @@
{'info': [],
'data': []},
'cell': \
- {'info': ["initial-value",
- "rate-of-change",
- "rate-start-time",
- "change-in-value",
- "change-start-time",
- ],
+ {'info': [],
'data': []}}
return
@@ -88,6 +85,15 @@
self.quadrature(self.bcQuadrature)
self.createSubMesh(mesh)
self.output.preinitialize(self)
+
+ fields = []
+ if not isinstance(self.inventory.dbInitial, NullComponent):
+ fields += ["initial_value"]
+ if not isinstance(self.inventory.dbRate, NullComponent):
+ fields += ["rate_of_change", "rate_start_time"]
+ if not isinstance(self.inventory.dbChange, NullComponent):
+ fields += ["change_in_value", "change_start_time"]
+ self.availableFields['cell']['info'] += fields
return
Modified: short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/FaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/FaultCohesiveDyn.py 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/FaultCohesiveDyn.py 2012-05-10 19:49:30 UTC (rev 20066)
@@ -48,7 +48,7 @@
fault opens.
\b Facilities
- @li \b db_initial_tractions Spatial database for initial tractions.
+ @li \b tract_perturbation Prescribed perturbation in fault tractions.
@li \b friction Fault constitutive model.
@li \b output Output manager associated with fault data.
@@ -68,9 +68,9 @@
"the fault opens, otherwise use initial tractions even when the " \
"fault opens."
- db = pyre.inventory.facility("db_initial_tractions", family="spatial_database",
+ tract = pyre.inventory.facility("traction_perturbation", family="traction_perturbation",
factory=NullComponent)
- db.meta['tip'] = "Spatial database for initial tractions."
+ tract.meta['tip'] = "Prescribed perturbation in fault tractions."
from pylith.friction.StaticFriction import StaticFriction
friction = pyre.inventory.facility("friction", family="friction_model",
@@ -115,6 +115,7 @@
self._info.log("Pre-initializing fault '%s'." % self.label())
FaultCohesive.preinitialize(self, mesh)
Integrator.preinitialize(self, mesh)
+ self.tract.preinitialize(mesh)
ModuleFaultCohesiveDyn.quadrature(self, self.faultQuadrature)
@@ -124,8 +125,9 @@
self.availableFields['vertex']['info'] += ["strike_dir",
"dip_dir"]
- if not isinstance(self.inventory.db, NullComponent):
- self.availableFields['vertex']['info'] += ["initial_traction"]
+ if not isinstance(self.tract, NullComponent):
+ print "HERE",self.tract.availableFields
+ self.availableFields['vertex']['info'] += self.tract.availableFields['vertex']['info']
self.availableFields['vertex']['info'] += \
self.friction.availableFields['vertex']['info']
@@ -209,8 +211,8 @@
Setup members using inventory.
"""
FaultCohesive._configure(self)
- if not isinstance(self.inventory.db, NullComponent):
- ModuleFaultCohesiveDyn.dbInitialTract(self, self.inventory.db)
+ if not isinstance(self.inventory.tract, NullComponent):
+ ModuleFaultCohesiveDyn.tractPerturbation(self, self.inventory.tract)
ModuleFaultCohesiveDyn.frictionModel(self, self.inventory.friction)
ModuleFaultCohesiveDyn.zeroTolerance(self, self.inventory.zeroTolerance)
ModuleFaultCohesiveDyn.openFreeSurf(self, self.inventory.openFreeSurf)
Copied: short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/TractPerturbation.py (from rev 20065, short/3D/PyLith/branches/v1.7-trunk/pylith/faults/TractPerturbation.py)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/TractPerturbation.py (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/pylith/faults/TractPerturbation.py 2012-05-10 19:49:30 UTC (rev 20066)
@@ -0,0 +1,94 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/faults/TractPerturbation.py
+##
+
+## @brief Python object for managing parameters for a kinematic
+## earthquake sources.
+##
+## TractPerturbation is responsible for providing the value of
+## specified traction at time t over a fault surface.
+##
+## Factory: eq_kinematic_src
+
+from pylith.bc.TimeDependent import TimeDependent
+from faults import TractPerturbation as ModuleTractPerturbation
+
+from pylith.utils.NullComponent import NullComponent
+
+# TractPerturbation class
+class TractPerturbation(TimeDependent, ModuleTractPerturbation):
+ """
+ Python object for managing specified tractions on a fault surface.
+
+ Factory: traction_perturbation
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="tractperturbation"):
+ """
+ Constructor.
+ """
+ TimeDependent.__init__(self, name)
+ self._createModuleObj()
+ self._loggingPrefix = "TrPt "
+ return
+
+
+ def preinitialize(self, mesh):
+ """
+ Do pre-initialization setup.
+ """
+ fields = []
+ if not isinstance(self.inventory.dbInitial, NullComponent):
+ fields += ["traction_initial_value"]
+ if not isinstance(self.inventory.dbRate, NullComponent):
+ fields += ["traction_rate_of_change", "traction_rate_start_time"]
+ if not isinstance(self.inventory.dbChange, NullComponent):
+ fields += ["traction_change_in_value", "traction_change_start_time"]
+
+ self.availableFields = \
+ {'vertex': {'info': fields,
+ 'data': []},
+ 'cell': {'info': [],
+ 'data': []}}
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _createModuleObj(self):
+ """
+ Create handle to corresponding C++ object.
+ """
+ ModuleTractPerturbation.__init__(self)
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def traction_perturbation():
+ """
+ Factory associated with TractPerturbation.
+ """
+ return TractPerturbation()
+
+
+# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/dynamic_slipweakening.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/dynamic_slipweakening.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/dynamic_slipweakening.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -24,7 +24,7 @@
[pylithapp.timedependent.interfaces.fault.output]
writer.filename = output/slipweakening-fault.vtk
-vertex_info_fields = [strike_dir,normal_dir,initial_traction,static_coefficient,dynamic_coefficient,slip_weakening_parameter,cohesion]
+vertex_info_fields = [strike_dir,normal_dir,traction_initial_value,static_coefficient,dynamic_coefficient,slip_weakening_parameter,cohesion]
vertex_data_fields = [slip,traction,cumulative_slip,previous_slip]
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/dynamic_timeweakening.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/dynamic_timeweakening.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/dynamic_timeweakening.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -24,7 +24,7 @@
[pylithapp.timedependent.interfaces.fault.output]
writer.filename = output/timeweakening-fault.vtk
-vertex_info_fields = [strike_dir,normal_dir,initial_traction,static_coefficient,dynamic_coefficient,time_weakening_parameter,cohesion]
+vertex_info_fields = [strike_dir,normal_dir,traction_initial_value,static_coefficient,dynamic_coefficient,time_weakening_parameter,cohesion]
vertex_data_fields = [slip,traction,elapsed_time]
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/pylithapp.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/faultstrip/pylithapp.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -102,12 +102,15 @@
quadrature.cell.dimension = 1
quadrature.cell.quad_order = 2
-db_initial_tractions = spatialdata.spatialdb.SimpleDB
-db_initial_tractions.label = Initial fault tractions
-db_initial_tractions.iohandler.filename = tractions.spatialdb
-db_initial_tractions.query_type = linear
+traction_perturbation = pylith.faults.TractPerturbation
+[pylithapp.timedependent.interfaces.fault.traction_perturbation]
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.label = Initial fault tractions
+db_initial.iohandler.filename = tractions.spatialdb
+db_initial.query_type = linear
+
# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/ratestate.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/ratestate.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/ratestate.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -41,11 +41,14 @@
# faults
# ----------------------------------------------------------------------
[pylithapp.timedependent.interfaces.fault]
-db_initial_tractions = spatialdata.spatialdb.UniformDB
-db_initial_tractions.label = Initial fault tractions
-db_initial_tractions.values = [traction-shear,traction-normal]
-db_initial_tractions.data = [0.0*MPa, -10.0*MPa]
+traction_perturbation = pylith.faults.TractPerturbation
+[pylithapp.timedependent.interfaces.fault.traction_perturbation]
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Initial fault tractions
+db_initial.values = [traction-shear,traction-normal]
+db_initial.data = [0.0*MPa, -10.0*MPa]
+
# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/tension.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/tension.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/2d/frictionslide/tension.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -97,11 +97,14 @@
friction.db_properties.values = [static-coefficient,dynamic-coefficient,slip-weakening-parameter,cohesion]
friction.db_properties.data = [0.6,0.2,1.0*mm,0.0*Pa]
-db_initial_tractions = spatialdata.spatialdb.UniformDB
-db_initial_tractions.label = Initial fault tractions
-db_initial_tractions.values = [traction-shear,traction-normal]
-db_initial_tractions.data = [0.0*MPa, -0.0*MPa]
+traction_perturbation = pylith.faults.TractPerturbation
+[pylithapp.timedependent.interfaces.fault.traction_perturbation]
+db_initial = spatialdata.spatialdb.UniformDB
+db_initial.label = Initial fault tractions
+db_initial.values = [traction-shear,traction-normal]
+db_initial.data = [0.0*MPa, -0.0*MPa]
+
[pylithapp.timedependent.interfaces.fault.output]
vertex_data_fields=[slip,slip_rate,traction]
Modified: short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/fieldsplit.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/fieldsplit.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/fieldsplit.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -4,15 +4,11 @@
use_custom_constraint_pc = True
[pylithapp.petsc]
-ksp_gmres_restart = 500
+ksp_gmres_restart = 100
fs_pc_type = fieldsplit
-fs_pc_fieldsplit_real_diagonal =
+fs_pc_fieldsplit_real_diagonal = true
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_1_pc_type = jacobi
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/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/pylithapp.cfg 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/tests/3d/cyclicfriction/pylithapp.cfg 2012-05-10 19:49:30 UTC (rev 20066)
@@ -128,11 +128,14 @@
friction.db_properties.values = [static-coefficient,dynamic-coefficient,slip-weakening-parameter,cohesion]
friction.db_properties.data = [0.6,0.1,1.0*mm,0.0*Pa]
-db_initial_tractions = spatialdata.spatialdb.SimpleDB
-db_initial_tractions.iohandler.filename = initial_tractions.spatialdb
-db_initial_tractions.label = Initial fault tractions
-db_initial_tractions.query_type = linear
+traction_perturbation = pylith.faults.TractPerturbation
+[pylithapp.timedependent.interfaces.fault.traction_perturbation]
+db_initial = spatialdata.spatialdb.SimpleDB
+db_initial.iohandler.filename = initial_tractions.spatialdb
+db_initial.label = Initial fault tractions
+db_initial.query_type = linear
+
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
@@ -150,7 +153,7 @@
writer.filename = output/cyclic-groundsurf.h5
[pylithapp.problem.interfaces.fault.output]
-vertex_info_fields = [initial_traction]
+vertex_info_fields = [traction_initial_value]
writer = pylith.meshio.DataWriterHDF5SubSubMesh
writer.filename = output/cyclic-fault.h5
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/Makefile.am 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/Makefile.am 2012-05-10 19:49:30 UTC (rev 20066)
@@ -57,6 +57,7 @@
TestFaultCohesiveDynQuad4.cc \
TestFaultCohesiveDynTet4.cc \
TestFaultCohesiveDynHex8.cc \
+ TestTractPerturbation.cc \
TestFaultCohesiveImpulses.cc \
TestFaultCohesiveImpulsesCases.cc \
test_faults.cc
@@ -93,6 +94,7 @@
TestFaultCohesiveDynQuad4.hh \
TestFaultCohesiveDynTet4.hh \
TestFaultCohesiveDynHex8.hh \
+ TestTractPerturbation.hh \
TestFaultCohesiveImpulses.hh \
TestFaultCohesiveImpulsesCases.hh
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -21,6 +21,7 @@
#include "TestFaultCohesiveDyn.hh" // Implementation of class methods
#include "pylith/faults/FaultCohesiveDyn.hh" // USES FaultCohesiveDyn
+#include "pylith/faults/TractPerturbation.hh" // USES TractPerturbation
#include "data/CohesiveDynData.hh" // USES CohesiveDynData
@@ -54,7 +55,7 @@
_data = 0;
_quadrature = new feassemble::Quadrature<topology::SubMesh>();
CPPUNIT_ASSERT(0 != _quadrature);
- _dbInitialTract = 0;
+ _tractPerturbation = 0;
_friction = 0;
_dbFriction = 0;
_flipFault = false;
@@ -67,7 +68,7 @@
{ // tearDown
delete _data; _data = 0;
delete _quadrature; _quadrature = 0;
- delete _dbInitialTract; _dbInitialTract = 0;
+ delete _tractPerturbation; _tractPerturbation = 0;
delete _friction; _friction = 0;
delete _dbFriction; _dbFriction = 0;
} // tearDown
@@ -81,19 +82,18 @@
} // testConstructor
// ----------------------------------------------------------------------
-// Test dbInitialTract().
+// Test tractPerturbation().
void
-pylith::faults::TestFaultCohesiveDyn::testDBInitialTract(void)
-{ // testDBInitialTract
+pylith::faults::TestFaultCohesiveDyn::testTractPerturbation(void)
+{ // testTractPerturbation
FaultCohesiveDyn fault;
const std::string& label = "test database";
- spatialdata::spatialdb::SimpleDB db;
- db.label(label.c_str());
- fault.dbInitialTract(&db);
- CPPUNIT_ASSERT(0 != fault._dbInitialTract);
- CPPUNIT_ASSERT_EQUAL(label, std::string(fault._dbInitialTract->label()));
- } // testDBInitialTract
+ TractPerturbation tract;
+ tract.label(label.c_str());
+ fault.tractPerturbation(&tract);
+ CPPUNIT_ASSERT(fault._tractPerturbation);
+ } // testTractPerturbation
// ----------------------------------------------------------------------
// Test zeroTolerance().
@@ -128,7 +128,7 @@
void
pylith::faults::TestFaultCohesiveDyn::testInitialize(void)
{ // testInitialize
- CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(_data);
topology::Mesh mesh;
FaultCohesiveDyn fault;
@@ -181,12 +181,15 @@
} // for
} // for
- // Initial tractions
- if (0 != fault._dbInitialTract) {
- //fault._fields->get("initial traction").view("INITIAL TRACTIONS"); // DEBUGGING
+ // Prescribed traction perturbation
+ if (fault._tractPerturbation) {
+ // :KLUDGE: Only check initial value
const ALE::Obj<RealSection>& initialTractionsSection =
- fault._fields->get("initial traction").section();
+ fault.vertexField("traction_initial_value").section();
CPPUNIT_ASSERT(!initialTractionsSection.isNull());
+
+ //initialTractionsSection->view("INITIAL TRACTIONS"); // DEBUGGING
+
const int spaceDim = _data->spaceDim;
iVertex = 0;
for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin;
@@ -719,15 +722,17 @@
_data->quadWts, _data->numQuadPts,
_data->spaceDim);
- // Setup initial tractions
- spatialdata::spatialdb::SimpleDB* db =
- new spatialdata::spatialdb::SimpleDB("initial tractions");
- CPPUNIT_ASSERT(0 != db);
+ // Setup prescribed traction perturbation
+ delete _tractPerturbation; _tractPerturbation = new TractPerturbation();
+ _tractPerturbation->label("traction perturbation");
+ spatialdata::spatialdb::SimpleDB* db = new spatialdata::spatialdb::SimpleDB("initial tractions");
+ CPPUNIT_ASSERT(db);
spatialdata::spatialdb::SimpleIOAscii ioInitialTract;
ioInitialTract.filename(_data->initialTractFilename);
db->ioHandler(&ioInitialTract);
delete _dbInitialTract; _dbInitialTract = db;
- fault->dbInitialTract(db);
+ _tractPerturbation->dbInitial(db);
+ fault->tractPerturbation(_tractPerturbation);
// Setup friction
spatialdata::spatialdb::SimpleDB* dbFriction =
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2012-05-10 19:49:30 UTC (rev 20066)
@@ -52,7 +52,7 @@
CPPUNIT_TEST_SUITE( TestFaultCohesiveDyn );
CPPUNIT_TEST( testConstructor );
- CPPUNIT_TEST( testDBInitialTract );
+ CPPUNIT_TEST( testTractPerturbation );
CPPUNIT_TEST( testZeroTolerance );
CPPUNIT_TEST( testOpenFreeSurf );
@@ -71,6 +71,7 @@
CohesiveDynData* _data; ///< Data for testing
feassemble::Quadrature<topology::SubMesh>* _quadrature; ///< Fault quad.
+ TractPerturbation* _tractPerturbation; ///< Initial tractions.
spatialdata::spatialdb::SpatialDB* _dbInitialTract; ///< Initial tractions.
friction::FrictionModel* _friction; ///< Friction model
spatialdata::spatialdb::SpatialDB* _dbFriction; ///< Friction parameters.
@@ -88,8 +89,8 @@
/// Test constructor.
void testConstructor(void);
- /// Test dbInitialTract().
- void testDBInitialTract(void);
+ /// Test tractPerturbation().
+ void testTractPerturbation(void);
/// Test zeroTolerance().
void testZeroTolerance(void);
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestTractPerturbation.cc (from rev 20065, short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestTractPerturbation.cc)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestTractPerturbation.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestTractPerturbation.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -0,0 +1,348 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestTractPerturbation.hh" // Implementation of class methods
+
+#include "pylith/faults/TractPerturbation.hh" // USES TractPerturbation
+
+#include "pylith/faults/CohesiveTopology.hh" // USES CohesiveTopology
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+#include "spatialdata/spatialdb/SimpleDB.hh" // USES SimpleDB
+#include "spatialdata/spatialdb/SimpleIOAscii.hh" // USES SimpleIOAscii
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::faults::TestTractPerturbation );
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RealUniformSection RealUniformSection;
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::faults::TestTractPerturbation::testConstructor(void)
+{ // testConstructor
+ TractPerturbation eqsrc;
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test label().
+void
+pylith::faults::TestTractPerturbation::testLabel(void)
+{ // testLabel
+ const std::string& label = "nucleation";
+
+ TractPerturbation tract;
+ tract.label(label.c_str());
+ CPPUNIT_ASSERT_EQUAL(label, tract._label);
+} // testLabel
+
+// ----------------------------------------------------------------------
+// Test hasParameter().
+void
+pylith::faults::TestTractPerturbation::testHasParameter(void)
+{ // testHasParameter
+ spatialdata::spatialdb::SimpleDB db;
+
+ TractPerturbation tract;
+
+ // no values
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_initial_value"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_rate_of_change"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_change_in_value"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_rate_start_time"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_change_start_time"));
+
+ // initial value
+ tract.dbInitial(&db);
+ CPPUNIT_ASSERT_EQUAL(true, tract.hasParameter("traction_initial_value"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_rate_of_change"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_change_in_value"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_rate_start_time"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_change_start_time"));
+
+ // change value
+ tract.dbChange(&db);
+ CPPUNIT_ASSERT_EQUAL(true, tract.hasParameter("traction_initial_value"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_rate_of_change"));
+ CPPUNIT_ASSERT_EQUAL(true, tract.hasParameter("traction_change_in_value"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_rate_start_time"));
+ CPPUNIT_ASSERT_EQUAL(true, tract.hasParameter("traction_change_start_time"));
+
+ // rate value, remove change
+ tract.dbRate(&db);
+ tract.dbChange(0);
+ CPPUNIT_ASSERT_EQUAL(true, tract.hasParameter("traction_initial_value"));
+ CPPUNIT_ASSERT_EQUAL(true, tract.hasParameter("traction_rate_of_change"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_change_in_value"));
+ CPPUNIT_ASSERT_EQUAL(true, tract.hasParameter("traction_rate_start_time"));
+ CPPUNIT_ASSERT_EQUAL(false, tract.hasParameter("traction_change_start_time"));
+} // testHasParameter
+
+// ----------------------------------------------------------------------
+// Test initialize() using 2-D mesh.
+void
+pylith::faults::TestTractPerturbation::testInitialize(void)
+{ // testInitialize
+ topology::Mesh mesh;
+ topology::SubMesh faultMesh;
+ TractPerturbation tract;
+ _initialize(&mesh, &faultMesh, &tract);
+
+ // Rely on testTraction() for verification of results.
+} // testInitialize
+
+// ----------------------------------------------------------------------
+// Test calculate() using 2-D mesh().
+void
+pylith::faults::TestTractPerturbation::testCalculate(void)
+{ // testCalculate
+ const PylithScalar tractionE[4] = {
+ -1.0*(-2.0+1.0), -1.0*(1.0-0.5), // initial + change
+ -1.0*(-2.1), -1.0*(1.1), // initial
+ };
+
+ topology::Mesh mesh;
+ topology::SubMesh faultMesh;
+ TractPerturbation tract;
+ _initialize(&mesh, &faultMesh, &tract);
+
+ const PylithScalar t = 2.134;
+ tract.calculate(t);
+
+
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
+ CPPUNIT_ASSERT(cs);
+ const int spaceDim = cs->spaceDim();
+
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = faultSieveMesh->depthStratum(0);
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ //traction.view("TRACTION"); // DEBUGGING
+
+ const PylithScalar tolerance = 1.0e-06;
+ int iPoint = 0;
+
+ CPPUNIT_ASSERT(tract._parameters);
+ const ALE::Obj<RealUniformSection>& paramsSection = tract._parameters->section();
+ CPPUNIT_ASSERT(!paramsSection.isNull());
+ const int fiberDim = tract._parameters->fiberDim();
+ const int valueIndex = tract._parameters->sectionIndex("value");
+ const int valueFiberDim = tract._parameters->sectionFiberDim("value");
+ CPPUNIT_ASSERT_EQUAL(spaceDim, valueFiberDim);
+ CPPUNIT_ASSERT(valueIndex + valueFiberDim < fiberDim);
+
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin(); v_iter != verticesEnd; ++v_iter, ++iPoint) {
+ CPPUNIT_ASSERT_EQUAL(fiberDim, paramsSection->getFiberDimension(*v_iter));
+ const PylithScalar* vals = paramsSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(vals);
+
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const PylithScalar valueE = tractionE[iPoint*spaceDim+iDim];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, vals[valueIndex+iDim], tolerance);
+ } // for
+ } // for
+} // testCalculate
+
+// ----------------------------------------------------------------------
+// Test parameterFields() using 2-D mesh.
+void
+pylith::faults::TestTractPerturbation::testParameterFields(void)
+{ // testParameterFields
+ const int fiberDimE = 7;
+ const PylithScalar parametersE[2*fiberDimE] = {
+ 0.0, 0.0, 2.0, -1.0, -1.0, 0.5, 1.5,
+ 0.0, 0.0, 2.1, -1.1, -0.8, 0.7, 2.5,
+ };
+
+ topology::Mesh mesh;
+ topology::SubMesh faultMesh;
+ TractPerturbation tract;
+ _initialize(&mesh, &faultMesh, &tract);
+
+ const topology::FieldsNew<topology::SubMesh>* parameters = tract.parameterFields();
+ const ALE::Obj<RealUniformSection>& parametersSection = parameters->section();
+ CPPUNIT_ASSERT(!parametersSection.isNull());
+
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = faultSieveMesh->depthStratum(0);
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ const PylithScalar tolerance = 1.0e-06;
+ int iPoint = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin(); v_iter != verticesEnd; ++v_iter, ++iPoint) {
+ const int fiberDim = parametersSection->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const PylithScalar* vals = parametersSection->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(vals);
+
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ const PylithScalar valueE = parametersE[iPoint*fiberDim+iDim];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, vals[iDim], tolerance);
+ } // for
+ } // for
+} // testParameterFields
+
+// ----------------------------------------------------------------------
+// Test vertexField() using 2-D mesh.
+void
+pylith::faults::TestTractPerturbation::testVertexField(void)
+{ // testVertexField
+ const int fiberDimE = 1;
+ const PylithScalar fieldE[2*fiberDimE] = {
+ 1.5,
+ 2.5,
+ };
+ const char* label = "traction_change_start_time";
+
+ topology::Mesh mesh;
+ topology::SubMesh faultMesh;
+ TractPerturbation tract;
+ _initialize(&mesh, &faultMesh, &tract);
+
+ const topology::Field<topology::SubMesh>& field = tract.vertexField(label);
+ const ALE::Obj<RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = faultSieveMesh->depthStratum(0);
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ const PylithScalar tolerance = 1.0e-06;
+ int iPoint = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin(); v_iter != verticesEnd; ++v_iter, ++iPoint) {
+ const int fiberDim = section->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const PylithScalar* vals = section->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(vals);
+
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ const PylithScalar valueE = fieldE[iPoint*fiberDim+iDim];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, vals[iDim], tolerance);
+ } // for
+ } // for
+} // testVertexField
+
+// ----------------------------------------------------------------------
+// Initialize TractPerturbation.
+void
+pylith::faults::TestTractPerturbation::_initialize(topology::Mesh* mesh,
+ topology::SubMesh* faultMesh,
+ TractPerturbation* tract)
+{ // _initialize
+ CPPUNIT_ASSERT(mesh);
+ CPPUNIT_ASSERT(faultMesh);
+ CPPUNIT_ASSERT(tract);
+
+ const char* meshFilename = "data/tri3.mesh";
+ const char* faultLabel = "fault";
+ const int faultId = 2;
+ const char* initialFilename = "data/tri3_initialtract.spatialdb";
+ const char* changeFilename = "data/tri3_changetract.spatialdb";
+ const PylithScalar orientationVertex[4] = {
+ 0.0, -1.0, -1.0, 0.0,
+ };
+
+ meshio::MeshIOAscii meshIO;
+ meshIO.filename(meshFilename);
+ meshIO.debug(false);
+ meshIO.interpolate(false);
+ meshIO.read(mesh);
+
+ // Set up coordinates
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(mesh->dimension());
+ cs.initialize();
+ mesh->coordsys(&cs);
+ const int spaceDim = cs.spaceDim();
+
+ // Create fault mesh
+ int firstFaultVertex = 0;
+ int firstLagrangeVertex = mesh->sieveMesh()->getIntSection(faultLabel)->size();
+ int firstFaultCell = mesh->sieveMesh()->getIntSection(faultLabel)->size();
+ const bool useLagrangeConstraints = true;
+ if (useLagrangeConstraints) {
+ firstFaultCell += mesh->sieveMesh()->getIntSection(faultLabel)->size();
+ } // if
+ ALE::Obj<SieveFlexMesh> faultBoundary = 0;
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ CohesiveTopology::createFault(faultMesh, faultBoundary,
+ *mesh, sieveMesh->getIntSection(faultLabel));
+ CohesiveTopology::create(mesh, *faultMesh, faultBoundary,
+ sieveMesh->getIntSection(faultLabel),
+ faultId,
+ firstFaultVertex, firstLagrangeVertex, firstFaultCell,
+ useLagrangeConstraints);
+ // Need to copy coordinates from mesh to fault mesh since we are
+ // using create() instead of createParallel().
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
+ CPPUNIT_ASSERT(!faultSieveMesh.isNull());
+ faultSieveMesh->setRealSection("coordinates",
+ sieveMesh->getRealSection("coordinates"));
+
+ // Setup databases
+ spatialdata::spatialdb::SimpleDB dbInitial("initial traction");
+ spatialdata::spatialdb::SimpleIOAscii ioInitial;
+ ioInitial.filename(initialFilename);
+ dbInitial.ioHandler(&ioInitial);
+
+ spatialdata::spatialdb::SimpleDB dbChange("traction change");
+ spatialdata::spatialdb::SimpleIOAscii ioChange;
+ ioChange.filename(changeFilename);
+ dbChange.ioHandler(&ioChange);
+
+ // Setup fault orientation
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = faultSieveMesh->depthStratum(0);
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ topology::Field<topology::SubMesh> faultOrientation(*faultMesh);
+ faultOrientation.newSection(vertices, spaceDim*spaceDim);
+ faultOrientation.allocate();
+ const ALE::Obj<RealSection>& orientationSection = faultOrientation.section();
+ CPPUNIT_ASSERT(!orientationSection.isNull());
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin(); v_iter != verticesEnd; ++v_iter) {
+ CPPUNIT_ASSERT_EQUAL(spaceDim*spaceDim, orientationSection->getFiberDimension(*v_iter));
+ orientationSection->updatePoint(*v_iter, orientationVertex);
+ } // for
+
+ spatialdata::units::Nondimensional normalizer;
+
+ // setup TractPerturbation
+ tract->dbInitial(&dbInitial);
+ tract->dbChange(&dbChange);
+
+ tract->initialize(*faultMesh, faultOrientation, normalizer);
+} // _initialize
+
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestTractPerturbation.hh (from rev 20065, short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/TestTractPerturbation.hh)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestTractPerturbation.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/TestTractPerturbation.hh 2012-05-10 19:49:30 UTC (rev 20066)
@@ -0,0 +1,102 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010-2012 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/faults/TestTractPerturbation.hh
+ *
+ * @brief C++ TestTractPerturbation object
+ *
+ * C++ unit testing for TractPerturbation.
+ */
+
+#if !defined(pylith_faults_testtractperturbation_hh)
+#define pylith_faults_testtractperturbation_hh
+
+#include "pylith/faults/faultsfwd.hh" // USES TractPerturbation
+#include "pylith/topology/topologyfwd.hh" // USES Mesh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace faults {
+ class TestTractPerturbation;
+ } // faults
+} // pylith
+
+/// C++ unit testing for TractPerturbation
+class pylith::faults::TestTractPerturbation : public CppUnit::TestFixture
+{ // class TestTractPerturbation
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestTractPerturbation );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testLabel );
+ CPPUNIT_TEST( testHasParameter );
+ CPPUNIT_TEST( testInitialize );
+ CPPUNIT_TEST( testCalculate );
+ CPPUNIT_TEST( testParameterFields );
+ CPPUNIT_TEST( testVertexField );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test label().
+ void testLabel(void);
+
+ /// Test hasParameter().
+ void testHasParameter(void);
+
+ /// Test initialize() with 2-D mesh.
+ void testInitialize(void);
+
+ /// Test calculate() with 2-D mesh.
+ void testCalculate(void);
+
+ /// Test parameterFields() with 2-D mesh.
+ void testParameterFields(void);
+
+ /// Test VertexField() with 2-D mesh.
+ void testVertexField(void);
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Initialize TractPerturbation.
+ *
+ * @param mesh Finite-element mesh of domain.
+ * @param faultMesh Finite-element mesh of fault.
+ * @param tract Traction perturbation.
+ */
+ static
+ void _initialize(topology::Mesh* mesh,
+ topology::SubMesh* faultMesh,
+ TractPerturbation* tract);
+
+}; // class TestTractPerturbation
+
+#endif // pylith_faults_testtractperturbation_hh
+
+
+// End of file
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-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataHex8.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -1340,10 +1340,11 @@
};
const PylithScalar pylith::faults::CohesiveDynDataHex8::_initialTractions[] = {
- +3.0, -1.0, +2.0,
- +3.1, -1.1, +2.1,
- +3.2, -1.2, +2.2,
- +3.3, -1.3, +2.3,
+ // Fault coordinate frame
+ +1.0, +2.0, -3.0,
+ +1.1, +2.1, -3.1,
+ +1.2, +2.2, -3.2,
+ +1.3, +2.3, -3.3,
};
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-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -281,8 +281,9 @@
};
const PylithScalar pylith::faults::CohesiveDynDataQuad4::_initialTractions[] = {
- 2.0, -1.0,
- 2.1, -1.1,
+ // Fault coordinate frame
+ 1.0, -2.0,
+ 1.1, -2.1,
};
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-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -485,9 +485,10 @@
};
const PylithScalar pylith::faults::CohesiveDynDataTet4::_initialTractions[] = {
- +3.0, -1.0, +2.0,
- +3.1, -1.1, +2.1,
- +3.2, -1.2, +2.2,
+ // Fault coordinate frame
+ +1.0, +2.0, -3.0,
+ +1.1, +2.1, -3.1,
+ +1.2, +2.2, -3.2,
};
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-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -251,8 +251,9 @@
};
const PylithScalar pylith::faults::CohesiveDynDataTri3::_initialTractions[] = {
- 2.0, -1.0,
- 2.1, -1.1,
+ // Fault coordinate frame
+ 1.0, -2.0,
+ 1.1, -2.1,
};
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-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc 2012-05-10 19:49:30 UTC (rev 20066)
@@ -437,9 +437,10 @@
};
const PylithScalar pylith::faults::CohesiveDynDataTri3d::_initialTractions[] = {
- 3.0*0.70710678118654757, 0.70710678118654757,
- 2.1, -1.1,
- 1.2, 2.2,
+ // Fault coordinate frame
+ 1.0, -2.0,
+ 1.1, -2.1,
+ 1.2, -2.2,
};
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/Makefile.am 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/Makefile.am 2012-05-10 19:49:30 UTC (rev 20066)
@@ -40,6 +40,7 @@
tri3d_sliptime.spatialdb \
tri3d_risetime.spatialdb \
tri3_initialtract.spatialdb \
+ tri3_changetract.spatialdb \
tri3d_initialtract.spatialdb \
tri3_impulses.spatialdb \
quad4.mesh \
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/tri3_changetract.spatialdb (from rev 20065, short/3D/PyLith/branches/v1.7-trunk/unittests/libtests/faults/data/tri3_changetract.spatialdb)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/tri3_changetract.spatialdb (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/libtests/faults/data/tri3_changetract.spatialdb 2012-05-10 19:49:30 UTC (rev 20066)
@@ -0,0 +1,15 @@
+#SPATIAL.ascii 1
+SimpleDB {
+ num-values = 3
+ value-names = traction-shear traction-normal change-start-time
+ value-units = Pa Pa s
+ num-locs = 2
+ data-dim = 1
+ space-dim = 2
+ cs-data = cartesian {
+ to-meters = 1.0
+ space-dim = 2
+ }
+}
+0.0 +1.0 -0.5 +1.0 1.5
+0.0 -1.0 -0.7 +0.8 2.5
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/TestFaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/TestFaultCohesiveDyn.py 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/TestFaultCohesiveDyn.py 2012-05-10 19:49:30 UTC (rev 20066)
@@ -280,7 +280,7 @@
quadrature.inventory.cell = cell
quadrature._configure()
- # Setup earthquake source
+ # Setup rupture info
from spatialdata.spatialdb.SimpleDB import SimpleDB
from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
ioTractions = SimpleIOAscii()
@@ -290,7 +290,11 @@
dbTractions.inventory.iohandler = ioTractions
dbTractions.inventory.label = "initial tractions"
dbTractions._configure()
-
+ from pylith.faults.TractPerturbation import TractPerturbation
+ tract = TractPerturbation()
+ tract.inventory.dbInitial = dbTractions
+ tract._configure()
+
ioFriction = SimpleIOAscii()
ioFriction.inventory.filename = "data/tri3_staticfriction.spatialdb"
ioFriction._configure()
@@ -313,6 +317,7 @@
fault.inventory.faultLabel = "fault"
fault.inventory.upDir = [0, 0, 1]
fault.inventory.faultQuadrature = quadrature
+ fault.inventory.tract = tract
fault.inventory.friction = friction
fault._configure()
Copied: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/TestTractPerturbation.py (from rev 20065, short/3D/PyLith/branches/v1.7-trunk/unittests/pytests/faults/TestTractPerturbation.py)
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/TestTractPerturbation.py (rev 0)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/TestTractPerturbation.py 2012-05-10 19:49:30 UTC (rev 20066)
@@ -0,0 +1,81 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2012 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/faults/TestTractPerturbation.py
+
+## @brief Unit testing of TractPerturbation object.
+
+import unittest
+
+from pylith.faults.TractPerturbation import TractPerturbation
+
+# ----------------------------------------------------------------------
+class TestTractPerturbation(unittest.TestCase):
+ """
+ Unit testing of TractPerturbation object.
+ """
+
+ def test_constructor(self):
+ """
+ Test constructor.
+ """
+ tract = TractPerturbation()
+ return
+
+
+ def test_configure(self):
+ """
+ Test initialize().
+ """
+ from spatialdata.spatialdb.SimpleDB import SimpleDB
+ from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+ from pyre.units.time import second
+
+ ioInitial = SimpleIOAscii()
+ ioInitial.inventory.filename = "tri3_initialtractions.spatialdb"
+ ioInitial._configure()
+ dbInitial = SimpleDB()
+ dbInitial.inventory.iohandler = ioInitial
+ dbInitial.inventory.label = "initial tractions"
+ dbInitial._configure()
+
+ ioChange = SimpleIOAscii()
+ ioChange.inventory.filename = "tri3_changetractions.spatialdb"
+ ioChange._configure()
+ dbChange = SimpleDB()
+ dbChange.inventory.iohandler = ioChange
+ dbChange.inventory.label = "traction change"
+ dbChange._configure()
+
+ tract = TractPerturbation()
+ tract.inventory.dbInitial = dbInitial
+ tract.inventory.dbChange = dbChange
+ tract._configure()
+ return
+
+
+ def test_factory(self):
+ """
+ Test factory method.
+ """
+ from pylith.faults.TractPerturbation import traction_perturbation
+ fn = traction_perturbation()
+ return
+
+
+# End of file
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/data/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/data/Makefile.am 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/data/Makefile.am 2012-05-10 19:49:30 UTC (rev 20066)
@@ -23,6 +23,7 @@
tri3_sliptime.spatialdb \
tri3_peakrate.spatialdb \
tri3_initialtractions.spatialdb \
+ tri3_changetractions.spatialdb \
tri3_staticfriction.spatialdb \
tri3_impulses.spatialdb \
slipfn.timedb
Modified: short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/testfaults.py
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/testfaults.py 2012-05-10 19:21:33 UTC (rev 20065)
+++ short/3D/PyLith/branches/pylith-scecdynrup/unittests/pytests/faults/testfaults.py 2012-05-10 19:49:30 UTC (rev 20066)
@@ -80,6 +80,9 @@
from TestEqKinSrc import TestEqKinSrc
suite.addTest(unittest.makeSuite(TestEqKinSrc))
+ from TestTractPerturbation import TestTractPerturbation
+ suite.addTest(unittest.makeSuite(TestTractPerturbation))
+
from TestFaultCohesiveKin import TestFaultCohesiveKin
suite.addTest(unittest.makeSuite(TestFaultCohesiveKin))
More information about the CIG-COMMITS
mailing list