[cig-commits] r15030 - in short/3D/PyLith/trunk: . libsrc/bc libsrc/faults libsrc/feassemble libsrc/materials libsrc/topology modulesrc/feassemble modulesrc/materials modulesrc/topology unittests/libtests/bc unittests/libtests/feassemble unittests/libtests/materials unittests/pytests/faults
brad at geodynamics.org
brad at geodynamics.org
Thu May 21 20:10:51 PDT 2009
Author: brad
Date: 2009-05-21 20:10:47 -0700 (Thu, 21 May 2009)
New Revision: 15030
Added:
short/3D/PyLith/trunk/libsrc/bc/PointForce.cc
short/3D/PyLith/trunk/libsrc/bc/PointForce.hh
short/3D/PyLith/trunk/libsrc/bc/PointForce.icc
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh
short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.cc
short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.hh
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
short/3D/PyLith/trunk/libsrc/bc/bcfwd.hh
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh
short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.cc
short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh
short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.cc
short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.hh
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/materials/Material.hh
short/3D/PyLith/trunk/libsrc/topology/Fields.hh
short/3D/PyLith/trunk/libsrc/topology/Fields.icc
short/3D/PyLith/trunk/modulesrc/feassemble/Quadrature.i
short/3D/PyLith/trunk/modulesrc/materials/Material.i
short/3D/PyLith/trunk/modulesrc/topology/Fields.i
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
Log:
Consolidated Field data members into Fields (makes memory logging easier).
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/TODO 2009-05-22 03:10:47 UTC (rev 15030)
@@ -7,12 +7,14 @@
field split
Brad
- Test precompute geometry
+ Test precompute geometry [DONE]
Nondimensionalization
cleanup
full-scale testing
test uniform refinement
+ time-dependent BC
+
Charles
3-D Power-law rheology
C++
@@ -46,8 +48,6 @@
2. Nondimensionalization
- Move nondimensionalization of coordinates to after distribution.
-
Add Quadrature::computeGeometry(coordinatesCell)
Use #define PRECOMPUTE_GEOMETRY [TEST]
@@ -94,6 +94,9 @@
libtests/topology/TestMesh::testNondimensionalize()
+ libtests/topology/TestFields::testHasField()
+ pytests/topology/TestFields.test_hasField()
+
libtests/topology/Field add constraints to field in unit tests
copy
+=
@@ -171,6 +174,7 @@
after assert(0) to insure error is trapped.
Use Fields object to hold Field
+ Material
Eliminate use of Inventory class.
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -42,7 +42,7 @@
// Default constructor.
pylith::bc::AbsorbingDampers::AbsorbingDampers(void) :
_boundaryMesh(0),
- _dampingConsts(0)
+ _parameters(0)
{ // constructor
} // constructor
@@ -51,7 +51,7 @@
pylith::bc::AbsorbingDampers::~AbsorbingDampers(void)
{ // destructor
delete _boundaryMesh; _boundaryMesh = 0;
- delete _dampingConsts; _dampingConsts = 0;
+ delete _parameters; _parameters = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -65,7 +65,7 @@
assert(0 != _db);
delete _boundaryMesh; _boundaryMesh = 0;
- delete _dampingConsts; _dampingConsts = 0;
+ delete _parameters; _parameters = 0;
_boundaryMesh = new topology::SubMesh(mesh, _label.c_str());
assert(0 != _boundaryMesh);
@@ -120,11 +120,14 @@
const int spaceDim = cellGeometry.spaceDim();
const int fiberDim = numQuadPts * spaceDim;
- _dampingConsts = new topology::Field<topology::SubMesh>(*_boundaryMesh);
- _dampingConsts->label("damping constants");
- assert(0 != _dampingConsts);
- _dampingConsts->newSection(cells, fiberDim);
- _dampingConsts->allocate();
+ _parameters =
+ new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
+ assert(0 != _parameters);
+ _parameters->add("damping constants", "damping_constants");
+ topology::Field<topology::SubMesh>& dampingConsts =
+ _parameters->get("damping constants");
+ dampingConsts.newSection(cells, fiberDim);
+ dampingConsts.allocate();
// Containers for orientation information
const int orientationSize = spaceDim * spaceDim;
@@ -170,7 +173,7 @@
const double velocityScale =
_normalizer->lengthScale() / _normalizer->timeScale();
- const ALE::Obj<SubRealSection>& dampersSection = _dampingConsts->section();
+ const ALE::Obj<SubRealSection>& dampersSection = dampingConsts.section();
assert(!dampersSection.isNull());
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
@@ -267,7 +270,7 @@
{ // integrateResidual
assert(0 != _quadrature);
assert(0 != _boundaryMesh);
- assert(0 != _dampingConsts);
+ assert(0 != _parameters);
assert(0 != fields);
// Get cell geometry information that doesn't depend on cell
@@ -291,7 +294,8 @@
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<SubRealSection>& dampersSection = _dampingConsts->section();
+ const ALE::Obj<SubRealSection>& dampersSection =
+ _parameters->get("damping constants").section();
assert(!dampersSection.isNull());
const ALE::Obj<RealSection>& residualSection = residual.section();
@@ -410,7 +414,8 @@
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<SubRealSection>& dampersSection = _dampingConsts->section();
+ const ALE::Obj<SubRealSection>& dampersSection =
+ _parameters->get("damping constants").section();
assert(!dampersSection.isNull());
const topology::Field<topology::Mesh>& solution = fields->solution();
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -125,8 +125,8 @@
/// Mesh of absorbing boundary
topology::SubMesh* _boundaryMesh;
- /// Damping constants in global coordinates at integration points.
- topology::Field<topology::SubMesh>* _dampingConsts;
+ /// Parameters for damping constants.
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
}; // class AbsorbingDampers
Modified: short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -17,6 +17,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/Fields.hh" // USES Fields
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -29,7 +30,7 @@
// Default constructor.
pylith::bc::DirichletBoundary::DirichletBoundary(void) :
_boundaryMesh(0),
- _tmpField(0)
+ _fields(0)
{ // constructor
} // constructor
@@ -38,7 +39,7 @@
pylith::bc::DirichletBoundary::~DirichletBoundary(void)
{ // destructor
delete _boundaryMesh; _boundaryMesh = 0;
- delete _tmpField; _tmpField = 0;
+ delete _fields; _fields = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -87,20 +88,27 @@
const int numPoints = _points.size();
const int numFixedDOF = _fixedDOF.size();
- if (0 == _tmpField) {
- _tmpField = new topology::Field<topology::SubMesh>(*_boundaryMesh);
- assert(0 != _tmpField);
- _tmpField->newSection(vertices, fiberDim);
- _tmpField->allocate();
+ if (0 == _fields) {
+ _fields =
+ new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
+ assert(0 != _fields);
+ _fields->add("output buffer", "temporary_buffer");
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("output buffer");
+ buffer.newSection(vertices, fiberDim);
+ buffer.allocate();
+ buffer.addDimensionOkay(true);
+ buffer.vectorFieldType(topology::Field<topology::SubMesh>::VECTOR);
} // if
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("output buffer");
if (0 == strcasecmp(name, "initial")) {
- _tmpField->label("displacement");
- _tmpField->vectorFieldType(topology::Field<topology::SubMesh>::VECTOR);
- _tmpField->scale(_normalizer->lengthScale());
- _tmpField->addDimensionOkay(true);
- _tmpField->zero();
- const ALE::Obj<RealSection>& section = _tmpField->section();
+ buffer.label("bc_displacement");
+ buffer.vectorFieldType(topology::Field<topology::SubMesh>::VECTOR);
+ buffer.scale(_normalizer->lengthScale());
+ buffer.zero();
+ const ALE::Obj<RealSection>& section = buffer.section();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
@@ -110,12 +118,14 @@
section->updatePointAll(_points[iPoint], &values[0]);
} // for
} else if (0 == strcasecmp(name, "rate-of-change")) {
- _tmpField->label("velocity");
- _tmpField->vectorFieldType(topology::Field<topology::SubMesh>::VECTOR);
- _tmpField->scale(_normalizer->lengthScale());
- _tmpField->addDimensionOkay(true);
- _tmpField->zero();
- const ALE::Obj<RealSection>& section = _tmpField->section();
+ buffer.label("bc_velocity");
+ buffer.vectorFieldType(topology::Field<topology::SubMesh>::VECTOR);
+ assert(0 != _normalizer->timeScale());
+ const double velocityScale =
+ _normalizer->lengthScale() / _normalizer->timeScale();
+ buffer.scale(velocityScale);
+ buffer.zero();
+ const ALE::Obj<RealSection>& section = buffer.section();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
@@ -132,7 +142,7 @@
throw std::runtime_error(msg.str());
} // else
- return *_tmpField;
+ return buffer;
} // getVertexField
Modified: short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/bc/DirichletBoundary.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -74,8 +74,10 @@
private :
topology::SubMesh* _boundaryMesh; ///< Boundary mesh.
- topology::Field<topology::SubMesh>* _tmpField; ///< Temporary field for output.
+ /// Fields manager (holds temporary field for output).
+ topology::Fields<topology::Field<topology::SubMesh> >* _fields;
+
}; // class DirichletBoundary
#include "DirichletBoundary.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -39,7 +39,7 @@
// Default constructor.
pylith::bc::Neumann::Neumann(void) :
_boundaryMesh(0),
- _tractions(0)
+ _parameters(0)
{ // constructor
} // constructor
@@ -48,7 +48,7 @@
pylith::bc::Neumann::~Neumann(void)
{ // destructor
delete _boundaryMesh; _boundaryMesh = 0;
- delete _tractions; _tractions = 0;
+ delete _parameters; _parameters = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -62,7 +62,7 @@
assert(0 != _db);
delete _boundaryMesh; _boundaryMesh = 0;
- delete _tractions; _tractions = 0;
+ delete _parameters; _parameters = 0;
_boundaryMesh = new topology::SubMesh(mesh, _label.c_str());
assert(0 != _boundaryMesh);
@@ -116,11 +116,13 @@
const int spaceDim = cellGeometry.spaceDim();
const int fiberDim = spaceDim * numQuadPts;
- _tractions = new topology::Field<topology::SubMesh>(*_boundaryMesh);
- _tractions->label("tractions");
- assert(0 != _tractions);
- _tractions->newSection(cells, fiberDim);
- _tractions->allocate();
+ _parameters =
+ new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
+ assert(0 != _parameters);
+ _parameters->add("traction", "traction");
+ topology::Field<topology::SubMesh>& traction = _parameters->get("traction");
+ traction.newSection(cells, fiberDim);
+ traction.allocate();
// Containers for orientation information
const int orientationSize = spaceDim * spaceDim;
@@ -177,8 +179,9 @@
coordinatesCell.size(),
&coordinatesCell[0]);
- const ALE::Obj<SubRealSection>& tractSection = _tractions->section();
- assert(!tractSection.isNull());
+ const ALE::Obj<SubRealSection>& tractionSection =
+ _parameters->get("traction").section();
+ assert(!tractionSection.isNull());
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
@@ -252,7 +255,7 @@
} // for
// Update tractionsGlobal
- tractSection->updatePoint(*c_iter, &cellTractionsGlobal[0]);
+ tractionSection->updatePoint(*c_iter, &cellTractionsGlobal[0]);
} // for
_db->close();
@@ -268,7 +271,7 @@
{ // integrateResidual
assert(0 != _quadrature);
assert(0 != _boundaryMesh);
- assert(0 != _tractions);
+ assert(0 != _parameters);
// Get cell geometry information that doesn't depend on cell
const int numQuadPts = _quadrature->numQuadPts();
@@ -290,8 +293,9 @@
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<SubRealSection>& tractSection = _tractions->section();
- assert(!tractSection.isNull());
+ const ALE::Obj<SubRealSection>& tractionSection =
+ _parameters->get("traction").section();
+ assert(!tractionSection.isNull());
const ALE::Obj<RealSection>& residualSection = residual.section();
topology::SubMesh::UpdateAddVisitor residualVisitor(*residualSection,
&_cellVector[0]);
@@ -320,7 +324,7 @@
_resetCellVector();
// Restrict tractions to cell
- tractSection->restrictPoint(*c_iter,
+ tractionSection->restrictPoint(*c_iter,
&tractionsCell[0], tractionsCell.size());
// Get cell geometry information that depends on cell
@@ -382,11 +386,11 @@
pylith::bc::Neumann::cellField(const char* name,
topology::SolutionFields* const fields)
{ // cellField
- assert(0 != _tractions);
+ assert(0 != _parameters);
assert(0 != name);
if (0 == strcasecmp(name, "tractions")) {
- return *_tractions;
+ return _parameters->get("traction");;
} else {
std::ostringstream msg;
msg << "Unknown field '" << name << "' requested for Neumann BC '"
@@ -394,7 +398,7 @@
throw std::runtime_error(msg.str());
} // else
- return *_tractions;
+ return _parameters->get("traction"); // Satisfy method definition
} // cellField
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -113,8 +113,9 @@
/// Mesh over which tractions are applied
topology::SubMesh* _boundaryMesh;
- /// Traction vector in global coordinates at integration points.
- topology::Field<topology::SubMesh>* _tractions;
+ /// Parameters for tractions vector in global coordinates at
+ /// integration points.
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
}; // class Neumann
Added: short/3D/PyLith/trunk/libsrc/bc/PointForce.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/PointForce.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/PointForce.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -0,0 +1,244 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "PointForce.hh" // implementation of object methods
+
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include <cstring> // USES strcpy()
+#include <cassert> // USES assert()
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::bc::PointForce::PointForce(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::bc::PointForce::~PointForce(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set indices of vertices with point forces.
+void
+pylith::bc::PointForce::forceDOF(const int* flags,
+ const int size)
+{ // forceDOF
+ if (size > 0)
+ assert(0 != flags);
+
+ _forceDOF.resize(size);
+ for (int i=0; i < size; ++i)
+ _forceDOF[i] = flags[i];
+} // forceDOF
+
+// ----------------------------------------------------------------------
+// Initialize boundary condition.
+void
+pylith::bc::PointForce::initialize(const topology::Mesh& mesh,
+ const double upDir[3])
+{ // initialize
+ const int numForceDOF = _forceDOF.size();
+ if (0 == numForceDOF)
+ return;
+
+ _getPoints(mesh);
+ _setupQueryDatabases();
+ _queryDatabases(mesh);
+} // initialize
+
+// ----------------------------------------------------------------------
+// Integrate contributions to residual term (r) for operator that
+// do not require assembly over cells, vertices, or processors.
+void
+integrateResidualAssembled(const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
+{ // integrateResidualAssembled
+ const int numPoints = _points.size();
+ const int numForceDOF = _forceDOF.size();
+
+ const ALE::Obj<RealSection>& amplitudeSection =
+ _parameters->get("amplitude").section();
+ assert(!amplitudeSection.isNull());
+
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ assert(!residualSection.isNull());
+
+ const ALE::Obj<RealSection>& startTimeSection = (0 == _dbStartTime) ?
+ 0 : _parameters->get("start time").section();
+
+ double_array forcesVertex(spaceDim);
+ double scale = 1.0;
+ for (int iPoint=0; iPoint < numPoints; ++iPoint) {
+ const int p_force = _points[iPoint];
+ assert(numForceDOF == amplitudeSection->getFiberDimension(p_force));
+ const double* amplitudeVertex = amplitudeSection->restrictPoint(p_force);
+ assert(0 != amplitudeVertex);
+
+ if (0 != _dbStartTime) {
+ assert(!startTimeSection.isNull()); // Expect section with start times
+ assert(0 != _dbTimeAmp); // Expect database with temporal evolution
+
+ assert(1 == startTimeSection->getFiberDimension(p_force));
+ const double tRef = *startTimeSection->restrictPoint(p_force);
+ const double tRel = t - tRef;
+ if (tRel > 0.0) {
+ _normalizer->dimensionalize(&tRel, 1, timeScale);
+ scale = _dbTimeAmp->query(tRel);
+ } // if
+ } // if
+
+ for (int iDOF=0; iDOF < numForceDOF; ++iDOF)
+ forcesVertex[_forceDOF[iDOF]] = amplitudeVertex[iDOF]*scale;
+ } // for
+} // interateResidualAssembled
+
+// ----------------------------------------------------------------------
+// Verify configuration is acceptable.
+void
+verifyConfiguration(const topology::Mesh& mesh) const
+{ // verifyConfiguration
+} // verifyConfiguration
+
+// ----------------------------------------------------------------------
+// Get mesh labels for points associated with point forces.
+void
+pylith::bc::PointForce::_getPoints(const topology::Mesh& mesh)
+{ // _getPoints
+ typedef topology::Mesh::IntSection::chart_type chart_type;
+
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ const ALE::Obj<SieveMesh::int_section_type>& groupField =
+ sieveMesh->getIntSection(_label);
+ if (groupField.isNull()) {
+ std::ostringstream msg;
+ msg << "Could not find group of points '" << _label << "' in mesh.";
+ throw std::runtime_error(msg.str());
+ } // if
+ assert(!groupField.isNull());
+ const chart_type& chart = groupField->getChart();
+ const chart_type::const_iterator& chartEnd = chart.end();
+ const int numPoints = groupField->size();
+ _points.resize(numPoints);
+ int i = 0;
+ for(chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter)
+ if (groupField->getFiberDimension(*c_iter))
+ _points[i++] = *c_iter;
+} // _getPoints
+
+// ----------------------------------------------------------------------
+// Setup initial and rate of change databases for querying.
+void
+pylith::bc::PointForce::_setupQueryDatabases(void)
+{ // _setupQueryDatabases
+ assert(0 != _db);
+
+ const int numForceDOF = _forceDOF.size();
+ char** valueNames = (numForceDOF > 0) ? new char*[numForceDOF] : 0;
+ for (int i=0; i < numForceDOF; ++i) {
+ std::ostringstream name;
+ name << "dof-" << _forceDOF[i];
+ const int size = 1 + name.str().length();
+ valueNames[i] = new char[size];
+ strcpy(valueNames[i], name.str().c_str());
+ } // for
+
+ // Setup initial database.
+ _db->open();
+ _db->queryVals(const_cast<const char**>(valueNames), numForceDOF);
+
+ // Setup rate database, if provided.
+ if (0 != _dbRate) {
+ _dbRate->open();
+ _dbRate->queryVals((const char**) valueNames, numForceDOF);
+ } // if
+ for (int i=0; i < numForceDOF; ++i) {
+ delete[] valueNames[i]; valueNames[i] = 0;
+ } // for
+ delete[] valueNames; valueNames = 0;
+} // _setupQueryDatabases
+
+// ----------------------------------------------------------------------
+// Query initial and rate of change databases for values.
+void
+pylith::bc::PointForce::_queryDatabases(const topology::Mesh& mesh)
+{ // _queryDatabases
+ assert(0 != _db);
+
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ const ALE::Obj<RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ assert(0 != cs);
+ const int spaceDim = cs->spaceDim();
+
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+ const double velocityScale =
+ _normalizer->lengthScale() / _normalizer->timeScale();
+
+ const int numPoints = _points.size();
+ const int numForceDOF = _forceDOF.size();
+ _valuesInitial.resize(numPoints*numForceDOF);
+ if (0 != _dbRate)
+ _valuesRate.resize(numPoints*numForceDOF);
+
+ double_array queryValues(numForceDOF);
+ double_array vCoordsGlobal(spaceDim);
+ for (int iPoint=0; iPoint < numPoints; ++iPoint) {
+ // Get coordinates of vertex
+ coordinates->restrictPoint(_points[iPoint],
+ &vCoordsGlobal[0], vCoordsGlobal.size());
+ _normalizer->dimensionalize(&vCoordsGlobal[0], vCoordsGlobal.size(),
+ lengthScale);
+ int err = _db->query(&queryValues[0], numForceDOF,
+ &vCoordsGlobal[0], vCoordsGlobal.size(), cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find initial values at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << vCoordsGlobal[i];
+ msg << ") using spatial database " << _db->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ for (int iDOF=0; iDOF < numForceDOF; ++iDOF)
+ _valuesInitial[numForceDOF*iPoint+iDOF] =
+ _normalizer->nondimensionalize(queryValues[iDOF], lengthScale);
+ } // for
+ _db->close();
+} // _queryDatabases
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/bc/PointForce.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/PointForce.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/PointForce.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -0,0 +1,130 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/bc/PointForce.hh
+ *
+ * @brief C++ implementation of point force on vertices.
+ */
+
+#if !defined(pylith_bc_pointforce_hh)
+#define pylith_bc_pointforce_hh
+
+// Include directives ---------------------------------------------------
+#include "BoundaryCondition.hh" // ISA BoundaryCondition
+#include "pylith/feassemble/Integrator.hh" // ISA Integrator
+
+#include "pylith/topology/topologyfwd.hh" // USES Fields<Mesh>
+#include "pylith/utils/array.hh" // HASA int_array
+
+// PointForce ------------------------------------------------------
+class pylith::bc::PointForce : public BoundaryCondition,
+ public feassemble::Integrator
+{ // class PointForce
+ friend class TestPointForce; // unit testing
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ PointForce(void);
+
+ /// Destructor.
+ ~PointForce(void);
+
+ /** Set indices of vertices with point forces.
+ *
+ * Note: forces at all points are applied to the same degrees of freedom.
+ *
+ * Example: [0, 1] to apply forces to x and y degrees of freedom in
+ * Cartesian system.
+ *
+ * @param flags Array of indices for degrees of freedom for forces.
+ * @param size Size of array
+ */
+ void forceDOF(const int* flags,
+ const int size);
+
+ /** Initialize boundary condition.
+ *
+ * @param mesh PETSc mesh
+ * @param upDir Vertical direction (somtimes used in 3-D problems).
+ */
+ void initialize(const topology::Mesh& mesh,
+ const double upDir[3]);
+
+ /** Integrate contributions to residual term (r) for operator.
+ *
+ * @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);
+
+ /** Verify configuration is acceptable.
+ *
+ * @param mesh Finite-element mesh
+ */
+ void verifyConfiguration(const topology::Mesh& mesh) const;
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Get mesh labels for points associated with applied forces.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ void _getPoints(const topology::Mesh& mesh);
+
+ /// Setup databases for querying for point forces.
+ void _setupQueryDatabases(void);
+
+ /** Query databases for point forces.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ void _queryDatabases(const topology::Mesh& mesh);
+
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ /// Not implemented
+ PointForce(const PointForce& m);
+
+ /// Not implemented
+ const PointForce& operator=(const PointForce& m);
+
+ // PROTECTED MEMBERS //////////////////////////////////////////////////
+protected :
+
+ /// Parameters for point forces.
+ topology::Fields<topology::Field<topology::Mesh> >* _parameters;
+
+ /// Start time for point forces.
+ spatialdata::spatialdb::SpatialDB* _dbStartTime;
+
+ /// Temporal evolution for amplitude of point forces.
+ spatialdata::spatialdb::SpatialDB* _dbTimeAmp;
+
+
+ int_array _points; ///< Points for forces.
+ int_array _forceDOF; ///< Indices of degrees of freedom with forces.
+
+}; // class PointForce
+
+#include "PointForce.icc" // inline methods
+
+#endif // pylith_bc_pointforce_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/bc/PointForce.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/PointForce.icc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/bc/PointForce.icc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -0,0 +1,18 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#if !defined(pylith_bc_pointforce_hh)
+#error "PointForce.icc can only be included from PointForce.hh"
+#endif
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/bc/bcfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/bcfwd.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/bc/bcfwd.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -29,6 +29,7 @@
class DirichletBoundary;
class Neumann;
class AbsorbingDampers;
+ class PointForce;
} // bc
} // pylith
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -36,7 +36,6 @@
pylith::faults::BruneSlipFn::BruneSlipFn(void) :
_slipTimeVertex(0),
_riseTimeVertex(0),
- _parameters(0),
_dbFinalSlip(0),
_dbSlipTime(0),
_dbRiseTime(0)
@@ -47,7 +46,6 @@
// Destructor.
pylith::faults::BruneSlipFn::~BruneSlipFn(void)
{ // destructor
- delete _parameters; _parameters = 0;
_dbFinalSlip = 0; // :TODO: Use shared pointer.
_dbSlipTime = 0; // :TODO: Use shared pointer.
_dbRiseTime = 0; // :TODO: Use shared pointer.
Modified: short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/BruneSlipFn.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -137,10 +137,6 @@
double _riseTimeVertex; ///< Rise time at a vertex.
double_array _slipVertex; ///< Slip at a vertex.
- /// Parameters for Brune slip time function, final slip (vector),
- /// rise time (scalar), slip time (scalar).
- topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
-
/// Spatial database for final slip.
spatialdata::spatialdb::SpatialDB* _dbFinalSlip;
Modified: short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -35,7 +35,6 @@
// Default constructor.
pylith::faults::ConstRateSlipFn::ConstRateSlipFn(void) :
_slipTimeVertex(0),
- _parameters(0),
_dbSlipRate(0),
_dbSlipTime(0)
{ // constructor
@@ -45,7 +44,6 @@
// Destructor.
pylith::faults::ConstRateSlipFn::~ConstRateSlipFn(void)
{ // destructor
- delete _parameters; _parameters = 0;
_dbSlipRate = 0; // :TODO: Use shared pointer.
_dbSlipTime = 0; // :TODO: Use shared pointer.
} // destructor
Modified: short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/ConstRateSlipFn.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -113,10 +113,6 @@
double _slipTimeVertex; ///< Slip time at a vertex.
double_array _slipRateVertex; ///< Slip rate at a vertex.
- /// Parameters for constant slip rate slip time function, slip rate
- /// (vector) and slip time (scalar).
- topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
-
/// Spatial database for slip rate.
spatialdata::spatialdb::SpatialDB* _dbSlipRate;
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesive.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -70,12 +70,6 @@
virtual
bool _useLagrangeConstraints(void) const = 0;
- // NOT IMPLEMENTED ////////////////////////////////////////////////////
-private :
-
- FaultCohesive(const FaultCohesive&); ///< Not implemented
- const FaultCohesive& operator=(const FaultCohesive&); ///< Not implemented
-
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
@@ -85,6 +79,12 @@
std::string _faultMeshFilename; /// Filename for fault mesh UCD file.
+ // NOT IMPLEMENTED ////////////////////////////////////////////////////
+private :
+
+ FaultCohesive(const FaultCohesive&); ///< Not implemented
+ const FaultCohesive& operator=(const FaultCohesive&); ///< Not implemented
+
}; // class FaultCohesive
#endif // pylith_faults_faultcohesive_hh
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -135,10 +135,6 @@
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
- /// Orientation of fault surface at vertices (fiber dimension is
- /// nonzero only at constraint vertices)
- topology::Field<topology::SubMesh>* _orientation;
-
}; // class FaultCohesiveDyn
#include "FaultCohesiveDyn.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -49,9 +49,7 @@
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveKin::FaultCohesiveKin(void) :
- _fields(0),
- _bufferVectorField(0),
- _bufferScalarField(0)
+ _fields(0)
{ // constructor
} // constructor
@@ -60,8 +58,6 @@
pylith::faults::FaultCohesiveKin::~FaultCohesiveKin(void)
{ // destructor
delete _fields; _fields = 0;
- delete _bufferVectorField; _bufferVectorField = 0;
- delete _bufferScalarField; _bufferScalarField = 0;
// :TODO: Use shared pointers for earthquake sources
} // destructor
@@ -715,10 +711,12 @@
orientationSection->getFibration(0);
assert(!dirSection.isNull());
_allocateBufferVectorField();
- assert(0 != _bufferVectorField);
- _bufferVectorField->copy(dirSection);
- _bufferVectorField->label("strike_dir");
- return *_bufferVectorField;
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ buffer.copy(dirSection);
+ buffer.label("strike_dir");
+ buffer.scale(1.0);
+ return buffer;
} else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
const ALE::Obj<RealSection>& orientationSection =
@@ -727,10 +725,12 @@
const ALE::Obj<RealSection>& dirSection =
orientationSection->getFibration(1);
_allocateBufferVectorField();
- assert(0 != _bufferVectorField);
- _bufferVectorField->copy(dirSection);
- _bufferVectorField->label("dip_dir");
- return *_bufferVectorField;
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ buffer.copy(dirSection);
+ buffer.label("dip_dir");
+ buffer.scale(1.0);
+ return buffer;
} else if (0 == strcasecmp("normal_dir", name)) {
const ALE::Obj<RealSection>& orientationSection =
@@ -742,10 +742,12 @@
orientationSection->getFibration(space);
assert(!dirSection.isNull());
_allocateBufferVectorField();
- assert(0 != _bufferVectorField);
- _bufferVectorField->copy(dirSection);
- _bufferVectorField->label("normal_dir");
- return *_bufferVectorField;
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ buffer.copy(dirSection);
+ buffer.label("normal_dir");
+ buffer.scale(1.0);
+ return buffer;
} else if (0 == strncasecmp("final_slip_X", name, slipStrLen)) {
const std::string value = std::string(name).substr(slipStrLen+1);
@@ -764,9 +766,10 @@
assert(0 != fields);
const topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
_allocateBufferVectorField();
- _calcTractionsChange(_bufferVectorField, dispT);
- _bufferVectorField->label("traction_change");
- return *_bufferVectorField;
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ _calcTractionsChange(&buffer, dispT);
+ return buffer;
} else {
std::ostringstream msg;
@@ -775,8 +778,12 @@
throw std::runtime_error(msg.str());
} // else
- assert(0 != _bufferScalarField);
- return *_bufferScalarField;
+
+ // Satisfy return values
+ assert(0 != _fields);
+ const topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ return buffer;
} // vertexField
// ----------------------------------------------------------------------
@@ -792,9 +799,11 @@
<< "' for fault '" << label() << ".";
throw std::runtime_error(msg.str());
- // Return generic section to satisfy member function definition.
- assert(0 != _bufferScalarField);
- return *_bufferScalarField;
+ // Satisfy return values
+ assert(0 != _fields);
+ const topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
+ return buffer;
} // cellField
// ----------------------------------------------------------------------
@@ -1154,7 +1163,11 @@
assert(0 != tractions);
assert(0 != _faultMesh);
assert(0 != _fields);
+ assert(0 != _normalizer);
+ tractions->scale(_normalizer->pressureScale());
+ tractions->label("traction");
+
// Get vertices from mesh of domain.
const ALE::Obj<SieveMesh>& sieveMesh = dispT.mesh().sieveMesh();
assert(!sieveMesh.isNull());
@@ -1265,19 +1278,20 @@
void
pylith::faults::FaultCohesiveKin::_allocateBufferVectorField(void)
{ // _allocateBufferVectorField
- if (0 != _bufferVectorField)
+ assert(0 != _fields);
+ if (_fields->hasField("buffer (vector)"))
return;
// Create vector field; use same shape/chart as cumulative slip field.
assert(0 != _faultMesh);
- assert(0 != _fields);
- _bufferVectorField = new topology::Field<topology::SubMesh>(*_faultMesh);
- _bufferVectorField->label("vector field buffer");
+ _fields->add("buffer (vector)", "buffer");
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (vector)");
const topology::Field<topology::SubMesh>& slip =
_fields->get("cumulative slip");
- _bufferVectorField->newSection(slip);
- _bufferVectorField->allocate();
- _bufferVectorField->zero();
+ buffer.newSection(slip);
+ buffer.allocate();
+ buffer.zero();
} // _allocateBufferVectorField
// ----------------------------------------------------------------------
@@ -1285,18 +1299,19 @@
void
pylith::faults::FaultCohesiveKin::_allocateBufferScalarField(void)
{ // _allocateBufferScalarField
- if (0 != _bufferScalarField)
+ assert(0 != _fields);
+ if (_fields->hasField("buffer (scalar)"))
return;
// Create vector field; use same shape/chart as area field.
assert(0 != _faultMesh);
- assert(0 != _fields);
- _bufferScalarField = new topology::Field<topology::SubMesh>(*_faultMesh);
- _bufferScalarField->label("scalar field buffer");
+ _fields->add("buffer (scalar)", "buffer");
+ topology::Field<topology::SubMesh>& buffer =
+ _fields->get("buffer (scalar)");
const topology::Field<topology::SubMesh>& area = _fields->get("area");
- _bufferScalarField->newSection(area);
- _bufferScalarField->allocate();
- _bufferScalarField->zero();
+ buffer.newSection(area);
+ buffer.allocate();
+ buffer.zero();
} // _allocateBufferScalarField
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -237,12 +237,6 @@
/// Fields for fault information.
topology::Fields<topology::Field<topology::SubMesh> >* _fields;
- /// Buffer for vector field over fault vertices.
- topology::Field<topology::SubMesh>* _bufferVectorField;
-
- /// Buffer for scalar field over fault vertices.
- topology::Field<topology::SubMesh>* _bufferScalarField;
-
/// Map label of cohesive cell to label of cells in fault mesh.
std::map<topology::Mesh::SieveMesh::point_type,
topology::SubMesh::SieveMesh::point_type> _cohesiveToFault;
Modified: short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -36,7 +36,6 @@
pylith::faults::LiuCosSlipFn::LiuCosSlipFn(void) :
_slipTimeVertex(0),
_riseTimeVertex(0),
- _parameters(0),
_dbFinalSlip(0),
_dbSlipTime(0),
_dbRiseTime(0)
@@ -47,7 +46,6 @@
// Destructor.
pylith::faults::LiuCosSlipFn::~LiuCosSlipFn(void)
{ // destructor
- delete _parameters; _parameters = 0;
_dbFinalSlip = 0; // :TODO: Use shared pointer.
_dbSlipTime = 0; // :TODO: Use shared pointer.
_dbRiseTime = 0; // :TODO: Use shared pointer.
Modified: short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/LiuCosSlipFn.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -137,10 +137,6 @@
double _riseTimeVertex; ///< Rise time at a vertex.
double_array _slipVertex; ///< Slip at a vertex.
- /// Parameters for Liu cosine/sine slip time function, final slip
- /// (vector), slip time (scalar), rise time (scalar).
- topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
-
/// Spatial database for final slip.
spatialdata::spatialdb::SpatialDB* _dbFinalSlip;
Modified: short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -14,9 +14,14 @@
#include "SlipTimeFn.hh" // implementation of object methods
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/Field.hh" // USES Field
+
// ----------------------------------------------------------------------
// Default constructor.
-pylith::faults::SlipTimeFn::SlipTimeFn(void)
+pylith::faults::SlipTimeFn::SlipTimeFn(void) :
+ _parameters(0)
{ // constructor
} // constructor
@@ -24,6 +29,7 @@
// Destructor.
pylith::faults::SlipTimeFn::~SlipTimeFn(void)
{ // destructor
+ delete _parameters; _parameters = 0;
} // destructor
Modified: short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/SlipTimeFn.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -23,7 +23,7 @@
// Include directives ---------------------------------------------------
#include "faultsfwd.hh" // forward declarations
-#include "pylith/topology/topologyfwd.hh" // USES Field<SubMesh>
+#include "pylith/topology/topologyfwd.hh" // USES Fields<SubMesh>
#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
#include "spatialdata/spatialdb/spatialdbfwd.hh" // USES SpatialDB
@@ -93,6 +93,12 @@
virtual
const topology::Field<topology::SubMesh>& slipTime(void) = 0;
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ /// Parameters for slip time function.
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -35,7 +35,6 @@
// Default constructor.
pylith::faults::StepSlipFn::StepSlipFn(void) :
_slipTimeVertex(0),
- _parameters(0),
_dbFinalSlip(0),
_dbSlipTime(0)
{ // constructor
@@ -45,7 +44,6 @@
// Destructor.
pylith::faults::StepSlipFn::~StepSlipFn(void)
{ // destructor
- delete _parameters; _parameters = 0;
_dbFinalSlip = 0; // :TODO: Use shared pointer
_dbSlipTime = 0; // :TODO: Use shared pointer
} // destructor
Modified: short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/faults/StepSlipFn.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -25,8 +25,6 @@
// Include directives ---------------------------------------------------
#include "SlipTimeFn.hh"
-#include "pylith/topology/topologyfwd.hh" // USES Fields<Field<SubMesh> >
-
#include "pylith/utils/array.hh" // HASA double_array
// StepSlipFn -----------------------------------------------------------
@@ -112,10 +110,6 @@
double _slipTimeVertex; ///< Slip time at a vertex.
double_array _slipVertex; ///< Final slip at a vertex.
- /// Parameters for step slip time function, final slip (vector) and
- /// slip time (scalar).
- topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
-
/// Spatial database for final slip
spatialdata::spatialdb::SpatialDB* _dbFinalSlip;
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -19,6 +19,7 @@
#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/utils/EventLogger.hh" // USES EventLogger
@@ -41,8 +42,7 @@
// Constructor
pylith::feassemble::IntegratorElasticity::IntegratorElasticity(void) :
_material(0),
- _bufferFieldTensor(0),
- _bufferFieldOther(0)
+ _outputFields(0)
{ // constructor
} // constructor
@@ -50,9 +50,8 @@
// Destructor
pylith::feassemble::IntegratorElasticity::~IntegratorElasticity(void)
{ // destructor
+ delete _outputFields; _outputFields = 0;
_material = 0; // Don't manage memory for material
- delete _bufferFieldTensor; _bufferFieldTensor = 0;
- delete _bufferFieldOther; _bufferFieldOther = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -319,32 +318,43 @@
// We assume the material stores the total_strain field if
// hasStateVars() is TRUE.
+ if (0 == _outputFields)
+ _outputFields =
+ new topology::Fields<topology::Field<topology::Mesh> >(mesh);
+
if (!_material->hasStateVars() &&
(0 == strcasecmp(name, "total_strain") ||
0 == strcasecmp(name, "stress") )) {
assert(0 != fields);
_allocateTensorField(mesh);
- _calcStrainStressField(_bufferFieldTensor, name, fields);
- _bufferFieldTensor->label(name);
- return *_bufferFieldTensor;
+ topology::Field<topology::Mesh>& buffer =
+ _outputFields->get("buffer (tensor)");
+ _calcStrainStressField(&buffer, name, fields);
+ buffer.label(name);
+ return buffer;
} else if (0 == strcasecmp(name, "stress")) {
assert(0 != fields);
_allocateTensorField(mesh);
- _material->getField(_bufferFieldTensor, "total_strain");
- _calcStressFromStrain(_bufferFieldTensor);
- _bufferFieldTensor->label(name);
- return *_bufferFieldTensor;
+ topology::Field<topology::Mesh>& buffer =
+ _outputFields->get("buffer (tensor)");
+ _material->getField(&buffer, "total_strain");
+ _calcStressFromStrain(&buffer);
+ buffer.label(name);
+ return buffer;
} else {
- if (0 == _bufferFieldOther)
- _bufferFieldOther = new topology::Field<topology::Mesh>(mesh);
- _bufferFieldOther->label("scalar field buffer");
- _material->getField(_bufferFieldOther, name);
- return *_bufferFieldOther;
+ if (!_outputFields->hasField("buffer (other)"))
+ _outputFields->add("buffer (other)", "buffer");
+ topology::Field<topology::Mesh>& buffer =
+ _outputFields->get("buffer (other)");
+ _material->getField(&buffer, name);
+ return buffer;
} // if/else
// Return tensor section to satisfy member function definition. Code
// should never get here.
- return *_bufferFieldTensor;
+ topology::Field<topology::Mesh>& buffer =
+ _outputFields->get("buffer (tensor)");
+ return buffer;
} // cellField
// ----------------------------------------------------------------------
@@ -355,6 +365,7 @@
{ // _allocateTensorField
assert(0 != _quadrature);
assert(0 != _material);
+ assert(0 != _outputFields);
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
assert(!sieveMesh.isNull());
@@ -368,13 +379,14 @@
const int spaceDim = _quadrature->spaceDim();
const int tensorSize = _material->tensorSize();
- if (0 == _bufferFieldTensor) {
- _bufferFieldTensor = new topology::Field<topology::Mesh>(mesh);
- _bufferFieldTensor->label("tensor field buffer");
- assert(0 != _bufferFieldTensor);
- _bufferFieldTensor->newSection(cells, numQuadPts*tensorSize);
- _bufferFieldTensor->allocate();
- _bufferFieldTensor->vectorFieldType(topology::FieldBase::MULTI_TENSOR);
+ if (!_outputFields->hasField("buffer (tensor)")) {
+ _outputFields->add("buffer (tensor)", "buffer");
+ topology::Field<topology::Mesh>& buffer =
+ _outputFields->get("buffer (tensor)");
+ buffer.newSection(cells, numQuadPts*tensorSize);
+ buffer.allocate();
+ buffer.vectorFieldType(topology::FieldBase::MULTI_TENSOR);
+ buffer.addDimensionOkay(true);
} // if
} // _allocateTensorField
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -223,12 +223,9 @@
/// Elastic material associated with integrator
materials::ElasticMaterial* _material;
- /// Buffer for storing cell tensor field.
- topology::Field<topology::Mesh>* _bufferFieldTensor;
+ /// Buffers for output.
+ topology::Fields<topology::Field<topology::Mesh> >* _outputFields;
- /// Buffer for storing cell state-variable field.
- topology::Field<topology::Mesh>* _bufferFieldOther;
-
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -23,6 +23,7 @@
#include "Quadrature2Din3D.hh"
#include "Quadrature3D.hh"
+#include "pylith/topology/Fields.hh" // HOLDSA Fields
#include "pylith/topology/Field.hh" // HOLDSA Field
#include <cstring> // USES memcpy()
@@ -35,10 +36,7 @@
template<typename mesh_type>
pylith::feassemble::Quadrature<mesh_type>::Quadrature(void) :
_engine(0),
- _quadPtsField(0),
- _jacobianField(0),
- _jacobianDetField(0),
- _basisDerivField(0),
+ _geometryFields(0),
_checkConditioning(false)
{ // constructor
} // constructor
@@ -49,10 +47,7 @@
pylith::feassemble::Quadrature<mesh_type>::~Quadrature(void)
{ // destructor
delete _engine; _engine = 0;
- delete _quadPtsField; _quadPtsField = 0;
- delete _jacobianField; _jacobianField = 0;
- delete _jacobianDetField; _jacobianDetField = 0;
- delete _basisDerivField; _basisDerivField = 0;
+ delete _geometryFields; _geometryFields = 0;
} // destructor
// ----------------------------------------------------------------------
@@ -61,10 +56,7 @@
pylith::feassemble::Quadrature<mesh_type>::Quadrature(const Quadrature& q) :
QuadratureRefCell(q),
_engine(0),
- _quadPtsField(0),
- _jacobianField(0),
- _jacobianDetField(0),
- _basisDerivField(0),
+ _geometryFields(0),
_checkConditioning(q._checkConditioning)
{ // copy constructor
if (0 != q._engine)
@@ -152,50 +144,53 @@
logger.setDebug(1);
logger.stagePush(loggingStage);
+ delete _geometryFields;
+ _geometryFields = new topology::Fields<topology::Field<mesh_type> >;
+
// Allocate field and cell buffer for quadrature points
+ _geometryFields->add("quadrature points", "quadrature_points");
+ topology::Field<mesh_type>& quadPtsField =
+ _geometryFields->get("quadrature points");
int fiberDim = _numQuadPts * _spaceDim;
- std::cout << "Quadrature points: " << _numQuadPts << " space dim: " << _spaceDim << " numCells: " << cells->size() << std::endl;
- _quadPtsField = new topology::Field<mesh_type>(mesh);
- _quadPtsField->label("quadPoints");
- assert(0 != _quadPtsField);
- _quadPtsField->newSection(cells, fiberDim);
- _quadPtsField->allocate();
+ quadPtsField.newSection(cells, fiberDim);
+ quadPtsField.allocate();
// Get chart for reuse in other fields
- const ALE::Obj<RealSection>& section = _quadPtsField->section();
+ const ALE::Obj<RealSection>& section = quadPtsField->section();
assert(!section.isNull());
const typename RealSection::chart_type& chart = section->getChart();
// Allocate field and cell buffer for Jacobian at quadrature points
logger.setDebug(2);
std::cout << "Jacobian: cell dim: " << _cellDim << std::endl;
+ _geometryFields->add("jacobian", "jacobian");
+ topology::Field<mesh_type>& jacobianField =
+ _geometryFields->get("jacobian");
fiberDim = (_cellDim > 0) ?
_numQuadPts * _cellDim * _spaceDim :
_numQuadPts * 1 * _spaceDim;
- _jacobianField = new topology::Field<mesh_type>(mesh);
- _jacobianField->label("jacobian");
- assert(0 != _jacobianField);
- _jacobianField->newSection(chart, fiberDim);
- _jacobianField->allocate();
+ jacobianField->newSection(chart, fiberDim);
+ jacobianField->allocate();
logger.setDebug(1);
// Allocate field and cell buffer for determinant of Jacobian at quad pts
std::cout << "Jacobian det:" << std::endl;
+ _geometryFields->add("determinant(jacobian)", "determinant_jacobian");
+ topology::Field<mesh_type>& jacobianDetField =
+ _geometryFields->get("determinant(jacobian)");
fiberDim = _numQuadPts;
- _jacobianDetField = new topology::Field<mesh_type>(mesh);
- _jacobianDetField->label("jacobianDet");
- assert(0 != _jacobianDetField);
- _jacobianDetField->newSection(chart, fiberDim);
- _jacobianDetField->allocate();
+ jacobianDetField.newSection(chart, fiberDim);
+ jacobianDetField.allocate();
// Allocate field for derivatives of basis functions at quad pts
std::cout << "Basis derivatives: num basis: " << _numBasis << std::endl;
+ _geometryFields->add("derivative basis functions",
+ "derivative_basis_functions");
+ topology::Field<mesh_type>& basisDerivField =
+ _geometryFields->get("jacobian");
fiberDim = _numQuadPts * _numBasis * _spaceDim;
- _basisDerivField = new topology::Field<mesh_type>(mesh);
- _basisDerivField->label("basis derivatives");
- assert(0 != _basisDerivField);
- _basisDerivField->newSection(chart, fiberDim);
- _basisDerivField->allocate();
+ basisDerivField.newSection(chart, fiberDim);
+ basisDerivField.allocate();
logger.stagePop();
logger.setDebug(0);
@@ -228,11 +223,11 @@
RestrictVisitor coordsVisitor(*coordinates,
coordinatesCell.size(), &coordinatesCell[0]);
- const ALE::Obj<RealSection>& quadPtsSection = _quadPtsField->section();
- const ALE::Obj<RealSection>& jacobianSection = _jacobianField->section();
+ const ALE::Obj<RealSection>& quadPtsSection = quadPtsField->section();
+ const ALE::Obj<RealSection>& jacobianSection = jacobianField->section();
const ALE::Obj<RealSection>& jacobianDetSection =
- _jacobianDetField->section();
- const ALE::Obj<RealSection>& basisDerivSection = _basisDerivField->section();
+ jacobianDetField->section();
+ const ALE::Obj<RealSection>& basisDerivSection = basisDerivField->section();
const double_array& quadPts = _engine->quadPts();
const double_array& jacobian = _engine->jacobian();
@@ -259,33 +254,33 @@
void
pylith::feassemble::Quadrature<mesh_type>::retrieveGeometry(const typename mesh_type::SieveMesh::point_type& cell)
{ // retrieveGeometry
- typedef typename mesh_type::RealSection RealSection;
-
- assert(0 != _quadPtsField);
- assert(0 != _jacobianField);
- assert(0 != _jacobianDetField);
- assert(0 != _basisDerivField);
+ assert(0 != _geometryFields);
assert(0 != _engine);
+ typedef typename mesh_type::RealSection RealSection;
+
const double_array& quadPts = _engine->quadPts();
const double_array& jacobian = _engine->jacobian();
const double_array& jacobianDet = _engine->jacobianDet();
const double_array& basisDeriv = _engine->basisDeriv();
- const ALE::Obj<RealSection>& quadPtsSection = _quadPtsField->section();
+ const ALE::Obj<RealSection>& quadPtsSection =
+ _geometryFields->get("quadrature points").section();
quadPtsSection->restrictPoint(cell, const_cast<double*>(&quadPts[0]),
quadPts.size());
- const ALE::Obj<RealSection>& jacobianSection = _jacobianField->section();
+ const ALE::Obj<RealSection>& jacobianSection =
+ _geometryFields->get("jacobian").section();
jacobianSection->restrictPoint(cell, const_cast<double*>(&jacobian[0]),
jacobian.size());
const ALE::Obj<RealSection>& jacobianDetSection =
- _jacobianDetField->section();
+ _geometryFields->get("determinant(jacobian)").section();
jacobianDetSection->restrictPoint(cell, const_cast<double*>(&jacobianDet[0]),
jacobianDet.size());
- const ALE::Obj<RealSection>& basisDerivSection = _basisDerivField->section();
+ const ALE::Obj<RealSection>& basisDerivSection =
+ _geometryFields->get("determinant basisfunctions").section();
basisDerivSection->restrictPoint(cell, const_cast<double*>(&basisDeriv[0]),
basisDeriv.size());
} // retrieveGeometry
@@ -298,11 +293,8 @@
{ // clear
delete _engine; _engine = 0;
- // Clear storage for fields
- delete _quadPtsField; _quadPtsField = 0;
- delete _jacobianField; _jacobianField = 0;
- delete _jacobianDetField; _jacobianDetField = 0;
- delete _basisDerivField; _basisDerivField = 0;
+ // Clear storage of precomputed geometry.
+ delete _geometryFields; _geometryFields = 0;
} // clear
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -97,30 +97,12 @@
*/
const double_array& jacobianDet(void) const;
- /** Get precomputed coordinates of quadrature points
- *
- * @returns Array of coordinates of quadrature points in cell
- */
- const pylith::topology::Field<mesh_type>& quadPtsPrecomp(void) const;
-
- /** Get precomputed derivatives of basis fns evaluated at quadrature points.
- *
- * @returns Array of derivatives of basis fns evaluated at
- * quadrature points
- */
- const pylith::topology::Field<mesh_type>& basisDerivPrecomp(void) const;
-
- /** Get precomputed Jacobians evaluated at quadrature points.
- *
- * @returns Array of Jacobian inverses evaluated at quadrature points.
- */
- const pylith::topology::Field<mesh_type>& jacobianPrecomp(void) const;
-
/** Get precomputed determinants of Jacobian evaluated at quadrature points.
*
* @returns Array of determinants of Jacobian evaluated at quadrature pts
*/
- const pylith::topology::Field<mesh_type>& jacobianDetPrecomp(void) const;
+ const topology::Fields<topology::Field<mesh_type> >&
+ geometryFields(void) const;
/** Compute geometric quantities for each cell.
*
@@ -153,16 +135,9 @@
QuadratureEngine* _engine; ///< Quadrature geometry engine.
- /** Fields and visitors for precomputing geometry information for
- * cells associated with this quadrature.
- */
- topology::Field<mesh_type>* _quadPtsField; ///< Coordinates of quad pts.
- topology::Field<mesh_type>* _jacobianField; ///< Jacobian at quad pts.
- topology::Field<mesh_type>* _jacobianDetField; ///< |J| at quad pts.
+ /// Fields for precomputing geometry information.
+ topology::Fields<topology::Field<mesh_type> >* _geometryFields;
- /// Derivatives of basis fns at quad pts.
- topology::Field<mesh_type>* _basisDerivField;
-
bool _checkConditioning; ///< True if checking for ill-conditioning.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -69,42 +69,15 @@
return _engine->jacobianDet();
}
-// Get precomputed coordinates of quadrature points
+// Get precomputed geometry fields.
template<typename mesh_type>
inline
-const pylith::topology::Field<mesh_type>&
-pylith::feassemble::Quadrature<mesh_type>::quadPtsPrecomp(void) const {
- assert(0 != _quadPtsField);
- return *_quadPtsField;
+const pylith::topology::Fields<pylith::topology::Field<mesh_type> >&
+pylith::feassemble::Quadrature<mesh_type>::geometryFields(void) const {
+ assert(0 != _geometryFields);
+ return *_geometryFields;
}
-// Get precomputed derivatives of basis fns evaluated at quadrature points.
-template<typename mesh_type>
-inline
-const pylith::topology::Field<mesh_type>&
-pylith::feassemble::Quadrature<mesh_type>::basisDerivPrecomp(void) const {
- assert(0 != _basisDerivField);
- return *_basisDerivField;
-}
-
-// Get precomputed Jacobians evaluated at quadrature points.
-template<typename mesh_type>
-inline
-const pylith::topology::Field<mesh_type>&
-pylith::feassemble::Quadrature<mesh_type>::jacobianPrecomp(void) const {
- assert(0 != _jacobianField);
- return *_jacobianField;
-}
-
-// Get precomputed determinants of Jacobian evaluated at quadrature points.
-template<typename mesh_type>
-inline
-const pylith::topology::Field<mesh_type>&
-pylith::feassemble::Quadrature<mesh_type>::jacobianDetPrecomp(void) const {
- assert(0 != _jacobianDetField);
- return *_jacobianDetField;
-}
-
// Precompute geometric quantities for each cell.
template<typename mesh_type>
inline
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -43,8 +43,7 @@
Material(dimension, tensorSize, metadata),
_dbInitialStress(0),
_dbInitialStrain(0),
- _initialStress(0),
- _initialStrain(0),
+ _initialFields(0),
_numQuadPts(0),
_numElasticConsts(numElasticConsts)
{ // constructor
@@ -54,8 +53,7 @@
// Destructor.
pylith::materials::ElasticMaterial::~ElasticMaterial(void)
{ // destructor
- delete _initialStress; _initialStress = 0;
- delete _initialStrain; _initialStrain = 0;
+ delete _initialFields; _initialFields = 0;
// Python db object owns databases, so just set point to null
// :KLUDGE: Should use shared pointer
@@ -76,6 +74,11 @@
assert(0 != quadrature);
_numQuadPts = quadrature->numQuadPts();
+ if (0 != _dbInitialStress || 0 != _dbInitialStrain) {
+ delete _initialFields;
+ _initialFields =
+ new topology::Fields<topology::Field<topology::Mesh> >(mesh);
+ } // if
_initializeInitialStress(mesh, quadrature);
_initializeInitialStrain(mesh, quadrature);
_allocateCellArrays();
@@ -101,22 +104,24 @@
_stateVarsCell.size());
} // if
- if (0 == _initialStress)
- _initialStressCell = 0.0;
- else {
- const ALE::Obj<RealSection>& stressSection = _initialStress->section();
- assert(!stressSection.isNull());
- stressSection->restrictPoint(cell, &_initialStressCell[0],
- _initialStressCell.size());
- } // if/else
- if (0 == _initialStrain)
- _initialStrainCell = 0.0;
- else {
- const ALE::Obj<RealSection>& strainSection = _initialStrain->section();
- assert(!strainSection.isNull());
- strainSection->restrictPoint(cell, &_initialStrainCell[0],
- _initialStrainCell.size());
- } // if/else
+ _initialStressCell = 0.0;
+ _initialStrainCell = 0.0;
+ if (0 != _initialFields) {
+ if (_initialFields->hasField("initial stress")) {
+ const ALE::Obj<RealSection>& stressSection =
+ _initialFields->get("initial stress").section();
+ assert(!stressSection.isNull());
+ stressSection->restrictPoint(cell, &_initialStressCell[0],
+ _initialStressCell.size());
+ } // if
+ if (_initialFields->hasField("initial strain")) {
+ const ALE::Obj<RealSection>& strainSection =
+ _initialFields->get("initial strain").section();
+ assert(!strainSection.isNull());
+ strainSection->restrictPoint(cell, &_initialStrainCell[0],
+ _initialStrainCell.size());
+ } // if
+ } // if
} // retrievePropsAndVars
// ----------------------------------------------------------------------
@@ -281,11 +286,18 @@
const pylith::topology::Field<pylith::topology::Mesh>&
pylith::materials::ElasticMaterial::initialStressField(void) const
{ // initialStressField
- if (0 == _initialStress)
+ if (0 == _initialFields)
+ throw std::runtime_error("Request for initial stress field, but no "
+ "initial fields exist.");
+ assert(0 != _initialFields);
+ if (!_initialFields->hasField("initial stress"))
throw std::runtime_error("Request for initial stress field, but field "
"does not exist.");
+
+ const topology::Field<topology::Mesh>& initialStress =
+ _initialFields->get("initial stress");
- return *_initialStress;
+ return initialStress;
} // initialStressField
// ----------------------------------------------------------------------
@@ -293,11 +305,18 @@
const pylith::topology::Field<pylith::topology::Mesh>&
pylith::materials::ElasticMaterial::initialStrainField(void) const
{ // initialStrainField
- if (0 == _initialStrain)
+ if (0 == _initialFields)
+ throw std::runtime_error("Request for initial strain field, but no "
+ "initial fields exist.");
+ assert(0 != _initialFields);
+ if (!_initialFields->hasField("initial strain"))
throw std::runtime_error("Request for initial strain field, but field "
"does not exist.");
+
+ const topology::Field<topology::Mesh>& initialStrain =
+ _initialFields->get("initial strain");
- return *_initialStrain;
+ return initialStrain;
} // initialStrainField
// ----------------------------------------------------------------------
@@ -327,10 +346,14 @@
const topology::Mesh& mesh,
feassemble::Quadrature<topology::Mesh>* quadrature)
{ // _initializeInitialStress
- delete _initialStress; _initialStress = 0;
if (0 == _dbInitialStress)
return;
+ assert(0 != _initialFields);
+ _initialFields->add("initial stress", "initial_stress");
+ topology::Field<topology::Mesh>& initialStress =
+ _initialFields->get("initial stress");
+
assert(0 != _dbInitialStress);
assert(0 != quadrature);
@@ -365,17 +388,12 @@
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
// Create field to hold initial stress state.
- delete _initialStress; _initialStress = 0;
- _initialStress = new topology::Field<topology::Mesh>(mesh);
- _initialStress->label("initial stress");
- assert(0 != _initialStress);
const int fiberDim = numQuadPts * tensorSize;
assert(fiberDim > 0);
- _initialStress->newSection(cells, fiberDim);
- _initialStress->allocate();
- _initialStress->zero();
- const ALE::Obj<RealSection>& initialStressSection =
- _initialStress->section();
+ initialStress.newSection(cells, fiberDim);
+ initialStress.allocate();
+ initialStress.zero();
+ const ALE::Obj<RealSection>& initialStressSection = initialStress.section();
// Setup databases for querying
_dbInitialStress->open();
@@ -415,9 +433,6 @@
const double lengthScale = _normalizer->pressureScale();
const double pressureScale = _normalizer->pressureScale();
- const ALE::Obj<RealSection>& stressSection = _initialStress->section();
- assert(!stressSection.isNull());
-
for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
@@ -457,7 +472,7 @@
_normalizer->nondimensionalize(&stressCell[0], stressCell.size(),
pressureScale);
- stressSection->updatePoint(*c_iter, &stressCell[0]);
+ initialStressSection->updatePoint(*c_iter, &stressCell[0]);
} // for
// Close databases
@@ -471,10 +486,14 @@
const topology::Mesh& mesh,
feassemble::Quadrature<topology::Mesh>* quadrature)
{ // _initializeInitialStrain
- delete _initialStrain; _initialStrain = 0;
if (0 == _dbInitialStrain)
return;
+ assert(0 != _initialFields);
+ _initialFields->add("initial strain", "initial_strain");
+ topology::Field<topology::Mesh>& initialStrain =
+ _initialFields->get("initial strain");
+
assert(0 != _dbInitialStrain);
assert(0 != quadrature);
@@ -509,17 +528,12 @@
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
// Create field to hold initial strain state.
- delete _initialStrain; _initialStrain = 0;
- _initialStrain = new topology::Field<topology::Mesh>(mesh);
- _initialStrain->label("initial stress");
- assert(0 != _initialStrain);
const int fiberDim = numQuadPts * tensorSize;
assert(fiberDim > 0);
- _initialStrain->newSection(cells, fiberDim);
- _initialStrain->allocate();
- _initialStrain->zero();
- const ALE::Obj<RealSection>& initialStrainSection =
- _initialStrain->section();
+ initialStrain.newSection(cells, fiberDim);
+ initialStrain.allocate();
+ initialStrain.zero();
+ const ALE::Obj<RealSection>& initialStrainSection = initialStrain.section();
// Setup databases for querying
_dbInitialStrain->open();
@@ -559,9 +573,6 @@
const double lengthScale = _normalizer->pressureScale();
const double pressureScale = _normalizer->pressureScale();
- const ALE::Obj<RealSection>& strainSection = _initialStrain->section();
- assert(!strainSection.isNull());
-
for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter) {
@@ -601,7 +612,7 @@
_normalizer->nondimensionalize(&strainCell[0], strainCell.size(),
pressureScale);
- strainSection->updatePoint(*c_iter, &strainCell[0]);
+ initialStrainSection->updatePoint(*c_iter, &strainCell[0]);
} // for
// Close databases
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -347,12 +347,9 @@
/// Database for initial strain tensor;
spatialdata::spatialdb::SpatialDB* _dbInitialStrain;
- /// Initial stress field.
- topology::Field<topology::Mesh>* _initialStress;
+ /// Initial stress/strain fields.
+ topology::Fields<topology::Field<topology::Mesh> >* _initialFields;
- /// Initial strain field.
- topology::Field<topology::Mesh>* _initialStrain;
-
/** Properties at quadrature points for current cell.
*
* size = numQuadPts * numPropsQuadPt
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -267,7 +267,7 @@
// ----------------------------------------------------------------------
// Get properties
-const pylith::materials::Material::field_type&
+const pylith::topology::Field<pylith::topology::Mesh>&
pylith::materials::Material::getProperties() const
{ // getProperties
return *_properties;
@@ -275,7 +275,7 @@
// ----------------------------------------------------------------------
// Get state variables
-const pylith::materials::Material::field_type&
+const pylith::topology::Field<pylith::topology::Mesh>&
pylith::materials::Material::getStateVars() const
{ // getStateVars
return *_stateVars;
@@ -284,7 +284,7 @@
// ----------------------------------------------------------------------
// Get physical property or state variable field.
void
-pylith::materials::Material::getField(field_type *field, const char* name) const
+pylith::materials::Material::getField(topology::Field<topology::Mesh> *field, const char* name) const
{ // getField
assert(0 != field);
assert(0 != _properties);
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -38,8 +38,6 @@
class pylith::materials::Material
{ // class Material
friend class TestMaterial; // unit testing
-public:
- typedef pylith::topology::Field<pylith::topology::Mesh> field_type;
// PUBLIC METHODS /////////////////////////////////////////////////////
public :
@@ -155,19 +153,19 @@
* @param field Field over material cells.
* @param name Name of field to retrieve.
*/
- void getField(field_type *field, const char* name) const;
+ void getField(topology::Field<topology::Mesh> *field, const char* name) const;
/** Get the properties.
*
* @returns the properties
*/
- const field_type& getProperties() const;
+ const topology::Field<topology::Mesh>& getProperties() const;
/** Get the state variables.
*
* @returns the state variables
*/
- const field_type& getStateVars() const;
+ const topology::Field<topology::Mesh>& getStateVars() const;
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -234,10 +232,10 @@
double _dt; ///< Current time step
/// Field containing physical properties of material.
- field_type *_properties;
+ topology::Field<topology::Mesh> *_properties;
/// Field containing the state variables for the material.
- field_type *_stateVars;
+ topology::Field<topology::Mesh> *_stateVars;
spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer
Modified: short/3D/PyLith/trunk/libsrc/topology/Fields.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Fields.hh 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/topology/Fields.hh 2009-05-22 03:10:47 UTC (rev 15030)
@@ -48,6 +48,13 @@
/// Deallocate PETSc and local data structures.
void deallocate(void);
+ /** Check if fields contains a given field.
+ *
+ * @param name Name of field.
+ * @return True if fields contains field, false otherwise.
+ */
+ bool hasField(const char* name) const;
+
/** Add field.
*
* @param name Name of field.
Modified: short/3D/PyLith/trunk/libsrc/topology/Fields.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Fields.icc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/libsrc/topology/Fields.icc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -47,14 +47,23 @@
} // deallocate
// ----------------------------------------------------------------------
+// Check if fields contains a given field.
+template<typename field_type>
+bool
+pylith::topology::Fields<field_type>::hasField(const char* name) const
+{ // hasField
+ typename map_type::const_iterator iter = _fields.find(name);
+ return iter != _fields.end();
+} // hasField
+
+// ----------------------------------------------------------------------
// Add field.
template<typename field_type>
void
pylith::topology::Fields<field_type>::add(const char* name,
const char* label)
{ // add
- typename map_type::iterator iter = _fields.find(name);
- if (iter != _fields.end()) {
+ if (hasField(name)) {
std::ostringstream msg;
msg << "Could not add field '" << name
<< "' to fields manager, because it already exists.";
@@ -75,8 +84,7 @@
const pylith::topology::FieldBase::DomainEnum domain,
const int fiberDim)
{ // add
- typename map_type::iterator iter = _fields.find(name);
- if (iter != _fields.end()) {
+ if (hasField(name)) {
std::ostringstream msg;
msg << "Could not add field '" << name
<< "' to fields manager, because it already exists.";
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/Quadrature.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/Quadrature.i 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/Quadrature.i 2009-05-22 03:10:47 UTC (rev 15030)
@@ -59,27 +59,8 @@
*
* @returns Array of coordinates of quadrature points in cell
*/
- const pylith::topology::Field<mesh_type>& quadPtsPrecomp(void) const;
+ const pylith::topology::Fields<pylith::topology::Field<mesh_type> >& geometryFields(void) const;
- /** Get precomputed derivatives of basis fns evaluated at quadrature points.
- *
- * @returns Array of derivatives of basis fns evaluated at
- * quadrature points
- */
- const pylith::topology::Field<mesh_type>& basisDerivPrecomp(void) const;
-
- /** Get precomputed Jacobians evaluated at quadrature points.
- *
- * @returns Array of Jacobian inverses evaluated at quadrature points.
- */
- const pylith::topology::Field<mesh_type>& jacobianPrecomp(void) const;
-
- /** Get precomputed determinants of Jacobian evaluated at quadrature points.
- *
- * @returns Array of determinants of Jacobian evaluated at quadrature pts
- */
- const pylith::topology::Field<mesh_type>& jacobianDetPrecomp(void) const;
-
}; // Quadrature
} // feassemble
Modified: short/3D/PyLith/trunk/modulesrc/materials/Material.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/materials/Material.i 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/modulesrc/materials/Material.i 2009-05-22 03:10:47 UTC (rev 15030)
@@ -21,8 +21,6 @@
class Material
{ // class Material
- public:
- typedef pylith::topology::Field<pylith::topology::Mesh> field_type;
// PUBLIC METHODS /////////////////////////////////////////////////
public :
@@ -124,19 +122,20 @@
* @param field Field over material cells.
* @param name Name of field to retrieve.
*/
- void getField(field_type *field, const char* name) const;
+ void getField(pylith::topology::Field<pylith::topology::Mesh>* field,
+ const char* name) const;
/** Get the properties.
*
* @returns the properties
*/
- const field_type& getProperties() const;
+ const pylith::topology::Field<pylith::topology::Mesh>& getProperties() const;
/** Get the state variables.
*
* @returns the state variables
*/
- const field_type& getStateVars() const;
+ const pylith::topology::Field<pylith::topology::Mesh>& getStateVars() const;
// PROTECTED METHODS //////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/modulesrc/topology/Fields.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/Fields.i 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/modulesrc/topology/Fields.i 2009-05-22 03:10:47 UTC (rev 15030)
@@ -38,6 +38,13 @@
/// Deallocate PETSc and local data structures.
void deallocate(void);
+ /** Check if fields contains a given field.
+ *
+ * @param name Name of field.
+ * @return True if fields contains field, false otherwise.
+ */
+ bool hasField(const char* name) const;
+
/** Add field.
*
* @param name Name of field.
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -120,7 +120,9 @@
const int fiberDim = numQuadPts * spaceDim;
double_array dampersCell(fiberDim);
int index = 0;
- const ALE::Obj<SubRealSection>& dampersSection = bc._dampingConsts->section();
+ CPPUNIT_ASSERT(0 != bc._parameters);
+ const ALE::Obj<SubRealSection>& dampersSection =
+ bc._parameters->get("damping constants").section();
const double tolerance = 1.0e-06;
for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -142,7 +142,9 @@
const int fiberDim = numQuadPts * spaceDim;
double_array tractionsCell(fiberDim);
int index = 0;
- const ALE::Obj<SubRealSection>& tractionSection = bc._tractions->section();
+ CPPUNIT_ASSERT(0 != bc._parameters);
+ const ALE::Obj<SubRealSection>& tractionSection =
+ bc._parameters->get("traction").section();
for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cells->end();
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -225,8 +225,18 @@
const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
CPPUNIT_ASSERT(!cells.isNull());
quadrature.initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
quadrature.computeGeometry(mesh, cells);
-
+#else
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<topology::Mesh::RealSection>& coordinates =
+ sieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+#endif
+
size_t size = 0;
// Check values from computeGeometry()
@@ -235,7 +245,13 @@
for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
+#if defined(PRECOMPUTE_GEOMETRY)
quadrature.retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ quadrature.computeGeometry(coordinatesCell, *c_iter);
+#endif
const double_array& quadPts = quadrature.quadPts();
size = numQuadPts * spaceDim;
@@ -265,10 +281,7 @@
// Check clear()
quadrature.clear();
- CPPUNIT_ASSERT(0 == quadrature._quadPtsField);
- CPPUNIT_ASSERT(0 == quadrature._jacobianField);
- CPPUNIT_ASSERT(0 == quadrature._jacobianDetField);
- CPPUNIT_ASSERT(0 == quadrature._basisDerivField);
+ CPPUNIT_ASSERT(0 == quadrature._geometryFields);
CPPUNIT_ASSERT(0 == quadrature._engine);
// Make sure caling clear without data doesn't generate errors
@@ -373,10 +386,7 @@
quadrature.clear();
- CPPUNIT_ASSERT(0 == quadrature._quadPtsField);
- CPPUNIT_ASSERT(0 == quadrature._jacobianField);
- CPPUNIT_ASSERT(0 == quadrature._jacobianDetField);
- CPPUNIT_ASSERT(0 == quadrature._basisDerivField);
+ CPPUNIT_ASSERT(0 == quadrature._geometryFields);
CPPUNIT_ASSERT(0 == quadrature._engine);
} // testComputeGeometryCell
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2009-05-22 03:10:47 UTC (rev 15030)
@@ -98,7 +98,7 @@
const int numQuadPts = data.numLocs;
// Test initialStress field
- CPPUNIT_ASSERT(0 != material._initialStress);
+ CPPUNIT_ASSERT(0 != material._initialFields);
const ALE::Obj<RealSection>& stressSection =
material.initialStressField().section();
CPPUNIT_ASSERT(!stressSection.isNull());
@@ -113,7 +113,7 @@
tolerance);
// Test initialStrain field
- CPPUNIT_ASSERT(0 != material._initialStrain);
+ CPPUNIT_ASSERT(0 != material._initialFields);
const ALE::Obj<RealSection>& strainSection =
material.initialStrainField().section();
CPPUNIT_ASSERT(!strainSection.isNull());
Modified: short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py 2009-05-21 23:38:23 UTC (rev 15029)
+++ short/3D/PyLith/trunk/unittests/pytests/faults/TestFaultCohesiveKin.py 2009-05-22 03:10:47 UTC (rev 15030)
@@ -326,7 +326,7 @@
fault.inventory.faultLabel = "fault"
fault.inventory.upDir = [0, 0, 1]
fault.inventory.normalDir = [1, 0, 0]
- fault.inventory.quadrature = quadrature
+ fault.inventory.faultQuadrature = quadrature
fault.inventory.matDB = dbMat
fault._configure()
eqsrc = fault.eqsrcs.components()[0]
More information about the CIG-COMMITS
mailing list