[cig-commits] r16393 - short/3D/PyLith/trunk/libsrc/faults
brad at geodynamics.org
brad at geodynamics.org
Tue Mar 9 11:34:24 PST 2010
Author: brad
Date: 2010-03-09 11:34:23 -0800 (Tue, 09 Mar 2010)
New Revision: 16393
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
Log:
More work on imposing initial fault tractions in integrateResidualAssembled().
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-03-09 01:34:39 UTC (rev 16392)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-03-09 19:34:23 UTC (rev 16393)
@@ -106,7 +106,7 @@
FaultCohesiveLagrange::initialize(mesh, upDir, normalDir);
// Get initial tractions using a spatial database.
- _getInitialTractions();
+ _setupInitialTractions();
// Setup fault constitutive model.
assert(0 != _friction);
@@ -155,7 +155,7 @@
// Get sections
double_array forcesInitialVertex(spaceDim);
const ALE::Obj<RealSection>& forcesInitialSection =
- _fields->get("initial force").section();
+ _fields->get("initial forces").section();
assert(!forcesInitialSection.isNull());
double_array residualVertex(spaceDim);
@@ -164,12 +164,12 @@
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
- const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
- // Get initial forces at fault vertex
+ // Get initial forces at fault vertex. Forces are in the global
+ // coordinate system so no rotation is necessary.
forcesInitialSection->restrictPoint(v_fault,
&forcesInitialVertex[0],
forcesInitialVertex.size());
@@ -186,6 +186,8 @@
} // for
PetscLogFlops(numVertices*spaceDim);
+
+ residual.view("RESIDUAL FAULT");
} // integrateResidualAssembled
// ----------------------------------------------------------------------
@@ -1015,12 +1017,11 @@
return buffer;
} else if (0 == strncasecmp("initial_traction", name, slipStrLen)) {
- // :TODO: Need to use scratch buffer to convert initial forces to
- // initial tractions.
assert(0 != _dbInitialTract);
- const topology::Field<topology::SubMesh>& initialTraction = _fields->get(
- "initial traction");
- return initialTraction;
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ _calcInitialTractions(&buffer);
+ return buffer;
} else if (0 == strcasecmp("traction", name)) {
assert(0 != fields);
@@ -1066,8 +1067,8 @@
// ----------------------------------------------------------------------
void
-pylith::faults::FaultCohesiveDyn::_getInitialTractions(void)
-{ // _getInitialTractions
+pylith::faults::FaultCohesiveDyn::_setupInitialTractions(void)
+{ // _setupInitialTractions
assert(0 != _normalizer);
assert(0 != _quadrature);
@@ -1209,17 +1210,50 @@
forcesInitial.complete(); // Assemble contributions
- //intialForces.view("INITIAL FORCES"); // DEBUGGING
-} // _getInitialTractions
+ // Rotate forces from fault coordinate system to global coordinate system
+ const int orientationSize = spaceDim * spaceDim;
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
+ double_array forcesInitialVertexFault(spaceDim);
+ double_array forcesInitialVertexGlobal(spaceDim);
+
+ const int numVertices = _cohesiveVertices.size();
+ for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_fault = _cohesiveVertices[iVertex].fault;
+
+ assert(orientationSize == orientationSection->getFiberDimension(v_fault));
+ assert(spaceDim == forcesInitialSection->getFiberDimension(v_fault));
+
+ const double* orientationVertex =
+ orientationSection->restrictPoint(v_fault);
+
+ forcesInitialSection->restrictPoint(v_fault,
+ &forcesInitialVertexFault[0],
+ forcesInitialVertexFault.size());
+
+ forcesInitialVertexGlobal = 0.0;
+ for (int iDim = 0; iDim < spaceDim; ++iDim)
+ for (int kDim = 0; kDim < spaceDim; ++kDim)
+ forcesInitialVertexGlobal[iDim] +=
+ forcesInitialVertexFault[kDim] *
+ orientationVertex[iDim*spaceDim+kDim];
+
+ assert(forcesInitialVertexGlobal.size() ==
+ forcesInitialSection->getFiberDimension(v_fault));
+ forcesInitialSection->updatePoint(v_fault, &forcesInitialVertexGlobal[0]);
+ } // for
+
+ forcesInitial.view("INITIAL FORCES"); // DEBUGGING
+} // _setupInitialTractions
+
// ----------------------------------------------------------------------
-// Compute change in tractions on fault surface using solution.
-// NOTE: We must convert vertex labels to fault vertex labels
+// Compute tractions on fault surface using solution.
void
pylith::faults::FaultCohesiveDyn::_calcTractions(
topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& dispT)
-{ // _calcTractionsChange
+{ // _calcTractions
assert(0 != tractions);
assert(0 != _faultMesh);
assert(0 != _fields);
@@ -1242,8 +1276,7 @@
//logger.stagePush("Fault");
const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- tractions->newSection(slip, fiberDim);
- tractions->allocate();
+ tractions->cloneSection(slip);
//logger.stagePop();
} // if
@@ -1282,6 +1315,88 @@
} // _calcTractions
// ----------------------------------------------------------------------
+// Compute initial tractions on fault surface.
+void
+pylith::faults::FaultCohesiveDyn::_calcInitialTractions(
+ topology::Field<topology::SubMesh>* tractions)
+{ // _calcInitialTractions
+ assert(0 != tractions);
+ assert(0 != _faultMesh);
+ assert(0 != _fields);
+ assert(0 != _normalizer);
+
+ // Fiber dimension of tractions matches spatial dimension.
+ const int spaceDim = _quadrature->spaceDim();
+ double_array tractionsVertexGlobal(spaceDim);
+ double_array tractionsVertexFault(spaceDim);
+
+ // Get sections.
+ const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
+ assert(!areaSection.isNull());
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
+ assert(!orientationSection.isNull());
+ const ALE::Obj<RealSection>& forcesInitialSection =
+ _fields->get("initial forces").section();
+ assert(!forcesInitialSection.isNull());
+
+ // Allocate buffer for tractions field (if necessary).
+ const ALE::Obj<RealSection>& tractionsSection = tractions->section();
+ if (tractionsSection.isNull()) {
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ //logger.stagePush("Fault");
+
+ const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
+ tractions->cloneSection(slip);
+
+ //logger.stagePop();
+ } // if
+ const double pressureScale = _normalizer->pressureScale();
+ tractions->label("initial_traction");
+ tractions->scale(pressureScale);
+ tractions->zero();
+
+ const int numVertices = _cohesiveVertices.size();
+ for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_fault = _cohesiveVertices[iVertex].fault;
+
+ assert(spaceDim == forcesInitialSection->getFiberDimension(v_fault));
+ assert(spaceDim == tractionsSection->getFiberDimension(v_fault));
+ assert(1 == areaSection->getFiberDimension(v_fault));
+
+ const double* forcesInitialVertex =
+ forcesInitialSection->restrictPoint(v_fault);
+ assert(0 != forcesInitialVertex);
+ const double* areaVertex = areaSection->restrictPoint(v_fault);
+ assert(0 != areaVertex);
+ const double* orientationVertex =
+ orientationSection->restrictPoint(v_fault);
+ assert(0 != orientationVertex);
+
+ for (int i = 0; i < spaceDim; ++i)
+ tractionsVertexGlobal[i] = forcesInitialVertex[i] / areaVertex[0];
+
+ // Rotate from global coordinate system to local coordinate system
+ tractionsVertexFault = 0.0;
+ for (int iDim = 0; iDim < spaceDim; ++iDim)
+ for (int kDim = 0; kDim < spaceDim; ++kDim)
+ tractionsVertexFault[iDim] +=
+ tractionsVertexGlobal[kDim] * orientationVertex[kDim*spaceDim+iDim];
+
+ assert(tractionsVertexFault.size() ==
+ tractionsSection->getFiberDimension(v_fault));
+ tractionsSection->updateAddPoint(v_fault, &tractionsVertexFault[0]);
+ } // for
+
+ PetscLogFlops(numVertices * (1 + spaceDim) );
+
+#if 0 // DEBUGGING
+ tractions->view("TRACTIONS");
+#endif
+
+} // _calcInitialTractions
+
+// ----------------------------------------------------------------------
// Update slip rate associated with Lagrange vertex k corresponding
// to diffential velocity between conventional vertices i and j.
void
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2010-03-09 01:34:39 UTC (rev 16392)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2010-03-09 19:34:23 UTC (rev 16393)
@@ -143,9 +143,9 @@
/** Get initial tractions using a spatial database.
*/
- void _getInitialTractions(void);
+ void _setupInitialTractions(void);
- /** Compute change in tractions on fault surface using solution.
+ /** Compute tractions on fault surface using solution.
*
* @param tractions Field for tractions.
* @param solution Solution over domain
@@ -153,6 +153,12 @@
void _calcTractions(topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& solution);
+ /** Compute initial tractions on fault surface.
+ *
+ * @param tractions Field for tractions.
+ */
+ void _calcInitialTractions(topology::Field<topology::SubMesh>* tractions);
+
/** Update slip rate associated with Lagrange vertex k corresponding
* to diffential velocity between conventional vertices i and j.
*
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-03-09 01:34:39 UTC (rev 16392)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-03-09 19:34:23 UTC (rev 16393)
@@ -1253,7 +1253,6 @@
// ----------------------------------------------------------------------
// Compute change in tractions on fault surface using solution.
-// NOTE: We must convert vertex labels to fault vertex labels
void
pylith::faults::FaultCohesiveLagrange::_calcTractionsChange(
topology::Field<topology::SubMesh>* tractions,
@@ -1268,8 +1267,8 @@
tractions->scale(_normalizer->pressureScale());
// Fiber dimension of tractions matches spatial dimension.
- const int fiberDim = _quadrature->spaceDim();
- double_array tractionsVertex(fiberDim);
+ const int spaceDim = _quadrature->spaceDim();
+ double_array tractionsVertex(spaceDim);
// Get sections.
const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
@@ -1284,8 +1283,7 @@
//logger.stagePush("Fault");
const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- tractions->newSection(slip, fiberDim);
- tractions->allocate();
+ tractions->cloneSection(slip);
//logger.stagePop();
} // if
@@ -1299,8 +1297,8 @@
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
- assert(fiberDim == dispTSection->getFiberDimension(v_lagrange));
- assert(fiberDim == tractionsSection->getFiberDimension(v_fault));
+ assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
+ assert(spaceDim == tractionsSection->getFiberDimension(v_fault));
assert(1 == areaSection->getFiberDimension(v_fault));
const double* dispTVertex = dispTSection->restrictPoint(v_lagrange);
@@ -1308,14 +1306,14 @@
const double* areaVertex = areaSection->restrictPoint(v_fault);
assert(0 != areaVertex);
- for (int i = 0; i < fiberDim; ++i)
+ for (int i = 0; i < spaceDim; ++i)
tractionsVertex[i] = dispTVertex[i] / areaVertex[0];
assert(tractionsVertex.size() == tractionsSection->getFiberDimension(v_fault));
tractionsSection->updatePoint(v_fault, &tractionsVertex[0]);
} // for
- PetscLogFlops(numVertices * (1 + fiberDim) );
+ PetscLogFlops(numVertices * (1 + spaceDim) );
#if 0 // DEBUGGING
tractions->view("TRACTIONS");
More information about the CIG-COMMITS
mailing list