[cig-commits] r14252 - in short/3D/PyLith/branches/pylith-swig: libsrc libsrc/feassemble libsrc/topology libsrc/utils pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Sat Mar 7 21:05:24 PST 2009
Author: brad
Date: 2009-03-07 21:05:24 -0800 (Sat, 07 Mar 2009)
New Revision: 14252
Added:
short/3D/PyLith/branches/pylith-swig/libsrc/utils/utilsfwd.hh
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.hh
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityImplicit.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/topology/Mesh.hh
short/3D/PyLith/branches/pylith-swig/libsrc/utils/Makefile.am
short/3D/PyLith/branches/pylith-swig/pylith/problems/Explicit.py
Log:
Worked on updating elasticity integrators.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-03-08 05:05:24 UTC (rev 14252)
@@ -51,6 +51,8 @@
feassemble/Quadrature2Din3D.cc \
feassemble/Quadrature3D.cc \
feassemble/IntegratorElasticity.cc \
+ feassemble/ElasticityImplicit.cc \
+ feassemble/ElasticityExplicit.cc \
materials/Metadata.cc \
materials/Material.cc \
materials/ElasticMaterial.cc \
@@ -90,8 +92,6 @@
# faults/LiuCosSlipFn.cc \
# faults/SlipTimeFn.cc \
# faults/StepSlipFn.cc \
-# feassemble/ElasticityExplicit.cc \
-# feassemble/ElasticityImplicit.cc \
# materials/GenMaxwellIsotropic3D.cc \
# meshio/CellFilter.cc \
# meshio/CellFilterAvg.cc \
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.cc 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.cc 2009-03-08 05:05:24 UTC (rev 14252)
@@ -18,17 +18,25 @@
#include "CellGeometry.hh" // USES CellGeometry
#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
-#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
-#include "pylith/utils/array.hh" // USES double_array
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+
+#include "pylith/utils/array.hh" // USES double_array
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+#include "pylith/utils/lapack.h" // USES LAPACKdgesvd
#include "petscmat.h" // USES PetscMat
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimendional
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityExplicit::ElasticityExplicit(void) :
_dtm1(-1.0)
@@ -72,11 +80,9 @@
// Integrate constributions to residual term (r) for operator.
void
pylith::feassemble::ElasticityExplicit::integrateResidual(
- const ALE::Obj<real_section_type>& residual,
- const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs)
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
{ // integrateResidual
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityExplicit::*elasticityResidual_fn_type)
@@ -84,31 +90,53 @@
assert(0 != _quadrature);
assert(0 != _material);
- assert(!residual.isNull());
+ assert(0 != _logger);
assert(0 != fields);
- assert(!mesh.isNull());
- PetscErrorCode err = 0;
+ const int setupEvent = _logger->eventId("ElIR setup");
+ const int geometryEvent = _logger->eventId("ElIR geometry");
+ const int computeEvent = _logger->eventId("ElIR compute");
+ const int restrictEvent = _logger->eventId("ElIR restrict");
+ const int stateVarsEvent = _logger->eventId("ElIR stateVars");
+ const int stressEvent = _logger->eventId("ElIR stress");
+ const int updateEvent = _logger->eventId("ElIR update");
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const int tensorSize = _material->tensorSize();
+ /** :TODO:
+ *
+ * If cellDim and spaceDim are different, we need to transform
+ * displacements into cellDim, compute action, and transform result
+ * back into spaceDim. We get this information from the Jacobian and
+ * inverse of the Jacobian.
+ */
+ if (cellDim != spaceDim)
+ throw std::logic_error("Integration for cells with spatial dimensions "
+ "different than the spatial dimension of the "
+ "domain not implemented yet.");
+
// Set variables dependent on dimension of cell
- const int cellDim = _quadrature->cellDim();
- int tensorSize = 0;
totalStrain_fn_type calcTotalStrainFn;
elasticityResidual_fn_type elasticityResidualFn;
if (1 == cellDim) {
- tensorSize = 1;
elasticityResidualFn =
&pylith::feassemble::ElasticityExplicit::_elasticityResidual1D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
} else if (2 == cellDim) {
- tensorSize = 3;
elasticityResidualFn =
&pylith::feassemble::ElasticityExplicit::_elasticityResidual2D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
} else if (3 == cellDim) {
- tensorSize = 6;
elasticityResidualFn =
&pylith::feassemble::ElasticityExplicit::_elasticityResidual3D;
calcTotalStrainFn =
@@ -116,72 +144,76 @@
} else
assert(0);
+ // Allocate vectors for cell values.
+ double_array dispTCell(numBasis*spaceDim);
+ double_array dispTmdtCell(numBasis*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+ double_array gravVec(spaceDim);
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+
// Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
const int materialId = _material->id();
- const ALE::Obj<Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", materialId);
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- const ALE::Obj<real_section_type>& dispT = fields->getFieldByHistory(1);
- assert(!dispT.isNull());
- const ALE::Obj<real_section_type>& dispTmdt = fields->getFieldByHistory(2);
- assert(!dispTmdt.isNull());
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp t").section();
+ assert(!dispTSection.isNull());
+ topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
+ numBasis*spaceDim,
+ &dispTCell[0]);
+ const ALE::Obj<RealSection>& dispTmdtSection =
+ fields->get("disp t-dt").section();
+ assert(!dispTmdtSection.isNull());
+ topology::Mesh::RestrictVisitor dispTmdtVisitor(*dispTmdtSection,
+ numBasis*spaceDim,
+ &dispTmdtCell[0]);
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
+ &_cellVector[0]);
+ assert(0 != _normalizer);
+ const double lengthScale = _normalizer->lengthScale();
+ const double gravityScale =
+ _normalizer->pressureScale() / (_normalizer->lengthScale() *
+ _normalizer->densityScale());
+
// Get parameters used in integration.
const double dt = _dt;
const double dt2 = dt*dt;
assert(dt > 0);
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const double_array& quadWts = _quadrature->quadWts();
- assert(quadWts.size() == numQuadPts);
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
+ _logger->eventEnd(setupEvent);
- /** :TODO:
- *
- * If cellDim and spaceDim are different, we need to transform
- * displacements into cellDim, compute action, and transform result
- * back into spaceDim. We get this information from the Jacobian and
- * inverse of the Jacobian.
- */
- if (cellDim != spaceDim)
- throw std::logic_error("Not implemented yet.");
-
- // Precompute the geometric and function space information
- _quadrature->precomputeGeometry(mesh, coordinates, cells);
-
- // Allocate vectors for cell values.
- _initCellVector();
- const int cellVecSize = numBasis*spaceDim;
- double_array dispTCell(cellVecSize);
- double_array dispTmdtCell(cellVecSize);
-
- // Allocate vector for total strain
- double_array totalStrain(numQuadPts*tensorSize);
- totalStrain = 0.0;
-
- int c_index = 0;
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter, ++c_index) {
+ ++c_iter) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
+ _logger->eventBegin(geometryEvent);
+ _quadrature->retrieveGeometry(*c_iter);
+ _logger->eventEnd(geometryEvent);
// Get state variables for cell.
- _material->getPropertiesCell(*c_iter, numQuadPts);
+ _logger->eventBegin(stateVarsEvent);
+ _material->retrievePropsAndVars(*c_iter);
+ _logger->eventEnd(stateVarsEvent);
// Reset element vector to zero
_resetCellVector();
- mesh->restrictClosure(dispT, *c_iter, &dispTCell[0], cellVecSize);
- mesh->restrictClosure(dispTmdt, *c_iter, &dispTmdtCell[0], cellVecSize);
+ // Restrict input fields to cell
+ _logger->eventBegin(restrictEvent);
+ dispTVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTVisitor);
+ dispTmdtVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTmdtVisitor);
+ _logger->eventEnd(restrictEvent);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -189,6 +221,7 @@
const double_array& jacobianDet = _quadrature->jacobianDet();
// Compute action for inertial terms
+ _logger->eventBegin(computeEvent);
const double_array& density = _material->calcDensity();
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
@@ -205,15 +238,24 @@
} // for
} // for
PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(6*spaceDim))));
+ _logger->eventEnd(computeEvent);
// Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&totalStrain, basisDeriv, dispTCell,
+ _logger->eventBegin(stressEvent);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTCell,
numBasis, numQuadPts);
- const double_array& stress = _material->calcStress(totalStrain, true);
- CALL_MEMBER_FN(*this, elasticityResidualFn)(stress);
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+ _logger->eventEnd(stressEvent);
+ _logger->eventBegin(computeEvent);
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell);
+ _logger->eventEnd(computeEvent);
+
// Assemble cell contribution into field
- mesh->updateAdd(residual, *c_iter, _cellVector);
+ _logger->eventBegin(updateEvent);
+ residualVisitor.clear();
+ sieveMesh->updateAdd(*c_iter, residualVisitor);
+ _logger->eventEnd(updateEvent);
} // for
} // integrateResidual
@@ -223,63 +265,82 @@
pylith::feassemble::ElasticityExplicit::integrateJacobian(
PetscMat* jacobian,
const double t,
- topology::FieldsManager* fields,
- const ALE::Obj<Mesh>& mesh)
+ topology::SolutionFields* fields)
{ // integrateJacobian
assert(0 != _quadrature);
assert(0 != _material);
assert(0 != jacobian);
assert(0 != fields);
- assert(!mesh.isNull());
- typedef ALE::ISieveVisitor::IndicesVisitor<Mesh::real_section_type,Mesh::order_type,PetscInt> visitor_type;
- // Get cell information
- const ALE::Obj<Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", _material->id());
- assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
+ const int setupEvent = _logger->eventId("ElIJ setup");
+ const int geometryEvent = _logger->eventId("ElIJ geometry");
+ const int computeEvent = _logger->eventId("ElIJ compute");
+ const int restrictEvent = _logger->eventId("ElIJ restrict");
+ const int stateVarsEvent = _logger->eventId("ElIJ stateVars");
+ const int updateEvent = _logger->eventId("ElIJ update");
- // Get sections
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- const ALE::Obj<real_section_type>& dispT = fields->getFieldByHistory(1);
- assert(!dispT.isNull());
+ _logger->eventBegin(setupEvent);
- // Get parameters used in integration.
- const double dt = _dt;
- const double dt2 = dt*dt;
- assert(dt > 0);
-
// Get cell geometry information that doesn't depend on cell
const int numQuadPts = _quadrature->numQuadPts();
const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
+ const int tensorSize = _material->tensorSize();
+ if (cellDim != spaceDim)
+ throw std::logic_error("Don't know how to integrate elasticity " \
+ "contribution to Jacobian matrix for cells with " \
+ "different dimensions than the spatial dimension.");
- // Precompute the geometric and function space information
- _quadrature->precomputeGeometry(mesh, coordinates, cells);
+ // Allocate vectors for cell data.
+ double_array dispTCell(numBasis*spaceDim);
- // Allocate vector for cell values (if necessary)
- _initCellMatrix();
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const int materialId = _material->id();
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+ assert(!cells.isNull());
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- const ALE::Obj<Mesh::order_type>& globalOrder = mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
+ // Get sections
+ const ALE::Obj<RealSection>& dispTSection =
+ fields->get("disp t").section();
+ assert(!dispTSection.isNull());
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+ assert(dt > 0);
+
+ const ALE::Obj<SieveMesh::order_type>& globalOrder =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispTSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- visitor_type iV(*dispT, *globalOrder, (int) pow(mesh->getSieve()->getMaxConeSize(), mesh->depth())*spaceDim);
+ topology::Mesh::IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
- int c_index = 0;
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ _logger->eventEnd(setupEvent);
+
+ // Loop over cells
+ for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter, ++c_index) {
+ ++c_iter) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
+ _logger->eventBegin(geometryEvent);
+ _quadrature->retrieveGeometry(*c_iter);
+ _logger->eventEnd(geometryEvent);
// Get state variables for cell.
- _material->getPropertiesCell(*c_iter, numQuadPts);
+ _logger->eventBegin(stateVarsEvent);
+ _material->retrievePropsAndVars(*c_iter);
+ _logger->eventEnd(stateVarsEvent);
- // Reset element vector to zero
+ // Reset element matrix to zero
_resetCellMatrix();
// Get cell geometry information that depends on cell
@@ -290,6 +351,7 @@
const double_array& density = _material->calcDensity();
// Compute Jacobian for inertial terms
+ _logger->eventBegin(computeEvent);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
@@ -306,12 +368,17 @@
} // for
} // for
PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
+ _logger->eventEnd(computeEvent);
- // Assemble cell contribution into PETSc Matrix
- PetscErrorCode err = updateOperator(*jacobian, *mesh->getSieve(), iV, *c_iter, _cellMatrix, ADD_VALUES);
+ // Assemble cell contribution into PETSc matrix.
+ _logger->eventBegin(updateEvent);
+ jacobianVisitor.clear();
+ PetscErrorCode err = updateOperator(*jacobian, *sieveMesh->getSieve(),
+ jacobianVisitor, *c_iter,
+ &_cellMatrix[0], ADD_VALUES);
if (err)
throw std::runtime_error("Update to PETSc Mat failed.");
- iV.clear();
+ _logger->eventEnd(updateEvent);
} // for
_needNewJacobian = false;
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.hh 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityExplicit.hh 2009-03-08 05:05:24 UTC (rev 14252)
@@ -55,24 +55,10 @@
#if !defined(pylith_feassemble_elasticityexplicit_hh)
#define pylith_feassemble_elasticityexplicit_hh
+// Include directives ---------------------------------------------------
#include "IntegratorElasticity.hh" // ISA IntegratorElasticity
-#include "pylith/utils/array.hh" // USES std::vector, double_array
-namespace pylith {
- namespace feassemble {
- class ElasticityExplicit;
- } // feassemble
-} // pylith
-
-namespace spatialdata {
- namespace spatialdb {
- class SpatialDB; // USES SpatialDB
- } // spatialdb
- namespace geocoords {
- class CoordSys; // USES CoordSys
- } // geocoords
-} // spatialdata
-
+// ElasticityExplicit ---------------------------------------------------
class pylith::feassemble::ElasticityExplicit : public IntegratorElasticity
{ // ElasticityExplicit
friend class TestElasticityExplicit; // unit testing
@@ -104,32 +90,27 @@
* @param residual Field containing values for residual
* @param t Current time
* @param fields Solution fields
- * @param mesh Finite-element mesh
*/
- void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ void integrateResidual(const topology::Field<topology::Mesh>& residual,
const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs);
+ topology::SolutionFields* const fields);
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
*
- * @param jacobian Sparse matrix to hold Jacobian of operator.
+ * @param jacobian Sparse matrix for Jacobian of system.
* @param t Current time
- * @param fields Solution fields.
- * @param mesh Finite-element mesh.
+ * @param fields Solution fields
*/
void integrateJacobian(PetscMat* jacobian,
const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ topology::SolutionFields* const fields);
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
/// Not implemented.
- ElasticityExplicit(const ElasticityExplicit& i);
+ ElasticityExplicit(const ElasticityExplicit&);
/// Not implemented
const ElasticityExplicit& operator=(const ElasticityExplicit&);
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityImplicit.cc 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityImplicit.cc 2009-03-08 05:05:24 UTC (rev 14252)
@@ -18,7 +18,10 @@
#include "CellGeometry.hh" // USES CellGeometry
#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
-#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+
+#include "pylith/utils/EventLogger.hh" // USES EventLogger
#include "pylith/utils/array.hh" // USES double_array
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
#include "pylith/utils/lapack.h" // USES LAPACKdgesvd
@@ -32,6 +35,10 @@
#include <stdexcept> // USES std::runtime_error
// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityImplicit::ElasticityImplicit(void) :
_dtm1(-1.0)
@@ -82,11 +89,9 @@
// Integrate constributions to residual term (r) for operator.
void
pylith::feassemble::ElasticityImplicit::integrateResidual(
- const ALE::Obj<real_section_type>& residual,
- const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs)
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields)
{ // integrateResidual
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityImplicit::*elasticityResidual_fn_type)
@@ -94,49 +99,46 @@
assert(0 != _quadrature);
assert(0 != _material);
- assert(!residual.isNull());
+ assert(0 != _logger);
assert(0 != fields);
- assert(!mesh.isNull());
- static PetscLogEvent setupEvent = 0, cellGeomEvent = 0, stateVarsEvent = 0, restrictEvent = 0, computeEvent = 0, updateEvent = 0, stressEvent;
+ const int setupEvent = _logger->eventId("ElIR setup");
+ const int geometryEvent = _logger->eventId("ElIR geometry");
+ const int computeEvent = _logger->eventId("ElIR compute");
+ const int restrictEvent = _logger->eventId("ElIR restrict");
+ const int stateVarsEvent = _logger->eventId("ElIR stateVars");
+ const int stressEvent = _logger->eventId("ElIR stress");
+ const int updateEvent = _logger->eventId("ElIR update");
- if (!setupEvent)
- PetscLogEventRegister("IRSetup", 0, &setupEvent);
- if (!cellGeomEvent)
- PetscLogEventRegister("IRCellGeom", 0, &cellGeomEvent);
- if (!stateVarsEvent)
- PetscLogEventRegister("IRProperties", 0, &stateVarsEvent);
- if (!restrictEvent)
- PetscLogEventRegister("IRRestrict", 0, &restrictEvent);
- if (!computeEvent)
- PetscLogEventRegister("IRCompute", 0, &computeEvent);
- if (!updateEvent)
- PetscLogEventRegister("IRUpdate", 0, &updateEvent);
- if (!stressEvent)
- PetscLogEventRegister("IRMaterialStress", 0, &stressEvent);
+ _logger->eventBegin(setupEvent);
- const Obj<sieve_type>& sieve = mesh->getSieve();
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+ const int tensorSize = _material->tensorSize();
+ if (cellDim != spaceDim)
+ throw std::logic_error("Integration for cells with spatial dimensions "
+ "different than the spatial dimension of the "
+ "domain not implemented yet.");
- PetscLogEventBegin(setupEvent,0,0,0,0);
// Set variables dependent on dimension of cell
- const int cellDim = _quadrature->cellDim();
- int tensorSize = 0;
totalStrain_fn_type calcTotalStrainFn;
elasticityResidual_fn_type elasticityResidualFn;
if (1 == cellDim) {
- tensorSize = 1;
elasticityResidualFn =
&pylith::feassemble::ElasticityImplicit::_elasticityResidual1D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
} else if (2 == cellDim) {
- tensorSize = 3;
elasticityResidualFn =
&pylith::feassemble::ElasticityImplicit::_elasticityResidual2D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
} else if (3 == cellDim) {
- tensorSize = 6;
elasticityResidualFn =
&pylith::feassemble::ElasticityImplicit::_elasticityResidual3D;
calcTotalStrainFn =
@@ -144,93 +146,63 @@
} else
assert(0);
+ // Allocate vectors for cell values.
+ double_array dispTBctpdtCell(numBasis*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+ double_array gravVec(spaceDim);
+ double_array quadPtsGlobal(numQuadPts*spaceDim);
+
// Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
const int materialId = _material->id();
- const ALE::Obj<Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", materialId);
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- const ALE::Obj<real_section_type>& dispTBctpdt =
- fields->getReal("dispTBctpdt");
- assert(!dispTBctpdt.isNull());
+ const ALE::Obj<RealSection>& dispTBctpdtSection =
+ fields->get("dispTBctpdt").section();
+ assert(!dispTBctpdtSection.isNull());
+ topology::Mesh::RestrictVisitor dispTBctpdtVisitor(*dispTBctpdtSection,
+ numBasis*spaceDim,
+ &dispTBctpdtCell[0]);
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
+ &_cellVector[0]);
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const double_array& quadWts = _quadrature->quadWts();
- assert(quadWts.size() == numQuadPts);
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
-
- // Precompute the geometric and function space information
- _quadrature->precomputeGeometry(mesh, coordinates, cells);
-
- // Allocate vectors for cell values.
- _initCellVector();
- const int cellVecSize = numBasis*spaceDim;
- double_array dispTBctpdtCell(cellVecSize);
- double_array totalStrain(numQuadPts*tensorSize);
- totalStrain = 0.0;
- double_array gravVec(spaceDim);
- double_array quadPtsGlobal(numQuadPts*spaceDim);
-
- // Set up gravity field database for querying
- if (0 != _gravityField) {
- _gravityField->open();
- if (1 == spaceDim){
- const char* queryNames[] = { "x"};
- _gravityField->queryVals(queryNames, spaceDim);
- } else if (2 == spaceDim){
- const char* queryNames[] = { "x", "y"};
- _gravityField->queryVals(queryNames, spaceDim);
- } else if (3 == spaceDim){
- const char* queryNames[] = { "x", "y", "z"};
- _gravityField->queryVals(queryNames, spaceDim);
- } else {
- assert(0);
- } // else
- } // if
-
- PetscLogEventEnd(setupEvent,0,0,0,0);
-
- ALE::ISieveVisitor::RestrictVisitor<real_section_type> rV(*dispTBctpdt, cellVecSize, &dispTBctpdtCell[0]);
- if (mesh->depth() > 1) {
- //ISieveVisitor::PointRetriever<sieve_type,ISieveVisitor::RestrictVisitor<Section> > pV((int) pow((double) mesh->getSieve()->getMaxConeSize(), this->depth())+1, rV, true);
- throw ALE::Exception("Need to reorganize to use a different visitor class");
- } // if
-
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
const double gravityScale =
_normalizer->pressureScale() / (_normalizer->lengthScale() *
_normalizer->densityScale());
+ _logger->eventEnd(setupEvent);
+
// Loop over cells
- int c_index = 0;
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter, ++c_index) {
+ ++c_iter) {
// Compute geometry information for current cell
- PetscLogEventBegin(cellGeomEvent,0,0,0,0);
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
- PetscLogEventEnd(cellGeomEvent,0,0,0,0);
+ _logger->eventBegin(geometryEvent);
+ _quadrature->retrieveGeometry(*c_iter);
+ _logger->eventEnd(geometryEvent);
// Get state variables for cell.
- PetscLogEventBegin(stateVarsEvent,0,0,0,0);
- _material->getPropertiesCell(*c_iter, numQuadPts);
- PetscLogEventEnd(stateVarsEvent,0,0,0,0);
+ _logger->eventBegin(stateVarsEvent);
+ _material->retrievePropsAndVars(*c_iter);
+ _logger->eventEnd(stateVarsEvent);
// Reset element vector to zero
_resetCellVector();
// Restrict input fields to cell
- PetscLogEventBegin(restrictEvent,0,0,0,0);
- mesh->restrictClosure(*c_iter, rV);
- PetscLogEventEnd(restrictEvent,0,0,0,0);
+ _logger->eventBegin(restrictEvent);
+ dispTBctpdtVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTBctpdtVisitor);
+ _logger->eventBegin(restrictEvent);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -238,13 +210,12 @@
const double_array& jacobianDet = _quadrature->jacobianDet();
const double_array& quadPtsNondim = _quadrature->quadPts();
- if (cellDim != spaceDim)
- throw std::logic_error("Integration for cells with spatial dimensions "
- "different than the spatial dimension of the "
- "domain not implemented yet.");
-
// Compute body force vector if gravity is being used.
if (0 != _gravityField) {
+ _logger->eventBegin(computeEvent);
+ const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
+ assert(0 != cs);
+
// Get density at quadrature points for this cell
const double_array& density = _material->calcDensity();
@@ -270,19 +241,19 @@
} // for
} // for
PetscLogFlops(numQuadPts*(2+numBasis*(1+2*spaceDim)));
- _gravityField->close();
+ _logger->eventEnd(computeEvent);
} // if
// Compute B(transpose) * sigma, first computing strains
- PetscLogEventBegin(stressEvent,0,0,0,0);
- calcTotalStrainFn(&totalStrain, basisDeriv, dispTBctpdtCell,
+ _logger->eventBegin(stressEvent);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTBctpdtCell,
numBasis, numQuadPts);
- const double_array& stress = _material->calcStress(totalStrain, true);
- PetscLogEventEnd(stressEvent,0,0,0,0);
+ const double_array& stressCell = _material->calcStress(strainCell, true);
+ _logger->eventEnd(stressEvent);
- PetscLogEventBegin(computeEvent,0,0,0,0);
- CALL_MEMBER_FN(*this, elasticityResidualFn)(stress);
- PetscLogEventEnd(computeEvent,0,0,0,0);
+ _logger->eventBegin(computeEvent);
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell);
+ _logger->eventEnd(computeEvent);
#if 0
std::cout << "Updating residual for cell " << *c_iter << std::endl;
@@ -291,10 +262,10 @@
}
#endif
// Assemble cell contribution into field
- PetscLogEventBegin(updateEvent,0,0,0,0);
- mesh->updateAdd(residual, *c_iter, _cellVector);
- PetscLogEventEnd(updateEvent,0,0,0,0);
- rV.clear();
+ _logger->eventBegin(updateEvent);
+ residualVisitor.clear();
+ sieveMesh->updateAdd(*c_iter, residualVisitor);
+ _logger->eventEnd(updateEvent);
} // for
} // integrateResidual
@@ -302,10 +273,9 @@
// Compute stiffness matrix.
void
pylith::feassemble::ElasticityImplicit::integrateJacobian(
- PetscMat* mat,
+ PetscMat* jacobian,
const double t,
- topology::FieldsManager* fields,
- const ALE::Obj<Mesh>& mesh)
+ topology::SolutionFields* fields)
{ // integrateJacobian
/// Member prototype for _elasticityJacobianXD()
typedef void (pylith::feassemble::ElasticityImplicit::*elasticityJacobian_fn_type)
@@ -313,30 +283,46 @@
assert(0 != _quadrature);
assert(0 != _material);
- assert(0 != mat);
+ assert(0 != _logger);
+ assert(0 != jacobian);
assert(0 != fields);
- assert(!mesh.isNull());
- typedef ALE::ISieveVisitor::IndicesVisitor<Mesh::real_section_type,Mesh::order_type,PetscInt> visitor_type;
- // Set variables dependent on dimension of cell
+ const int setupEvent = _logger->eventId("ElIJ setup");
+ const int geometryEvent = _logger->eventId("ElIJ geometry");
+ const int computeEvent = _logger->eventId("ElIJ compute");
+ const int restrictEvent = _logger->eventId("ElIJ restrict");
+ const int stateVarsEvent = _logger->eventId("ElIJ stateVars");
+ const int updateEvent = _logger->eventId("ElIJ update");
+
+ _logger->eventBegin(setupEvent);
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
- int tensorSize = 0;
+ const int tensorSize = _material->tensorSize();
+ if (cellDim != spaceDim)
+ throw std::logic_error("Don't know how to integrate elasticity " \
+ "contribution to Jacobian matrix for cells with " \
+ "different dimensions than the spatial dimension.");
+
+ // Set variables dependent on dimension of cell
totalStrain_fn_type calcTotalStrainFn;
elasticityJacobian_fn_type elasticityJacobianFn;
if (1 == cellDim) {
- tensorSize = 1;
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian1D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
} else if (2 == cellDim) {
- tensorSize = 3;
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian2D;
calcTotalStrainFn =
&pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
} else if (3 == cellDim) {
- tensorSize = 6;
elasticityJacobianFn =
&pylith::feassemble::ElasticityImplicit::_elasticityJacobian3D;
calcTotalStrainFn =
@@ -344,85 +330,83 @@
} else
assert(0);
+ // Allocate vector for total strain
+ double_array dispTBctpdtCell(numBasis*spaceDim);
+ double_array strainCell(numQuadPts*tensorSize);
+ strainCell = 0.0;
+
// Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
const int materialId = _material->id();
- const ALE::Obj<Mesh::label_sequence>& cells =
- mesh->getLabelStratum("material-id", materialId);
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
assert(!cells.isNull());
- const Mesh::label_sequence::iterator cellsEnd = cells->end();
+ const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
// Get sections
- const ALE::Obj<real_section_type>& coordinates =
- mesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- const ALE::Obj<real_section_type>& dispTBctpdt =
- fields->getReal("dispTBctpdt");
- assert(!dispTBctpdt.isNull());
+ const ALE::Obj<RealSection>& dispTBctpdtSection =
+ fields->get("dispTBctpdt").section();
+ assert(!dispTBctpdtSection.isNull());
+ topology::Mesh::RestrictVisitor dispTBctpdtVisitor(*dispTBctpdtSection,
+ numBasis*spaceDim,
+ &dispTBctpdtCell[0]);
// Get parameters used in integration.
const double dt = _dt;
assert(dt > 0);
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const double_array& quadWts = _quadrature->quadWts();
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
-
- if (cellDim != spaceDim)
- throw std::logic_error("Don't know how to integrate elasticity " \
- "contribution to Jacobian matrix for cells with " \
- "different dimensions than the spatial dimension.");
-
- // Precompute the geometric and function space information
- _quadrature->precomputeGeometry(mesh, coordinates, cells);
-
- // Allocate matrix and vectors for cell values.
- _initCellMatrix();
- const int cellVecSize = numBasis*spaceDim;
- double_array dispTBctpdtCell(cellVecSize);
-
- // Allocate vector for total strain
- double_array totalStrain(numQuadPts*tensorSize);
- totalStrain = 0.0;
-
- const ALE::Obj<Mesh::order_type>& globalOrder = mesh->getFactory()->getGlobalOrder(mesh, "default", dispTBctpdt);
+ const ALE::Obj<SieveMesh::order_type>& globalOrder =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
+ dispTBctpdtSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- visitor_type iV(*dispTBctpdt, *globalOrder, (int) pow(mesh->getSieve()->getMaxConeSize(), mesh->depth())*spaceDim);
- ALE::ISieveVisitor::RestrictVisitor<real_section_type> rV(*dispTBctpdt, cellVecSize, &dispTBctpdtCell[0]);
+ topology::Mesh::IndicesVisitor jacobianVisitor(*dispTBctpdtSection,
+ *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
+ _logger->eventEnd(setupEvent);
+
// Loop over cells
- int c_index = 0;
- for (Mesh::label_sequence::iterator c_iter=cells->begin();
+ for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter, ++c_index) {
+ ++c_iter) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
+ _logger->eventBegin(geometryEvent);
+ _quadrature->retrieveGeometry(*c_iter);
+ _logger->eventEnd(geometryEvent);
// Get state variables for cell.
- _material->getPropertiesCell(*c_iter, numQuadPts);
+ _logger->eventBegin(stateVarsEvent);
+ _material->retrievePropsAndVars(*c_iter);
+ _logger->eventEnd(stateVarsEvent);
// Reset element matrix to zero
_resetCellMatrix();
// Restrict input fields to cell
- mesh->restrictClosure(*c_iter, rV);
+ _logger->eventBegin(restrictEvent);
+ dispTBctpdtVisitor.clear();
+ sieveMesh->restrictClosure(*c_iter, dispTBctpdtVisitor);
+ _logger->eventBegin(restrictEvent);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
const double_array& basisDeriv = _quadrature->basisDeriv();
const double_array& jacobianDet = _quadrature->jacobianDet();
+ _logger->eventBegin(computeEvent);
// Compute strains
- calcTotalStrainFn(&totalStrain, basisDeriv, dispTBctpdtCell,
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTBctpdtCell,
numBasis, numQuadPts);
// Get "elasticity" matrix at quadrature points for this cell
const double_array& elasticConsts =
- _material->calcDerivElastic(totalStrain);
+ _material->calcDerivElastic(strainCell);
CALL_MEMBER_FN(*this, elasticityJacobianFn)(elasticConsts);
+ _logger->eventEnd(computeEvent);
if (_quadrature->checkConditioning()) {
int n = numBasis*spaceDim;
@@ -446,7 +430,7 @@
throw std::runtime_error("Lapack SVD failed");
minSV = svalues[n-7];
maxSV = svalues[0];
- std::cout << "Element " << c_index << std::endl;
+ std::cout << "Element " << *c_iter << std::endl;
for(int i = 0; i < n; ++i)
std::cout << " sV["<<i<<"] = " << svalues[i] << std::endl;
std::cout << " kappa(elemMat) = " << maxSV/minSV << std::endl;
@@ -455,13 +439,15 @@
delete [] work;
} // if
- // Assemble cell contribution into field. Not sure if this is correct for
- // global stiffness matrix.
- PetscErrorCode err = updateOperator(*mat, *mesh->getSieve(), iV, *c_iter, _cellMatrix, ADD_VALUES);
+ // Assemble cell contribution into PETSc matrix.
+ _logger->eventBegin(updateEvent);
+ jacobianVisitor.clear();
+ PetscErrorCode err = updateOperator(*jacobian, *sieveMesh->getSieve(),
+ jacobianVisitor, *c_iter,
+ &_cellMatrix[0], ADD_VALUES);
if (err)
throw std::runtime_error("Update to PETSc Mat failed.");
- iV.clear();
- rV.clear();
+ _logger->eventEnd(updateEvent);
} // for
_needNewJacobian = false;
_material->resetNeedNewJacobian();
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityImplicit.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityImplicit.hh 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/ElasticityImplicit.hh 2009-03-08 05:05:24 UTC (rev 14252)
@@ -95,35 +95,30 @@
* external loads due to body forces plus the
* element internal forces for the current stress state.
*
- * @param residual Residual field (output)
+ * @param residual Field containing values for residual
* @param t Current time
* @param fields Solution fields
- * @param mesh Mesh object
- * @param cs Coordinate system
*/
- void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ void integrateResidual(const topology::Field<topology::Mesh>& residual,
const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs);
+ topology::SolutionFields* const fields);
- /** Compute Jacobian matrix (A) associated with operator.
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
*
- * @param mat Sparse matrix
+ * @param jacobian Sparse matrix for Jacobian of system.
* @param t Current time
* @param fields Solution fields
- * @param mesh Mesh object
*/
- void integrateJacobian(PetscMat* mat,
+ void integrateJacobian(PetscMat* jacobian,
const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ topology::SolutionFields* const fields);
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
/// Not implemented.
- ElasticityImplicit(const ElasticityImplicit& i);
+ ElasticityImplicit(const ElasticityImplicit&);
/// Not implemented
const ElasticityImplicit& operator=(const ElasticityImplicit&);
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.cc 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.cc 2009-03-08 05:05:24 UTC (rev 14252)
@@ -13,6 +13,7 @@
#include <portinfo>
#include "Quadrature.hh" // USES Quadrature
+#include "pylith/utils/EventLogger.hh" // USES EventLogger
#include "pylith/utils/constdefs.h" // USES MAXDOUBLE
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -28,6 +29,7 @@
_quadrature(0),
_normalizer(new spatialdata::units::Nondimensional),
_gravityField(0),
+ _logger(0),
_needNewJacobian(true),
_useSolnIncr(false)
{ // constructor
@@ -40,6 +42,8 @@
{ // destructor
delete _quadrature; _quadrature = 0;
delete _normalizer; _normalizer = 0;
+ delete _logger; _logger = 0;
+ _gravityField = 0; /// Memory managed elsewhere :TODO: use shared pointer
} // destructor
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.hh 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.hh 2009-03-08 05:05:24 UTC (rev 14252)
@@ -28,6 +28,8 @@
#include "feassemblefwd.hh" // forward declarations
#include "pylith/topology/topologyfwd.hh" // USES Mesh, Field, SolutionFields
+#include "pylith/utils/utilsfwd.hh" // HOLDSA EventLogger
+
#include "spatialdata/spatialdb/spatialdbfwd.hh" // USES GravityField
#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
@@ -100,6 +102,13 @@
virtual
void useSolnIncr(const bool flag);
+ /** Initialize integrator.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ virtual
+ void initialize(const topology::Mesh& mesh);
+
/** Integrate contributions to residual term (r) for operator.
*
* @param residual Field containing values for residual
@@ -190,6 +199,8 @@
spatialdata::units::Nondimensional* _normalizer; ///< Nondimensionalizer.
spatialdata::spatialdb::GravityField* _gravityField; ///< Gravity field.
+ utils::EventLogger* _logger; ///< Event logger.
+
/// Vector local to cell containing result of integration action
double_array _cellVector;
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.icc 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/Integrator.icc 2009-03-08 05:05:24 UTC (rev 14252)
@@ -38,6 +38,13 @@
_useSolnIncr = flag;
} // useSolnIncr
+// Initialize integrator.
+template<typename quadrature_type>
+inline
+void
+pylith::feassemble::Integrator<quadrature_type>::initialize(const topology::Mesh& mesh) {
+} // initialize
+
// Integrate contributions to residual term (r) for operator.
template<typename quadrature_type>
inline
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.cc 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.cc 2009-03-08 05:05:24 UTC (rev 14252)
@@ -20,6 +20,7 @@
#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/utils/EventLogger.hh" // USES EventLogger
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -83,6 +84,70 @@
} // useSolnIncr
// ----------------------------------------------------------------------
+// Initialize integrator.
+void
+pylith::feassemble::IntegratorElasticity::initialize(const topology::Mesh& mesh)
+{ // initialize
+ assert(0 != _quadrature);
+ assert(0 != _material);
+
+ // Get cell information
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const int materialId = _material->id();
+ const ALE::Obj<SieveMesh::label_sequence>& cells =
+ sieveMesh->getLabelStratum("material-id", materialId);
+
+ // Compute geometry for quadrature operations.
+ _quadrature->computeGeometry(mesh, cells);
+
+ // Initialize material.
+ _material->initialize(mesh, _quadrature);
+
+ // Allocate vectors and matrices for cell values.
+ _initCellVector();
+ _initCellMatrix();
+
+ // Setup event logger.
+ _logger = new utils::EventLogger;
+ assert(0 != _logger);
+ _logger->className("ElasticityIntegrator");
+ _logger->initialize();
+ _logger->registerEvent("ElIR setup");
+ _logger->registerEvent("ElIR geometry");
+ _logger->registerEvent("ElIR compute");
+ _logger->registerEvent("ElIR restrict");
+ _logger->registerEvent("ElIR stateVars");
+ _logger->registerEvent("ElIR stress");
+ _logger->registerEvent("ElIR update");
+
+ _logger->registerEvent("ElIJ setup");
+ _logger->registerEvent("ElIJ geometry");
+ _logger->registerEvent("ElIJ compute");
+ _logger->registerEvent("ElIJ restrict");
+ _logger->registerEvent("ElIJ stateVars");
+ _logger->registerEvent("ElIJ update");
+
+ // Set up gravity field database for querying
+ if (0 != _gravityField) {
+ const int spaceDim = _quadrature->spaceDim();
+ _gravityField->open();
+ if (1 == spaceDim){
+ const char* queryNames[] = { "x"};
+ _gravityField->queryVals(queryNames, spaceDim);
+ } else if (2 == spaceDim){
+ const char* queryNames[] = { "x", "y"};
+ _gravityField->queryVals(queryNames, spaceDim);
+ } else if (3 == spaceDim){
+ const char* queryNames[] = { "x", "y", "z"};
+ _gravityField->queryVals(queryNames, spaceDim);
+ } else {
+ assert(0);
+ } // else
+ } // if
+} // initialize
+
+// ----------------------------------------------------------------------
// Update state variables as needed.
void
pylith::feassemble::IntegratorElasticity::updateStateVars(
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.hh 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/feassemble/IntegratorElasticity.hh 2009-03-08 05:05:24 UTC (rev 14252)
@@ -75,6 +75,12 @@
*/
void useSolnIncr(const bool flag);
+ /** Initialize integrator.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ void initialize(const topology::Mesh& mesh);
+
/** Update state variables as needed.
*
* @param t Current time
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh 2009-03-08 05:05:24 UTC (rev 14252)
@@ -51,6 +51,7 @@
typedef ALE::IMesh<ALE::LabelSifter<int, SieveMesh::point_type> > SieveSubMesh;
typedef ALE::ISieveVisitor::RestrictVisitor<RealSection> RestrictVisitor;
typedef ALE::ISieveVisitor::UpdateAddVisitor<RealSection> UpdateAddVisitor;
+ typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveMesh::order_type,PetscInt> IndicesVisitor;
// PUBLIC METHODS ///////////////////////////////////////////////////////
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/utils/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/utils/Makefile.am 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/utils/Makefile.am 2009-03-08 05:05:24 UTC (rev 14252)
@@ -24,7 +24,8 @@
petscfwd.h \
sievefwd.hh \
sievetypes.hh \
- vectorfields.hh
+ vectorfields.hh \
+ utilsfwd.hh
noinst_HEADERS =
Added: short/3D/PyLith/branches/pylith-swig/libsrc/utils/utilsfwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/utils/utilsfwd.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/utils/utilsfwd.hh 2009-03-08 05:05:24 UTC (rev 14252)
@@ -0,0 +1,36 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/** @file libsrc/utils/utilsfwd.hh
+ *
+ * @brief Forward declarations for PyLith utils objects.
+ *
+ * Including this header file eliminates the need to use separate
+ * forward declarations.
+ */
+
+#if !defined(pylith_utils_utilsfwd_hh)
+#define pylith_utils_utilsfwd_hh
+
+namespace pylith {
+ namespace utils {
+
+ class EventLogger;
+
+ } // utils
+} // pylith
+
+
+#endif // pylith_utils_utilsfwd_hh
+
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-swig/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/problems/Explicit.py 2009-03-08 03:22:52 UTC (rev 14251)
+++ short/3D/PyLith/branches/pylith-swig/pylith/problems/Explicit.py 2009-03-08 05:05:24 UTC (rev 14252)
@@ -101,6 +101,7 @@
needNewJacobian = False
for integrator in self.integrators:
+ integrator.timeStep(dt)
if integrator.needNewJacobian():
needNewJacobian = True
if needNewJacobian:
More information about the CIG-COMMITS
mailing list