[cig-commits] r21724 - short/3D/PyLith/trunk/libsrc/pylith/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Apr 4 16:35:43 PDT 2013
Author: brad
Date: 2013-04-04 16:35:43 -0700 (Thu, 04 Apr 2013)
New Revision: 21724
Modified:
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.hh
Log:
Code cleanup (feassemble explicit).
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2013-04-04 21:47:27 UTC (rev 21723)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2013-04-04 23:35:43 UTC (rev 21724)
@@ -33,7 +33,6 @@
#include "pylith/utils/array.hh" // USES scalar_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
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2013-04-04 21:47:27 UTC (rev 21723)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2013-04-04 23:35:43 UTC (rev 21724)
@@ -27,10 +27,12 @@
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
#include "pylith/utils/array.hh" // USES scalar_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
@@ -70,6 +72,8 @@
void
pylith::feassemble::ElasticityExplicitLgDeform::timeStep(const PylithScalar dt)
{ // timeStep
+ PYLITH_METHOD_BEGIN;
+
if (_dt != -1.0)
_dtm1 = _dt;
else
@@ -78,6 +82,8 @@
assert(_dt == _dtm1); // For now, don't allow variable time step
if (0 != _material)
_material->timeStep(_dt);
+
+ PYLITH_METHOD_END;
} // timeStep
// ----------------------------------------------------------------------
@@ -85,8 +91,10 @@
PylithScalar
pylith::feassemble::ElasticityExplicitLgDeform::stableTimeStep(const topology::Mesh& mesh) const
{ // stableTimeStep
+ PYLITH_METHOD_BEGIN;
+
assert(_material);
- return _material->stableTimeStepExplicit(mesh, _quadrature);
+ PYLITH_METHOD_RETURN(_material->stableTimeStepExplicit(mesh, _quadrature));
} // stableTimeStep
// ----------------------------------------------------------------------
@@ -94,6 +102,8 @@
void
pylith::feassemble::ElasticityExplicitLgDeform::normViscosity(const PylithScalar viscosity)
{ // normViscosity
+ PYLITH_METHOD_BEGIN;
+
if (viscosity < 0.0) {
std::ostringstream msg;
msg << "Normalized viscosity (" << viscosity << ") must be nonnegative.";
@@ -101,24 +111,27 @@
} // if
_normViscosity = viscosity;
+
+ PYLITH_METHOD_END;
} // normViscosity
// ----------------------------------------------------------------------
// Integrate constributions to residual term (r) for operator.
void
-pylith::feassemble::ElasticityExplicitLgDeform::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::feassemble::ElasticityExplicitLgDeform::integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateResidualLumped
+ PYLITH_METHOD_BEGIN;
+
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityExplicitLgDeform::*elasticityResidual_fn_type)
(const scalar_array&, const scalar_array&);
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != _logger);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_material);
+ assert(_logger);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIR setup");
const int geometryEvent = _logger->eventId("ElIR geometry");
@@ -178,54 +191,38 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- DM dmMesh = fields->mesh().dmMesh();
- IS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- PetscErrorCode err;
+ PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
+ topology::StratumIS materialIS(dmMesh, "material-id", _material->id());
+ const PetscInt* cells = materialIS.points();
+ const PetscInt numCells = materialIS.size();
- assert(dmMesh);
- err = DMPlexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
- err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
- err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ // Setup field visitors.
+ topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
+ PetscScalar* accCell = NULL;
+ PetscInt accSize = 0;
- // Get sections
- topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
- PetscSection accSection = acc.petscSection();
- Vec accVec = acc.localVector();
- assert(accSection);assert(accVec);
+ topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
+ PetscScalar* velCell = NULL;
+ PetscInt velSize = 0;
scalar_array dispAdjCell(numBasis*spaceDim);
- topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
- PetscSection velSection = vel.petscSection();
- Vec velVec = vel.localVector();
- assert(velSection);assert(velVec);
+ topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- PetscSection dispTSection = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
- assert(dispTSection);assert(dispTVec);
+ topology::VecVisitorMesh residualVisitor(residual);
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
-
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar *coordsCell = NULL;
+ PetscInt coordsSize = 0;
- assert(0 != _normalizer);
+ assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
- const PylithScalar gravityScale =
- _normalizer->pressureScale() / (_normalizer->lengthScale() *
- _normalizer->densityScale());
+ const PylithScalar gravityScale = _normalizer->pressureScale() / (_normalizer->lengthScale() * _normalizer->densityScale());
- const PylithScalar dt = _dt;
- assert(_normViscosity >= 0.0);
- assert(dt > 0);
- const PylithScalar viscosity = dt*_normViscosity;
+ const PylithScalar dt = _dt;assert(dt > 0);
+ const PylithScalar viscosity = dt*_normViscosity;assert(_normViscosity >= 0.0);
// Get parameters used in integration.
scalar_array valuesIJ(numBasis);
@@ -240,12 +237,10 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(cell);
#else
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- _quadrature->computeGeometry(coordinatesCell, cell);
- err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsCell[i];} // :TODO: Remove copy
+ _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
#endif
// Get state variables for cell.
@@ -255,13 +250,11 @@
_resetCellVector();
// Restrict input fields to cell
- PetscScalar *accArray, *velArray, *dispTArray;
- PetscInt accSize, velSize, dispTSize;
- err = DMPlexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- assert(velSize == accSize);
- err = DMPlexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- assert(dispTSize == accSize);
+ accVisitor.getClosure(&accCell, &accSize, cell);
+ velVisitor.getClosure(&velCell, &velSize, cell);
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);
+ assert(velSize == accSize);
+ assert(dispSize == accSize);
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -270,29 +263,25 @@
const scalar_array& quadPtsNondim = _quadrature->quadPts();
// Compute body force vector if gravity is being used.
- if (0 != _gravityField) {
- const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
- assert(0 != cs);
+ if (_gravityField) {
+ const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();assert(cs);
// Get density at quadrature points for this cell
const scalar_array& density = _material->calcDensity();
quadPtsGlobal = quadPtsNondim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
// Compute action for element body forces
spatialdata::spatialdb::SpatialDB* db = _gravityField;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const int err = db->query(&gravVec[0], gravVec.size(),
- &quadPtsGlobal[0], spaceDim, cs);
- if (err)
+ const int err = db->query(&gravVec[0], gravVec.size(), &quadPtsGlobal[0], spaceDim, cs);
+ if (err) {
throw std::runtime_error("Unable to get gravity vector for point.");
- _normalizer->nondimensionalize(&gravVec[0], gravVec.size(),
- gravityScale);
+ } // if
+ _normalizer->nondimensionalize(&gravVec[0], gravVec.size(), gravityScale);
const PylithScalar wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
- for (int iBasis=0, iQ=iQuad*numBasis;
- iBasis < numBasis; ++iBasis) {
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const PylithScalar valI = wt*basis[iQ+iBasis];
for (int iDim=0; iDim < spaceDim; ++iDim) {
_cellVector[iBasis*spaceDim+iDim] += valI*gravVec[iDim];
@@ -309,28 +298,30 @@
const PylithScalar wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
const int iQ = iQuad * numBasis;
PylithScalar valJ = 0.0;
- for (int jBasis = 0; jBasis < numBasis; ++jBasis)
+ for (int jBasis = 0; jBasis < numBasis; ++jBasis) {
valJ += basis[iQ + jBasis];
+ } // for
valJ *= wt;
- for (int iBasis = 0; iBasis < numBasis; ++iBasis)
+ for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
valuesIJ[iBasis] += basis[iQ + iBasis] * valJ;
+ } // for
} // for
- for (int iBasis = 0; iBasis < numBasis; ++iBasis)
- for (int iDim = 0; iDim < spaceDim; ++iDim)
- _cellVector[iBasis*spaceDim+iDim] -= valuesIJ[iBasis] *
- accArray[iBasis*spaceDim+iDim];
+ for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
+ for (int iDim = 0; iDim < spaceDim; ++iDim) {
+ _cellVector[iBasis*spaceDim+iDim] -= valuesIJ[iBasis] * accCell[iBasis*spaceDim+iDim];
+ } // for
+ } // for
PetscLogFlops(numQuadPts*(4+numBasis*3));
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- for (PetscInt i = 0; i < dispTSize; ++i) {
- dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];
+ for (PetscInt i = 0; i < dispSize; ++i) {
+ dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
} // for
- err = DMPlexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ accVisitor.restoreClosure(&accCell, &accSize, cell);
+ velVisitor.restoreClosure(&velCell, &velSize, cell);
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
- err = DMPlexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
-
// Compute B(transpose) * sigma, first computing strains
_calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispAdjCell, numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
@@ -339,37 +330,41 @@
CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispAdjCell);
// Assemble cell contribution into field
- err = DMPlexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
} // for
- err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_logger->eventEnd(computeEvent);
+
+ PYLITH_METHOD_END;
} // integrateResidual
// ----------------------------------------------------------------------
// Compute matrix associated with operator.
void
-pylith::feassemble::ElasticityExplicitLgDeform::integrateJacobian(
- topology::Jacobian* jacobian,
- const PylithScalar t,
- topology::SolutionFields* fields)
+pylith::feassemble::ElasticityExplicitLgDeform::integrateJacobian(topology::Jacobian* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* fields)
{ // integrateJacobian
+ PYLITH_METHOD_BEGIN;
+
throw std::logic_error("ElasticityExplicit::integrateJacobian() not implemented. Use integrateJacobian(lumped) instead.");
+
+ PYLITH_METHOD_END;
} // integrateJacobian
// ----------------------------------------------------------------------
// Compute matrix associated with operator.
void
-pylith::feassemble::ElasticityExplicitLgDeform::integrateJacobian(
- topology::Field<topology::Mesh>* jacobian,
- const PylithScalar t,
- topology::SolutionFields* fields)
+pylith::feassemble::ElasticityExplicitLgDeform::integrateJacobian(topology::Field<topology::Mesh>* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* fields)
{ // integrateJacobian
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != jacobian);
- assert(0 != fields);
+ PYLITH_METHOD_BEGIN;
+
+ assert(_quadrature);
+ assert(_material);
+ assert(jacobian);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIJ setup");
const int geometryEvent = _logger->eventId("ElIJ geometry");
@@ -394,32 +389,24 @@
"different dimensions than the spatial dimension.");
// Get cell information
- DM dmMesh = fields->mesh().dmMesh();
- IS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- PetscErrorCode err;
+ PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
+ topology::StratumIS materialIS(dmMesh, "material-id", _material->id());
+ const PetscInt* cells = materialIS.points();
+ const PetscInt numCells = materialIS.size();
- assert(dmMesh);
- err = DMPlexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
- err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
- err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-
// Get parameters used in integration.
const PylithScalar dt = _dt;
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
- // Get sections
- PetscSection jacSection = jacobian->petscSection();
- Vec jacVec = jacobian->localVector();
+ // Setup field visitors.
+ topology::VecVisitorMesh jacobianVisitor(*jacobian);
+ PetscScalar* jacobianCell = NULL;
+ PetscInt jacobianSize = 0;
- scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
@@ -431,12 +418,9 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(cell);
#else
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- _quadrature->computeGeometry(coordinatesCell, cell);
- err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
#endif
// Get state variables for cell.
@@ -472,15 +456,15 @@
_lumpCellMatrix();
// Assemble cell contribution into lumped matrix.
- err = DMPlexVecSetClosure(dmMesh, jacSection, jacVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ jacobianVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
} // for
- err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_needNewJacobian = false;
_material->resetNeedNewJacobian();
_logger->eventEnd(computeEvent);
+
+ PYLITH_METHOD_END;
} // integrateJacobian
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc 2013-04-04 21:47:27 UTC (rev 21723)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc 2013-04-04 23:35:43 UTC (rev 21724)
@@ -26,10 +26,12 @@
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
#include "pylith/utils/array.hh" // USES scalar_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
@@ -54,7 +56,6 @@
_dtm1(-1.0),
_normViscosity(0.1)
{ // constructor
- _basisDerivArray.resize(_numQuadPts*_numBasis*_spaceDim);
} // constructor
// ----------------------------------------------------------------------
@@ -77,14 +78,18 @@
void
pylith::feassemble::ElasticityExplicitTet4::timeStep(const PylithScalar dt)
{ // timeStep
+ PYLITH_METHOD_BEGIN;
+
if (_dt != -1.0)
_dtm1 = _dt;
else
_dtm1 = dt;
_dt = dt;
assert(_dt == _dtm1); // For now, don't allow variable time step
- if (0 != _material)
+ if (_material)
_material->timeStep(_dt);
+
+ PYLITH_METHOD_END;
} // timeStep
// ----------------------------------------------------------------------
@@ -92,8 +97,10 @@
PylithScalar
pylith::feassemble::ElasticityExplicitTet4::stableTimeStep(const topology::Mesh& mesh) const
{ // stableTimeStep
+ PYLITH_METHOD_BEGIN;
+
assert(_material);
- return _material->stableTimeStepExplicit(mesh, _quadrature);
+ PYLITH_METHOD_RETURN(_material->stableTimeStepExplicit(mesh, _quadrature));
} // stableTimeStep
// ----------------------------------------------------------------------
@@ -101,6 +108,8 @@
void
pylith::feassemble::ElasticityExplicitTet4::normViscosity(const PylithScalar viscosity)
{ // normViscosity
+ PYLITH_METHOD_BEGIN;
+
if (viscosity < 0.0) {
std::ostringstream msg;
msg << "Normalized viscosity (" << viscosity << ") must be nonnegative.";
@@ -108,6 +117,8 @@
} // if
_normViscosity = viscosity;
+
+ PYLITH_METHOD_END;
} // normViscosity
// ----------------------------------------------------------------------
@@ -118,14 +129,16 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateResidual
+ PYLITH_METHOD_BEGIN;
+
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityExplicitTet4::*elasticityResidual_fn_type)
(const scalar_array&);
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != _logger);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_material);
+ assert(_logger);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIR setup");
const int geometryEvent = _logger->eventId("ElIR geometry");
@@ -148,6 +161,7 @@
const int tensorSize = _tensorSize;
const int numBasis = _numBasis;
const int numQuadPts = _numQuadPts;
+ const int cellVectorSize = numBasis*spaceDim;
/** :TODO:
*
* If cellDim and spaceDim are different, we need to transform
@@ -167,54 +181,37 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- DM dmMesh = fields->mesh().dmMesh();
- IS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- PetscErrorCode err;
+ PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
+ topology::StratumIS materialIS(dmMesh, "material-id", _material->id());
+ const PetscInt* cells = materialIS.points();
+ const PetscInt numCells = materialIS.size();
- assert(dmMesh);
- err = DMPlexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
- err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
- err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ // Setup field visitors.
+ topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
+ PetscScalar* accCell = NULL;
+ PetscInt accSize = 0;
- // Get sections
- topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
- PetscSection accSection = acc.petscSection();
- Vec accVec = acc.localVector();
- assert(accSection);assert(accVec);
+ topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
+ PetscScalar* velCell = NULL;
+ PetscInt velSize = 0;
- topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
- PetscSection velSection = vel.petscSection();
- Vec velVec = vel.localVector();
- assert(velSection);assert(velVec);
-
scalar_array dispAdjCell(numBasis*spaceDim);
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- PetscSection dispTSection = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
- assert(dispTSection);assert(dispTVec);
+ topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
+ topology::VecVisitorMesh residualVisitor(residual);
- scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar *coordsCell = NULL;
+ PetscInt coordsSize = 0;
- assert(0 != _normalizer);
+ assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
- const PylithScalar gravityScale =
- _normalizer->pressureScale() / (_normalizer->lengthScale() *
- _normalizer->densityScale());
+ const PylithScalar gravityScale = _normalizer->pressureScale() / (_normalizer->lengthScale() * _normalizer->densityScale());
- const PylithScalar dt = _dt;
- assert(_normViscosity > 0.0);
- assert(dt > 0);
- const PylithScalar viscosity = dt*_normViscosity;
+ const PylithScalar dt = _dt;assert(dt > 0);
+ const PylithScalar viscosity = dt*_normViscosity;assert(_normViscosity >= 0.0);
// Get parameters used in integration.
scalar_array valuesIJ(numBasis);
@@ -232,13 +229,11 @@
#endif
// Restrict input fields to cell
- PetscScalar *accArray, *velArray, *dispTArray;
- PetscInt accSize, velSize, dispTSize;
- err = DMPlexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- assert(velSize == accSize);
- assert(dispTSize == accSize);
+ accVisitor.getClosure(&accCell, &accSize, cell);
+ velVisitor.getClosure(&velCell, &velSize, cell);
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);
+ assert(velSize == accSize);
+ assert(dispSize == accSize);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@ -246,13 +241,8 @@
#endif
// Compute geometry information for current cell
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- const PylithScalar volume = _volume(coordinatesCell);
- assert(volume > 0.0);
- err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ const PylithScalar volume = _volume(coordsCell, coordsSize);assert(volume > 0.0);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -274,36 +264,39 @@
_resetCellVector();
// Compute body force vector if gravity is being used.
- if (0 != _gravityField) {
- const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
- assert(0 != cs);
+ if (_gravityField) {
+ const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();assert(cs);
quadPtsGlobal = 0.0;
- for (int iBasis=0; iBasis < numBasis; ++iBasis)
- for (int iDim=0; iDim < spaceDim; ++iDim)
- quadPtsGlobal[iDim] +=
- coordinatesCell[iBasis*spaceDim+iDim] / numBasis;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ quadPtsGlobal[iDim] += coordsCell[iBasis*spaceDim+iDim] / numBasis;
+ } // for
+ } // for
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
// Compute action for element body forces
spatialdata::spatialdb::SpatialDB* db = _gravityField;
- const int err = db->query(&gravVec[0], gravVec.size(),
- &quadPtsGlobal[0], spaceDim, cs);
- if (err)
+ const int err = db->query(&gravVec[0], gravVec.size(), &quadPtsGlobal[0], spaceDim, cs);
+ if (err) {
throw std::runtime_error("Unable to get gravity vector for point.");
- _normalizer->nondimensionalize(&gravVec[0], gravVec.size(),
- gravityScale);
+ } // if
+ _normalizer->nondimensionalize(&gravVec[0], gravVec.size(), gravityScale);
const PylithScalar wtVertex = density[0] * volume / 4.0;
- for (int iBasis=0; iBasis < numBasis; ++iBasis)
- for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
_cellVector[iBasis * spaceDim + iDim] += wtVertex * gravVec[iDim];
+ } // for
+ } // for
PetscLogFlops(numBasis*spaceDim*2 + numBasis*spaceDim*2);
} // if
// Compute action for inertial terms
const PylithScalar wtVertex = density[0] * volume / 4.0;
- for(PetscInt i = 0; i < dispTSize; ++i) {_cellVector[i] -= wtVertex * accArray[i];}
+ assert(cellVectorSize == dispSize);
+ for(PetscInt i = 0; i < cellVectorSize; ++i) {
+ _cellVector[i] -= wtVertex * accCell[i];
+ } // for
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(2 + numBasis*spaceDim*2);
@@ -313,27 +306,29 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
- err = DMPlexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < cellVectorSize; ++i) {
+ dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
+ } // for
+ accVisitor.restoreClosure(&accCell, &accSize, cell);
+ velVisitor.restoreClosure(&velCell, &velSize, cell);
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
// Compute B(transpose) * sigma, first computing strains
- const PylithScalar x0 = coordinatesCell[0];
- const PylithScalar y0 = coordinatesCell[1];
- const PylithScalar z0 = coordinatesCell[2];
+ const PylithScalar x0 = coordsCell[0];
+ const PylithScalar y0 = coordsCell[1];
+ const PylithScalar z0 = coordsCell[2];
- const PylithScalar x1 = coordinatesCell[3];
- const PylithScalar y1 = coordinatesCell[4];
- const PylithScalar z1 = coordinatesCell[5];
+ const PylithScalar x1 = coordsCell[3];
+ const PylithScalar y1 = coordsCell[4];
+ const PylithScalar z1 = coordsCell[5];
- const PylithScalar x2 = coordinatesCell[6];
- const PylithScalar y2 = coordinatesCell[7];
- const PylithScalar z2 = coordinatesCell[8];
+ const PylithScalar x2 = coordsCell[6];
+ const PylithScalar y2 = coordsCell[7];
+ const PylithScalar z2 = coordsCell[8];
- const PylithScalar x3 = coordinatesCell[9];
- const PylithScalar y3 = coordinatesCell[10];
- const PylithScalar z3 = coordinatesCell[11];
+ const PylithScalar x3 = coordsCell[9];
+ const PylithScalar y3 = coordsCell[10];
+ const PylithScalar z3 = coordsCell[11];
const PylithScalar scaleB = 6.0 * volume;
const PylithScalar b1 = (y1*(z3-z2)-y2*z3+y3*z2-(y3-y2)*z1) / scaleB;
@@ -388,30 +383,18 @@
assert(_cellVector.size() == 12);
assert(stressCell.size() == 6);
- _cellVector[0] -=
- (d1*stressCell[5]+c1*stressCell[3]+b1*stressCell[0]) * volume;
- _cellVector[1] -=
- (d1*stressCell[4]+b1*stressCell[3]+c1*stressCell[1]) * volume;
- _cellVector[2] -=
- (b1*stressCell[5]+c1*stressCell[4]+d1*stressCell[2]) * volume;
- _cellVector[3] -=
- (d2*stressCell[5]+c2*stressCell[3]+b2*stressCell[0]) * volume;
- _cellVector[4] -=
- (d2*stressCell[4]+b2*stressCell[3]+c2*stressCell[1]) * volume;
- _cellVector[5] -=
- (b2*stressCell[5]+c2*stressCell[4]+d2*stressCell[2]) * volume;
- _cellVector[6] -=
- (d3*stressCell[5]+c3*stressCell[3]+b3*stressCell[0]) * volume;
- _cellVector[7] -=
- (d3*stressCell[4]+b3*stressCell[3]+c3*stressCell[1]) * volume;
- _cellVector[8] -=
- (b3*stressCell[5]+c3*stressCell[4]+d3*stressCell[2]) * volume;
- _cellVector[9] -=
- (d4*stressCell[5]+c4*stressCell[3]+b4*stressCell[0]) * volume;
- _cellVector[10] -=
- (d4*stressCell[4]+b4*stressCell[3]+c4*stressCell[1]) * volume;
- _cellVector[11] -=
- (b4*stressCell[5]+c4*stressCell[4]+d4*stressCell[2]) * volume;
+ _cellVector[ 0] -= (d1*stressCell[5]+c1*stressCell[3]+b1*stressCell[0]) * volume;
+ _cellVector[ 1] -= (d1*stressCell[4]+b1*stressCell[3]+c1*stressCell[1]) * volume;
+ _cellVector[ 2] -= (b1*stressCell[5]+c1*stressCell[4]+d1*stressCell[2]) * volume;
+ _cellVector[ 3] -= (d2*stressCell[5]+c2*stressCell[3]+b2*stressCell[0]) * volume;
+ _cellVector[ 4] -= (d2*stressCell[4]+b2*stressCell[3]+c2*stressCell[1]) * volume;
+ _cellVector[ 5] -= (b2*stressCell[5]+c2*stressCell[4]+d2*stressCell[2]) * volume;
+ _cellVector[ 6] -= (d3*stressCell[5]+c3*stressCell[3]+b3*stressCell[0]) * volume;
+ _cellVector[ 7] -= (d3*stressCell[4]+b3*stressCell[3]+c3*stressCell[1]) * volume;
+ _cellVector[ 8] -= (b3*stressCell[5]+c3*stressCell[4]+d3*stressCell[2]) * volume;
+ _cellVector[ 9] -= (d4*stressCell[5]+c4*stressCell[3]+b4*stressCell[0]) * volume;
+ _cellVector[10] -= (d4*stressCell[4]+b4*stressCell[3]+c4*stressCell[1]) * volume;
+ _cellVector[11] -= (b4*stressCell[5]+c4*stressCell[4]+d4*stressCell[2]) * volume;
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(84);
@@ -419,31 +402,36 @@
_logger->eventBegin(updateEvent);
#endif
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+
// Assemble cell contribution into field
- err = DMPlexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
- err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numCells*(2 + numBasis*spaceDim*2 + 196+84));
_logger->eventEnd(computeEvent);
#endif
+
+ PYLITH_METHOD_END;
} // integrateResidual
// ----------------------------------------------------------------------
// Compute matrix associated with operator.
void
-pylith::feassemble::ElasticityExplicitTet4::integrateJacobian(
- topology::Jacobian* jacobian,
- const PylithScalar t,
- topology::SolutionFields* fields)
+pylith::feassemble::ElasticityExplicitTet4::integrateJacobian(topology::Jacobian* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* fields)
{ // integrateJacobian
+ PYLITH_METHOD_BEGIN;
+
throw std::logic_error("ElasticityExplicit::integrateJacobian() not implemented. Use integrateJacobian(lumped) instead.");
+
+ PYLITH_METHOD_END;
} // integrateJacobian
// ----------------------------------------------------------------------
@@ -453,10 +441,12 @@
const PylithScalar t,
topology::SolutionFields* fields)
{ // integrateJacobian
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != jacobian);
- assert(0 != fields);
+ PYLITH_METHOD_BEGIN;
+
+ assert(_quadrature);
+ assert(_material);
+ assert(jacobian);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIJ setup");
const int geometryEvent = _logger->eventId("ElIJ geometry");
@@ -481,32 +471,24 @@
"different dimensions than the spatial dimension.");
// Get cell information
- DM dmMesh = fields->mesh().dmMesh();
- IS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- PetscErrorCode err;
+ PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
+ topology::StratumIS materialIS(dmMesh, "material-id", _material->id());
+ const PetscInt* cells = materialIS.points();
+ const PetscInt numCells = materialIS.size();
- assert(dmMesh);
- err = DMPlexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
- err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
- err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-
// Get parameters used in integration.
const PylithScalar dt = _dt;
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
- // Get sections
- PetscSection jacSection = jacobian->petscSection();
- Vec jacVec = jacobian->localVector();
+ // Setup visitors.
+ topology::VecVisitorMesh jacobianVisitor(*jacobian);
+ PetscScalar* jacobianCell = NULL;
+ PetscInt jacobianSize = 0;
- scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -519,13 +501,9 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- const PylithScalar volume = _volume(coordinatesCell);
- assert(volume > 0.0);
- err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ const PylithScalar volume = _volume(coordsCell, coordsSize);assert(volume > 0.0);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -542,7 +520,6 @@
// Compute Jacobian for inertial terms
const scalar_array& density = _material->calcDensity();
- assert(volume > 0.0);
_cellVector = density[0] * volume / (4.0 * dt2);
#if defined(DETAILED_EVENT_LOGGING)
@@ -552,14 +529,12 @@
#endif
// Assemble cell contribution into lumped matrix.
- err = DMPlexVecSetClosure(dmMesh, jacSection, jacVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ jacobianVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
- err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numCells*3);
@@ -568,6 +543,8 @@
_needNewJacobian = false;
_material->resetNeedNewJacobian();
+
+ PYLITH_METHOD_END;
} // integrateJacobian
// ----------------------------------------------------------------------
@@ -575,6 +552,8 @@
void
pylith::feassemble::ElasticityExplicitTet4::verifyConfiguration(const topology::Mesh& mesh) const
{ // verifyConfiguration
+ PYLITH_METHOD_BEGIN;
+
IntegratorElasticity::verifyConfiguration(mesh);
assert(_quadrature);
@@ -588,14 +567,17 @@
<< " # quad points: " << _numQuadPts << " (code), " << _quadrature->numQuadPts() << " (user)";
throw std::runtime_error(msg.str());
} // if
+
+ PYLITH_METHOD_END;
} // verifyConfiguration
// ----------------------------------------------------------------------
// Compute volume of tetrahedral cell.
PylithScalar
-pylith::feassemble::ElasticityExplicitTet4::_volume(const scalar_array& coordinatesCell) const
+pylith::feassemble::ElasticityExplicitTet4::_volume(const PylithScalar coordinatesCell[],
+ const int coordinatesSize) const
{ // __volume
- assert(12 == coordinatesCell.size());
+ assert(12 == coordinatesSize);
const PylithScalar x0 = coordinatesCell[0];
const PylithScalar y0 = coordinatesCell[1];
@@ -622,7 +604,7 @@
const PylithScalar volume = det / 6.0;
PetscLogFlops(48);
- return volume;
+ return volume;
} // _volume
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.hh 2013-04-04 21:47:27 UTC (rev 21723)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.hh 2013-04-04 23:35:43 UTC (rev 21724)
@@ -138,22 +138,15 @@
/** Compute volume of tetrahedral cell.
*
* @param coordinatesCell Coordinates of vertices of cell.
+ * @param coordinatesSize Size of array.
* @returns Volume of cell.
*/
- PylithScalar _volume(const scalar_array& coordinatesCell) const;
+ PylithScalar _volume(const PylithScalar coordinatesCell[],
+ const int coordinatesSize) const;
- /** Compute derivatives of basis functions of tetrahedral cell.
- *
- * @param coordinatesCell Coordinates of vertices of cell.
- * @returns Derivatives of basis functions.
- */
- const scalar_array& _basisDeriv(const scalar_array& coordinatesCell) const;
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
- scalar_array _basisDerivArray; ///< Array of basis derivatives
-
PylithScalar _dtm1; ///< Time step for t-dt1 -> t
PylithScalar _normViscosity; ///< Normalized viscosity for numerical damping.
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc 2013-04-04 21:47:27 UTC (rev 21723)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc 2013-04-04 23:35:43 UTC (rev 21724)
@@ -26,10 +26,12 @@
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
+#include "pylith/topology/Stratum.hh" // USES Stratum
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
#include "pylith/utils/array.hh" // USES scalar_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
@@ -54,7 +56,6 @@
_dtm1(-1.0),
_normViscosity(0.1)
{ // constructor
- _basisDerivArray.resize(_numQuadPts*_numBasis*_spaceDim);
} // constructor
// ----------------------------------------------------------------------
@@ -77,6 +78,8 @@
void
pylith::feassemble::ElasticityExplicitTri3::timeStep(const PylithScalar dt)
{ // timeStep
+ PYLITH_METHOD_BEGIN;
+
if (_dt != -1.0)
_dtm1 = _dt;
else
@@ -85,6 +88,8 @@
assert(_dt == _dtm1); // For now, don't allow variable time step
if (0 != _material)
_material->timeStep(_dt);
+
+ PYLITH_METHOD_END;
} // timeStep
// ----------------------------------------------------------------------
@@ -92,8 +97,10 @@
PylithScalar
pylith::feassemble::ElasticityExplicitTri3::stableTimeStep(const topology::Mesh& mesh) const
{ // stableTimeStep
+ PYLITH_METHOD_BEGIN;
+
assert(_material);
- return _material->stableTimeStepExplicit(mesh, _quadrature);
+ PYLITH_METHOD_RETURN(_material->stableTimeStepExplicit(mesh, _quadrature));
} // stableTimeStep
// ----------------------------------------------------------------------
@@ -101,6 +108,8 @@
void
pylith::feassemble::ElasticityExplicitTri3::normViscosity(const PylithScalar viscosity)
{ // normViscosity
+ PYLITH_METHOD_BEGIN;
+
if (viscosity < 0.0) {
std::ostringstream msg;
msg << "Normalized viscosity (" << viscosity << ") must be nonnegative.";
@@ -108,6 +117,8 @@
} // if
_normViscosity = viscosity;
+
+ PYLITH_METHOD_END;
} // normViscosity
// ----------------------------------------------------------------------
@@ -117,14 +128,16 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateResidual
+ PYLITH_METHOD_BEGIN;
+
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityExplicitTri3::*elasticityResidual_fn_type)
(const scalar_array&);
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != _logger);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_material);
+ assert(_logger);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIR setup");
const int geometryEvent = _logger->eventId("ElIR geometry");
@@ -147,6 +160,7 @@
const int tensorSize = _tensorSize;
const int numBasis = _numBasis;
const int numQuadPts = _numQuadPts;
+ const int cellVectorSize = numBasis*spaceDim;
/** :TODO:
*
* If cellDim and spaceDim are different, we need to transform
@@ -166,54 +180,37 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- DM dmMesh = fields->mesh().dmMesh();
- IS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- PetscErrorCode err;
+ PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
+ topology::StratumIS materialIS(dmMesh, "material-id", _material->id());
+ const PetscInt* cells = materialIS.points();
+ const PetscInt numCells = materialIS.size();
- assert(dmMesh);
- err = DMPlexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
- err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
- err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ // Setup field visitors.
+ topology::VecVisitorMesh accVisitor(fields->get("acceleration(t)"));
+ PetscScalar* accCell = NULL;
+ PetscInt accSize = 0;
- // Get sections
- topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
- PetscSection accSection = acc.petscSection();
- Vec accVec = acc.localVector();
- assert(accSection);assert(accVec);
+ topology::VecVisitorMesh velVisitor(fields->get("velocity(t)"));
+ PetscScalar* velCell = NULL;
+ PetscInt velSize = 0;
- topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
- PetscSection velSection = vel.petscSection();
- Vec velVec = vel.localVector();
- assert(velSection);assert(velVec);
-
scalar_array dispAdjCell(numBasis*spaceDim);
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- PetscSection dispTSection = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
- assert(dispTSection);assert(dispTVec);
+ topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
+ topology::VecVisitorMesh residualVisitor(residual);
- scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar *coordsCell = NULL;
+ PetscInt coordsSize = 0;
- assert(0 != _normalizer);
+ assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
- const PylithScalar gravityScale =
- _normalizer->pressureScale() / (_normalizer->lengthScale() *
- _normalizer->densityScale());
+ const PylithScalar gravityScale = _normalizer->pressureScale() / (_normalizer->lengthScale() * _normalizer->densityScale());
- const PylithScalar dt = _dt;
- assert(_normViscosity > 0.0);
- assert(dt > 0);
- const PylithScalar viscosity = dt*_normViscosity;
+ const PylithScalar dt = _dt;assert(dt > 0);
+ const PylithScalar viscosity = dt*_normViscosity;assert(_normViscosity >= 0.0);
// Get parameters used in integration.
scalar_array valuesIJ(numBasis);
@@ -231,13 +228,11 @@
#endif
// Restrict input fields to cell
- PetscScalar *accArray, *velArray, *dispTArray;
- PetscInt accSize, velSize, dispTSize;
- err = DMPlexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- assert(velSize == accSize);
- assert(dispTSize == accSize);
+ accVisitor.getClosure(&accCell, &accSize, cell);
+ velVisitor.getClosure(&velCell, &velSize, cell);
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);
+ assert(velSize == accSize);
+ assert(dispSize == accSize);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@ -245,13 +240,8 @@
#endif
// Compute geometry information for current cell
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- const PylithScalar area = _area(coordinatesCell);
- assert(area > 0.0);
- err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ const PylithScalar area = _area(coordsCell, coordsSize);assert(area > 0.0);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -273,36 +263,39 @@
_resetCellVector();
// Compute body force vector if gravity is being used.
- if (0 != _gravityField) {
- const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
- assert(0 != cs);
+ if (_gravityField) {
+ const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();assert(cs);
quadPtsGlobal = 0.0;
- for (int iBasis=0; iBasis < numBasis; ++iBasis)
- for (int iDim=0; iDim < spaceDim; ++iDim)
- quadPtsGlobal[iDim] +=
- coordinatesCell[iBasis*spaceDim+iDim] / numBasis;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ quadPtsGlobal[iDim] += coordsCell[iBasis*spaceDim+iDim] / numBasis;
+ } // for
+ } // for
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
// Compute action for element body forces
spatialdata::spatialdb::SpatialDB* db = _gravityField;
- const int err = db->query(&gravVec[0], gravVec.size(),
- &quadPtsGlobal[0], spaceDim, cs);
- if (err)
+ const int err = db->query(&gravVec[0], gravVec.size(), &quadPtsGlobal[0], spaceDim, cs);
+ if (err) {
throw std::runtime_error("Unable to get gravity vector for point.");
- _normalizer->nondimensionalize(&gravVec[0], gravVec.size(),
- gravityScale);
+ } // if
+ _normalizer->nondimensionalize(&gravVec[0], gravVec.size(), gravityScale);
const PylithScalar wtVertex = density[0] * area / 3.0;
- for (int iBasis=0; iBasis < numBasis; ++iBasis)
- for (int iDim=0; iDim < spaceDim; ++iDim)
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
_cellVector[iBasis * spaceDim + iDim] += wtVertex * gravVec[iDim];
+ } // for
+ } // for
PetscLogFlops(numBasis*spaceDim*2 + numBasis*spaceDim*2);
} // if
// Compute action for inertial terms
const PylithScalar wtVertex = density[0] * area / 3.0;
- for(PetscInt i = 0; i < dispTSize; ++i) {_cellVector[i] -= wtVertex * accArray[i];}
+ assert(cellVectorSize == dispSize);
+ for(PetscInt i = 0; i < cellVectorSize; ++i) {
+ _cellVector[i] -= wtVertex * accCell[i];
+ } // for
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(2 + numBasis*spaceDim*2);
@@ -312,20 +305,22 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
- err = DMPlexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < cellVectorSize; ++i) {
+ dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
+ } // for
+ accVisitor.restoreClosure(&accCell, &accSize, cell);
+ velVisitor.restoreClosure(&velCell, &velSize, cell);
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
// Compute B(transpose) * sigma, first computing strains
- const PylithScalar x0 = coordinatesCell[0];
- const PylithScalar y0 = coordinatesCell[1];
+ const PylithScalar x0 = coordsCell[0];
+ const PylithScalar y0 = coordsCell[1];
- const PylithScalar x1 = coordinatesCell[2];
- const PylithScalar y1 = coordinatesCell[3];
+ const PylithScalar x1 = coordsCell[2];
+ const PylithScalar y1 = coordsCell[3];
- const PylithScalar x2 = coordinatesCell[4];
- const PylithScalar y2 = coordinatesCell[5];
+ const PylithScalar x2 = coordsCell[4];
+ const PylithScalar y2 = coordsCell[5];
const PylithScalar scaleB = 2.0 * area;
const PylithScalar b0 = (y1 - y2) / scaleB;
@@ -368,20 +363,22 @@
_logger->eventBegin(updateEvent);
#endif
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+
// Assemble cell contribution into field
- err = DMPlexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
- err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numCells*(2 + numBasis*spaceDim*2 + 34+30));
_logger->eventEnd(computeEvent);
#endif
+
+ PYLITH_METHOD_END;
} // integrateResidual
// ----------------------------------------------------------------------
@@ -391,7 +388,11 @@
const PylithScalar t,
topology::SolutionFields* fields)
{ // integrateJacobian
+ PYLITH_METHOD_BEGIN;
+
throw std::logic_error("ElasticityExplicit::integrateJacobian() not implemented. Use integrateJacobian(lumped) instead.");
+
+ PYLITH_METHOD_END;
} // integrateJacobian
// ----------------------------------------------------------------------
@@ -402,11 +403,13 @@
const PylithScalar t,
topology::SolutionFields* fields)
{ // integrateJacobian
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != jacobian);
- assert(0 != fields);
+ PYLITH_METHOD_BEGIN;
+ assert(_quadrature);
+ assert(_material);
+ assert(jacobian);
+ assert(fields);
+
const int setupEvent = _logger->eventId("ElIJ setup");
const int geometryEvent = _logger->eventId("ElIJ geometry");
const int computeEvent = _logger->eventId("ElIJ compute");
@@ -430,32 +433,24 @@
"different dimensions than the spatial dimension.");
// Get cell information
- DM dmMesh = fields->mesh().dmMesh();
- IS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- PetscErrorCode err;
+ PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
+ topology::StratumIS materialIS(dmMesh, "material-id", _material->id());
+ const PetscInt* cells = materialIS.points();
+ const PetscInt numCells = materialIS.size();
- assert(dmMesh);
- err = DMPlexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
- err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
- err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
-
// Get parameters used in integration.
const PylithScalar dt = _dt;
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
- // Get sections
- PetscSection jacSection = jacobian->petscSection();
- Vec jacVec = jacobian->localVector();
+ // Setup visitors.
+ topology::VecVisitorMesh jacobianVisitor(*jacobian);
+ PetscScalar* jacobianCell = NULL;
+ PetscInt jacobianSize = 0;
- scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -468,13 +463,9 @@
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- const PylithScalar area = _area(coordinatesCell);
- assert(area > 0.0);
- err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);
+ const PylithScalar area = _area(coordsCell, coordsSize);assert(area > 0.0);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -500,14 +491,12 @@
#endif
// Assemble cell contribution into lumped matrix.
- err = DMPlexVecSetClosure(dmMesh, jacSection, jacVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ jacobianVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
- err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numCells*3);
@@ -516,6 +505,8 @@
_needNewJacobian = false;
_material->resetNeedNewJacobian();
+
+ PYLITH_METHOD_END;
} // integrateJacobian
// ----------------------------------------------------------------------
@@ -523,6 +514,8 @@
void
pylith::feassemble::ElasticityExplicitTri3::verifyConfiguration(const topology::Mesh& mesh) const
{ // verifyConfiguration
+ PYLITH_METHOD_BEGIN;
+
IntegratorElasticity::verifyConfiguration(mesh);
assert(_quadrature);
@@ -541,10 +534,10 @@
// ----------------------------------------------------------------------
// Compute area of triangular cell.
PylithScalar
-pylith::feassemble::ElasticityExplicitTri3::_area(
- const scalar_array& coordinatesCell) const
+pylith::feassemble::ElasticityExplicitTri3::_area(const PylithScalar coordinatesCell[],
+ const int coordinatesSize) const
{ // __area
- assert(6 == coordinatesCell.size());
+ assert(6 == coordinatesSize);
const PylithScalar x0 = coordinatesCell[0];
const PylithScalar y0 = coordinatesCell[1];
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.hh 2013-04-04 21:47:27 UTC (rev 21723)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.hh 2013-04-04 23:35:43 UTC (rev 21724)
@@ -138,22 +138,15 @@
/** Compute area of triangular cell.
*
* @param coordinatesCell Coordinates of vertices of cell.
- * @returns Area of cell.
+ * @param coordinatesSize Size of array.
+ * @returns Volume of cell.
*/
- PylithScalar _area(const scalar_array& coordinatesCell) const;
+ PylithScalar _area(const PylithScalar coordinatesCell[],
+ const int coordinatesSize) const;
- /** Compute derivatives of basis functions of triangular cell.
- *
- * @param coordinatesCell Coordinates of vertices of cell.
- * @returns Derivatives of basis functions.
- */
- const scalar_array& _basisDeriv(const scalar_array& coordinatesCell) const;
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
- scalar_array _basisDerivArray; ///< Array of basis derivatives
-
PylithScalar _dtm1; ///< Time step for t-dt1 -> t
PylithScalar _normViscosity; ///< Normalized viscosity for numerical damping.
More information about the CIG-COMMITS
mailing list