[cig-commits] r21721 - short/3D/PyLith/trunk/libsrc/pylith/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Apr 4 13:01:52 PDT 2013
Author: brad
Date: 2013-04-04 13:01:51 -0700 (Thu, 04 Apr 2013)
New Revision: 21721
Modified:
short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
Log:
Code cleanup (feassemble).
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.cc 2013-04-04 18:29:46 UTC (rev 21720)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/CellGeometry.cc 2013-04-04 20:01:51 UTC (rev 21721)
@@ -176,7 +176,7 @@
const PylithScalar jacobianDet,
const scalar_array& upDir)
{ // _orient0D
- assert(0 != orientation);
+ assert(orientation);
assert(1 == orientation->size());
(*orientation) = 1.0;
} // _orient0D
@@ -190,7 +190,7 @@
const scalar_array& upDir)
{ // _orient1D
const int orientSize = 4;
- assert(0 != orientation);
+ assert(orientation);
assert(orientSize == orientation->size());
const int jacobianSize = 2;
assert(jacobianSize == jacobian.size());
@@ -216,9 +216,9 @@
const scalar_array& upDir)
{ // _orient2D
const int orientSize = 9;
- assert(0 != orientation);
+ const int jacobianSize = 6;
+ assert(orientation);
assert(orientSize == orientation->size());
- const int jacobianSize = 6;
assert(jacobianSize == jacobian.size());
assert(3 == upDir.size());
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2013-04-04 18:29:46 UTC (rev 21720)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2013-04-04 20:01:51 UTC (rev 21721)
@@ -27,6 +27,9 @@
#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
@@ -72,14 +75,20 @@
void
pylith::feassemble::ElasticityExplicit::timeStep(const PylithScalar dt)
{ // timeStep
- if (_dt != -1.0)
+ PYLITH_METHOD_BEGIN;
+
+ if (_dt != -1.0) {
_dtm1 = _dt;
- else
+ } else {
_dtm1 = dt;
+ } // if/else
_dt = dt;
assert(_dt == _dtm1); // For now, don't allow variable time step
- if (0 != _material)
+ if (_material) {
_material->timeStep(_dt);
+ } // if
+
+ PYLITH_METHOD_END;
} // timeStep
// ----------------------------------------------------------------------
@@ -87,8 +96,10 @@
PylithScalar
pylith::feassemble::ElasticityExplicit::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
// ----------------------------------------------------------------------
@@ -96,6 +107,8 @@
void
pylith::feassemble::ElasticityExplicit::normViscosity(const PylithScalar viscosity)
{ // normViscosity
+ PYLITH_METHOD_BEGIN;
+
if (viscosity < 0.0) {
std::ostringstream msg;
msg << "Normalized viscosity (" << viscosity << ") must be nonnegative.";
@@ -103,6 +116,8 @@
} // if
_normViscosity = viscosity;
+
+ PYLITH_METHOD_END;
} // normViscosity
// ----------------------------------------------------------------------
@@ -112,14 +127,16 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateResidual
+ PYLITH_METHOD_BEGIN;
+
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityExplicit::*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");
@@ -181,54 +198,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);
@@ -248,12 +248,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
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -272,13 +269,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);
- 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);
@@ -292,26 +287,23 @@
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) {
const PylithScalar valI = wt * basis[iQ + iBasis];
@@ -330,33 +322,37 @@
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
+ // Numerical damping. Compute displacements adjusted by velocity
+ // times normalized viscosity.
+ for(PetscInt i = 0; i < dispSize; ++i) {
+ dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
+ } // for
+ accVisitor.restoreClosure(&accCell, &accSize, cell);
+ velVisitor.restoreClosure(&velCell, &velSize, cell);
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(4+numBasis*3));
_logger->eventEnd(computeEvent);
_logger->eventBegin(stressEvent);
#endif
- // 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);
-
// Compute B(transpose) * sigma, first computing strains
calcTotalStrainFn(&strainCell, basisDeriv, dispAdjCell, numBasis, numQuadPts);
-
const scalar_array& stressCell = _material->calcStress(strainCell, false);
#if defined(DETAILED_EVENT_LOGGING)
@@ -372,44 +368,48 @@
#endif
// 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*numQuadPts*(4+numBasis*3));
_logger->eventEnd(computeEvent);
#endif
+
+ PYLITH_METHOD_END;
} // integrateResidualLumped
// ----------------------------------------------------------------------
// Compute matrix associated with operator.
void
-pylith::feassemble::ElasticityExplicit::integrateJacobian(
- topology::Jacobian* jacobian,
- const PylithScalar t,
- topology::SolutionFields* fields)
+pylith::feassemble::ElasticityExplicit::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::ElasticityExplicit::integrateJacobian(
- topology::Field<topology::Mesh>* jacobian,
- const PylithScalar t,
- topology::SolutionFields* fields)
+pylith::feassemble::ElasticityExplicit::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");
@@ -433,33 +433,26 @@
"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);
+
+ // Setup field visitors.
scalar_array valuesIJ(numBasis);
- // Get sections
- PetscSection jacSection = jacobian->petscSection();
- Vec jacVec = jacobian->localVector();
+ 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,6 +461,7 @@
// Loop over cells
for(PetscInt c = 0; c < numCells; ++c) {
const PetscInt cell = cells[c];
+
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
@@ -475,12 +469,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
#if defined(DETAILED_EVENT_LOGGING)
@@ -509,20 +500,22 @@
// Compute Jacobian for inertial terms
valuesIJ = 0.0;
for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
- const PylithScalar wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad]
- / dt2;
+ const PylithScalar wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
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) {
valuesIJ[iBasis] += basis[iQ + iBasis] * valJ;
} // for
} // for
- 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] += valuesIJ[iBasis];
+ } // for
+ } // for
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(4 + numBasis*3) + numBasis*spaceDim);
@@ -531,14 +524,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*(numQuadPts*(4 + numBasis*3) + numBasis*spaceDim));
@@ -547,6 +538,8 @@
_needNewJacobian = false;
_material->resetNeedNewJacobian();
+
+ PYLITH_METHOD_END;
} // integrateJacobian
More information about the CIG-COMMITS
mailing list