[cig-commits] r15565 - in short/3D/PyLith/branches/pylith-friction: libsrc/faults modulesrc/faults pylith pylith/faults pylith/meshio
surendra at geodynamics.org
surendra at geodynamics.org
Thu Aug 20 12:43:54 PDT 2009
Author: surendra
Date: 2009-08-20 12:43:53 -0700 (Thu, 20 Aug 2009)
New Revision: 15565
Added:
short/3D/PyLith/branches/pylith-friction/pylith/faults/FaultCohesiveDyn.py
short/3D/PyLith/branches/pylith-friction/pylith/meshio/OutputFaultDyn.py
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/pylith-friction/modulesrc/faults/FaultCohesiveDyn.i
short/3D/PyLith/branches/pylith-friction/pylith/Makefile.am
Log:
Fixed some bugs in dynamic fault interface.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc 2009-08-20 16:24:13 UTC (rev 15564)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc 2009-08-20 19:43:53 UTC (rev 15565)
@@ -13,6 +13,7 @@
#include <portinfo>
#include "FaultCohesiveDyn.hh" // implementation of object methods
+#include "CohesiveTopology.hh" // USES CohesiveTopology
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include <cassert> // USES assert()
#include <sstream> // USES std::ostringstream
@@ -66,6 +67,10 @@
assert(0 != normalDir);
assert(0 != _quadrature);
+ delete _faultMesh; _faultMesh = new topology::SubMesh();
+ CohesiveTopology::createFaultParallel(_faultMesh, &_cohesiveToFault,
+ mesh, id(), useLagrangeConstraints());
+
// Reset fields.
delete _fields;
_fields =
@@ -212,19 +217,19 @@
// Use basis functions to compute displacement at quadrature point and
// then difference displacements to get slip.
// ADD STUFF HERE
- double_array dispQuadPt(numQuadPts*spaceDim*2);
+ double_array dispQuadPt(spaceDim*2);
double_array slipGlobal(spaceDim);
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- dispQuadPt[iQuad*spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
+ dispQuadPt[iSpace] += basis[iQuad*numBasis+iBasis]
*dispCell[iBasis*spaceDim+iSpace];
- dispQuadPt[(iQuad+numQuadPts)*spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
+ dispQuadPt[spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
*dispCell[(iBasis+numBasis)*spaceDim+iSpace];
}
}
for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- slipGlobal[iSpace] = dispQuadPt[(iQuad+numQuadPts)*spaceDim+iSpace]
- -dispQuadPt[iQuad*spaceDim+iSpace];
+ slipGlobal[iSpace] = dispQuadPt[spaceDim+iSpace]
+ -dispQuadPt[iSpace];
}
// Compute slip in fault orientation.
@@ -248,19 +253,19 @@
// Use basis functions to compute forces at quadrature point, then difference to get relative forces,
// and divide by area to get traction vector.
// ADD STUFF HERE
- double_array forcesQuadPt(numQuadPts*spaceDim*2);
+ double_array forcesQuadPt(2*spaceDim);
double_array tractionCurrentGlobal(spaceDim);
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- forcesQuadPt[iQuad*spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
+ forcesQuadPt[iSpace] += basis[iQuad*numBasis+iBasis]
*forcesCurrentCell[iBasis*spaceDim+iSpace];
- forcesQuadPt[(iQuad+numQuadPts)*spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
+ forcesQuadPt[spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
*forcesCurrentCell[(iBasis+numBasis)*spaceDim+iSpace];
}
}
for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- tractionCurrentGlobal[iSpace] = forcesQuadPt[(iQuad+numQuadPts)*spaceDim+iSpace]
- -forcesQuadPt[iQuad*spaceDim+iSpace];
+ tractionCurrentGlobal[iSpace] = forcesQuadPt[spaceDim+iSpace]
+ -forcesQuadPt[iSpace];
tractionCurrentGlobal[iSpace] /= wt;
}
@@ -296,7 +301,7 @@
// Out: frictionFault [spaceDim]
// BEGIN TEMPORARY
// Simple fault constitutive model with static friction.
- const double mu = 0.6;
+ const double mu = 0.7;
// ADD STUFF HERE
double_array frictionFault(spaceDim);
frictionFault = 0.0;
@@ -304,22 +309,30 @@
{ // switch
case 1 : {
if (tractionTotalFault[0] < 0) {
- frictionFault[0] = tractionCurrentFault[0];
+ frictionFault[0] = tractionCurrentFault[0]/2;
}
break;
} // case 1
case 2 : {
if (tractionTotalFault[1] < 0) {
- frictionFault[1] = tractionCurrentFault[1];
+ frictionFault[1] = tractionCurrentFault[1]/2;
frictionFault[0] = -mu * tractionTotalFault[1];
+ if (frictionFault[0] > tractionCurrentFault[0])
+ frictionFault[0] = tractionCurrentFault[0];
}
break;
} // case 2
case 3 : {
if (tractionTotalFault[2] < 0) {
- frictionFault[2] = tractionCurrentFault[2];
+ frictionFault[2] = tractionCurrentFault[2]/2;
frictionFault[1] = -mu * tractionTotalFault[2] * tractionTotalFault[1] / sqrt(pow(tractionTotalFault[1],2) +pow(tractionTotalFault[0],2));
frictionFault[0] = -mu * tractionTotalFault[2] * tractionTotalFault[0] / sqrt(pow(tractionTotalFault[1],2) +pow(tractionTotalFault[0],2));
+
+ if (frictionFault[0] > tractionCurrentFault[0])
+ frictionFault[0] = tractionCurrentFault[0];
+
+ if (frictionFault[1] > tractionCurrentFault[1])
+ frictionFault[1] = tractionCurrentFault[1];
}
break;
} // case 3
@@ -348,16 +361,33 @@
}
}
+ std::cout << "wt: " << wt
+ << ", dispQuadPt (-): (" << dispQuadPt[0] << "," << dispQuadPt[1] << ")\n"
+ << ", dispQuadPt (+): (" << dispQuadPt[spaceDim+0] << "," << dispQuadPt[spaceDim+1] << ")\n"
+ << ", slipGlobal: (" << slipGlobal[0] << "," << slipGlobal[1] << ")\n"
+ << ", slipFault: (" << slipFault[0] << "," << slipFault[1] << ")\n"
+ << ", forcesQuadPt (-): (" << forcesQuadPt[0] << "," << forcesQuadPt[1] << ")\n"
+ << ", forcesQuadPt (+): (" << forcesQuadPt[spaceDim+0] << "," << forcesQuadPt[spaceDim+1] << ")\n"
+ << ", tractionCurrentGlobal: (" << tractionCurrentGlobal[0] << "," << tractionCurrentGlobal[1] << ")\n"
+ << ", tractionCurrentFault: (" << tractionCurrentFault[0] << "," << tractionCurrentFault[1] << ")\n"
+ << ", tractionTotalFault: (" << tractionTotalFault[0] << "," << tractionTotalFault[1] << ")\n"
+ << ", frictionFault: (" << frictionFault[0] << "," << frictionFault[1] << ")\n"
+ << ", tractionCell: (" << tractionCell[iQuad*spaceDim+0] << "," << tractionCell[iQuad*spaceDim+1] << ")\n"
+ << std::endl;
+
// Compute action for dynamic fault term
for (int iBasis=0; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQuad*numBasis+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
const double valIJ = valI * basis[iQuad*numBasis+jBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
// :TODO: action for each side of the fault
residualCell[iBasis*spaceDim+iDim] +=
tractionCell[iQuad*spaceDim+iDim] * valIJ;
+ residualCell[(iBasis+numBasis)*spaceDim+iDim] +=
+ -tractionCell[iQuad*spaceDim+iDim] * valIJ;
+ } // for
} // for
} // for
} // for
@@ -410,7 +440,8 @@
<< "'.";
throw std::runtime_error(msg.str());
} // if
- const int numCorners = _quadrature->numBasis();
+
+ const int numCorners = _quadrature->refGeometry().numCorners();
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->getLabelStratum("material-id", id());
assert(!cells.isNull());
@@ -420,7 +451,7 @@
c_iter != cellsEnd;
++c_iter) {
const int cellNumCorners = sieveMesh->getNumCellCorners(*c_iter);
- if (3*numCorners != cellNumCorners) {
+ if (2*numCorners != cellNumCorners) {
std::ostringstream msg;
msg << "Number of vertices in reference cell (" << numCorners
<< ") is not compatible with number of vertices (" << cellNumCorners
@@ -667,6 +698,9 @@
} // for
_dbInitial->close();
+
+ // debugging
+ traction.view("INITIAL TRACTIONS");
} // if
} // _getInitialTractions
Modified: short/3D/PyLith/branches/pylith-friction/modulesrc/faults/FaultCohesiveDyn.i
===================================================================
--- short/3D/PyLith/branches/pylith-friction/modulesrc/faults/FaultCohesiveDyn.i 2009-08-20 16:24:13 UTC (rev 15564)
+++ short/3D/PyLith/branches/pylith-friction/modulesrc/faults/FaultCohesiveDyn.i 2009-08-20 19:43:53 UTC (rev 15565)
@@ -34,6 +34,12 @@
/// Deallocate PETSc and local data structures.
virtual
void deallocate(void);
+
+ /** Sets the spatial database for the inital tractions
+ * @param dbs spatial database for initial tractions
+ */
+ void dbInitial(spatialdata::spatialdb::SpatialDB* dbs);
+
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
Modified: short/3D/PyLith/branches/pylith-friction/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/pylith/Makefile.am 2009-08-20 16:24:13 UTC (rev 15564)
+++ short/3D/PyLith/branches/pylith-friction/pylith/Makefile.am 2009-08-20 19:43:53 UTC (rev 15565)
@@ -32,6 +32,7 @@
faults/Fault.py \
faults/FaultCohesive.py \
faults/FaultCohesiveKin.py \
+ faults/FaultCohesiveDyn.py \
faults/LiuCosSlipFn.py \
faults/SingleRupture.py \
faults/SlipTimeFn.py \
@@ -78,6 +79,7 @@
meshio/OutputManagerSubMesh.py \
meshio/OutputSoln.py \
meshio/OutputFaultKin.py \
+ meshio/OutputFaultDyn.py \
meshio/OutputMatElastic.py \
meshio/OutputNeumann.py \
meshio/OutputSolnSubset.py \
Added: short/3D/PyLith/branches/pylith-friction/pylith/faults/FaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/pylith/faults/FaultCohesiveDyn.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/pylith/faults/FaultCohesiveDyn.py 2009-08-20 19:43:53 UTC (rev 15565)
@@ -0,0 +1,220 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/faults/FaultCohesiveDyn.py
+##
+
+## @brief Python object for a fault surface with dynamic
+## (friction) fault implemented with cohesive elements.
+##
+## Factory: fault
+
+from FaultCohesive import FaultCohesive
+from pylith.feassemble.Integrator import Integrator
+from faults import FaultCohesiveDyn as ModuleFaultCohesiveDyn
+
+from pylith.utils.NullComponent import NullComponent
+
+# FaultCohesiveDyn class
+class FaultCohesiveDyn(FaultCohesive, Integrator, ModuleFaultCohesiveDyn):
+ """
+ Python object for a fault surface with dynamic (friction) fault
+ implemented with cohesive elements.
+
+ Inventory
+
+ @class Inventory
+ Python object for managing FaultCohesiveDyn facilities and properties.
+
+ \b Properties
+ @li None
+
+ \b Facilities
+ @li \b db_initial_tractions Spatial database for initial tractions.
+ @li \b output Output manager associated with fault data.
+
+ Factory: fault
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ import pyre.inventory
+
+ db = pyre.inventory.facility("db_initial_tractions", family="spatial_database",
+ factory=NullComponent)
+ db.meta['tip'] = "Spatial database for initial tractions."
+
+ from pylith.meshio.OutputFaultDyn import OutputFaultDyn
+ output = pyre.inventory.facility("output", family="output_manager",
+ factory=OutputFaultDyn)
+ output.meta['tip'] = "Output manager associated with fault data."
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="faultcohesivedyn"):
+ """
+ Initialize configuration.
+ """
+ FaultCohesive.__init__(self, name)
+ Integrator.__init__(self)
+ self._loggingPrefix = "CoDy "
+
+ self.availableFields = \
+ {'vertex': \
+ {'info': [],
+ 'data': []},
+ 'cell': \
+ {'info': ["normal_dir"],
+ 'data': ["slip",
+ "traction"]},
+}
+ return
+
+
+ def preinitialize(self, mesh):
+ """
+ Do pre-initialization setup.
+ """
+ self._info.log("Pre-initializing fault '%s'." % self.label())
+ FaultCohesive.preinitialize(self, mesh)
+ Integrator.preinitialize(self, mesh)
+
+ ModuleFaultCohesiveDyn.quadrature(self, self.faultQuadrature)
+
+ if mesh.dimension() == 2:
+ self.availableFields['cell']['info'] += ["strike_dir"]
+ elif mesh.dimension() == 3:
+ self.availableFields['cell']['info'] += ["strike_dir",
+ "dip_dir"]
+
+ if not isinstance(self.inventory.db, NullComponent):
+ self.availableFields['cell']['info'] += ["initial_traction"]
+ return
+
+
+ def verifyConfiguration(self):
+ """
+ Verify compatibility of configuration.
+ """
+ logEvent = "%sverify" % self._loggingPrefix
+ self._eventLogger.eventBegin(logEvent)
+
+ FaultCohesive.verifyConfiguration(self)
+ Integrator.verifyConfiguration(self)
+ ModuleFaultCohesiveDyn.verifyConfiguration(self, self.mesh)
+
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+ def initialize(self, totalTime, numTimeSteps, normalizer):
+ """
+ Initialize cohesive elements.
+ """
+ logEvent = "%sinit" % self._loggingPrefix
+ self._eventLogger.eventBegin(logEvent)
+ self._info.log("Initializing fault '%s'." % self.label())
+
+ Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
+
+ FaultCohesive.initialize(self, totalTime, numTimeSteps, normalizer)
+
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+ def poststep(self, t, dt, totalTime, fields):
+ """
+ Hook for doing stuff after advancing time step.
+ """
+ logEvent = "%spoststep" % self._loggingPrefix
+ self._eventLogger.eventBegin(logEvent)
+
+ Integrator.poststep(self, t, dt, totalTime, fields)
+ FaultCohesive.poststep(self, t, dt, totalTime, fields)
+
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+ def getVertexField(self, name, fields=None):
+ """
+ Get vertex field.
+ """
+ if None == fields:
+ field = ModuleFaultCohesiveDyn.vertexField(self, name)
+ else:
+ field = ModuleFaultCohesiveDyn.vertexField(self, name, fields)
+ return field
+
+
+ def getCellField(self, name, fields=None):
+ """
+ Get cell field.
+ """
+ if None == fields:
+ field = ModuleFaultCohesiveDyn.cellField(self, name)
+ else:
+ field = ModuleFaultCohesiveDyn.cellField(self, name, fields)
+ return field
+
+
+ def finalize(self):
+ """
+ Cleanup.
+ """
+ FaultCohesive.finalize(self)
+ Integrator.finalize(self)
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Setup members using inventory.
+ """
+ FaultCohesive._configure(self)
+ if not isinstance(self.inventory.db, NullComponent):
+ ModuleFaultCohesiveDyn.dbInitial(self, self.inventory.db)
+ self.output = self.inventory.output
+ return
+
+
+ def _createModuleObj(self):
+ """
+ Create handle to C++ FaultCohesiveDyn.
+ """
+ ModuleFaultCohesiveDyn.__init__(self)
+ return
+
+
+ def _modelMemoryUse(self):
+ """
+ Model memory allocation.
+ """
+ self.perfLogger.logFault("Fault", self)
+ #self.perfLogger.logFields("Fault", self.fields())
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def fault():
+ """
+ Factory associated with FaultCohesiveDyn.
+ """
+ return FaultCohesiveDyn()
+
+
+# End of file
Added: short/3D/PyLith/branches/pylith-friction/pylith/meshio/OutputFaultDyn.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/pylith/meshio/OutputFaultDyn.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/pylith/meshio/OutputFaultDyn.py 2009-08-20 19:43:53 UTC (rev 15565)
@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pyre/meshio/OutputFaultDyn.py
+##
+## @brief Python object for managing output of finite-element
+## information for faults with dynamic ruptures.
+##
+## Factory: output_manager
+
+from OutputManagerSubMesh import OutputManagerSubMesh
+
+# OutputFaultDyn class
+class OutputFaultDyn(OutputManagerSubMesh):
+ """
+ Python object for managing output of finite-element information for
+ faults with dynamic ruptures.
+
+ Inventory
+
+ @class Inventory
+ Python object for managing OutputFaultDyn facilities and properties.
+
+ \b Properties
+ @li \b cell_info_fields Names of cell info fields to output.
+ @li \b cell_data_fields Names of cell data fields to output.
+
+ \b Facilities
+ @li None
+
+ Factory: output_manager
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ import pyre.inventory
+
+ cellInfoFields = pyre.inventory.list("cell_info_fields",
+ default=["normal_dir"])
+ cellInfoFields.meta['tip'] = "Names of cell info fields to output."
+
+ cellDataFields = pyre.inventory.list("cell_data_fields",
+ default=["slip",
+ "traction"])
+ cellDataFields.meta['tip'] = "Names of cell data fields to output."
+
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="outputfaultdyn"):
+ """
+ Constructor.
+ """
+ OutputManagerSubMesh.__init__(self, name)
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based using inventory.
+ """
+ OutputManagerSubMesh._configure(self)
+ self.cellInfoFields = self.inventory.cellInfoFields
+ self.cellDataFields = self.inventory.cellDataFields
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def output_manager():
+ """
+ Factory associated with OutputManager.
+ """
+ return OutputFaultDyn()
+
+
+# End of file
More information about the CIG-COMMITS
mailing list