[cig-commits] r21527 - in short/3D/PyLith/trunk: libsrc/pylith/bc libsrc/pylith/faults libsrc/pylith/topology unittests/libtests/bc
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 13 13:29:41 PDT 2013
Author: brad
Date: 2013-03-13 13:29:41 -0700 (Wed, 13 Mar 2013)
New Revision: 21527
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.hh
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.icc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
short/3D/PyLith/trunk/libsrc/pylith/topology/Makefile.am
short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc
short/3D/PyLith/trunk/libsrc/pylith/topology/Stratum.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.icc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
Log:
Small fixes to visitors. More use of visitors.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -25,17 +25,21 @@
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
#include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/VisitorSubMesh.hh" // USES VecVisitorSubMesh, MatVisitorSubMesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
#include <cstring> // USES memcpy()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
+//#define PRECOMPUTE_GEOMETRY
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
@@ -82,12 +86,11 @@
const int numCorners = _quadrature->numBasis();
// Get 'surface' cells (1 dimension lower than top-level cells)
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
+ const PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
- PetscInt cStart, cEnd;
- PetscErrorCode err;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
-
// Get damping constants at each quadrature point and rotate to
// global coordinate frame using orientation information
const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
@@ -135,13 +138,12 @@
// Container for damping constants for current cell
scalar_array dampingConstsLocal(spaceDim);
topology::Field<topology::SubMesh>& dampingConsts = _parameters->get("damping constants");
+ topology::VecVisitorMesh dampingConstsVisitor(dampingConsts);
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ PetscScalar *coordsCell = NULL;
+ PetscInt coordsSize = 0;
+ topology::CoordsVisitor coordsVisitor(dmSubMesh);
assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -154,20 +156,19 @@
// Compute quadrature information
_quadrature->initializeGeometry();
- PetscScalar* dampingConstsArray = dampingConsts.getLocalArray();
+ PetscScalar* dampingConstsArray = dampingConstsVisitor.localArray();
+
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for (PetscInt i=0; i < coordsSize; ++i) {
- coordinatesCell[i] = coords[i];
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+ for (PetscInt i=0; i < coordsSize; ++i) { // :TODO: Remove copy.
+ coordinatesCell[i] = coordsCell[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
- const PetscInt doff = dampingConsts.sectionOffset(c);
- assert(fiberDim == dampingConsts.sectionDof(c));
+ const PetscInt doff = dampingConstsVisitor.sectionOffset(c);
+ assert(fiberDim == dampingConstsVisitor.sectionDof(c));
const scalar_array& quadPtsNondim = _quadrature->quadPts();
const scalar_array& quadPtsRef = _quadrature->quadPtsRef();
@@ -217,7 +218,6 @@
} // for
} // for
} // for
- dampingConsts.restoreLocalArray(&dampingConstsArray);
_db->close();
} // initialize
@@ -253,47 +253,42 @@
// Allocate vectors for cell values.
_initCellVector();
-
- // Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
- PetscIS subpointIS = NULL;
- PetscInt cStart, cEnd;
- PetscErrorCode err;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
-
// Get sections
topology::Field<topology::SubMesh>& dampingConsts = _parameters->get("damping constants");
+ topology::VecVisitorMesh dampingConstsVisitor(dampingConsts);
+ PetscScalar* dampingConstsArray = dampingConstsVisitor.localArray();
+ // Get subsections
+ const PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::SubMeshIS submeshIS(*_boundaryMesh);
+
// Use _cellVector for cell residual.
- PetscVec residualVec = residual.localVector();assert(residualVec);
- PetscSection residualSection = residual.petscSection();assert(residualSection);
- PetscSection residualSubsection = NULL;
- err = PetscSectionCreateSubmeshSection(residualSection, subpointIS, &residualSubsection);
+ topology::VecVisitorSubMesh residualVisitor(residual, submeshIS);
- topology::Field<topology::Mesh>& velocity = fields->get("velocity(t)");
- PetscSection velSection = velocity.petscSection();assert(velSection);
- PetscVec velVec = velocity.localVector();assert(velVec);
- PetscSection velSubsection = NULL;
- err = PetscSectionCreateSubmeshSection(velSection, subpointIS, &velSubsection);
- err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
+ PetscScalar *velocityArray = NULL;
+ PetscInt velocitySize = 0;
+ topology::VecVisitorSubMesh velocityVisitor(fields->get("velocity(t)"), submeshIS);
+
+ submeshIS.deallocate();
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);assert(coordSection);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);assert(coordVec);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
+topology::CoordsVisitor coordsVisitor(dmSubMesh);
#endif
+ // Get 'surface' cells (1 dimension lower than top-level cells)
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
+
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
#endif
- PetscScalar* dampingConstsArray = dampingConsts.getLocalArray();
-
- for(PetscInt c = cStart; c < cEnd; ++c) {
+ for (PetscInt c = cStart; c < cEnd; ++c) {
// Get geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
@@ -301,14 +296,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
- PetscScalar *coordsArray;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
- for (PetscInt i=0; i < coordsSize; ++i) {
- coordinatesCell[i] = coordsArray[i];
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+ for (PetscInt i=0; i < coordsSize; ++i) { // :TODO: Remove copy.
+ coordinatesCell[i] = coordsCell[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -320,13 +313,11 @@
_resetCellVector();
// Restrict input fields to cell
- PetscScalar *velArray = NULL;
- PetscInt velSize;
- err = DMPlexVecGetClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- assert(velSize == numBasis*spaceDim);
+ velocityVisitor.getClosure(&velocityArray, &velocitySize, c);
+ assert(velocitySize == numBasis*spaceDim);
- const PetscInt doff = dampingConsts.sectionOffset(c);
- assert(numQuadPts*spaceDim == dampingConsts.sectionDof(c));
+ const PetscInt doff = dampingConstsVisitor.sectionOffset(c);
+ assert(numQuadPts*spaceDim == dampingConstsVisitor.sectionDof(c));
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@ -348,13 +339,14 @@
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -=
dampingConstsArray[doff+iQuad*spaceDim+iDim] *
- valIJ * velArray[jBasis*spaceDim+iDim];
+ valIJ * velocityArray[jBasis*spaceDim+iDim];
} // for
} // for
} // for
- err = DMPlexVecRestoreClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecSetClosure(subMesh, residualSubsection, residualVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ velocityVisitor.restoreClosure(&velocityArray, &velocitySize, c);
+ residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), c, ADD_VALUES);
+
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
_logger->eventEnd(computeEvent);
@@ -365,9 +357,6 @@
_logger->eventEnd(updateEvent);
#endif
} // for
- dampingConsts.restoreLocalArray(&dampingConstsArray);
- err = PetscSectionDestroy(&residualSubsection);CHECK_PETSC_ERROR(err);
- err = PetscSectionDestroy(&velSubsection);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops((cEnd-cStart)*numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
@@ -406,36 +395,33 @@
// Allocate vectors for cell values.
_initCellVector();
- // Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
- PetscIS subpointIS = NULL;
- PetscInt cStart, cEnd;
- PetscErrorCode err = 0;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
+ // Get 'surface' cells (1 dimension lower than top-level cells)
+ const PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
// Get sections
topology::Field<topology::SubMesh>& dampingConsts = _parameters->get("damping constants");
+ topology::VecVisitorMesh dampingConstsVisitor(dampingConsts);
+ PetscScalar* dampingConstsArray = dampingConstsVisitor.localArray();
+ // Get subsections
+ topology::SubMeshIS submeshIS(*_boundaryMesh);
// Use _cellVector for cell values.
- PetscSection residualSection = residual.petscSection();assert(residualSection);
- PetscVec residualVec = residual.localVector();assert(residualVec);
- PetscSection residualSubsection = NULL;
- err = PetscSectionCreateSubmeshSection(residualSection, subpointIS, &residualSubsection);
+ topology::VecVisitorSubMesh residualVisitor(residual, submeshIS);
- PetscSection velSection = fields->get("velocity(t)").petscSection();assert(velSection);
- PetscVec velVec = fields->get("velocity(t)").localVector();assert(velVec);
- PetscSection velSubsection = NULL;
- err = PetscSectionCreateSubmeshSection(velSection, subpointIS, &velSubsection);
- err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
+ PetscScalar *velocityArray = NULL;
+ PetscInt velocitySize;
+ topology::VecVisitorSubMesh velocityVisitor(fields->get("velocity(t)"), submeshIS);
+ submeshIS.deallocate();
+
#if !defined(PRECOMPUTE_GEOMETRY)
+ PetscScalar *coordsCell = NULL;
+ PetscInt coordsSize = 0;
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmSubMesh);
#endif
_logger->eventEnd(setupEvent);
@@ -443,8 +429,6 @@
_logger->eventBegin(computeEvent);
#endif
- PetscScalar* dampingConstsArray = dampingConsts.getLocalArray();
-
for (PetscInt c=cStart; c < cEnd; ++c) {
// Get geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
@@ -453,14 +437,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
- PetscScalar *coordsArray;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
- for (PetscInt i=0; i < coordsSize; ++i) {
- coordinatesCell[i] = coordsArray[i];
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+ for (PetscInt i=0; i < coordsSize; ++i) { // :TODO: Remove copy.
+ coordinatesCell[i] = coordsCell[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -472,10 +454,8 @@
_resetCellVector();
// Restrict input fields to cell
- PetscScalar *velArray = NULL;
- PetscInt velSize;
- err = DMPlexVecGetClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- assert(velSize == numBasis*spaceDim);
+ velocityVisitor.getClosure(&velocityArray, &velocitySize, c);
+ assert(velocitySize == numBasis*spaceDim);
const PetscInt doff = dampingConsts.sectionOffset(c);
assert(numQuadPts*spaceDim == dampingConsts.sectionDof(c));
@@ -502,11 +482,13 @@
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -=
dampingConstsArray[doff+iQuad*spaceDim+iDim] *
- valIJ * velArray[iBasis*spaceDim+iDim];
+ valIJ * velocityArray[iBasis*spaceDim+iDim];
} // for
} // for
- err = DMPlexVecRestoreClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
- err = DMPlexVecSetClosure(subMesh, residualSubsection, residualVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ velocityVisitor.restoreClosure(&velocityArray, &velocitySize, c);
+
+ residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), c, ADD_VALUES);
+
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3)));
_logger->eventEnd(computeEvent);
@@ -517,9 +499,6 @@
_logger->eventEnd(updateEvent);
#endif
} // for
- dampingConsts.restoreLocalArray(&dampingConstsArray);
- err = PetscSectionDestroy(&residualSubsection);CHECK_PETSC_ERROR(err);
- err = PetscSectionDestroy(&velSubsection);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops((cEnd-cStart)*numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3)));
@@ -530,10 +509,9 @@
// ----------------------------------------------------------------------
// Integrate contributions to Jacobian matrix (A) associated with
void
-pylith::bc::AbsorbingDampers::integrateJacobian(
- topology::Jacobian* jacobian,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::bc::AbsorbingDampers::integrateJacobian(topology::Jacobian* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateJacobian
assert(_quadrature);
assert(_boundaryMesh);
@@ -556,31 +534,23 @@
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
- // Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
- PetscIS subpointIS = NULL;
- PetscInt cStart, cEnd;
- PetscErrorCode err;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
+ // Get 'surface' cells (1 dimension lower than top-level cells)
+ const PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
// Get sections
topology::Field<topology::SubMesh>& dampingConsts = _parameters->get("damping constants");
+ topology::VecVisitorMesh dampingConstsVisitor(dampingConsts);
+ PetscScalar* dampingConstsArray = dampingConstsVisitor.localArray();
- const topology::Field<topology::Mesh>& solution = fields->solution();
- PetscSection solutionSection = solution.petscSection();assert(solutionSection);
- PetscVec solutionVec = solution.localVector();assert(solutionVec);
- PetscSection solutionGlobalSection, solutionSubsection, solutionGlobalSubsection;
- PetscSF sf = NULL;
-
- err = PetscSectionCreateSubmeshSection(solutionSection, subpointIS, &solutionSubsection);CHECK_PETSC_ERROR(err);
- err = DMGetPointSF(solution.dmMesh(), &sf);CHECK_PETSC_ERROR(err);
- err = PetscSectionCreateGlobalSection(solutionSection, sf, PETSC_FALSE, &solutionGlobalSection);CHECK_PETSC_ERROR(err);
- err = PetscSectionCreateSubmeshSection(solutionGlobalSection, subpointIS, &solutionGlobalSubsection);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
-
// Get sparse matrix
+ const topology::Field<topology::Mesh>& solution = fields->solution();
+ topology::SubMeshIS submeshIS(*_boundaryMesh);
const PetscMat jacobianMat = jacobian->matrix();assert(jacobianMat);
+ topology::MatVisitorSubMesh jacobianVisitor(jacobianMat, solution, submeshIS);
+ submeshIS.deallocate();
// Get parameters used in integration.
const PylithScalar dt = _dt;
@@ -590,12 +560,10 @@
_initCellMatrix();
#if !defined(PRECOMPUTE_GEOMETRY)
+ PetscScalar *coordsCell = NULL;
+ PetscInt coordsSize = 0;
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ topology::CoordsVisitor coordsVisitor(dmSubMesh);
#endif
_logger->eventEnd(setupEvent);
@@ -603,8 +571,6 @@
_logger->eventBegin(computeEvent);
#endif
- PetscScalar* dampingConstsArray = dampingConsts.getLocalArray();
-
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
@@ -613,14 +579,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
#error("Code for PRECOMPUTE_GEOMETRY not implemented")
#else
- PetscScalar *coordsArray;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
- for (PetscInt i = 0; i < coordsSize; ++i) {
- coordinatesCell[i] = coordsArray[i];
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+ for (PetscInt i = 0; i < coordsSize; ++i) { // :TODO: Remove copy.
+ coordinatesCell[i] = coordsCell[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -628,8 +592,8 @@
#endif
// Get damping constants
- const PetscInt doff = dampingConsts.sectionOffset(c);
- assert(numQuadPts*spaceDim == dampingConsts.sectionDof(c));
+ const PetscInt doff = dampingConstsVisitor.sectionOffset(c);
+ assert(numQuadPts*spaceDim == dampingConstsVisitor.sectionDof(c));
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@ -666,17 +630,12 @@
#endif
// Assemble cell contribution into PETSc Matrix
- err = DMPlexMatSetClosure(subMesh, solutionSubsection, solutionGlobalSubsection, jacobianMat, c, &_cellMatrix[0], ADD_VALUES);
- CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
+ jacobianVisitor.setClosure(&_cellMatrix[0], _cellMatrix.size(), c, ADD_VALUES);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
- dampingConsts.restoreLocalArray(&dampingConstsArray);
- err = PetscSectionDestroy(&solutionSubsection);CHECK_PETSC_ERROR(err);
- err = PetscSectionDestroy(&solutionGlobalSection);CHECK_PETSC_ERROR(err);
- err = PetscSectionDestroy(&solutionGlobalSubsection);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops((cEnd-cStart)*numQuadPts*(3+numBasis*(1+numBasis*(1+2*spaceDim))));
@@ -714,13 +673,11 @@
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
- // Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
- PetscIS subpointIS = NULL;
- PetscInt cStart, cEnd;
- PetscErrorCode err;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
+ // Get 'surface' cells (1 dimension lower than top-level cells)
+ const PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
// Get parameters used in integration.
const PylithScalar dt = _dt;
@@ -732,20 +689,19 @@
// Get sections
topology::Field<topology::SubMesh>& dampingConsts = _parameters->get("damping constants");
-
- PetscSection jacobianSection = jacobian->petscSection();assert(jacobianSection);
- PetscVec jacobianVec = jacobian->localVector();assert(jacobianVec);
- PetscSection jacobianSubsection;
- err = PetscSectionCreateSubmeshSection(jacobianSection, subpointIS, &jacobianSubsection);
- err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
+ topology::VecVisitorMesh dampingConstsVisitor(dampingConsts);
+ PetscScalar* dampingConstsArray = dampingConstsVisitor.localArray();
+ topology::SubMeshIS submeshIS(*_boundaryMesh);
+ topology::VecVisitorSubMesh jacobianVisitor(*jacobian, submeshIS);
+
+ submeshIS.deallocate();
+
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
+ topology::CoordsVisitor coordsVisitor(dmSubMesh);
#endif
_logger->eventEnd(setupEvent);
@@ -753,8 +709,6 @@
_logger->eventBegin(computeEvent);
#endif
- PetscScalar* dampingConstsArray = dampingConsts.getLocalArray();
-
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
@@ -763,14 +717,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
#error("Code for PRECOMPUTE_GEOMETRY not implemented");
#else
- PetscScalar *coordsArray;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
- for (PetscInt i = 0; i < coordsSize; ++i) {
- coordinatesCell[i] = coordsArray[i];
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+ for (PetscInt i = 0; i < coordsSize; ++i) { // :TODO: Remove copy.
+ coordinatesCell[i] = coordsCell[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -778,8 +730,8 @@
#endif
// Get damping constants
- const PetscInt doff = dampingConsts.sectionOffset(c);
- assert(numQuadPts*spaceDim == dampingConsts.sectionDof(c));
+ const PetscInt doff = dampingConstsVisitor.sectionOffset(c);
+ assert(numQuadPts*spaceDim == dampingConstsVisitor.sectionDof(c));
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
@@ -808,7 +760,9 @@
* dampingConstsArray[doff+iQuad * spaceDim + iDim];
} // for
} // for
- err = DMPlexVecSetClosure(subMesh, jacobianSubsection, jacobianVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+
+ jacobianVisitor.setClosure(&_cellVector[0], _cellVector.size(), c, ADD_VALUES);
+
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(4+numBasis+numBasis*(1+spaceDim*2)));
_logger->eventEnd(computeEvent);
@@ -818,8 +772,6 @@
_logger->eventEnd(updateEvent);
#endif
} // for
- dampingConsts.restoreLocalArray(&dampingConstsArray);
- err = PetscSectionDestroy(&jacobianSubsection);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
PetscLogFlops((cEnd-cStart)*numQuadPts*(4+numBasis+numBasis*(1+spaceDim*2)));
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.hh 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.hh 2013-03-13 20:29:41 UTC (rev 21527)
@@ -109,8 +109,8 @@
* @param fields Solution fields
*/
void integrateResidualLumped(const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields);
+ const PylithScalar t,
+ topology::SolutionFields* const fields);
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -23,6 +23,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -184,28 +185,28 @@
// Calculate spatial and temporal variation of value for BC.
_calculateValue(t);
+ // Get sections
assert(_parameters);
topology::Field<topology::Mesh>& valueField = _parameters->get("value");
- PetscSection fieldSection = field.petscSection();assert(fieldSection);
+ topology::VecVisitorMesh valueVisitor(valueField);
+ PetscScalar* valueArray = valueVisitor.localArray();
- PetscScalar* valueArray = valueField.getLocalArray();
- PetscScalar* fieldArray = field.getLocalArray();
+ topology::VecVisitorMesh fieldVisitor(field);
+ PetscScalar* fieldArray = fieldVisitor.localArray();
const int numPoints = _points.size();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
PetscInt p_bc = _points[iPoint];
- const PetscInt off = field.sectionOffset(p_bc);
- const PetscInt voff = valueField.sectionOffset(p_bc);
- assert(numFixedDOF == valueField.sectionDof(p_bc));
+ const PetscInt off = fieldVisitor.sectionOffset(p_bc);
+ const PetscInt voff = valueVisitor.sectionOffset(p_bc);
+ assert(numFixedDOF == valueVisitor.sectionDof(p_bc));
for(PetscInt iDOF = 0; iDOF < numFixedDOF; ++iDOF) {
- assert(_bcDOF[iDOF] < field.sectionDof(p_bc));
+ assert(_bcDOF[iDOF] < fieldVisitor.sectionDof(p_bc));
fieldArray[_bcDOF[iDOF]+off] = valueArray[voff+iDOF];
} // for
} // for
- valueField.restoreLocalArray(&valueArray);
- field.restoreLocalArray(&fieldArray);
} // setField
// ----------------------------------------------------------------------
@@ -222,27 +223,27 @@
// Calculate spatial and temporal variation of value for BC.
_calculateValueIncr(t0, t1);
+ // Get sections
assert(_parameters);
topology::Field<topology::Mesh>& valueField = _parameters->get("value");
- PetscSection fieldSection = field.petscSection();assert(fieldSection);
+ topology::VecVisitorMesh valueVisitor(valueField);
+ PetscScalar* valueArray = valueVisitor.localArray();
- PetscScalar* valueArray = valueField.getLocalArray();
- PetscScalar* fieldArray = field.getLocalArray();
+ topology::VecVisitorMesh fieldVisitor(field);
+ PetscScalar* fieldArray = fieldVisitor.localArray();
const int numPoints = _points.size();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
PetscInt p_bc = _points[iPoint];
- const PetscInt off = field.sectionOffset(p_bc);
- const PetscInt voff = valueField.sectionOffset(p_bc);
- assert(numFixedDOF == valueField.sectionDof(p_bc));
+ const PetscInt off = fieldVisitor.sectionOffset(p_bc);
+ const PetscInt voff = valueVisitor.sectionOffset(p_bc);
+ assert(numFixedDOF == valueVisitor.sectionDof(p_bc));
for(PetscInt iDOF=0; iDOF < numFixedDOF; ++iDOF) {
- assert(_bcDOF[iDOF] < field.sectionDof(p_bc));
+ assert(_bcDOF[iDOF] < fieldVisitor.sectionDof(p_bc));
fieldArray[_bcDOF[iDOF]+off] = valueArray[voff+iDOF];
} // for
} // for
- valueField.restoreLocalArray(&valueArray);
- field.restoreLocalArray(&fieldArray);
} // setFieldIncr
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -24,6 +24,7 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/Fields.hh" // USES Fields
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -89,17 +90,21 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("BoundaryConditions");
- if (0 == _outputFields)
+ if (!_outputFields) {
_outputFields = new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
+ } // if
assert(_outputFields);
_outputFields->add("buffer (vector)", "buffer_vector", topology::FieldBase::CELLS_FIELD, spaceDim);
- _outputFields->get("buffer (vector)").vectorFieldType(topology::FieldBase::VECTOR);
- _outputFields->get("buffer (vector)").scale(lengthScale);
- _outputFields->get("buffer (vector)").allocate();
+ topology::Field<topology::SubMesh>& bufferVector = _outputFields->get("buffer (vector)");
+ bufferVector.vectorFieldType(topology::FieldBase::VECTOR);
+ bufferVector.scale(lengthScale);
+ bufferVector.allocate();
+
_outputFields->add("buffer (scalar)", "buffer_scalar", topology::FieldBase::CELLS_FIELD, 1);
- _outputFields->get("buffer (scalar)").vectorFieldType(topology::FieldBase::SCALAR);
- _outputFields->get("buffer (scalar)").scale(timeScale);
- _outputFields->get("buffer (scalar)").allocate();
+ topology::Field<topology::SubMesh>& bufferScalar = _outputFields->get("buffer (scalar)");
+ bufferScalar.vectorFieldType(topology::FieldBase::SCALAR);
+ bufferScalar.scale(timeScale);
+ bufferScalar.allocate();
logger.stagePop();
@@ -132,11 +137,6 @@
const char* label,
const PylithScalar scale)
{ // _bufferVector
- typedef topology::SubMesh::SieveMesh SieveMesh;
- typedef topology::Mesh::RealUniformSection RealUniformSection;
- typedef topology::SubMesh::RealSection SubRealSection;
- typedef topology::SubMesh::RealUniformSection SubRealUniformSection;
-
assert(_boundaryMesh);
assert(_parameters);
assert(_outputFields);
@@ -148,43 +148,33 @@
throw std::runtime_error(msg.str());
} // if
+ assert(_outputFields->hasField("buffer (vector)"));
+ topology::Field<topology::SubMesh>& buffer = _outputFields->get("buffer (vector)");
+ topology::VecVisitorMesh bufferVisitor(buffer);
+ PetscScalar* bufferArray = bufferVisitor.localArray();
+
+ topology::Field<topology::Mesh>& field = _parameters->get(name);
+ topology::VecVisitorMesh fieldVisitor(field);
+ PetscScalar* fieldArray = fieldVisitor.localArray();
+
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
assert(cs);
const int spaceDim = cs->spaceDim();
const int numPoints = _points.size();
const int numFixedDOF = _bcDOF.size();
- PetscErrorCode err = 0;
-
- assert(_outputFields->hasField("buffer (vector)"));
- PetscSection outputSection = _outputFields->get("buffer (vector)").petscSection();assert(outputSection);
- PetscVec outputVec = _outputFields->get("buffer (vector)").localVector();assert(outputVec);
- PetscScalar *outputArray = NULL;
- err = VecGetArray(outputVec, &outputArray);CHECK_PETSC_ERROR(err);
-
- PetscSection fieldSection = _parameters->get(name).petscSection();assert(fieldSection);
- PetscVec fieldVec = _parameters->get(name).localVector();assert(fieldVec);
- PetscScalar *fieldArray = NULL;
- err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
-
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const SieveMesh::point_type point = _points[iPoint];
+ const PetscInt point = _points[iPoint];
- PetscInt odof, ooff, fdof, foff;
- err = PetscSectionGetDof(outputSection, point, &odof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(outputSection, point, &ooff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(fieldSection, point, &fdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(fieldSection, point, &foff);CHECK_PETSC_ERROR(err);
- assert(spaceDim == odof);
- assert(fdof == numFixedDOF);
+ const PetscInt boff = bufferVisitor.sectionOffset(point);
+ const PetscInt foff = fieldVisitor.sectionOffset(point);
+ assert(spaceDim == bufferVisitor.sectionDof(point));
+ assert(numFixedDOF == fieldVisitor.sectionDof(point));
for(PetscInt iDOF=0; iDOF < numFixedDOF; ++iDOF)
- outputArray[ooff+_bcDOF[iDOF]] = fieldArray[foff+iDOF];
+ bufferArray[boff+_bcDOF[iDOF]] = fieldArray[foff+iDOF];
} // for
- err = VecRestoreArray(outputVec, &outputArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
- topology::Field<topology::SubMesh>& buffer = _outputFields->get("buffer (vector)");
buffer.label(label);
buffer.scale(scale);
@@ -198,11 +188,6 @@
const char* label,
const PylithScalar scale)
{ // _bufferScalar
- typedef topology::SubMesh::SieveMesh SieveMesh;
- typedef topology::Mesh::RealUniformSection RealUniformSection;
- typedef topology::SubMesh::RealSection SubRealSection;
- typedef topology::SubMesh::RealUniformSection SubRealUniformSection;
-
assert(_boundaryMesh);
assert(_parameters);
assert(_outputFields);
@@ -214,35 +199,27 @@
throw std::runtime_error(msg.str());
} // if
- const int numPoints = _points.size();
- PetscErrorCode err = 0;
-
assert(_outputFields->hasField("buffer (scalar)"));
- PetscSection outputSection = _outputFields->get("buffer (scalar)").petscSection();assert(outputSection);
- PetscVec outputVec = _outputFields->get("buffer (scalar)").localVector();assert(outputVec);
- PetscScalar *outputArray = NULL;
- err = VecGetArray(outputVec, &outputArray);CHECK_PETSC_ERROR(err);
-
- PetscSection fieldSection = _parameters->get(name).petscSection();assert(fieldSection);
- PetscVec fieldVec = _parameters->get(name).localVector();assert(fieldVec);
- PetscScalar *fieldArray = NULL;
- err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::SubMesh>& buffer = _outputFields->get("buffer (scalar)");
+ topology::VecVisitorMesh bufferVisitor(buffer);
+ PetscScalar* bufferArray = bufferVisitor.localArray();
+ topology::Field<topology::Mesh>& field = _parameters->get(name);
+ topology::VecVisitorMesh fieldVisitor(field);
+ PetscScalar* fieldArray = fieldVisitor.localArray();
+
+ const int numPoints = _points.size();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const SieveMesh::point_type point = _points[iPoint];
+ const PetscInt point = _points[iPoint];
- PetscInt odof, ooff, fdof, foff;
- err = PetscSectionGetDof(outputSection, point, &odof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(outputSection, point, &ooff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(fieldSection, point, &fdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(fieldSection, point, &foff);CHECK_PETSC_ERROR(err);
- assert(1 == odof);
- assert(fdof == 1);
+ const PetscInt boff = bufferVisitor.sectionOffset(point);
+ const PetscInt foff = fieldVisitor.sectionOffset(point);
+ assert(1 == bufferVisitor.sectionDof(point));
+ assert(1 == fieldVisitor.sectionDof(point));
- outputArray[ooff] = fieldArray[foff];
+ bufferArray[boff] = fieldArray[foff];
} // for
- topology::Field<topology::SubMesh>& buffer = _outputFields->get("buffer (scalar)");
buffer.label(label);
buffer.scale(scale);
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -23,6 +23,11 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Fields.hh" // HOLDSA Fields
#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/CoordsVisitor.hh" // USES CoordsVisitor
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
+#include "pylith/topology/VisitorSubMesh.hh" // USES VecVisitorSubMesh
+#include "pylith/topology/Stratum.hh" // USES Stratum
+
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -97,55 +102,48 @@
scalar_array tractionsCell(numQuadPts*spaceDim);
// Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
- PetscIS subpointIS;
- PetscInt cStart, cEnd;
- PetscErrorCode err = 0;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
+ PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
// Get sections
_calculateValue(t);
topology::Field<topology::SubMesh>& valueField = _parameters->get("value");
-
- PetscSection residualSection = residual.petscSection();assert(residualSection);
- PetscVec residualVec = residual.localVector();assert(residualVec);
- PetscSection residualSubsection = NULL;
- err = PetscSectionCreateSubmeshSection(residualSection, subpointIS, &residualSubsection);CHECK_PETSC_ERROR(err);
- err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
+ topology::VecVisitorMesh valueVisitor(valueField);
+ PetscScalar* valueArray = valueVisitor.localArray();
+ // Get subsections
+ topology::SubMeshIS submeshIS(*_boundaryMesh);
+ topology::VecVisitorSubMesh residualVisitor(residual, submeshIS);
+ submeshIS.deallocate();
+
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
+ topology::CoordsVisitor coordsVisitor(dmSubMesh);
#endif
- PetscScalar* valueArray = valueField.getLocalArray();
-
// Loop over faces and integrate contribution from each face
for(PetscInt c = cStart; c < cEnd; ++c) {
#if defined(PRECOMPUTE_GEOMETRY)
#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for (PetscInt i=0; i < coordsSize; ++i) {
- coordinatesCell[i] = coords[i];
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+ for (PetscInt i=0; i < coordsSize; ++i) { // :TODO: Remove copy.
+ coordinatesCell[i] = coordsCell[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
#endif
// Reset element vector to zero
_resetCellVector();
// Restrict tractions to cell
- const PetscInt voff = valueField.sectionOffset(c);
- assert(numQuadPts*spaceDim == valueField.sectionDof(c));
+ const PetscInt voff = valueVisitor.sectionOffset(c);
+ assert(numQuadPts*spaceDim == valueVisitor.sectionDof(c));
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -163,10 +161,11 @@
} // for
} // for
} // for
- err = DMPlexVecSetClosure(subMesh, residualSubsection, residualVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+
+ residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), c, ADD_VALUES);
+
PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*(1+2*spaceDim))));
} // for
- valueField.restoreLocalArray(&valueArray);
} // integrateResidual
// ----------------------------------------------------------------------
@@ -390,10 +389,10 @@
assert(_parameters);
// Get 'surface' cells (1 dimension lower than top-level cells)
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
- PetscInt cStart, cEnd;
- PetscErrorCode err = 0;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
const int numBasis = _quadrature->numBasis();
@@ -406,15 +405,16 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get sections.
- scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
-
topology::Field<topology::SubMesh>& valueField = _parameters->get(name);
+ topology::VecVisitorMesh valueVisitor(valueField);
+ PetscScalar* valueArray = valueVisitor.localArray();
+ // Get coordinates
+ scalar_array coordinatesCell(numBasis*spaceDim);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
+ topology::CoordsVisitor coordsVisitor(dmSubMesh);
+
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
assert(cs);
@@ -426,8 +426,6 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->computeGeometry(*_boundaryMesh, cells);
#endif
-
- PetscScalar* valueArray = valueField.getLocalArray();
// Loop over cells in boundary mesh and perform queries.
for(PetscInt c = cStart; c < cEnd; ++c) {
@@ -435,14 +433,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for (PetscInt i = 0; i < coordsSize; ++i) {
- coordinatesCell[i] = coords[i];
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+ for (PetscInt i = 0; i < coordsSize; ++i) { // :TODO: Remove copy.
+ coordinatesCell[i] = coordsCell[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
#endif
const scalar_array& quadPtsNondim = _quadrature->quadPts();
quadPtsGlobal = quadPtsNondim;
@@ -472,7 +468,6 @@
for(PetscInt d = 0; d < vdof; ++d)
valueArray[voff+d] = valuesCell[d];
} // for
- valueField.restoreLocalArray(&valueArray);
} // _queryDB
// ----------------------------------------------------------------------
@@ -490,10 +485,10 @@
up[i] = upDir[i];
// Get 'surface' cells (1 dimension lower than top-level cells)
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
- PetscInt cStart, cEnd;
- PetscErrorCode err = 0;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
// Quadrature related values.
const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
@@ -511,30 +506,27 @@
PylithScalar jacobianDet = 0;
scalar_array orientation(orientationSize);
- // Get sections.
+ // Get coordinates.
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection = NULL;
- PetscVec coordVec = NULL;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);assert(coordSection);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);assert(coordVec);
+ PetscScalar* coordsCell = NULL;
+ PetscInt coordsSize = 0;
+ topology::CoordsVisitor coordsVisitor(dmSubMesh);
- scalar_array parametersCellLocal(spaceDim);
-
- topology::Field<topology::SubMesh>& valueField = _parameters->get("value");
- PetscScalar* valueArray = valueField.getLocalArray();
-
+ // Get sections
+ scalar_array tmpLocal(spaceDim);
+
topology::Field<topology::SubMesh>* initialField = (_dbInitial) ? &_parameters->get("initial") : 0;
+ topology::VecVisitorMesh* initialVisitor = (initialField) ? new topology::VecVisitorMesh(*initialField) : 0;
+ PetscScalar* initialArray = (initialVisitor) ? initialVisitor->localArray() : NULL;
+
topology::Field<topology::SubMesh>* rateField = (_dbRate) ? &_parameters->get("rate") : 0;
- topology::Field<topology::SubMesh>* rateTimeField = (_dbRate) ? &_parameters->get("rate time") : 0;
+ topology::VecVisitorMesh* rateVisitor = (rateField) ? new topology::VecVisitorMesh(*rateField) : 0;
+ PetscScalar* rateArray = (rateVisitor) ? rateVisitor->localArray() : NULL;
+
topology::Field<topology::SubMesh>* changeField = (_dbChange) ? &_parameters->get("change") : 0;
- topology::Field<topology::SubMesh>* changeTimeField = (_dbChange) ? &_parameters->get("change time") : 0;
+ topology::VecVisitorMesh* changeVisitor = (changeField) ? new topology::VecVisitorMesh(*changeField) : 0;
+ PetscScalar* changeArray = (changeVisitor) ? changeVisitor->localArray() : NULL;
- PetscScalar* initialArray = (_dbInitial) ? initialField->getLocalArray() : NULL;
- PetscScalar* rateArray = (_dbRate) ? rateField->getLocalArray() : NULL;
- PetscScalar* rateTimeArray = (_dbRate) ? rateTimeField->getLocalArray() : NULL;
- PetscScalar* changeArray = (_dbChange) ? changeField->getLocalArray() : NULL;
- PetscScalar* changeTimeArray = (_dbChange) ? changeTimeField->getLocalArray() : NULL;
-
// Loop over cells in boundary mesh, compute orientations, and then
// rotate corresponding traction vector from local to global coordinates.
for(PetscInt c = cStart; c < cEnd; ++c) {
@@ -542,14 +534,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
- PetscScalar *coords;
- PetscInt coordsSize;
- err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for (PetscInt i=0; i < coordsSize; ++i) {
- coordinatesCell[i] = coords[i];
+ coordsVisitor.getClosure(&coordsCell, &coordsSize, c);
+ for (PetscInt i=0; i < coordsSize; ++i) { // :TODO: Remove copy.
+ coordinatesCell[i] = coordsCell[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ coordsVisitor.restoreClosure(&coordsCell, &coordsSize, c);
#endif
for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
// Compute Jacobian and determinant at quadrature point, then get orientation.
@@ -564,75 +554,57 @@
if (_dbInitial) {
// Rotate traction vector from local coordinate system to global coordinate system
- PylithScalar *initialLocal = ¶metersCellLocal[0];
- assert(initialField);
- const PetscInt ioff = initialField->sectionOffset(c);
- assert(numQuadPts*spaceDim == initialField->sectionDof(c));
+ assert(initialVisitor);
+ const PetscInt ioff = initialVisitor->sectionOffset(c);
+ assert(numQuadPts*spaceDim == initialVisitor->sectionDof(c));
for(int iDim = 0; iDim < spaceDim; ++iDim) {
- initialLocal[iDim] = initialArray[ioff+iSpace+iDim];
+ tmpLocal[iDim] = initialArray[ioff+iSpace+iDim];
} // for
for(int iDim = 0; iDim < spaceDim; ++iDim) {
initialArray[ioff+iSpace+iDim] = 0.0;
for(int jDim = 0; jDim < spaceDim; ++jDim) {
- initialArray[ioff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * initialLocal[jDim];
+ initialArray[ioff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * tmpLocal[jDim];
} // for
} // for
} // if
if (_dbRate) {
// Rotate traction vector from local coordinate system to global coordinate system
- PylithScalar *rateLocal = ¶metersCellLocal[0];
- assert(rateField);
- const PetscInt roff = rateField->sectionOffset(c);
- assert(numQuadPts*spaceDim == rateField->sectionDof(c));
+ assert(rateVisitor);
+ const PetscInt roff = rateVisitor->sectionOffset(c);
+ assert(numQuadPts*spaceDim == rateVisitor->sectionDof(c));
for(int iDim = 0; iDim < spaceDim; ++iDim) {
- rateLocal[iDim] = rateArray[roff+iSpace+iDim];
+ tmpLocal[iDim] = rateArray[roff+iSpace+iDim];
} // for
for(int iDim = 0; iDim < spaceDim; ++iDim) {
rateArray[roff+iSpace+iDim] = 0.0;
for(int jDim = 0; jDim < spaceDim; ++jDim) {
- rateArray[roff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * rateLocal[jDim];
+ rateArray[roff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * tmpLocal[jDim];
} // for
} // for
} // if
if (_dbChange) {
// Rotate traction vector from local coordinate system to global coordinate system
- PylithScalar *changeLocal = ¶metersCellLocal[0];
- assert(changeField);
- const PetscInt coff = changeField->sectionOffset(c);
- assert(numQuadPts*spaceDim == changeField->sectionDof(c));
+ assert(changeVisitor);
+ const PetscInt coff = changeVisitor->sectionOffset(c);
+ assert(numQuadPts*spaceDim == changeVisitor->sectionDof(c));
for (int iDim = 0; iDim < spaceDim; ++iDim) {
- changeLocal[iDim] = changeArray[coff+iSpace+iDim];
+ tmpLocal[iDim] = changeArray[coff+iSpace+iDim];
} // for
for(int iDim = 0; iDim < spaceDim; ++iDim) {
changeArray[coff+iSpace+iDim] = 0.0;
for(int jDim = 0; jDim < spaceDim; ++jDim) {
- changeArray[coff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * changeLocal[jDim];
+ changeArray[coff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * tmpLocal[jDim];
} // for
} // for
} // if
} // for
} // for
- valueField.restoreLocalArray(&valueArray);
- if (_dbInitial) {
- assert(initialField);
- initialField->restoreLocalArray(&initialArray);
- } // if
- if (_dbRate) {
- assert(rateField);
- rateField->restoreLocalArray(&rateArray);
- assert(rateTimeField);
- rateTimeField->restoreLocalArray(&rateTimeArray);
- } // if
- if (_dbChange) {
- assert(changeField);
- changeField->restoreLocalArray(&changeArray);
- assert(changeTimeField);
- changeTimeField->restoreLocalArray(&changeTimeArray);
- } // if
-
+ delete initialVisitor; initialVisitor = 0;
+ delete rateVisitor; rateVisitor = 0;
+ delete changeVisitor; changeVisitor = 0;
} // paramsLocalToGlobal
// ----------------------------------------------------------------------
@@ -647,33 +619,42 @@
const PylithScalar timeScale = _getNormalizer().timeScale();
// Get 'surface' cells (1 dimension lower than top-level cells)
- PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
- PetscInt cStart, cEnd;
- PetscErrorCode err = 0;
- err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ PetscDM dmSubMesh = _boundaryMesh->dmMesh();assert(dmSubMesh);
+ topology::Stratum heightStratum(dmSubMesh, topology::Stratum::HEIGHT, 1);
+ const PetscInt cStart = heightStratum.begin();
+ const PetscInt cEnd = heightStratum.end();
const int spaceDim = _quadrature->spaceDim();
const int numQuadPts = _quadrature->numQuadPts();
-
+ // Get sections
topology::Field<topology::SubMesh>& valueField = _parameters->get("value");
- PetscScalar* valueArray = valueField.getLocalArray();
-
+ topology::VecVisitorMesh valueVisitor(valueField);
+ PetscScalar* valueArray = valueVisitor.localArray();
+
topology::Field<topology::SubMesh>* initialField = (_dbInitial) ? &_parameters->get("initial") : 0;
+ topology::VecVisitorMesh* initialVisitor = (initialField) ? new topology::VecVisitorMesh(*initialField) : 0;
+ PetscScalar* initialArray = (initialVisitor) ? initialVisitor->localArray() : NULL;
+
topology::Field<topology::SubMesh>* rateField = (_dbRate) ? &_parameters->get("rate") : 0;
+ topology::VecVisitorMesh* rateVisitor = (rateField) ? new topology::VecVisitorMesh(*rateField) : 0;
+ PetscScalar* rateArray = (rateVisitor) ? rateVisitor->localArray() : NULL;
+
topology::Field<topology::SubMesh>* rateTimeField = (_dbRate) ? &_parameters->get("rate time") : 0;
+ topology::VecVisitorMesh* rateTimeVisitor = (rateTimeField) ? new topology::VecVisitorMesh(*rateTimeField) : 0;
+ PetscScalar* rateTimeArray = (rateTimeVisitor) ? rateTimeVisitor->localArray() : NULL;
+
topology::Field<topology::SubMesh>* changeField = (_dbChange) ? &_parameters->get("change") : 0;
+ topology::VecVisitorMesh* changeVisitor = (changeField) ? new topology::VecVisitorMesh(*changeField) : 0;
+ PetscScalar* changeArray = (changeVisitor) ? changeVisitor->localArray() : NULL;
+
topology::Field<topology::SubMesh>* changeTimeField = (_dbChange) ? &_parameters->get("change time") : 0;
+ topology::VecVisitorMesh* changeTimeVisitor = (changeTimeField) ? new topology::VecVisitorMesh(*changeTimeField) : 0;
+ PetscScalar* changeTimeArray = (changeTimeVisitor) ? changeTimeVisitor->localArray() : NULL;
- PetscScalar* initialArray = (_dbInitial) ? initialField->getLocalArray() : NULL;
- PetscScalar* rateArray = (_dbRate) ? rateField->getLocalArray() : NULL;
- PetscScalar* rateTimeArray = (_dbRate) ? rateTimeField->getLocalArray() : NULL;
- PetscScalar* changeArray = (_dbChange) ? changeField->getLocalArray() : NULL;
- PetscScalar* changeTimeArray = (_dbChange) ? changeTimeField->getLocalArray() : NULL;
-
for(PetscInt c = cStart; c < cEnd; ++c) {
- const PetscInt voff = valueField.sectionOffset(c);
- const PetscInt vdof = valueField.sectionDof(c);
+ const PetscInt voff = valueVisitor.sectionOffset(c);
+ const PetscInt vdof = valueVisitor.sectionDof(c);
assert(numQuadPts*spaceDim == vdof);
for (PetscInt d = 0; d < vdof; ++d) {
valueArray[voff+d] = 0.0;
@@ -681,9 +662,9 @@
// Contribution from initial value
if (_dbInitial) {
- assert(initialField);
- const PetscInt ioff = initialField->sectionOffset(c);
- const PetscInt idof = initialField->sectionDof(c);
+ assert(initialVisitor);
+ const PetscInt ioff = initialVisitor->sectionOffset(c);
+ const PetscInt idof = initialVisitor->sectionDof(c);
assert(numQuadPts*spaceDim == idof);
for (PetscInt d = 0; d < idof; ++d) {
valueArray[voff+d] += initialArray[ioff+d];
@@ -692,13 +673,13 @@
// Contribution from rate of change of value
if (_dbRate) {
- assert(rateField);
- const PetscInt roff = rateField->sectionOffset(c);
- const PetscInt rdof = rateField->sectionDof(c);
+ assert(rateVisitor);
+ const PetscInt roff = rateVisitor->sectionOffset(c);
+ const PetscInt rdof = rateVisitor->sectionDof(c);
assert(numQuadPts*spaceDim == rdof);
- assert(rateTimeField);
- const PetscInt rtoff = rateTimeField->sectionOffset(c);
- assert(numQuadPts == rateTimeField->sectionDof(c));
+ assert(rateTimeVisitor);
+ const PetscInt rtoff = rateTimeVisitor->sectionOffset(c);
+ assert(numQuadPts == rateTimeVisitor->sectionDof(c));
for(int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
const PylithScalar tRel = t - rateTimeArray[rtoff+iQuad];
@@ -711,13 +692,12 @@
// Contribution from change of value
if (_dbChange) {
- assert(changeField);
- const PetscInt coff = changeField->sectionOffset(c);
- const PetscInt cdof = changeField->sectionDof(c);
- assert(numQuadPts*spaceDim == cdof);
+ assert(changeVisitor);
+ const PetscInt coff = changeVisitor->sectionOffset(c);
+ const PetscInt cdof = changeVisitor->sectionDof(c);
assert(changeTimeField);
- const PetscInt ctoff = changeTimeField->sectionOffset(c);
- assert(numQuadPts == changeTimeField->sectionDof(c));
+ const PetscInt ctoff = changeTimeVisitor->sectionOffset(c);
+ assert(numQuadPts == changeTimeVisitor->sectionDof(c));
for(int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
const PylithScalar tRel = t - changeTimeArray[ctoff+iQuad];
@@ -736,29 +716,18 @@
} // if
} // if
for (int iDim = 0; iDim < spaceDim; ++iDim) {
- valueArray[voff+iQuad*spaceDim+iDim] += changeArray[coff+iQuad*spaceDim+iDim] * scale;
+ valueArray[voff+iQuad*spaceDim+iDim] += changeArray[coff+iQuad*spaceDim+iDim]*scale;
} // for
} // if
} // for
} // if
} // for
- valueField.restoreLocalArray(&valueArray);
- if (_dbInitial) {
- assert(initialField);
- initialField->restoreLocalArray(&initialArray);
- } // if
- if (_dbRate) {
- assert(rateField);
- rateField->restoreLocalArray(&rateArray);
- assert(rateTimeField);
- rateTimeField->restoreLocalArray(&rateTimeArray);
- } // if
- if (_dbChange) {
- assert(changeField);
- changeField->restoreLocalArray(&changeArray);
- assert(changeTimeField);
- changeTimeField->restoreLocalArray(&changeTimeArray);
- } // if
+
+ delete initialVisitor; initialVisitor = 0;
+ delete rateVisitor; rateVisitor = 0;
+ delete rateTimeVisitor; rateTimeVisitor = 0;
+ delete changeVisitor; changeVisitor = 0;
+ delete changeTimeVisitor; changeTimeVisitor = 0;
} // _calculateValue
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -23,6 +23,7 @@
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/VisitorMesh.hh" // USES VecVisitorMesh
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -85,52 +86,43 @@
// Calculate spatial and temporal variation of value for BC.
_calculateValue(t);
- const int numPoints = _points.size();
- const int numBCDOF = _bcDOF.size();
-
const topology::Mesh& mesh = residual.mesh();
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
assert(cs);
const int spaceDim = cs->spaceDim();
- PetscSection residualSection = residual.petscSection();assert(residualSection);
- PetscVec residualVec = residual.localVector();assert(residualVec);
- PetscScalar *residualArray;
-
+ topology::VecVisitorMesh residualVisitor(residual);
+ PetscScalar* residualArray = residualVisitor.localArray();
+
+ topology::Field<topology::Mesh>& valueField = _parameters->get("value");
+ topology::VecVisitorMesh valueVisitor(valueField);
+ PetscScalar* valueArray = valueVisitor.localArray();
+
// Get global order
PetscDM dmMesh = residual.dmMesh();assert(dmMesh);
PetscSection globalSection = NULL;
PetscErrorCode err;
err = DMGetDefaultGlobalSection(dmMesh, &globalSection);CHECK_PETSC_ERROR(err);
- PetscSection valueSection = _parameters->get("value").petscSection();assert(valueSection);
- PetscVec valueVec = _parameters->get("value").localVector();assert(valueVec);
- PetscScalar* valueArray = NULL;
- err = VecGetArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
-
+ const int numPoints = _points.size();
+ const int numBCDOF = _bcDOF.size();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const PetscInt p_bc = _points[iPoint]; // Get point label.
- PetscInt goff;
// Contribute to residual only if point is local.
+ PetscInt goff = -1;
err = PetscSectionGetOffset(globalSection, p_bc, &goff);CHECK_PETSC_ERROR(err);
if (goff < 0) continue;
- PetscInt vdof, voff, rdof, roff;
+ const PetscInt roff = residualVisitor.sectionOffset(p_bc);
+ const PetscInt voff = valueVisitor.sectionOffset(p_bc);
+ assert(numBCDOF == valueVisitor.sectionDof(p_bc));
+ assert(spaceDim == residualVisitor.sectionDof(p_bc));
- err = PetscSectionGetDof(valueSection, p_bc, &vdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(valueSection, p_bc, &voff);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetDof(residualSection, p_bc, &rdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(residualSection, p_bc, &roff);CHECK_PETSC_ERROR(err);
- assert(vdof == numBCDOF);
- assert(rdof == spaceDim);
-
- for (int iDOF=0; iDOF < numBCDOF; ++iDOF)
+ for (int iDOF=0; iDOF < numBCDOF; ++iDOF) {
residualArray[roff+_bcDOF[iDOF]] += valueArray[voff+iDOF];
+ } // for
} // for
- err = VecRestoreArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
} // integrateResidualAssembled
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -219,9 +219,9 @@
// Query database for values.
void
pylith::bc::TimeDependentPoints::_queryDB(const char* name,
- spatialdata::spatialdb::SpatialDB* const db,
- const int querySize,
- const PylithScalar scale)
+ spatialdata::spatialdb::SpatialDB* const db,
+ const int querySize,
+ const PylithScalar scale)
{ // _queryDB
assert(name);
assert(db);
@@ -244,17 +244,14 @@
err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
- PetscSection parametersSection = _parameters->get(name).petscSection();assert(parametersSection);
- PetscVec parametersVec = _parameters->get(name).localVector();assert(parametersVec);
- PetscScalar *parametersArray = NULL;
- err = VecGetArray(parametersVec, ¶metersArray);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::Mesh>& parametersField = _parameters->get(name);
+ PetscScalar* parametersArray = parametersField.getLocalArray();
scalar_array valueVertex(querySize);
const int numPoints = _points.size();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- PetscInt cdof, coff;
-
// Get dimensionalized coordinates of vertex
+ PetscInt cdof, coff;
err = PetscSectionGetDof(coordSection, _points[iPoint], &cdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(coordSection, _points[iPoint], &coff);CHECK_PETSC_ERROR(err);
assert(cdof == spaceDim);
@@ -274,15 +271,12 @@
_getNormalizer().nondimensionalize(&valueVertex[0], valueVertex.size(), scale);
// Update section
- PetscInt dof, off;
-
- err = PetscSectionGetDof(parametersSection, _points[iPoint], &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(parametersSection, _points[iPoint], &off);CHECK_PETSC_ERROR(err);
- //assert(parametersFiberDim == dof);
+ const PetscInt off = parametersField.sectionOffset(_points[iPoint]);
+ const PetscInt dof = parametersField.sectionDof(_points[iPoint]);
for(int i = 0; i < dof; ++i)
parametersArray[off+i] = valueVertex[i];
} // for
- err = VecRestoreArray(parametersVec, ¶metersArray);CHECK_PETSC_ERROR(err);
+ parametersField.restoreLocalArray(¶metersArray);
err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
} // _queryDB
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -990,7 +990,7 @@
err = DMPlexGetLabel(dm, labelName, &label);CHECK_PETSC_ERROR(err);
// Completes the set of cells scheduled to be replaced
// Have to do internal fault vertices before fault boundary vertices, and this is the only thing I use faultBoundary for
- err = DMLabelCohesiveComplete(dm, label);CHECK_PETSC_ERROR(err);
+ err = DMPlexLabelCohesiveComplete(dm, label);CHECK_PETSC_ERROR(err);
err = DMPlexConstructCohesiveCells(dm, label, &sdm);CHECK_PETSC_ERROR(err);
mesh->setDMMesh(sdm);
} // createInterpolated
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.icc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/CoordsVisitor.icc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -58,6 +58,7 @@
// ----------------------------------------------------------------------
// Clear cached data.
+inline
void
pylith::topology::CoordsVisitor::clear(void)
{ // clear
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2013-03-13 20:29:41 UTC (rev 21527)
@@ -203,6 +203,8 @@
/** Get fiber dimension for point.
*
+ * @preq Must call cachePetscSection().
+ *
* @param point Point in mesh.
* @returns Fiber dimension.
*/
@@ -210,6 +212,8 @@
/** Get offset into array for point.
*
+ * @preq Must call cachePetscSection().
+ *
* @param point Point in mesh.
* @returns Offset.
*/
@@ -457,7 +461,7 @@
/// Data structures used in scattering to/from PETSc Vecs.
struct ScatterInfo {
- DM dm; ///< PETSc DM defining the communication pattern
+ PetscDM dm; ///< PETSc DM defining the communication pattern
PetscVec vector; ///< PETSc vector associated with field.
// Deprecated
PetscVecScatter scatter; ///< PETSc scatter associated with field.
@@ -506,14 +510,16 @@
private :
map_type _metadata;
+
/* Old construction */
const mesh_type& _mesh; ///< Mesh associated with section.
scatter_map_type _scatters; ///< Collection of scatters.
+
/* New construction */
- DM _dm; /* Holds the PetscSection */
- Vec _globalVec;
- Vec _localVec;
- std::map<std::string, int> _tmpFields;
+ PetscDM _dm; ///< Manages the PetscSection
+ PetscVec _globalVec; ///< Global PETSc vector
+ PetscVec _localVec; ///< Local PETSc vector
+ std::map<std::string, int> _tmpFields; ///< Map of fields to bundle together.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -79,7 +79,7 @@
PetscInt
pylith::topology::Field<mesh_type>::sectionDof(const PetscInt point) const
{ // sectionDof
- const PetscSection& section = petscSection();assert(section);
+ PetscSection section = petscSection();assert(section);
PetscInt dof;
PetscErrorCode err = PetscSectionGetDof(section, point, &dof);CHECK_PETSC_ERROR(err);
return dof;
@@ -91,7 +91,7 @@
PetscInt
pylith::topology::Field<mesh_type>::sectionOffset(const PetscInt point) const
{ // sectionOffset
- const PetscSection& section = petscSection();assert(section);
+ PetscSection section = petscSection();assert(section);
PetscInt offset;
PetscErrorCode err = PetscSectionGetOffset(section, point, &offset);CHECK_PETSC_ERROR(err);
return offset;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Makefile.am 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Makefile.am 2013-03-13 20:29:41 UTC (rev 21527)
@@ -36,8 +36,6 @@
Mesh.hh \
Mesh.icc \
MeshOps.hh \
- MeshVisitor.hh \
- MeshVisitor.icc \
RefineUniform.hh \
ReverseCuthillMcKee.hh \
SolutionFields.hh \
@@ -46,8 +44,10 @@
SubMesh.hh \
SubMesh.icc \
SubMesh.cc \
- SubMeshVisitor.hh \
- SubMeshVisitor.icc \
+ VisitorMesh.hh \
+ VisitorMesh.icc \
+ VisitorSubMesh.hh \
+ VisitorSubMesh.icc \
ISectionSpaces.hh \
ISectionSpaces.cc \
topologyfwd.hh
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.hh 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.hh 2013-03-13 20:29:41 UTC (rev 21527)
@@ -32,6 +32,7 @@
#include "ISectionSpaces.hh" // USES ISectionSpaces
+#include "pylith/utils/petscfwd.h" // HASA PetscDM
#include "pylith/utils/sievetypes.hh" // HASA pylith::SieveMesh
// Mesh -----------------------------------------------------------------
@@ -117,13 +118,13 @@
*
* @returns DMPlex mesh.
*/
- DM dmMesh(void) const;
+ PetscDM dmMesh(void) const;
/** Set DMPlex mesh.
*
* @param DMPlex mesh.
*/
- void setDMMesh(DM dm);
+ void setDMMesh(PetscDM dm);
/** Get sizes for all point types.
*
@@ -211,6 +212,8 @@
*/
int commRank(void) const;
+ /** Get
+
/** Print mesh to stdout.
*
* @param label Label for mesh.
@@ -241,11 +244,13 @@
private :
ALE::Obj<SieveMesh> _mesh; ///< Sieve mesh.
- DM _newMesh;
+ PetscDM _newMesh;
+
/* The old-style point numbering: [normalCells, normalVertices, shadowVertices, lagrangeVertices, cohesiveCells]
The new-style point numbering: [normalCells, cohesiveCells, normalVertices, shadowVertices, lagrangeVertices]
*/
PetscInt _numNormalCells, _numCohesiveCells, _numNormalVertices, _numShadowVertices, _numLagrangeVertices;
+
spatialdata::geocoords::CoordSys* _coordsys; ///< Coordinate system.
MPI_Comm _comm; ///< MPI communicator for mesh.
bool _debug; ///< Debugging flag for mesh.
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Mesh.icc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -179,7 +179,6 @@
PetscErrorCode err = DMView(_newMesh, PETSC_VIEWER_STDOUT_WORLD);CHECK_PETSC_ERROR(err);
}
-
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Stratum.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Stratum.hh 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Stratum.hh 2013-03-13 20:29:41 UTC (rev 21527)
@@ -88,7 +88,7 @@
PetscInt _begin; ///< Starting point.
PetscInt _end; ///< End point.
-}; // HeightStratum
+}; // Stratum
// StratumIS ------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -20,12 +20,13 @@
#include "SubMesh.hh" // implementation of class methods
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "pylith/utils/petscfwd.h" // USES PetscVec
+#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
#include <Selection.hh> // USES ALE::Selection
-#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
-
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
#include <cassert> // USES assert()
@@ -196,5 +197,4 @@
_coordsys->initialize();
} // initialize
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh 2013-03-13 20:29:41 UTC (rev 21527)
@@ -100,7 +100,7 @@
*
* @returns DMPlex mesh.
*/
- DM dmMesh(void) const;
+ PetscDM dmMesh(void) const;
/** Set DMPlex mesh.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.icc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/VisitorMesh.icc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -20,6 +20,10 @@
#error "VisitorMesh.icc must be included only from VisitorMesh.hh"
#else
+#include "Mesh.hh" // USES Mesh
+#include "SubMesh.hh" // USES SubMesh
+#include "Field.hh" // USES Field
+
#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -232,11 +232,11 @@
bc.setConstraints(field);
PetscSection fieldSection = field.petscSection();
- PetscVec fieldVec = field.localVector();
+ PetscVec fieldVec = field.localVector();
CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
const PetscInt numCells = cEnd - cStart;
- const PetscInt offset = numCells;
+ const PetscInt offset = numCells;
int iConstraint = 0;
for(PetscInt v = vStart; v < vEnd; ++v) {
const PetscInt *cInd, *fcInd;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2013-03-13 18:23:22 UTC (rev 21526)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2013-03-13 20:29:41 UTC (rev 21527)
@@ -330,8 +330,7 @@
// Check change start time.
const topology::Field<topology::SubMesh>& changeTime = bc._parameters->get("change time");
- _TestNeumann::_checkValues(_TestNeumann::changeTime,
- numQuadPts, changeTime);
+ _TestNeumann::_checkValues(_TestNeumann::changeTime, numQuadPts, changeTime);
th.close();
} // test_queryDatabases
@@ -441,7 +440,7 @@
const PylithScalar tolerance = 1.0e-06;
const int spaceDim = _TestNeumann::spaceDim;
const int numQuadPts = _TestNeumann::numQuadPts;
- CPPUNIT_ASSERT(0 != bc._parameters);
+ CPPUNIT_ASSERT(bc._parameters);
// Check values.
const topology::Field<topology::SubMesh>& value = bc._parameters->get("value");
@@ -717,8 +716,8 @@
// Check values in section against expected values.
void
pylith::bc::_TestNeumann::_checkValues(const PylithScalar* valuesE,
- const int fiberDimE,
- const topology::Field<topology::SubMesh>& field)
+ const int fiberDimE,
+ const topology::Field<topology::SubMesh>& field)
{ // _checkValues
CPPUNIT_ASSERT(valuesE);
@@ -746,8 +745,13 @@
err = PetscSectionGetDof(fieldSection, c, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(fieldSection, c, &off);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
- for (int iDim=0; iDim < fiberDimE; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[icell*fiberDimE+iDim]/scale, fieldArray[off+iDim], tolerance);
+ for (int iDim=0; iDim < fiberDimE; ++iDim) {
+ if (valuesE[icell*fiberDimE+iDim] != 0.0) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, fieldArray[off+iDim]/valuesE[icell*fiberDimE+iDim]*scale, tolerance);
+ } else {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[icell*fiberDimE+iDim], fieldArray[off+iDim]*scale, tolerance);
+ } // if/else
+ } // for
} // for
err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
} // _checkValues
More information about the CIG-COMMITS
mailing list