[cig-commits] r20886 - in short/3D/PyLith/trunk: libsrc/pylith/bc libsrc/pylith/topology unittests/libtests/bc
knepley at geodynamics.org
knepley at geodynamics.org
Tue Oct 23 07:46:46 PDT 2012
Author: knepley
Date: 2012-10-23 07:46:46 -0700 (Tue, 23 Oct 2012)
New Revision: 20886
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.hh
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/BoundaryConditionPoints.hh
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.hh
short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/FieldBase.hh
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/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependentPoints.cc
Log:
Updated BC to DMComplex, tests fail with value problems now, added FACES_FIELD for submeshes
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -20,7 +20,7 @@
#include "AbsorbingDampers.hh" // implementation of object methods
-#include "pylith/topology/FieldsNew.hh" // HOLDSA FieldsNew
+#include "pylith/topology/Fields.hh" // HOLDSA Fields
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
@@ -91,14 +91,13 @@
const int numCorners = _quadrature->numBasis();
// Get 'surface' cells (1 dimension lower than top-level cells)
- const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
- assert(!sieveSubMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- sieveSubMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
// Get damping constants at each quadrature point and rotate to
// global coordinate frame using orientation information
const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
@@ -111,11 +110,10 @@
logger.stagePush("BoundaryConditions");
delete _parameters;
- _parameters =
- new topology::FieldsNew<topology::SubMesh>(*_boundaryMesh);
+ _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
assert(0 != _parameters);
- _parameters->add("damping constants", "damping_constants", fiberDim);
- _parameters->allocate(cells);
+ _parameters->add("damping constants", "damping_constants", topology::FieldBase::FACES_FIELD, fiberDim);
+ _parameters->get("damping constants").allocate();
logger.stagePop();
@@ -145,15 +143,14 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Container for damping constants for current cell
- scalar_array dampingConstsLocal(fiberDim);
- scalar_array dampingConstsGlobal(fiberDim);
+ scalar_array dampingConstsLocal(spaceDim);
- scalar_array coordinatesCell(numCorners*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveSubMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ scalar_array coordinatesCell(numBasis*spaceDim);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -162,9 +159,10 @@
const PylithScalar velocityScale =
_normalizer->lengthScale() / _normalizer->timeScale();
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
+ PetscSection valueSection = _parameters->get("damping constants").petscSection();
+ Vec valueVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingArray;
+ assert(valueSection);assert(valueVec);
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
@@ -174,58 +172,58 @@
_quadrature->computeGeometry(*_boundaryMesh, cells);
#endif
- for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ err = VecGetArray(valueVec, &dampingArray);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
- coordsVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
+ PetscInt ddof, doff;
+ err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ assert(ddof == fiberDim);
+
const scalar_array& quadPtsNondim = _quadrature->quadPts();
const scalar_array& quadPtsRef = _quadrature->quadPtsRef();
quadPtsGlobal = quadPtsNondim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
- dampingConstsGlobal = 0.0;
for(int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
// Compute damping constants in normal/tangential coordinates
const int err = _db->query(&queryData[0], numValues,
- &quadPtsGlobal[iQuad*spaceDim], spaceDim, cs);
+ &quadPtsGlobal[iQuad*spaceDim], spaceDim, cs);
if (err) {
- std::ostringstream msg;
- msg << "Could not find parameters for physical properties at \n"
- << "(";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPtsGlobal[iQuad*spaceDim+i];
- msg << ") for absorbing boundary condition " << _label << "\n"
- << "using spatial database " << _db->label() << ".";
- throw std::runtime_error(msg.str());
+ std::ostringstream msg;
+ msg << "Could not find parameters for physical properties at \n"
+ << "(";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << quadPtsGlobal[iQuad*spaceDim+i];
+ msg << ") for absorbing boundary condition " << _label << "\n"
+ << "using spatial database " << _db->label() << ".";
+ throw std::runtime_error(msg.str());
} // if
// Nondimensionalize damping constants
- const PylithScalar densityN =
- _normalizer->nondimensionalize(queryData[0], densityScale);
- const PylithScalar vpN =
- _normalizer->nondimensionalize(queryData[1], velocityScale);
- const PylithScalar vsN = (3 == numValues) ?
- _normalizer->nondimensionalize(queryData[2], velocityScale) :
- 0.0;
+ const PylithScalar densityN = _normalizer->nondimensionalize(queryData[0], densityScale);
+ const PylithScalar vpN = _normalizer->nondimensionalize(queryData[1], velocityScale);
+ const PylithScalar vsN = (3 == numValues) ? _normalizer->nondimensionalize(queryData[2], velocityScale) : 0.0;
const PylithScalar constTangential = densityN * vsN;
const PylithScalar constNormal = densityN * vpN;
const int numTangential = spaceDim-1;
for (int iDim=0; iDim < numTangential; ++iDim)
- dampingConstsLocal[iDim] = constTangential;
+ dampingConstsLocal[iDim] = constTangential;
dampingConstsLocal[spaceDim-1] = constNormal;
// Compute normal/tangential orientation
- memcpy(&quadPtRef[0], &quadPtsRef[iQuad*cellDim],
- cellDim*sizeof(PylithScalar));
+ memcpy(&quadPtRef[0], &quadPtsRef[iQuad*cellDim], cellDim*sizeof(PylithScalar));
#if defined(PRECOMPUTE_GEOMETRY)
coordsVisitor.clear();
sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
@@ -236,16 +234,15 @@
orientation /= jacobianDet;
for (int iDim=0; iDim < spaceDim; ++iDim) {
- for (int jDim=0; jDim < spaceDim; ++jDim)
- dampingConstsGlobal[iQuad*spaceDim+iDim] +=
- dampingConstsLocal[jDim]*orientation[jDim*spaceDim+iDim];
- // Ensure damping constants are positive
- dampingConstsGlobal[iQuad*spaceDim+iDim] =
- fabs(dampingConstsGlobal[iQuad*spaceDim+iDim]);
+ dampingArray[doff+iQuad*spaceDim+iDim] = 0.0;
+ for (int jDim=0; jDim < spaceDim; ++jDim)
+ dampingArray[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]);
} // for
} // for
- parametersSection->updatePoint(*c_iter, &dampingConstsGlobal[0]);
} // for
+ err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
_db->close();
} // initialize
@@ -283,37 +280,39 @@
_initCellVector();
// Get cell information
- const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
- assert(!sieveSubMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- sieveSubMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ IS subpointMap;
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(subMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+
// Get sections
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
+ PetscSection valueSection = _parameters->get("damping constants").petscSection();
+ Vec valueVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingArray;
+ assert(valueSection);assert(valueVec);
// Use _cellVector for cell residual.
- PetscSection residualSection = residual.petscSection();
+ PetscSection residualSection = residual.petscSection(), residualSubsection;
Vec residualVec = residual.localVector();
assert(residualSection);assert(residualVec);
- ///UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+ err = PetscSectionCreateSubmeshSection(residualSection, subpointMap, &residualSubsection);
- scalar_array velCell(numBasis*spaceDim);
- PetscSection velSection = fields->get("velocity(t)").petscSection();
+ PetscSection velSection = fields->get("velocity(t)").petscSection(), velSubsection;
Vec velVec = fields->get("velocity(t)").localVector();
assert(velSection);assert(velVec);
- ///RestrictVisitor velVisitor(*velSection, velCell.size(), &velCell[0]);
+ err = PetscSectionCreateSubmeshSection(velSection, subpointMap, &velSubsection);
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveSubMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
#endif
_logger->eventEnd(setupEvent);
@@ -321,9 +320,8 @@
_logger->eventBegin(computeEvent);
#endif
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
// Get geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
@@ -331,9 +329,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- coordsVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -345,12 +346,16 @@
_resetCellVector();
// Restrict input fields to cell
- ///velVisitor.clear();
- ///sieveSubMesh->restrictClosure(*c_iter, velVisitor);
- assert(numQuadPts*spaceDim ==
- parametersSection->getFiberDimension(*c_iter));
- const PylithScalar* dampersCell = parametersSection->restrictPoint(*c_iter);
+ const PetscScalar *velArray = PETSC_NULL;
+ PetscInt velSize;
+ PetscInt ddof, doff;
+ err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ assert(ddof == numQuadPts*spaceDim);
+ err = DMComplexVecGetClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == numBasis*spaceDim);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -370,11 +375,13 @@
const PylithScalar valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -=
- dampersCell[iQuad*spaceDim+iDim] *
- valIJ * velCell[jBasis*spaceDim+iDim];
+ dampingArray[doff+iQuad*spaceDim+iDim] *
+ valIJ * velArray[jBasis*spaceDim+iDim];
} // for
} // for
} // for
+ err = DMComplexVecRestoreClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(subMesh, residualSubsection, residualVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
@@ -382,17 +389,16 @@
_logger->eventBegin(updateEvent);
#endif
- // Assemble cell contribution into field
- ///residualVisitor.clear();
- ///sieveSubMesh->updateClosure(*c_iter, residualVisitor);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = PetscSectionDestroy(&residualSubsection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionDestroy(&velSubsection);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
+ PetscLogFlops((cEnd-cStart)*numQuadPts*(1+numBasis*(1+numBasis*(3*spaceDim))));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -430,36 +436,39 @@
_initCellVector();
// Get cell information
- const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
- assert(!sieveSubMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- sieveSubMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ IS subpointMap;
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(subMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+
// Get sections
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
+ PetscSection valueSection = _parameters->get("damping constants").petscSection();
+ Vec valueVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingArray;
+ assert(valueSection);assert(valueVec);
// Use _cellVector for cell values.
- const ALE::Obj<RealSection>& residualSection = residual.section();
- assert(!residualSection.isNull());
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+ PetscSection residualSection = residual.petscSection(), residualSubsection;
+ Vec residualVec = residual.localVector();
+ assert(residualSection);assert(residualVec);
+ err = PetscSectionCreateSubmeshSection(residualSection, subpointMap, &residualSubsection);
- scalar_array velCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& velSection =
- fields->get("velocity(t)").section();
- assert(!velSection.isNull());
- RestrictVisitor velVisitor(*velSection, velCell.size(), &velCell[0]);
+ PetscSection velSection = fields->get("velocity(t)").petscSection(), velSubsection;
+ Vec velVec = fields->get("velocity(t)").localVector();
+ assert(velSection);assert(velVec);
+ err = PetscSectionCreateSubmeshSection(velSection, subpointMap, &velSubsection);
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveSubMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
#endif
_logger->eventEnd(setupEvent);
@@ -467,9 +476,8 @@
_logger->eventBegin(computeEvent);
#endif
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
// Get geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
@@ -477,9 +485,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- coordsVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -491,12 +502,16 @@
_resetCellVector();
// Restrict input fields to cell
- velVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, velVisitor);
- assert(numQuadPts*spaceDim ==
- parametersSection->getFiberDimension(*c_iter));
- const PylithScalar* dampersCell = parametersSection->restrictPoint(*c_iter);
+ const PetscScalar *velArray = PETSC_NULL;
+ PetscInt velSize;
+ PetscInt ddof, doff;
+ err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ assert(ddof == numQuadPts*spaceDim);
+ err = DMComplexVecGetClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == numBasis*spaceDim);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -517,26 +532,29 @@
for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
const PylithScalar valIJ = basis[iQ + iBasis] * valJ;
for (int iDim = 0; iDim < spaceDim; ++iDim)
- _cellVector[iBasis*spaceDim+iDim] -= valIJ *
- dampersCell[iQuad*spaceDim+iDim] * velCell[iBasis*spaceDim+iDim];
+ _cellVector[iBasis*spaceDim+iDim] -=
+ dampingArray[doff+iQuad*spaceDim+iDim] *
+ valIJ * velArray[iBasis*spaceDim+iDim];
} // for
} // for
+ err = DMComplexVecRestoreClosure(subMesh, velSubsection, velVec, c, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(subMesh, residualSubsection, residualVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3)));
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
#endif
- // Assemble cell contribution into field
- residualVisitor.clear();
- sieveSubMesh->updateClosure(*c_iter, residualVisitor);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = PetscSectionDestroy(&residualSubsection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionDestroy(&velSubsection);CHECK_PETSC_ERROR(err);
+
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3)));
+ PetscLogFlops((cEnd-cStart)*numQuadPts*(1+numBasis+numBasis*(1+spaceDim*3)));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidualLumped
@@ -571,35 +589,27 @@
const int spaceDim = _quadrature->spaceDim();
// Get cell information
- const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
- assert(!sieveSubMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- sieveSubMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ IS subpointMap;
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(subMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+
// Get sections
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
+ PetscSection valueSection = _parameters->get("damping constants").petscSection();
+ Vec valueVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingArray;
+ assert(valueSection);assert(valueVec);
const topology::Field<topology::Mesh>& solution = fields->solution();
- const ALE::Obj<SieveMesh>& sieveMesh = solution.mesh().sieveMesh();
- const ALE::Obj<RealSection>& solutionSection = solution.section();
- assert(!solutionSection.isNull());
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- solutionSection);
- assert(!globalOrder.isNull());
+ PetscSection solutionSection = solution.petscSection(), solutionSubsection;
+ Vec solutionVec = solution.localVector();
+ assert(solutionSection);assert(solutionVec);
+ err = PetscSectionCreateSubmeshSection(solutionSection, subpointMap, &solutionSubsection);
- const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
- assert(!sieve.isNull());
- const int closureSize =
- int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
- IndicesVisitor jacobianVisitor(*solutionSection, *globalOrder,
- closureSize*spaceDim);
-
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
assert(0 != jacobianMat);
@@ -613,10 +623,11 @@
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveSubMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
#endif
_logger->eventEnd(setupEvent);
@@ -624,9 +635,8 @@
_logger->eventBegin(computeEvent);
#endif
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
@@ -634,9 +644,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- coordsVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -644,9 +657,12 @@
#endif
// Get damping constants
- assert(numQuadPts*spaceDim == parametersSection->getFiberDimension(*c_iter));
- const PylithScalar* dampingConstsCell = parametersSection->restrictPoint(*c_iter);
+ PetscInt ddof, doff;
+ err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ assert(ddof == numQuadPts*spaceDim);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -662,7 +678,7 @@
// Compute Jacobian for absorbing bc terms
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const PylithScalar wt =
- quadWts[iQuad] * jacobianDet[iQuad] / (2.0 * dt);
+ quadWts[iQuad] * jacobianDet[iQuad] / (2.0 * dt);
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const PylithScalar valI = wt*basis[iQ+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
@@ -670,8 +686,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 * dampingConstsCell[iQuad*spaceDim+iDim];
+ _cellMatrix[iBlock+jBlock] += valIJ * dampingArray[doff+iQuad*spaceDim+iDim];
} // for
} // for
} // for
@@ -683,19 +698,18 @@
#endif
// Assemble cell contribution into PETSc Matrix
- jacobianVisitor.clear();
- PetscErrorCode err = updateOperator(jacobianMat, *sieveSubMesh->getSieve(),
- jacobianVisitor, *c_iter,
- &_cellMatrix[0], ADD_VALUES);
+ err = DMComplexMatSetClosure(subMesh, solutionSubsection, PETSC_NULL, jacobianMat, c, &_cellMatrix[0], ADD_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = PetscSectionDestroy(&solutionSubsection);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(1+2*spaceDim))));
+ PetscLogFlops((cEnd-cStart)*numQuadPts*(3+numBasis*(1+numBasis*(1+2*spaceDim))));
_logger->eventEnd(computeEvent);
#endif
@@ -732,14 +746,15 @@
const int spaceDim = _quadrature->spaceDim();
// Get cell information
- const ALE::Obj<SieveSubMesh>& sieveSubMesh = _boundaryMesh->sieveMesh();
- assert(!sieveSubMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- sieveSubMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ IS subpointMap;
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(subMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+
// Get parameters used in integration.
const PylithScalar dt = _dt;
assert(dt > 0);
@@ -749,26 +764,23 @@
_initCellVector();
// Get sections
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
+ PetscSection valueSection = _parameters->get("damping constants").petscSection();
+ Vec valueVec = _parameters->get("damping constants").localVector();
+ PetscScalar *dampingArray;
+ assert(valueSection);assert(valueVec);
- const topology::Field<topology::Mesh>& solution = fields->solution();
- const ALE::Obj<SieveMesh>& sieveMesh = solution.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<RealSection>& solutionSection = solution.section();
- assert(!solutionSection.isNull());
+ PetscSection jacobianSection = jacobian->petscSection(), jacobianSubsection;
+ Vec jacobianVec = jacobian->localVector();
+ assert(jacobianSection);assert(jacobianVec);
+ err = PetscSectionCreateSubmeshSection(jacobianSection, subpointMap, &jacobianSubsection);
- const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
- assert(!jacobianSection.isNull());
- UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
-
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveSubMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
#endif
_logger->eventEnd(setupEvent);
@@ -776,9 +788,8 @@
_logger->eventBegin(computeEvent);
#endif
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ err = VecGetArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
@@ -786,9 +797,12 @@
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- coordsVisitor.clear();
- sieveSubMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -796,9 +810,12 @@
#endif
// Get damping constants
- assert(numQuadPts*spaceDim == parametersSection->getFiberDimension(*c_iter));
- const PylithScalar* dampingConstsCell = parametersSection->restrictPoint(*c_iter);
+ PetscInt ddof, doff;
+ err = PetscSectionGetDof(valueSection, c, &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, c, &doff);CHECK_PETSC_ERROR(err);
+ assert(ddof == numQuadPts*spaceDim);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -823,25 +840,24 @@
const PylithScalar valIJ = basis[iQ + iBasis] * valJ;
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis * spaceDim + iDim] += valIJ
- * dampingConstsCell[iQuad * spaceDim + iDim];
+ * dampingArray[doff+iQuad * spaceDim + iDim];
} // for
} // for
+ err = DMComplexVecSetClosure(subMesh, jacobianSubsection, jacobianVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(4+numBasis+numBasis*(1+spaceDim*2)));
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
#endif
-
- // Assemble cell contribution into lumped matrix.
- jacobianVisitor.clear();
- sieveSubMesh->updateClosure(*c_iter, jacobianVisitor);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = VecRestoreArray(valueVec, &dampingArray);CHECK_PETSC_ERROR(err);
+ err = PetscSectionDestroy(&jacobianSubsection);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(4+numBasis+numBasis*(1+spaceDim*2)));
+ PetscLogFlops((cEnd-cStart)*numQuadPts*(4+numBasis+numBasis*(1+spaceDim*2)));
_logger->eventEnd(computeEvent);
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -21,7 +21,6 @@
#include "BCIntegratorSubMesh.hh" // implementation of object methods
#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
// ----------------------------------------------------------------------
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.hh 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.hh 2012-10-23 14:46:46 UTC (rev 20886)
@@ -63,7 +63,7 @@
*
* @returns Parameter fields.
*/
- const topology::FieldsNew<topology::SubMesh>*
+ const topology::Fields<topology::Field<topology::SubMesh> >*
parameterFields(void) const;
/** Get boundary mesh.
@@ -90,7 +90,7 @@
topology::SubMesh* _boundaryMesh; ///< Boundary mesh.
/// Parameters for boundary condition.
- topology::FieldsNew<topology::SubMesh>* _parameters;
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.icc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BCIntegratorSubMesh.icc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -24,7 +24,7 @@
// Get parameter fields.
inline
-const pylith::topology::FieldsNew<pylith::topology::SubMesh>*
+const pylith::topology::Fields<pylith::topology::Field<pylith::topology::SubMesh> >*
pylith::bc::BCIntegratorSubMesh::parameterFields(void) const {
return _parameters;
} // parameterFields
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -21,7 +21,7 @@
#include "BoundaryConditionPoints.hh" // implementation of object methods
#include "pylith/topology/Mesh.hh" // USES Mesh
-#include "pylith/topology/FieldsNew.hh" // HOLDSA FieldsNew
+#include "pylith/topology/Fields.hh" // HOLDSA Fields
#include "pylith/topology/Field.hh" // USES Field
#include <stdexcept> // USES std::runtime_error()
@@ -54,7 +54,7 @@
// ----------------------------------------------------------------------
// Get parameter fields.
-const pylith::topology::FieldsNew<pylith::topology::Mesh>*
+const pylith::topology::Fields<pylith::topology::Field<pylith::topology::Mesh> >*
pylith::bc::BoundaryConditionPoints::parameterFields(void) const
{ // parameterFields
return _parameters;
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.hh 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.hh 2012-10-23 14:46:46 UTC (rev 20886)
@@ -30,7 +30,7 @@
// Include directives ---------------------------------------------------
#include "BoundaryCondition.hh" // ISA BoundaryCondition
-#include "pylith/topology/topologyfwd.hh" // HOLDSA FieldsNew
+#include "pylith/topology/topologyfwd.hh" // HOLDSA Fields
#include "pylith/utils/array.hh" // HASA int_array
// BoundaryConditionPoints ----------------------------------------------
@@ -59,7 +59,7 @@
*
* @returns Parameter fields.
*/
- const topology::FieldsNew<topology::Mesh>* parameterFields(void) const;
+ const topology::Fields<topology::Field<topology::Mesh> >* parameterFields(void) const;
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -74,7 +74,7 @@
protected :
/// Parameters for boundary condition.
- topology::FieldsNew<topology::Mesh>* _parameters;
+ topology::Fields<topology::Field<topology::Mesh> >* _parameters;
int_array _points; ///< Points for boundary condition.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -22,7 +22,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -193,35 +193,33 @@
const int numPoints = _points.size();
assert(_parameters);
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(valueFiberDim == numFixedDOF);
- assert(valueIndex+valueFiberDim <= parametersFiberDim);
+ PetscSection valueSection = _parameters->get("value").petscSection();
+ Vec valueVec = _parameters->get("value").localVector();
+ PetscScalar *valueArray;
+ assert(valueSection);assert(valueVec);
- PetscSection fieldSectionP = field.petscSection();
- Vec localVec = field.localVector();
+ PetscSection fieldSection = field.petscSection();
+ Vec localVec = field.localVector();
PetscScalar *array;
PetscErrorCode err;
- assert(fieldSectionP);
+ assert(fieldSection);assert(localVec);
err = VecGetArray(localVec, &array);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
PetscInt p_bc = _points[iPoint];
- PetscInt dof, off;
+ PetscInt dof, off, vdof, voff;
- assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
- err = PetscSectionGetDof(fieldSectionP, p_bc, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(fieldSectionP, p_bc, &off);CHECK_PETSC_ERROR(err);
- const PylithScalar *parametersVertex = parametersSection->restrictPoint(p_bc);
-
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- array[_bcDOF[iDOF]+off] = parametersVertex[valueIndex+iDOF];
+ 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];
} // for
err = VecRestoreArray(localVec, &array);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
} // setField
// ----------------------------------------------------------------------
@@ -241,35 +239,33 @@
const int numPoints = _points.size();
assert(_parameters);
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(valueFiberDim == numFixedDOF);
- assert(valueIndex+valueFiberDim <= parametersFiberDim);
+ PetscSection valueSection = _parameters->get("value").petscSection();
+ Vec valueVec = _parameters->get("value").localVector();
+ PetscScalar *valueArray;
+ assert(valueSection);assert(valueVec);
- PetscSection fieldSectionP = field.petscSection();
- Vec localVec = field.localVector();
+ PetscSection fieldSection = field.petscSection();
+ Vec localVec = field.localVector();
PetscScalar *array;
PetscErrorCode err;
- assert(fieldSectionP);
+ assert(fieldSection);assert(localVec);
err = VecGetArray(localVec, &array);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
PetscInt p_bc = _points[iPoint];
- PetscInt dof, off;
+ PetscInt dof, off, vdof, voff;
- assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
- err = PetscSectionGetDof(fieldSectionP, p_bc, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(fieldSectionP, p_bc, &off);CHECK_PETSC_ERROR(err);
- const PylithScalar* parametersVertex = parametersSection->restrictPoint(p_bc);
-
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- array[_bcDOF[iDOF]+off] = parametersVertex[valueIndex+iDOF];
+ 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];
} // for
err = VecRestoreArray(localVec, &array);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 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -23,7 +23,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -90,17 +90,16 @@
logger.stagePush("BoundaryConditions");
if (0 == _outputFields)
- _outputFields = new topology::FieldsNew<topology::SubMesh>(*_boundaryMesh);
+ _outputFields = new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
assert(0 != _outputFields);
- _outputFields->add("buffer (vector)", "buffer_vector",
- spaceDim,
- topology::FieldBase::VECTOR,
- lengthScale);
- _outputFields->add("buffer (scalar)", "buffer_scalar",
- 1,
- topology::FieldBase::SCALAR,
- timeScale);
- _outputFields->allocate(topology::FieldBase::CELLS_FIELD, 1);
+ _outputFields->add("buffer (vector)", "buffer_vector", topology::FieldBase::CELLS_FIELD, spaceDim);
+ _outputFields->get("buffer (vector)").vectorFieldType(topology::FieldBase::VECTOR);
+ _outputFields->get("buffer (vector)").scale(lengthScale);
+ _outputFields->get("buffer (vector)").allocate();
+ _outputFields->add("buffer (scalar)", "buffer_scalar", topology::FieldBase::CELLS_FIELD, 1);
+ _outputFields->get("buffer (scalar)").vectorFieldType(topology::FieldBase::SCALAR);
+ _outputFields->get("buffer (scalar)").scale(timeScale);
+ _outputFields->get("buffer (scalar)").allocate();
logger.stagePop();
@@ -149,9 +148,6 @@
throw std::runtime_error(msg.str());
} // if
- const ALE::Obj<SieveMesh>& sieveMesh = _boundaryMesh->sieveMesh();
- assert(!sieveMesh.isNull());
-
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
assert(0 != cs);
const int spaceDim = cs->spaceDim();
@@ -160,42 +156,36 @@
const int numFixedDOF = _bcDOF.size();
assert(_outputFields->hasField("buffer (vector)"));
- const ALE::Obj<SubRealUniformSection>& outputSection =
- _outputFields->section();
- assert(!outputSection.isNull());
- const int outputFiberDim = _outputFields->fiberDim();
- scalar_array outputVertex(outputFiberDim);
- const int bufferIndex = _outputFields->sectionIndex("buffer (vector)");
- const int bufferFiberDim = _outputFields->sectionFiberDim("buffer (vector)");
- assert(bufferIndex + bufferFiberDim <= outputFiberDim);
- assert(spaceDim == bufferFiberDim);
+ 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);
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- const int fieldIndex = _parameters->sectionIndex(name);
- const int fieldFiberDim = _parameters->sectionFiberDim(name);
- assert(fieldIndex + fieldFiberDim <= parametersFiberDim);
- assert(fieldFiberDim == numFixedDOF);
-
+ err = VecGetArray(outputVec, &outputArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
- outputVertex = 0.0;
+ PetscInt odof, ooff, fdof, foff;
- assert(parametersFiberDim == parametersSection->getFiberDimension(point));
- const PylithScalar* parametersVertex = parametersSection->restrictPoint(point);
- assert(parametersVertex);
-
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- outputVertex[bufferIndex+_bcDOF[iDOF]] =
- parametersVertex[fieldIndex+iDOF];
- assert(outputFiberDim == outputSection->getFiberDimension(point));
- outputSection->updatePointAll(point, &outputVertex[0]);
+ 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
+ err = VecRestoreArray(outputVec, &outputArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
- topology::Field<topology::SubMesh>& buffer =
- _outputFields->get("buffer (vector)");
+ topology::Field<topology::SubMesh>& buffer = _outputFields->get("buffer (vector)");
buffer.label(label);
buffer.scale(scale);
@@ -225,46 +215,36 @@
throw std::runtime_error(msg.str());
} // if
- const ALE::Obj<SieveMesh>& sieveMesh = _boundaryMesh->sieveMesh();
- assert(!sieveMesh.isNull());
-
const int numPoints = _points.size();
assert(_outputFields->hasField("buffer (scalar)"));
- const ALE::Obj<SubRealUniformSection>& outputSection =
- _outputFields->section();
- assert(!outputSection.isNull());
- const int outputFiberDim = _outputFields->fiberDim();
- scalar_array outputVertex(outputFiberDim);
- const int bufferIndex = _outputFields->sectionIndex("buffer (vector)");
- const int bufferFiberDim = _outputFields->sectionFiberDim("buffer (vector)");
- assert(bufferIndex + bufferFiberDim <= outputFiberDim);
- assert(1 == bufferFiberDim);
+ PetscSection outputSection = _outputFields->get("buffer (scalar)").petscSection();
+ Vec outputVec = _outputFields->get("buffer (scalar)").localVector();
+ PetscScalar *outputArray;
+ PetscErrorCode err;
+ assert(outputSection);assert(outputVec);
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- const int fieldIndex = _parameters->sectionIndex(name);
- const int fieldFiberDim = _parameters->sectionFiberDim(name);
- assert(fieldIndex + fieldFiberDim <= parametersFiberDim);
- assert(1 == fieldFiberDim);
+ PetscSection fieldSection = _parameters->get(name).petscSection();
+ Vec fieldVec = _parameters->get(name).localVector();
+ PetscScalar *fieldArray;
+ assert(fieldSection);assert(fieldVec);
+ err = VecGetArray(outputVec, &outputArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
- outputVertex = 0.0;
+ PetscInt odof, ooff, fdof, foff;
- assert(parametersFiberDim == parametersSection->getFiberDimension(point));
- const PylithScalar* parametersVertex = parametersSection->restrictPoint(point);
- assert(parametersVertex);
-
- outputVertex[bufferIndex] = parametersVertex[fieldIndex];
- assert(outputFiberDim == outputSection->getFiberDimension(point));
- outputSection->updatePointAll(point, &outputVertex[0]);
+ 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
- topology::Field<topology::SubMesh>& buffer =
- _outputFields->get("buffer (scalar)");
+ topology::Field<topology::SubMesh>& buffer = _outputFields->get("buffer (scalar)");
buffer.label(label);
buffer.scale(scale);
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.hh 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBoundary.hh 2012-10-23 14:46:46 UTC (rev 20886)
@@ -107,7 +107,7 @@
topology::SubMesh* _boundaryMesh; ///< Boundary mesh.
/// Fields manager (holds temporary field for output).
- topology::FieldsNew<topology::SubMesh>* _outputFields;
+ topology::Fields<topology::Field<topology::SubMesh> >* _outputFields;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -21,7 +21,7 @@
#include "Neumann.hh" // implementation of object methods
#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/topology/FieldsNew.hh" // HOLDSA FieldsNew
+#include "pylith/topology/Fields.hh" // HOLDSA Fields
#include "pylith/topology/Field.hh" // USES Field
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
@@ -106,58 +106,60 @@
scalar_array tractionsCell(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveSubMesh>& subSieveMesh = _boundaryMesh->sieveMesh();
- assert(!subSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- subSieveMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ IS subpointMap;
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(subMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+
// Get sections
_calculateValue(t);
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(valueFiberDim == tractionsCell.size());
- assert(valueIndex+valueFiberDim <= parametersFiberDim);
+ PetscSection valueSection = _parameters->get("value").petscSection();
+ Vec valueVec = _parameters->get("value").localVector();
+ PetscScalar *tractionsArray;
+ assert(valueSection);assert(valueVec);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+ PetscSection residualSection = residual.petscSection(), residualSubsection;
+ Vec residualVec = residual.localVector();
+ assert(residualSection);assert(residualVec);
+ err = PetscSectionCreateSubmeshSection(residualSection, subpointMap, &residualSubsection);
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- subSieveMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
#endif
// Loop over faces and integrate contribution from each face
- for (SieveSubMesh::label_sequence::iterator c_iter=cells->begin();
- c_iter != cellsEnd;
- ++c_iter) {
+ err = VecGetArray(valueVec, &tractionsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(c);
#else
- coordsVisitor.clear();
- subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Reset element vector to zero
_resetCellVector();
// Restrict tractions to cell
- assert(parametersFiberDim ==
- parametersSection->getFiberDimension(*c_iter));
- const PylithScalar* parametersCell = parametersSection->restrictPoint(*c_iter);
- assert(parametersCell);
- const PylithScalar* tractionsCell = ¶metersCell[valueIndex];
- assert(tractionsCell);
+ PetscInt vdof, voff;
+ err = PetscSectionGetDof(valueSection, c, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, c, &voff);CHECK_PETSC_ERROR(err);
+ assert(vdof == numQuadPts*spaceDim);
+
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
const scalar_array& jacobianDet = _quadrature->jacobianDet();
@@ -171,16 +173,14 @@
const PylithScalar valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] +=
- tractionsCell[iQuad*spaceDim+iDim] * valIJ;
+ tractionsArray[voff+iQuad*spaceDim+iDim] * valIJ;
} // for
} // for
} // for
- // Assemble cell contribution into field
- residualVisitor.clear();
- subSieveMesh->updateClosure(*c_iter, residualVisitor);
-
+ err = DMComplexVecSetClosure(subMesh, residualSubsection, residualVec, c, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
PetscLogFlops(numQuadPts*(1+numBasis*(1+numBasis*(1+2*spaceDim))));
} // for
+ err = VecRestoreArray(valueVec, &tractionsArray);CHECK_PETSC_ERROR(err);
} // integrateResidual
// ----------------------------------------------------------------------
@@ -245,37 +245,39 @@
const int numQuadPts = _quadrature->numQuadPts();
delete _parameters;
- _parameters =
- new topology::FieldsNew<topology::SubMesh>(*_boundaryMesh);
+ _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(*_boundaryMesh);
// Create section to hold time dependent values
- _parameters->add("value", "traction",
- numQuadPts*spaceDim, topology::FieldBase::MULTI_VECTOR,
- pressureScale);
- if (_dbInitial)
- _parameters->add("initial", "initial_traction",
- numQuadPts*spaceDim, topology::FieldBase::MULTI_VECTOR,
- pressureScale);
+ _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();
+ 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();
+ }
if (_dbRate) {
- _parameters->add("rate", "traction_rate",
- numQuadPts*spaceDim, topology::FieldBase::MULTI_VECTOR,
- rateScale);
- _parameters->add("rate time", "traction_rate_time",
- numQuadPts, topology::FieldBase::MULTI_SCALAR,
- timeScale);
+ _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();
+ _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();
} // if
if (_dbChange) {
- _parameters->add("change", "change_traction",
- numQuadPts*spaceDim, topology::FieldBase::MULTI_VECTOR,
- pressureScale);
- _parameters->add("change time", "change_traction_time",
- numQuadPts, topology::FieldBase::MULTI_SCALAR,
- timeScale);
+ _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();
+ _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();
} // if
- _parameters->allocate(topology::FieldBase::CELLS_FIELD, 1);
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
if (0 != _dbInitial) { // Setup initial values, if provided.
_dbInitial->open();
@@ -396,14 +398,13 @@
assert(_parameters);
// Get 'surface' cells (1 dimension lower than top-level cells)
- const ALE::Obj<SieveSubMesh>& subSieveMesh = _boundaryMesh->sieveMesh();
- assert(!subSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- subSieveMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
const int numBasis = _quadrature->numBasis();
const int numQuadPts = _quadrature->numQuadPts();
@@ -416,20 +417,16 @@
// Get sections.
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- subSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- const int valueIndex = _parameters->sectionIndex(name);
- const int valueFiberDim = _parameters->sectionFiberDim(name);
- assert(valueIndex+valueFiberDim <= parametersFiberDim);
- scalar_array parametersCell(parametersFiberDim);
+ PetscSection valueSection = _parameters->get(name).petscSection();
+ Vec valueVec = _parameters->get(name).localVector();
+ PetscScalar *valueArray;
+ assert(valueSection);assert(valueVec);
const spatialdata::geocoords::CoordSys* cs = _boundaryMesh->coordsys();
assert(cs);
@@ -444,52 +441,49 @@
#endif
// Loop over cells in boundary mesh and perform queries.
- for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ err = VecGetArray(valueVec, &valueArray);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
- coordsVisitor.clear();
- subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
const scalar_array& quadPtsNondim = _quadrature->quadPts();
quadPtsGlobal = quadPtsNondim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
valuesCell = 0.0;
- for (int iQuad=0, iSpace=0;
- iQuad < numQuadPts;
- ++iQuad, iSpace+=spaceDim) {
+ for (int iQuad=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iSpace+=spaceDim) {
const int err = db->query(&valuesCell[iQuad*querySize], querySize,
- &quadPtsGlobal[iSpace], spaceDim, cs);
+ &quadPtsGlobal[iSpace], spaceDim, cs);
if (err) {
- std::ostringstream msg;
- msg << "Could not find values at (";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPtsGlobal[i+iSpace];
- msg << ") for traction boundary condition " << _label << "\n"
- << "using spatial database " << db->label() << ".";
- throw std::runtime_error(msg.str());
+ std::ostringstream msg;
+ msg << "Could not find values at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << quadPtsGlobal[i+iSpace];
+ msg << ") for traction boundary condition " << _label << "\n"
+ << "using spatial database " << db->label() << ".";
+ throw std::runtime_error(msg.str());
} // if
-
} // for
- _normalizer->nondimensionalize(&valuesCell[0], valuesCell.size(),
- scale);
+ _normalizer->nondimensionalize(&valuesCell[0], valuesCell.size(), scale);
// Update section
- assert(parametersFiberDim == parametersSection->getFiberDimension(*c_iter));
- parametersSection->restrictPoint(*c_iter,
- ¶metersCell[0], parametersCell.size());
- for (int i=0; i < valueFiberDim; ++i)
- parametersCell[valueIndex+i] = valuesCell[i];
-
- parametersSection->updatePoint(*c_iter, ¶metersCell[0]);
+ PetscInt vdof, voff;
+
+ err = PetscSectionGetDof(valueSection, c, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, c, &voff);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < vdof; ++d)
+ valueArray[voff+d] = valuesCell[d];
} // for
+ err = VecRestoreArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
} // _queryDB
// ----------------------------------------------------------------------
@@ -507,14 +501,13 @@
up[i] = upDir[i];
// Get 'surface' cells (1 dimension lower than top-level cells)
- const ALE::Obj<SieveSubMesh>& subSieveMesh = _boundaryMesh->sieveMesh();
- assert(!subSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- subSieveMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
// Quadrature related values.
const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
@@ -533,118 +526,139 @@
// Get sections.
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- subSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(subMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(subMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- scalar_array parametersCellLocal(parametersFiberDim);
- scalar_array parametersCellGlobal(parametersFiberDim);
+ scalar_array parametersCellLocal(spaceDim);
+ PetscSection initialSection, rateSection, rateTimeSection, changeSection, changeTimeSection;
+ Vec initialVec, rateVec, rateTimeVec, changeVec, changeTimeVec;
+ PetscScalar *valuesArray, *initialArray, *rateArray, *rateTimeArray, *changeArray, *changeTimeArray;
- const int initialIndex =
- (_dbInitial) ? _parameters->sectionIndex("initial") : -1;
- const int initialFiberDim =
- (_dbInitial) ? _parameters->sectionFiberDim("initial") : 0;
+ PetscSection valuesSection = _parameters->get("value").petscSection();
+ Vec valuesVec = _parameters->get("value").localVector();
+ assert(valuesSection);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);
+ }
+ 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);
+ }
+ 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);
+ err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
- const int rateIndex = (_dbRate) ? _parameters->sectionIndex("rate") : -1;
- const int rateFiberDim =
- (_dbRate) ? _parameters->sectionFiberDim("rate") : 0;
-
- const int changeIndex = (_dbChange) ? _parameters->sectionIndex("change") : -1;
- const int changeFiberDim =
- (_dbChange) ? _parameters->sectionFiberDim("change") : 0;
-
// Loop over cells in boundary mesh, compute orientations, and then
// rotate corresponding traction vector from local to global coordinates.
- for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- coordsVisitor.clear();
- subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(subMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
- // Get parameters
- parametersSection->restrictPoint(*c_iter,
- ¶metersCellLocal[0],
- parametersCellLocal.size());
- // Copy parameter values to global (not all are tranformed)
- parametersCellGlobal = parametersCellLocal;
-
- for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts;
- ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
- // Compute Jacobian and determinant at quadrature point, then get
- // orientation.
+ for(int iQuad=0, iRef=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
+ // Compute Jacobian and determinant at quadrature point, then get orientation.
memcpy(&quadPtRef[0], &quadPtsRef[iRef], cellDim*sizeof(PylithScalar));
#if defined(PRECOMPUTE_GEOMETRY)
coordsVisitor.clear();
subSieveMesh->restrictClosure(*c_iter, coordsVisitor);
#endif
- cellGeometry.jacobian(&jacobian, &jacobianDet,
- coordinatesCell, quadPtRef);
+ cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell, quadPtRef);
cellGeometry.orientation(&orientation, jacobian, jacobianDet, up);
assert(jacobianDet > 0.0);
orientation /= jacobianDet;
if (0 != _dbInitial) {
- // Rotate traction vector from local coordinate system to global
- // coordinate system
- assert(initialIndex >= 0);
- assert(initialFiberDim == numQuadPts*spaceDim);
- PylithScalar* initialGlobal = ¶metersCellGlobal[initialIndex+iSpace];
- PylithScalar* initialLocal = ¶metersCellLocal[initialIndex+iSpace];
- for(int iDim = 0; iDim < spaceDim; ++iDim) {
- initialGlobal[iDim] = 0.0;
- for(int jDim = 0; jDim < spaceDim; ++jDim)
- initialGlobal[iDim] +=
- orientation[jDim*spaceDim+iDim] * initialLocal[jDim];
- } // for
+ // Rotate traction vector from local coordinate system to global coordinate system
+ PylithScalar *initialLocal = ¶metersCellLocal[0];
+ 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(int iDim = 0; iDim < spaceDim; ++iDim) {
+ initialArray[ioff+iSpace+iDim] = 0.0;
+ for(int jDim = 0; jDim < spaceDim; ++jDim)
+ initialArray[ioff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * initialLocal[jDim];
+ } // for
} // if
if (0 != _dbRate) {
- // Rotate traction vector from local coordinate system to global
- // coordinate system
- assert(rateIndex >= 0);
- assert(rateFiberDim == numQuadPts*spaceDim);
- PylithScalar* rateGlobal = ¶metersCellGlobal[rateIndex+iSpace];
- PylithScalar* rateLocal = ¶metersCellLocal[rateIndex+iSpace];
- for(int iDim = 0; iDim < spaceDim; ++iDim) {
- rateGlobal[iDim] = 0.0;
- for(int jDim = 0; jDim < spaceDim; ++jDim)
- rateGlobal[iDim] +=
- orientation[jDim*spaceDim+iDim] * rateLocal[jDim];
- } // for
+ // Rotate traction vector from local coordinate system to global coordinate system
+ PylithScalar *rateLocal = ¶metersCellLocal[0];
+ 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(int iDim = 0; iDim < spaceDim; ++iDim) {
+ rateArray[roff+iSpace+iDim] = 0.0;
+ for(int jDim = 0; jDim < spaceDim; ++jDim)
+ rateArray[roff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * rateLocal[jDim];
+ } // for
} // if
if (0 != _dbChange) {
- // Rotate traction vector from local coordinate system to global
- // coordinate system
- assert(changeIndex >= 0);
- assert(changeFiberDim == numQuadPts*spaceDim);
- PylithScalar* changeGlobal = ¶metersCellGlobal[changeIndex+iSpace];
- PylithScalar* changeLocal = ¶metersCellLocal[changeIndex+iSpace];
- for(int iDim = 0; iDim < spaceDim; ++iDim) {
- changeGlobal[iDim] = 0.0;
- for(int jDim = 0; jDim < spaceDim; ++jDim)
- changeGlobal[iDim] +=
- orientation[jDim*spaceDim+iDim] * changeLocal[jDim];
- } // for
+ // Rotate traction vector from local coordinate system to global coordinate system
+ PylithScalar *changeLocal = ¶metersCellLocal[0];
+ 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) {
+ changeLocal[iDim] = changeArray[coff+iSpace+iDim];
+ }
+ for(int iDim = 0; iDim < spaceDim; ++iDim) {
+ changeArray[coff+iSpace+iDim] = 0.0;
+ for(int jDim = 0; jDim < spaceDim; ++jDim)
+ changeArray[coff+iSpace+iDim] += orientation[jDim*spaceDim+iDim] * changeLocal[jDim];
+ } // for
} // if
-
} // for
-
- // Update sections
- assert(parametersFiberDim == parametersSection->getFiberDimension(*c_iter));
- parametersSection->updatePoint(*c_iter, ¶metersCellGlobal[0]);
} // for
+ err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ if (_dbInitial) {
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbRate) {
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbChange) {
+ err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
} // paramsLocalToGlobal
// ----------------------------------------------------------------------
@@ -659,116 +673,132 @@
const PylithScalar timeScale = _getNormalizer().timeScale();
// Get 'surface' cells (1 dimension lower than top-level cells)
- const ALE::Obj<SieveSubMesh>& subSieveMesh = _boundaryMesh->sieveMesh();
- assert(!subSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- subSieveMesh->heightStratum(1);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM subMesh = _boundaryMesh->dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(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;
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- scalar_array parametersCell(parametersFiberDim);
-
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(numQuadPts*spaceDim == valueFiberDim);
+ PetscSection valuesSection = _parameters->get("value").petscSection();
+ Vec valuesVec = _parameters->get("value").localVector();
+ assert(valuesSection);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);
+ }
+ 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);
+ }
+ 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);
+ err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
- const int initialIndex =
- (_dbInitial) ? _parameters->sectionIndex("initial") : -1;
- const int initialFiberDim =
- (_dbInitial) ? _parameters->sectionFiberDim("initial") : 0;
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ PetscInt vdof, voff;
- const int rateIndex =
- (_dbRate) ? _parameters->sectionIndex("rate") : -1;
- const int rateFiberDim =
- (_dbRate) ? _parameters->sectionFiberDim("rate") : 0;
- const int rateTimeIndex =
- (_dbRate) ? _parameters->sectionIndex("rate time") : -1;
- const int rateTimeFiberDim =
- (_dbRate) ? _parameters->sectionFiberDim("rate time") : 0;
+ 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)
+ valuesArray[voff+d] = 0.0;
- const int changeIndex =
- (_dbChange) ? _parameters->sectionIndex("change") : -1;
- const int changeFiberDim =
- (_dbChange) ? _parameters->sectionFiberDim("change") : 0;
- const int changeTimeIndex =
- (_dbChange) ? _parameters->sectionIndex("change time") : -1;
- const int changeTimeFiberDim =
- (_dbChange) ? _parameters->sectionFiberDim("change time") : 0;
-
- for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
-
- assert(parametersFiberDim == parametersSection->getFiberDimension(*c_iter));
- parametersSection->restrictPoint(*c_iter,
- ¶metersCell[0], parametersCell.size());
- for (int i=0; i < valueFiberDim; ++i)
- parametersCell[valueIndex+i] = 0.0;
-
// Contribution from initial value
if (0 != _dbInitial) {
- assert(initialIndex >= 0);
- assert(initialFiberDim == valueFiberDim);
- for (int i=0; i < initialFiberDim; ++i)
- parametersCell[valueIndex+i] += parametersCell[initialIndex+i];
+ 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)
+ valuesArray[voff+d] += initialArray[ioff+d];
} // if
// Contribution from rate of change of value
if (0 != _dbRate) {
- assert(rateIndex >= 0);
- assert(rateFiberDim == valueFiberDim);
- assert(rateTimeIndex >= 0);
- assert(rateTimeFiberDim == numQuadPts);
-
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const PylithScalar tRel = t - parametersCell[rateTimeIndex+iQuad];
- if (tRel > 0.0) // rate of change integrated over time
- for (int iDim=0; iDim < spaceDim; ++iDim)
- parametersCell[valueIndex+iQuad*spaceDim+iDim] +=
- parametersCell[rateIndex+iQuad*spaceDim+iDim] * tRel;
+ 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);
+ err = PetscSectionGetOffset(rateTimeSection, c, &rtoff);CHECK_PETSC_ERROR(err);
+ assert(rdof == vdof);
+ assert(rtdof == numQuadPts);
+
+ 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)
+ valuesArray[voff+iQuad*spaceDim+iDim] += rateArray[roff+iQuad*spaceDim+iDim] * tRel;
} // for
} // if
// Contribution from change of value
if (0 != _dbChange) {
- assert(changeIndex >= 0);
- assert(changeFiberDim == valueFiberDim);
- assert(changeTimeIndex >= 0);
- assert(changeTimeFiberDim == numQuadPts);
+ PetscInt cdof, coff, ctdof, ctoff;
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const PylithScalar tRel = t - parametersCell[changeTimeIndex+iQuad];
- if (tRel >= 0) { // change in value over time
- PylithScalar scale = 1.0;
- if (0 != _dbTimeHistory) {
- PylithScalar tDim = tRel;
- _getNormalizer().dimensionalize(&tDim, 1, timeScale);
- const int err = _dbTimeHistory->query(&scale, tDim);
- if (0 != err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- } // if
- for (int iDim=0; iDim < spaceDim; ++iDim)
- parametersCell[valueIndex+iQuad*spaceDim+iDim] +=
- parametersCell[changeIndex+iQuad*spaceDim+iDim] * scale;
- } // if
+ 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);
+ err = PetscSectionGetOffset(changeTimeSection, c, &ctoff);CHECK_PETSC_ERROR(err);
+ assert(cdof == vdof);
+ assert(ctdof == numQuadPts);
+
+ for(int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
+ const PylithScalar tRel = t - changeTimeArray[ctoff+iQuad];
+ if (tRel >= 0) { // change in value over time
+ PylithScalar scale = 1.0;
+ if (0 != _dbTimeHistory) {
+ PylithScalar tDim = tRel;
+ _getNormalizer().dimensionalize(&tDim, 1, timeScale);
+ const int err = _dbTimeHistory->query(&scale, tDim);
+ if (0 != err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ for(int iDim = 0; iDim < spaceDim; ++iDim)
+ valuesArray[voff+iQuad*spaceDim+iDim] += changeArray[coff+iQuad*spaceDim+iDim] * scale;
+ } // if
} // for
} // if
-
- parametersSection->updatePoint(*c_iter, ¶metersCell[0]);
} // for
+ err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ if (_dbInitial) {
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbRate) {
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbChange) {
+ err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
} // _calculateValue
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/PointForce.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -21,7 +21,7 @@
#include "PointForce.hh" // implementation of object methods
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
@@ -99,46 +99,48 @@
assert(0 != cs);
const int spaceDim = cs->spaceDim();
- scalar_array residualVertex(spaceDim);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- assert(!residualSection.isNull());
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *residualArray;
+ assert(residualSection);assert(residualVec);
// Get global order
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- residualSection);
- assert(!globalOrder.isNull());
+ DM dmMesh = residual.dmMesh();
+ PetscSection globalSection;
+ PetscErrorCode err;
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(valueIndex+valueFiberDim <= parametersFiberDim);
- assert(valueFiberDim == numBCDOF);
+ 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);
+
+ err = VecGetArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const int p_bc = _points[iPoint]; // Get point label.
+ const PetscInt p_bc = _points[iPoint]; // Get point label.
+ PetscInt goff;
// Contribute to residual only if point is local.
- if (!globalOrder->isLocal(p_bc))
- continue;
+ err = PetscSectionGetOffset(globalSection, p_bc, &goff);CHECK_PETSC_ERROR(err);
+ if (goff < 0) continue;
- residualVertex *= 0.0; // Reset residual contribution to zero.
-
- assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
- const PylithScalar* parametersVertex = parametersSection->restrictPoint(p_bc);
- assert(parametersVertex);
+ PetscInt vdof, voff, rdof, roff;
+ err = PetscSectionGetDof(valueSection, p_bc, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, p_bc, &voff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(residualSection, p_bc, &rdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(residualSection, p_bc, &roff);CHECK_PETSC_ERROR(err);
+ assert(vdof == numBCDOF);
+ assert(rdof == spaceDim);
+
for (int iDOF=0; iDOF < numBCDOF; ++iDOF)
- residualVertex[_bcDOF[iDOF]] += parametersVertex[valueIndex+iDOF];
-
- assert(residualVertex.size() == residualSection->getFiberDimension(p_bc));
- residualSection->updateAddPoint(p_bc, &residualVertex[0]);
+ residualArray[roff+_bcDOF[iDOF]] += valueArray[voff+iDOF];
} // for
+ err = VecRestoreArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
} // integrateResidualAssembled
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/TimeDependentPoints.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -22,7 +22,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
#include "spatialdata/spatialdb/TimeHistory.hh" // USES TimeHistory
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -83,6 +83,8 @@
{ // _queryDatabases
const PylithScalar timeScale = _getNormalizer().timeScale();
const PylithScalar rateScale = valueScale / timeScale;
+ Vec v;
+ PetscErrorCode err;
const int numPoints = _points.size();
const int numBCDOF = _bcDOF.size();
@@ -119,46 +121,56 @@
logger.stagePush("BoundaryConditions");
delete _parameters;
- _parameters = new topology::FieldsNew<topology::Mesh>(mesh);
+ _parameters = new topology::Fields<topology::Field<topology::Mesh> >(mesh);
- _parameters->add("value", "value", numBCDOF, topology::FieldBase::OTHER,
- valueScale);
+ _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);
if (_dbInitial) {
- const std::string& fieldLabel =
- std::string("initial_") + std::string(fieldName);
- _parameters->add("initial", fieldLabel.c_str(),
- numBCDOF, topology::FieldBase::OTHER,
- valueScale);
+ 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);
} // if
if (_dbRate) {
- const std::string& fieldLabel =
- std::string("rate_") + std::string(fieldName);
- _parameters->add("rate", fieldLabel.c_str(),
- numBCDOF, topology::FieldBase::OTHER,
- rateScale);
- const std::string& timeLabel =
- std::string("rate_time_") + std::string(fieldName);
- _parameters->add("rate time", timeLabel.c_str(),
- 1, topology::FieldBase::SCALAR,
- timeScale);
+ 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);
+ 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);
} // if
if (_dbChange) {
- const std::string& fieldLabel =
- std::string("change_") + std::string(fieldName);
- _parameters->add("change", fieldLabel.c_str(),
- numBCDOF, topology::FieldBase::OTHER,
- valueScale);
- const std::string& timeLabel =
- std::string("change_time_") + std::string(fieldName);
- _parameters->add("change time", timeLabel.c_str(),
- 1, topology::FieldBase::SCALAR,
- timeScale);
+ 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);
+ 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);
} // if
- _parameters->allocate(_points);
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
if (_dbInitial) { // Setup initial values, if provided.
_dbInitial->open();
@@ -229,18 +241,23 @@
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
+ 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
scalar_array valueVertex(querySize);
const int numPoints = _points.size();
+ PetscScalar *array;
+ PetscErrorCode err;
+ err = VecGetArray(parametersVec, &array);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
// Get dimensionalized coordinates of vertex
coordinates->restrictPoint(_points[iPoint],
@@ -261,15 +278,15 @@
scale);
// Update section
- assert(parametersFiberDim ==
- parametersSection->getFiberDimension(_points[iPoint]));
- parametersSection->restrictPoint(_points[iPoint], ¶metersVertex[0],
- parametersVertex.size());
- for (int i=0; i < valueFiberDim; ++i)
- parametersVertex[valueIndex+i] = valueVertex[i];
-
- parametersSection->updatePoint(_points[iPoint], ¶metersVertex[0]);
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(parametersSection, _points[iPoint], &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(parametersSection, _points[iPoint], &off);CHECK_PETSC_ERROR(err);
+ //assert(parametersFiberDim == dof);
+ for(int i = 0; i < dof; ++i)
+ array[off+i] = valueVertex[i];
} // for
+ err = VecRestoreArray(parametersVec, &array);CHECK_PETSC_ERROR(err);
} // _queryDB
// ----------------------------------------------------------------------
@@ -282,101 +299,115 @@
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;
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- scalar_array parametersVertex(parametersFiberDim);
-
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(numBCDOF == valueFiberDim);
-
- const int initialIndex =
- (_dbInitial) ? _parameters->sectionIndex("initial") : -1;
- const int initialFiberDim =
- (_dbInitial) ? _parameters->sectionFiberDim("initial") : 0;
+ PetscSection valuesSection = _parameters->get("value").petscSection();
+ Vec valuesVec = _parameters->get("value").localVector();
+ assert(valuesSection);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);
+ }
+ 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);
+ }
+ 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);
+ err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
- const int rateIndex =
- (_dbRate) ? _parameters->sectionIndex("rate") : -1;
- const int rateFiberDim =
- (_dbRate) ? _parameters->sectionFiberDim("rate") : 0;
- const int rateTimeIndex =
- (_dbRate) ? _parameters->sectionIndex("rate time") : -1;
- const int rateTimeFiberDim =
- (_dbRate) ? _parameters->sectionFiberDim("rate time") : 0;
-
- const int changeIndex =
- (_dbChange) ? _parameters->sectionIndex("change") : -1;
- const int changeFiberDim =
- (_dbChange) ? _parameters->sectionFiberDim("change") : 0;
- const int changeTimeIndex =
- (_dbChange) ? _parameters->sectionIndex("change time") : -1;
- const int changeTimeFiberDim =
- (_dbChange) ? _parameters->sectionFiberDim("change time") : 0;
-
- for (int iPoint=0; iPoint < numPoints; ++iPoint) {
+ for(int iPoint=0; iPoint < numPoints; ++iPoint) {
const int p_bc = _points[iPoint]; // Get point label
-
- assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
- parametersSection->restrictPoint(p_bc, ¶metersVertex[0],
- parametersVertex.size());
- for (int i=0; i < valueFiberDim; ++i)
- parametersVertex[valueIndex+i] = 0.0;
+ PetscInt vdof, voff;
+ err = PetscSectionGetDof(valuesSection, p_bc, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valuesSection, p_bc, &voff);CHECK_PETSC_ERROR(err);
// Contribution from initial value
if (_dbInitial) {
- assert(initialIndex >= 0);
- assert(initialFiberDim == valueFiberDim);
- for (int i=0; i < initialFiberDim; ++i)
- parametersVertex[valueIndex+i] += parametersVertex[initialIndex+i];
+ PetscInt idof, ioff;
+
+ err = PetscSectionGetDof(initialSection, p_bc, &idof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(initialSection, p_bc, &ioff);CHECK_PETSC_ERROR(err);
+ assert(idof == vdof);
+ for(PetscInt d = 0; d < vdof; ++d)
+ valuesArray[voff+d] += initialArray[ioff+d];
} // if
// Contribution from rate of change of value
if (_dbRate) {
- assert(rateIndex >= 0);
- assert(rateFiberDim == valueFiberDim);
- assert(rateTimeIndex >= 0);
- assert(rateTimeFiberDim == 1);
-
- const PylithScalar tRel = t - parametersVertex[rateTimeIndex];
+ PetscInt rdof, roff, rtdof, rtoff;
+
+ err = PetscSectionGetDof(rateSection, p_bc, &rdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(rateSection, p_bc, &roff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(rateTimeSection, p_bc, &rtdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(rateTimeSection, p_bc, &rtoff);CHECK_PETSC_ERROR(err);
+ assert(rdof == vdof);
+ assert(rtdof == 1);
+ const PylithScalar tRel = t - rateTimeArray[rtoff];
if (tRel > 0.0) // rate of change integrated over time
- for (int iDim=0; iDim < numBCDOF; ++iDim)
- parametersVertex[valueIndex+iDim] +=
- parametersVertex[rateIndex+iDim] * tRel;
+ for(PetscInt d = 0; d < vdof; ++d)
+ valuesArray[voff+d] += rateArray[roff+d] * tRel;
} // if
// Contribution from change of value
if (_dbChange) {
- assert(changeIndex >= 0);
- assert(changeFiberDim == valueFiberDim);
- assert(changeTimeIndex >= 0);
- assert(changeTimeFiberDim == 1);
+ PetscInt cdof, coff, ctdof, ctoff;
- const PylithScalar tRel = t - parametersVertex[changeTimeIndex];
+ err = PetscSectionGetDof(changeSection, p_bc, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(changeSection, p_bc, &coff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(changeTimeSection, p_bc, &ctdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(changeTimeSection, p_bc, &ctoff);CHECK_PETSC_ERROR(err);
+ assert(cdof == vdof);
+ assert(ctdof == 1);
+ const PylithScalar tRel = t - changeTimeArray[ctoff];
if (tRel >= 0) { // change in value over time
- PylithScalar scale = 1.0;
- if (0 != _dbTimeHistory) {
- PylithScalar tDim = tRel;
- _getNormalizer().dimensionalize(&tDim, 1, timeScale);
- const int err = _dbTimeHistory->query(&scale, tDim);
- if (0 != err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- } // if
- for (int iDim=0; iDim < numBCDOF; ++iDim)
- parametersVertex[valueIndex+iDim] +=
- parametersVertex[changeIndex+iDim] * scale;
+ PylithScalar scale = 1.0;
+ if (0 != _dbTimeHistory) {
+ PylithScalar tDim = tRel;
+ _getNormalizer().dimensionalize(&tDim, 1, timeScale);
+ const int err = _dbTimeHistory->query(&scale, tDim);
+ if (0 != err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ for(PetscInt d = 0; d < vdof; ++d)
+ valuesArray[voff+d] += changeArray[coff+d] * scale;
} // if
} // if
-
- parametersSection->updatePoint(p_bc, ¶metersVertex[0]);
} // for
+ err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ if (_dbInitial) {
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbRate) {
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbChange) {
+ err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
} // _calculateValue
// ----------------------------------------------------------------------
@@ -391,127 +422,142 @@
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;
- const ALE::Obj<RealUniformSection>& parametersSection =
- _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- scalar_array parametersVertex(parametersFiberDim);
-
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(numBCDOF == valueFiberDim);
-
- const int rateIndex =
- (_dbRate) ? _parameters->sectionIndex("rate") : -1;
- const int rateFiberDim =
- (_dbRate) ? _parameters->sectionFiberDim("rate") : 0;
- const int rateTimeIndex =
- (_dbRate) ? _parameters->sectionIndex("rate time") : -1;
- const int rateTimeFiberDim =
- (_dbRate) ? _parameters->sectionFiberDim("rate time") : 0;
+ PetscSection valuesSection = _parameters->get("value").petscSection();
+ Vec valuesVec = _parameters->get("value").localVector();
+ assert(valuesSection);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);
+ }
+ 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);
+ }
+ 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);
+ err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
- const int changeIndex =
- (_dbChange) ? _parameters->sectionIndex("change") : -1;
- const int changeFiberDim =
- (_dbChange) ? _parameters->sectionFiberDim("change") : 0;
- const int changeTimeIndex =
- (_dbChange) ? _parameters->sectionIndex("change time") : -1;
- const int changeTimeFiberDim =
- (_dbChange) ? _parameters->sectionFiberDim("change time") : 0;
-
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const int p_bc = _points[iPoint]; // Get point label
-
- assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
- parametersSection->restrictPoint(p_bc, ¶metersVertex[0],
- parametersVertex.size());
- for (int i=0; i < valueFiberDim; ++i)
- parametersVertex[valueIndex+i] = 0.0;
+ PetscInt vdof, voff;
+ err = PetscSectionGetDof(valuesSection, p_bc, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valuesSection, p_bc, &voff);CHECK_PETSC_ERROR(err);
// No contribution from initial value
// Contribution from rate of change of value
if (_dbRate) {
- assert(rateIndex >= 0);
- assert(rateFiberDim == valueFiberDim);
- assert(rateTimeIndex >= 0);
- assert(rateTimeFiberDim == 1);
-
+ PetscInt rdof, roff, rtdof, rtoff;
+
+ err = PetscSectionGetDof(rateSection, p_bc, &rdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(rateSection, p_bc, &roff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(rateTimeSection, p_bc, &rtdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(rateTimeSection, p_bc, &rtoff);CHECK_PETSC_ERROR(err);
+ assert(rdof == vdof);
+ assert(rtdof == 1);
// Account for when rate dependence begins.
- const PylithScalar tRate = parametersVertex[rateTimeIndex];
+ const PylithScalar tRate = rateTimeArray[rtoff];
PylithScalar tIncr = 0.0;
if (t0 > tRate) // rate dependence for t0 to t1
- tIncr = t1 - t0;
+ tIncr = t1 - t0;
else if (t1 > tRate) // rate dependence for tRef to t1
- tIncr = t1 - tRate;
+ tIncr = t1 - tRate;
else
- tIncr = 0.0; // no rate dependence for t0 to t1
+ tIncr = 0.0; // no rate dependence for t0 to t1
if (tIncr > 0.0) // rate of change integrated over time
- for (int iDim=0; iDim < numBCDOF; ++iDim)
- parametersVertex[valueIndex+iDim] +=
- parametersVertex[rateIndex+iDim] * tIncr;
+ for(PetscInt d = 0; d < vdof; ++d)
+ valuesArray[voff+d] += rateArray[roff+d] * tIncr;
} // if
// Contribution from change of value
if (_dbChange) {
- assert(changeIndex >= 0);
- assert(changeFiberDim == valueFiberDim);
- assert(changeTimeIndex >= 0);
- assert(changeTimeFiberDim == 1);
+ PetscInt cdof, coff, ctdof, ctoff;
- const PylithScalar tChange = parametersVertex[changeTimeIndex];
+ err = PetscSectionGetDof(changeSection, p_bc, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(changeSection, p_bc, &coff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(changeTimeSection, p_bc, &ctdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(changeTimeSection, p_bc, &ctoff);CHECK_PETSC_ERROR(err);
+ assert(cdof == vdof);
+ assert(ctdof == 1);
+ const PylithScalar tChange = changeTimeArray[ctoff];
if (t0 >= tChange) { // increment is after change starts
- PylithScalar scale0 = 1.0;
- PylithScalar scale1 = 1.0;
- if (0 != _dbTimeHistory) {
- PylithScalar tDim = t0 - tChange;
- _getNormalizer().dimensionalize(&tDim, 1, timeScale);
- int err = _dbTimeHistory->query(&scale0, tDim);
- if (0 != err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- tDim = t1 - tChange;
- _getNormalizer().dimensionalize(&tDim, 1, timeScale);
- err = _dbTimeHistory->query(&scale1, tDim);
- if (0 != err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- } // if
- for (int iDim=0; iDim < numBCDOF; ++iDim)
- parametersVertex[valueIndex+iDim] +=
- parametersVertex[changeIndex+iDim] * (scale1 - scale0);
+ PylithScalar scale0 = 1.0;
+ PylithScalar scale1 = 1.0;
+ if (0 != _dbTimeHistory) {
+ PylithScalar tDim = t0 - tChange;
+ _getNormalizer().dimensionalize(&tDim, 1, timeScale);
+ int err = _dbTimeHistory->query(&scale0, tDim);
+ if (0 != err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ tDim = t1 - tChange;
+ _getNormalizer().dimensionalize(&tDim, 1, timeScale);
+ err = _dbTimeHistory->query(&scale1, tDim);
+ if (0 != err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ for(PetscInt d = 0; d < vdof; ++d)
+ valuesArray[voff+d] += changeArray[coff+d] * (scale1 - scale0);
} else if (t1 >= tChange) { // increment spans when change starts
- PylithScalar scale1 = 1.0;
- if (0 != _dbTimeHistory) {
- PylithScalar tDim = t1 - tChange;
- _getNormalizer().dimensionalize(&tDim, 1, timeScale);
- int err = _dbTimeHistory->query(&scale1, tDim);
- if (0 != err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- } // if
- for (int iDim=0; iDim < numBCDOF; ++iDim)
- parametersVertex[valueIndex+iDim] +=
- parametersVertex[changeIndex+iDim] * scale1;
+ PylithScalar scale1 = 1.0;
+ if (0 != _dbTimeHistory) {
+ PylithScalar tDim = t1 - tChange;
+ _getNormalizer().dimensionalize(&tDim, 1, timeScale);
+ int err = _dbTimeHistory->query(&scale1, tDim);
+ if (0 != err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ for(PetscInt d = 0; d < vdof; ++d)
+ valuesArray[voff+d] += changeArray[coff+d] * scale1;
} // if/else
} // if
-
- parametersSection->updatePoint(p_bc, ¶metersVertex[0]);
} // for
+ err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ if (_dbInitial) {
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbRate) {
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbChange) {
+ err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
} // _calculateValueIncr
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -388,6 +388,9 @@
case CELLS_FIELD:
err = DMComplexGetHeightStratum(_dm, stratum, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
break;
+ case FACES_FIELD:
+ err = DMComplexGetHeightStratum(_dm, stratum+1, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
case POINTS_FIELD:
err = DMComplexGetChart(_dm, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
break;
@@ -410,6 +413,8 @@
points = sieveMesh->depthStratum(stratum);
else if (CELLS_FIELD == domain)
points = sieveMesh->heightStratum(stratum);
+ else if (FACES_FIELD == domain)
+ points = sieveMesh->heightStratum(stratum+1);
else {
std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
assert(0);
@@ -1808,6 +1813,9 @@
case CELLS_FIELD:
err = DMComplexGetHeightStratum(_dm, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
break;
+ case FACES_FIELD:
+ err = DMComplexGetHeightStratum(_dm, 1, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
case POINTS_FIELD:
err = DMComplexGetChart(_dm, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
break;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/FieldBase.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/FieldBase.hh 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/FieldBase.hh 2012-10-23 14:46:46 UTC (rev 20886)
@@ -56,6 +56,7 @@
VERTICES_FIELD=0, ///< FieldBase over vertices.
CELLS_FIELD=1, ///< FieldBase over cells.
POINTS_FIELD=2, ///< FieldBase over all points.
+ FACES_FIELD=3, ///< FieldBase over faces.
}; // DomainEnum
// PUBLIC STRUCTS ///////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -27,7 +27,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
@@ -101,64 +101,66 @@
CPPUNIT_ASSERT(0 != _data);
const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
+ DM subMesh = boundaryMesh.dmMesh();
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
// Check boundary mesh
- CPPUNIT_ASSERT(!submesh.isNull());
+ CPPUNIT_ASSERT(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(subMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
const int cellDim = boundaryMesh.dimension();
const int numCorners = _data->numCorners;
const int spaceDim = _data->spaceDim;
- const ALE::Obj<SieveSubMesh::label_sequence>& cells = submesh->heightStratum(1);
- const int numVertices = submesh->depthStratum(0)->size();
- const int numCells = cells->size();
- const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
+ const int numVertices = vEnd-vStart;
+ const int numCells = cEnd-cStart;
+ //const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
CPPUNIT_ASSERT_EQUAL(_data->cellDim, cellDim);
CPPUNIT_ASSERT_EQUAL(_data->numVertices, numVertices);
CPPUNIT_ASSERT_EQUAL(_data->numCells, numCells);
- const ALE::Obj<SieveMesh::sieve_type>& sieve = submesh->getSieve();
- ALE::ISieveVisitor::PointRetriever<SieveSubMesh::sieve_type> pV(sieve->getMaxConeSize());
- int dp = 0;
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter) {
- const int numCorners = submesh->getNumCellCorners(*c_iter, boundaryDepth);
+ PetscInt dp = 0;
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ const PetscInt *cone;
+ PetscInt numCorners;
+
+ // Assume non-interpolated mesh
+ err = DMComplexGetConeSize(subMesh, c, &numCorners);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
-
- sieve->cone(*c_iter, pV);
- const SieveSubMesh::point_type *cone = pV.getPoints();
- for(int p = 0; p < pV.getSize(); ++p, ++dp) {
+ err = DMComplexGetCone(subMesh, c, &cone);CHECK_PETSC_ERROR(err);
+ for(PetscInt p = 0; p < numCorners; ++p, ++dp) {
CPPUNIT_ASSERT_EQUAL(_data->cells[dp], cone[p]);
}
- pV.clear();
} // for
// Check damping constants
const int numQuadPts = _data->numQuadPts;
const int fiberDim = numQuadPts * spaceDim;
- scalar_array dampersCell(fiberDim);
- int index = 0;
+ PetscInt index = 0;
CPPUNIT_ASSERT(0 != bc._parameters);
- const ALE::Obj<SubRealUniformSection>& dampersSection =
- bc._parameters->section();
+ PetscSection dampersSection = bc._parameters->get("damping constants").petscSection();
+ Vec dampersVec = bc._parameters->get("damping constants").localVector();
+ CPPUNIT_ASSERT(dampersSection);CPPUNIT_ASSERT(dampersVec);
const PylithScalar tolerance = 1.0e-06;
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter) {
- dampersSection->restrictPoint(*c_iter,
- &dampersCell[0], dampersCell.size());
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
- for (int iDim =0; iDim < spaceDim; ++iDim) {
- const PylithScalar dampersCellData = _data->dampingConsts[index];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0,
- dampersCell[iQuad*spaceDim+iDim]/dampersCellData,
- tolerance);
- ++index;
+ PetscScalar *dampersValues;
+ err = VecGetArray(dampersVec, &dampersValues);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);
+ 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);
+ ++index;
} // for
} // for
+ err = VecRestoreArray(dampersVec, &dampersValues);CHECK_PETSC_ERROR(err);
} // testInitialize
// ----------------------------------------------------------------------
@@ -174,7 +176,8 @@
_initialize(&mesh, &bc, &fields);
const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
+ DM subMesh = boundaryMesh.dmMesh();
+ PetscErrorCode err;
topology::Field<topology::Mesh>& residual = fields.get("residual");
const PylithScalar t = 0.0;
@@ -183,7 +186,6 @@
DM dmMesh = mesh.dmMesh();
PetscInt vStart, vEnd;
const PylithScalar* valsE = _data->valsResidual;
- PetscErrorCode err;
CPPUNIT_ASSERT(dmMesh);
err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
@@ -218,13 +220,14 @@
{ // testIntegrateJacobian
CPPUNIT_ASSERT(0 != _data);
+ CPPUNIT_ASSERT(0);
topology::Mesh mesh;
AbsorbingDampers bc;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &bc, &fields);
const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
+ DM subMesh = boundaryMesh.dmMesh();
topology::Field<topology::Mesh>& solution = fields.solution();
topology::Jacobian jacobian(solution);
@@ -302,8 +305,8 @@
jacobian.allocate();
const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
- CPPUNIT_ASSERT(!submesh.isNull());
+ DM subMesh = boundaryMesh.dmMesh();
+ CPPUNIT_ASSERT(subMesh);
const PylithScalar t = 1.0;
bc.integrateJacobian(&jacobian, t, &fields);
@@ -420,8 +423,6 @@
fields->add("velocity(t)", "velocity");
fields->solutionName("dispIncr(t->t+dt)");
- const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
-
topology::Field<topology::Mesh>& residual = fields->get("residual");
DM dmMesh = mesh->dmMesh();
PetscInt vStart, vEnd;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -26,7 +26,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
@@ -94,45 +94,43 @@
if (numFixedDOF > 0) {
// Check values
CPPUNIT_ASSERT(0 != bc._parameters);
- const ALE::Obj<RealUniformSection>& parametersSection =
- bc._parameters->section();
- CPPUNIT_ASSERT(!parametersSection.isNull());
- const int parametersFiberDim = bc._parameters->fiberDim();
- const int initialIndex = bc._parameters->sectionIndex("initial");
- const int initialFiberDim = bc._parameters->sectionFiberDim("initial");
- CPPUNIT_ASSERT_EQUAL(numFixedDOF, initialFiberDim);
+ PetscSection initialSection = bc._parameters->get("initial").petscSection();
+ Vec initialVec = bc._parameters->get("initial").localVector();
+ PetscScalar *initialArray;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(initialSection);CPPUNIT_ASSERT(initialVec);
const PylithScalar tolerance = 1.0e-06;
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
for (int i=0; i < numPoints; ++i) {
- const int p_value = _data->constrainedPoints[i]+offset;
- CPPUNIT_ASSERT_EQUAL(parametersFiberDim,
- parametersSection->getFiberDimension(p_value));
- const PylithScalar* parametersVertex =
- parametersSection->restrictPoint(p_value);
- CPPUNIT_ASSERT(parametersVertex);
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->valuesInitial[i*numFixedDOF+iDOF],
- parametersVertex[initialIndex+iDOF],
- tolerance);
+ const PetscInt p_value = _data->constrainedPoints[i]+offset;
+ 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
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
// Check rate of change
- const int rateIndex = bc._parameters->sectionIndex("rate");
- const int rateFiberDim = bc._parameters->sectionFiberDim("rate");
- CPPUNIT_ASSERT_EQUAL(numFixedDOF, rateFiberDim);
+ PetscSection rateSection = bc._parameters->get("rate").petscSection();
+ Vec rateVec = bc._parameters->get("rate").localVector();
+ PetscScalar *rateArray;
+ err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
for (int i=0; i < numPoints; ++i) {
- const int p_value = _data->constrainedPoints[i]+offset;
- CPPUNIT_ASSERT_EQUAL(parametersFiberDim,
- parametersSection->getFiberDimension(p_value));
- const PylithScalar* parametersVertex =
- parametersSection->restrictPoint(p_value);
- CPPUNIT_ASSERT(parametersVertex);
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->valueRate,
- parametersVertex[rateIndex+iDOF],
- tolerance);
+ const PetscInt p_value = _data->constrainedPoints[i]+offset;
+ 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);
} // for
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
} // if
} // testInitialize
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -27,7 +27,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
@@ -159,73 +159,64 @@
CPPUNIT_ASSERT(0 != _data);
const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
+ DM subMesh = boundaryMesh.dmMesh();
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
- // Check boundary mesh
- CPPUNIT_ASSERT(!submesh.isNull());
+ assert(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(subMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int cellDim = boundaryMesh.dimension();
const int numCorners = _data->numCorners;
const int spaceDim = _data->spaceDim;
- const ALE::Obj<SieveSubMesh::label_sequence>& cells = submesh->heightStratum(1);
- const int numVertices = submesh->depthStratum(0)->size();
- const int numCells = cells->size();
- const int boundaryDepth = submesh->depth()-1; // depth of boundary cells
+ const int numVertices = vEnd-vStart;
+ const int numCells = cEnd-cStart;
CPPUNIT_ASSERT_EQUAL(_data->cellDim, cellDim);
CPPUNIT_ASSERT_EQUAL(_data->numVertices, numVertices);
CPPUNIT_ASSERT_EQUAL(_data->numCells, numCells);
- const ALE::Obj<SieveMesh::sieve_type>& sieve = submesh->getSieve();
- ALE::ISieveVisitor::PointRetriever<SieveSubMesh::sieve_type> pV(sieve->getMaxConeSize());
- int dp = 0;
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter) {
- const int numCorners = submesh->getNumCellCorners(*c_iter, boundaryDepth);
- CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
+ PetscInt dp = 0;
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ const PetscInt *cone;
+ PetscInt numCorners;
- sieve->cone(*c_iter, pV);
- const SieveSubMesh::point_type *cone = pV.getPoints();
- for(int p = 0; p < pV.getSize(); ++p, ++dp) {
+ err = DMComplexGetConeSize(subMesh, c, &numCorners);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCone(subMesh, c, &cone);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(_data->numCorners, numCorners);
+ for(PetscInt p = 0; p < numCorners; ++p, ++dp) {
CPPUNIT_ASSERT_EQUAL(_data->cells[dp], cone[p]);
}
- pV.clear();
} // for
// Check traction values
const int numQuadPts = _data->numQuadPts;
const int fiberDim = numQuadPts * spaceDim;
scalar_array tractionsCell(fiberDim);
- int index = 0;
+ PetscInt index = 0;
CPPUNIT_ASSERT(0 != bc._parameters);
- const ALE::Obj<SubRealUniformSection>& parametersSection =
- bc._parameters->section();
- CPPUNIT_ASSERT(!parametersSection.isNull());
- const int parametersFiberDim = bc._parameters->fiberDim();
- const int initialIndex = bc._parameters->sectionIndex("initial");
- const int initialFiberDim = bc._parameters->sectionFiberDim("initial");
+ PetscSection initialSection = bc._parameters->get("initial").petscSection();
+ Vec initialVec = bc._parameters->get("initial").localVector();
+ PetscScalar *initialArray;
+ CPPUNIT_ASSERT(initialSection);CPPUNIT_ASSERT(initialVec);
const PylithScalar tolerance = 1.0e-06;
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter) {
- CPPUNIT_ASSERT_EQUAL(parametersFiberDim,
- parametersSection->getFiberDimension(*c_iter));
- const PylithScalar* parametersCell = parametersSection->restrictPoint(*c_iter);
- CPPUNIT_ASSERT(parametersCell);
- CPPUNIT_ASSERT(initialIndex + initialFiberDim <= parametersFiberDim);
- const PylithScalar* tractionsCell = ¶metersCell[initialIndex];
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(initialSection, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(initialSection, c, &off);CHECK_PETSC_ERROR(err);
+ 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,
- tractionsCell[iQuad*spaceDim+iDim],
- tolerance);
- ++index;
+ const PylithScalar tractionsCellData = _data->tractionsCell[index];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionsCellData, initialArray[off+iQuad*spaceDim+iDim], tolerance);
+ ++index;
} // for
} // for
-
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
} // testInitialize
// ----------------------------------------------------------------------
@@ -244,31 +235,35 @@
const PylithScalar t = 0.0;
bc.integrateResidual(residual, t, &fields);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- CPPUNIT_ASSERT(!sieveMesh->depthStratum(0).isNull());
-
+ DM dmMesh = mesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
const PylithScalar* valsE = _data->valsResidual;
- const int totalNumVertices = sieveMesh->depthStratum(0)->size();
+
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ const int totalNumVertices = vEnd - vStart;
const int sizeE = _data->spaceDim * totalNumVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ 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;
- // std::cout << "computed residuals: " << std::endl;
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
for (int i=0; i < size; ++i)
- // std::cout << " " << vals[i] << std::endl;
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);
+ err = VecRestoreArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateResidual
// ----------------------------------------------------------------------
@@ -759,33 +754,36 @@
assert(0 != valuesE);
const topology::SubMesh& boundaryMesh = field.mesh();
- const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
- CPPUNIT_ASSERT(!submesh.isNull());
- const ALE::Obj<RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- submesh->heightStratum(1);
+ DM subMesh = boundaryMesh.dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
- const PylithScalar scale = field.scale();
+ CPPUNIT_ASSERT(subMesh);
+ err = DMComplexGetHeightStratum(subMesh, 1, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- const size_t ncells = _TestNeumann::ncells;
- CPPUNIT_ASSERT_EQUAL(ncells, cells->size());
+ const PylithScalar scale = field.scale();
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ PetscScalar *array;
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+ const PetscInt ncells = _TestNeumann::ncells;
+ CPPUNIT_ASSERT_EQUAL(ncells, cEnd-cStart);
+
// Check values associated with BC.
int icell = 0;
const PylithScalar tolerance = 1.0e-06;
- for(SieveSubMesh::label_sequence::iterator c_iter = cells->begin();
- c_iter != cells->end();
- ++c_iter, ++icell) {
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ PetscInt dof, off;
- const int fiberDim = section->getFiberDimension(*c_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
-
- const PylithScalar* values = section->restrictPoint(*c_iter);
+ err = PetscSectionGetDof(section, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(section, 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,
- values[iDim], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[icell*fiberDimE+iDim]/scale, array[off+iDim], tolerance);
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // _checkValues
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestPointForce.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -26,7 +26,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
@@ -92,10 +92,13 @@
_initialize(&mesh, &bc);
CPPUNIT_ASSERT(0 != _data);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
-
- const int numCells = sieveMesh->heightStratum(0)->size();
+ DM dmMesh = mesh.dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ const int numCells = cEnd-cStart;
const int numForceDOF = _data->numForceDOF;
const size_t numPoints = _data->numForcePts;
@@ -108,50 +111,44 @@
} // if
CPPUNIT_ASSERT(0 != bc._parameters);
- const ALE::Obj<RealUniformSection>& parametersSection =
- bc._parameters->section();
- CPPUNIT_ASSERT(!parametersSection.isNull());
- const int parametersFiberDim = bc._parameters->fiberDim();
+ PetscSection initialSection = bc._parameters->get("initial").petscSection();
+ Vec initialVec = bc._parameters->get("initial").localVector();
+ PetscScalar *initialArray;
+ CPPUNIT_ASSERT(initialSection);CPPUNIT_ASSERT(initialVec);
- const PylithScalar tolerance = 1.0e-06;
-
// Check values
- const int initialIndex = bc._parameters->sectionIndex("initial");
- const int initialFiberDim = bc._parameters->sectionFiberDim("initial");
- CPPUNIT_ASSERT_EQUAL(numForceDOF, initialFiberDim);
-
+ const PylithScalar tolerance = 1.0e-06;
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
for (int i=0; i < numPoints; ++i) {
const int p_force = _data->forcePoints[i]+offset;
+ PetscInt dof, off;
- CPPUNIT_ASSERT_EQUAL(parametersFiberDim,
- parametersSection->getFiberDimension(p_force));
- const PylithScalar* parametersVertex = parametersSection->restrictPoint(p_force);
- CPPUNIT_ASSERT(parametersVertex);
-
+ 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],
- parametersVertex[initialIndex+iDOF],
- tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->forceInitial[i*numForceDOF+iDOF], initialArray[off+iDOF], tolerance);
} // for
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
// Check rate of change
- const int rateIndex = bc._parameters->sectionIndex("rate");
- const int rateFiberDim = bc._parameters->sectionFiberDim("rate");
- CPPUNIT_ASSERT_EQUAL(numForceDOF, rateFiberDim);
+ PetscSection rateSection = bc._parameters->get("rate").petscSection();
+ Vec rateVec = bc._parameters->get("rate").localVector();
+ PetscScalar *rateArray;
+ CPPUNIT_ASSERT(rateSection);CPPUNIT_ASSERT(rateVec);
+ 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;
- CPPUNIT_ASSERT_EQUAL(parametersFiberDim,
- parametersSection->getFiberDimension(p_force));
- const PylithScalar* parametersVertex = parametersSection->restrictPoint(p_force);
- CPPUNIT_ASSERT(parametersVertex);
-
+ 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,
- parametersVertex[rateIndex+iDOF],
- tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->forceRate, rateArray[off+iDOF], tolerance);
} // for
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
} // testInitialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependentPoints.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependentPoints.cc 2012-10-23 03:13:10 UTC (rev 20885)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestTimeDependentPoints.cc 2012-10-23 14:46:46 UTC (rev 20886)
@@ -24,7 +24,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
@@ -109,7 +109,8 @@
static
void _checkValues(const PylithScalar* valuesE,
const int fiberDimE,
- const ALE::Obj<RealSection>& section,
+ PetscSection section,
+ Vec vec,
const PylithScalar scale);
} // _TestTimeDependentPoints
} // bc
@@ -228,39 +229,34 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check initial values.
- const ALE::Obj<RealSection>& initialSection =
- _bc->_parameters->get("initial").section();
- CPPUNIT_ASSERT(!initialSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::initial,
- numBCDOF, initialSection, forceScale);
+ PetscSection initialSection = _bc->_parameters->get("initial").petscSection();
+ Vec initialVec = _bc->_parameters->get("initial").localVector();
+ CPPUNIT_ASSERT(initialSection);CPPUNIT_ASSERT(initialVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::initial, numBCDOF, initialSection, initialVec, forceScale);
// Check rate values.
- const ALE::Obj<RealSection>& rateSection =
- _bc->_parameters->get("rate").section();
- CPPUNIT_ASSERT(!rateSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::rate,
- numBCDOF, rateSection, forceScale/timeScale);
+ PetscSection rateSection = _bc->_parameters->get("rate").petscSection();
+ Vec rateVec = _bc->_parameters->get("rate").localVector();
+ CPPUNIT_ASSERT(rateSection);CPPUNIT_ASSERT(rateVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::rate, numBCDOF, rateSection, rateVec, forceScale/timeScale);
// Check rate start time.
- const ALE::Obj<RealSection>& rateTimeSection =
- _bc->_parameters->get("rate time").section();
- CPPUNIT_ASSERT(!rateTimeSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::rateTime,
- 1, rateTimeSection, timeScale);
+ PetscSection rateTimeSection = _bc->_parameters->get("rate time").petscSection();
+ Vec rateTimeVec = _bc->_parameters->get("rate time").localVector();
+ CPPUNIT_ASSERT(rateTimeSection);CPPUNIT_ASSERT(rateTimeVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::rateTime, 1, rateTimeSection, rateTimeVec, timeScale);
// Check change values.
- const ALE::Obj<RealSection>& changeSection =
- _bc->_parameters->get("change").section();
- CPPUNIT_ASSERT(!changeSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::change,
- numBCDOF, changeSection, forceScale);
+ PetscSection changeSection = _bc->_parameters->get("change").petscSection();
+ Vec changeVec = _bc->_parameters->get("change").localVector();
+ CPPUNIT_ASSERT(changeSection);CPPUNIT_ASSERT(changeVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::change, numBCDOF, changeSection, changeVec, forceScale);
// Check change start time.
- const ALE::Obj<RealSection>& changeTimeSection =
- _bc->_parameters->get("change time").section();
- CPPUNIT_ASSERT(!changeTimeSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::changeTime,
- 1, changeTimeSection, timeScale);
+ PetscSection changeTimeSection = _bc->_parameters->get("change time").petscSection();
+ Vec changeTimeVec = _bc->_parameters->get("change time").localVector();
+ CPPUNIT_ASSERT(changeTimeSection);CPPUNIT_ASSERT(changeTimeVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::changeTime, 1, changeTimeSection, changeTimeVec, timeScale);
th.close();
} // testQueryDatabases
@@ -293,11 +289,10 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check values.
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::initial,
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::initial, numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueInitial
// ----------------------------------------------------------------------
@@ -329,11 +324,10 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check values.
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesRate,
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesRate, numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueRate
// ----------------------------------------------------------------------
@@ -365,11 +359,10 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check values.
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesChange,
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesChange, numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueChange
// ----------------------------------------------------------------------
@@ -404,11 +397,10 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check values.
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesChangeTH,
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesChangeTH, numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueChangeTH
// ----------------------------------------------------------------------
@@ -466,11 +458,10 @@
_TestTimeDependentPoints::valuesRate[i] +
_TestTimeDependentPoints::valuesChangeTH[i];
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(&valuesE[0],
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(&valuesE[0], numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueAll
// ----------------------------------------------------------------------
@@ -504,11 +495,10 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check values.
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesIncrInitial,
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesIncrInitial, numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueIncrInitial
// ----------------------------------------------------------------------
@@ -542,11 +532,10 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check values.
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesIncrRate,
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesIncrRate, numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueIncrRate
// ----------------------------------------------------------------------
@@ -580,11 +569,10 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check values.
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesIncrChange,
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesIncrChange, numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueIncrChange
// ----------------------------------------------------------------------
@@ -621,11 +609,10 @@
CPPUNIT_ASSERT(_bc->_parameters);
// Check values.
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesIncrChangeTH,
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(_TestTimeDependentPoints::valuesIncrChangeTH, numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueIncrChangeTH
// ----------------------------------------------------------------------
@@ -685,11 +672,10 @@
_TestTimeDependentPoints::valuesIncrRate[i] +
_TestTimeDependentPoints::valuesIncrChangeTH[i];
- const ALE::Obj<RealSection>& valueSection =
- _bc->_parameters->get("value").section();
- CPPUNIT_ASSERT(!valueSection.isNull());
- _TestTimeDependentPoints::_checkValues(&valuesE[0],
- numBCDOF, valueSection, forceScale);
+ PetscSection valueSection = _bc->_parameters->get("value").petscSection();
+ Vec valueVec = _bc->_parameters->get("value").localVector();
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ _TestTimeDependentPoints::_checkValues(&valuesE[0], numBCDOF, valueSection, valueVec, forceScale);
} // testCalculateValueIncrAll
// ----------------------------------------------------------------------
@@ -697,26 +683,30 @@
void
pylith::bc::_TestTimeDependentPoints::_checkValues(const PylithScalar* valuesE,
const int fiberDimE,
- const ALE::Obj<RealSection>& section,
+ PetscSection section,
+ Vec vec,
const PylithScalar scale)
{ // _checkValues
- CPPUNIT_ASSERT(!section.isNull());
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+ PetscScalar *array;
+ PetscErrorCode err;
const PylithScalar tolerance = 1.0e-06;
// Check values at points associated with BC.
const int npointsIn = _TestTimeDependentPoints::npointsIn;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
for (int i=0; i < npointsIn; ++i) {
+ PetscInt dof, off;
const int p_bc = _TestTimeDependentPoints::pointsIn[i];
- const int fiberDim = section->getFiberDimension(p_bc);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* values = section->restrictPoint(p_bc);
- CPPUNIT_ASSERT(values);
+ err = PetscSectionGetDof(section, p_bc, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(section, p_bc, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
for (int iDim=0; iDim < fiberDimE; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i*fiberDimE+iDim]/scale,
- values[iDim], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i*fiberDimE+iDim]/scale, array[off+iDim], tolerance);
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // _checkValues
More information about the CIG-COMMITS
mailing list