[cig-commits] r21456 - in short/3D/PyLith/trunk: libsrc/pylith/bc unittests/libtests/bc unittests/libtests/bc/data
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 6 16:43:39 PST 2013
Author: brad
Date: 2013-03-06 16:43:38 -0800 (Wed, 06 Mar 2013)
New Revision: 21456
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
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/data/AbsorbingDampersData.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.hh
short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletData.cc
short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletData.hh
Log:
Added scales to unit tests for DirichletDB. Code cleanup.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-03-07 00:43:38 UTC (rev 21456)
@@ -36,14 +36,9 @@
#include <stdexcept> // USES std::runtime_error
#include <sstream> // USES std::ostringstream
-//#define PRECOMPUTE_GEOMETRY
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
-typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::AbsorbingDampers::AbsorbingDampers(void) :
_db(0)
@@ -73,20 +68,21 @@
pylith::bc::AbsorbingDampers::initialize(const topology::Mesh& mesh,
const PylithScalar upDir[3])
{ // initialize
- assert(0 != _boundaryMesh);
- assert(0 != _quadrature);
- assert(0 != _db);
+ assert(_boundaryMesh);
+ assert(_quadrature);
+ assert(_db);
_initializeLogger();
scalar_array up(3);
- for (int i=0; i < 3; ++i)
+ for (int i=0; i < 3; ++i) {
up[i] = upDir[i];
+ } // for
const int numCorners = _quadrature->numBasis();
// Get 'surface' cells (1 dimension lower than top-level cells)
- DM subMesh = _boundaryMesh->dmMesh();
+ PetscDM subMesh = _boundaryMesh->dmMesh();
PetscInt cStart, cEnd;
PetscErrorCode err;
@@ -106,7 +102,7 @@
delete _parameters;
_parameters = new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
- assert(0 != _parameters);
+ assert(_parameters);
_parameters->add("damping constants", "damping_constants", topology::FieldBase::FACES_FIELD, fiberDim);
_parameters->get("damping constants").allocate();
@@ -142,21 +138,20 @@
scalar_array coordinatesCell(numBasis*spaceDim);
PetscSection coordSection;
- Vec coordVec;
+ PetscVec coordVec;
err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
assert(coordSection);assert(coordVec);
- assert(0 != _normalizer);
+ assert(_normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
const PylithScalar densityScale = _normalizer->densityScale();
assert(_normalizer->timeScale() > 0);
- const PylithScalar velocityScale =
- _normalizer->lengthScale() / _normalizer->timeScale();
+ const PylithScalar velocityScale = _normalizer->lengthScale() / _normalizer->timeScale();
PetscSection valueSection = _parameters->get("damping constants").petscSection();
- Vec valueVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingArray;
+ PetscVec valueVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingConstsArray;
assert(valueSection);assert(valueVec);
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
@@ -167,16 +162,18 @@
_quadrature->computeGeometry(*_boundaryMesh, cells);
#endif
- err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ 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;
+ 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
@@ -219,25 +216,21 @@
// Compute normal/tangential orientation
memcpy(&quadPtRef[0], &quadPtsRef[iQuad*cellDim], cellDim*sizeof(PylithScalar));
-#if defined(PRECOMPUTE_GEOMETRY)
- coordsVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
-#endif
cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell, quadPtRef);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
assert(jacobianDet > 0.0);
orientation /= jacobianDet;
for (int iDim=0; iDim < spaceDim; ++iDim) {
- dampingArray[doff+iQuad*spaceDim+iDim] = 0.0;
+ dampingConstsArray[doff+iQuad*spaceDim+iDim] = 0.0;
for (int jDim=0; jDim < spaceDim; ++jDim)
- dampingArray[doff+iQuad*spaceDim+iDim] += dampingConstsLocal[jDim]*orientation[jDim*spaceDim+iDim];
+ dampingConstsArray[doff+iQuad*spaceDim+iDim] += dampingConstsLocal[jDim]*orientation[jDim*spaceDim+iDim];
// Ensure damping constants are positive
- dampingArray[doff+iQuad*spaceDim+iDim] = fabs(dampingArray[doff+iQuad*spaceDim+iDim]);
+ dampingConstsArray[doff+iQuad*spaceDim+iDim] = fabs(dampingConstsArray[doff+iQuad*spaceDim+iDim]);
} // for
} // for
} // for
- err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(valueVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
_db->close();
} // initialize
@@ -250,11 +243,11 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateResidual
- assert(0 != _quadrature);
- assert(0 != _boundaryMesh);
- assert(0 != _parameters);
- assert(0 != fields);
- assert(0 != _logger);
+ assert(_quadrature);
+ assert(_boundaryMesh);
+ assert(_parameters);
+ assert(fields);
+ assert(_logger);
const int setupEvent = _logger->eventId("AdIR setup");
const int geometryEvent = _logger->eventId("AdIR geometry");
@@ -275,8 +268,8 @@
_initCellVector();
// Get cell information
- DM subMesh = _boundaryMesh->dmMesh();
- IS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();
+ PetscIS subpointIS;
PetscInt cStart, cEnd;
PetscErrorCode err;
@@ -285,19 +278,19 @@
err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
// Get sections
- PetscSection valueSection = _parameters->get("damping constants").petscSection();
- Vec valueVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingArray;
- assert(valueSection);assert(valueVec);
+ PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();
+ PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingConstsArray;
+ assert(dampingConstsSection);assert(dampingConstsVec);
// Use _cellVector for cell residual.
PetscSection residualSection = residual.petscSection(), residualSubsection;
- Vec residualVec = residual.localVector();
+ PetscVec residualVec = residual.localVector();
assert(residualSection);assert(residualVec);
err = PetscSectionCreateSubmeshSection(residualSection, subpointIS, &residualSubsection);
PetscSection velSection = fields->get("velocity(t)").petscSection(), velSubsection;
- Vec velVec = fields->get("velocity(t)").localVector();
+ PetscVec velVec = fields->get("velocity(t)").localVector();
assert(velSection);assert(velVec);
err = PetscSectionCreateSubmeshSection(velSection, subpointIS, &velSubsection);
err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
@@ -305,7 +298,7 @@
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
PetscSection coordSection;
- Vec coordVec;
+ PetscVec coordVec;
err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
assert(coordSection);assert(coordVec);
@@ -316,7 +309,7 @@
_logger->eventBegin(computeEvent);
#endif
- err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
for(PetscInt c = cStart; c < cEnd; ++c) {
// Get geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
@@ -325,12 +318,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- const PetscScalar *coords;
+ const PetscScalar *coordsArray;
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];}
+ err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsArray[i];}
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -342,12 +335,12 @@
_resetCellVector();
// Restrict input fields to cell
- const PetscScalar *velArray = PETSC_NULL;
+ const PetscScalar *velArray = NULL;
PetscInt velSize;
PetscInt ddof, doff;
- err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(dampingConstsSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dampingConstsSection, c, &doff);CHECK_PETSC_ERROR(err);
assert(ddof == numQuadPts*spaceDim);
err = DMPlexVecGetClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
assert(velSize == numBasis*spaceDim);
@@ -371,7 +364,7 @@
const PylithScalar valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -=
- dampingArray[doff+iQuad*spaceDim+iDim] *
+ dampingConstsArray[doff+iQuad*spaceDim+iDim] *
valIJ * velArray[jBasis*spaceDim+iDim];
} // for
} // for
@@ -389,7 +382,7 @@
_logger->eventEnd(updateEvent);
#endif
} // for
- err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
err = PetscSectionDestroy(&residualSubsection);CHECK_PETSC_ERROR(err);
err = PetscSectionDestroy(&velSubsection);CHECK_PETSC_ERROR(err);
@@ -407,11 +400,11 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateResidualLumped
- assert(0 != _quadrature);
- assert(0 != _boundaryMesh);
- assert(0 != _parameters);
- assert(0 != fields);
- assert(0 != _logger);
+ assert(_quadrature);
+ assert(_boundaryMesh);
+ assert(_parameters);
+ assert(fields);
+ assert(_logger);
const int setupEvent = _logger->eventId("AdIR setup");
const int geometryEvent = _logger->eventId("AdIR geometry");
@@ -432,8 +425,8 @@
_initCellVector();
// Get cell information
- DM subMesh = _boundaryMesh->dmMesh();
- IS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();
+ PetscIS subpointIS;
PetscInt cStart, cEnd;
PetscErrorCode err;
@@ -442,19 +435,19 @@
err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
// Get sections
- PetscSection valueSection = _parameters->get("damping constants").petscSection();
- Vec valueVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingArray;
- assert(valueSection);assert(valueVec);
+ PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();
+ PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingConstsArray;
+ assert(dampingConstsSection);assert(dampingConstsVec);
// Use _cellVector for cell values.
PetscSection residualSection = residual.petscSection(), residualSubsection;
- Vec residualVec = residual.localVector();
+ PetscVec residualVec = residual.localVector();
assert(residualSection);assert(residualVec);
err = PetscSectionCreateSubmeshSection(residualSection, subpointIS, &residualSubsection);
PetscSection velSection = fields->get("velocity(t)").petscSection(), velSubsection;
- Vec velVec = fields->get("velocity(t)").localVector();
+ PetscVec velVec = fields->get("velocity(t)").localVector();
assert(velSection);assert(velVec);
err = PetscSectionCreateSubmeshSection(velSection, subpointIS, &velSubsection);
err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
@@ -462,7 +455,7 @@
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
PetscSection coordSection;
- Vec coordVec;
+ PetscVec coordVec;
err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
assert(coordSection);assert(coordVec);
@@ -473,7 +466,7 @@
_logger->eventBegin(computeEvent);
#endif
- err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
for(PetscInt c = cStart; c < cEnd; ++c) {
// Get geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
@@ -482,12 +475,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- const PetscScalar *coords;
+ const PetscScalar *coordsArray;
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];}
+ err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsArray[i];}
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -499,12 +492,12 @@
_resetCellVector();
// Restrict input fields to cell
- const PetscScalar *velArray = PETSC_NULL;
+ const PetscScalar *velArray = NULL;
PetscInt velSize;
PetscInt ddof, doff;
- err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(dampingConstsSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dampingConstsSection, c, &doff);CHECK_PETSC_ERROR(err);
assert(ddof == numQuadPts*spaceDim);
err = DMPlexVecGetClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
assert(velSize == numBasis*spaceDim);
@@ -530,7 +523,7 @@
const PylithScalar valIJ = basis[iQ + iBasis] * valJ;
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -=
- dampingArray[doff+iQuad*spaceDim+iDim] *
+ dampingConstsArray[doff+iQuad*spaceDim+iDim] *
valIJ * velArray[iBasis*spaceDim+iDim];
} // for
} // for
@@ -546,7 +539,7 @@
_logger->eventEnd(updateEvent);
#endif
} // for
- err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
err = PetscSectionDestroy(&residualSubsection);CHECK_PETSC_ERROR(err);
err = PetscSectionDestroy(&velSubsection);CHECK_PETSC_ERROR(err);
@@ -564,11 +557,11 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateJacobian
- assert(0 != _quadrature);
- assert(0 != _boundaryMesh);
- assert(0 != _logger);
- assert(0 != jacobian);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_boundaryMesh);
+ assert(_logger);
+ assert(jacobian);
+ assert(fields);
const int setupEvent = _logger->eventId("AdIJ setup");
const int geometryEvent = _logger->eventId("AdIJ geometry");
@@ -586,24 +579,24 @@
const int spaceDim = _quadrature->spaceDim();
// Get cell information
- DM subMesh = _boundaryMesh->dmMesh();
- IS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();
+ PetscIS subpointIS;
PetscInt cStart, cEnd;
- PetscErrorCode err;
+ PetscErrorCode err = 0;
assert(subMesh);
err = DMPlexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
err = DMPlexCreateSubpointIS(subMesh, &subpointIS);CHECK_PETSC_ERROR(err);
// Get sections
- PetscSection valueSection = _parameters->get("damping constants").petscSection();
- Vec valueVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingArray;
- assert(valueSection);assert(valueVec);
+ PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();
+ PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingConstsArray;
+ assert(dampingConstsSection);assert(dampingConstsVec);
const topology::Field<topology::Mesh>& solution = fields->solution();
PetscSection solutionSection = solution.petscSection(), solutionGlobalSection, solutionSubsection, solutionGlobalSubsection;
- Vec solutionVec = solution.localVector();
+ PetscVec solutionVec = solution.localVector();
PetscSF sf;
assert(solutionSection);assert(solutionVec);
err = PetscSectionCreateSubmeshSection(solutionSection, subpointIS, &solutionSubsection);CHECK_PETSC_ERROR(err);
@@ -614,7 +607,7 @@
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
- assert(0 != jacobianMat);
+ assert(jacobianMat);
// Get parameters used in integration.
const PylithScalar dt = _dt;
@@ -626,7 +619,7 @@
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
PetscSection coordSection;
- Vec coordVec;
+ PetscVec coordVec;
err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
assert(coordSection);assert(coordVec);
@@ -637,7 +630,7 @@
_logger->eventBegin(computeEvent);
#endif
- err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
@@ -646,12 +639,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- const PetscScalar *coords;
+ const PetscScalar *coordsArray;
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];}
+ err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsArray[i];}
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -661,8 +654,8 @@
// Get damping constants
PetscInt ddof, doff;
- err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(dampingConstsSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dampingConstsSection, c, &doff);CHECK_PETSC_ERROR(err);
assert(ddof == numQuadPts*spaceDim);
#if defined(DETAILED_EVENT_LOGGING)
@@ -688,7 +681,7 @@
for (int iDim=0; iDim < spaceDim; ++iDim) {
const int iBlock = (iBasis*spaceDim + iDim) * (numBasis*spaceDim);
const int jBlock = (jBasis*spaceDim + iDim);
- _cellMatrix[iBlock+jBlock] += valIJ * dampingArray[doff+iQuad*spaceDim+iDim];
+ _cellMatrix[iBlock+jBlock] += valIJ * dampingConstsArray[doff+iQuad*spaceDim+iDim];
} // for
} // for
} // for
@@ -707,7 +700,7 @@
_logger->eventEnd(updateEvent);
#endif
} // for
- err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
err = PetscSectionDestroy(&solutionSubsection);CHECK_PETSC_ERROR(err);
err = PetscSectionDestroy(&solutionGlobalSection);CHECK_PETSC_ERROR(err);
err = PetscSectionDestroy(&solutionGlobalSubsection);CHECK_PETSC_ERROR(err);
@@ -728,11 +721,11 @@
const PylithScalar t,
topology::SolutionFields* const fields)
{ // integrateJacobian
- assert(0 != _quadrature);
- assert(0 != _boundaryMesh);
- assert(0 != _logger);
- assert(0 != jacobian);
- assert(0 != fields);
+ assert(_quadrature);
+ assert(_boundaryMesh);
+ assert(_logger);
+ assert(jacobian);
+ assert(fields);
const int setupEvent = _logger->eventId("AdIJ setup");
const int geometryEvent = _logger->eventId("AdIJ geometry");
@@ -750,8 +743,8 @@
const int spaceDim = _quadrature->spaceDim();
// Get cell information
- DM subMesh = _boundaryMesh->dmMesh();
- IS subpointIS;
+ PetscDM subMesh = _boundaryMesh->dmMesh();
+ PetscIS subpointIS;
PetscInt cStart, cEnd;
PetscErrorCode err;
@@ -768,13 +761,13 @@
_initCellVector();
// Get sections
- PetscSection valueSection = _parameters->get("damping constants").petscSection();
- Vec valueVec = _parameters->get("damping constants").localVector();
- PetscScalar *dampingArray;
- assert(valueSection);assert(valueVec);
+ PetscSection dampingConstsSection = _parameters->get("damping constants").petscSection();
+ PetscVec dampingConstsVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingConstsArray;
+ assert(dampingConstsSection);assert(dampingConstsVec);
PetscSection jacobianSection = jacobian->petscSection(), jacobianSubsection;
- Vec jacobianVec = jacobian->localVector();
+ PetscVec jacobianVec = jacobian->localVector();
assert(jacobianSection);assert(jacobianVec);
err = PetscSectionCreateSubmeshSection(jacobianSection, subpointIS, &jacobianSubsection);
err = ISDestroy(&subpointIS);CHECK_PETSC_ERROR(err);
@@ -782,7 +775,7 @@
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
PetscSection coordSection;
- Vec coordVec;
+ PetscVec coordVec;
err = DMPlexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
assert(coordSection);assert(coordVec);
@@ -793,7 +786,7 @@
_logger->eventBegin(computeEvent);
#endif
- err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
@@ -802,12 +795,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- const PetscScalar *coords;
+ const PetscScalar *coordsArray;
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];}
+ err = DMPlexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coordsArray[i];}
_quadrature->computeGeometry(coordinatesCell, c);
- err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ err = DMPlexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coordsArray);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -817,8 +810,8 @@
// Get damping constants
PetscInt ddof, doff;
- err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(dampingConstsSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dampingConstsSection, c, &doff);CHECK_PETSC_ERROR(err);
assert(ddof == numQuadPts*spaceDim);
#if defined(DETAILED_EVENT_LOGGING)
@@ -845,7 +838,7 @@
const PylithScalar valIJ = basis[iQ + iBasis] * valJ;
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis * spaceDim + iDim] += valIJ
- * dampingArray[doff+iQuad * spaceDim + iDim];
+ * dampingConstsArray[doff+iQuad * spaceDim + iDim];
} // for
} // for
err = DMPlexVecSetClosure(subMesh, jacobianSubsection, jacobianVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
@@ -858,7 +851,7 @@
_logger->eventEnd(updateEvent);
#endif
} // for
- err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
err = PetscSectionDestroy(&jacobianSubsection);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -883,7 +876,7 @@
pylith::bc::AbsorbingDampers::_initializeLogger(void)
{ // initializeLogger
delete _logger; _logger = new utils::EventLogger;
- assert(0 != _logger);
+ assert(_logger);
_logger->className("AbsorbingDampers");
_logger->initialize();
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2013-03-07 00:43:38 UTC (rev 21456)
@@ -83,8 +83,8 @@
return;
PetscSection sectionP = field.petscSection();
- PetscInt numFields;
- PetscErrorCode err;
+ PetscInt numFields;
+ PetscErrorCode err = 0;
assert(sectionP);
// Set constraints in field
@@ -93,7 +93,7 @@
err = PetscSectionGetNumFields(sectionP, &numFields);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const PetscInt point = _points[iPoint];
- PetscInt dof, cdof;
+ PetscInt dof, cdof;
err = PetscSectionGetDof(sectionP, point, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetConstraintDof(sectionP, point, &cdof);CHECK_PETSC_ERROR(err);
@@ -122,9 +122,9 @@
if (0 == numFixedDOF)
return;
- PetscSection sectionP = field.petscSection();
- PetscInt numFields;
- PetscErrorCode err;
+ PetscSection sectionP = field.petscSection();
+ PetscInt numFields;
+ PetscErrorCode err = 0;
assert(sectionP);
const int numPoints = _points.size();
@@ -194,14 +194,14 @@
assert(_parameters);
PetscSection valueSection = _parameters->get("value").petscSection();
- Vec valueVec = _parameters->get("value").localVector();
+ PetscVec valueVec = _parameters->get("value").localVector();
PetscScalar *valueArray;
assert(valueSection);assert(valueVec);
- PetscSection fieldSection = field.petscSection();
- Vec localVec = field.localVector();
- PetscScalar *array;
- PetscErrorCode err;
+ 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);
@@ -240,14 +240,14 @@
assert(_parameters);
PetscSection valueSection = _parameters->get("value").petscSection();
- Vec valueVec = _parameters->get("value").localVector();
+ PetscVec valueVec = _parameters->get("value").localVector();
PetscScalar *valueArray;
assert(valueSection);assert(valueVec);
- PetscSection fieldSection = field.petscSection();
- Vec localVec = field.localVector();
- PetscScalar *array;
- PetscErrorCode err;
+ 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);
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2013-03-07 00:43:38 UTC (rev 21456)
@@ -34,11 +34,6 @@
#include <sstream> // USES std::ostringstream
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealUniformSection RealUniformSection;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::TimeDependentPoints::TimeDependentPoints(void)
{ // constructor
@@ -83,8 +78,7 @@
{ // _queryDatabases
const PylithScalar timeScale = _getNormalizer().timeScale();
const PylithScalar rateScale = valueScale / timeScale;
- Vec v;
- PetscErrorCode err;
+ PetscErrorCode err = 0;
const int numPoints = _points.size();
const int numBCDOF = _bcDOF.size();
@@ -124,52 +118,58 @@
_parameters = new topology::Fields<topology::Field<topology::Mesh> >(mesh);
_parameters->add("value", "value", topology::FieldBase::VERTICES_FIELD, numBCDOF);
- _parameters->get("value").vectorFieldType(topology::FieldBase::OTHER);
- _parameters->get("value").scale(valueScale);
- _parameters->get("value").allocate();
- v = _parameters->get("value").localVector();assert(v);
- err = VecSet(v, 0.0);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::Mesh>& value = _parameters->get("value");
+ value.vectorFieldType(topology::FieldBase::OTHER);
+ value.scale(valueScale);
+ value.allocate();
+ PetscVec valueVec = value.localVector();assert(valueVec);
+ err = VecSet(valueVec, 0.0);CHECK_PETSC_ERROR(err);
if (_dbInitial) {
const std::string& fieldLabel = std::string("initial_") + std::string(fieldName);
_parameters->add("initial", fieldLabel.c_str(), topology::FieldBase::VERTICES_FIELD, numBCDOF);
- _parameters->get("initial").vectorFieldType(topology::FieldBase::OTHER);
- _parameters->get("initial").scale(valueScale);
- _parameters->get("initial").allocate();
- v = _parameters->get("initial").localVector();assert(v);
- err = VecSet(v, 0.0);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::Mesh>& initial = _parameters->get("initial");
+ initial.vectorFieldType(topology::FieldBase::OTHER);
+ initial.scale(valueScale);
+ initial.allocate();
+ PetscVec initialVec = initial.localVector();assert(initialVec);
+ err = VecSet(initialVec, 0.0);CHECK_PETSC_ERROR(err);
} // if
if (_dbRate) {
const std::string& fieldLabel = std::string("rate_") + std::string(fieldName);
_parameters->add("rate", fieldLabel.c_str(), topology::FieldBase::VERTICES_FIELD, numBCDOF);
- _parameters->get("rate").vectorFieldType(topology::FieldBase::OTHER);
- _parameters->get("rate").scale(rateScale);
- _parameters->get("rate").allocate();
- v = _parameters->get("rate").localVector();assert(v);
- err = VecSet(v, 0.0);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::Mesh>& rate = _parameters->get("rate");
+ rate.vectorFieldType(topology::FieldBase::OTHER);
+ rate.scale(rateScale);
+ rate.allocate();
+ PetscVec rateVec = rate.localVector();assert(rateVec);
+ err = VecSet(rateVec, 0.0);CHECK_PETSC_ERROR(err);
const std::string& timeLabel = std::string("rate_time_") + std::string(fieldName);
_parameters->add("rate time", timeLabel.c_str(), topology::FieldBase::VERTICES_FIELD, 1);
- _parameters->get("rate time").vectorFieldType(topology::FieldBase::SCALAR);
- _parameters->get("rate time").scale(timeScale);
- _parameters->get("rate time").allocate();
- v = _parameters->get("rate time").localVector();assert(v);
- err = VecSet(v, 0.0);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::Mesh>& rateTime = _parameters->get("rate time");
+ rateTime.vectorFieldType(topology::FieldBase::SCALAR);
+ rateTime.scale(timeScale);
+ rateTime.allocate();
+ PetscVec rateTimeVec = rateTime.localVector();assert(rateTimeVec);
+ err = VecSet(rateTimeVec, 0.0);CHECK_PETSC_ERROR(err);
} // if
if (_dbChange) {
const std::string& fieldLabel = std::string("change_") + std::string(fieldName);
_parameters->add("change", fieldLabel.c_str(), topology::FieldBase::VERTICES_FIELD, numBCDOF);
- _parameters->get("change").vectorFieldType(topology::FieldBase::OTHER);
- _parameters->get("change").scale(valueScale);
- _parameters->get("change").allocate();
- v = _parameters->get("change").localVector();assert(v);
- err = VecSet(v, 0.0);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::Mesh>& change = _parameters->get("change");
+ change.vectorFieldType(topology::FieldBase::OTHER);
+ change.scale(valueScale);
+ change.allocate();
+ PetscVec changeVec = change.localVector();assert(changeVec);
+ err = VecSet(changeVec, 0.0);CHECK_PETSC_ERROR(err);
const std::string& timeLabel = std::string("change_time_") + std::string(fieldName);
_parameters->add("change time", timeLabel.c_str(), topology::FieldBase::VERTICES_FIELD, 1);
- _parameters->get("change time").vectorFieldType(topology::FieldBase::SCALAR);
- _parameters->get("change time").scale(timeScale);
- _parameters->get("change time").allocate();
- v = _parameters->get("change time").localVector();assert(v);
- err = VecSet(v, 0.0);CHECK_PETSC_ERROR(err);
+ topology::Field<topology::Mesh>& changeTime = _parameters->get("change time");
+ changeTime.vectorFieldType(topology::FieldBase::SCALAR);
+ changeTime.scale(timeScale);
+ changeTime.allocate();
+ PetscVec changeTimeVec = _parameters->get("change time").localVector();assert(changeTimeVec);
+ err = VecSet(changeTimeVec, 0.0);CHECK_PETSC_ERROR(err);
} // if
if (_dbInitial) { // Setup initial values, if provided.
@@ -280,8 +280,7 @@
msg << ") using spatial database " << db->label() << ".";
throw std::runtime_error(msg.str());
} // if
- _getNormalizer().nondimensionalize(&valueVertex[0], valueVertex.size(),
- scale);
+ _getNormalizer().nondimensionalize(&valueVertex[0], valueVertex.size(), scale);
// Update section
PetscInt dof, off;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2013-03-07 00:43:38 UTC (rev 21456)
@@ -124,7 +124,7 @@
PetscInt dp = 0;
for(PetscInt c = cStart; c < cEnd; ++c) {
- PetscInt *closure = PETSC_NULL;
+ PetscInt *closure = NULL;
PetscInt closureSize, numCorners = 0;
err = DMPlexGetTransitiveClosure(subMesh, c, PETSC_TRUE, &closureSize, &closure);CHECK_PETSC_ERROR(err);
@@ -146,26 +146,28 @@
const int fiberDim = numQuadPts * spaceDim;
PetscInt index = 0;
CPPUNIT_ASSERT(0 != bc._parameters);
- PetscSection dampersSection = bc._parameters->get("damping constants").petscSection();
- Vec dampersVec = bc._parameters->get("damping constants").localVector();
- CPPUNIT_ASSERT(dampersSection);CPPUNIT_ASSERT(dampersVec);
+ PetscSection dampingConstsSection = bc._parameters->get("damping constants").petscSection();
+ PetscVec dampingConstsVec = bc._parameters->get("damping constants").localVector();
+ CPPUNIT_ASSERT(dampingConstsSection);CPPUNIT_ASSERT(dampingConstsVec);
const PylithScalar tolerance = 1.0e-06;
- PetscScalar *dampersValues;
- err = VecGetArray(dampersVec, &dampersValues);CHECK_PETSC_ERROR(err);
+ const PylithScalar* dampingConstsE = _data->dampingConsts;
+ const PylithScalar dampingConstsScale = _data->densityScale * _data->lengthScale / _data->timeScale;
+ PetscScalar* dampingConstsArray;
+ err = VecGetArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
for(PetscInt c = cStart; c < cEnd; ++c) {
PetscInt dof, off;
- err = PetscSectionGetDof(dampersSection, c, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(dampersSection, c, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(dampingConstsSection, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dampingConstsSection, c, &off);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
for(PetscInt iQuad=0; iQuad < numQuadPts; ++iQuad)
for(PetscInt iDim =0; iDim < spaceDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, dampersValues[off+iQuad*spaceDim+iDim]/_data->dampingConsts[index], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, dampingConstsArray[off+iQuad*spaceDim+iDim]/dampingConstsE[index]*dampingConstsScale, tolerance);
++index;
} // for
} // for
- err = VecRestoreArray(dampersVec, &dampersValues);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dampingConstsVec, &dampingConstsArray);CHECK_PETSC_ERROR(err);
} // testInitialize
// ----------------------------------------------------------------------
@@ -191,6 +193,9 @@
DM dmMesh = mesh.dmMesh();
PetscInt vStart, vEnd;
const PylithScalar* valsE = _data->valsResidual;
+ const PylithScalar dampingConstsScale = _data->densityScale * _data->lengthScale / _data->timeScale;
+ const PylithScalar velocityScale = 1.0; // Input velocity is nondimensional.
+ const PylithScalar residualScale = dampingConstsScale*velocityScale*pow(_data->lengthScale, _data->spaceDim-1);
CPPUNIT_ASSERT(dmMesh);
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
@@ -198,7 +203,7 @@
const int sizeE = _data->spaceDim * totalNumVertices;
PetscSection residualSection = residual.petscSection();
- Vec residualVec = residual.localVector();
+ PetscVec residualVec = residual.localVector();
PetscScalar *vals;
PetscInt size;
@@ -212,9 +217,9 @@
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);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i]*residualScale, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i]/residualScale, vals[i], tolerance);
err = VecRestoreArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateResidual
@@ -248,10 +253,12 @@
CPPUNIT_ASSERT(dmMesh);
err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const PylithScalar* valsE = _data->valsJacobian;
+ const PylithScalar dampingConstsScale = _data->densityScale * _data->lengthScale / _data->timeScale;
+ const PylithScalar jacobianScale = dampingConstsScale*pow(_data->lengthScale, _data->spaceDim-1);
+
const int totalNumVertices = vEnd - vStart;
const int nrowsE = totalNumVertices * _data->spaceDim;
const int ncolsE = totalNumVertices * _data->spaceDim;
-
const PetscMat jacobianMat = jacobian.matrix();
int nrows = 0;
int ncols = 0;
@@ -283,9 +290,9 @@
for (int iCol=0; iCol < ncols; ++iCol) {
const int index = ncols*iRow+iCol;
if (fabs(valsE[index]) > 1.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[index]/valsE[index], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[index]/valsE[index]*jacobianScale, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[index], vals[index], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[index]/jacobianScale, vals[index], tolerance);
} // for
err = MatDestroy(&jDense);CHECK_PETSC_ERROR(err);
} // testIntegrateJacobian
@@ -327,6 +334,9 @@
const int totalNumVertices = vEnd - vStart;
const int sizeE = totalNumVertices * _data->spaceDim;
scalar_array valsE(sizeE);
+ const PylithScalar dampingConstsScale = _data->densityScale * _data->lengthScale / _data->timeScale;
+ const PylithScalar jacobianScale = dampingConstsScale*pow(_data->lengthScale, _data->spaceDim-1);
+
const int spaceDim = _data->spaceDim;
for (int iVertex=0; iVertex < totalNumVertices; ++iVertex)
for (int iDim=0; iDim < spaceDim; ++iDim) {
@@ -349,7 +359,7 @@
#endif // DEBUGGING
PetscSection jacobianSection = jacobian.petscSection();
- Vec jacobianVec = jacobian.localVector();
+ PetscVec jacobianVec = jacobian.localVector();
PetscScalar *vals;
PetscInt size;
@@ -360,9 +370,9 @@
err = VecGetArray(jacobianVec, &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);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i]*jacobianScale, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i]/jacobianScale, vals[i], tolerance);
err = VecRestoreArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateJacobianLumped
@@ -387,6 +397,10 @@
// Set coordinate system
spatialdata::geocoords::CSCart 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);
@@ -413,6 +427,7 @@
bc->timeStep(_data->dt);
bc->label(_data->label);
bc->db(&db);
+ bc->normalizer(normalizer);
bc->createSubMesh(*mesh);
bc->initialize(*mesh, upDir);
@@ -437,18 +452,19 @@
residual.newSection(pylith::topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
residual.allocate();
residual.zero();
+ residual.scale(normalizer.lengthScale());
fields->copyLayout("residual");
const int totalNumVertices = vEnd - vStart;
const int fieldSize = _data->spaceDim * totalNumVertices;
PetscSection dispTIncrSection = fields->get("dispIncr(t->t+dt)").petscSection();
- Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
- PetscSection dispTSection = fields->get("disp(t)").petscSection();
- Vec dispTVec = fields->get("disp(t)").localVector();
+ PetscVec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ PetscSection dispTSection = fields->get("disp(t)").petscSection();
+ PetscVec dispTVec = fields->get("disp(t)").localVector();
PetscSection dispTmdtSection = fields->get("disp(t-dt)").petscSection();
- Vec dispTmdtVec = fields->get("disp(t-dt)").localVector();
- PetscSection velSection = fields->get("velocity(t)").petscSection();
- Vec velVec = fields->get("velocity(t)").localVector();
+ PetscVec dispTmdtVec = fields->get("disp(t-dt)").localVector();
+ PetscSection velSection = fields->get("velocity(t)").petscSection();
+ PetscVec velVec = fields->get("velocity(t)").localVector();
PetscScalar *dispTIncrArray, *dispTArray, *dispTmdtArray, *velArray;
const int spaceDim = _data->spaceDim;
const PylithScalar dt = _data->dt;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2013-03-07 00:43:38 UTC (rev 21456)
@@ -38,11 +38,6 @@
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::bc::TestDirichletBC );
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RealUniformSection RealUniformSection;
-
-// ----------------------------------------------------------------------
// Setup testing data.
void
pylith::bc::TestDirichletBC::setUp(void)
@@ -95,40 +90,43 @@
// Check values
CPPUNIT_ASSERT(0 != bc._parameters);
PetscSection initialSection = bc._parameters->get("initial").petscSection();
- Vec initialVec = bc._parameters->get("initial").localVector();
- PetscScalar *initialArray;
- PetscErrorCode err;
+ PetscVec initialVec = bc._parameters->get("initial").localVector();
+ PetscScalar* initialArray;
+ PetscErrorCode err = 0;
CPPUNIT_ASSERT(initialSection);CPPUNIT_ASSERT(initialVec);
const PylithScalar tolerance = 1.0e-06;
+ const PylithScalar dispScale = _data->lengthScale;
+ const PylithScalar velocityScale = _data->lengthScale / _data->timeScale;
err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
for (int i=0; i < numPoints; ++i) {
const PetscInt p_value = _data->constrainedPoints[i]+offset;
- PetscInt dof, off;
+ PetscInt dof, off;
err = PetscSectionGetDof(initialSection, p_value, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(initialSection, p_value, &off);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(numFixedDOF, dof);
- for(int iDOF = 0; iDOF < numFixedDOF; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->valuesInitial[i*numFixedDOF+iDOF], initialArray[off+iDOF], tolerance);
+ for(int iDOF = 0; iDOF < numFixedDOF; ++iDOF) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->valuesInitial[i*numFixedDOF+iDOF]/dispScale, initialArray[off+iDOF], 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();
+ PetscVec rateVec = bc._parameters->get("rate").localVector();
PetscScalar *rateArray;
err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
for (int i=0; i < numPoints; ++i) {
const PetscInt p_value = _data->constrainedPoints[i]+offset;
- PetscInt dof, off;
+ PetscInt dof, off;
err = PetscSectionGetDof(rateSection, p_value, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(rateSection, p_value, &off);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(numFixedDOF, dof);
for(int iDOF = 0; iDOF < numFixedDOF; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->valueRate, rateArray[off+iDOF], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->valueRate/velocityScale, rateArray[off+iDOF], tolerance);
} // for
err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
} // if
@@ -157,10 +155,10 @@
_initialize(&mesh, &bc);
CPPUNIT_ASSERT(0 != _data);
- DM dmMesh = mesh.dmMesh();
+ PetscDM dmMesh = mesh.dmMesh();
CPPUNIT_ASSERT(dmMesh);
- PetscInt cStart, cEnd, vStart, vEnd;
- PetscErrorCode err;
+ 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);
@@ -176,7 +174,7 @@
field.allocate();
PetscSection fieldSection = field.petscSection();
- Vec fieldVec = field.localVector();
+ PetscVec fieldVec = field.localVector();
CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
const PetscInt numCells = cEnd - cStart;
@@ -214,9 +212,9 @@
_initialize(&mesh, &bc);
CPPUNIT_ASSERT(0 != _data);
- DM dmMesh = mesh.dmMesh();
+ PetscDM dmMesh = mesh.dmMesh();
CPPUNIT_ASSERT(dmMesh);
- PetscInt cStart, cEnd, vStart, vEnd;
+ PetscInt cStart, cEnd, vStart, vEnd;
PetscErrorCode err;
err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
@@ -234,7 +232,7 @@
bc.setConstraints(field);
PetscSection fieldSection = field.petscSection();
- Vec fieldVec = field.localVector();
+ PetscVec fieldVec = field.localVector();
CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
const PetscInt numCells = cEnd - cStart;
@@ -242,7 +240,7 @@
int iConstraint = 0;
for(PetscInt v = vStart; v < vEnd; ++v) {
const PetscInt *cInd, *fcInd;
- PetscInt dof, cdof, fdof, fcdof;
+ PetscInt dof, cdof, fdof, fcdof;
err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
@@ -276,10 +274,10 @@
_initialize(&mesh, &bc);
CPPUNIT_ASSERT(0 != _data);
- DM dmMesh = mesh.dmMesh();
+ PetscDM dmMesh = mesh.dmMesh();
CPPUNIT_ASSERT(dmMesh);
- PetscInt cStart, cEnd, vStart, vEnd;
- PetscErrorCode err;
+ 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);
@@ -295,9 +293,12 @@
bc.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 dispScale = _data->lengthScale;
+ const PylithScalar velocityScale = _data->lengthScale / _data->timeScale;
+ const PylithScalar timeScale = _data->timeScale;
// All values should be zero.
PetscScalar *values;
@@ -331,8 +332,8 @@
} // for
assert(index == numFreeDOF);
- const PetscInt numCells = cEnd - cStart;
- const PetscInt offset = numCells;
+ const PetscInt numCells = cEnd - cStart;
+ const PetscInt offset = numCells;
const PetscInt numFixedDOF = _data->numFixedDOF;
int iConstraint = 0;
@@ -357,9 +358,9 @@
// check constrained DOF
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
const int index = iConstraint * numFixedDOF + iDOF;
- const PylithScalar valueE = (t > _data->tRef) ?
- _data->valuesInitial[index] + (t-_data->tRef)*_data->valueRate :
- _data->valuesInitial[index];
+ const PylithScalar valueE = (t > _data->tRef/timeScale) ?
+ _data->valuesInitial[index]/dispScale + (t-_data->tRef/timeScale)*_data->valueRate/velocityScale :
+ _data->valuesInitial[index]/dispScale;
CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[off+_data->fixedDOF[iDOF]], tolerance);
} // for
++iConstraint;
@@ -378,10 +379,10 @@
_initialize(&mesh, &bc);
CPPUNIT_ASSERT(0 != _data);
- DM dmMesh = mesh.dmMesh();
+ PetscDM dmMesh = mesh.dmMesh();
CPPUNIT_ASSERT(dmMesh);
- PetscInt cStart, cEnd, vStart, vEnd;
- PetscErrorCode err;
+ 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);
@@ -397,9 +398,12 @@
bc.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 dispScale = _data->lengthScale;
+ const PylithScalar velocityScale = _data->lengthScale / _data->timeScale;
+ const PylithScalar timeScale = _data->timeScale;
// All values should be zero.
PetscScalar *values;
@@ -435,8 +439,8 @@
} // for
assert(index == numFreeDOF);
- const PetscInt numCells = cEnd - cStart;
- const PetscInt offset = numCells;
+ const PetscInt numCells = cEnd - cStart;
+ const PetscInt offset = numCells;
const PetscInt numFixedDOF = _data->numFixedDOF;
int iConstraint = 0;
@@ -460,8 +464,8 @@
// check constrained DOF
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
- const PylithScalar valueE = (t0 > _data->tRef) ? (t1-t0)*_data->valueRate :
- (t1 > _data->tRef) ? (t1-_data->tRef)*_data->valueRate : 0.0;
+ 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;
CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[off+_data->fixedDOF[iDOF]], tolerance);
} // for
++iConstraint;
@@ -483,6 +487,10 @@
spatialdata::geocoords::CSCart 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);
@@ -520,6 +528,7 @@
bc->dbInitial(&db);
bc->dbRate(&dbRate);
bc->bcDOF(_data->fixedDOF, _data->numFixedDOF);
+ bc->normalizer(normalizer);
bc->initialize(*mesh, upDir);
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.cc 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.cc 2013-03-07 00:43:38 UTC (rev 21456)
@@ -22,6 +22,10 @@
// Constructor
pylith::bc::AbsorbingDampersData::AbsorbingDampersData(void) :
meshFilename(0),
+ lengthScale(1.0e+3),
+ pressureScale(2.25e+10),
+ densityScale(1.0),
+ timeScale(2.0),
numBasis(0),
numQuadPts(0),
quadPts(0),
@@ -45,6 +49,8 @@
valsResidual(0),
valsJacobian(0)
{ // constructor
+ const PylithScalar velScale = lengthScale / timeScale;
+ densityScale = pressureScale / (velScale*velScale);
} // constructor
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.hh 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/AbsorbingDampersData.hh 2013-03-07 00:43:38 UTC (rev 21456)
@@ -44,6 +44,14 @@
char* meshFilename; ///< Name of file with input mesh
+ /// @name Scales information for nondimensionalization.
+ //@{
+ PylithScalar lengthScale; ///< Length scale.
+ PylithScalar pressureScale; ///< Pressure scale.
+ PylithScalar timeScale; ///< Time scale.
+ PylithScalar densityScale; ///< Density scale.
+ //@}
+
/// @name Quadrature information
//@{
int numBasis; ///< Number of basis functions for cell
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletData.cc 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletData.cc 2013-03-07 00:43:38 UTC (rev 21456)
@@ -32,8 +32,14 @@
constrainedPoints(0),
valuesInitial(0),
meshFilename(0),
- dbFilename(0)
+ dbFilename(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/DirichletData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletData.hh 2013-03-06 23:57:22 UTC (rev 21455)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/data/DirichletData.hh 2013-03-07 00:43:38 UTC (rev 21456)
@@ -58,6 +58,15 @@
char* meshFilename; ///< Filename for input mesh.
char* dbFilename; ///< Filename of simple spatial database.
+
+ /// @name Scales information for nondimensionalization.
+ //@{
+ PylithScalar lengthScale; ///< Length scale.
+ PylithScalar pressureScale; ///< Pressure scale.
+ PylithScalar timeScale; ///< Time scale.
+ PylithScalar densityScale; ///< Density scale.
+ //@}
+
};
#endif // pylith_bc_dirichletdata_hh
More information about the CIG-COMMITS
mailing list