[cig-commits] r16279 - in short/3D/PyLith/trunk: libsrc libsrc/faults modulesrc/faults pylith pylith/faults unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Thu Feb 18 13:19:27 PST 2010
Author: brad
Date: 2010-02-18 13:19:26 -0800 (Thu, 18 Feb 2010)
New Revision: 16279
Added:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.icc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.hh
short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveTract.i
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveTract.py
Removed:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc
short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py
Modified:
short/3D/PyLith/trunk/libsrc/Makefile.am
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.icc
short/3D/PyLith/trunk/libsrc/faults/Makefile.am
short/3D/PyLith/trunk/libsrc/faults/faultsfwd.hh
short/3D/PyLith/trunk/modulesrc/faults/Makefile.am
short/3D/PyLith/trunk/modulesrc/faults/faults.i
short/3D/PyLith/trunk/pylith/Makefile.am
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
Log:
Renamed FaultCohesiveDyn to FaultCohesiveTract (will become object associated with fault traction boundary condition; no Lagrange multiplier constraints).
Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am 2010-02-18 21:19:26 UTC (rev 16279)
@@ -45,7 +45,7 @@
faults/TopologyOps.cc \
faults/CohesiveTopology.cc \
faults/FaultCohesive.cc \
- faults/FaultCohesiveDyn.cc \
+ faults/FaultCohesiveTract.cc \
faults/FaultCohesiveDynL.cc \
faults/FaultCohesiveKin.cc \
feassemble/CellGeometry.cc \
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2010-02-18 21:19:26 UTC (rev 16279)
@@ -93,8 +93,7 @@
* @returns True if implementation using Lagrange multiplier
* constraints, false otherwise.
*/
- virtual
- bool useLagrangeConstraints(void) const = 0;
+ bool useLagrangeConstraints(void) const;
// PROTECTED MEMBERS //////////////////////////////////////////////////
protected :
@@ -106,6 +105,8 @@
std::map<topology::Mesh::SieveMesh::point_type,
topology::SubMesh::SieveMesh::point_type> _cohesiveToFault;
+ bool _useLagrangeConstraints; ///< True if uses Lagrange multipliers.
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
@@ -123,6 +124,8 @@
}; // class FaultCohesive
+#include "FaultCohesive.icc" // inline methods
+
#endif // pylith_faults_faultcohesive_hh
Added: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.icc 2010-02-18 21:19:26 UTC (rev 16279)
@@ -0,0 +1,25 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_faults_faultcohesive_hh)
+#error "FaultCohesive.icc can only be included from FaultCohesive.hh"
+#endif
+
+// Cohesive cells use Lagrange multiplier constraints?
+inline
+bool
+pylith::faults::FaultCohesive::useLagrangeConstraints(void) const {
+ return _useLagrangeConstraints;
+} // useLagrangeConstraints
+
+
+// End of file
Deleted: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-02-18 21:19:26 UTC (rev 16279)
@@ -1,738 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#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
-#include <stdexcept> // USES std::runtime_error
-// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::RealSection SubRealSection;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-
-// ----------------------------------------------------------------------
-// Default constructor.
-pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(void) :
- _dbInitial(0)
-{ // constructor
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor.
-pylith::faults::FaultCohesiveDyn::~FaultCohesiveDyn(void)
-{ // destructor
- deallocate();
-} // destructor
-
-// ----------------------------------------------------------------------
-// Deallocate PETSc and local data structures.
-void
-pylith::faults::FaultCohesiveDyn::deallocate(void)
-{ // deallocate
- FaultCohesive::deallocate();
-
- _dbInitial = 0; // :TODO: Use shared pointer
-} // deallocate
-
-// ----------------------------------------------------------------------
-// Sets the spatial database for the inital tractions
-void pylith::faults::FaultCohesiveDyn::dbInitial(spatialdata::spatialdb::SpatialDB* dbs)
-{ // dbInitial
- _dbInitial = dbs;
-} // dbInitial
-
-// ----------------------------------------------------------------------
-// Initialize fault. Determine orientation and setup boundary
-void
-pylith::faults::FaultCohesiveDyn::initialize(const topology::Mesh& mesh,
- const double upDir[3],
- const double normalDir[3])
-{ // initialize
- assert(0 != upDir);
- assert(0 != normalDir);
- assert(0 != _quadrature);
-
- delete _faultMesh; _faultMesh = new topology::SubMesh();
- CohesiveTopology::createFaultParallel(_faultMesh, &_cohesiveToFault,
- mesh, id(), useLagrangeConstraints());
-
- // Reset fields.
- delete _fields;
- _fields =
- new topology::Fields<topology::Field<topology::SubMesh> >(*_faultMesh);
-
- // Initialize quadrature geometry.
- _quadrature->initializeGeometry();
-
- // Compute orientation at quadrature points in fault mesh.
- _calcOrientation(upDir, normalDir);
-
- // Get initial tractions using a spatial database.
- _getInitialTractions();
-
- // Setup fault constitutive model.
- _initConstitutiveModel();
-} // initialize
-
-// ----------------------------------------------------------------------
-// Integrate contribution of cohesive cells to residual term.
-void
-pylith::faults::FaultCohesiveDyn::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const double t,
- topology::SolutionFields* const fields)
-{ // integrateResidual
- assert(0 != _quadrature);
- assert(0 != _faultMesh);
- assert(0 != _fields);
-
- // debugging
- residual.view("RESIDUAL BEFORE FRICTION");
-
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const double_array& quadWts = _quadrature->quadWts();
- assert(quadWts.size() == numQuadPts);
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
- const int numCornersCohesive = 2*numBasis;
- const int orientationSize = spaceDim*spaceDim;
-
- // Get cohesive cells.
- const ALE::Obj<SieveMesh>& sieveMesh = residual.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
- sieveMesh->getLabelStratum("material-id", id());
- assert(!cellsCohesive.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cellsCohesive->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cellsCohesive->end();
-
- // Get fault mesh.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
-
- // Get sections.
-
- // Get residual.
- const ALE::Obj<RealSection>& residualSection = residual.section();
- assert(!residualSection.isNull());
- double_array residualCell(numCornersCohesive*spaceDim);
- topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
- &residualCell[0]);
- double_array forcesCurrentCell(numCornersCohesive*spaceDim);
- RestrictVisitor forcesCurrentVisitor(*residualSection,
- forcesCurrentCell.size(), &forcesCurrentCell[0]);
-
- // Get displacements.
- double_array dispCell(numCornersCohesive*spaceDim);
- const ALE::Obj<RealSection>& dispSection = fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection,
- dispCell.size(), &dispCell[0]);
- double_array dispIncrCell(numCornersCohesive*spaceDim);
- const ALE::Obj<RealSection>& dispIncrSection = fields->get("dispIncr(t->t+dt)").section();
- assert(!dispIncrSection.isNull());
- RestrictVisitor dispIncrVisitor(*dispIncrSection,
- dispIncrCell.size(), &dispIncrCell[0]);
- double_array dispTpdtCell(numCornersCohesive*spaceDim);
-
- // Get initial tractions (if exist).
- double_array tractionInitialFault(numQuadPts*spaceDim);
- const ALE::Obj<RealSection>& tractionInitialSection = (0 != _dbInitial) ?
- _fields->get("initial traction").section() : 0;
-
- // Get orientation.
- double_array orientationCell(numQuadPts*orientationSize);
- const ALE::Obj<RealSection>& orientationSection = _fields->get("orientation").section();
- assert(!orientationSection.isNull());
-
- // Values defined at quadrature points.
- double_array tractionsCurrentCell(numQuadPts*spaceDim);
- double_array tractionCell(numQuadPts*spaceDim); // friction
- double_array slip(spaceDim); // slip at quad pt
- double_array tractionsCurrentFault(spaceDim); // current tractions at quad pt
- double_array tractionFault(spaceDim); // friction at quad pt
-
-#if !defined(PRECOMPUTE_GEOMETRY)
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
-#endif
-
- // Loop over faces and integrate contribution from each face
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
- residualCell = 0.0;
-
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
-#else
- coordsVisitor.clear();
- faultSieveMesh->restrictClosure(c_fault, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, c_fault);
-#endif
-
- // Get cell geometry information that depends on cell
- const double_array& basis = _quadrature->basis();
- const double_array& jacobianDet = _quadrature->jacobianDet();
-
- // Get displacements at vertices of cohesive cell.
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
- dispIncrVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispIncrVisitor);
-
- // Compute current estimate of displacement at time t+dt using
- // solution increment.
- dispTpdtCell = dispCell + dispIncrCell;
-
- // Get current forces at vertices of cohesive cell (current residual).
- forcesCurrentVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, forcesCurrentVisitor);
-
- // Get initial tractions for fault cell.
- if (0 != _dbInitial) {
- assert(!tractionInitialSection.isNull());
- tractionInitialSection->restrictPoint(c_fault, &tractionInitialFault[0],
- tractionInitialFault.size());
- } // if
-
- // Get fault orientation at quadrature points.
- orientationSection->restrictPoint(c_fault, &orientationCell[0],
- orientationCell.size());
-
- tractionCell = 0.0;
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- // wt is also area associated with quadrature point
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
-
- // Compute slip at current quad pt in global coordinate system.
- // In: dispTpdtCell [2*numBasis*spaceDim] (negative side then positive side??)
- // basis [numQuadpts*numBasis]
- // Out: slipGlobal [spaceDim]
- // Use basis functions to compute displacement at quadrature point and
- // then difference displacements to get slip.
- // ADD STUFF HERE
- double_array dispQuadPt(spaceDim*2);
- double_array slipGlobal(spaceDim);
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- dispQuadPt[iSpace] += basis[iQuad*numBasis+iBasis]
- *dispTpdtCell[iBasis*spaceDim+iSpace];
- dispQuadPt[spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
- *dispTpdtCell[(iBasis+numBasis)*spaceDim+iSpace];
- }
- }
- for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- slipGlobal[iSpace] = dispQuadPt[spaceDim+iSpace]
- -dispQuadPt[iSpace];
- }
-
- // Compute slip in fault orientation.
- // In: slipGlobal [spaceDim]
- // orientationCell [numQuadPts*spaceDim*spaceDim, iQuadPt*spaceDim*spaceDim+iDim*spaceDim+jDim]
- // Out: slipFault [spaceDim]
- // Use orientation to rotate from global to fault orientation.
- // ADD STUFF HERE
- double_array slipFault(spaceDim);
- for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- for (int jSpace=0; jSpace < spaceDim; ++jSpace) {
- slipFault[iSpace] += slipGlobal[jSpace]
- *orientationCell[iQuad*spaceDim*spaceDim+iSpace*spaceDim+jSpace];
- }
- }
-
-
- // Compute traction from current deformation in global coordinate system. (tractionCurrentGlobal)
- // In: forcesCurrentCell [2*numBasis*spaceDim] (negative side then positive side)
- // Out: tractionCurrentGlobal [spaceDim]
- // 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(2*spaceDim);
- double_array tractionCurrentGlobal(spaceDim);
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- forcesQuadPt[iSpace] += basis[iQuad*numBasis+iBasis]
- *forcesCurrentCell[iBasis*spaceDim+iSpace];
- forcesQuadPt[spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
- *forcesCurrentCell[(iBasis+numBasis)*spaceDim+iSpace];
- }
- }
- for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- tractionCurrentGlobal[iSpace] = forcesQuadPt[spaceDim+iSpace]
- -forcesQuadPt[iSpace];
- tractionCurrentGlobal[iSpace] /= wt;
- }
-
-
- // Compute traction in fault orientation.
- // In: tractionCurrentGlobal [spaceDim]
- // Out: tractionCurrentFault [spaceDim]
- // Use orientation to rotate from global to fault orientation.
- // ADD STUFF HERE
- double_array tractionCurrentFault(spaceDim);
- for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- for (int jSpace=0; jSpace < spaceDim; ++jSpace) {
- tractionCurrentFault[iSpace] += tractionCurrentGlobal[jSpace]
- *orientationCell[iQuad*spaceDim*spaceDim+iSpace*spaceDim+jSpace];
- }
- }
-
- // Compute total traction (initial + current).
- // In: tractionCurrentFault [spaceDim]
- // tractionInitialFault [numQuadPts*spaceDim]
- // Out: tractionTotalFault [spaceDim]
- // ADD STUFF HERE
- double_array tractionTotalFault(spaceDim);
- for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- tractionTotalFault[iSpace] += tractionCurrentFault[iSpace]
- +tractionInitialFault[iQuad*spaceDim+iSpace];
- }
-
- // Compute traction ("friction") using fault constitutive model in fault orientation.
- // In: slipFault [spaceDim]
- // tractionCurrentFault [spaceDim]
- // tractionTotalFault [spaceDim]
- // Out: frictionFault [spaceDim]
- // BEGIN TEMPORARY
- // Simple fault constitutive model with static friction.
- const double mu = 0.7;
- // ADD STUFF HERE
- double_array frictionFault(spaceDim);
- frictionFault = 0.0;
- switch (spaceDim)
- { // switch
- case 1 : {
- if (tractionTotalFault[0] < 0) {
- frictionFault[0] = tractionCurrentFault[0];
- }
- break;
- } // case 1
- case 2 : {
- if (tractionTotalFault[1] < 0) {
- frictionFault[1] = tractionCurrentFault[1];
- 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[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
- default :
- std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
- assert(0);
- throw std::logic_error("Bad spatial dimension in Friction.");
- } // switch
-
-
- // END TEMPORARY
-
- // If normal traction is negative (compression), prevent
- // interpenetration by setting traction to exactly counteract
- // current forces acting normal to fault.
-
- // Compute traction associated with "friction" in global coordinate system. (tractionCell)
- // In: frictionFault [spaceDim]
- // Out: tractionCell [numQuadPts*spaceDim]
- // Use orientation to rotate from global to fault orientation.
- // ADD STUFF HERE
- for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
- for (int jSpace=0; jSpace < spaceDim; ++jSpace) {
- tractionCell[iQuad*spaceDim+iSpace] += frictionFault[jSpace]
- *orientationCell[iQuad*spaceDim*spaceDim+iSpace*spaceDim+jSpace];
- }
- }
- tractionCell /= 2.0;
-
- std::cout << " wt: " << wt
- << " dispTpdtCell (-): (" << dispTpdtCell[0*spaceDim+0] << "," << dispTpdtCell[0*spaceDim+1] << ")\n"
- << " (" << dispTpdtCell[1*spaceDim+0] << "," << dispTpdtCell[1*spaceDim+1] << ")\n"
- << " dispTpdtCell (+): (" << dispTpdtCell[2*spaceDim+0] << "," << dispTpdtCell[2*spaceDim+1] << ")\n"
- << " (" << dispTpdtCell[3*spaceDim+0] << "," << dispTpdtCell[3*spaceDim+1] << ")\n"
- << " 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) {
- // :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
-
- // Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
-
- PetscLogFlops(numQuadPts*(0)); // :TODO: Count number of operations
- } // for
-
- // debugging
- residual.view("RESIDUAL AFTER FRICTION");
-} // integrateResidual
-
-// ----------------------------------------------------------------------
-// Compute Jacobian matrix (A) associated with operator.
-void
-pylith::faults::FaultCohesiveDyn::integrateJacobian(
- topology::Jacobian* jacobian,
- const double t,
- topology::SolutionFields* const fields)
-{ // integrateJacobian
- _needNewJacobian = false;
-} // integrateJacobian
-
-// ----------------------------------------------------------------------
-// Verify configuration is acceptable.
-void
-pylith::faults::FaultCohesiveDyn::verifyConfiguration(
- const topology::Mesh& mesh) const
-{ // verifyConfiguration
- assert(0 != _quadrature);
-
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
-
- if (!sieveMesh->hasIntSection(label())) {
- std::ostringstream msg;
- msg << "Mesh missing group of vertices '" << label()
- << " for boundary condition.";
- throw std::runtime_error(msg.str());
- } // if
-
- // check compatibility of mesh and quadrature scheme
- const int dimension = mesh.dimension()-1;
- if (_quadrature->cellDim() != dimension) {
- std::ostringstream msg;
- msg << "Dimension of reference cell in quadrature scheme ("
- << _quadrature->cellDim()
- << ") does not match dimension of cells in mesh ("
- << dimension << ") for fault '" << label()
- << "'.";
- throw std::runtime_error(msg.str());
- } // if
-
- const int numCorners = _quadrature->refGeometry().numCorners();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- const int cellNumCorners = sieveMesh->getNumCellCorners(*c_iter);
- if (2*numCorners != cellNumCorners) {
- std::ostringstream msg;
- msg << "Number of vertices in reference cell (" << numCorners
- << ") is not compatible with number of vertices (" << cellNumCorners
- << ") in cohesive cell " << *c_iter << " for fault '"
- << label() << "'.";
- throw std::runtime_error(msg.str());
- } // if
- } // for
-} // verifyConfiguration
-
-// ----------------------------------------------------------------------
-// Get vertex field associated with integrator.
-const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::faults::FaultCohesiveDyn::vertexField(
- const char* name,
- const topology::SolutionFields* fields)
-{ // vertexField
- throw std::logic_error("FaultCohesiveDyn::vertexField() not implemented.");
-} // vertexField
-
-// ----------------------------------------------------------------------
-// Get cell field associated with integrator.
-const pylith::topology::Field<pylith::topology::SubMesh>&
-pylith::faults::FaultCohesiveDyn::cellField(
- const char* name,
- const topology::SolutionFields* fields)
-{ // cellField
- throw std::logic_error("FaultCohesiveDyn::cellField() not implemented.");
-} // cellField
-
-// ----------------------------------------------------------------------
-// Calculate orientation at fault vertices.
-void
-pylith::faults::FaultCohesiveDyn::_calcOrientation(const double upDir[3],
- const double normalDir[3])
-{ // _calcOrientation
- assert(0 != _fields);
- assert(0 != _quadrature);
-
- double_array up(upDir, 3);
-
- // Get 'fault' cells.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
-
- // Quadrature related values.
- const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
- const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
- const int numBasis = _quadrature->numBasis();
- const int numQuadPts = _quadrature->numQuadPts();
- const int spaceDim = cellGeometry.spaceDim();
- double_array quadPtRef(cellDim);
- const double_array& quadPtsRef = _quadrature->quadPtsRef();
-
- // Containers for orientation information
- const int orientationSize = spaceDim * spaceDim;
- const int fiberDim = numQuadPts * orientationSize;
- const int jacobianSize = spaceDim * cellDim;
- double_array jacobian(jacobianSize);
- double jacobianDet = 0;
- double_array orientationQuadPt(orientationSize);
- double_array orientationCell(fiberDim);
-
- // Get sections.
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
-
- // :TODO: Use spaces to create subsections like in FaultCohesiveKin.
- _fields->add("orientation", "orientation",
- topology::FieldBase::CELLS_FIELD, fiberDim);
- topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
- orientation.allocate();
- const ALE::Obj<RealSection>& orientationSection = orientation.section();
- assert(!orientationSection.isNull());
-
- // Loop over cells in fault mesh and compute orientations.
- for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- // Compute geometry information for current cell
- coordsVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
-
- // Reset orientation to zero.
- orientationCell = 0.0;
-
- // Compute orientation at each quadrature point of current cell.
- for (int iQuad=0, iRef=0, iSpace=0;
- iQuad < numQuadPts;
- ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
- // Reset orientation at quad pt to zero.
- orientationQuadPt = 0.0;
-
- // Compute Jacobian and determinant at quadrature point, then get
- // orientation.
- memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(double));
- cellGeometry.jacobian(&jacobian, &jacobianDet,
- coordinatesCell, quadPtRef);
- cellGeometry.orientation(&orientationQuadPt, jacobian, jacobianDet, up);
- assert(jacobianDet > 0.0);
- orientationQuadPt /= jacobianDet;
-
- memcpy(&orientationCell[iQuad*orientationSize],
- &orientationQuadPt[0], orientationSize*sizeof(double));
- } // for
-
- orientationSection->updatePoint(*c_iter, &orientationCell[0]);
- } // for
-
- // debugging
- orientation.view("FAULT ORIENTATION");
-} // _calcOrientation
-
-// ----------------------------------------------------------------------
-void
-pylith::faults::FaultCohesiveDyn::_getInitialTractions(void)
-{ // _getInitialTractions
- assert(0 != _normalizer);
- assert(0 != _quadrature);
-
- const double pressureScale = _normalizer->pressureScale();
- const double lengthScale = _normalizer->lengthScale();
-
- const int spaceDim = _quadrature->spaceDim();
- const int numQuadPts = _quadrature->numQuadPts();
-
- if (0 != _dbInitial) { // Setup initial values, if provided.
- // Create section to hold initial tractions.
- _fields->add("initial traction", "initial_traction");
- topology::Field<topology::SubMesh>& traction = _fields->get("initial traction");
- traction.scale(pressureScale);
- traction.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- traction.newSection(topology::FieldBase::CELLS_FIELD, numQuadPts*spaceDim);
- traction.allocate();
- const ALE::Obj<RealSection>& tractionSection = traction.section();
- assert(!tractionSection.isNull());
-
- _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-leftlateral",
- "traction-shear-updip",
- "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 Neumann.");
- } // switch
-
- // Get 'fault' cells.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
-
- const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
- const int numBasis = _quadrature->numBasis();
- const int numQuadPts = _quadrature->numQuadPts();
- const int spaceDim = _quadrature->spaceDim();
-
- // Containers for database query results and quadrature coordinates in
- // reference geometry.
- double_array tractionCell(numQuadPts*spaceDim);
- double_array quadPtsGlobal(numQuadPts*spaceDim);
-
- // Get sections.
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
-
- const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
-
- // Compute quadrature information
-
- // Loop over cells in boundary mesh and perform queries.
- for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- // Compute geometry information for current cell
- coordsVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
-
- const double_array& quadPtsNondim = _quadrature->quadPts();
- quadPtsGlobal = quadPtsNondim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
-
- tractionCell = 0.0;
- for (int iQuad=0, iSpace=0;
- iQuad < numQuadPts;
- ++iQuad, iSpace+=spaceDim) {
- const int err = _dbInitial->query(&tractionCell[iQuad*spaceDim], spaceDim,
- &quadPtsGlobal[iSpace], spaceDim, cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find initial tractions at (";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPtsGlobal[i+iSpace];
- msg << ") for dynamic fault interface " << label() << "\n"
- << "using spatial database " << _dbInitial->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
-
- } // for
- _normalizer->nondimensionalize(&tractionCell[0], tractionCell.size(),
- pressureScale);
-
- // Update section
- assert(tractionCell.size() == tractionSection->getFiberDimension(*c_iter));
- tractionSection->updatePoint(*c_iter, &tractionCell[0]);
- } // for
-
- _dbInitial->close();
-
- // debugging
- traction.view("INITIAL TRACTIONS");
- } // if
-} // _getInitialTractions
-
-// ----------------------------------------------------------------------
-void
-pylith::faults::FaultCohesiveDyn::_initConstitutiveModel(void)
-{ // _initConstitutiveModel
- // :TODO: ADD STUFF HERE
-} // _initConstitutiveModel
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2010-02-18 21:19:26 UTC (rev 16279)
@@ -1,174 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-/** @file libsrc/faults/FaultCohesiveDyn.hh
- *
- * @brief C++ implementation for a fault surface with spontaneous
- * (dynamic) slip implemented with cohesive elements.
- */
-
-#if !defined(pylith_faults_faultcohesivedyn_hh)
-#define pylith_faults_faultcohesivedyn_hh
-
-// Include directives ---------------------------------------------------
-#include "FaultCohesive.hh" // ISA FaultCohesive
-
-// FaultCohesiveDyn -----------------------------------------------------
-/**
- * @brief C++ implementation for a fault surface with spontaneous
- * (dynamic) slip implemented with cohesive elements.
- *
- * The ordering of vertices in a cohesive cell is the vertices on the
- * negative side of the fault and then the corresponding entries on
- * the positive side of the fault.
- */
-class pylith::faults::FaultCohesiveDyn : public FaultCohesive
-{ // class FaultCohesiveDyn
- friend class TestFaultCohesiveDyn; // unit testing
-
- // PUBLIC METHODS /////////////////////////////////////////////////////
-public :
-
- /// Default constructor.
- FaultCohesiveDyn(void);
-
- /// Destructor.
- virtual
- ~FaultCohesiveDyn(void);
-
- /// 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.
- *
- * @param mesh Finite-element mesh.
- * @param upDir Direction perpendicular to along-strike direction that is
- * not collinear with fault normal (usually "up" direction but could
- * be up-dip direction; only applies to fault surfaces in a 3-D domain).
- * @param normalDir General preferred direction for fault normal
- * (used to pick which of two possible normal directions for
- * interface; only applies to fault surfaces in a 3-D domain).
- */
- void initialize(const topology::Mesh& mesh,
- const double upDir[3],
- const double normalDir[3]);
-
- /** Integrate contribution of cohesive cells to residual term.
- *
- * @param residual Field containing values for residual
- * @param t Current time
- * @param fields Solution fields
- */
- void integrateResidual(const topology::Field<topology::Mesh>& residual,
- const double t,
- topology::SolutionFields* const fields);
-
- /** Integrate contributions to Jacobian matrix (A) associated with
- * operator.
- *
- * @param jacobian Sparse matrix for Jacobian of system.
- * @param t Current time
- * @param fields Solution fields
- */
- void integrateJacobian(topology::Jacobian* jacobian,
- const double t,
- topology::SolutionFields* const fields);
-
- /** Verify configuration is acceptable.
- *
- * @param mesh Finite-element mesh
- */
- void verifyConfiguration(const topology::Mesh& mesh) const;
-
- /** Get vertex field associated with integrator.
- *
- * @param name Name of vertex field.
- * @param fields Solution fields.
- *
- * @returns Vertex field.
- */
- const topology::Field<topology::SubMesh>&
- vertexField(const char* name,
- const topology::SolutionFields* fields =0);
-
- /** Get cell field associated with integrator.
- *
- * @param name Name of cell field.
- * @param fields Solution fields.
- *
- * @returns Cell field.
- */
- const topology::Field<topology::SubMesh>&
- cellField(const char* name,
- const topology::SolutionFields* fields =0);
-
- /** Cohesive cells use Lagrange multiplier constraints?
- *
- * @returns True if implementation using Lagrange multiplier
- * constraints, false otherwise.
- */
- bool useLagrangeConstraints(void) const;
-
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
-
- /** Calculate orientation at quadrature points.
- *
- * @param upDir Direction perpendicular to along-strike direction that is
- * not collinear with fault normal (usually "up" direction but could
- * be up-dip direction; only applies to fault surfaces in a 3-D domain).
- * @param normalDir General preferred direction for fault normal
- * (used to pick which of two possible normal directions for
- * interface; only applies to fault surfaces in a 3-D domain).
- */
- void _calcOrientation(const double upDir[3],
- const double normalDir[3]);
-
- /** Get initial tractions using a spatial database.
- */
- void _getInitialTractions(void);
-
- /** Setup fault constitutive model.
- */
- void _initConstitutiveModel(void);
-
- // PRIVATE MEMBERS ////////////////////////////////////////////////////
-private :
-
- /// Database for initial tractions.
- spatialdata::spatialdb::SpatialDB* _dbInitial;
-
-
- // NOT IMPLEMENTED ////////////////////////////////////////////////////
-private :
-
- /// Not implemented
- FaultCohesiveDyn(const FaultCohesiveDyn&);
-
- /// Not implemented
- const FaultCohesiveDyn& operator=(const FaultCohesiveDyn&);
-
-}; // class FaultCohesiveDyn
-
-#include "FaultCohesiveDyn.icc" // inline methods
-
-#endif // pylith_faults_faultcohesivedyn_hh
-
-
-// End of file
Deleted: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.icc 2010-02-18 21:19:26 UTC (rev 16279)
@@ -1,25 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-#if !defined(pylith_faults_faultcohesivedyn_hh)
-#error "FaultCohesiveDyn.icc can only be included from FaultCohesiveDyn.hh"
-#endif
-
-// Cohesive cells use Lagrange multiplier constraints?
-inline
-bool
-pylith::faults::FaultCohesiveDyn::useLagrangeConstraints(void) const {
- return false;
-} // useLagrangeConstraints
-
-
-// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDynL.cc 2010-02-18 21:19:26 UTC (rev 16279)
@@ -53,6 +53,7 @@
_dbInitialTract(0),
_friction(0)
{ // constructor
+ _useLagrangeConstraints = true;
_needJacobianDiag = true;
_needVelocity = true;
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-02-18 21:19:26 UTC (rev 16279)
@@ -43,18 +43,23 @@
// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void) { // constructor
+pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void)
+{ // constructor
+ _useLagrangeConstraints = true;
} // constructor
// ----------------------------------------------------------------------
// Destructor.
-pylith::faults::FaultCohesiveKin::~FaultCohesiveKin(void) { // destructor
+pylith::faults::FaultCohesiveKin::~FaultCohesiveKin(void)
+{ // destructor
deallocate();
} // destructor
// ----------------------------------------------------------------------
// Deallocate PETSc and local data structures.
-void pylith::faults::FaultCohesiveKin::deallocate(void) { // deallocate
+void
+pylith::faults::FaultCohesiveKin::deallocate(void)
+{ // deallocate
FaultCohesive::deallocate();
// :TODO: Use shared pointers for earthquake sources
@@ -62,10 +67,12 @@
// ----------------------------------------------------------------------
// Set kinematic earthquake source.
-void pylith::faults::FaultCohesiveKin::eqsrcs(const char* const * names,
- const int numNames,
- EqKinSrc** sources,
- const int numSources) { // eqsrcs
+void
+pylith::faults::FaultCohesiveKin::eqsrcs(const char* const * names,
+ const int numNames,
+ EqKinSrc** sources,
+ const int numSources)
+{ // eqsrcs
assert(numNames == numSources);
// :TODO: Use shared pointers for earthquake sources
@@ -79,9 +86,11 @@
// ----------------------------------------------------------------------
// Initialize fault. Determine orientation and setup boundary
-void pylith::faults::FaultCohesiveKin::initialize(const topology::Mesh& mesh,
- const double upDir[3],
- const double normalDir[3]) { // initialize
+void
+pylith::faults::FaultCohesiveKin::initialize(const topology::Mesh& mesh,
+ const double upDir[3],
+ const double normalDir[3])
+{ // initialize
assert(0 != upDir);
assert(0 != normalDir);
assert(0 != _quadrature);
@@ -145,7 +154,8 @@
} // initialize
// ----------------------------------------------------------------------
-void pylith::faults::FaultCohesiveKin::splitField(topology::Field<
+void
+pylith::faults::FaultCohesiveKin::splitField(topology::Field<
topology::Mesh>* field)
{ // splitField
assert(0 != field);
@@ -173,7 +183,8 @@
// ----------------------------------------------------------------------
// Integrate contribution of cohesive cells to residual term that do
// not require assembly across cells, vertices, or processors.
-void pylith::faults::FaultCohesiveKin::integrateResidualAssembled(const topology::Field<
+void
+pylith::faults::FaultCohesiveKin::integrateResidualAssembled(const topology::Field<
topology::Mesh>& residual,
const double t,
topology::SolutionFields* const fields) { // integrateResidualAssembled
@@ -334,9 +345,11 @@
// ----------------------------------------------------------------------
// Compute Jacobian matrix (A) associated with operator that do not
// require assembly across cells, vertices, or processors.
-void pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(topology::Jacobian* jacobian,
+void
+pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(topology::Jacobian* jacobian,
const double t,
- topology::SolutionFields* const fields) { // integrateJacobianAssembled
+ topology::SolutionFields* const fields)
+{ // integrateJacobianAssembled
assert(0 != jacobian);
assert(0 != fields);
assert(0 != _fields);
@@ -511,10 +524,12 @@
// ----------------------------------------------------------------------
// Compute Jacobian matrix (A) associated with operator that do not
// require assembly across cells, vertices, or processors.
-void pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(topology::Field<
+void
+pylith::faults::FaultCohesiveKin::integrateJacobianAssembled(topology::Field<
topology::Mesh>* jacobian,
const double t,
- topology::SolutionFields* const fields) { // integrateJacobianAssembled
+ topology::SolutionFields* const fields)
+{ // integrateJacobianAssembled
assert(0 != jacobian);
assert(0 != fields);
assert(0 != _fields);
@@ -569,8 +584,10 @@
// ----------------------------------------------------------------------
// Update state variables as needed.
-void pylith::faults::FaultCohesiveKin::updateStateVars(const double t,
- topology::SolutionFields* const fields) { // updateStateVars
+void
+pylith::faults::FaultCohesiveKin::updateStateVars(const double t,
+ topology::SolutionFields* const fields)
+{ // updateStateVars
assert(0 != fields);
assert(0 != _fields);
@@ -579,9 +596,11 @@
// ----------------------------------------------------------------------
// Adjust solution from solver with lumped Jacobian to match Lagrange
// multiplier constraints.
-void pylith::faults::FaultCohesiveKin::adjustSolnLumped(topology::SolutionFields* const fields,
+void
+pylith::faults::FaultCohesiveKin::adjustSolnLumped(topology::SolutionFields* const fields,
const topology::Field<
- topology::Mesh>& jacobian) { // adjustSolnLumped
+ topology::Mesh>& jacobian)
+{ // adjustSolnLumped
assert(0 != fields);
assert(0 != _quadrature);
@@ -873,7 +892,9 @@
// ----------------------------------------------------------------------
// Verify configuration is acceptable.
-void pylith::faults::FaultCohesiveKin::verifyConfiguration(const topology::Mesh& mesh) const { // verifyConfiguration
+void
+pylith::faults::FaultCohesiveKin::verifyConfiguration(const topology::Mesh& mesh) const
+{ // verifyConfiguration
assert(0 != _quadrature);
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
@@ -921,7 +942,8 @@
// Get vertex field associated with integrator.
const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::FaultCohesiveKin::vertexField(const char* name,
- const topology::SolutionFields* fields) { // vertexField
+ const topology::SolutionFields* fields)
+{ // vertexField
assert(0 != _faultMesh);
assert(0 != _quadrature);
assert(0 != _normalizer);
@@ -1025,7 +1047,8 @@
// Get cell field associated with integrator.
const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::faults::FaultCohesiveKin::cellField(const char* name,
- const topology::SolutionFields* fields) { // cellField
+ const topology::SolutionFields* fields)
+{ // cellField
// Should not reach this point if requested field was found
std::ostringstream msg;
msg << "Request for unknown cell field '" << name << "' for fault '"
@@ -1041,7 +1064,8 @@
// ----------------------------------------------------------------------
// Initialize auxiliary cohesive cell information.
-void pylith::faults::FaultCohesiveKin::_initializeCohesiveInfo(const topology::Mesh& mesh) { // _initializeCohesiveInfo
+void pylith::faults::FaultCohesiveKin::_initializeCohesiveInfo(const topology::Mesh& mesh)
+{ // _initializeCohesiveInfo
assert(0 != _quadrature);
// Get cohesive cells
@@ -1120,7 +1144,9 @@
// ----------------------------------------------------------------------
// Initialize logger.
-void pylith::faults::FaultCohesiveKin::_initializeLogger(void) { // initializeLogger
+void
+pylith::faults::FaultCohesiveKin::_initializeLogger(void)
+{ // initializeLogger
delete _logger;
_logger = new utils::EventLogger;
assert(0 != _logger);
@@ -1148,8 +1174,10 @@
// ----------------------------------------------------------------------
// Calculate orientation at fault vertices.
-void pylith::faults::FaultCohesiveKin::_calcOrientation(const double upDir[3],
- const double normalDir[3]) { // _calcOrientation
+void
+pylith::faults::FaultCohesiveKin::_calcOrientation(const double upDir[3],
+ const double normalDir[3])
+{ // _calcOrientation
assert(0 != upDir);
assert(0 != normalDir);
assert(0 != _faultMesh);
@@ -1332,7 +1360,9 @@
} // _calcOrientation
// ----------------------------------------------------------------------
-void pylith::faults::FaultCohesiveKin::_calcArea(void) { // _calcArea
+void
+pylith::faults::FaultCohesiveKin::_calcArea(void)
+{ // _calcArea
assert(0 != _faultMesh);
assert(0 != _fields);
@@ -1496,7 +1526,9 @@
// ----------------------------------------------------------------------
// Allocate buffer for vector field.
-void pylith::faults::FaultCohesiveKin::_allocateBufferVectorField(void) { // _allocateBufferVectorField
+void
+pylith::faults::FaultCohesiveKin::_allocateBufferVectorField(void)
+{ // _allocateBufferVectorField
assert(0 != _fields);
if (_fields->hasField("buffer (vector)"))
return;
@@ -1517,7 +1549,9 @@
// ----------------------------------------------------------------------
// Allocate buffer for scalar field.
-void pylith::faults::FaultCohesiveKin::_allocateBufferScalarField(void) { // _allocateBufferScalarField
+void
+pylith::faults::FaultCohesiveKin::_allocateBufferScalarField(void)
+{ // _allocateBufferScalarField
assert(0 != _fields);
if (_fields->hasField("buffer (scalar)"))
return;
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2010-02-18 21:19:26 UTC (rev 16279)
@@ -218,13 +218,6 @@
cellField(const char* name,
const topology::SolutionFields* fields =0);
- /** Cohesive cells use Lagrange multiplier constraints?
- *
- * @returns True if implementation using Lagrange multiplier
- * constraints, false otherwise.
- */
- bool useLagrangeConstraints(void) const;
-
/** Get fields associated with fault.
*
* @returns Fields associated with fault.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.icc 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.icc 2010-02-18 21:19:26 UTC (rev 16279)
@@ -14,13 +14,6 @@
#error "FaultCohesiveKin.icc can only be included from FaultCohesiveKin.hh"
#endif
-// Cohesive cells use Lagrange multiplier constraints?
-inline
-bool
-pylith::faults::FaultCohesiveKin::useLagrangeConstraints(void) const {
- return true;
-} // useLagrangeConstraints
-
// Get fields associated with fault.
inline
const pylith::topology::Fields<pylith::topology::Field<pylith::topology::SubMesh> >*
Copied: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.cc (from rev 16277, short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc)
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.cc 2010-02-18 21:19:26 UTC (rev 16279)
@@ -0,0 +1,740 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "FaultCohesiveTract.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
+#include <stdexcept> // USES std::runtime_error
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::RealSection SubRealSection;
+typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::faults::FaultCohesiveTract::FaultCohesiveTract(void) :
+ _dbInitial(0)
+{ // constructor
+ _useLagrangeConstraints = false;
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::faults::FaultCohesiveTract::~FaultCohesiveTract(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate PETSc and local data structures.
+void
+pylith::faults::FaultCohesiveTract::deallocate(void)
+{ // deallocate
+ FaultCohesive::deallocate();
+
+ _dbInitial = 0; // :TODO: Use shared pointer
+} // deallocate
+
+// ----------------------------------------------------------------------
+// Sets the spatial database for the inital tractions
+void pylith::faults::FaultCohesiveTract::dbInitial(spatialdata::spatialdb::SpatialDB* dbs)
+{ // dbInitial
+ _dbInitial = dbs;
+} // dbInitial
+
+// ----------------------------------------------------------------------
+// Initialize fault. Determine orientation and setup boundary
+void
+pylith::faults::FaultCohesiveTract::initialize(const topology::Mesh& mesh,
+ const double upDir[3],
+ const double normalDir[3])
+{ // initialize
+ assert(0 != upDir);
+ assert(0 != normalDir);
+ assert(0 != _quadrature);
+
+ delete _faultMesh; _faultMesh = new topology::SubMesh();
+ CohesiveTopology::createFaultParallel(_faultMesh, &_cohesiveToFault,
+ mesh, id(), useLagrangeConstraints());
+
+ // Reset fields.
+ delete _fields;
+ _fields =
+ new topology::Fields<topology::Field<topology::SubMesh> >(*_faultMesh);
+
+ // Initialize quadrature geometry.
+ _quadrature->initializeGeometry();
+
+ // Compute orientation at quadrature points in fault mesh.
+ _calcOrientation(upDir, normalDir);
+
+ // Get initial tractions using a spatial database.
+ _getInitialTractions();
+
+ // Setup fault constitutive model.
+ _initConstitutiveModel();
+} // initialize
+
+// ----------------------------------------------------------------------
+// Integrate contribution of cohesive cells to residual term.
+void
+pylith::faults::FaultCohesiveTract::integrateResidual(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidual
+ assert(0 != _quadrature);
+ assert(0 != _faultMesh);
+ assert(0 != _fields);
+
+ // debugging
+ residual.view("RESIDUAL BEFORE FRICTION");
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int numCornersCohesive = 2*numBasis;
+ const int orientationSize = spaceDim*spaceDim;
+
+ // Get cohesive cells.
+ const ALE::Obj<SieveMesh>& sieveMesh = residual.mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive =
+ sieveMesh->getLabelStratum("material-id", id());
+ assert(!cellsCohesive.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cellsCohesive->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cellsCohesive->end();
+
+ // Get fault mesh.
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+
+ // Get sections.
+
+ // Get residual.
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ assert(!residualSection.isNull());
+ double_array residualCell(numCornersCohesive*spaceDim);
+ topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
+ &residualCell[0]);
+ double_array forcesCurrentCell(numCornersCohesive*spaceDim);
+ RestrictVisitor forcesCurrentVisitor(*residualSection,
+ forcesCurrentCell.size(), &forcesCurrentCell[0]);
+
+ // Get displacements.
+ double_array dispCell(numCornersCohesive*spaceDim);
+ const ALE::Obj<RealSection>& dispSection = fields->get("disp(t)").section();
+ assert(!dispSection.isNull());
+ RestrictVisitor dispVisitor(*dispSection,
+ dispCell.size(), &dispCell[0]);
+ double_array dispIncrCell(numCornersCohesive*spaceDim);
+ const ALE::Obj<RealSection>& dispIncrSection = fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispIncrSection.isNull());
+ RestrictVisitor dispIncrVisitor(*dispIncrSection,
+ dispIncrCell.size(), &dispIncrCell[0]);
+ double_array dispTpdtCell(numCornersCohesive*spaceDim);
+
+ // Get initial tractions (if exist).
+ double_array tractionInitialFault(numQuadPts*spaceDim);
+ const ALE::Obj<RealSection>& tractionInitialSection = (0 != _dbInitial) ?
+ _fields->get("initial traction").section() : 0;
+
+ // Get orientation.
+ double_array orientationCell(numQuadPts*orientationSize);
+ const ALE::Obj<RealSection>& orientationSection = _fields->get("orientation").section();
+ assert(!orientationSection.isNull());
+
+ // Values defined at quadrature points.
+ double_array tractionsCurrentCell(numQuadPts*spaceDim);
+ double_array tractionCell(numQuadPts*spaceDim); // friction
+ double_array slip(spaceDim); // slip at quad pt
+ double_array tractionsCurrentFault(spaceDim); // current tractions at quad pt
+ double_array tractionFault(spaceDim); // friction at quad pt
+
+#if !defined(PRECOMPUTE_GEOMETRY)
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
+#endif
+
+ // Loop over faces and integrate contribution from each face
+ for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
+ residualCell = 0.0;
+
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ faultSieveMesh->restrictClosure(c_fault, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, c_fault);
+#endif
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Get displacements at vertices of cohesive cell.
+ dispVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispVisitor);
+ dispIncrVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispIncrVisitor);
+
+ // Compute current estimate of displacement at time t+dt using
+ // solution increment.
+ dispTpdtCell = dispCell + dispIncrCell;
+
+ // Get current forces at vertices of cohesive cell (current residual).
+ forcesCurrentVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, forcesCurrentVisitor);
+
+ // Get initial tractions for fault cell.
+ if (0 != _dbInitial) {
+ assert(!tractionInitialSection.isNull());
+ tractionInitialSection->restrictPoint(c_fault, &tractionInitialFault[0],
+ tractionInitialFault.size());
+ } // if
+
+ // Get fault orientation at quadrature points.
+ orientationSection->restrictPoint(c_fault, &orientationCell[0],
+ orientationCell.size());
+
+ tractionCell = 0.0;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ // wt is also area associated with quadrature point
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+
+ // Compute slip at current quad pt in global coordinate system.
+ // In: dispTpdtCell [2*numBasis*spaceDim] (negative side then positive side??)
+ // basis [numQuadpts*numBasis]
+ // Out: slipGlobal [spaceDim]
+ // Use basis functions to compute displacement at quadrature point and
+ // then difference displacements to get slip.
+ // ADD STUFF HERE
+ double_array dispQuadPt(spaceDim*2);
+ double_array slipGlobal(spaceDim);
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
+ dispQuadPt[iSpace] += basis[iQuad*numBasis+iBasis]
+ *dispTpdtCell[iBasis*spaceDim+iSpace];
+ dispQuadPt[spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
+ *dispTpdtCell[(iBasis+numBasis)*spaceDim+iSpace];
+ }
+ }
+ for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
+ slipGlobal[iSpace] = dispQuadPt[spaceDim+iSpace]
+ -dispQuadPt[iSpace];
+ }
+
+ // Compute slip in fault orientation.
+ // In: slipGlobal [spaceDim]
+ // orientationCell [numQuadPts*spaceDim*spaceDim, iQuadPt*spaceDim*spaceDim+iDim*spaceDim+jDim]
+ // Out: slipFault [spaceDim]
+ // Use orientation to rotate from global to fault orientation.
+ // ADD STUFF HERE
+ double_array slipFault(spaceDim);
+ for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
+ for (int jSpace=0; jSpace < spaceDim; ++jSpace) {
+ slipFault[iSpace] += slipGlobal[jSpace]
+ *orientationCell[iQuad*spaceDim*spaceDim+iSpace*spaceDim+jSpace];
+ }
+ }
+
+
+ // Compute traction from current deformation in global coordinate system. (tractionCurrentGlobal)
+ // In: forcesCurrentCell [2*numBasis*spaceDim] (negative side then positive side)
+ // Out: tractionCurrentGlobal [spaceDim]
+ // 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(2*spaceDim);
+ double_array tractionCurrentGlobal(spaceDim);
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
+ forcesQuadPt[iSpace] += basis[iQuad*numBasis+iBasis]
+ *forcesCurrentCell[iBasis*spaceDim+iSpace];
+ forcesQuadPt[spaceDim+iSpace] += basis[iQuad*numBasis+iBasis]
+ *forcesCurrentCell[(iBasis+numBasis)*spaceDim+iSpace];
+ }
+ }
+ for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
+ tractionCurrentGlobal[iSpace] = forcesQuadPt[spaceDim+iSpace]
+ -forcesQuadPt[iSpace];
+ tractionCurrentGlobal[iSpace] /= wt;
+ }
+
+
+ // Compute traction in fault orientation.
+ // In: tractionCurrentGlobal [spaceDim]
+ // Out: tractionCurrentFault [spaceDim]
+ // Use orientation to rotate from global to fault orientation.
+ // ADD STUFF HERE
+ double_array tractionCurrentFault(spaceDim);
+ for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
+ for (int jSpace=0; jSpace < spaceDim; ++jSpace) {
+ tractionCurrentFault[iSpace] += tractionCurrentGlobal[jSpace]
+ *orientationCell[iQuad*spaceDim*spaceDim+iSpace*spaceDim+jSpace];
+ }
+ }
+
+ // Compute total traction (initial + current).
+ // In: tractionCurrentFault [spaceDim]
+ // tractionInitialFault [numQuadPts*spaceDim]
+ // Out: tractionTotalFault [spaceDim]
+ // ADD STUFF HERE
+ double_array tractionTotalFault(spaceDim);
+ for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
+ tractionTotalFault[iSpace] += tractionCurrentFault[iSpace]
+ +tractionInitialFault[iQuad*spaceDim+iSpace];
+ }
+
+ // Compute traction ("friction") using fault constitutive model in fault orientation.
+ // In: slipFault [spaceDim]
+ // tractionCurrentFault [spaceDim]
+ // tractionTotalFault [spaceDim]
+ // Out: frictionFault [spaceDim]
+ // BEGIN TEMPORARY
+ // Simple fault constitutive model with static friction.
+ const double mu = 0.7;
+ // ADD STUFF HERE
+ double_array frictionFault(spaceDim);
+ frictionFault = 0.0;
+ switch (spaceDim)
+ { // switch
+ case 1 : {
+ if (tractionTotalFault[0] < 0) {
+ frictionFault[0] = tractionCurrentFault[0];
+ }
+ break;
+ } // case 1
+ case 2 : {
+ if (tractionTotalFault[1] < 0) {
+ frictionFault[1] = tractionCurrentFault[1];
+ 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[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
+ default :
+ std::cerr << "Bad spatial dimension '" << spaceDim << "'." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad spatial dimension in Friction.");
+ } // switch
+
+
+ // END TEMPORARY
+
+ // If normal traction is negative (compression), prevent
+ // interpenetration by setting traction to exactly counteract
+ // current forces acting normal to fault.
+
+ // Compute traction associated with "friction" in global coordinate system. (tractionCell)
+ // In: frictionFault [spaceDim]
+ // Out: tractionCell [numQuadPts*spaceDim]
+ // Use orientation to rotate from global to fault orientation.
+ // ADD STUFF HERE
+ for (int iSpace=0; iSpace < spaceDim; ++iSpace) {
+ for (int jSpace=0; jSpace < spaceDim; ++jSpace) {
+ tractionCell[iQuad*spaceDim+iSpace] += frictionFault[jSpace]
+ *orientationCell[iQuad*spaceDim*spaceDim+iSpace*spaceDim+jSpace];
+ }
+ }
+ tractionCell /= 2.0;
+
+ std::cout << " wt: " << wt
+ << " dispTpdtCell (-): (" << dispTpdtCell[0*spaceDim+0] << "," << dispTpdtCell[0*spaceDim+1] << ")\n"
+ << " (" << dispTpdtCell[1*spaceDim+0] << "," << dispTpdtCell[1*spaceDim+1] << ")\n"
+ << " dispTpdtCell (+): (" << dispTpdtCell[2*spaceDim+0] << "," << dispTpdtCell[2*spaceDim+1] << ")\n"
+ << " (" << dispTpdtCell[3*spaceDim+0] << "," << dispTpdtCell[3*spaceDim+1] << ")\n"
+ << " 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) {
+ // :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
+
+ // Assemble cell contribution into field
+ residualVisitor.clear();
+ sieveMesh->updateClosure(*c_iter, residualVisitor);
+
+ PetscLogFlops(numQuadPts*(0)); // :TODO: Count number of operations
+ } // for
+
+ // debugging
+ residual.view("RESIDUAL AFTER FRICTION");
+} // integrateResidual
+
+// ----------------------------------------------------------------------
+// Compute Jacobian matrix (A) associated with operator.
+void
+pylith::faults::FaultCohesiveTract::integrateJacobian(
+ topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateJacobian
+ _needNewJacobian = false;
+} // integrateJacobian
+
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+pylith::faults::FaultCohesiveTract::verifyConfiguration(
+ const topology::Mesh& mesh) const
+{ // verifyConfiguration
+ assert(0 != _quadrature);
+
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ if (!sieveMesh->hasIntSection(label())) {
+ std::ostringstream msg;
+ msg << "Mesh missing group of vertices '" << label()
+ << " for boundary condition.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ // check compatibility of mesh and quadrature scheme
+ const int dimension = mesh.dimension()-1;
+ if (_quadrature->cellDim() != dimension) {
+ std::ostringstream msg;
+ msg << "Dimension of reference cell in quadrature scheme ("
+ << _quadrature->cellDim()
+ << ") does not match dimension of cells in mesh ("
+ << dimension << ") for fault '" << label()
+ << "'.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ const int numCorners = _quadrature->refGeometry().numCorners();
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", id());
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ const int cellNumCorners = sieveMesh->getNumCellCorners(*c_iter);
+ if (2*numCorners != cellNumCorners) {
+ std::ostringstream msg;
+ msg << "Number of vertices in reference cell (" << numCorners
+ << ") is not compatible with number of vertices (" << cellNumCorners
+ << ") in cohesive cell " << *c_iter << " for fault '"
+ << label() << "'.";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // for
+} // verifyConfiguration
+
+// ----------------------------------------------------------------------
+// Get vertex field associated with integrator.
+const pylith::topology::Field<pylith::topology::SubMesh>&
+pylith::faults::FaultCohesiveTract::vertexField(
+ const char* name,
+ const topology::SolutionFields* fields)
+{ // vertexField
+ throw std::logic_error("FaultCohesiveTract::vertexField() not implemented.");
+} // vertexField
+
+// ----------------------------------------------------------------------
+// Get cell field associated with integrator.
+const pylith::topology::Field<pylith::topology::SubMesh>&
+pylith::faults::FaultCohesiveTract::cellField(
+ const char* name,
+ const topology::SolutionFields* fields)
+{ // cellField
+ throw std::logic_error("FaultCohesiveTract::cellField() not implemented.");
+} // cellField
+
+// ----------------------------------------------------------------------
+// Calculate orientation at fault vertices.
+void
+pylith::faults::FaultCohesiveTract::_calcOrientation(const double upDir[3],
+ const double normalDir[3])
+{ // _calcOrientation
+ assert(0 != _fields);
+ assert(0 != _quadrature);
+
+ double_array up(upDir, 3);
+
+ // Get 'fault' cells.
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ faultSieveMesh->heightStratum(0);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Quadrature related values.
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
+ const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
+ const int numBasis = _quadrature->numBasis();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int spaceDim = cellGeometry.spaceDim();
+ double_array quadPtRef(cellDim);
+ const double_array& quadPtsRef = _quadrature->quadPtsRef();
+
+ // Containers for orientation information
+ const int orientationSize = spaceDim * spaceDim;
+ const int fiberDim = numQuadPts * orientationSize;
+ const int jacobianSize = spaceDim * cellDim;
+ double_array jacobian(jacobianSize);
+ double jacobianDet = 0;
+ double_array orientationQuadPt(orientationSize);
+ double_array orientationCell(fiberDim);
+
+ // Get sections.
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ // :TODO: Use spaces to create subsections like in FaultCohesiveKin.
+ _fields->add("orientation", "orientation",
+ topology::FieldBase::CELLS_FIELD, fiberDim);
+ topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+ orientation.allocate();
+ const ALE::Obj<RealSection>& orientationSection = orientation.section();
+ assert(!orientationSection.isNull());
+
+ // Loop over cells in fault mesh and compute orientations.
+ for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ coordsVisitor.clear();
+ faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+
+ // Reset orientation to zero.
+ orientationCell = 0.0;
+
+ // Compute orientation at each quadrature point of current cell.
+ for (int iQuad=0, iRef=0, iSpace=0;
+ iQuad < numQuadPts;
+ ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
+ // Reset orientation at quad pt to zero.
+ orientationQuadPt = 0.0;
+
+ // Compute Jacobian and determinant at quadrature point, then get
+ // orientation.
+ memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(double));
+ cellGeometry.jacobian(&jacobian, &jacobianDet,
+ coordinatesCell, quadPtRef);
+ cellGeometry.orientation(&orientationQuadPt, jacobian, jacobianDet, up);
+ assert(jacobianDet > 0.0);
+ orientationQuadPt /= jacobianDet;
+
+ memcpy(&orientationCell[iQuad*orientationSize],
+ &orientationQuadPt[0], orientationSize*sizeof(double));
+ } // for
+
+ orientationSection->updatePoint(*c_iter, &orientationCell[0]);
+ } // for
+
+ // debugging
+ orientation.view("FAULT ORIENTATION");
+} // _calcOrientation
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::FaultCohesiveTract::_getInitialTractions(void)
+{ // _getInitialTractions
+ assert(0 != _normalizer);
+ assert(0 != _quadrature);
+
+ const double pressureScale = _normalizer->pressureScale();
+ const double lengthScale = _normalizer->lengthScale();
+
+ const int spaceDim = _quadrature->spaceDim();
+ const int numQuadPts = _quadrature->numQuadPts();
+
+ if (0 != _dbInitial) { // Setup initial values, if provided.
+ // Create section to hold initial tractions.
+ _fields->add("initial traction", "initial_traction");
+ topology::Field<topology::SubMesh>& traction = _fields->get("initial traction");
+ traction.scale(pressureScale);
+ traction.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ traction.newSection(topology::FieldBase::CELLS_FIELD, numQuadPts*spaceDim);
+ traction.allocate();
+ const ALE::Obj<RealSection>& tractionSection = traction.section();
+ assert(!tractionSection.isNull());
+
+ _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-leftlateral",
+ "traction-shear-updip",
+ "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 Neumann.");
+ } // switch
+
+ // Get 'fault' cells.
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ faultSieveMesh->heightStratum(0);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
+ const int numBasis = _quadrature->numBasis();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int spaceDim = _quadrature->spaceDim();
+
+ // Containers for database query results and quadrature coordinates in
+ // reference geometry.
+ double_array tractionCell(numQuadPts*spaceDim);
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+
+ // Get sections.
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
+
+ // Compute quadrature information
+
+ // Loop over cells in boundary mesh and perform queries.
+ for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Compute geometry information for current cell
+ coordsVisitor.clear();
+ faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+
+ const double_array& quadPtsNondim = _quadrature->quadPts();
+ quadPtsGlobal = quadPtsNondim;
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
+ lengthScale);
+
+ tractionCell = 0.0;
+ for (int iQuad=0, iSpace=0;
+ iQuad < numQuadPts;
+ ++iQuad, iSpace+=spaceDim) {
+ const int err = _dbInitial->query(&tractionCell[iQuad*spaceDim], spaceDim,
+ &quadPtsGlobal[iSpace], spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find initial tractions at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << quadPtsGlobal[i+iSpace];
+ msg << ") for dynamic fault interface " << label() << "\n"
+ << "using spatial database " << _dbInitial->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ } // for
+ _normalizer->nondimensionalize(&tractionCell[0], tractionCell.size(),
+ pressureScale);
+
+ // Update section
+ assert(tractionCell.size() == tractionSection->getFiberDimension(*c_iter));
+ tractionSection->updatePoint(*c_iter, &tractionCell[0]);
+ } // for
+
+ _dbInitial->close();
+
+ // debugging
+ traction.view("INITIAL TRACTIONS");
+ } // if
+} // _getInitialTractions
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::FaultCohesiveTract::_initConstitutiveModel(void)
+{ // _initConstitutiveModel
+ // :TODO: ADD STUFF HERE
+} // _initConstitutiveModel
+
+
+// End of file
Copied: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.hh (from rev 16277, short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh)
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.hh 2010-02-18 21:19:26 UTC (rev 16279)
@@ -0,0 +1,161 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/faults/FaultCohesiveTract.hh
+ *
+ * @brief C++ implementation for a fault surface with tractions
+ * applied to the fault surface using cohesive cells.
+ */
+
+#if !defined(pylith_faults_faultcohesivetract_hh)
+#define pylith_faults_faultcohesivetract_hh
+
+// Include directives ---------------------------------------------------
+#include "FaultCohesive.hh" // ISA FaultCohesive
+
+// FaultCohesiveTract -----------------------------------------------------
+/**
+ * @brief C++ implementation for a fault surface with tractions
+ * applied to the fault surface using cohesive cells.
+ */
+class pylith::faults::FaultCohesiveTract : public FaultCohesive
+{ // class FaultCohesiveTract
+ friend class TestFaultCohesiveTract; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ FaultCohesiveTract(void);
+
+ /// Destructor.
+ virtual
+ ~FaultCohesiveTract(void);
+
+ /// 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.
+ *
+ * @param mesh Finite-element mesh.
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction; only applies to fault surfaces in a 3-D domain).
+ * @param normalDir General preferred direction for fault normal
+ * (used to pick which of two possible normal directions for
+ * interface; only applies to fault surfaces in a 3-D domain).
+ */
+ void initialize(const topology::Mesh& mesh,
+ const double upDir[3],
+ const double normalDir[3]);
+
+ /** Integrate contribution of cohesive cells to residual term.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param jacobian Sparse matrix for Jacobian of system.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateJacobian(topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Verify configuration is acceptable.
+ *
+ * @param mesh Finite-element mesh
+ */
+ void verifyConfiguration(const topology::Mesh& mesh) const;
+
+ /** Get vertex field associated with integrator.
+ *
+ * @param name Name of vertex field.
+ * @param fields Solution fields.
+ *
+ * @returns Vertex field.
+ */
+ const topology::Field<topology::SubMesh>&
+ vertexField(const char* name,
+ const topology::SolutionFields* fields =0);
+
+ /** Get cell field associated with integrator.
+ *
+ * @param name Name of cell field.
+ * @param fields Solution fields.
+ *
+ * @returns Cell field.
+ */
+ const topology::Field<topology::SubMesh>&
+ cellField(const char* name,
+ const topology::SolutionFields* fields =0);
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Calculate orientation at quadrature points.
+ *
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction; only applies to fault surfaces in a 3-D domain).
+ * @param normalDir General preferred direction for fault normal
+ * (used to pick which of two possible normal directions for
+ * interface; only applies to fault surfaces in a 3-D domain).
+ */
+ void _calcOrientation(const double upDir[3],
+ const double normalDir[3]);
+
+ /** Get initial tractions using a spatial database.
+ */
+ void _getInitialTractions(void);
+
+ /** Setup fault constitutive model.
+ */
+ void _initConstitutiveModel(void);
+
+ // PRIVATE MEMBERS ////////////////////////////////////////////////////
+private :
+
+ /// Database for initial tractions.
+ spatialdata::spatialdb::SpatialDB* _dbInitial;
+
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ FaultCohesiveTract(const FaultCohesiveTract&);
+
+ /// Not implemented
+ const FaultCohesiveTract& operator=(const FaultCohesiveTract&);
+
+}; // class FaultCohesiveTract
+
+#endif // pylith_faults_faultcohesivetract_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/Makefile.am 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/Makefile.am 2010-02-18 21:19:26 UTC (rev 16279)
@@ -23,8 +23,8 @@
Fault.hh \
Fault.icc \
FaultCohesive.hh \
- FaultCohesiveDyn.hh \
- FaultCohesiveDyn.icc \
+ FaultCohesive.icc \
+ FaultCohesiveTract.hh \
FaultCohesiveDynL.hh \
FaultCohesiveDynL.icc \
FaultCohesiveKin.hh \
Modified: short/3D/PyLith/trunk/libsrc/faults/faultsfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/faultsfwd.hh 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/libsrc/faults/faultsfwd.hh 2010-02-18 21:19:26 UTC (rev 16279)
@@ -28,9 +28,9 @@
class Fault;
class FaultCohesive;
- class FaultCohesiveDyn;
class FaultCohesiveDynL;
class FaultCohesiveKin;
+ class FaultCohesiveTract;
class EqKinSrc;
class SlipTimeFn;
Deleted: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i 2010-02-18 21:19:26 UTC (rev 16279)
@@ -1,121 +0,0 @@
-// -*- C++ -*-
-//
-// ----------------------------------------------------------------------
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ----------------------------------------------------------------------
-//
-
-/** @file modulesrc/faults/FaultCohesiveDyn.i
- *
- * @brief Python interface to C++ FaultCohesiveDyn object.
- */
-
-namespace pylith {
- namespace faults {
-
- class FaultCohesiveDyn : public FaultCohesive
- { // class FaultCohesiveDyn
-
- // PUBLIC METHODS /////////////////////////////////////////////////
- public :
-
- /// Default constructor.
- FaultCohesiveDyn(void);
-
- /// Destructor.
- virtual
- ~FaultCohesiveDyn(void);
-
- /// 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.
- *
- * @param mesh Finite-element mesh.
- * @param upDir Direction perpendicular to along-strike direction that is
- * not collinear with fault normal (usually "up" direction but could
- * be up-dip direction; only applies to fault surfaces in a 3-D domain).
- * @param normalDir General preferred direction for fault normal
- * (used to pick which of two possible normal directions for
- * interface; only applies to fault surfaces in a 3-D domain).
- */
- void initialize(const pylith::topology::Mesh& mesh,
- const double upDir[3],
- const double normalDir[3]);
-
- /** Integrate contribution of cohesive cells to residual term.
- *
- * @param residual Field containing values for residual
- * @param t Current time
- * @param fields Solution fields
- */
- void integrateResidual(const pylith::topology::Field<pylith::topology::Mesh>& residual,
- const double t,
- pylith::topology::SolutionFields* const fields);
-
- /** Integrate contributions to Jacobian matrix (A) associated with
- * operator.
- *
- * @param jacobian Sparse matrix for Jacobian of system.
- * @param t Current time
- * @param fields Solution fields
- */
- void integrateJacobian(pylith::topology::Jacobian* jacobian,
- const double t,
- pylith::topology::SolutionFields* const fields);
-
- /** Verify configuration is acceptable.
- *
- * @param mesh Finite-element mesh
- */
- void verifyConfiguration(const pylith::topology::Mesh& mesh) const;
-
- /** Get vertex field associated with integrator.
- *
- * @param name Name of vertex field.
- * @param fields Solution fields.
- *
- * @returns Vertex field.
- */
- const pylith::topology::Field<pylith::topology::SubMesh>&
- vertexField(const char* name,
- const pylith::topology::SolutionFields* fields =0);
-
- /** Get cell field associated with integrator.
- *
- * @param name Name of cell field.
- * @param fields Solution fields.
- *
- * @returns Cell field.
- */
- const pylith::topology::Field<pylith::topology::SubMesh>&
- cellField(const char* name,
- const pylith::topology::SolutionFields* fields =0);
-
- /** Cohesive cells use Lagrange multiplier constraints?
- *
- * @returns True if implementation using Lagrange multiplier
- * constraints, false otherwise.
- */
- bool useLagrangeConstraints(void) const;
-
- }; // class FaultCohesiveDyn
-
- } // faults
-} // pylith
-
-
-// End of file
Copied: short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveTract.i (from rev 16277, short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveDyn.i)
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveTract.i (rev 0)
+++ short/3D/PyLith/trunk/modulesrc/faults/FaultCohesiveTract.i 2010-02-18 21:19:26 UTC (rev 16279)
@@ -0,0 +1,114 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file modulesrc/faults/FaultCohesiveTract.i
+ *
+ * @brief Python interface to C++ FaultCohesiveTract object.
+ */
+
+namespace pylith {
+ namespace faults {
+
+ class FaultCohesiveTract : public FaultCohesive
+ { // class FaultCohesiveTract
+
+ // PUBLIC METHODS /////////////////////////////////////////////////
+ public :
+
+ /// Default constructor.
+ FaultCohesiveTract(void);
+
+ /// Destructor.
+ virtual
+ ~FaultCohesiveTract(void);
+
+ /// 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.
+ *
+ * @param mesh Finite-element mesh.
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction; only applies to fault surfaces in a 3-D domain).
+ * @param normalDir General preferred direction for fault normal
+ * (used to pick which of two possible normal directions for
+ * interface; only applies to fault surfaces in a 3-D domain).
+ */
+ void initialize(const pylith::topology::Mesh& mesh,
+ const double upDir[3],
+ const double normalDir[3]);
+
+ /** Integrate contribution of cohesive cells to residual term.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateResidual(const pylith::topology::Field<pylith::topology::Mesh>& residual,
+ const double t,
+ pylith::topology::SolutionFields* const fields);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param jacobian Sparse matrix for Jacobian of system.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ void integrateJacobian(pylith::topology::Jacobian* jacobian,
+ const double t,
+ pylith::topology::SolutionFields* const fields);
+
+ /** Verify configuration is acceptable.
+ *
+ * @param mesh Finite-element mesh
+ */
+ void verifyConfiguration(const pylith::topology::Mesh& mesh) const;
+
+ /** Get vertex field associated with integrator.
+ *
+ * @param name Name of vertex field.
+ * @param fields Solution fields.
+ *
+ * @returns Vertex field.
+ */
+ const pylith::topology::Field<pylith::topology::SubMesh>&
+ vertexField(const char* name,
+ const pylith::topology::SolutionFields* fields =0);
+
+ /** Get cell field associated with integrator.
+ *
+ * @param name Name of cell field.
+ * @param fields Solution fields.
+ *
+ * @returns Cell field.
+ */
+ const pylith::topology::Field<pylith::topology::SubMesh>&
+ cellField(const char* name,
+ const pylith::topology::SolutionFields* fields =0);
+
+ }; // class FaultCohesiveTract
+
+ } // faults
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/trunk/modulesrc/faults/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/Makefile.am 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/modulesrc/faults/Makefile.am 2010-02-18 21:19:26 UTC (rev 16279)
@@ -27,9 +27,9 @@
EqKinSrc.i \
Fault.i \
FaultCohesive.i \
- FaultCohesiveDyn.i \
FaultCohesiveDynL.i \
FaultCohesiveKin.i \
+ FaultCohesiveTract.i \
../topology/SubMesh.i \
../feassemble/Quadrature.i \
../feassemble/Integrator.i
Modified: short/3D/PyLith/trunk/modulesrc/faults/faults.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/faults.i 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/modulesrc/faults/faults.i 2010-02-18 21:19:26 UTC (rev 16279)
@@ -24,9 +24,9 @@
#include "pylith/faults/EqKinSrc.hh"
#include "pylith/faults/Fault.hh"
#include "pylith/faults/FaultCohesive.hh"
-#include "pylith/faults/FaultCohesiveDyn.hh"
#include "pylith/faults/FaultCohesiveDynL.hh"
#include "pylith/faults/FaultCohesiveKin.hh"
+#include "pylith/faults/FaultCohesiveTract.hh"
#include "pylith/topology/SubMesh.hh"
#include "pylith/feassemble/Quadrature.hh"
@@ -73,9 +73,9 @@
%include "EqKinSrc.i"
%include "Fault.i"
%include "FaultCohesive.i"
-%include "FaultCohesiveDyn.i"
%include "FaultCohesiveDynL.i"
%include "FaultCohesiveKin.i"
+%include "FaultCohesiveTract.i"
// End of file
Modified: short/3D/PyLith/trunk/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/pylith/Makefile.am 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/pylith/Makefile.am 2010-02-18 21:19:26 UTC (rev 16279)
@@ -33,8 +33,8 @@
faults/Fault.py \
faults/FaultCohesive.py \
faults/FaultCohesiveKin.py \
- faults/FaultCohesiveDyn.py \
faults/FaultCohesiveDynL.py \
+ faults/FaultCohesiveTract.py \
faults/LiuCosSlipFn.py \
faults/SingleRupture.py \
faults/SlipTimeFn.py \
Deleted: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py 2010-02-18 21:19:26 UTC (rev 16279)
@@ -1,220 +0,0 @@
-#!/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
Copied: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveTract.py (from rev 16277, short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py)
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveTract.py (rev 0)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveTract.py 2010-02-18 21:19:26 UTC (rev 16279)
@@ -0,0 +1,220 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/faults/FaultCohesiveTract.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 FaultCohesiveTract as ModuleFaultCohesiveTract
+
+from pylith.utils.NullComponent import NullComponent
+
+# FaultCohesiveTract class
+class FaultCohesiveTract(FaultCohesive, Integrator, ModuleFaultCohesiveTract):
+ """
+ Python object for a fault surface with dynamic (friction) fault
+ implemented with cohesive elements.
+
+ Inventory
+
+ @class Inventory
+ Python object for managing FaultCohesiveTract 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.OutputFaultTract import OutputFaultTract
+ #output = pyre.inventory.facility("output", family="output_manager",
+ # factory=OutputFaultTract)
+ #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 = "CoTr "
+
+ 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)
+
+ ModuleFaultCohesiveTract.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)
+ ModuleFaultCohesiveTract.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 = ModuleFaultCohesiveTract.vertexField(self, name)
+ else:
+ field = ModuleFaultCohesiveTract.vertexField(self, name, fields)
+ return field
+
+
+ def getCellField(self, name, fields=None):
+ """
+ Get cell field.
+ """
+ if None == fields:
+ field = ModuleFaultCohesiveTract.cellField(self, name)
+ else:
+ field = ModuleFaultCohesiveTract.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):
+ ModuleFaultCohesiveTract.dbInitial(self, self.inventory.db)
+ #self.output = self.inventory.output
+ return
+
+
+ def _createModuleObj(self):
+ """
+ Create handle to C++ FaultCohesiveTract.
+ """
+ ModuleFaultCohesiveTract.__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 FaultCohesiveTract.
+ """
+ return FaultCohesiveTract()
+
+
+# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2010-02-18 18:22:21 UTC (rev 16278)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesive.cc 2010-02-18 21:19:26 UTC (rev 16279)
@@ -15,7 +15,7 @@
#include "TestFaultCohesive.hh" // Implementation of class methods
#include "pylith/faults/FaultCohesiveKin.hh" // USES FaultsCohesiveKin
-#include "pylith/faults/FaultCohesiveDyn.hh" // USES FaultsCohesiveDyn
+#include "pylith/faults/FaultCohesiveTract.hh" // USES FaultsCohesiveTract
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/utils/array.hh" // USES int_array, double_array
@@ -71,7 +71,7 @@
void
pylith::faults::TestFaultCohesive::testUseFaultMesh(void)
{ // testUseFaultMesh
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
CPPUNIT_ASSERT(!fault._useFaultMesh);
fault.useFaultMesh(true);
@@ -83,7 +83,7 @@
void
pylith::faults::TestFaultCohesive::testFaultMeshFilename(void)
{ // testFaultMeshFilename
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
CPPUNIT_ASSERT_EQUAL(std::string("fault.inp"), fault._faultMeshFilename);
const std::string filename = "SanAndreas.inp";
@@ -97,7 +97,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyLine2(void)
{ // testAdjustTopologyLine2
CohesiveDataLine2 data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyLine2
@@ -107,7 +107,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTri3(void)
{ // testAdjustTopologyTri3
CohesiveDataTri3 data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyTri3
@@ -117,7 +117,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTri3b(void)
{ // testAdjustTopologyTri3b
CohesiveDataTri3b data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyTri3b
@@ -127,7 +127,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTri3c(void)
{ // testAdjustTopologyTri3c
CohesiveDataTri3c data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyTri3c
@@ -137,7 +137,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTri3d(void)
{ // testAdjustTopologyTri3d
CohesiveDataTri3d data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyTri3d
@@ -147,7 +147,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTri3e(void)
{ // testAdjustTopologyTri3e
CohesiveDataTri3e data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyTri3e
@@ -157,7 +157,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTri3f(void)
{ // testAdjustTopologyTri3f
CohesiveDataTri3f data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyTri3f
@@ -167,7 +167,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4(void)
{ // testAdjustTopologyQuad4
CohesiveDataQuad4 data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyQuad4
@@ -177,7 +177,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4b(void)
{ // testAdjustTopologyQuad4b
CohesiveDataQuad4b data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyQuad4b
@@ -187,7 +187,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4c(void)
{ // testAdjustTopologyQuad4c
CohesiveDataQuad4c data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyQuad4c
@@ -197,7 +197,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4d(void)
{ // testAdjustTopologyQuad4d
CohesiveDataQuad4d data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyQuad4d
@@ -207,7 +207,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4e(void)
{ // testAdjustTopologyQuad4e
CohesiveDataQuad4e data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyQuad4e
@@ -217,7 +217,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4f(void)
{ // testAdjustTopologyQuad4f
CohesiveDataQuad4f data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyQuad4f
@@ -227,7 +227,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4g(void)
{ // testAdjustTopologyQuad4g
CohesiveDataQuad4g data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyQuad4g
@@ -237,8 +237,8 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyQuad4h(void)
{ // testAdjustTopologyQuad4h
CohesiveDataQuad4h data;
- FaultCohesiveDyn faultA;
- FaultCohesiveDyn faultB;
+ FaultCohesiveTract faultA;
+ FaultCohesiveTract faultB;
_testAdjustTopology(&faultA, &faultB, data, true, true);
} // testAdjustTopologyQuad4h
@@ -248,7 +248,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4(void)
{ // testAdjustTopologyTet4
CohesiveDataTet4 data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyTet4
@@ -258,7 +258,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4b(void)
{ // testAdjustTopologyTet4b
CohesiveDataTet4b data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyTet4b
@@ -268,7 +268,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4c(void)
{ // testAdjustTopologyTet4c
CohesiveDataTet4c data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyTet4c
@@ -278,7 +278,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4d(void)
{ // testAdjustTopologyTet4d
CohesiveDataTet4d data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyTet4d
@@ -288,7 +288,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4f(void)
{ // testAdjustTopologyTet4f
CohesiveDataTet4f data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyTet4f
@@ -298,7 +298,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4g(void)
{ // testAdjustTopologyTet4g
CohesiveDataTet4g data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyTet4g
@@ -308,7 +308,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4h(void)
{ // testAdjustTopologyTet4h
CohesiveDataTet4h data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyTet4h
@@ -318,7 +318,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4i(void)
{ // testAdjustTopologyTet4i
CohesiveDataTet4i data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyTet4i
@@ -328,7 +328,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyTet4j(void)
{ // testAdjustTopologyTet4j
CohesiveDataTet4j data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyTet4j
@@ -338,7 +338,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8(void)
{ // testAdjustTopologyHex8
CohesiveDataHex8 data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyHex8
@@ -348,7 +348,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8b(void)
{ // testAdjustTopologyHex8b
CohesiveDataHex8b data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyHex8b
@@ -358,7 +358,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8c(void)
{ // testAdjustTopologyHex8c
CohesiveDataHex8c data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyHex8c
@@ -368,7 +368,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8d(void)
{ // testAdjustTopologyHex8d
CohesiveDataHex8d data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyHex8d
@@ -378,7 +378,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8e(void)
{ // testAdjustTopologyHex8e
CohesiveDataHex8e data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyHex8e
@@ -388,7 +388,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8f(void)
{ // testAdjustTopologyHex8f
CohesiveDataHex8f data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyHex8f
@@ -398,7 +398,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8g(void)
{ // testAdjustTopologyHex8g
CohesiveDataHex8g data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, true);
} // testAdjustTopologyHex8g
@@ -408,7 +408,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8h(void)
{ // testAdjustTopologyHex8h
CohesiveDataHex8h data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyHex8h
@@ -418,7 +418,7 @@
pylith::faults::TestFaultCohesive::testAdjustTopologyHex8i(void)
{ // testAdjustTopologyHex8i
CohesiveDataHex8i data;
- FaultCohesiveDyn fault;
+ FaultCohesiveTract fault;
_testAdjustTopology(&fault, data, false);
} // testAdjustTopologyHex8i
More information about the CIG-COMMITS
mailing list