[cig-commits] r20044 - in short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith: . bc faults
brad at geodynamics.org
brad at geodynamics.org
Fri May 4 17:13:07 PDT 2012
Author: brad
Date: 2012-05-04 17:13:07 -0700 (Fri, 04 May 2012)
New Revision: 20044
Removed:
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Nucleator.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Nucleator.hh
Modified:
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/Neumann.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/Neumann.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/TimeDependent.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/TimeDependentPoints.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveKin.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.hh
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Makefile.am
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/TractPerturbation.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/TractPerturbation.hh
Log:
Switched TractPerturbation to time-dependent formulation used in BC.
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/Makefile.am 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/Makefile.am 2012-05-05 00:13:07 UTC (rev 20044)
@@ -69,9 +69,7 @@
faults/BruneSlipFn.cc \
faults/TimeHistorySlipFn.cc \
faults/LiuCosSlipFn.cc \
- faults/Nucleator.cc \
faults/TractPerturbation.cc \
- faults/StaticPerturbation.cc \
feassemble/CellGeometry.cc \
feassemble/Constraint.cc \
feassemble/GeometryPoint1D.cc \
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/Neumann.cc 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/Neumann.cc 2012-05-05 00:13:07 UTC (rev 20044)
@@ -23,7 +23,6 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/FieldsNew.hh" // HOLDSA FieldsNew
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/Fields.hh" // USES Fields
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -63,6 +62,7 @@
void
pylith::bc::Neumann::deallocate(void)
{ // deallocate
+ TimeDependent::deallocate();
} // deallocate
// ----------------------------------------------------------------------
@@ -259,7 +259,7 @@
_parameters->add("rate", "traction_rate",
numQuadPts*spaceDim, topology::FieldBase::MULTI_VECTOR,
rateScale);
- _parameters->add("rate time", "traction_rate__time",
+ _parameters->add("rate time", "traction_rate_time",
numQuadPts, topology::FieldBase::MULTI_SCALAR,
timeScale);
} // if
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/Neumann.hh 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/Neumann.hh 2012-05-05 00:13:07 UTC (rev 20044)
@@ -76,9 +76,7 @@
/** Get cell field with BC information.
*
- * @param fieldType Type of field.
* @param name Name of field.
- * @param mesh Finite-element mesh.
* @param fields Solution fields.
*
* @returns Traction vector at integration points.
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/TimeDependent.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/TimeDependent.hh 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/TimeDependent.hh 2012-05-05 00:13:07 UTC (rev 20044)
@@ -93,13 +93,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;
-
// PROTECTED MEMBERS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/TimeDependentPoints.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/TimeDependentPoints.hh 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/bc/TimeDependentPoints.hh 2012-05-05 00:13:07 UTC (rev 20044)
@@ -73,6 +73,13 @@
*/
const char* _getLabel(void) const;
+ /** Get manager of scales used to nondimensionalize problem.
+ *
+ * @returns Nondimensionalizer.
+ */
+ virtual
+ const spatialdata::units::Nondimensional& _getNormalizer(void) const = 0;
+
/** Query databases for time dependent parameters.
*
* @param mesh Finite-element mesh.
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-05-05 00:13:07 UTC (rev 20044)
@@ -1484,6 +1484,8 @@
const int cohesiveDim = _faultMesh->dimension();
const int spaceDim = _quadrature->spaceDim();
+ const topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+
PylithScalar scale = 0.0;
int fiberDim = 0;
if (0 == strcasecmp("slip", name)) {
@@ -1494,7 +1496,7 @@
_fields->get("buffer (vector)");
buffer.copy(dispRel);
buffer.label("slip");
- _globalToFault(&buffer);
+ FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
} else if (0 == strcasecmp("slip_rate", name)) {
@@ -1505,7 +1507,7 @@
_fields->get("buffer (vector)");
buffer.copy(velRel);
buffer.label("slip_rate");
- _globalToFault(&buffer);
+ FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
@@ -1561,7 +1563,7 @@
topology::Field<topology::SubMesh>& tractions =
_fields->get("initial traction");
buffer.copy(tractions);
- _globalToFault(&buffer);
+ FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
} else if (0 == strcasecmp("traction", name)) {
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc 2012-05-05 00:13:07 UTC (rev 20044)
@@ -168,7 +168,8 @@
// Transform slip from local (fault) coordinate system to relative
// displacement field in global coordinate system
- _faultToGlobal(&dispRel);
+ const topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+ FaultCohesiveLagrange::faultToGlobal(&dispRel, orientation);
_logger->eventEnd(setupEvent);
@@ -190,6 +191,8 @@
const int cohesiveDim = _faultMesh->dimension();
const int spaceDim = _quadrature->spaceDim();
+ const topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+
PylithScalar scale = 0.0;
int fiberDim = 0;
if (0 == strcasecmp("slip", name)) {
@@ -200,7 +203,7 @@
_fields->get("buffer (vector)");
buffer.copy(dispRel);
buffer.label("slip");
- _globalToFault(&buffer);
+ FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveKin.cc 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveKin.cc 2012-05-05 00:13:07 UTC (rev 20044)
@@ -142,7 +142,8 @@
// Transform slip from local (fault) coordinate system to relative
// displacement field in global coordinate system
- _faultToGlobal(&dispRel);
+ const topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+ FaultCohesiveLagrange::faultToGlobal(&dispRel, orientation);
_logger->eventEnd(setupEvent);
@@ -164,6 +165,8 @@
const int cohesiveDim = _faultMesh->dimension();
const int spaceDim = _quadrature->spaceDim();
+ const topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+
const int slipStrLen = strlen("final_slip");
const int timeStrLen = strlen("slip_time");
@@ -177,7 +180,7 @@
_fields->get("buffer (vector)");
buffer.copy(dispRel);
buffer.label("slip");
- _globalToFault(&buffer);
+ FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-05-05 00:13:07 UTC (rev 20044)
@@ -1163,6 +1163,118 @@
} // initializeLogger
// ----------------------------------------------------------------------
+// Transform field from local (fault) coordinate system to
+// global coordinate system.
+void
+pylith::faults::FaultCohesiveLagrange::faultToGlobal(topology::Field<topology::SubMesh>* field,
+ const topology::Field<topology::SubMesh>& faultOrientation)
+{ // faultToGlobal
+ assert(field);
+
+ // Fiber dimension of vector field matches spatial dimension.
+ const spatialdata::geocoords::CoordSys* cs = field->mesh().coordsys();
+ assert(cs);
+ const int spaceDim = cs->spaceDim();
+ scalar_array fieldVertexGlobal(spaceDim);
+
+ // Get sections.
+ const ALE::Obj<RealSection>& fieldSection = field->section();
+ assert(!fieldSection.isNull());
+
+ const ALE::Obj<RealSection>& orientationSection = faultOrientation.section();
+ assert(!orientationSection.isNull());
+
+ const ALE::Obj<SieveSubMesh>& sieveMesh = field->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();
+
+ for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin; v_iter != verticesEnd; ++v_iter) {
+ assert(spaceDim == fieldSection->getFiberDimension(*v_iter));
+ const PylithScalar* fieldVertexFault = fieldSection->restrictPoint(*v_iter);
+ assert(fieldVertexFault);
+
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(*v_iter));
+ const PylithScalar* orientationVertex = orientationSection->restrictPoint(*v_iter);
+ assert(orientationVertex);
+
+ // Rotate from fault to global coordinate system (transpose orientation)
+ fieldVertexGlobal = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < spaceDim; ++jDim)
+ fieldVertexGlobal[iDim] +=
+ orientationVertex[jDim*spaceDim+iDim] * fieldVertexFault[jDim];
+
+ assert(fieldVertexGlobal.size() == fieldSection->getFiberDimension(*v_iter));
+ fieldSection->updatePoint(*v_iter, &fieldVertexGlobal[0]);
+ } // for
+
+ PetscLogFlops(vertices->size() * (2*spaceDim*spaceDim) );
+
+#if 0 // DEBUGGING
+ field->view("FIELD (GLOBAL)");
+#endif
+} // faultToGlobal
+
+// ----------------------------------------------------------------------
+// Transform field from global coordinate system to local (fault)
+// coordinate system.
+void
+pylith::faults::FaultCohesiveLagrange::globalToFault(topology::Field<topology::SubMesh>* field,
+ const topology::Field<topology::SubMesh>& faultOrientation)
+{ // globalToFault
+ assert(field);
+
+ // Fiber dimension of vector field matches spatial dimension.
+ const spatialdata::geocoords::CoordSys* cs = field->mesh().coordsys();
+ assert(cs);
+ const int spaceDim = cs->spaceDim();
+ scalar_array fieldVertexFault(spaceDim);
+
+ // Get sections.
+ const ALE::Obj<RealSection>& fieldSection = field->section();
+ assert(!fieldSection.isNull());
+
+ const ALE::Obj<RealSection>& orientationSection = faultOrientation.section();
+ assert(!orientationSection.isNull());
+
+ const ALE::Obj<SieveSubMesh>& sieveMesh = field->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();
+
+ for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin; v_iter != verticesEnd; ++v_iter) {
+ assert(spaceDim == fieldSection->getFiberDimension(*v_iter));
+ const PylithScalar* fieldVertexGlobal = fieldSection->restrictPoint(*v_iter);
+ assert(fieldVertexGlobal);
+
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(*v_iter));
+ const PylithScalar* orientationVertex = orientationSection->restrictPoint(*v_iter);
+ assert(orientationVertex);
+
+ // Rotate from global coordinate system to fault (orientation)
+ fieldVertexFault = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int jDim=0; jDim < spaceDim; ++jDim)
+ fieldVertexFault[iDim] +=
+ orientationVertex[iDim*spaceDim+jDim] * fieldVertexGlobal[jDim];
+
+ assert(fieldVertexFault.size() == fieldSection->getFiberDimension(*v_iter));
+ fieldSection->updatePoint(*v_iter, &fieldVertexFault[0]);
+ } // for
+
+ PetscLogFlops(vertices->size() * (2*spaceDim*spaceDim) );
+
+#if 0 // DEBUGGING
+ field->view("FIELD (FAULT)");
+#endif
+} // faultToGlobal
+
+// ----------------------------------------------------------------------
// Calculate orientation at fault vertices.
void
pylith::faults::FaultCohesiveLagrange::_calcOrientation(const PylithScalar upDir[3])
@@ -1178,10 +1290,8 @@
// Get vertices in fault mesh.
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveSubMesh::label_sequence::iterator verticesBegin =
- vertices->begin();
+ const ALE::Obj<SieveSubMesh::label_sequence>& vertices = faultSieveMesh->depthStratum(0);
+ const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
// Containers for orientation information.
@@ -1599,112 +1709,6 @@
} // _calcTractionsChange
// ----------------------------------------------------------------------
-// Transform field from local (fault) coordinate system to
-// global coordinate system.
-void
-pylith::faults::FaultCohesiveLagrange::_faultToGlobal(topology::Field<topology::SubMesh>* field)
-{ // _faultToGlobal
- assert(field);
- assert(0 != _faultMesh);
- assert(0 != _fields);
-
- // Fiber dimension of vector field matches spatial dimension.
- const int spaceDim = _quadrature->spaceDim();
- scalar_array fieldVertexGlobal(spaceDim);
-
- // Get sections.
- const ALE::Obj<RealSection>& fieldSection = field->section();
- assert(!fieldSection.isNull());
-
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
- assert(!orientationSection.isNull());
-
- const int numVertices = _cohesiveVertices.size();
- for (int iVertex=0; iVertex < numVertices; ++iVertex) {
- const int v_fault = _cohesiveVertices[iVertex].fault;
-
- assert(spaceDim == fieldSection->getFiberDimension(v_fault));
- const PylithScalar* fieldVertexFault = fieldSection->restrictPoint(v_fault);
- assert(fieldVertexFault);
-
- assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
- const PylithScalar* orientationVertex = orientationSection->restrictPoint(v_fault);
- assert(orientationVertex);
-
- // Rotate from fault to global coordinate system (transpose orientation)
- fieldVertexGlobal = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim)
- for (int jDim=0; jDim < spaceDim; ++jDim)
- fieldVertexGlobal[iDim] +=
- orientationVertex[jDim*spaceDim+iDim] * fieldVertexFault[jDim];
-
- assert(fieldVertexGlobal.size() ==
- fieldSection->getFiberDimension(v_fault));
- fieldSection->updatePoint(v_fault, &fieldVertexGlobal[0]);
- } // for
-
- PetscLogFlops(numVertices * (2*spaceDim*spaceDim) );
-
-#if 0 // DEBUGGING
- field->view("FIELD (GLOBAL)");
-#endif
-} // _faultToGlobal
-
-// ----------------------------------------------------------------------
-// Transform field from global coordinate system to local (fault)
-// coordinate system.
-void
-pylith::faults::FaultCohesiveLagrange::_globalToFault(topology::Field<topology::SubMesh>* field)
-{ // _globalToFault
- assert(field);
- assert(0 != _faultMesh);
- assert(0 != _fields);
-
- // Fiber dimension of vector field matches spatial dimension.
- const int spaceDim = _quadrature->spaceDim();
- scalar_array fieldVertexFault(spaceDim);
-
- // Get sections.
- const ALE::Obj<RealSection>& fieldSection = field->section();
- assert(!fieldSection.isNull());
-
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
- assert(!orientationSection.isNull());
-
- const int numVertices = _cohesiveVertices.size();
- for (int iVertex=0; iVertex < numVertices; ++iVertex) {
- const int v_fault = _cohesiveVertices[iVertex].fault;
-
- assert(spaceDim == fieldSection->getFiberDimension(v_fault));
- const PylithScalar* fieldVertexGlobal = fieldSection->restrictPoint(v_fault);
- assert(fieldVertexGlobal);
-
- assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
- const PylithScalar* orientationVertex = orientationSection->restrictPoint(v_fault);
- assert(orientationVertex);
-
- // Rotate from global coordinate system to fault (orientation)
- fieldVertexFault = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim)
- for (int jDim=0; jDim < spaceDim; ++jDim)
- fieldVertexFault[iDim] +=
- orientationVertex[iDim*spaceDim+jDim] * fieldVertexGlobal[jDim];
-
- assert(fieldVertexFault.size() ==
- fieldSection->getFiberDimension(v_fault));
- fieldSection->updatePoint(v_fault, &fieldVertexFault[0]);
- } // for
-
- PetscLogFlops(numVertices * (2*spaceDim*spaceDim) );
-
-#if 0 // DEBUGGING
- field->view("FIELD (FAULT)");
-#endif
-} // _faultToGlobal
-
-// ----------------------------------------------------------------------
// Allocate buffer for vector field.
void
pylith::faults::FaultCohesiveLagrange::_allocateBufferVectorField(void)
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.hh 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveLagrange.hh 2012-05-05 00:13:07 UTC (rev 20044)
@@ -206,6 +206,26 @@
cellField(const char* name,
const topology::SolutionFields* fields =0);
+ /** Transform field from local (fault) coordinate system to
+ * global coordinate system.
+ *
+ * @param field Field to transform.
+ * @param faultOrientation Orientation of vertices on fault.
+ */
+ static
+ void faultToGlobal(topology::Field<topology::SubMesh>* field,
+ const topology::Field<topology::SubMesh>& faultOrientation);
+
+ /** Transform field from global coordinate system to local (fault)
+ * coordinate system.
+ *
+ * @param field Field to transform.
+ * @param faultOrientation Orientation of vertices on fault.
+ */
+ static
+ void globalToFault(topology::Field<topology::SubMesh>* field,
+ const topology::Field<topology::SubMesh>& faultOrientation);
+
// PROTECTED STRUCTS //////////////////////////////////////////////////
protected :
@@ -236,20 +256,6 @@
void _calcTractionsChange(topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& solution);
- /** Transform field from local (fault) coordinate system to
- * global coordinate system.
- *
- * @param field Field to transform.
- */
- void _faultToGlobal(topology::Field<topology::SubMesh>* field);
-
- /** Transform field from global coordinate system to local (fault)
- * coordinate system.
- *
- * @param field Field to transform.
- */
- void _globalToFault(topology::Field<topology::SubMesh>* field);
-
/// Allocate buffer for vector field.
void _allocateBufferVectorField(void);
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Makefile.am 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Makefile.am 2012-05-05 00:13:07 UTC (rev 20044)
@@ -42,9 +42,7 @@
StepSlipFn.icc \
TimeHistorySlipFn.hh \
TimeHistorySlipFn.icc \
- Nucleator.hh \
TractPerturbation.hh \
- StaticPerturbation.hh \
faultsfwd.hh
noinst_HEADERS = \
Deleted: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Nucleator.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Nucleator.cc 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Nucleator.cc 2012-05-05 00:13:07 UTC (rev 20044)
@@ -1,110 +0,0 @@
-// -*- 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 "Nucleator.hh" // implementation of object methods
-
-#include "TractPerturbation.hh" // USES SlipTimeFn
-
-#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-
-#include <cassert> // USES assert()
-
-// ----------------------------------------------------------------------
-// Default constructor.
-pylith::faults::Nucleator::Nucleator(void) :
- _originTime(0.0),
- _tractfn(0)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor.
-pylith::faults::Nucleator::~Nucleator(void)
-{ // destructor
- deallocate();
-} // destructor
-
-// ----------------------------------------------------------------------
-// Deallocate PETSc and local data structures.
-void
-pylith::faults::Nucleator::deallocate(void)
-{ // deallocate
- _tractfn = 0; // :TODO: Use shared pointer.
-} // deallocate
-
-// ----------------------------------------------------------------------
-// Set origin time for earthquake source.
-void
-pylith::faults::Nucleator::originTime(const PylithScalar value)
-{ // originTime
- _originTime = value;
-} // originTime
-
-// ----------------------------------------------------------------------
-// Get origin time for earthquake source.
-PylithScalar
-pylith::faults::Nucleator::originTime(void) const
-{ // originTime
- return _originTime;
-} // originTime
-
-// ----------------------------------------------------------------------
-// Set slip time function.
-void
-pylith::faults::Nucleator::perturbationFn(TractPerturbation* tractfn)
-{ // perturbationFn
- _tractfn = tractfn; // :TODO: Use shared pointer.
-} // perturbationFn
-
-// ----------------------------------------------------------------------
-// Initialize slip time function.
-void
-pylith::faults::Nucleator::initialize(
- const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer)
-{ // initialize
- // :TODO: Normalize origin time in Python?
- normalizer.nondimensionalize(&_originTime, 1, normalizer.timeScale());
- assert(_tractfn);
- _tractfn->initialize(faultMesh, normalizer);
-} // initialize
-
-// ----------------------------------------------------------------------
-// Get slip on fault surface at time t.
-void
-pylith::faults::Nucleator::traction(
- topology::Field<topology::SubMesh>* const tractionField,
- const PylithScalar t)
-{ // slip
- assert(_tractfn);
- _tractfn->traction(tractionField, t-_originTime);
-} // slip
-
-// ----------------------------------------------------------------------
-// Get amplitude of spatial variation of traction.
-const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::faults::Nucleator::tractionAmp(void) const
-{ // finalSlip
- assert(_tractfn);
- return _tractfn->amplitude();
-} // tractionAmp
-
-
-// End of file
Deleted: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Nucleator.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Nucleator.hh 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/Nucleator.hh 2012-05-05 00:13:07 UTC (rev 20044)
@@ -1,115 +0,0 @@
-// -*- 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 libsrc/faults/Nucleator.hh
- *
- * @brief C++ object for managing traction perturbations on faults
- * with spontaneous rupture.
- */
-
-#if !defined(pylith_faults_nucleator_hh)
-#define pylith_faults_nucleator_hh
-
-// Include directives ---------------------------------------------------
-#include "faultsfwd.hh" // forward declarations
-
-#include "pylith/topology/topologyfwd.hh" // USES Field<SubMesh>
-
-#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
-
-// Nucleator -------------------------------------------------------------
-/** @brief Traction perturbations on faults with spontaneous rupture.
- *
- * Nucleator is responsible for providing perturbations in the
- * traction field over a fault surface at time t.
- */
-class pylith::faults::Nucleator
-{ // class Nucleator
- friend class TestNucleator; // unit testing
-
- // PUBLIC METHODS /////////////////////////////////////////////////////
-public :
-
- /// Default constructor.
- Nucleator(void);
-
- /// Destructor.
- ~Nucleator(void);
-
- /// Deallocate PETSc and local data structures.
- virtual
- void deallocate(void);
-
- /** Set origin time for perturbation.
- *
- * @param value Origin time for perturbation.
- */
- void originTime(const PylithScalar value);
-
- /** Get origin time for perturbation.
- *
- * @returns Origin time for perturbation.
- */
- PylithScalar originTime(void) const;
-
- /** Set traction perturbation.
- *
- * @param tractfn Traction perturbation function.
- */
- void perturbationFn(TractPerturbation* tractfn);
-
- /** Initialize traction perturbation function.
- *
- * @param faultMesh Finite-element mesh of fault.
- * @param normalizer Nondimensionalization of scales.
- */
- void initialize(const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer);
-
- /** Get traction perturbation on fault surface at time t.
- *
- * @param tractionField Traction field over fault mesh.
- * @param t Time t.
- */
- void traction(topology::Field<topology::SubMesh>* const tractionField,
- const PylithScalar t);
-
- /** Get amplitude of spatial variation of traction.
- *
- * @returns Spatial variation of traction amplitude.
- */
- const topology::Field<topology::SubMesh>& tractionAmp(void) const;
-
- // NOT IMPLEMENTED ////////////////////////////////////////////////////
-private :
-
- Nucleator(const Nucleator&); ///< Not implemented
- const Nucleator& operator=(const Nucleator&); ///< Not implemented
-
- // PRIVATE MEMBERS ////////////////////////////////////////////////////
-private :
-
- PylithScalar _originTime; ///< Origin time for traction perturbation.
- TractPerturbation* _tractfn; ///< Traction perturbation function.
-
-}; // class Nucleator
-
-#endif // pylith_faults_nucleator_hh
-
-
-// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/TractPerturbation.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/TractPerturbation.cc 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/TractPerturbation.cc 2012-05-05 00:13:07 UTC (rev 20044)
@@ -20,14 +20,32 @@
#include "TractPerturbation.hh" // implementation of object methods
+#include "FaultCohesiveLagrange.hh" // USES faultToGlobal()
+
#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/FieldsNew.hh" // HOLDSA FieldsNew
#include "pylith/topology/Field.hh" // USES Field
+#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
+#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include <strings.h> // USES strcasecmp()
+#include <cassert> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+
// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::RealUniformSection SubRealUniformSection;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::TractPerturbation::TractPerturbation(void) :
- _parameters(0)
+ _parameters(0),
+ _timeScale(1.0)
{ // constructor
} // constructor
@@ -43,16 +61,444 @@
void
pylith::faults::TractPerturbation::deallocate(void)
{ // deallocate
+ TimeDependent::deallocate();
+
delete _parameters; _parameters = 0;
} // deallocate
// ----------------------------------------------------------------------
+// Set label for traction perturbation.
+void
+pylith::faults::TractPerturbation::label(const char* value)
+{ // label
+ _label = value;
+} // label
+
+// ----------------------------------------------------------------------
+// Initialize traction perturbation function.
+void
+pylith::faults::TractPerturbation::initialize(const topology::SubMesh& faultMesh,
+ const topology::Field<topology::SubMesh>& faultOrientation,
+ const spatialdata::units::Nondimensional& normalizer)
+{ // initialize
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Fault");
+
+ const PylithScalar pressureScale = normalizer.pressureScale();
+ const PylithScalar timeScale = normalizer.timeScale();
+ const PylithScalar rateScale = pressureScale / timeScale;
+ _timeScale = timeScale;
+
+ const spatialdata::geocoords::CoordSys* cs = faultMesh.coordsys();
+ assert(cs);
+ const int spaceDim = cs->spaceDim();
+
+ delete _parameters;
+ _parameters = new topology::FieldsNew<topology::SubMesh>(faultMesh);
+
+ // 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);
+ if (_dbRate) {
+ _parameters->add("rate", "traction_rate", spaceDim, topology::FieldBase::VECTOR, rateScale);
+ _parameters->add("rate time", "traction_rate_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);
+ } // if
+ _parameters->allocate(topology::FieldBase::VERTICES_FIELD, 1);
+ const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
+ assert(!parametersSection.isNull());
+
+ if (_dbInitial) { // Setup initial values, if provided.
+ _dbInitial->open();
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ const char* valueNames[] = {"traction-normal"};
+ _dbInitial->queryVals(valueNames, 1);
+ break;
+ } // case 1
+ case 2 : {
+ const char* valueNames[] = {"traction-shear", "traction-normal"};
+ _dbInitial->queryVals(valueNames, 2);
+ break;
+ } // case 2
+ case 3 : {
+ const char* valueNames[] = {"traction-shear-horiz",
+ "traction-shear-vert",
+ "traction-normal"};
+ _dbInitial->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 TractPerturbation.");
+ } // switch
+ _queryDB("initial", _dbInitial, spaceDim, pressureScale, normalizer);
+ _dbInitial->close();
+ pylith::topology::Field<pylith::topology::SubMesh>& initial = _parameters->get("initial");
+ FaultCohesiveLagrange::faultToGlobal(&initial, faultOrientation);
+ } // if
+
+ if (_dbRate) { // Setup rate of change of values, if provided.
+ _dbRate->open();
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ const char* valueNames[] = {"traction-rate-normal"};
+ _dbRate->queryVals(valueNames, 1);
+ break;
+ } // case 1
+ case 2 : {
+ const char* valueNames[] = {"traction-rate-shear",
+ "traction-rate-normal"};
+ _dbRate->queryVals(valueNames, 2);
+ break;
+ } // case 2
+ case 3 : {
+ const char* valueNames[] = {"traction-rate-shear-horiz",
+ "traction-rate-shear-vert",
+ "traction-rate-normal"};
+ _dbRate->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 TractPerturbation.");
+ } // switch
+ _queryDB("rate", _dbRate, spaceDim, rateScale, normalizer);
+
+ const char* timeNames[1] = { "rate-start-time" };
+ _dbRate->queryVals(timeNames, 1);
+ _queryDB("rate time", _dbRate, 1, timeScale, normalizer);
+ _dbRate->close();
+ pylith::topology::Field<pylith::topology::SubMesh>& rate = _parameters->get("rate");
+ FaultCohesiveLagrange::faultToGlobal(&rate, faultOrientation);
+ } // if
+
+ if (_dbChange) { // Setup change of values, if provided.
+ _dbChange->open();
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ const char* valueNames[] = {"traction-normal"};
+ _dbChange->queryVals(valueNames, 1);
+ break;
+ } // case 1
+ case 2 : {
+ const char* valueNames[] = {"traction-shear", "traction-normal"};
+ _dbChange->queryVals(valueNames, 2);
+ break;
+ } // case 2
+ case 3 : {
+ const char* valueNames[] = {"traction-shear-horiz",
+ "traction-shear-vert",
+ "traction-normal"};
+ _dbChange->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 TractPerturbation.");
+ } // switch
+ _queryDB("change", _dbChange, spaceDim, pressureScale, normalizer);
+
+ const char* timeNames[1] = { "change-start-time" };
+ _dbChange->queryVals(timeNames, 1);
+ _queryDB("change time", _dbChange, 1, timeScale, normalizer);
+ _dbChange->close();
+ pylith::topology::Field<pylith::topology::SubMesh>& change = _parameters->get("change");
+ FaultCohesiveLagrange::faultToGlobal(&change, faultOrientation);
+
+ if (_dbTimeHistory)
+ _dbTimeHistory->open();
+ } // if
+
+ logger.stagePop();
+} // initialize
+
+// ----------------------------------------------------------------------
+// Get traction perturbation on fault surface at time t.
+void
+pylith::faults::TractPerturbation::traction(topology::Field<topology::SubMesh>* const tractionField,
+ const PylithScalar t)
+{ // traction
+ assert(tractionField);
+ assert(_parameters);
+
+ _calculateValue(t);
+
+ 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();
+ 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();
+ assert(!parametersSection.isNull());
+ const int parametersFiberDim = _parameters->fiberDim();
+ const int valueIndex = _parameters->sectionIndex("value");
+ const int valueFiberDim = _parameters->sectionFiberDim("value");
+ assert(valueFiberDim == tractionsVertex.size());
+ assert(valueIndex+valueFiberDim <= parametersFiberDim);
+
+ const ALE::Obj<RealSection>& tractionSection = tractionField->section();
+ assert(!tractionSection.isNull());
+
+ 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);
+
+ // Update field
+ assert(spaceDim == tractionSection->getFiberDimension(*v_iter));
+ tractionSection->updateAddPoint(*v_iter, &tractionsVertex[0]);
+ } // for
+
+} // traction
+
+// ----------------------------------------------------------------------
// Get parameter fields.
-const pylith::topology::Fields<pylith::topology::Field<pylith::topology::SubMesh> >*
+const pylith::topology::FieldsNew<pylith::topology::SubMesh>*
pylith::faults::TractPerturbation::parameterFields(void) const
{ // parameterFields
return _parameters;
} // parameterFields
+// ----------------------------------------------------------------------
+// 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)
+{ // vertexField
+ assert(_parameters);
+ assert(name);
-// End of file
+ if (0 == strcasecmp(name, "initial-value"))
+ return _parameters->get("initial");
+
+ else if (0 == strcasecmp(name, "rate-of-change"))
+ return _parameters->get("rate");
+
+ else if (0 == strcasecmp(name, "change-in-value"))
+ return _parameters->get("change");
+
+ else if (0 == strcasecmp(name, "rate-start-time"))
+ return _parameters->get("rate time");
+
+ else if (0 == strcasecmp(name, "change-start-time"))
+ return _parameters->get("change time");
+
+ else {
+ std::ostringstream msg;
+ msg << "Unknown field '" << name << "' requested for fault traction perturbation '"
+ << _label << "'.";
+ throw std::runtime_error(msg.str());
+ } // else
+
+ return _parameters->get("traction"); // Satisfy method definition
+} // vertexField
+
+// ----------------------------------------------------------------------
+// Get label of boundary condition surface.
+const char*
+pylith::faults::TractPerturbation::_getLabel(void) const
+{ // _getLabel
+ return _label.c_str();
+} // _getLabel
+
+// ----------------------------------------------------------------------
+// Query database for values.
+void
+pylith::faults::TractPerturbation::_queryDB(const char* name,
+ spatialdata::spatialdb::SpatialDB* const db,
+ const int querySize,
+ const PylithScalar scale,
+ const spatialdata::units::Nondimensional& normalizer)
+{ // _queryDB
+ assert(name);
+ assert(db);
+ assert(_parameters);
+
+ // 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 PylithScalar lengthScale = normalizer.lengthScale();
+
+ // Containers for database query results and quadrature coordinates in
+ // reference geometry.
+ scalar_array valuesVertex(querySize);
+ scalar_array coordsVertexGlobal(spaceDim);
+
+ // Get sections.
+ const ALE::Obj<RealSection>& coordinates = sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
+ assert(!parametersSection.isNull());
+ const int parametersFiberDim = _parameters->fiberDim();
+ const int valueIndex = _parameters->sectionIndex(name);
+ const int valueFiberDim = _parameters->sectionFiberDim(name);
+ assert(valueIndex+valueFiberDim <= parametersFiberDim);
+ scalar_array parametersVertex(parametersFiberDim);
+
+ // Loop over cells in boundary mesh and perform queries.
+ for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter != verticesEnd; ++v_iter) {
+ normalizer.dimensionalize(&coordsVertexGlobal[0], coordsVertexGlobal.size(), lengthScale);
+
+ valuesVertex = 0.0;
+ const int err = db->query(&valuesVertex[0], querySize, &coordsVertexGlobal[0], spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find values at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << coordsVertexGlobal[i];
+ msg << ") for traction boundary condition " << _label << "\n"
+ << "using spatial database " << db->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ normalizer.nondimensionalize(&valuesVertex[0], valuesVertex.size(), scale);
+
+ // Update section
+ assert(parametersFiberDim == parametersSection->getFiberDimension(*v_iter));
+ parametersSection->restrictPoint(*v_iter, ¶metersVertex[0], parametersVertex.size());
+ for (int i=0; i < valueFiberDim; ++i)
+ parametersVertex[valueIndex+i] = valuesVertex[i];
+
+ parametersSection->updatePoint(*v_iter, ¶metersVertex[0]);
+ } // 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/v1.7-trunk/libsrc/pylith/faults/TractPerturbation.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/TractPerturbation.hh 2012-05-04 21:29:45 UTC (rev 20043)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/TractPerturbation.hh 2012-05-05 00:13:07 UTC (rev 20044)
@@ -18,7 +18,8 @@
/** @file libsrc/faults/TractPerturbation.hh
*
- * @brief C++ abstract base class for traction perturbation function.
+ * @brief C++ implementation of a spatial and temporal perturbation in
+ * tractions.
*/
#if !defined(pylith_faults_tractperturbation_hh)
@@ -26,18 +27,18 @@
// Include directives ---------------------------------------------------
#include "faultsfwd.hh" // forward declarations
+#include "pylith/bc/TimeDependent.hh" // ISA TimeDependent
-#include "pylith/topology/topologyfwd.hh" // USES Fields<SubMesh>
+#include "pylith/topology/topologyfwd.hh" // HOLSA FieldsNew
#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
// TractPerturbation -----------------------------------------------------------
/**
- * @brief Abstract base class for traction perturbation function.
- *
- * Interface definition for traction perturbation function.
+ * @brief C++ implementation of a spatial and temporal perturbation in
+ * tractions.
*/
-class pylith::faults::TractPerturbation
+class pylith::faults::TractPerturbation : public pylith::bc::TimeDependent
{ // class TractPerturbation
friend class TestTractPerturbation; // unit testing
@@ -55,46 +56,88 @@
virtual
void deallocate(void);
+ /** Set label for traction perturbation.
+ *
+ * @param value Label.
+ */
+ void label(const char* value);
+
/** Initialize slip time function.
*
* @param faultMesh Finite-element mesh of fault.
- * @param cs Coordinate system for mesh
+ * @param faultOrientation Orientation of fault.
* @param normalizer Nondimensionalization of scales.
- * @param originTime Origin time for earthquake source.
*/
- virtual
void initialize(const topology::SubMesh& faultMesh,
- const spatialdata::units::Nondimensional& normalizer) = 0;
+ const topology::Field<topology::SubMesh>& faultOrientation,
+ const spatialdata::units::Nondimensional& normalizer);
/** Get traction perturbation on fault surface at time t.
*
- * @param tractionField Traction field over fault surface.
+ * @param tractionField Traction field over fault surface [output].
* @param t Time t.
*/
- virtual
void traction(topology::Field<topology::SubMesh>* const tractionField,
- const PylithScalar t) = 0;
+ const PylithScalar t);
- /** Get amplitude of traction perturbation.
- *
- * @returns Traction field.
- */
- virtual
- const topology::Field<topology::SubMesh>& amplitude(void) = 0;
-
/** Get parameter fields.
*
* @returns Parameter fields.
*/
- const topology::Fields<topology::Field<topology::SubMesh> >*
- parameterFields(void) const;
+ const topology::FieldsNew<topology::SubMesh>* parameterFields(void) const;
+
+ /** Get vertex field with traction perturbation information.
+ *
+ * @param name Name of field.
+ * @param fields Solution fields.
+ *
+ * @returns Traction vector field.
+ */
+ const topology::Field<topology::SubMesh>&
+ vertexField(const char* name,
+ topology::SolutionFields* const fields =0);
-// PROTECTED MEMBERS ////////////////////////////////////////////////////
+ // PROTECTED METHODS //////////////////////////////////////////////////
protected :
- /// Parameters for slip time function.
- topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
+ /** Get label of boundary condition surface.
+ *
+ * @returns Label of surface (from mesh generator).
+ */
+ const char* _getLabel(void) const;
+ /** Query database for values.
+ *
+ * @param name Name of field associated with database.
+ * @param db Spatial database with values.
+ * @param querySize Number of values at each location.
+ * @param scale Dimension scale associated with values.
+ * @param normalizer Nondimensionalization of scales.
+ */
+ void _queryDB(const char* name,
+ spatialdata::spatialdb::SpatialDB* const db,
+ const int querySize,
+ 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 :
+
+ /// Parameters for perturbations.
+ topology::FieldsNew<topology::SubMesh>* _parameters;
+
+ /// Time scale for current time.
+ PylithScalar _timeScale;
+
+ /// Label for traction perturbation.
+ std::string _label;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
More information about the CIG-COMMITS
mailing list