[cig-commits] r21725 - short/3D/PyLith/trunk/libsrc/pylith/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Apr 4 17:23:00 PDT 2013
Author: brad
Date: 2013-04-04 17:22:59 -0700 (Thu, 04 Apr 2013)
New Revision: 21725
Modified:
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
Log:
Code cleanup (feassemble implicit).
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2013-04-04 23:35:43 UTC (rev 21724)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2013-04-05 00:22:59 UTC (rev 21725)
@@ -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/EventLogger.hh" // USES EventLogger
#include "pylith/utils/array.hh" // USES scalar_array
@@ -45,9 +48,6 @@
//#define PRECOMPUTE_GEOMETRY
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityImplicit::ElasticityImplicit(void) :
_dtm1(-1.0)
@@ -74,13 +74,17 @@
void
pylith::feassemble::ElasticityImplicit::timeStep(const PylithScalar dt)
{ // timeStep
+ PYLITH_METHOD_BEGIN;
+
if (_dt != -1.0)
_dtm1 = _dt;
else
_dtm1 = dt;
_dt = dt;
- if (0 != _material)
+ if (_material)
_material->timeStep(_dt);
+
+ PYLITH_METHOD_END;
} // timeStep
// ----------------------------------------------------------------------
@@ -88,26 +92,28 @@
PylithScalar
pylith::feassemble::ElasticityImplicit::stableTimeStep(const topology::Mesh& mesh) const
{ // stableTimeStep
- assert(0 != _material);
- return _material->stableTimeStepImplicit(mesh);
+ PYLITH_METHOD_BEGIN;
+
+ assert(_material);
+ PYLITH_METHOD_RETURN(_material->stableTimeStepImplicit(mesh));
} // stableTimeStep
// ----------------------------------------------------------------------
void
-pylith::feassemble::ElasticityImplicit::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::feassemble::ElasticityImplicit::integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateResidual
+ PYLITH_METHOD_BEGIN;
+
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityImplicit::*elasticityResidual_fn_type)
(const scalar_array&);
- PetscErrorCode err;
- 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");
@@ -161,43 +167,29 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- DM dmMesh = fields->mesh().dmMesh();
- IS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- 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);
+ 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();
- // Get sections
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- PetscSection dispTSection = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
- assert(dispTSection);assert(dispTVec);
+ // Setup field visitors.
+ topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
- topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
- PetscSection dispTIncrSection = dispTIncr.petscSection();
- Vec dispTIncrVec = dispTIncr.localVector();
- assert(dispTIncrSection);assert(dispTIncrVec);
+ topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
+ PetscScalar* dispIncrCell = NULL;
+ PetscInt dispIncrSize = 0;
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
+ topology::VecVisitorMesh residualVisitor(residual);
-#if !defined(PRECOMPUTE_GEOMETRY)
- 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);
-#endif
+ 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());
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
@@ -209,12 +201,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.
@@ -224,11 +213,9 @@
_resetCellVector();
// Restrict input fields to cell
- PetscScalar *dispTArray, *dispTIncrArray;
- PetscInt dispTSize, dispTIncrSize;
- err = DMPlexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
- assert(dispTSize == dispTIncrSize);
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);
+ dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);
+ assert(dispSize == dispIncrSize);
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -237,14 +224,15 @@
const scalar_array& quadPtsNondim = _quadrature->quadPts();
// Compute current estimate of displacement at time t+dt using solution increment.
- for(PetscInt i = 0; i < dispTSize; ++i) {dispTpdtCell[i] = dispTArray[i] + dispTIncrArray[i];}
- err = DMPlexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < dispSize; ++i) {
+ dispTpdtCell[i] = dispCell[i] + dispIncrCell[i];
+ } // for
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+ dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// 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();
@@ -256,7 +244,9 @@
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) throw std::runtime_error("Unable to get gravity vector for point.");
+ if (err) {
+ throw std::runtime_error("Unable to get gravity vector for point.");
+ } // 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) {
@@ -267,7 +257,7 @@
} // for
} // for
PetscLogFlops(numQuadPts * (2 + numBasis * (1 + 2 * spaceDim)));
- }
+ } // if
// residualSection->view("After gravity contribution");
// Compute B(transpose) * sigma, first computing strains
@@ -280,34 +270,35 @@
std::cout << "Updating residual for cell " << cell << std::endl;
for(PetscInt i = 0; i < spaceDim*numBasis; ++i) {
std::cout << " v["<<i<<"]: " << _cellVector[i] << std::endl;
- }
+ } // for
#endif
// Assemble cell contribution into field
- err = DMPlexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
- }
- err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+ residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
+ } // for
+
_logger->eventEnd(computeEvent);
+
+ PYLITH_METHOD_END;
} // integrateResidual
// ----------------------------------------------------------------------
// Compute stiffness matrix.
void
-pylith::feassemble::ElasticityImplicit::integrateJacobian(
- topology::Jacobian* jacobian,
- const PylithScalar t,
- topology::SolutionFields* fields)
+pylith::feassemble::ElasticityImplicit::integrateJacobian(topology::Jacobian* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* fields)
{ // integrateJacobian
+ PYLITH_METHOD_BEGIN;
+
/// Member prototype for _elasticityJacobianXD()
typedef void (pylith::feassemble::ElasticityImplicit::*elasticityJacobian_fn_type)
(const scalar_array&);
- PetscErrorCode err;
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != _logger);
- assert(0 != jacobian);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_material);
+ assert(_logger);
+ assert(jacobian);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIJ setup");
const int geometryEvent = _logger->eventId("ElIJ geometry");
@@ -358,42 +349,32 @@
strainCell = 0.0;
// Get cell information
- DM dmMesh = fields->mesh().dmMesh();
- IS cellIS;
- const PetscInt *cells;
- PetscInt numCells;
- assert(dmMesh);
- const int materialId = _material->id();
- err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
- err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
- err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(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();
- // Get sections
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- PetscSection dispTSection = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
- assert(dispTSection);assert(dispTVec);
+ // Setup field visitors.
+ topology::VecVisitorMesh dispVisitor(fields->get("disp(t)"));
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
- topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
- PetscSection dispTIncrSection = dispTIncr.petscSection();
- Vec dispTIncrVec = dispTIncr.localVector();
- assert(dispTIncrSection);assert(dispTIncrVec);
+ topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
+ PetscScalar* dispIncrCell = NULL;
+ PetscInt dispIncrSize = 0;
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
+
// Get sparse matrix
- const PetscMat jacobianMat = jacobian->matrix();
- assert(0 != jacobianMat);
+ const PetscMat jacobianMat = jacobian->matrix();assert(jacobianMat);
+ topology::MatVisitorMesh jacobianVisitor(jacobianMat, fields->get("disp(t)"));
// Get parameters used in integration.
const PylithScalar dt = _dt;
assert(dt > 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);
-
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
@@ -404,12 +385,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 physical properties and state variables for cell.
@@ -419,11 +397,9 @@
_resetCellMatrix();
// Restrict input fields to cell
- PetscScalar *dispTArray, *dispTIncrArray;
- PetscInt dispTSize, dispTIncrSize;
- err = DMPlexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
- assert(dispTSize == dispTIncrSize);
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);
+ dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);
+ assert(dispSize == dispIncrSize);
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -431,9 +407,11 @@
const scalar_array& jacobianDet = _quadrature->jacobianDet();
// Compute current estimate of displacement at time t+dt using solution increment.
- for(PetscInt i = 0; i < dispTSize; ++i) {dispTpdtCell[i] = dispTArray[i] + dispTIncrArray[i];}
- err = DMPlexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < dispSize; ++i) {
+ dispTpdtCell[i] = dispCell[i] + dispIncrCell[i];
+ } // for
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+ dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// Compute strains
calcTotalStrainFn(&strainCell, basisDeriv, dispTpdtCell, numBasis, numQuadPts);
@@ -476,15 +454,15 @@
// Assemble cell contribution into PETSc matrix.
// Notice that we are using the default sections
- err = DMPlexMatSetClosure(dmMesh, dispTSection, PETSC_NULL, jacobianMat, cell, &_cellMatrix[0], ADD_VALUES);
- CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
+ jacobianVisitor.setClosure(&_cellMatrix[0], _cellMatrix.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/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2013-04-04 23:35:43 UTC (rev 21724)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2013-04-05 00:22:59 UTC (rev 21725)
@@ -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/EventLogger.hh" // USES EventLogger
#include "pylith/utils/array.hh" // USES scalar_array
@@ -71,12 +74,14 @@
void
pylith::feassemble::ElasticityImplicitLgDeform::timeStep(const PylithScalar dt)
{ // timeStep
+ PYLITH_METHOD_BEGIN;
+
if (_dt != -1.0)
_dtm1 = _dt;
else
_dtm1 = dt;
_dt = dt;
- if (0 != _material)
+ if (_material)
_material->timeStep(_dt);
} // timeStep
@@ -85,8 +90,10 @@
PylithScalar
pylith::feassemble::ElasticityImplicitLgDeform::stableTimeStep(const topology::Mesh& mesh) const
{ // stableTimeStep
- assert(0 != _material);
- return _material->stableTimeStepImplicit(mesh);
+ PYLITH_METHOD_BEGIN;
+
+ assert(_material);
+ PYLITH_METHOD_RETURN(_material->stableTimeStepImplicit(mesh));
} // stableTimeStep
// ----------------------------------------------------------------------
@@ -97,14 +104,16 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateResidual
+ PYLITH_METHOD_BEGIN;
+
/// Member prototype for _elasticityResidualXD()
typedef void (pylith::feassemble::ElasticityImplicitLgDeform::*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");
@@ -160,39 +169,28 @@
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 dispVisitor(fields->get("disp(t)"));
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
- // Get sections
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- PetscSection dispTSection = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
- assert(dispTSection);assert(dispTVec);
+ topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
+ PetscScalar* dispIncrCell = NULL;
+ PetscInt dispIncrSize = 0;
- topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
- PetscSection dispTIncrSection = dispTIncr.petscSection();
- Vec dispTIncrVec = dispTIncr.localVector();
- assert(dispTIncrSection);assert(dispTIncrVec);
+ 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() *
@@ -208,12 +206,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);
+ _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsCell[i];} // :TODO: Remove copy
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
#endif
// Get state variables for cell.
@@ -223,11 +219,9 @@
_resetCellVector();
// Restrict input fields to cell
- PetscScalar *dispTArray, *dispTIncrArray;
- PetscInt dispTSize, dispTIncrSize;
- err = DMPlexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
- assert(dispTSize == dispTIncrSize);
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);
+ dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);
+ assert(dispSize == dispIncrSize);
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -237,34 +231,32 @@
// Compute current estimate of displacement at time t+dt using
// solution increment.
- for(PetscInt i = 0; i < dispTSize; ++i) {dispTpdtCell[i] = dispTArray[i] + dispTIncrArray[i];}
- err = DMPlexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < dispSize; ++i) {
+ dispTpdtCell[i] = dispCell[i] + dispIncrCell[i];
+ } // for
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+ dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// 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];
@@ -276,8 +268,7 @@
// Compute B(transpose) * sigma, first computing deformation
// tensor and strains
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTpdtCell,
- numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTpdtCell, numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, true);
@@ -290,31 +281,32 @@
}
#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);
} // for
- err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_logger->eventEnd(computeEvent);
+
+ PYLITH_METHOD_END;
} // integrateResidual
// ----------------------------------------------------------------------
// Compute stiffness matrix.
void
-pylith::feassemble::ElasticityImplicitLgDeform::integrateJacobian(
- topology::Jacobian* jacobian,
- const PylithScalar t,
- topology::SolutionFields* fields)
+pylith::feassemble::ElasticityImplicitLgDeform::integrateJacobian(topology::Jacobian* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* fields)
{ // integrateJacobian
+ PYLITH_METHOD_BEGIN;
+
/// Member prototype for _elasticityJacobianXD()
typedef void (pylith::feassemble::ElasticityImplicitLgDeform::*elasticityJacobian_fn_type)
(const scalar_array&, const scalar_array&, const scalar_array&);
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != _logger);
- assert(0 != jacobian);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_material);
+ assert(_logger);
+ assert(jacobian);
+ assert(fields);
const int setupEvent = _logger->eventId("ElIJ setup");
const int geometryEvent = _logger->eventId("ElIJ geometry");
@@ -368,44 +360,33 @@
scalar_array stressCell(numQuadPts*tensorSize);
// 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);
- const int materialId = _material->id();
- err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &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 dispVisitor(fields->get("disp(t)"));
+ PetscScalar* dispCell = NULL;
+ PetscInt dispSize = 0;
- // Get sections
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- PetscSection dispTSection = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
- assert(dispTSection);assert(dispTVec);
+ topology::VecVisitorMesh dispIncrVisitor(fields->get("dispIncr(t->t+dt)"));
+ PetscScalar* dispIncrCell = NULL;
+ PetscInt dispIncrSize = 0;
- topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
- PetscSection dispTIncrSection = dispTIncr.petscSection();
- Vec dispTIncrVec = dispTIncr.localVector();
- assert(dispTIncrSection);assert(dispTIncrVec);
+ scalar_array coordinatesCell(numBasis*spaceDim);
+ topology::CoordsVisitor coordsVisitor(dmMesh);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
// Get sparse matrix
- const PetscMat jacobianMat = jacobian->matrix();
- assert(0 != jacobianMat);
+ const PetscMat jacobianMat = jacobian->matrix();assert(jacobianMat);
+ topology::MatVisitorMesh jacobianVisitor(jacobianMat, fields->get("disp(t)"));
// Get parameters used in integration.
const PylithScalar dt = _dt;
assert(dt > 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);
-
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
@@ -416,12 +397,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);
+ _quadrature->computeGeometry(coordsCell, coordsSize, cell);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsCell[i];} // :TODO: Remove copy.
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
#endif
// Get state variables for cell.
@@ -431,11 +410,9 @@
_resetCellMatrix();
// Restrict input fields to cell
- PetscScalar *dispTArray, *dispTIncrArray;
- PetscInt dispTSize, dispTIncrSize;
- err = DMPlexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecGetClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
- assert(dispTSize == dispTIncrSize);
+ dispVisitor.getClosure(&dispCell, &dispSize, cell);
+ dispIncrVisitor.getClosure(&dispIncrCell, &dispIncrSize, cell);
+ assert(dispSize == dispIncrSize);
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -444,13 +421,14 @@
// Compute current estimate of displacement at time t+dt using
// solution increment.
- for(PetscInt i = 0; i < dispTSize; ++i) {dispTpdtCell[i] = dispTArray[i] + dispTIncrArray[i];}
- err = DMPlexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecRestoreClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < dispSize; ++i) {
+ dispTpdtCell[i] = dispCell[i] + dispIncrCell[i];
+ } // for
+ dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+ dispIncrVisitor.restoreClosure(&dispIncrCell, &dispIncrSize, cell);
// Compute deformation tensor, strains, and stresses
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTpdtCell,
- numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTpdtCell, numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
// Get "elasticity" matrix at quadrature points for this cell
@@ -460,8 +438,7 @@
// Get Second Priola-Kirchoff stress tensor
const scalar_array& stressCell = _material->calcStress(strainCell, true);
- CALL_MEMBER_FN(*this, elasticityJacobianFn)(elasticConsts, stressCell,
- dispTpdtCell);
+ CALL_MEMBER_FN(*this, elasticityJacobianFn)(elasticConsts, stressCell, dispTpdtCell);
if (_quadrature->checkConditioning()) {
int n = numBasis*spaceDim;
@@ -495,14 +472,15 @@
} // if
// Assemble cell contribution into PETSc matrix.
- err = DMPlexMatSetClosure(dmMesh, dispTSection, PETSC_NULL, jacobianMat, cell, &_cellMatrix[0], ADD_VALUES);
- CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
+ jacobianVisitor.setClosure(&_cellMatrix[0], _cellMatrix.size(), cell, ADD_VALUES);
} // for
_logger->eventEnd(computeEvent);
_needNewJacobian = false;
_material->resetNeedNewJacobian();
+
+ PYLITH_METHOD_END;
} // integrateJacobian
More information about the CIG-COMMITS
mailing list