[cig-commits] r21467 - in short/3D/PyLith/trunk: libsrc/pylith/bc libsrc/pylith/utils unittests/libtests/bc unittests/libtests/bc/data
brad at geodynamics.org
brad at geodynamics.org
Thu Mar 7 15:38:15 PST 2013
Author: brad
Date: 2013-03-07 15:38:14 -0800 (Thu, 07 Mar 2013)
New Revision: 21467
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.icc
short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc
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/TimeDependent.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
short/3D/PyLith/trunk/libsrc/pylith/utils/petscfwd.h
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletDataMulti.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletDataMulti.hh
Log:
Added use of scales in Neumann and PointForce unit tests. Code cleanup.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -82,11 +82,10 @@
const int numCorners = _quadrature->numBasis();
// Get 'surface' cells (1 dimension lower than top-level cells)
- PetscDM subMesh = _boundaryMesh->dmMesh();
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
+
PetscInt cStart, cEnd;
PetscErrorCode err;
-
- assert(subMesh);
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
// Get damping constants at each quadrature point and rotate to
@@ -137,8 +136,8 @@
scalar_array dampingConstsLocal(spaceDim);
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- PetscVec coordVec;
+ 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);
@@ -151,32 +150,25 @@
PetscSection valueSection = _parameters->get("damping constants").petscSection();
PetscVec valueVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingConstsArray;
+ PetscScalar *dampingConstsArray = NULL;
assert(valueSection);assert(valueVec);
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
// Compute quadrature information
_quadrature->initializeGeometry();
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->computeGeometry(*_boundaryMesh, cells);
-#endif
err = VecGetArray(valueVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
-#else
const PetscScalar *coords;
PetscInt coordsSize;
err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
- for (PetscInt i = 0; i < coordsSize; ++i) {
+ for (PetscInt i=0; i < coordsSize; ++i) {
coordinatesCell[i] = coords[i];
} // for
_quadrature->computeGeometry(coordinatesCell, c);
err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
-#endif
PetscInt ddof, doff;
err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
@@ -210,8 +202,9 @@
const PylithScalar constTangential = densityN * vsN;
const PylithScalar constNormal = densityN * vpN;
const int numTangential = spaceDim-1;
- for (int iDim=0; iDim < numTangential; ++iDim)
+ for (int iDim=0; iDim < numTangential; ++iDim) {
dampingConstsLocal[iDim] = constTangential;
+ } // for
dampingConstsLocal[spaceDim-1] = constNormal;
// Compute normal/tangential orientation
@@ -223,8 +216,9 @@
for (int iDim=0; iDim < spaceDim; ++iDim) {
dampingConstsArray[doff+iQuad*spaceDim+iDim] = 0.0;
- for (int jDim=0; jDim < spaceDim; ++jDim)
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
dampingConstsArray[doff+iQuad*spaceDim+iDim] += dampingConstsLocal[jDim]*orientation[jDim*spaceDim+iDim];
+ } // for
// Ensure damping constants are positive
dampingConstsArray[doff+iQuad*spaceDim+iDim] = fabs(dampingConstsArray[doff+iQuad*spaceDim+iDim]);
} // for
@@ -238,10 +232,9 @@
// ----------------------------------------------------------------------
// Integrate contributions to residual term (r) for operator.
void
-pylith::bc::AbsorbingDampers::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::bc::AbsorbingDampers::integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateResidual
assert(_quadrature);
assert(_boundaryMesh);
@@ -267,38 +260,37 @@
// Allocate vectors for cell values.
_initCellVector();
+
// Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();
- PetscIS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
+ PetscIS subpointIS = NULL;
PetscInt cStart, cEnd;
PetscErrorCode err;
-
- assert(subMesh);
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
// Get sections
PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();
PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingConstsArray;
+ PetscScalar *dampingConstsArray = NULL;
assert(dampingConstsSection);assert(dampingConstsVec);
// Use _cellVector for cell residual.
- PetscSection residualSection = residual.petscSection(), residualSubsection;
- PetscVec residualVec = residual.localVector();
- assert(residualSection);assert(residualVec);
+ PetscSection residualSection = residual.petscSection();assert(residualSection);
+ PetscVec residualVec = residual.localVector();assert(residualVec);
+ PetscSection residualSubsection = NULL;
err = PetscSectionCreateSubmeshSection(residualSection, subpointIS, &residualSubsection);
- PetscSection velSection = fields->get("velocity(t)").petscSection(), velSubsection;
- PetscVec velVec = fields->get("velocity(t)").localVector();
- assert(velSection);assert(velVec);
+ 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);
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- PetscVec coordVec;
+ 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);
@@ -316,12 +308,14 @@
_logger->eventBegin(geometryEvent);
#endif
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
const PetscScalar *coordsArray;
- PetscInt coordsSize;
+ 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];}
+ for (PetscInt i=0; i < coordsSize; ++i) {
+ coordinatesCell[i] = coordsArray[i];
+ } // for
_quadrature->computeGeometry(coordinatesCell, c);
err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
@@ -395,10 +389,9 @@
// ----------------------------------------------------------------------
// Integrate contributions to residual term (r) for operator.
void
-pylith::bc::AbsorbingDampers::integrateResidualLumped(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::bc::AbsorbingDampers::integrateResidualLumped(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateResidualLumped
assert(_quadrature);
assert(_boundaryMesh);
@@ -425,37 +418,35 @@
_initCellVector();
// Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();
- PetscIS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
+ PetscIS subpointIS = NULL;
PetscInt cStart, cEnd;
- PetscErrorCode err;
-
- assert(subMesh);
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
// Get sections
PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();
PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingConstsArray;
+ PetscScalar *dampingConstsArray = NULL;
assert(dampingConstsSection);assert(dampingConstsVec);
// Use _cellVector for cell values.
- PetscSection residualSection = residual.petscSection(), residualSubsection;
- PetscVec residualVec = residual.localVector();
- assert(residualSection);assert(residualVec);
+ PetscSection residualSection = residual.petscSection();assert(residualSection);
+ PetscVec residualVec = residual.localVector();assert(residualVec);
+ PetscSection residualSubsection = NULL;
err = PetscSectionCreateSubmeshSection(residualSection, subpointIS, &residualSubsection);
- PetscSection velSection = fields->get("velocity(t)").petscSection(), velSubsection;
- PetscVec velVec = fields->get("velocity(t)").localVector();
- assert(velSection);assert(velVec);
+ 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);
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- PetscVec coordVec;
+ 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);
@@ -467,18 +458,20 @@
#endif
err = VecGetArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
- 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);
#endif
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
const PetscScalar *coordsArray;
- PetscInt coordsSize;
+ 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];}
+ for (PetscInt i=0; i < coordsSize; ++i) {
+ coordinatesCell[i] = coordsArray[i];
+ } // for
_quadrature->computeGeometry(coordinatesCell, c);
err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
@@ -495,7 +488,6 @@
const PetscScalar *velArray = NULL;
PetscInt velSize;
PetscInt ddof, doff;
-
err = PetscSectionGetDof(dampingConstsSection, c, &ddof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dampingConstsSection, c, &doff);CHECK_PETSC_ERROR(err);
assert(ddof == numQuadPts*spaceDim);
@@ -579,26 +571,24 @@
const int spaceDim = _quadrature->spaceDim();
// Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();
- PetscIS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
+ PetscIS subpointIS = NULL;
PetscInt cStart, cEnd;
- PetscErrorCode err = 0;
-
- assert(subMesh);
+ PetscErrorCode err;
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
// Get sections
- PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();
- PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingConstsArray;
- assert(dampingConstsSection);assert(dampingConstsVec);
+ PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();assert(dampingConstsSection);
+ PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();assert(dampingConstsVec);
+ PetscScalar *dampingConstsArray = NULL;
const topology::Field<topology::Mesh>& solution = fields->solution();
- PetscSection solutionSection = solution.petscSection(), solutionGlobalSection, solutionSubsection, solutionGlobalSubsection;
- PetscVec solutionVec = solution.localVector();
- PetscSF sf;
- assert(solutionSection);assert(solutionVec);
+ 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);
@@ -606,8 +596,7 @@
err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
// Get sparse matrix
- const PetscMat jacobianMat = jacobian->matrix();
- assert(jacobianMat);
+ const PetscMat jacobianMat = jacobian->matrix();assert(jacobianMat);
// Get parameters used in integration.
const PylithScalar dt = _dt;
@@ -618,8 +607,8 @@
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- PetscVec coordVec;
+ 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);
@@ -637,12 +626,14 @@
_logger->eventBegin(geometryEvent);
#endif
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+#error("Code for PRECOMPUTE_GEOMETRY not implemented")
#else
const PetscScalar *coordsArray;
- PetscInt coordsSize;
+ 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];}
+ for (PetscInt i = 0; i < coordsSize; ++i) {
+ coordinatesCell[i] = coordsArray[i];
+ } // for
_quadrature->computeGeometry(coordinatesCell, c);
err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
@@ -716,10 +707,9 @@
// ----------------------------------------------------------------------
// Integrate contributions to Jacobian matrix (A) associated with
void
-pylith::bc::AbsorbingDampers::integrateJacobian(
- topology::Field<topology::Mesh>* jacobian,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::bc::AbsorbingDampers::integrateJacobian(topology::Field<topology::Mesh>* jacobian,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateJacobian
assert(_quadrature);
assert(_boundaryMesh);
@@ -743,12 +733,10 @@
const int spaceDim = _quadrature->spaceDim();
// Get cell information
- PetscDM subMesh = _boundaryMesh->dmMesh();
- PetscIS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
+ PetscIS subpointIS = NULL;
PetscInt cStart, cEnd;
PetscErrorCode err;
-
- assert(subMesh);
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
@@ -761,21 +749,21 @@
_initCellVector();
// Get sections
- PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();
- PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();
+ PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();assert(dampingConstsSection);
+ PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();assert(dampingConstsVec);
PetscScalar *dampingConstsArray;
- assert(dampingConstsSection);assert(dampingConstsVec);
-
- PetscSection jacobianSection = jacobian->petscSection(), jacobianSubsection;
- PetscVec jacobianVec = jacobian->localVector();
- assert(jacobianSection);assert(jacobianVec);
+
+ 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);
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- PetscVec coordVec;
+ 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);
@@ -793,12 +781,14 @@
_logger->eventBegin(geometryEvent);
#endif
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+#error("Code for PRECOMPUTE_GEOMETRY not implemented");
#else
const PetscScalar *coordsArray;
- PetscInt coordsSize;
+ 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];}
+ for (PetscInt i = 0; i < coordsSize; ++i) {
+ coordinatesCell[i] = coordsArray[i];
+ } // for
_quadrature->computeGeometry(coordinatesCell, c);
err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
@@ -809,7 +799,6 @@
// Get damping constants
PetscInt ddof, doff;
-
err = PetscSectionGetDof(dampingConstsSection, c, &ddof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(dampingConstsSection, c, &doff);CHECK_PETSC_ERROR(err);
assert(ddof == numQuadPts*spaceDim);
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -61,7 +61,7 @@
delete _parameters; _parameters = 0;
_boundaryMesh = new topology::SubMesh(mesh, _label.c_str());
- assert(0 != _boundaryMesh);
+ assert(_boundaryMesh);
} // createSubMesh
// ----------------------------------------------------------------------
@@ -69,8 +69,8 @@
void
pylith::bc::BCIntegratorSubMesh::verifyConfiguration(const topology::Mesh& mesh) const
{ // verifyConfiguration
- assert(0 != _quadrature);
- assert(0 != _boundaryMesh);
+ assert(_quadrature);
+ assert(_boundaryMesh);
BoundaryCondition::verifyConfiguration(mesh);
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.icc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.icc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -33,7 +33,7 @@
inline
const pylith::topology::SubMesh&
pylith::bc::BCIntegratorSubMesh::boundaryMesh(void) const {
- assert(0 != _boundaryMesh);
+ assert(_boundaryMesh);
return *_boundaryMesh;
}
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -70,17 +70,15 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("BoundaryConditions");
- DM dmMesh = mesh.dmMesh();
- DMLabel label;
- IS pointIS;
+ PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
+ PetscDMLabel label = NULL;
+ PetscIS pointIS = NULL;
const PetscInt *points;
- PetscInt numPoints;
- PetscBool has;
- PetscErrorCode err;
-
- assert(dmMesh);
- err = DMPlexHasLabel(dmMesh, _label.c_str(), &has);CHECK_PETSC_ERROR(err);
- if (!has) {
+ PetscInt numPoints;
+ PetscBool hasLabel;
+ PetscErrorCode err;
+ err = DMPlexHasLabel(dmMesh, _label.c_str(), &hasLabel);CHECK_PETSC_ERROR(err);
+ if (!hasLabel) {
std::ostringstream msg;
msg << "Could not find group of points '" << _label << "' in mesh.";
throw std::runtime_error(msg.str());
@@ -94,6 +92,7 @@
for(PetscInt p = 0; p < numPoints; ++p) {_points[p] = points[p];}
err = ISRestoreIndices(pointIS, &points);CHECK_PETSC_ERROR(err);
err = ISDestroy(&pointIS);CHECK_PETSC_ERROR(err);
+
logger.stagePop();
} // _getPoints
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -68,7 +68,7 @@
_getPoints(mesh);
- assert(0 != _normalizer);
+ assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
_queryDatabases(mesh, lengthScale, "displacement");
} // initialize
@@ -82,15 +82,15 @@
if (0 == numFixedDOF)
return;
- PetscSection sectionP = field.petscSection();
+ PetscSection sectionP = field.petscSection();assert(sectionP);
PetscInt numFields;
PetscErrorCode err = 0;
- assert(sectionP);
// Set constraints in field
const int numPoints = _points.size();
_offsetLocal.resize(numPoints);
err = PetscSectionGetNumFields(sectionP, &numFields);CHECK_PETSC_ERROR(err);
+
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const PetscInt point = _points[iPoint];
PetscInt dof, cdof;
@@ -122,20 +122,19 @@
if (0 == numFixedDOF)
return;
- PetscSection sectionP = field.petscSection();
+ PetscSection sectionP = field.petscSection();assert(sectionP);
PetscInt numFields;
PetscErrorCode err = 0;
- assert(sectionP);
const int numPoints = _points.size();
err = PetscSectionGetNumFields(sectionP, &numFields);CHECK_PETSC_ERROR(err);
+
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
// Get list of currently constrained DOF
PetscInt cdof;
const PetscInt *cInd;
-
err = PetscSectionGetConstraintDof(sectionP, point, &cdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetConstraintIndices(sectionP, point, &cInd);CHECK_PETSC_ERROR(err);
@@ -191,34 +190,33 @@
_calculateValue(t);
const int numPoints = _points.size();
+ PetscErrorCode err = 0;
assert(_parameters);
- PetscSection valueSection = _parameters->get("value").petscSection();
- PetscVec valueVec = _parameters->get("value").localVector();
- PetscScalar *valueArray;
- assert(valueSection);assert(valueVec);
-
- PetscSection fieldSection = field.petscSection();
- PetscVec localVec = field.localVector();
- PetscScalar *array;
- PetscErrorCode err = 0;
- assert(fieldSection);assert(localVec);
-
- err = VecGetArray(localVec, &array);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);
+
+ PetscSection fieldSection = field.petscSection();assert(fieldSection);
+ PetscVec fieldVec = field.localVector();assert(fieldVec);
+ PetscScalar *fieldArray = NULL;
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
PetscInt p_bc = _points[iPoint];
- PetscInt dof, off, vdof, voff;
+ PetscInt dof, off, vdof, voff;
err = PetscSectionGetDof(fieldSection, p_bc, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(fieldSection, p_bc, &off);CHECK_PETSC_ERROR(err);
err = PetscSectionGetDof(valueSection, p_bc, &vdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(valueSection, p_bc, &voff);CHECK_PETSC_ERROR(err);
assert(vdof == numFixedDOF);
+
for(PetscInt iDOF = 0; iDOF < numFixedDOF; ++iDOF)
- array[_bcDOF[iDOF]+off] = valueArray[voff+iDOF];
+ fieldArray[_bcDOF[iDOF]+off] = valueArray[voff+iDOF];
} // for
- err = VecRestoreArray(localVec, &array);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
} // setField
@@ -237,21 +235,19 @@
_calculateValueIncr(t0, t1);
const int numPoints = _points.size();
+ PetscErrorCode err = 0;
assert(_parameters);
- PetscSection valueSection = _parameters->get("value").petscSection();
- PetscVec valueVec = _parameters->get("value").localVector();
- PetscScalar *valueArray;
- assert(valueSection);assert(valueVec);
-
- PetscSection fieldSection = field.petscSection();
- PetscVec localVec = field.localVector();
- PetscScalar *array;
- PetscErrorCode err = 0;
- assert(fieldSection);assert(localVec);
-
- err = VecGetArray(localVec, &array);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);
+
+ PetscSection fieldSection = field.petscSection();assert(fieldSection);
+ PetscVec fieldVec = field.localVector();assert(fieldVec);
+ PetscScalar *fieldArray = NULL;
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
PetscInt p_bc = _points[iPoint];
PetscInt dof, off, vdof, voff;
@@ -260,11 +256,12 @@
err = PetscSectionGetOffset(fieldSection, p_bc, &off);CHECK_PETSC_ERROR(err);
err = PetscSectionGetDof(valueSection, p_bc, &vdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(valueSection, p_bc, &voff);CHECK_PETSC_ERROR(err);
+
assert(vdof == numFixedDOF);
for(PetscInt iDOF=0; iDOF < numFixedDOF; ++iDOF)
- array[_bcDOF[iDOF]+off] = valueArray[voff+iDOF];
+ fieldArray[_bcDOF[iDOF]+off] = valueArray[voff+iDOF];
} // for
- err = VecRestoreArray(localVec, &array);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
} // setFieldIncr
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -67,7 +67,7 @@
DirichletBC::initialize(mesh, upDir);
_boundaryMesh = new topology::SubMesh(mesh, _label.c_str());
- assert(0 != _boundaryMesh);
+ assert(_boundaryMesh);
} // initialize
// ----------------------------------------------------------------------
@@ -83,7 +83,7 @@
assert(_boundaryMesh);
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
- assert(0 != cs);
+ assert(cs);
const int spaceDim = cs->spaceDim();
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
@@ -91,7 +91,7 @@
if (0 == _outputFields)
_outputFields = new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
- assert(0 != _outputFields);
+ 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);
@@ -149,36 +149,35 @@
} // if
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
- assert(0 != cs);
+ 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();
- Vec outputVec = _outputFields->get("buffer (vector)").localVector();
- PetscScalar *outputArray;
- PetscErrorCode err;
- assert(outputSection);assert(outputVec);
-
- PetscSection fieldSection = _parameters->get(name).petscSection();
- Vec fieldVec = _parameters->get(name).localVector();
- PetscScalar *fieldArray;
- assert(fieldSection);assert(fieldVec);
-
+ 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];
- PetscInt odof, ooff, fdof, foff;
+ 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);
+
for(PetscInt iDOF=0; iDOF < numFixedDOF; ++iDOF)
outputArray[ooff+_bcDOF[iDOF]] = fieldArray[foff+iDOF];
} // for
@@ -216,31 +215,30 @@
} // if
const int numPoints = _points.size();
+ PetscErrorCode err = 0;
assert(_outputFields->hasField("buffer (scalar)"));
- PetscSection outputSection = _outputFields->get("buffer (scalar)").petscSection();
- Vec outputVec = _outputFields->get("buffer (scalar)").localVector();
- PetscScalar *outputArray;
- PetscErrorCode err;
- assert(outputSection);assert(outputVec);
-
- PetscSection fieldSection = _parameters->get(name).petscSection();
- Vec fieldVec = _parameters->get(name).localVector();
- PetscScalar *fieldArray;
- assert(fieldSection);assert(fieldVec);
-
+ 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);
+
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
- PetscInt odof, ooff, fdof, foff;
+ 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);
+
outputArray[ooff] = fieldArray[foff];
} // for
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -37,11 +37,6 @@
//#define PRECOMPUTE_GEOMETRY
// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::RealUniformSection SubRealUniformSection;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::Neumann::Neumann(void)
{ // constructor
@@ -82,10 +77,9 @@
// ----------------------------------------------------------------------
// Integrate contributions to residual term (r) for operator.
void
-pylith::bc::Neumann::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::bc::Neumann::integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateResidual
assert(_quadrature);
assert(_boundaryMesh);
@@ -103,32 +97,29 @@
scalar_array tractionsCell(numQuadPts*spaceDim);
// Get cell information
- DM subMesh = _boundaryMesh->dmMesh();
- IS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
+ PetscIS subpointIS;
PetscInt cStart, cEnd;
- PetscErrorCode err;
-
- assert(subMesh);
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
// Get sections
_calculateValue(t);
- PetscSection valueSection = _parameters->get("value").petscSection();
- Vec valueVec = _parameters->get("value").localVector();
- PetscScalar *tractionsArray;
- assert(valueSection);assert(valueVec);
-
- PetscSection residualSection = residual.petscSection(), residualSubsection;
- Vec residualVec = residual.localVector();
- assert(residualSection);assert(residualVec);
+ PetscSection valueSection = _parameters->get("value").petscSection();assert(valueSection);
+ PetscVec valueVec = _parameters->get("value").localVector();assert(valueVec);
+ PetscScalar *tractionsArray = NULL;
+
+ 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);
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
+ 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);
@@ -138,12 +129,14 @@
err = VecGetArray(valueVec, &tractionsArray);CHECK_PETSC_ERROR(err);
for(PetscInt c = cStart; c < cEnd; ++c) {
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(c);
+#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
const PetscScalar *coords;
- PetscInt coordsSize;
+ 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];}
+ for (PetscInt i=0; i < coordsSize; ++i) {
+ coordinatesCell[i] = coords[i];
+ } // for
_quadrature->computeGeometry(coordinatesCell, c);
err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
@@ -247,34 +240,40 @@
// Create section to hold time dependent values
_parameters->add("value", "traction", topology::FieldBase::FACES_FIELD, numQuadPts*spaceDim);
- _parameters->get("value").vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- _parameters->get("value").scale(pressureScale);
- _parameters->get("value").allocate();
+ topology::Field<topology::SubMesh>& value = _parameters->get("value");
+ value.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ value.scale(pressureScale);
+ value.allocate();
if (_dbInitial) {
_parameters->add("initial", "initial_traction", topology::FieldBase::FACES_FIELD, numQuadPts*spaceDim);
- _parameters->get("initial").vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- _parameters->get("initial").scale(pressureScale);
- _parameters->get("initial").allocate();
+ topology::Field<topology::SubMesh>& initial = _parameters->get("initial");
+ initial.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ initial.scale(pressureScale);
+ initial.allocate();
}
if (_dbRate) {
_parameters->add("rate", "traction_rate", topology::FieldBase::FACES_FIELD, numQuadPts*spaceDim);
- _parameters->get("rate").vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- _parameters->get("rate").scale(rateScale);
- _parameters->get("rate").allocate();
+ topology::Field<topology::SubMesh>& rate = _parameters->get("rate");
+ rate.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ rate.scale(rateScale);
+ rate.allocate();
_parameters->add("rate time", "traction_rate_time", topology::FieldBase::FACES_FIELD, numQuadPts);
- _parameters->get("rate time").vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- _parameters->get("rate time").scale(timeScale);
- _parameters->get("rate time").allocate();
+ topology::Field<topology::SubMesh>& rateTime = _parameters->get("rate time");
+ rateTime.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ rateTime.scale(timeScale);
+ rateTime.allocate();
} // if
if (_dbChange) {
_parameters->add("change", "change_traction", topology::FieldBase::FACES_FIELD, numQuadPts*spaceDim);
- _parameters->get("change").vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- _parameters->get("change").scale(pressureScale);
- _parameters->get("change").allocate();
+ topology::Field<topology::SubMesh>& change = _parameters->get("change");
+ change.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ change.scale(pressureScale);
+ change.allocate();
_parameters->add("change time", "change_traction_time", topology::FieldBase::FACES_FIELD, numQuadPts);
- _parameters->get("change time").vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- _parameters->get("change time").scale(timeScale);
- _parameters->get("change time").allocate();
+ topology::Field<topology::SubMesh>& changeTime = _parameters->get("change time");
+ changeTime.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
+ changeTime.scale(timeScale);
+ changeTime.allocate();
} // if
if (0 != _dbInitial) { // Setup initial values, if provided.
@@ -396,11 +395,9 @@
assert(_parameters);
// Get 'surface' cells (1 dimension lower than top-level cells)
- DM subMesh = _boundaryMesh->dmMesh();
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
PetscInt cStart, cEnd;
- PetscErrorCode err;
-
- assert(subMesh);
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
@@ -415,16 +412,15 @@
// Get sections.
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
+ 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);
- PetscSection valueSection = _parameters->get(name).petscSection();
- Vec valueVec = _parameters->get(name).localVector();
- PetscScalar *valueArray;
- assert(valueSection);assert(valueVec);
+ PetscSection valueSection = _parameters->get(name).petscSection();assert(valueSection);
+ PetscVec valueVec = _parameters->get(name).localVector();assert(valueVec);
+ PetscScalar *valueArray = NULL;
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
assert(cs);
@@ -443,12 +439,14 @@
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
const PetscScalar *coords;
- PetscInt coordsSize;
+ 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];}
+ for (PetscInt i = 0; i < coordsSize; ++i) {
+ coordinatesCell[i] = coords[i];
+ } // for
_quadrature->computeGeometry(coordinatesCell, c);
err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
@@ -499,11 +497,9 @@
up[i] = upDir[i];
// Get 'surface' cells (1 dimension lower than top-level cells)
- DM subMesh = _boundaryMesh->dmMesh();
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
PetscInt cStart, cEnd;
- PetscErrorCode err;
-
- assert(subMesh);
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
// Quadrature related values.
@@ -524,43 +520,40 @@
// Get sections.
scalar_array coordinatesCell(numBasis*spaceDim);
- PetscSection coordSection;
- Vec coordVec;
- err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
- assert(coordSection);assert(coordVec);
+ 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);
- scalar_array parametersCellLocal(spaceDim);
- PetscSection initialSection, rateSection, rateTimeSection, changeSection, changeTimeSection;
- Vec initialVec, rateVec, rateTimeVec, changeVec, changeTimeVec;
- PetscScalar *valuesArray, *initialArray, *rateArray, *rateTimeArray, *changeArray, *changeTimeArray;
+ scalar_array parametersCellLocal(spaceDim);
+ PetscSection initialSection=NULL, rateSection=NULL, rateTimeSection=NULL, changeSection=NULL, changeTimeSection=NULL;
+ PetscVec initialVec=NULL, rateVec=NULL, rateTimeVec=NULL, changeVec=NULL, changeTimeVec=NULL;
+ PetscScalar *valuesArray=NULL, *initialArray=NULL, *rateArray=NULL, *rateTimeArray=NULL, *changeArray=NULL, *changeTimeArray=NULL;
- PetscSection valuesSection = _parameters->get("value").petscSection();
- Vec valuesVec = _parameters->get("value").localVector();
- assert(valuesSection);assert(valuesVec);
+ PetscSection valuesSection = _parameters->get("value").petscSection();assert(valuesSection);
+ PetscVec valuesVec = _parameters->get("value").localVector();assert(valuesVec);
err = VecGetArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
if (_dbInitial) {
- initialSection = _parameters->get("initial").petscSection();
- initialVec = _parameters->get("initial").localVector();
- assert(initialSection);assert(initialVec);
- err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ initialSection = _parameters->get("initial").petscSection();assert(initialSection);
+ initialVec = _parameters->get("initial").localVector();assert(initialVec);
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
}
if (_dbRate) {
- rateSection = _parameters->get("rate").petscSection();
- rateTimeSection = _parameters->get("rate time").petscSection();
- rateVec = _parameters->get("rate").localVector();
- rateTimeVec = _parameters->get("rate time").localVector();
- assert(rateSection);assert(rateTimeSection);assert(rateVec);assert(rateTimeVec);
- err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ rateSection = _parameters->get("rate").petscSection();assert(rateSection);
+ rateVec = _parameters->get("rate").localVector();assert(rateVec);
+ err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+
+ rateTimeSection = _parameters->get("rate time").petscSection();assert(rateTimeSection);
+ rateTimeVec = _parameters->get("rate time").localVector();assert(rateTimeVec);
+ err = VecGetArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
}
if (_dbChange) {
- changeSection = _parameters->get("change").petscSection();
- changeTimeSection = _parameters->get("change time").petscSection();
- changeVec = _parameters->get("change").localVector();
- changeTimeVec = _parameters->get("change time").localVector();
- assert(changeSection);assert(changeTimeSection);assert(changeVec);assert(changeTimeVec);
- err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ changeSection = _parameters->get("change").petscSection();assert(changeSection);
+ changeVec = _parameters->get("change").localVector();assert(changeVec);
+ err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+
+ changeTimeSection = _parameters->get("change time").petscSection();assert(changeTimeSection);
+ changeTimeVec = _parameters->get("change time").localVector();assert(changeTimeVec);
err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
}
@@ -569,12 +562,14 @@
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#else
const PetscScalar *coords;
- PetscInt coordsSize;
+ 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];}
+ for (PetscInt i=0; i < coordsSize; ++i) {
+ coordinatesCell[i] = coords[i];
+ } // for
_quadrature->computeGeometry(coordinatesCell, c);
err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
@@ -582,65 +577,67 @@
// Compute Jacobian and determinant at quadrature point, then get orientation.
memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(PylithScalar));
#if defined(PRECOMPUTE_GEOMETRY)
- coordsVisitor.clear();
- subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+#error("Code for PRECOMPUTE_GEOMETRY not implemented.")
#endif
cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell, quadPtRef);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
assert(jacobianDet > 0.0);
orientation /= jacobianDet;
- if (0 != _dbInitial) {
+ if (_dbInitial) {
// Rotate traction vector from local coordinate system to global coordinate system
PylithScalar *initialLocal = ¶metersCellLocal[0];
- PetscInt idof, ioff;
+ PetscInt idof, ioff;
err = PetscSectionGetDof(initialSection, c, &idof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(initialSection, c, &ioff);CHECK_PETSC_ERROR(err);
assert(idof == numQuadPts*spaceDim);
for(int iDim = 0; iDim < spaceDim; ++iDim) {
initialLocal[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)
+ for(int jDim = 0; jDim < spaceDim; ++jDim) {
initialArray[ioff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * initialLocal[jDim];
+ } // for
} // for
} // if
- if (0 != _dbRate) {
+ if (_dbRate) {
// Rotate traction vector from local coordinate system to global coordinate system
PylithScalar *rateLocal = ¶metersCellLocal[0];
- PetscInt rdof, roff;
+ PetscInt rdof, roff;
err = PetscSectionGetDof(rateSection, c, &rdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(rateSection, c, &roff);CHECK_PETSC_ERROR(err);
assert(rdof == numQuadPts*spaceDim);
for(int iDim = 0; iDim < spaceDim; ++iDim) {
rateLocal[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)
+ for(int jDim = 0; jDim < spaceDim; ++jDim) {
rateArray[roff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * rateLocal[jDim];
+ } // for
} // for
} // if
- if (0 != _dbChange) {
+ if (_dbChange) {
// Rotate traction vector from local coordinate system to global coordinate system
PylithScalar *changeLocal = ¶metersCellLocal[0];
- PetscInt cdof, coff;
+ PetscInt cdof, coff;
err = PetscSectionGetDof(changeSection, c, &cdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(changeSection, c, &coff);CHECK_PETSC_ERROR(err);
assert(cdof == numQuadPts*spaceDim);
- for(int iDim = 0; iDim < spaceDim; ++iDim) {
+ for (int iDim = 0; iDim < spaceDim; ++iDim) {
changeLocal[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)
+ for(int jDim = 0; jDim < spaceDim; ++jDim) {
changeArray[coff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * changeLocal[jDim];
+ } // for
} // for
} // if
} // for
@@ -648,15 +645,15 @@
err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
if (_dbInitial) {
err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
if (_dbRate) {
err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
if (_dbChange) {
err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
} // paramsLocalToGlobal
// ----------------------------------------------------------------------
@@ -671,72 +668,67 @@
const PylithScalar timeScale = _getNormalizer().timeScale();
// Get 'surface' cells (1 dimension lower than top-level cells)
- DM subMesh = _boundaryMesh->dmMesh();
+ PetscDM subMesh = _boundaryMesh->dmMesh();assert(subMesh);
PetscInt cStart, cEnd;
- PetscErrorCode err;
-
- assert(subMesh);
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
const int spaceDim = _quadrature->spaceDim();
const int numQuadPts = _quadrature->numQuadPts();
- PetscSection initialSection, rateSection, rateTimeSection, changeSection, changeTimeSection;
- Vec initialVec, rateVec, rateTimeVec, changeVec, changeTimeVec;
- PetscScalar *valuesArray, *initialArray, *rateArray, *rateTimeArray, *changeArray, *changeTimeArray;
+ PetscSection initialSection=NULL, rateSection=NULL, rateTimeSection=NULL, changeSection=NULL, changeTimeSection=NULL;
+ PetscVec initialVec=NULL, rateVec=NULL, rateTimeVec=NULL, changeVec=NULL, changeTimeVec=NULL;
+ PetscScalar *valuesArray=NULL, *initialArray=NULL, *rateArray=NULL, *rateTimeArray=NULL, *changeArray=NULL, *changeTimeArray=NULL;
- PetscSection valuesSection = _parameters->get("value").petscSection();
- Vec valuesVec = _parameters->get("value").localVector();
- assert(valuesSection);assert(valuesVec);
+ PetscSection valuesSection = _parameters->get("value").petscSection();assert(valuesSection);
+ PetscVec valuesVec = _parameters->get("value").localVector();assert(valuesVec);
err = VecGetArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
if (_dbInitial) {
- initialSection = _parameters->get("initial").petscSection();
- initialVec = _parameters->get("initial").localVector();
- assert(initialSection);assert(initialVec);
- err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- }
+ initialSection = _parameters->get("initial").petscSection();assert(initialSection);
+ initialVec = _parameters->get("initial").localVector();assert(initialVec);
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ } // if
if (_dbRate) {
- rateSection = _parameters->get("rate").petscSection();
- rateTimeSection = _parameters->get("rate time").petscSection();
- rateVec = _parameters->get("rate").localVector();
- rateTimeVec = _parameters->get("rate time").localVector();
- assert(rateSection);assert(rateTimeSection);assert(rateVec);assert(rateTimeVec);
- err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ rateSection = _parameters->get("rate").petscSection();assert(rateSection);
+ rateVec = _parameters->get("rate").localVector();assert(rateVec);
+ err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+
+ rateTimeSection = _parameters->get("rate time").petscSection();assert(rateTimeSection);
+ rateTimeVec = _parameters->get("rate time").localVector();assert(rateTimeVec);
err = VecGetArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
if (_dbChange) {
- changeSection = _parameters->get("change").petscSection();
- changeTimeSection = _parameters->get("change time").petscSection();
- changeVec = _parameters->get("change").localVector();
- changeTimeVec = _parameters->get("change time").localVector();
- assert(changeSection);assert(changeTimeSection);assert(changeVec);assert(changeTimeVec);
- err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ changeSection = _parameters->get("change").petscSection();assert(changeSection);
+ changeVec = _parameters->get("change").localVector();assert(changeVec);
+ err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+
+ changeTimeSection = _parameters->get("change time").petscSection();assert(changeTimeSection);
+ changeTimeVec = _parameters->get("change time").localVector();assert(changeTimeVec);
err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
for(PetscInt c = cStart; c < cEnd; ++c) {
PetscInt vdof, voff;
-
err = PetscSectionGetDof(valuesSection, c, &vdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(valuesSection, c, &voff);CHECK_PETSC_ERROR(err);
assert(vdof == numQuadPts*spaceDim);
- for(PetscInt d = 0; d < vdof; ++d)
+ for (PetscInt d = 0; d < vdof; ++d) {
valuesArray[voff+d] = 0.0;
+ } // for
// Contribution from initial value
- if (0 != _dbInitial) {
+ if (_dbInitial) {
PetscInt idof, ioff;
-
err = PetscSectionGetDof(initialSection, c, &idof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(initialSection, c, &ioff);CHECK_PETSC_ERROR(err);
assert(idof == vdof);
- for(PetscInt d = 0; d < idof; ++d)
+ for (PetscInt d = 0; d < idof; ++d) {
valuesArray[voff+d] += initialArray[ioff+d];
+ } // for
} // if
// Contribution from rate of change of value
- if (0 != _dbRate) {
+ if (_dbRate) {
PetscInt rdof, roff, rtdof, rtoff;
-
err = PetscSectionGetDof(rateSection, c, &rdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(rateSection, c, &roff);CHECK_PETSC_ERROR(err);
err = PetscSectionGetDof(rateTimeSection, c, &rtdof);CHECK_PETSC_ERROR(err);
@@ -747,15 +739,15 @@
for(int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
const PylithScalar tRel = t - rateTimeArray[rtoff+iQuad];
if (tRel > 0.0) // rate of change integrated over time
- for(int iDim = 0; iDim < spaceDim; ++iDim)
+ for(int iDim = 0; iDim < spaceDim; ++iDim) {
valuesArray[voff+iQuad*spaceDim+iDim] += rateArray[roff+iQuad*spaceDim+iDim] * tRel;
+ } // for
} // for
} // if
// Contribution from change of value
- if (0 != _dbChange) {
+ if (_dbChange) {
PetscInt cdof, coff, ctdof, ctoff;
-
err = PetscSectionGetDof(changeSection, c, &cdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(changeSection, c, &coff);CHECK_PETSC_ERROR(err);
err = PetscSectionGetDof(changeTimeSection, c, &ctdof);CHECK_PETSC_ERROR(err);
@@ -779,24 +771,25 @@
throw std::runtime_error(msg.str());
} // if
} // if
- for(int iDim = 0; iDim < spaceDim; ++iDim)
+ for (int iDim = 0; iDim < spaceDim; ++iDim) {
valuesArray[voff+iQuad*spaceDim+iDim] += changeArray[coff+iQuad*spaceDim+iDim] * scale;
+ } // for
} // if
} // for
} // if
} // for
- err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
if (_dbInitial) {
- err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- }
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ } // if
if (_dbRate) {
- err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
- }
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ } // if
if (_dbChange) {
- err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
} // _calculateValue
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -32,11 +32,6 @@
#include <sstream> // USES std::ostringstream
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RealUniformSection RealUniformSection;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::PointForce::PointForce(void)
{ // constructor
@@ -69,7 +64,7 @@
_getPoints(mesh);
- assert(0 != _normalizer);
+ assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
const PylithScalar pressureScale = _normalizer->pressureScale();
const PylithScalar forceScale = pressureScale * lengthScale * lengthScale;
@@ -80,10 +75,9 @@
// ----------------------------------------------------------------------
// Integrate contributions to residual term (r) for operator.
void
-pylith::bc::PointForce::integrateResidual(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
+pylith::bc::PointForce::integrateResidual(const topology::Field<topology::Mesh>& residual,
+ const PylithScalar t,
+ topology::SolutionFields* const fields)
{ // integrateResidualAssembled
assert(_parameters);
assert(_normalizer);
@@ -96,29 +90,25 @@
const topology::Mesh& mesh = residual.mesh();
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- assert(0 != cs);
+ assert(cs);
const int spaceDim = cs->spaceDim();
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
+ PetscSection residualSection = residual.petscSection();assert(residualSection);
+ PetscVec residualVec = residual.localVector();assert(residualVec);
PetscScalar *residualArray;
- assert(residualSection);assert(residualVec);
-
+
// Get global order
- DM dmMesh = residual.dmMesh();
- PetscSection globalSection;
+ PetscDM dmMesh = residual.dmMesh();assert(dmMesh);
+ PetscSection globalSection = NULL;
PetscErrorCode err;
-
- assert(dmMesh);
err = DMGetDefaultGlobalSection(dmMesh, &globalSection);CHECK_PETSC_ERROR(err);
- PetscSection valueSection = _parameters->get("value").petscSection();
- Vec valueVec = _parameters->get("value").localVector();
- PetscScalar *valueArray;
- assert(valueSection);assert(valueVec);
-
+ 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);
+
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const PetscInt p_bc = _points[iPoint]; // Get point label.
PetscInt goff;
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependent.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependent.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependent.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -60,7 +60,7 @@
void
pylith::bc::TimeDependent::verifyConfiguration(const topology::Mesh& mesh) const
{ // verifyConfiguration
- if (0 == _dbChange && 0 != _dbTimeHistory) {
+ if (!_dbChange && _dbTimeHistory) {
std::ostringstream msg;
msg << "Time dependent boundary condition '" << _getLabel() << "',\n has a "
<< "time history database but not change in value spatial database.";
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -62,7 +62,7 @@
const int size)
{ // bcDOF
if (size > 0)
- assert(0 != flags);
+ assert(flags);
_bcDOF.resize(size);
for (int i=0; i < size; ++i)
@@ -200,7 +200,7 @@
_queryDB("change time", _dbChange, 1, timeScale);
_dbChange->close();
- if (0 != _dbTimeHistory)
+ if (_dbTimeHistory)
_dbTimeHistory->open();
} // if
@@ -229,39 +229,28 @@
const topology::Mesh& mesh = _parameters->mesh();
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- assert(0 != cs);
+ assert(cs);
const int spaceDim = cs->spaceDim();
const PylithScalar lengthScale = _getNormalizer().lengthScale();
+ PetscErrorCode err = 0;
scalar_array coordsVertex(spaceDim);
- DM dmMesh = mesh.dmMesh();
- PetscSection coordSection;
- Vec coordVec;
- PetscScalar *coordArray;
- PetscErrorCode err;
-
- assert(dmMesh);
+ PetscDM dmMesh = mesh.dmMesh();assert(dmMesh);
+ PetscSection coordSection = NULL;
+ PetscVec coordVec = NULL;
+ PetscScalar *coordArray = NULL;
err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
- PetscSection parametersSection = _parameters->get(name).petscSection();
- Vec parametersVec = _parameters->get(name).localVector();
- assert(parametersSection);assert(parametersVec);
-#if 0
- const int parametersFiberDim = _parameters->fiberDim();
- const int valueIndex = _parameters->sectionIndex(name);
- const int valueFiberDim = _parameters->sectionFiberDim(name);
- assert(valueIndex+valueFiberDim <= parametersFiberDim);
- scalar_array parametersVertex(parametersFiberDim);
-#endif
+ 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);
scalar_array valueVertex(querySize);
-
const int numPoints = _points.size();
- PetscScalar *array;
- err = VecGetArray(parametersVec, &array);CHECK_PETSC_ERROR(err);
- err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
PetscInt cdof, coff;
@@ -269,7 +258,9 @@
err = PetscSectionGetDof(coordSection, _points[iPoint], &cdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(coordSection, _points[iPoint], &coff);CHECK_PETSC_ERROR(err);
assert(cdof == spaceDim);
- for (PetscInt d = 0; d < spaceDim; ++d) {coordsVertex[d] = coordArray[coff+d];}
+ for (PetscInt d = 0; d < spaceDim; ++d) {
+ coordsVertex[d] = coordArray[coff+d];
+ } // for
_getNormalizer().dimensionalize(&coordsVertex[0], coordsVertex.size(), lengthScale);
int err = db->query(&valueVertex[0], valueVertex.size(), &coordsVertex[0], coordsVertex.size(), cs);
if (err) {
@@ -289,9 +280,9 @@
err = PetscSectionGetOffset(parametersSection, _points[iPoint], &off);CHECK_PETSC_ERROR(err);
//assert(parametersFiberDim == dof);
for(int i = 0; i < dof; ++i)
- array[off+i] = valueVertex[i];
+ parametersArray[off+i] = valueVertex[i];
} // for
- err = VecRestoreArray(parametersVec, &array);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(parametersVec, ¶metersArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
} // _queryDB
@@ -305,39 +296,38 @@
const int numPoints = _points.size();
const int numBCDOF = _bcDOF.size();
const PylithScalar timeScale = _getNormalizer().timeScale();
- PetscSection initialSection, rateSection, rateTimeSection, changeSection, changeTimeSection;
- Vec initialVec, rateVec, rateTimeVec, changeVec, changeTimeVec;
- PetscScalar *valuesArray, *initialArray, *rateArray, *rateTimeArray, *changeArray, *changeTimeArray;
- PetscErrorCode err;
+ PetscSection initialSection=NULL, rateSection=NULL, rateTimeSection=NULL, changeSection=NULL, changeTimeSection=NULL;
+ PetscVec initialVec=NULL, rateVec=NULL, rateTimeVec=NULL, changeVec=NULL, changeTimeVec=NULL;
+ PetscScalar *valuesArray=NULL, *initialArray=NULL, *rateArray=NULL, *rateTimeArray=NULL, *changeArray=NULL, *changeTimeArray=NULL;
+ PetscErrorCode err = 0;
- PetscSection valuesSection = _parameters->get("value").petscSection();
- Vec valuesVec = _parameters->get("value").localVector();
- assert(valuesSection);assert(valuesVec);
+ PetscSection valuesSection = _parameters->get("value").petscSection();assert(valuesSection);
+ PetscVec valuesVec = _parameters->get("value").localVector();assert(valuesVec);
err = VecGetArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+
if (_dbInitial) {
- initialSection = _parameters->get("initial").petscSection();
- initialVec = _parameters->get("initial").localVector();
- assert(initialSection);assert(initialVec);
- err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- }
+ initialSection = _parameters->get("initial").petscSection();assert(initialSection);
+ initialVec = _parameters->get("initial").localVector();assert(initialVec);
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ } // if
if (_dbRate) {
- rateSection = _parameters->get("rate").petscSection();
- rateTimeSection = _parameters->get("rate time").petscSection();
- rateVec = _parameters->get("rate").localVector();
- rateTimeVec = _parameters->get("rate time").localVector();
- assert(rateSection);assert(rateTimeSection);assert(rateVec);assert(rateTimeVec);
+ rateSection = _parameters->get("rate").petscSection();assert(rateSection);
+ rateVec = _parameters->get("rate").localVector();assert(rateVec);
err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+
+ rateTimeSection = _parameters->get("rate time").petscSection();assert(rateTimeSection);
+ rateTimeVec = _parameters->get("rate time").localVector();assert(rateTimeVec);
err = VecGetArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
if (_dbChange) {
- changeSection = _parameters->get("change").petscSection();
- changeTimeSection = _parameters->get("change time").petscSection();
- changeVec = _parameters->get("change").localVector();
- changeTimeVec = _parameters->get("change time").localVector();
- assert(changeSection);assert(changeTimeSection);assert(changeVec);assert(changeTimeVec);
- err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ changeSection = _parameters->get("change").petscSection();assert(changeSection);
+ changeVec = _parameters->get("change").localVector();assert(changeVec);
+ err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+
+ changeTimeSection = _parameters->get("change time").petscSection();assert(changeTimeSection);
+ changeTimeVec = _parameters->get("change time").localVector();assert(changeTimeVec);
err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
for(int iPoint=0; iPoint < numPoints; ++iPoint) {
const int p_bc = _points[iPoint]; // Get point label
@@ -385,11 +375,11 @@
const PylithScalar tRel = t - changeTimeArray[ctoff];
if (tRel >= 0) { // change in value over time
PylithScalar scale = 1.0;
- if (0 != _dbTimeHistory) {
+ if (_dbTimeHistory) {
PylithScalar tDim = tRel;
_getNormalizer().dimensionalize(&tDim, 1, timeScale);
const int err = _dbTimeHistory->query(&scale, tDim);
- if (0 != err) {
+ if (err) {
std::ostringstream msg;
msg << "Error querying for time '" << tDim
<< "' in time history database "
@@ -428,39 +418,38 @@
const int numPoints = _points.size();
const int numBCDOF = _bcDOF.size();
const PylithScalar timeScale = _getNormalizer().timeScale();
- PetscSection initialSection, rateSection, rateTimeSection, changeSection, changeTimeSection;
- Vec initialVec, rateVec, rateTimeVec, changeVec, changeTimeVec;
- PetscScalar *valuesArray, *initialArray, *rateArray, *rateTimeArray, *changeArray, *changeTimeArray;
- PetscErrorCode err;
+ PetscSection initialSection=NULL, rateSection=NULL, rateTimeSection=NULL, changeSection=NULL, changeTimeSection=NULL;
+ PetscVec initialVec=NULL, rateVec=NULL, rateTimeVec=NULL, changeVec=NULL, changeTimeVec=NULL;
+ PetscScalar *valuesArray=NULL, *initialArray=NULL, *rateArray=NULL, *rateTimeArray=NULL, *changeArray=NULL, *changeTimeArray=NULL;
+ PetscErrorCode err = 0;
- PetscSection valuesSection = _parameters->get("value").petscSection();
- Vec valuesVec = _parameters->get("value").localVector();
- assert(valuesSection);assert(valuesVec);
+ PetscSection valuesSection = _parameters->get("value").petscSection();assert(valuesSection);
+ PetscVec valuesVec = _parameters->get("value").localVector();assert(valuesVec);
err = VecGetArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+
if (_dbInitial) {
- initialSection = _parameters->get("initial").petscSection();
- initialVec = _parameters->get("initial").localVector();
- assert(initialSection);assert(initialVec);
- err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- }
+ initialSection = _parameters->get("initial").petscSection();assert(initialSection);
+ initialVec = _parameters->get("initial").localVector();assert(initialVec);
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ } // if
if (_dbRate) {
- rateSection = _parameters->get("rate").petscSection();
- rateTimeSection = _parameters->get("rate time").petscSection();
- rateVec = _parameters->get("rate").localVector();
- rateTimeVec = _parameters->get("rate time").localVector();
- assert(rateSection);assert(rateTimeSection);assert(rateVec);assert(rateTimeVec);
+ rateSection = _parameters->get("rate").petscSection();assert(rateSection);
+ rateVec = _parameters->get("rate").localVector();assert(rateVec);
err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+
+ rateTimeSection = _parameters->get("rate time").petscSection();assert(rateTimeSection);
+ rateTimeVec = _parameters->get("rate time").localVector();assert(rateTimeVec);
err = VecGetArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
if (_dbChange) {
- changeSection = _parameters->get("change").petscSection();
- changeTimeSection = _parameters->get("change time").petscSection();
- changeVec = _parameters->get("change").localVector();
- changeTimeVec = _parameters->get("change time").localVector();
- assert(changeSection);assert(changeTimeSection);assert(changeVec);assert(changeTimeVec);
- err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ changeSection = _parameters->get("change").petscSection();assert(changeSection);
+ changeVec = _parameters->get("change").localVector();assert(changeVec);
+ err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+
+ changeTimeSection = _parameters->get("change time").petscSection();assert(changeTimeSection);
+ changeTimeVec = _parameters->get("change time").localVector();assert(changeTimeVec);
err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const int p_bc = _points[iPoint]; // Get point label
@@ -509,11 +498,11 @@
if (t0 >= tChange) { // increment is after change starts
PylithScalar scale0 = 1.0;
PylithScalar scale1 = 1.0;
- if (0 != _dbTimeHistory) {
+ if (_dbTimeHistory) {
PylithScalar tDim = t0 - tChange;
_getNormalizer().dimensionalize(&tDim, 1, timeScale);
int err = _dbTimeHistory->query(&scale0, tDim);
- if (0 != err) {
+ if (err) {
std::ostringstream msg;
msg << "Error querying for time '" << tDim
<< "' in time history database "
@@ -523,7 +512,7 @@
tDim = t1 - tChange;
_getNormalizer().dimensionalize(&tDim, 1, timeScale);
err = _dbTimeHistory->query(&scale1, tDim);
- if (0 != err) {
+ if (err) {
std::ostringstream msg;
msg << "Error querying for time '" << tDim
<< "' in time history database "
@@ -535,11 +524,11 @@
valuesArray[voff+d] += changeArray[coff+d] * (scale1 - scale0);
} else if (t1 >= tChange) { // increment spans when change starts
PylithScalar scale1 = 1.0;
- if (0 != _dbTimeHistory) {
+ if (_dbTimeHistory) {
PylithScalar tDim = t1 - tChange;
_getNormalizer().dimensionalize(&tDim, 1, timeScale);
int err = _dbTimeHistory->query(&scale1, tDim);
- if (0 != err) {
+ if (err) {
std::ostringstream msg;
msg << "Error querying for time '" << tDim
<< "' in time history database "
@@ -552,18 +541,18 @@
} // if/else
} // if
} // for
- err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
if (_dbInitial) {
- err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- }
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ } // if
if (_dbRate) {
- err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
- }
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ } // if
if (_dbChange) {
- err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
- }
+ } // if
} // _calculateValueIncr
Modified: short/3D/PyLith/trunk/libsrc/pylith/utils/petscfwd.h
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/utils/petscfwd.h 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/libsrc/pylith/utils/petscfwd.h 2013-03-07 23:38:14 UTC (rev 21467)
@@ -9,7 +9,7 @@
// This code was developed as part of the Computational Infrastructure
// for Geodynamics (http://geodynamics.org).
//
-// Copyright (c) 2010-2011 University of California, Davis
+// Copyright (c) 2010-2013 University of California, Davis
//
// See COPYING for license information.
//
@@ -49,10 +49,13 @@
/// forward declaration for PETSc PC
typedef struct _p_PC* PetscPC;
-/// forward declaration for PETSc Mat
+/// forward declaration for PETSc DM
typedef struct _p_DM* PetscDM;
-/// forward declaration for PETSc Mat
+/// forward declaration for PETSc DMLabel
+typedef struct _n_DMLabel* PetscDMLabel;
+
+/// forward declaration for PETSc IS
typedef struct _p_IS* PetscIS;
/// forward declaration for PETSc ISLocalToGlobalMapping
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -390,14 +390,15 @@
// Set coordinate system
spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(mesh->dimension());
+ cs.initialize();
+ mesh->coordsys(&cs);
+
spatialdata::units::Nondimensional normalizer;
normalizer.lengthScale(_data->lengthScale);
normalizer.pressureScale(_data->pressureScale);
normalizer.densityScale(_data->densityScale);
normalizer.timeScale(_data->timeScale);
- cs.setSpaceDim(mesh->dimension());
- cs.initialize();
- mesh->coordsys(&cs);
mesh->nondimensionalize(normalizer);
// Set up quadrature
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -315,7 +315,7 @@
err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
// Only unconstrained values should be zero.
- const PylithScalar t = 1.0;
+ const PylithScalar t = 1.0 / timeScale;
bc.setField(t, field);
// Create list of unconstrained DOF at constrained DOF
@@ -337,6 +337,8 @@
const PetscInt numFixedDOF = _data->numFixedDOF;
int iConstraint = 0;
+ const PylithScalar tRef = _data->tRef / timeScale;
+
err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
PetscInt dof, cdof, off;
@@ -358,8 +360,8 @@
// check constrained DOF
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
const int index = iConstraint * numFixedDOF + iDOF;
- const PylithScalar valueE = (t > _data->tRef/timeScale) ?
- _data->valuesInitial[index]/dispScale + (t-_data->tRef/timeScale)*_data->valueRate/velocityScale :
+ const PylithScalar valueE = (t > tRef) ?
+ _data->valuesInitial[index]/dispScale + (t-tRef)*_data->valueRate/velocityScale :
_data->valuesInitial[index]/dispScale;
CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[off+_data->fixedDOF[iDOF]], tolerance);
} // for
@@ -421,8 +423,8 @@
err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
// Only unconstrained values should be zero.
- const PylithScalar t0 = 1.0;
- const PylithScalar t1 = 2.0;
+ const PylithScalar t0 = 1.0 / timeScale;
+ const PylithScalar t1 = 2.0 / timeScale;
bc.setFieldIncr(t0, t1, field);
// Create list of unconstrained DOF at constrained DOF
@@ -444,6 +446,8 @@
const PetscInt numFixedDOF = _data->numFixedDOF;
int iConstraint = 0;
+ const PylithScalar tRef = _data->tRef / timeScale;
+
err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v) {
PetscInt dof, cdof, off;
@@ -464,8 +468,8 @@
// check constrained DOF
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
- const PylithScalar valueE = (t0 > _data->tRef/timeScale) ? (t1-t0)*_data->valueRate/velocityScale :
- (t1 > _data->tRef/timeScale) ? (t1-_data->tRef/timeScale)*_data->valueRate/velocityScale : 0.0;
+ const PylithScalar valueE = (t0 > tRef) ? (t1-t0)*_data->valueRate/velocityScale :
+ (t1 > tRef) ? (t1-tRef)*_data->valueRate/velocityScale : 0.0;
CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[off+_data->fixedDOF[iDOF]], tolerance);
} // for
++iConstraint;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -65,10 +65,9 @@
_initialize(&mesh, &bcA, &bcB, &bcC);
CPPUNIT_ASSERT(0 != _data);
- DM dmMesh = mesh.dmMesh();
- CPPUNIT_ASSERT(dmMesh);
- PetscInt cStart, cEnd, vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
@@ -84,12 +83,11 @@
bcC.setConstraintSizes(field);
field.allocate();
- PetscSection fieldSection = field.petscSection();
- Vec fieldVec = field.localVector();
- CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
+ PetscSection fieldSection = field.petscSection();CPPUNIT_ASSERT(fieldSection);
+ PetscVec fieldVec = field.localVector();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) {
PetscInt dof, cdof;
@@ -113,10 +111,9 @@
_initialize(&mesh, &bcA, &bcB, &bcC);
CPPUNIT_ASSERT(0 != _data);
- DM dmMesh = mesh.dmMesh();
- CPPUNIT_ASSERT(dmMesh);
- PetscInt cStart, cEnd, vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
@@ -135,18 +132,17 @@
bcB.setConstraints(field);
bcC.setConstraints(field);
- PetscSection fieldSection = field.petscSection();
- Vec fieldVec = field.localVector();
- CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
+ PetscSection fieldSection = field.petscSection();CPPUNIT_ASSERT(fieldSection);
+ PetscVec fieldVec = field.localVector();CPPUNIT_ASSERT(fieldVec);
const PetscInt numCells = cEnd - cStart;
- const PetscInt offset = numCells;
+ const PetscInt offset = numCells;
int index = 0;
int iConstraint = 0;
for(PetscInt v = vStart; v < vEnd; ++v) {
const int numConstrainedDOF = _data->constraintSizes[v-offset];
const PetscInt *cInd;
- PetscInt dof, cdof;
+ PetscInt dof, cdof;
err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
@@ -171,10 +167,9 @@
_initialize(&mesh, &bcA, &bcB, &bcC);
CPPUNIT_ASSERT(0 != _data);
- DM dmMesh = mesh.dmMesh();
- CPPUNIT_ASSERT(dmMesh);
- PetscInt cStart, cEnd, vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
@@ -193,10 +188,11 @@
bcB.setConstraints(field);
bcC.setConstraints(field);
- PetscSection fieldSection = field.petscSection();
- Vec fieldVec = field.localVector();
- CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
+ PetscSection fieldSection = field.petscSection();CPPUNIT_ASSERT(fieldSection);
+ PetscVec fieldVec = field.localVector();CPPUNIT_ASSERT(fieldVec);
+
const PylithScalar tolerance = 1.0e-06;
+ const PylithScalar valueScale = _data->lengthScale;
// All values should be zero.
PetscScalar *values;
@@ -214,7 +210,7 @@
// Only unconstrained values should be zero.
// Expected values set in _data->field
- const PylithScalar t = 10.0;
+ const PylithScalar t = 10.0/_data->timeScale;
bcA.setField(t, field);
bcB.setField(t, field);
bcC.setField(t, field);
@@ -227,7 +223,7 @@
err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
for(int iDOF = 0; iDOF < dof; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->field[i++], values[off+iDOF], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->field[i++], values[off+iDOF]*valueScale, tolerance);
} // for
err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
} // testSetField
@@ -244,7 +240,7 @@
_initialize(&mesh, &bcA, &bcB, &bcC);
CPPUNIT_ASSERT(0 != _data);
- DM dmMesh = mesh.dmMesh();
+ PetscDM dmMesh = mesh.dmMesh();
CPPUNIT_ASSERT(dmMesh);
PetscInt cStart, cEnd, vStart, vEnd;
PetscErrorCode err;
@@ -267,9 +263,11 @@
bcC.setConstraints(field);
PetscSection fieldSection = field.petscSection();
- Vec fieldVec = field.localVector();
+ PetscVec fieldVec = field.localVector();
CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
+
const PylithScalar tolerance = 1.0e-06;
+ const PylithScalar valueScale = _data->lengthScale;
// All values should be zero.
PetscScalar *values;
@@ -287,8 +285,8 @@
// Only unconstrained values should be zero.
// Expected values set in _data->field
- const PylithScalar t0 = 10.0;
- const PylithScalar t1 = 14.0;
+ const PylithScalar t0 = 10.0/_data->timeScale;
+ const PylithScalar t1 = 14.0/_data->timeScale;
bcA.setFieldIncr(t0, t1, field);
bcB.setFieldIncr(t0, t1, field);
bcC.setFieldIncr(t0, t1, field);
@@ -300,8 +298,9 @@
err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
- for(int iDOF = 0; iDOF < dof; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->fieldIncr[i++], values[off+iDOF], tolerance);
+ for(int iDOF = 0; iDOF < dof; ++iDOF) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->fieldIncr[i++], values[off+iDOF]*valueScale, tolerance);
+ } // for
} // for
err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
} // testSetFieldIncr
@@ -323,10 +322,15 @@
iohandler.read(mesh);
spatialdata::geocoords::CSCart cs;
- spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+
+ spatialdata::units::Nondimensional normalizer;
+ normalizer.lengthScale(_data->lengthScale);
+ normalizer.pressureScale(_data->pressureScale);
+ normalizer.densityScale(_data->densityScale);
+ normalizer.timeScale(_data->timeScale);
mesh->nondimensionalize(normalizer);
// Setup boundary condition A
@@ -346,6 +350,7 @@
bcA->dbInitial(&db);
bcA->dbRate(&dbRate);
bcA->bcDOF(_data->fixedDOFA, _data->numFixedDOFA);
+ bcA->normalizer(normalizer);
bcA->initialize(*mesh, upDir);
// Setup boundary condition B
@@ -359,6 +364,7 @@
bcB->dbInitial(&db);
bcB->dbRate(&dbRate);
bcB->bcDOF(_data->fixedDOFB, _data->numFixedDOFB);
+ bcB->normalizer(normalizer);
bcB->initialize(*mesh, upDir);
// Setup boundary condition C
@@ -373,6 +379,7 @@
bcC->dbInitial(&db);
bcC->dbRate(&dbRate);
bcC->bcDOF(_data->fixedDOFC, _data->numFixedDOFC);
+ bcC->normalizer(normalizer);
bcC->initialize(*mesh, upDir);
} // if
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -45,16 +45,9 @@
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::bc::TestNeumann );
// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-
-// ----------------------------------------------------------------------
namespace pylith {
namespace bc {
namespace _TestNeumann {
- const PylithScalar pressureScale = 4.0;
- const PylithScalar lengthScale = 1.0; // Mesh coordinates have scale=1.0
- const PylithScalar timeScale = 0.5;
const int ncells = 2;
const int numQuadPts = 2;
const int spaceDim = 2;
@@ -152,14 +145,12 @@
topology::SolutionFields fields(mesh);
_initialize(&mesh, &bc, &fields);
- CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(_data);
const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
- DM subMesh = boundaryMesh.dmMesh();
+ PetscDM subMesh = boundaryMesh.dmMesh();assert(subMesh);
PetscInt cStart, cEnd, vStart, vEnd;
- PetscErrorCode err;
-
- assert(subMesh);
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexGetDepthStratum(subMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
@@ -174,7 +165,7 @@
CPPUNIT_ASSERT_EQUAL(_data->numCells, numCells);
PetscInt dp = 0;
- for(PetscInt c = cStart; c < cEnd; ++c) {
+ for (PetscInt c = cStart; c < cEnd; ++c) {
PetscInt *closure = PETSC_NULL;
PetscInt closureSize, numCorners = 0;
@@ -183,12 +174,12 @@
const PetscInt point = closure[p];
if ((point >= vStart) && (point < vEnd)) {
closure[numCorners++] = point;
- }
- }
+ } // if
+ } // for
CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
for(PetscInt p = 0; p < numCorners; ++p, ++dp) {
CPPUNIT_ASSERT_EQUAL(_data->cells[dp], closure[p]);
- }
+ } // for
err = DMPlexRestoreTransitiveClosure(subMesh, c, PETSC_TRUE, &closureSize, &closure);CHECK_PETSC_ERROR(err);
} // for
@@ -199,13 +190,14 @@
PetscInt index = 0;
CPPUNIT_ASSERT(0 != bc._parameters);
PetscSection initialSection = bc._parameters->get("initial").petscSection();
- Vec initialVec = bc._parameters->get("initial").localVector();
+ PetscVec initialVec = bc._parameters->get("initial").localVector();
PetscScalar *initialArray;
CPPUNIT_ASSERT(initialSection);CPPUNIT_ASSERT(initialVec);
const PylithScalar tolerance = 1.0e-06;
+ const PylithScalar pressureScale = _data->pressureScale;
err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
- for(PetscInt c = cStart; c < cEnd; ++c) {
+ for (PetscInt c = cStart; c < cEnd; ++c) {
PetscInt dof, off;
err = PetscSectionGetDof(initialSection, c, &dof);CHECK_PETSC_ERROR(err);
@@ -213,8 +205,8 @@
CPPUNIT_ASSERT(dof == numQuadPts*spaceDim);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
for (int iDim =0; iDim < spaceDim; ++iDim) {
- const PylithScalar tractionsCellData = _data->tractionsCell[index];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionsCellData, initialArray[off+iQuad*spaceDim+iDim], tolerance);
+ const PylithScalar tractionE = _data->tractionsCell[index];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionE, initialArray[off+iQuad*spaceDim+iDim]*pressureScale, tolerance);
++index;
} // for
} // for
@@ -237,34 +229,35 @@
const PylithScalar t = 0.0;
bc.integrateResidual(residual, t, &fields);
- DM dmMesh = mesh.dmMesh();
- PetscInt vStart, vEnd;
- PetscErrorCode err;
+ PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ PetscErrorCode err = 0;
const PylithScalar* valsE = _data->valsResidual;
- CPPUNIT_ASSERT(dmMesh);
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int totalNumVertices = vEnd - vStart;
const int sizeE = _data->spaceDim * totalNumVertices;
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
+ PetscSection residualSection = residual.petscSection();CPPUNIT_ASSERT(residualSection);
+ PetscVec residualVec = residual.localVector();CPPUNIT_ASSERT(residualVec);
PetscScalar *vals;
- PetscInt size;
+ PetscInt size;
- CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
//residual.view("RESIDUAL");
const PylithScalar tolerance = 1.0e-06;
+ const PylithScalar residualScale = _data->pressureScale * pow(_data->lengthScale, _data->spaceDim-1);
+
err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
for (int i=0; i < size; ++i)
- if (fabs(valsE[i]) > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ if (fabs(valsE[i]) > 1.0) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i]*residualScale, tolerance);
+ } else {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i]*residualScale, tolerance);
+ } // if/else
err = VecRestoreArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateResidual
@@ -280,7 +273,7 @@
topology::Mesh mesh;
Neumann bc;
- _preinitialize(&mesh, &bc, true);
+ _preinitialize(&mesh, &bc);
spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDatabases");
spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
@@ -308,8 +301,8 @@
bc.dbChange(&dbChange);
bc.dbTimeHistory(&th);
- const PylithScalar pressureScale = _TestNeumann::pressureScale;
- const PylithScalar timeScale = _TestNeumann::timeScale;
+ const PylithScalar pressureScale = _data->pressureScale;
+ const PylithScalar timeScale = _data->timeScale;
bc._queryDatabases();
const PylithScalar tolerance = 1.0e-06;
@@ -320,31 +313,23 @@
// bc._parameters->view("PARAMETERS"); // DEBUGGING
// Check initial values.
- const topology::Field<topology::SubMesh>& initial =
- bc._parameters->get("initial");
- _TestNeumann::_checkValues(_TestNeumann::initial,
- numQuadPts*spaceDim, initial);
+ const topology::Field<topology::SubMesh>& initial = bc._parameters->get("initial");
+ _TestNeumann::_checkValues(_TestNeumann::initial, numQuadPts*spaceDim, initial);
// Check rate values.
const topology::Field<topology::SubMesh>& rate = bc._parameters->get("rate");
- _TestNeumann::_checkValues(_TestNeumann::rate,
- numQuadPts*spaceDim, rate);
+ _TestNeumann::_checkValues(_TestNeumann::rate, numQuadPts*spaceDim, rate);
// Check rate start time.
- const topology::Field<topology::SubMesh>& rateTime =
- bc._parameters->get("rate time");
- _TestNeumann::_checkValues(_TestNeumann::rateTime,
- numQuadPts, rateTime);
+ const topology::Field<topology::SubMesh>& rateTime = bc._parameters->get("rate time");
+ _TestNeumann::_checkValues(_TestNeumann::rateTime, numQuadPts, rateTime);
// Check change values.
- const topology::Field<topology::SubMesh>& change =
- bc._parameters->get("change");
- _TestNeumann::_checkValues(_TestNeumann::change,
- numQuadPts*spaceDim, change);
+ const topology::Field<topology::SubMesh>& change = bc._parameters->get("change");
+ _TestNeumann::_checkValues(_TestNeumann::change, numQuadPts*spaceDim, change);
// Check change start time.
- const topology::Field<topology::SubMesh>& changeTime =
- bc._parameters->get("change time");
+ const topology::Field<topology::SubMesh>& changeTime = bc._parameters->get("change time");
_TestNeumann::_checkValues(_TestNeumann::changeTime,
numQuadPts, changeTime);
th.close();
@@ -357,12 +342,12 @@
{ // test_paramsLocalToGlobal
_data = new NeumannDataQuad4();
feassemble::GeometryLine2D geometry;
- CPPUNIT_ASSERT(0 != _quadrature);
+ CPPUNIT_ASSERT(_quadrature);
_quadrature->refGeometry(&geometry);
topology::Mesh mesh;
Neumann bc;
- _preinitialize(&mesh, &bc, true);
+ _preinitialize(&mesh, &bc);
spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDatabases");
spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
@@ -386,8 +371,8 @@
bc.dbRate(&dbRate);
bc.dbChange(&dbChange);
- const PylithScalar pressureScale = _TestNeumann::pressureScale;
- const PylithScalar timeScale = _TestNeumann::timeScale;
+ const PylithScalar pressureScale = _data->pressureScale;
+ const PylithScalar timeScale = _data->timeScale;
bc._queryDatabases();
const PylithScalar upDir[3] = { 0.0, 0.0, 1.0 };
bc._paramsLocalToGlobal(upDir);
@@ -407,10 +392,8 @@
valuesE[i+0] = _TestNeumann::initial[i+0]; // x
valuesE[i+1] = -_TestNeumann::initial[i+1]; // y
} // for
- const topology::Field<topology::SubMesh>& initial =
- bc._parameters->get("initial");
- _TestNeumann::_checkValues(&valuesE[0],
- numQuadPts*spaceDim, initial);
+ const topology::Field<topology::SubMesh>& initial = bc._parameters->get("initial");
+ _TestNeumann::_checkValues(&valuesE[0], numQuadPts*spaceDim, initial);
// Check rate values.
for (int i=0; i < valuesE.size(); i+=spaceDim) {
@@ -418,18 +401,15 @@
valuesE[i+1] = -_TestNeumann::rate[i+1]; // y
} // for
const topology::Field<topology::SubMesh>& rate = bc._parameters->get("rate");
- _TestNeumann::_checkValues(&valuesE[0],
- numQuadPts*spaceDim, rate);
+ _TestNeumann::_checkValues(&valuesE[0], numQuadPts*spaceDim, rate);
// Check change values.
for (int i=0; i < valuesE.size(); i+=spaceDim) {
valuesE[i+0] = _TestNeumann::change[i+0]; // x
valuesE[i+1] = -_TestNeumann::change[i+1]; // y
} // for
- const topology::Field<topology::SubMesh>& change =
- bc._parameters->get("change");
- _TestNeumann::_checkValues(&valuesE[0],
- numQuadPts*spaceDim, change);
+ const topology::Field<topology::SubMesh>& change = bc._parameters->get("change");
+ _TestNeumann::_checkValues(&valuesE[0], numQuadPts*spaceDim, change);
} // test_paramsLocalToGlobal
// ----------------------------------------------------------------------
@@ -444,7 +424,7 @@
topology::Mesh mesh;
Neumann bc;
- _preinitialize(&mesh, &bc, true);
+ _preinitialize(&mesh, &bc);
spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDatabases");
spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
@@ -454,7 +434,7 @@
bc.dbInitial(&dbInitial);
- const PylithScalar timeScale = _TestNeumann::timeScale;
+ const PylithScalar timeScale = _data->timeScale;
bc._queryDatabases();
bc._calculateValue(_TestNeumann::tValue/timeScale);
@@ -464,10 +444,8 @@
CPPUNIT_ASSERT(0 != bc._parameters);
// Check values.
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann::_checkValues(_TestNeumann::initial,
- numQuadPts*spaceDim, value);
+ const topology::Field<topology::SubMesh>& value = bc._parameters->get("value");
+ _TestNeumann::_checkValues(_TestNeumann::initial, numQuadPts*spaceDim, value);
} // test_calculateValueInitial
// ----------------------------------------------------------------------
@@ -482,7 +460,7 @@
topology::Mesh mesh;
Neumann bc;
- _preinitialize(&mesh, &bc, true);
+ _preinitialize(&mesh, &bc);
spatialdata::spatialdb::SimpleDB dbRate("_TestNeumann _queryDatabases");
spatialdata::spatialdb::SimpleIOAscii dbRateIO;
@@ -492,7 +470,7 @@
bc.dbRate(&dbRate);
- const PylithScalar timeScale = _TestNeumann::timeScale;
+ const PylithScalar timeScale = _data->timeScale;
bc._queryDatabases();
bc._calculateValue(_TestNeumann::tValue/timeScale);
@@ -502,10 +480,8 @@
CPPUNIT_ASSERT(0 != bc._parameters);
// Check values.
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann::_checkValues(_TestNeumann::valuesRate,
- numQuadPts*spaceDim, value);
+ const topology::Field<topology::SubMesh>& value = bc._parameters->get("value");
+ _TestNeumann::_checkValues(_TestNeumann::valuesRate, numQuadPts*spaceDim, value);
} // test_calculateValueRate
// ----------------------------------------------------------------------
@@ -520,7 +496,7 @@
topology::Mesh mesh;
Neumann bc;
- _preinitialize(&mesh, &bc, true);
+ _preinitialize(&mesh, &bc);
spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann _queryDatabases");
spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
@@ -530,7 +506,7 @@
bc.dbChange(&dbChange);
- const PylithScalar timeScale = _TestNeumann::timeScale;
+ const PylithScalar timeScale = _data->timeScale;
bc._queryDatabases();
bc._calculateValue(_TestNeumann::tValue/timeScale);
@@ -540,10 +516,8 @@
CPPUNIT_ASSERT(0 != bc._parameters);
// Check values.
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann::_checkValues(_TestNeumann::valuesChange,
- numQuadPts*spaceDim, value);
+ const topology::Field<topology::SubMesh>& value = bc._parameters->get("value");
+ _TestNeumann::_checkValues(_TestNeumann::valuesChange, numQuadPts*spaceDim, value);
} // test_calculateValueChange
// ----------------------------------------------------------------------
@@ -558,9 +532,8 @@
topology::Mesh mesh;
Neumann bc;
- _preinitialize(&mesh, &bc, true);
+ _preinitialize(&mesh, &bc);
-
spatialdata::spatialdb::SimpleDB dbChange("_TestNeumann _queryDatabases");
spatialdata::spatialdb::SimpleIOAscii dbChangeIO;
dbChangeIO.filename("data/quad4_traction_change.spatialdb");
@@ -573,7 +546,7 @@
bc.dbChange(&dbChange);
bc.dbTimeHistory(&th);
- const PylithScalar timeScale = _TestNeumann::timeScale;
+ const PylithScalar timeScale = _data->timeScale;
bc._queryDatabases();
bc._calculateValue(_TestNeumann::tValue/timeScale);
@@ -583,10 +556,8 @@
CPPUNIT_ASSERT(0 != bc._parameters);
// Check values.
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
- _TestNeumann::_checkValues(_TestNeumann::valuesChangeTH,
- numQuadPts*spaceDim, value);
+ const topology::Field<topology::SubMesh>& value = bc._parameters->get("value");
+ _TestNeumann::_checkValues(_TestNeumann::valuesChangeTH, numQuadPts*spaceDim, value);
} // test_calculateValueChangeTH
// ----------------------------------------------------------------------
@@ -601,9 +572,8 @@
topology::Mesh mesh;
Neumann bc;
- _preinitialize(&mesh, &bc, true);
+ _preinitialize(&mesh, &bc);
-
spatialdata::spatialdb::SimpleDB dbInitial("_TestNeumann _queryDatabases");
spatialdata::spatialdb::SimpleIOAscii dbInitialIO;
dbInitialIO.filename("data/quad4_traction_initial.spatialdb");
@@ -630,7 +600,7 @@
bc.dbChange(&dbChange);
bc.dbTimeHistory(&th);
- const PylithScalar timeScale = _TestNeumann::timeScale;
+ const PylithScalar timeScale = _data->timeScale;
bc._queryDatabases();
bc._calculateValue(_TestNeumann::tValue/timeScale);
@@ -648,21 +618,19 @@
_TestNeumann::valuesRate[i] +
_TestNeumann::valuesChangeTH[i];
- const topology::Field<topology::SubMesh>& value =
- bc._parameters->get("value");
+ const topology::Field<topology::SubMesh>& value = bc._parameters->get("value");
_TestNeumann::_checkValues(&valuesE[0], numQuadPts*spaceDim, value);
} // test_calculateValueAll
// ----------------------------------------------------------------------
void
pylith::bc::TestNeumann::_preinitialize(topology::Mesh* mesh,
- Neumann* const bc,
- const bool useScales) const
+ Neumann* const bc) const
{ // _initialize
- CPPUNIT_ASSERT(0 != _data);
- CPPUNIT_ASSERT(0 != _quadrature);
- CPPUNIT_ASSERT(0 != mesh);
- CPPUNIT_ASSERT(0 != bc);
+ CPPUNIT_ASSERT(_data);
+ CPPUNIT_ASSERT(_quadrature);
+ CPPUNIT_ASSERT(mesh);
+ CPPUNIT_ASSERT(bc);
try {
// Set up mesh
@@ -674,15 +642,13 @@
spatialdata::geocoords::CSCart cs;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
+ mesh->coordsys(&cs);
spatialdata::units::Nondimensional normalizer;
- if (useScales) {
- normalizer.lengthScale(_TestNeumann::lengthScale);
- normalizer.pressureScale(_TestNeumann::pressureScale);
- normalizer.timeScale(_TestNeumann::timeScale);
- } // if
-
- mesh->coordsys(&cs);
+ normalizer.lengthScale(_data->lengthScale);
+ normalizer.pressureScale(_data->pressureScale);
+ normalizer.densityScale(_data->densityScale);
+ normalizer.timeScale(_data->timeScale);
mesh->nondimensionalize(normalizer);
// Set up quadrature
@@ -708,11 +674,11 @@
Neumann* const bc,
topology::SolutionFields* fields) const
{ // _initialize
- CPPUNIT_ASSERT(0 != _data);
- CPPUNIT_ASSERT(0 != mesh);
- CPPUNIT_ASSERT(0 != bc);
- CPPUNIT_ASSERT(0 != fields);
- CPPUNIT_ASSERT(0 != _quadrature);
+ CPPUNIT_ASSERT(_data);
+ CPPUNIT_ASSERT(mesh);
+ CPPUNIT_ASSERT(bc);
+ CPPUNIT_ASSERT(fields);
+ CPPUNIT_ASSERT(_quadrature);
try {
_preinitialize(mesh, bc);
@@ -738,6 +704,7 @@
topology::Field<topology::Mesh>& residual = fields->get("residual");
residual.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
residual.allocate();
+ residual.scale(_data->lengthScale);
residual.zero();
fields->copyLayout("residual");
@@ -753,21 +720,18 @@
const int fiberDimE,
const topology::Field<topology::SubMesh>& field)
{ // _checkValues
- assert(0 != valuesE);
+ CPPUNIT_ASSERT(valuesE);
const topology::SubMesh& boundaryMesh = field.mesh();
- DM subMesh = boundaryMesh.dmMesh();
+ PetscDM subMesh = boundaryMesh.dmMesh();CPPUNIT_ASSERT(subMesh);
PetscInt cStart, cEnd;
- PetscErrorCode err;
-
- CPPUNIT_ASSERT(subMesh);
+ PetscErrorCode err = 0;
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- const PylithScalar scale = field.scale();
- PetscSection section = field.petscSection();
- Vec vec = field.localVector();
- PetscScalar *array;
- CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+ const PylithScalar scale = field.scale();
+ PetscSection fieldSection = field.petscSection();CPPUNIT_ASSERT(fieldSection);
+ PetscVec fieldVec = field.localVector();assert(fieldVec);
+ PetscScalar *fieldArray;
const PetscInt ncells = _TestNeumann::ncells;
CPPUNIT_ASSERT_EQUAL(ncells, cEnd-cStart);
@@ -775,17 +739,17 @@
// Check values associated with BC.
int icell = 0;
const PylithScalar tolerance = 1.0e-06;
- err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
- for(PetscInt c = cStart; c < cEnd; ++c, ++icell) {
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ for (PetscInt c = cStart; c < cEnd; ++c, ++icell) {
PetscInt dof, off;
- err = PetscSectionGetDof(section, c, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(section, c, &off);CHECK_PETSC_ERROR(err);
+ 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, array[off+iDim], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[icell*fiberDimE+iDim]/scale, fieldArray[off+iDim], tolerance);
} // for
- err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
} // _checkValues
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.hh 2013-03-07 23:38:14 UTC (rev 21467)
@@ -119,11 +119,9 @@
*
* @param mesh Finite-element mesh to initialize
* @param bc Neumann boundary condition to initialize.
- * @param useScales Use scales provided by local constants.
*/
void _preinitialize(topology::Mesh* mesh,
- Neumann* const bc,
- const bool useScales =false) const;
+ Neumann* const bc) const;
/** Initialize Neumann boundary condition.
*
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -39,11 +39,6 @@
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::bc::TestPointForce );
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RealUniformSection RealUniformSection;
-
-// ----------------------------------------------------------------------
// Setup testing data.
void
pylith::bc::TestPointForce::setUp(void)
@@ -90,14 +85,13 @@
topology::Mesh mesh;
PointForce bc;
_initialize(&mesh, &bc);
- CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(_data);
- DM dmMesh = mesh.dmMesh();
- PetscInt cStart, cEnd;
+ PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd;
PetscErrorCode err;
+ err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- CPPUNIT_ASSERT(dmMesh);
- err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
const int numCells = cEnd-cStart;
const int numForceDOF = _data->numForceDOF;
const size_t numPoints = _data->numForcePts;
@@ -106,18 +100,21 @@
const int offset = numCells;
if (numForceDOF > 0) {
CPPUNIT_ASSERT_EQUAL(numPoints, bc._points.size());
- for (int i=0; i < numPoints; ++i)
+ for (int i=0; i < numPoints; ++i) {
CPPUNIT_ASSERT_EQUAL(_data->forcePoints[i]+offset, bc._points[i]);
+ } // for
} // if
- CPPUNIT_ASSERT(0 != bc._parameters);
- PetscSection initialSection = bc._parameters->get("initial").petscSection();
- Vec initialVec = bc._parameters->get("initial").localVector();
- PetscScalar *initialArray;
- CPPUNIT_ASSERT(initialSection);CPPUNIT_ASSERT(initialVec);
+ CPPUNIT_ASSERT(bc._parameters);
+ PetscSection initialSection = bc._parameters->get("initial").petscSection();CPPUNIT_ASSERT(initialSection);
+ PetscVec initialVec = bc._parameters->get("initial").localVector();CPPUNIT_ASSERT(initialVec);
+ PetscScalar *initialArray = NULL;
// Check values
const PylithScalar tolerance = 1.0e-06;
+ const PylithScalar forceScale = _data->pressureScale*pow(_data->lengthScale, 2);
+ const PylithScalar timeScale = _data->timeScale;
+
err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
for (int i=0; i < numPoints; ++i) {
const int p_force = _data->forcePoints[i]+offset;
@@ -126,18 +123,18 @@
err = PetscSectionGetDof(initialSection, p_force, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(initialSection, p_force, &off);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(numForceDOF, dof);
- for (int iDOF=0; iDOF < numForceDOF; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->forceInitial[i*numForceDOF+iDOF], initialArray[off+iDOF], tolerance);
+ for (int iDOF=0; iDOF < numForceDOF; ++iDOF) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->forceInitial[i*numForceDOF+iDOF], initialArray[off+iDOF]*forceScale, tolerance);
+ } // for
} // for
err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
// Check rate of change
- PetscSection rateSection = bc._parameters->get("rate").petscSection();
- Vec rateVec = bc._parameters->get("rate").localVector();
- PetscScalar *rateArray;
- CPPUNIT_ASSERT(rateSection);CPPUNIT_ASSERT(rateVec);
+ PetscSection rateSection = bc._parameters->get("rate").petscSection();CPPUNIT_ASSERT(rateSection);
+ PetscVec rateVec = bc._parameters->get("rate").localVector();CPPUNIT_ASSERT(rateVec);
+ PetscScalar *rateArray = NULL;
+ err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
- err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
for (int i=0; i < numPoints; ++i) {
const int p_force = _data->forcePoints[i]+offset;
PetscInt dof, off;
@@ -145,8 +142,9 @@
err = PetscSectionGetDof(rateSection, p_force, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(rateSection, p_force, &off);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(numForceDOF, dof);
- for (int iDOF=0; iDOF < numForceDOF; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->forceRate, rateArray[off+iDOF], tolerance);
+ for (int iDOF=0; iDOF < numForceDOF; ++iDOF) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->forceRate, rateArray[off+iDOF]*forceScale/timeScale, tolerance);
+ } // for
} // for
err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
} // testInitialize
@@ -162,7 +160,7 @@
topology::Field<topology::Mesh> residual(mesh);
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- CPPUNIT_ASSERT(0 != cs);
+ CPPUNIT_ASSERT(cs);
const int spaceDim = cs->spaceDim();
residual.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
residual.allocate();
@@ -170,38 +168,38 @@
topology::SolutionFields fields(mesh);
- const PylithScalar t = _data->tResidual;
+ const PylithScalar t = _data->tResidual/_data->timeScale;
bc.integrateResidual(residual, t, &fields);
- DM dmMesh = mesh.dmMesh();
- PetscInt vStart, vEnd;
+ PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
PetscErrorCode err;
-
- CPPUNIT_ASSERT(dmMesh);
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const PylithScalar* valsE = _data->residual;
const int totalNumVertices = vEnd-vStart;
const int sizeE = spaceDim * totalNumVertices;
- PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
- PetscScalar *vals;
- PetscInt size;
-
- CPPUNIT_ASSERT(residualSection);
+ PetscSection residualSection = residual.petscSection();CPPUNIT_ASSERT(residualSection);
+ PetscVec residualVec = residual.localVector();CPPUNIT_ASSERT(residualVec);
+ PetscScalar *vals = NULL;
+ PetscInt size;
err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
- //residual.view("RESIDUAL");
+ // residual.view("RESIDUAL"); // DEBUGGING
const PylithScalar tolerance = 1.0e-06;
+ const PylithScalar forceScale = _data->pressureScale*pow(_data->lengthScale, 2);
+ const PylithScalar residualScale = forceScale;
+
err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
for (int i=0; i < size; ++i)
- if (fabs(valsE[i]) > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
- else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ if (fabs(valsE[i]) > 1.0) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i]*residualScale, tolerance);
+ } else {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i]*residualScale, tolerance);
+ } // if/else
err = VecRestoreArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateResidual
@@ -222,18 +220,24 @@
pylith::bc::TestPointForce::_initialize(topology::Mesh* mesh,
PointForce* const bc) const
{ // _initialize
- CPPUNIT_ASSERT(0 != _data);
- CPPUNIT_ASSERT(0 != bc);
+ CPPUNIT_ASSERT(_data);
+ CPPUNIT_ASSERT(bc);
meshio::MeshIOAscii iohandler;
iohandler.filename(_data->meshFilename);
iohandler.read(mesh);
+ // Set coordinate system
spatialdata::geocoords::CSCart cs;
- spatialdata::units::Nondimensional normalizer;
cs.setSpaceDim(mesh->dimension());
cs.initialize();
mesh->coordsys(&cs);
+
+ spatialdata::units::Nondimensional normalizer;
+ normalizer.lengthScale(_data->lengthScale);
+ normalizer.pressureScale(_data->pressureScale);
+ normalizer.densityScale(_data->densityScale);
+ normalizer.timeScale(_data->timeScale);
mesh->nondimensionalize(normalizer);
spatialdata::spatialdb::SimpleDB dbInitial("TestPointForce initial");
@@ -270,6 +274,7 @@
bc->dbInitial(&dbInitial);
bc->dbRate(&dbRate);
bc->bcDOF(_data->forceDOF, _data->numForceDOF);
+ bc->normalizer(normalizer);
bc->initialize(*mesh, upDir);
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletDataMulti.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletDataMulti.cc 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletDataMulti.cc 2013-03-07 23:38:14 UTC (rev 21467)
@@ -53,8 +53,14 @@
fieldIncr(0), // General
constraintSizes(0),
constrainedDOF(0),
- meshFilename(0)
+ meshFilename(0),
+ lengthScale(1.0e+3),
+ pressureScale(2.25e+10),
+ densityScale(1.0),
+ timeScale(2.0)
{ // constructor
+ const PylithScalar velScale = lengthScale / timeScale;
+ densityScale = pressureScale / (velScale*velScale);
} // constructor
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletDataMulti.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletDataMulti.hh 2013-03-07 22:56:07 UTC (rev 21466)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletDataMulti.hh 2013-03-07 23:38:14 UTC (rev 21467)
@@ -86,6 +86,14 @@
int* constrainedDOF; ///< Indices of constrained DOF at each constrained vertex
char* meshFilename; ///< Filename for input mesh.
+
+ /// @name Scales information for nondimensionalization.
+ //@{
+ PylithScalar lengthScale; ///< Length scale.
+ PylithScalar pressureScale; ///< Pressure scale.
+ PylithScalar timeScale; ///< Time scale.
+ PylithScalar densityScale; ///< Density scale.
+ //@}
};
#endif // pylith_bc_dirichletdatamulti_hh
More information about the CIG-COMMITS
mailing list