[cig-commits] r14051 - in short/3D/PyLith/branches/pylith-swig: libsrc libsrc/bc libsrc/feassemble libsrc/topology unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri Feb 13 16:55:31 PST 2009
Author: brad
Date: 2009-02-13 16:55:30 -0800 (Fri, 13 Feb 2009)
New Revision: 14051
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.cc
short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.hh
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.hh
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.cc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.hh
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.icc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.hh
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh
short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestIntegrator.hh
Log:
Worked on updating integrator implementation.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-02-14 00:55:30 UTC (rev 14051)
@@ -25,6 +25,7 @@
bc/BoundaryCondition.cc \
bc/DirichletBC.cc \
bc/DirichletBoundary.cc \
+ bc/Neumann.cc \
feassemble/CellGeometry.cc \
feassemble/Constraint.cc \
feassemble/GeometryPoint1D.cc \
@@ -64,8 +65,6 @@
topology/MeshOps.cc \
utils/EventLogger.cc
-# feassemble/Integrator.cc \
-# bc/Neumann.cc \
# faults/BruneSlipFn.cc \
# faults/ConstRateSlipFn.cc \
# faults/CohesiveTopology.cc \
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.cc 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.cc 2009-02-14 00:55:30 UTC (rev 14051)
@@ -12,11 +12,12 @@
#include <portinfo>
+#include "pylith/topology/SubMesh.hh" // HOLDSA SubMesh
+#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+
#include "Neumann.hh" // implementation of object methods
-#include "pylith/topology/SubMesh.hh" // HOLDSA SubMesh
-#include "pylith/topology/FieldSubMesh.hh" // HOLDSA FieldSubMesh
-#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
+#include "pylith/topology/Field.hh" // HOLDSA Field
#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
@@ -29,6 +30,11 @@
#include <sstream> // USES std::ostringstream
// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::RealSection SubRealSection;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::Neumann::Neumann(void) :
_boundaryMesh(0),
@@ -60,6 +66,8 @@
_boundaryMesh = new topology::SubMesh(mesh, _label.c_str());
assert(0 != _boundaryMesh);
+ double_array up(upDir, 3);
+
// check compatibility of quadrature and boundary mesh
if (_quadrature->cellDim() != _boundaryMesh->dimension()) {
std::ostringstream msg;
@@ -76,10 +84,10 @@
const ALE::Obj<SieveSubMesh>& submesh = _boundaryMesh->sieveMesh();
assert(!submesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- _boundaryMesh->heightStratum(1);
+ submesh->heightStratum(1);
assert(!cells.isNull());
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
- const int boundaryDepth = _boundaryMesh->depth()-1; // depth of bndry cells
+ const int boundaryDepth = submesh->depth()-1; // depth of bndry cells
// Make sure surface cells are compatible with quadrature.
for (SieveSubMesh::label_sequence::iterator c_iter=cells->begin();
@@ -99,25 +107,25 @@
} // for
// Create section for traction vector in global coordinates
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
const int cellDim = _quadrature->cellDim();
const int numBasis = _quadrature->numBasis();
const int numQuadPts = _quadrature->numQuadPts();
- const int spaceDim = cs->spaceDim();
+ const int spaceDim = cellGeometry.spaceDim();
const int fiberDim = spaceDim * numQuadPts;
- _tractions = new FieldSubMesh(submesh);
+ _tractions = new topology::Field<topology::SubMesh>(*_boundaryMesh);
assert(0 != _tractions);
_tractions->newSection(cells, fiberDim);
_tractions->allocate();
// Containers for orientation information
const int orientationSize = spaceDim * spaceDim;
- const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
const int jacobianSize = (cellDim > 0) ? spaceDim * cellDim : 1;
double_array jacobian(jacobianSize);
double jacobianDet = 0;
double_array orientation(orientationSize);
- double_array cellVertices(numBasis*spaceDim);
+ double_array cellVertices(numCorners*spaceDim);
// Set names based on dimension of problem.
// 1-D problem = {'normal-traction'}
@@ -159,29 +167,38 @@
double_array cellTractionsGlobal(fiberDim);
// Get sections.
- const ALE::Obj<MeshRealSection>& coordinates =
+ const ALE::Obj<RealSection>& coordinates =
submesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- const ALE::Obj<SubMeshRealSection>& tractSection = _tractions->section();
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ numCorners*spaceDim);
+
+ const ALE::Obj<SubRealSection>& tractSection = _tractions->section();
assert(!tractSection.isNull());
+ const ALE::Obj<topology::Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
const double pressureScale = _normalizer->pressureScale();
+ // Compute quadrature information
+ _quadrature->computeGeometry(*_boundaryMesh, cells);
+
// Loop over cells in boundary mesh, compute orientations, and then
// compute corresponding traction vector in global coordinates
// (store values in _tractionGlobal).
for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
c_iter != cellsEnd;
++c_iter) {
- _quadrature->computeGeometry(submesh, coordinates, *c_iter);
+ _quadrature->retrieveGeometry(*c_iter);
const double_array& quadPtsNondim = _quadrature->quadPts();
quadPtsGlobal = quadPtsNondim;
_normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
lengthScale);
- submesh->restrictClosure(coordinates, *c_iter,
- &cellVertices[0], cellVertices.size());
+
+ sieveMesh->restrictClosure(*c_iter, coordsVisitor);
cellTractionsGlobal = 0.0;
for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts;
@@ -205,8 +222,10 @@
// Compute Jacobian and determinant at quadrature point, then get
// orientation.
memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(double));
+ memcpy(&cellVertices[0], coordsVisitor.getValues(),
+ cellVertices.size()*sizeof(double));
cellGeometry.jacobian(&jacobian, &jacobianDet, cellVertices, quadPtRef);
- cellGeometry.orientation(&orientation, jacobian, jacobianDet, upDir);
+ cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
orientation /= jacobianDet;
// Rotate traction vector from local coordinate system to global
@@ -228,14 +247,13 @@
// ----------------------------------------------------------------------
// Integrate contributions to residual term (r) for operator.
void
-pylith::bc::Neumann::integrateResidual(const topology::Field& residual,
+pylith::bc::Neumann::integrateResidual(const topology::Field<topology::Mesh>& residual,
const double t,
topology::SolutionFields* const fields)
{ // integrateResidual
assert(0 != _quadrature);
assert(0 != _boundaryMesh);
assert(0 != _tractions);
- assert(!residual.isNull());
assert(0 != fields);
PetscErrorCode err = 0;
@@ -249,11 +267,9 @@
const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<MeshRealSection>& coordinates =
- submesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- const ALE::Obj<SubMeshRealSection>& tractSection = _tractions->section();
+ const ALE::Obj<SubRealSection>& tractSection = _tractions->section();
assert(!tractSection.isNull());
+ const ALE::Obj<RealSection>& residualSection = residual.section();
// Get cell geometry information that doesn't depend on cell
const int numQuadPts = _quadrature->numQuadPts();
@@ -269,18 +285,18 @@
double_array tractionsCell(numQuadPts*spaceDim);
// Loop over faces and integrate contribution from each face
- for (SubMesh::label_sequence::iterator c_iter=cellsBegin;
+ for (SieveSubMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
- // Compute geometry information for current cell
- _quadrature->computeGeometry(submesh, coordinates, *c_iter);
+ // Get geometry information for current cell
+ _quadrature->retrieveGeometry(*c_iter);
// Reset element vector to zero
_resetCellVector();
// Restrict tractions to cell
- submesh->restrictClosure(tractSection, *c_iter,
- &tractionsCell[0], tractionsCell.size());
+ tractSection->restrictPoint(*c_iter,
+ &tractionsCell[0], tractionsCell.size());
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -300,7 +316,8 @@
} // for
} // for
// Assemble cell contribution into field
- submesh->updateAdd(residual, *c_iter, _cellVector);
+ // :TODO: Use an UpdateAddVisitor?
+ submesh->updateAdd(residualSection, *c_iter, &_cellVector[0]);
PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*(1+2*spaceDim))));
} // for
@@ -325,7 +342,7 @@
// ----------------------------------------------------------------------
// Get boundary mesh.
-const topology::SubMesh&
+const pylith::topology::SubMesh&
pylith::bc::Neumann::boundaryMesh(void) const
{ // dataMesh
assert(0 != _boundaryMesh);
@@ -335,7 +352,7 @@
// ----------------------------------------------------------------------
// Get cell field for tractions.
-const FieldSubMesh&
+const pylith::topology::Field<pylith::topology::SubMesh>&
pylith::bc::Neumann::cellField(const char* name,
topology::SolutionFields* const fields)
{ // cellField
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.hh 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/bc/Neumann.hh 2009-02-14 00:55:30 UTC (rev 14051)
@@ -27,7 +27,7 @@
// Neumann --------------------------------------------------------------
class pylith::bc::Neumann : public BoundaryCondition,
- public feassemble::Integrator
+ public feassemble::Integrator<feassemble::Quadrature<topology::SubMesh> >
{ // class Neumann
friend class TestNeumann; // unit testing
@@ -55,7 +55,7 @@
* @param t Current time.
* @param fields Solution fields.
*/
- void integrateResidual(const topology::Field& residual,
+ void integrateResidual(const topology::Field<topology::Mesh>& residual,
const double t,
topology::SolutionFields* const fields);
@@ -91,7 +91,7 @@
*
* @returns Traction vector at integration points.
*/
- const topology::Field&
+ const topology::Field<topology::SubMesh>&
cellField(const char* name,
topology::SolutionFields* const fields);
@@ -111,7 +111,7 @@
topology::SubMesh* _boundaryMesh;
/// Traction vector in global coordinates at integration points.
- topology::FieldSubMesh* _tractions;
+ topology::Field<topology::SubMesh>* _tractions;
}; // class Neumann
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.hh 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.hh 2009-02-14 00:55:30 UTC (rev 14051)
@@ -61,7 +61,6 @@
namespace pylith {
namespace feassemble {
class ElasticityExplicit;
- class TestElasticityExplicit;
} // feassemble
} // pylith
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.cc 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.cc 2009-02-14 00:55:30 UTC (rev 14051)
@@ -12,8 +12,6 @@
#include <portinfo>
-#include "Integrator.hh" // implementation of class methods
-
#include "Quadrature.hh" // USES Quadrature
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
@@ -24,7 +22,8 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Integrator::Integrator(void) :
+template<typename quadrature_type>
+pylith::feassemble::Integrator<quadrature_type>::Integrator(void) :
_dt(-1.0),
_quadrature(0),
_normalizer(new spatialdata::units::Nondimensional),
@@ -36,7 +35,8 @@
// ----------------------------------------------------------------------
// Destructor
-pylith::feassemble::Integrator::~Integrator(void)
+template<typename quadrature_type>
+pylith::feassemble::Integrator<quadrature_type>::~Integrator(void)
{ // destructor
delete _quadrature; _quadrature = 0;
delete _normalizer; _normalizer = 0;
@@ -44,8 +44,9 @@
// ----------------------------------------------------------------------
// Set quadrature for integrating finite-element quantities.
+template<typename quadrature_type>
void
-pylith::feassemble::Integrator::quadrature(const Quadrature* q)
+pylith::feassemble::Integrator<quadrature_type>::quadrature(const quadrature_type* q)
{ // quadrature
delete _quadrature;
_quadrature = (0 != q) ? q->clone() : 0;
@@ -57,8 +58,9 @@
// ----------------------------------------------------------------------
// Set manager of scales used to nondimensionalize problem.
+template<typename quadrature_type>
void
-pylith::feassemble::Integrator::normalizer(const spatialdata::units::Nondimensional& dim)
+pylith::feassemble::Integrator<quadrature_type>::normalizer(const spatialdata::units::Nondimensional& dim)
{ // normalizer
if (0 == _normalizer)
_normalizer = new spatialdata::units::Nondimensional(dim);
@@ -68,16 +70,18 @@
// ----------------------------------------------------------------------
// Set gravity field.
+template<typename quadrature_type>
void
-pylith::feassemble::Integrator::gravityField(spatialdata::spatialdb::GravityField* const gravityField)
+pylith::feassemble::Integrator<quadrature_type>::gravityField(spatialdata::spatialdb::GravityField* const gravityField)
{ // gravityField
_gravityField = gravityField;
} // gravityField
// ----------------------------------------------------------------------
// Get stable time step for advancing from time t to time t+dt.
+template<typename quadrature_type>
double
-pylith::feassemble::Integrator::stableTimeStep(void) const
+pylith::feassemble::Integrator<quadrature_type>::stableTimeStep(void) const
{ // stableTimeStep
// Assume any time step will work.
return pylith::PYLITH_MAXDOUBLE;
@@ -85,8 +89,9 @@
// ----------------------------------------------------------------------
// Initialize vector containing result of integration action for cell.
+template<typename quadrature_type>
void
-pylith::feassemble::Integrator::_initCellVector(void)
+pylith::feassemble::Integrator<quadrature_type>::_initCellVector(void)
{ // _initCellVector
assert(0 != _quadrature);
const int size = _quadrature->spaceDim() * _quadrature->numBasis();
@@ -96,16 +101,18 @@
// ----------------------------------------------------------------------
// Zero out vector containing result of integration actions for cell.
+template<typename quadrature_type>
void
-pylith::feassemble::Integrator::_resetCellVector(void)
+pylith::feassemble::Integrator<quadrature_type>::_resetCellVector(void)
{ // _resetCellVector
_cellVector = 0.0;
} // _resetCellVector
// ----------------------------------------------------------------------
// Initialize matrix containing result of integration for cell.
+template<typename quadrature_type>
void
-pylith::feassemble::Integrator::_initCellMatrix(void)
+pylith::feassemble::Integrator<quadrature_type>::_initCellMatrix(void)
{ // _initCellMatrix
assert(0 != _quadrature);
const int size =
@@ -117,8 +124,9 @@
// ----------------------------------------------------------------------
// Zero out matrix containing result of integration for cell.
+template<typename quadrature_type>
void
-pylith::feassemble::Integrator::_resetCellMatrix(void)
+pylith::feassemble::Integrator<quadrature_type>::_resetCellMatrix(void)
{ // _resetCellMatrix
_cellMatrix = 0.0;
} // _resetCellMatrix
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.hh 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.hh 2009-02-14 00:55:30 UTC (rev 14051)
@@ -35,6 +35,7 @@
#include "pylith/utils/array.hh" // HASA double_array
// Integrator -----------------------------------------------------------
+template<typename quadrature_type>
class pylith::feassemble::Integrator
{ // Integrator
friend class TestIntegrator; // unit testing
@@ -54,7 +55,7 @@
*
* @param q Quadrature for integrating.
*/
- void quadrature(const Quadrature* q);
+ void quadrature(const quadrature_type* q);
/** Set manager of scales used to nondimensionalize problem.
*
@@ -62,7 +63,7 @@
*/
void normalizer(const spatialdata::units::Nondimensional& dim);
- /** Set gravity field. Gravity Field should already be initialized.
+ /** Set gravity field.
*
* @param g Gravity field.
*/
@@ -179,21 +180,12 @@
/// Zero out matrix containing result of integration for cell.
void _resetCellMatrix(void);
-// PRIVATE METHODS //////////////////////////////////////////////////////
-private :
-
- // Not implemented.
- Integrator(const Integrator& i);
-
- /// Not implemented
- const Integrator& operator=(const Integrator&);
-
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
double _dt; ///< Time step for t -> t+dt
- Quadrature* _quadrature; ///< Quadrature for integrating finite-element
+ quadrature_type* _quadrature; ///< Quadrature for integrating finite-element
spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer.
spatialdata::spatialdb::GravityField* _gravityField; ///< Gravity field.
@@ -212,10 +204,18 @@
/// solution or an incremental field solution
bool _useSolnIncr;
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ Integrator(const Integrator& i); ///< Not implemented
+ const Integrator& operator=(const Integrator&); ///< Not implemented
+
}; // Integrator
#include "Integrator.icc" // inline methods
+#include "Integrator.cc" // template methods
#endif // pylith_feassemble_integrator_hh
+
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.icc 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.icc 2009-02-14 00:55:30 UTC (rev 14051)
@@ -15,72 +15,80 @@
#else
// Set time step for advancing from time t to time t+dt.
+template<typename quadrature_type>
inline
void
-pylith::feassemble::Integrator::timeStep(const double dt) {
+pylith::feassemble::Integrator<quadrature_type>::timeStep(const double dt) {
_dt = dt;
} // timeStep
// Check whether Jacobian needs to be recomputed.
+template<typename quadrature_type>
inline
bool
-pylith::feassemble::Integrator::needNewJacobian(void) const {
+pylith::feassemble::Integrator<quadrature_type>::needNewJacobian(void) const {
return _needNewJacobian;
} // needNewJacobian
// Set flag for setting constraints for total field solution or
+template<typename quadrature_type>
inline
void
-pylith::feassemble::Integrator::useSolnIncr(const bool flag) {
+pylith::feassemble::Integrator<quadrature_type>::useSolnIncr(const bool flag) {
_useSolnIncr = flag;
} // useSolnIncr
// Integrate contributions to residual term (r) for operator.
+template<typename quadrature_type>
inline
void
-pylith::feassemble::Integrator::integrateResidual(
- const topology::Field& residual,
- const double t,
- topology::SolutionFields* const fields) {
+pylith::feassemble::Integrator<quadrature_type>::integrateResidual(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields) {
} // integrateResidual
// Integrate contributions to Jacobian matrix (A) associated with
// operator.
+template<typename quadrature_type>
inline
void
-pylith::feassemble::Integrator::integrateJacobian(
- PetscMat* mat,
- const double t,
- topology::SolutionFields* const fields) {
+pylith::feassemble::Integrator<quadrature_type>::integrateJacobian(
+ PetscMat* mat,
+ const double t,
+ topology::SolutionFields* const fields) {
_needNewJacobian = false;
} // integrateJacobian
// Integrate contributions to residual term (r) for operator that
// do not require assembly over cells, vertices, or processors.
+template<typename quadrature_type>
inline
void
-pylith::feassemble::Integrator::integrateResidualAssembled(
- const topology::Field& residual,
- const double t,
- topology::SolutionFields* const fields) {
+pylith::feassemble::Integrator<quadrature_type>::integrateResidualAssembled(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields) {
_needNewJacobian = false;
} // integrateResidualAssembled
// Integrate contributions to Jacobian matrix (A) associated with
// operator that do not require assembly over cells, vertices, or
// processors
+template<typename quadrature_type>
inline
void
-pylith::feassemble::Integrator::integrateJacobianAssembled(
- PetscMat* mat,
- const double t,
- topology::SolutionFields* const fields) {
+pylith::feassemble::Integrator<quadrature_type>::integrateJacobianAssembled(
+ PetscMat* mat,
+ const double t,
+ topology::SolutionFields* const fields) {
} // integrateJacobianAssembled
// Update state variables as needed.
+template<typename quadrature_type>
inline
void
-pylith::feassemble::Integrator::updateState(
+pylith::feassemble::Integrator<quadrature_type>::updateState(
const double t,
topology::SolutionFields* const fields) {
} // updateState
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.cc 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.cc 2009-02-14 00:55:30 UTC (rev 14051)
@@ -79,9 +79,8 @@
// Update state variables as needed.
void
pylith::feassemble::IntegratorElasticity::updateState(
- const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh)
+ const double t,
+ topology::Solution* const fields)
{ // updateState
assert(0 != _quadrature);
assert(0 != _material);
@@ -112,7 +111,7 @@
// Get cell information
const int materialId = _material->id();
- const ALE::Obj<Mesh::label_sequence>& cells =
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
mesh->getLabelStratum("material-id", materialId);
assert(!cells.isNull());
const Mesh::label_sequence::iterator cellsEnd = cells->end();
@@ -135,7 +134,6 @@
totalStrain = 0.0;
const ALE::Obj<real_section_type>& disp = fields->getSolution();
- /// const int dispAtlasTag = fields->getSolutionAtlasTag(materialId);
// Loop over cells
int c_index = 0;
@@ -143,7 +141,7 @@
c_iter != cellsEnd;
++c_iter, ++c_index) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
+ _quadrature->retrieveGeometry(*c_iter);
// Restrict input fields to cell
mesh->restrictClosure(disp, *c_iter, &dispCell[0], cellVecSize);
@@ -163,7 +161,7 @@
// Verify configuration is acceptable.
void
pylith::feassemble::IntegratorElasticity::verifyConfiguration(
- const ALE::Obj<Mesh>& mesh) const
+ const topology::Mesh& mesh) const
{ // verifyConfiguration
assert(0 != _quadrature);
assert(0 != _material);
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.hh 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.hh 2009-02-14 00:55:30 UTC (rev 14051)
@@ -21,21 +21,17 @@
#define pylith_feassemble_integratorelasticity_hh
#include "Integrator.hh" // ISA Integrator
+#include "pylith/topology/Mesh.hh" // ISA Integrator<Mesh>
+
#include "pylith/utils/array.hh" // USES std::vector, double_array
-#include "pylith/utils/vectorfields.hh" // USES VectorFieldEnum
namespace pylith {
- namespace feassemble {
- class IntegratorElasticity;
- class TestIntegratorElasticity;
- } // feassemble
-
namespace materials {
class ElasticMaterial;
} // feassemble
} // pylith
-class pylith::feassemble::IntegratorElasticity : public Integrator
+class pylith::feassemble::IntegratorElasticity : public Integrator<topology::Mesh>
{ // IntegratorElasticity
friend class TestIntegratorElasticity; // unit testing
@@ -84,14 +80,13 @@
* @param mesh Finite-element mesh
*/
void updateState(const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ topology::SolutionFields* const fields);
/** Verify configuration is acceptable.
*
* @param mesh Finite-element mesh
*/
- void verifyConfiguration(const ALE::Obj<Mesh>& mesh) const;
+ void verifyConfiguration(const topology::Mesh& mesh) const;
/** Get cell field associated with integrator.
*
@@ -101,11 +96,9 @@
* @param fields Fields manager.
* @returns Cell field.
*/
- const ALE::Obj<real_section_type>&
- cellField(VectorFieldEnum* fieldType,
- const char* name,
- const ALE::Obj<Mesh>& mesh,
- topology::FieldsManager* const fields);
+ const topology::Field<topology::Mesh>&
+ cellField(const char* name,
+ topology::SolutionFields* const fields);
// PROTECTED METHODS ////////////////////////////////////////////////////
@@ -231,16 +224,16 @@
materials::ElasticMaterial* _material;
/// Buffer for storing scalar cell field.
- ALE::Obj<real_section_type> _bufferCellScalar;
+ Field<topology::Mesh>* _bufferCellScalar;
/// Buffer for storing vector cell field.
- ALE::Obj<real_section_type> _bufferCellVector;
+ Field<topology::Mesh>* _bufferCellVector;
/// Buffer for storing cell tensor field.
- ALE::Obj<real_section_type> _bufferCellTensor;
+ Field<topology::Mesh>* _bufferCellTensor;
/// Buffer for storing other cell fields.
- ALE::Obj<real_section_type> _bufferCellOther;
+ Field<topology::Mesh>* _bufferCellOther;
}; // IntegratorElasticity
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Makefile.am 2009-02-14 00:55:30 UTC (rev 14051)
@@ -22,6 +22,7 @@
ElasticityImplicit.hh \
Integrator.hh \
Integrator.icc \
+ Integrator.cc \
IntegratorElasticity.hh \
GeometryPoint1D.hh \
GeometryPoint2D.hh \
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh 2009-02-14 00:55:30 UTC (rev 14051)
@@ -41,6 +41,7 @@
typedef Mesh::SieveSubMesh SieveMesh;
typedef Mesh::RealSection RealSection;
typedef Mesh::IntSection IntSection;
+ typedef Mesh::RestrictVisitor RestrictVisitor;
// Sieve mesh for higher level domain (mesh, not submesh)
typedef Mesh::SieveMesh DomainSieveMesh;
Modified: short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestIntegrator.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestIntegrator.hh 2009-02-13 23:06:42 UTC (rev 14050)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/feassemble/TestIntegrator.hh 2009-02-14 00:55:30 UTC (rev 14051)
@@ -27,17 +27,9 @@
namespace pylith {
namespace feassemble {
class TestIntegrator;
-
- class Quadrature1D; // USES Quadrature1D
} // feassemble
} // pylith
-namespace spatialdata {
- namespace spatialdb {
- class GravityField; // USES GravityField
- } // spatialdb
-} // spatialdata
-
/// C++ unit testing for Integrator
class pylith::feassemble::TestIntegrator : public CppUnit::TestFixture
{ // class TestIntegrator
@@ -84,15 +76,6 @@
/// Test _resetCellMatrix()
void testResetCellMatrix(void);
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
-
- /** Initialize 1-D quadrature object.
- *
- * @param quadrature Quadrature object
- */
- void _initQuadrature(Quadrature1D* quadrature);
-
}; // class TestIntegrator
#endif // pylith_feassemble_testintegrator_hh
More information about the CIG-COMMITS
mailing list