[cig-commits] r20694 - in short/3D/PyLith/trunk: libsrc/pylith/feassemble libsrc/pylith/materials libsrc/pylith/topology unittests/libtests/feassemble
knepley at geodynamics.org
knepley at geodynamics.org
Thu Sep 6 08:12:20 PDT 2012
Author: knepley
Date: 2012-09-06 08:12:20 -0700 (Thu, 06 Sep 2012)
New Revision: 20694
Modified:
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitLgDeform.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicitLgDeform.cc
Log:
Disable field.section(), corrected assembly routines so that all feassemble tests pass
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -191,43 +191,43 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- scalar_array accCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
- RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+ topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
+ PetscSection accSection = acc.petscSection();
+ Vec accVec = acc.localVector();
+ assert(accSection);assert(accVec);
+
+ topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
+ PetscSection velSection = vel.petscSection();
+ Vec velVec = vel.localVector();
+ assert(velSection);assert(velVec);
- 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]);
-
- scalar_array dispCell(numBasis*spaceDim);
scalar_array dispAdjCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection = fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
-
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -246,19 +246,21 @@
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -266,7 +268,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -277,15 +279,14 @@
_resetCellVector();
// Restrict input fields to cell
- accVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, accVisitor);
+ const PetscScalar *accArray, *velArray, *dispTArray;
+ PetscInt accSize, velSize, dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == accSize);
+ assert(dispTSize == accSize);
- velVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, velVisitor);
-
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -306,8 +307,7 @@
const scalar_array& density = _material->calcDensity();
quadPtsGlobal = quadPtsNondim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
+ _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(), lengthScale);
// Compute action for element body forces
spatialdata::spatialdb::SpatialDB* db = _gravityField;
@@ -316,8 +316,7 @@
&quadPtsGlobal[0], spaceDim, cs);
if (err)
throw std::runtime_error("Unable to get gravity vector for point.");
- _normalizer->nondimensionalize(&gravVec[0], gravVec.size(),
- gravityScale);
+ _normalizer->nondimensionalize(&gravVec[0], gravVec.size(), gravityScale);
const PylithScalar wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
for (int iBasis = 0, iQ = iQuad * numBasis; iBasis < numBasis; ++iBasis) {
const PylithScalar valI = wt * basis[iQ + iBasis];
@@ -339,7 +338,7 @@
const PylithScalar valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -= valIJ *
- accCell[jBasis*spaceDim+iDim];
+ accArray[jBasis*spaceDim+iDim];
} // for
} // for
} // for
@@ -351,11 +350,13 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- dispAdjCell = dispCell + viscosity * velCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispAdjCell,
- numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispAdjCell, numBasis, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, false);
@@ -372,15 +373,16 @@
#endif
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(2*spaceDim))));
+ PetscLogFlops(numCells*numQuadPts*(3+numBasis*(1+numBasis*(2*spaceDim))));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -462,44 +464,43 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- scalar_array accCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
- RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+ topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
+ PetscSection accSection = acc.petscSection();
+ Vec accVec = acc.localVector();
+ assert(accSection);assert(accVec);
+
+ topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
+ PetscSection velSection = vel.petscSection();
+ Vec velVec = vel.localVector();
+ assert(velSection);assert(velVec);
- 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]);
-
- scalar_array dispCell(numBasis*spaceDim);
scalar_array dispAdjCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection =
- fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
-
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -521,19 +522,21 @@
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -541,7 +544,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -552,15 +555,14 @@
_resetCellVector();
// Restrict input fields to cell
- accVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, accVisitor);
+ const PetscScalar *accArray, *velArray, *dispTArray;
+ PetscInt accSize, velSize, dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == accSize);
+ assert(dispTSize == accSize);
- velVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, velVisitor);
-
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -620,7 +622,7 @@
for (int iBasis = 0; iBasis < numBasis; ++iBasis)
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -= valuesIJ[iBasis] *
- accCell[iBasis*spaceDim+iDim];
+ accArray[iBasis*spaceDim+iDim];
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(numQuadPts*(4+numBasis*3));
@@ -630,11 +632,14 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- dispAdjCell = dispCell + viscosity * velCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispAdjCell,
- numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispAdjCell, numBasis, numQuadPts);
+
const scalar_array& stressCell = _material->calcStress(strainCell, false);
#if defined(DETAILED_EVENT_LOGGING)
@@ -650,16 +655,17 @@
#endif
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(4+numBasis*3));
+ PetscLogFlops(numCells*numQuadPts*(4+numBasis*3));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidualLumped
@@ -699,19 +705,28 @@
"different dimensions than the spatial dimension.");
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- const ALE::Obj<RealSection>& solnSection = fields->solution().section();
- assert(!solnSection.isNull());
+ PetscSection solnSection = fields->solution().petscSection();
+ Vec solnVec = fields->solution().localVector();
+ scalar_array coordinatesCell(numBasis*spaceDim);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
+
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
assert(0 != jacobianMat);
@@ -721,44 +736,27 @@
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solnSection);
- assert(!globalOrder.isNull());
- // We would need to request unique points here if we had an interpolated mesh
- const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
- assert(!sieve.isNull());
- const int closureSize =
- int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
- assert(closureSize >= 0);
- IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
- closureSize*spaceDim);
-
- scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
-
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -767,7 +765,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -807,19 +805,18 @@
#endif
// Assemble cell contribution into PETSc matrix.
- jacobianVisitor.clear();
- PetscErrorCode err = updateOperator(jacobianMat, *sieveMesh->getSieve(),
- jacobianVisitor, *c_iter,
- &_cellMatrix[0], ADD_VALUES);
+ err = DMComplexMatSetClosure(dmMesh, solnSection, PETSC_NULL, jacobianMat, cell, &_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 = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
+ PetscLogFlops(numCells*numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
_logger->eventEnd(computeEvent);
#endif
@@ -862,15 +859,17 @@
"different dimensions than the spatial dimension.");
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get parameters used in integration.
const PylithScalar dt = _dt;
const PylithScalar dt2 = dt*dt;
@@ -878,35 +877,36 @@
scalar_array valuesIJ(numBasis);
// Get sections
- const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
- assert(!jacobianSection.isNull());
- UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
+ PetscSection jacSection = jacobian->petscSection();
+ Vec jacVec = jacobian->localVector();
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
#if defined(DETAILED_EVENT_LOGGING)
@@ -915,7 +915,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -957,16 +957,17 @@
#endif
// Assemble cell contribution into lumped matrix.
- jacobianVisitor.clear();
- sieveMesh->updateClosure(*c_iter, jacobianVisitor);
+ err = DMComplexVecSetClosure(dmMesh, jacSection, jacVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(numQuadPts*(4 + numBasis*3) + numBasis*spaceDim));
+ PetscLogFlops(numCells*(numQuadPts*(4 + numBasis*3) + numBasis*spaceDim));
_logger->eventEnd(computeEvent);
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -181,21 +181,22 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- scalar_array accCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
- RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+ topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
+ PetscSection accSection = acc.petscSection();
+ Vec accVec = acc.localVector();
+ assert(accSection);assert(accVec);
#if 0 // Numerical damping not yet implemented
scalar_array velCell(numBasis*spaceDim);
@@ -203,24 +204,25 @@
fields->get("velocity(t)").section();
assert(!velSection.isNull());
RestrictVisitor velVisitor(*velSection, velCell.size(), &velCell[0]);
-
scalar_array dispAdjCell(numBasis*spaceDim);
+#else
+ scalar_array dispTCell(numBasis*spaceDim);
#endif
- scalar_array dispCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection = fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -237,36 +239,37 @@
_logger->eventBegin(computeEvent);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Reset element vector to zero
_resetCellVector();
// Restrict input fields to cell
- accVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, accVisitor);
-
+ const PetscScalar *accArray, *velArray, *dispTArray;
+ PetscInt accSize, velSize, dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
#if 0 // Numerical damping not yet implemented.
- velVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, velVisitor);
+ err = DMComplexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == accSize);
#endif
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ assert(dispTSize == accSize);
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
-
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
const scalar_array& basisDeriv = _quadrature->basisDeriv();
@@ -317,7 +320,7 @@
const PylithScalar valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -=
- valIJ * accCell[jBasis*spaceDim+iDim];
+ valIJ * accArray[jBasis*spaceDim+iDim];
} // for
} // for
} // for
@@ -328,21 +331,26 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- dispAdjCell = dispCell + viscosity * velCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+#else
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispTCell[i] = dispTArray[i];}
#endif
+ err = DMComplexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Compute B(transpose) * sigma, first computing strains
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCell,
- numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTCell, numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, true);
- CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispCell);
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispTCell);
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_logger->eventEnd(computeEvent);
} // integrateResidual
@@ -422,46 +430,47 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- scalar_array accCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
- RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+ topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
+ PetscSection accSection = acc.petscSection();
+ Vec accVec = acc.localVector();
+ assert(accSection);assert(accVec);
#if 0 // Numerical damping not yet implemented
- 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]);
-
scalar_array dispAdjCell(numBasis*spaceDim);
+ topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
+ PetscSection velSection = vel.petscSection();
+ Vec velVec = vel.localVector();
+ assert(velSection);assert(velVec);
+#else
+ scalar_array dispTCell(numBasis*spaceDim);
#endif
+
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
- scalar_array dispCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection = fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
-
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -476,36 +485,37 @@
_logger->eventBegin(computeEvent);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Reset element vector to zero
_resetCellVector();
// Restrict input fields to cell
- accVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, accVisitor);
-
+ const PetscScalar *accArray, *velArray, *dispTArray;
+ PetscInt accSize, velSize, dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
#if 0 // Numerical damping not yet implemented.
- velVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, velVisitor);
+ err = DMComplexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == accSize);
#endif
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ assert(dispTSize == accSize);
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
-
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
const scalar_array& basisDeriv = _quadrature->basisDeriv();
@@ -561,7 +571,7 @@
for (int iBasis = 0; iBasis < numBasis; ++iBasis)
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -= valuesIJ[iBasis] *
- accCell[iBasis*spaceDim+iDim];
+ accArray[iBasis*spaceDim+iDim];
PetscLogFlops(numQuadPts*(4+numBasis*3));
#if 0 // Numerical damping not yet implemented. Is small strain
@@ -569,20 +579,26 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- dispAdjCell = dispCell + viscosity * velCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+#else
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispTCell[i] = dispTArray[i];}
#endif
+ err = DMComplexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+
// Compute B(transpose) * sigma, first computing strains
- _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCell,
- numBasis, numQuadPts, spaceDim);
+ _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTCell, numBasis, numQuadPts, spaceDim);
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
const scalar_array& stressCell = _material->calcStress(strainCell, true);
- CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispCell);
+ CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispTCell);
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_logger->eventEnd(computeEvent);
} // integrateResidualLumped
@@ -623,19 +639,28 @@
"different dimensions than the spatial dimension.");
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- const ALE::Obj<RealSection>& solnSection = fields->solution().section();
- assert(!solnSection.isNull());
+ PetscSection solnSection = fields->solution().petscSection();
+ Vec solnVec = fields->solution().localVector();
+ scalar_array coordinatesCell(numBasis*spaceDim);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
+
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
assert(0 != jacobianMat);
@@ -645,43 +670,26 @@
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solnSection);
- assert(!globalOrder.isNull());
- // We would need to request unique points here if we had an interpolated mesh
- const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
- assert(!sieve.isNull());
- const int closureSize =
- int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
- assert(closureSize >= 0);
- IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
- closureSize*spaceDim);
-
- scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
-
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Reset element matrix to zero
_resetCellMatrix();
@@ -712,12 +720,11 @@
PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
// Assemble cell contribution into PETSc matrix.
- jacobianVisitor.clear();
- PetscErrorCode err = updateOperator(jacobianMat, *sieveMesh->getSieve(),
- jacobianVisitor, *c_iter,
- &_cellMatrix[0], ADD_VALUES);
+ err = DMComplexMatSetClosure(dmMesh, solnSection, PETSC_NULL, jacobianMat, cell, &_cellMatrix[0], ADD_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_needNewJacobian = false;
_material->resetNeedNewJacobian();
@@ -761,50 +768,53 @@
"different dimensions than the spatial dimension.");
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get parameters used in integration.
const PylithScalar dt = _dt;
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
// Get sections
- const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
- assert(!jacobianSection.isNull());
- UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
+ PetscSection jacSection = jacobian->petscSection();
+ Vec jacVec = jacobian->localVector();
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Reset element matrix to zero
_resetCellMatrix();
@@ -836,9 +846,10 @@
_lumpCellMatrix();
// Assemble cell contribution into lumped matrix.
- jacobianVisitor.clear();
- sieveMesh->updateClosure(*c_iter, jacobianVisitor);
+ err = DMComplexVecSetClosure(dmMesh, jacSection, jacVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_needNewJacobian = false;
_material->resetNeedNewJacobian();
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTet4.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -169,44 +169,43 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- scalar_array accCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
- RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+ topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
+ PetscSection accSection = acc.petscSection();
+ Vec accVec = acc.localVector();
+ assert(accSection);assert(accVec);
- 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]);
+ topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
+ PetscSection velSection = vel.petscSection();
+ Vec velVec = vel.localVector();
+ assert(velSection);assert(velVec);
- scalar_array dispCell(numBasis*spaceDim);
scalar_array dispAdjCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection =
- fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -225,18 +224,20 @@
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
// Compute geometry information for current cell
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
const PylithScalar volume = _volume(coordinatesCell);
assert(volume > 0.0);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -244,7 +245,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -255,15 +256,14 @@
_resetCellVector();
// Restrict input fields to cell
- accVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, accVisitor);
+ const PetscScalar *accArray, *velArray, *dispTArray;
+ PetscInt accSize, velSize, dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == accSize);
+ assert(dispTSize == accSize);
- velVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, velVisitor);
-
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -306,7 +306,7 @@
for (int jBasis = 0; jBasis < numBasis; ++jBasis)
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -=
- wtVertex * accCell[jBasis*spaceDim+iDim];
+ wtVertex * accArray[jBasis*spaceDim+iDim];
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(3 + numBasis*numBasis*spaceDim*2);
@@ -316,7 +316,10 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- dispAdjCell = dispCell + viscosity * velCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Compute B(transpose) * sigma, first computing strains
const PylithScalar x0 = coordinatesCell[0];
@@ -408,15 +411,14 @@
#endif
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(3 + numBasis*numBasis*spaceDim*2 + 196+84));
+ PetscLogFlops(numCells*(3 + numBasis*numBasis*spaceDim*2 + 196+84));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -478,47 +480,43 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
- const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
- assert(!sieve.isNull());
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
// Get sections
- scalar_array accCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
- RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+ topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
+ PetscSection accSection = acc.petscSection();
+ Vec accVec = acc.localVector();
+ assert(accSection);assert(accVec);
- 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]);
-
- scalar_array dispCell(numBasis*spaceDim);
+ topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
+ PetscSection velSection = vel.petscSection();
+ Vec velVec = vel.localVector();
+ assert(velSection);assert(velVec);
+
scalar_array dispAdjCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection =
- fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
+
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
-
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -540,47 +538,34 @@
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(restrictEvent);
#endif
// Restrict input fields to cell
-#if 1
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const PetscScalar *accArray, *velArray, *dispTArray;
+ PetscInt accSize, velSize, dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == accSize);
+ assert(dispTSize == accSize);
- accVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, accVisitor);
-
- velVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, velVisitor);
-
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
-#else
- coordsVisitor.clear();
- sieve->orientedConeOpt(*c_iter, coordsVisitor, numBasis, spaceDim);
-
- accVisitor.clear();
- sieve->orientedConeOpt(*c_iter, accVisitor, numBasis, spaceDim);
-
- velVisitor.clear();
- sieve->orientedConeOpt(*c_iter, velVisitor, numBasis, spaceDim);
-
- dispVisitor.clear();
- sieve->orientedConeOpt(*c_iter, dispVisitor, numBasis, spaceDim);
-#endif
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(geometryEvent);
#endif
// Compute geometry information for current cell
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
const PylithScalar volume = _volume(coordinatesCell);
+ assert(volume > 0.0);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -588,7 +573,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Get density at quadrature points for this cell
const scalar_array& density = _material->calcDensity();
@@ -631,7 +616,7 @@
// Compute action for inertial terms
const PylithScalar wtVertex = density[0] * volume / 4.0;
- _cellVector -= wtVertex * accCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {_cellVector[i] -= wtVertex * accArray[i];}
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(2 + numBasis*spaceDim*2);
@@ -641,7 +626,10 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- dispAdjCell = dispCell + viscosity * velCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Compute B(transpose) * sigma, first computing strains
const PylithScalar x0 = coordinatesCell[0];
@@ -745,16 +733,17 @@
#endif
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(2 + numBasis*spaceDim*2 + 196+84));
+ PetscLogFlops(numCells*(2 + numBasis*spaceDim*2 + 196+84));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidualLumped
@@ -795,20 +784,22 @@
"different dimensions than the spatial dimension.");
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- scalar_array dispCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection =
- fields->get("disp(t)").section();
- assert(!dispSection.isNull());
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
@@ -819,20 +810,12 @@
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispSection);
- assert(!globalOrder.isNull());
- // We would need to request unique points here if we had an interpolated mesh
- IndicesVisitor jacobianVisitor(*dispSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
-
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -840,18 +823,20 @@
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
// Compute geometry information for current cell
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
const PylithScalar volume = _volume(coordinatesCell);
assert(volume > 0.0);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -859,7 +844,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -890,19 +875,18 @@
#endif
// Assemble cell contribution into PETSc matrix.
- jacobianVisitor.clear();
- PetscErrorCode err = updateOperator(jacobianMat, *sieveMesh->getSieve(),
- jacobianVisitor, *c_iter,
- &_cellMatrix[0], ADD_VALUES);
+ err = DMComplexMatSetClosure(dmMesh, dispTSection, PETSC_NULL, jacobianMat, cell, &_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 = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(3+numBasis*numBasis*spaceDim*1));
+ PetscLogFlops(numCells*(3+numBasis*numBasis*spaceDim*1));
_logger->eventEnd(computeEvent);
#endif
@@ -946,47 +930,51 @@
"different dimensions than the spatial dimension.");
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get parameters used in integration.
const PylithScalar dt = _dt;
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
// Get sections
- const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
- assert(!jacobianSection.isNull());
- UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
+ PetscSection jacSection = jacobian->petscSection();
+ Vec jacVec = jacobian->localVector();
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
const PylithScalar volume = _volume(coordinatesCell);
+ assert(volume > 0.0);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -994,7 +982,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -1013,16 +1001,17 @@
#endif
// Assemble cell contribution into lumped matrix.
- jacobianVisitor.clear();
- sieveMesh->updateClosure(*c_iter, jacobianVisitor);
+ err = DMComplexVecSetClosure(dmMesh, jacSection, jacVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*3);
+ PetscLogFlops(numCells*3);
_logger->eventEnd(computeEvent);
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -169,44 +169,43 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- scalar_array accCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
- RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+ topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
+ PetscSection accSection = acc.petscSection();
+ Vec accVec = acc.localVector();
+ assert(accSection);assert(accVec);
- 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]);
+ topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
+ PetscSection velSection = vel.petscSection();
+ Vec velVec = vel.localVector();
+ assert(velSection);assert(velVec);
- scalar_array dispCell(numBasis*spaceDim);
scalar_array dispAdjCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection =
- fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -225,18 +224,20 @@
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
// Compute geometry information for current cell
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
const PylithScalar area = _area(coordinatesCell);
assert(area > 0.0);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -244,7 +245,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -255,15 +256,14 @@
_resetCellVector();
// Restrict input fields to cell
- accVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, accVisitor);
+ const PetscScalar *accArray, *velArray, *dispTArray;
+ PetscInt accSize, velSize, dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == accSize);
+ assert(dispTSize == accSize);
- velVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, velVisitor);
-
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -306,7 +306,7 @@
for (int jBasis = 0; jBasis < numBasis; ++jBasis)
for (int iDim = 0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] -=
- wtVertex * accCell[jBasis*spaceDim+iDim];
+ wtVertex * accArray[jBasis*spaceDim+iDim];
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(3 + numBasis*numBasis*spaceDim*2);
@@ -316,7 +316,10 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- dispAdjCell = dispCell + viscosity * velCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Compute B(transpose) * sigma, first computing strains
const PylithScalar x0 = coordinatesCell[0];
@@ -371,15 +374,16 @@
#endif
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(3 + numBasis*numBasis*spaceDim*2 + 34+30));
+ PetscLogFlops(numCells*(3 + numBasis*numBasis*spaceDim*2 + 34+30));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidual
@@ -441,47 +445,43 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
- const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
- assert(!sieve.isNull());
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
// Get sections
- scalar_array accCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
- RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
+ topology::Field<topology::Mesh>& acc = fields->get("acceleration(t)");
+ PetscSection accSection = acc.petscSection();
+ Vec accVec = acc.localVector();
+ assert(accSection);assert(accVec);
- 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]);
-
- scalar_array dispCell(numBasis*spaceDim);
+ topology::Field<topology::Mesh>& vel = fields->get("velocity(t)");
+ PetscSection velSection = vel.petscSection();
+ Vec velVec = vel.localVector();
+ assert(velSection);assert(velVec);
+
scalar_array dispAdjCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection =
- fields->get("disp(t)").section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
+
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
-
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -503,45 +503,34 @@
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(restrictEvent);
#endif
// Restrict input fields to cell
-#if 1
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const PetscScalar *accArray, *velArray, *dispTArray;
+ PetscInt accSize, velSize, dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ assert(velSize == accSize);
+ assert(dispTSize == accSize);
- accVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, accVisitor);
-
- velVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, velVisitor);
-
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
-#else
- coordsVisitor.clear();
- sieve->orientedConeOpt(*c_iter, coordsVisitor, numBasis, spaceDim);
-
- accVisitor.clear();
- sieve->orientedConeOpt(*c_iter, accVisitor, numBasis, spaceDim);
-
- dispVisitor.clear();
- sieve->orientedConeOpt(*c_iter, dispVisitor, numBasis, spaceDim);
-#endif
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(geometryEvent);
#endif
// Compute geometry information for current cell
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
const PylithScalar area = _area(coordinatesCell);
assert(area > 0.0);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -549,7 +538,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Get density at quadrature points for this cell
const scalar_array& density = _material->calcDensity();
@@ -592,7 +581,7 @@
// Compute action for inertial terms
const PylithScalar wtVertex = density[0] * area / 3.0;
- _cellVector -= wtVertex * accCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {_cellVector[i] -= wtVertex * accArray[i];}
#if defined(DETAILED_EVENT_LOGGING)
PetscLogFlops(2 + numBasis*spaceDim*2);
@@ -602,7 +591,10 @@
// Numerical damping. Compute displacements adjusted by velocity
// times normalized viscosity.
- dispAdjCell = dispCell + viscosity * velCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispAdjCell[i] = dispTArray[i] + viscosity * velArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, accSection, accVec, cell, &accSize, &accArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, velSection, velVec, cell, &velSize, &velArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Compute B(transpose) * sigma, first computing strains
const PylithScalar x0 = coordinatesCell[0];
@@ -656,16 +648,17 @@
#endif
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(2 + numBasis*spaceDim*2 + 34+30));
+ PetscLogFlops(numCells*(2 + numBasis*spaceDim*2 + 34+30));
_logger->eventEnd(computeEvent);
#endif
} // integrateResidualLumped
@@ -706,20 +699,22 @@
"different dimensions than the spatial dimension.");
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- scalar_array dispCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispSection =
- fields->get("disp(t)").section();
- assert(!dispSection.isNull());
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
@@ -730,20 +725,12 @@
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispSection);
- assert(!globalOrder.isNull());
- // We would need to request unique points here if we had an interpolated mesh
- IndicesVisitor jacobianVisitor(*dispSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
-
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -751,18 +738,20 @@
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
// Compute geometry information for current cell
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
const PylithScalar area = _area(coordinatesCell);
assert(area > 0.0);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -770,7 +759,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -801,19 +790,18 @@
#endif
// Assemble cell contribution into PETSc matrix.
- jacobianVisitor.clear();
- PetscErrorCode err = updateOperator(jacobianMat, *sieveMesh->getSieve(),
- jacobianVisitor, *c_iter,
- &_cellMatrix[0], ADD_VALUES);
+ err = DMComplexMatSetClosure(dmMesh, dispTSection, PETSC_NULL, jacobianMat, cell, &_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 = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*(3+numBasis*numBasis*spaceDim*1));
+ PetscLogFlops(numCells*(3+numBasis*numBasis*spaceDim*1));
_logger->eventEnd(computeEvent);
#endif
@@ -857,48 +845,51 @@
"different dimensions than the spatial dimension.");
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get parameters used in integration.
const PylithScalar dt = _dt;
const PylithScalar dt2 = dt*dt;
assert(dt > 0);
// Get sections
- const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
- assert(!jacobianSection.isNull());
- UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
+ PetscSection jacSection = jacobian->petscSection();
+ Vec jacVec = jacobian->localVector();
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
#endif
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(geometryEvent);
#endif
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
const PylithScalar area = _area(coordinatesCell);
assert(area > 0.0);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(geometryEvent);
@@ -906,7 +897,7 @@
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(stateVarsEvent);
@@ -924,16 +915,17 @@
#endif
// Assemble cell contribution into lumped matrix.
- jacobianVisitor.clear();
- sieveMesh->updateClosure(*c_iter, jacobianVisitor);
+ err = DMComplexVecSetClosure(dmMesh, jacSection, jacVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
- PetscLogFlops(cells->size()*3);
+ PetscLogFlops(numCells*3);
_logger->eventEnd(computeEvent);
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -171,36 +171,37 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _material->id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get sections
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim, &dispTCell[0]);
- const ALE::Obj<RealSection>& dispTIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- assert(!dispTIncrSection.isNull());
- RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
- numBasis*spaceDim, &dispTIncrCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
+ topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
+ PetscSection dispTIncrSection = dispTIncr.petscSection();
+ Vec dispTIncrVec = dispTIncr.localVector();
+ assert(dispTIncrSection);assert(dispTIncrVec);
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -212,29 +213,32 @@
_logger->eventBegin(computeEvent);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Reset element vector to zero
_resetCellVector();
// Restrict input fields to cell
- dispTVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTVisitor);
- dispTIncrVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+ const PetscScalar *dispTArray, *dispTIncrArray;
+ PetscInt dispTSize, dispTIncrSize;
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ assert(dispTSize == dispTIncrSize);
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -244,7 +248,9 @@
// Compute current estimate of displacement at time t+dt using
// solution increment.
- dispTpdtCell = dispTCell + dispTIncrCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispTpdtCell[i] = dispTArray[i] + dispTIncrArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
// Compute body force vector if gravity is being used.
if (0 != _gravityField) {
@@ -295,9 +301,10 @@
}
#endif
// Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_logger->eventEnd(computeEvent);
} // integrateResidual
@@ -372,27 +379,29 @@
scalar_array stressCell(numQuadPts*tensorSize);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
+ DM dmMesh = fields->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+
+ assert(dmMesh);
const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ err = DMComplexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
// Get sections
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim, &dispTCell[0]);
- const ALE::Obj<RealSection>& dispTIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- assert(!dispTIncrSection.isNull());
- RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
- numBasis*spaceDim, &dispTIncrCell[0]);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
+ topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
+ PetscSection dispTIncrSection = dispTIncr.petscSection();
+ Vec dispTIncrVec = dispTIncr.localVector();
+ assert(dispTIncrSection);assert(dispTIncrVec);
+
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
assert(0 != jacobianMat);
@@ -401,53 +410,43 @@
const PylithScalar dt = _dt;
assert(dt > 0);
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- dispTSection);
- assert(!globalOrder.isNull());
- // We would need to request unique points here if we had an interpolated mesh
- const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
- assert(!sieve.isNull());
- const int closureSize =
- int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
- assert(closureSize >= 0);
- IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
- closureSize*spaceDim);
-
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Reset element matrix to zero
_resetCellMatrix();
// Restrict input fields to cell
- dispTVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTVisitor);
- dispTIncrVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispTIncrVisitor);
+ const PetscScalar *dispTArray, *dispTIncrArray;
+ PetscInt dispTSize, dispTIncrSize;
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecGetClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ assert(dispTSize == dispTIncrSize);
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -456,7 +455,9 @@
// Compute current estimate of displacement at time t+dt using
// solution increment.
- dispTpdtCell = dispTCell + dispTIncrCell;
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispTpdtCell[i] = dispTArray[i] + dispTIncrArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecRestoreClosure(dmMesh, dispTIncrSection, dispTIncrVec, cell, &dispTIncrSize, &dispTIncrArray);CHECK_PETSC_ERROR(err);
// Compute deformation tensor, strains, and stresses
_calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispTpdtCell,
@@ -495,7 +496,7 @@
throw std::runtime_error("Lapack SVD failed");
minSV = svalues[n-7];
maxSV = svalues[0];
- std::cout << "Element " << *c_iter << std::endl;
+ std::cout << "Element " << cell << std::endl;
for(int i = 0; i < n; ++i)
std::cout << " sV["<<i<<"] = " << svalues[i] << std::endl;
std::cout << " kappa(elemMat) = " << maxSV/minSV << std::endl;
@@ -505,10 +506,7 @@
} // if
// Assemble cell contribution into PETSc matrix.
- jacobianVisitor.clear();
- PetscErrorCode err = updateOperator(jacobianMat, *sieveMesh->getSieve(),
- jacobianVisitor, *c_iter,
- &_cellMatrix[0], ADD_VALUES);
+ err = DMComplexMatSetClosure(dmMesh, dispTSection, PETSC_NULL, jacobianMat, cell, &_cellMatrix[0], ADD_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
} // for
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -107,34 +107,68 @@
assert(0 != _properties);
assert(0 != _stateVars);
- const ALE::Obj<RealSection>& propertiesSection = _properties->section();
- assert(!propertiesSection.isNull());
- propertiesSection->restrictPoint(cell, &_propertiesCell[0],
- _propertiesCell.size());
+ PetscSection propertiesSection = _properties->petscSection();
+ Vec propertiesVec = _properties->localVector();
+ PetscScalar *propertiesArray;
+ PetscInt dof, off;
+ PetscErrorCode err;
+ assert(propertiesSection);assert(propertiesVec);
+ err = PetscSectionGetDof(propertiesSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(propertiesSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(_propertiesCell.size() == dof);
+ err = VecGetArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {_propertiesCell[d] = propertiesArray[off+d];}
+ err = VecRestoreArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
+
if (hasStateVars()) {
- const ALE::Obj<RealSection>& stateVarsSection = _stateVars->section();
- assert(!stateVarsSection.isNull());
- stateVarsSection->restrictPoint(cell, &_stateVarsCell[0],
- _stateVarsCell.size());
+ PetscSection stateVarsSection = _stateVars->petscSection();
+ Vec stateVarsVec = _stateVars->localVector();
+ PetscScalar *stateVarsArray;
+ PetscInt dof, off;
+ PetscErrorCode err;
+
+ assert(stateVarsSection);assert(stateVarsVec);
+ err = PetscSectionGetDof(stateVarsSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(stateVarsSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(_stateVarsCell.size() == dof);
+ err = VecGetArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {_stateVarsCell[d] = stateVarsArray[off+d];}
+ err = VecRestoreArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);
} // if
_initialStressCell = 0.0;
_initialStrainCell = 0.0;
if (0 != _initialFields) {
if (_initialFields->hasField("initial stress")) {
- const ALE::Obj<RealSection>& stressSection =
- _initialFields->get("initial stress").section();
- assert(!stressSection.isNull());
- stressSection->restrictPoint(cell, &_initialStressCell[0],
- _initialStressCell.size());
+ PetscSection stressSection = _initialFields->get("initial stress").petscSection();
+ Vec stressVec = _initialFields->get("initial stress").localVector();
+ PetscScalar *stressArray;
+ PetscInt dof, off;
+ PetscErrorCode err;
+
+ assert(stressSection);assert(stressVec);
+ err = PetscSectionGetDof(stressSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(stressSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(_initialStressCell.size() == dof);
+ err = VecGetArray(stressVec, &stressArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {_initialStressCell[d] = stressArray[off+d];}
+ err = VecRestoreArray(stressVec, &stressArray);CHECK_PETSC_ERROR(err);
} // if
if (_initialFields->hasField("initial strain")) {
- const ALE::Obj<RealSection>& strainSection =
- _initialFields->get("initial strain").section();
- assert(!strainSection.isNull());
- strainSection->restrictPoint(cell, &_initialStrainCell[0],
- _initialStrainCell.size());
+ PetscSection strainSection = _initialFields->get("initial strain").petscSection();
+ Vec strainVec = _initialFields->get("initial strain").localVector();
+ PetscScalar *strainArray;
+ PetscInt dof, off;
+ PetscErrorCode err;
+
+ assert(strainSection);assert(strainVec);
+ err = PetscSectionGetDof(strainSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(strainSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(_initialStrainCell.size() == dof);
+ err = VecGetArray(strainVec, &strainArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {_initialStrainCell[d] = strainArray[off+d];}
+ err = VecRestoreArray(strainVec, &strainArray);CHECK_PETSC_ERROR(err);
} // if
} // if
} // retrievePropsAndVars
@@ -246,29 +280,33 @@
PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
// Get cells associated with material
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- retrievePropsAndVars(*c_iter);
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
+
+ retrievePropsAndVars(cell);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const PylithScalar dt =
- _stableTimeStepImplicit(&_propertiesCell[iQuad*numPropsQuadPt],
- numPropsQuadPt,
- &_stateVarsCell[iQuad*numVarsQuadPt],
- numVarsQuadPt);
+ _stableTimeStepImplicit(&_propertiesCell[iQuad*numPropsQuadPt],
+ numPropsQuadPt,
+ &_stateVarsCell[iQuad*numVarsQuadPt],
+ numVarsQuadPt);
if (dt < dtStable)
- dtStable = dt;
+ dtStable = dt;
} // for
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
return dtStable;
} // stableTimeStepImplicit
@@ -319,15 +357,26 @@
const int spaceDim = quadrature->spaceDim();
const int numBasis = quadrature->numBasis();
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
+ // Get cells associated with material
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
#endif
// Create arrays for querying
@@ -335,21 +384,15 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
scalar_array stressCell(numQuadPts*tensorSize);
- // Get cells associated with material
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
-
// Create field to hold initial stress state.
const int fiberDim = numQuadPts * tensorSize;
assert(fiberDim > 0);
- initialStress.newSection(cells, fiberDim);
+ int_array cellsTmp(cells, numCells);
+ initialStress.newSection(cellsTmp, fiberDim);
initialStress.allocate();
initialStress.zero();
- const ALE::Obj<RealSection>& initialStressSection = initialStress.section();
+ PetscSection initialStressSection = initialStress.petscSection();
+ Vec initialStressVec = initialStress.localVector();
// Setup databases for querying
_dbInitialStress->open();
@@ -391,16 +434,18 @@
const PylithScalar lengthScale = _normalizer->lengthScale();
const PylithScalar pressureScale = _normalizer->pressureScale();
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
quadrature->retrieveGeometry(*c_iter);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Dimensionalize coordinates for querying
@@ -430,7 +475,7 @@
_normalizer->nondimensionalize(&stressCell[0], stressCell.size(),
pressureScale);
- initialStressSection->updatePoint(*c_iter, &stressCell[0]);
+ err = DMComplexVecSetClosure(dmMesh, initialStressSection, initialStressVec, cell, &stressCell[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // for
// Close databases
@@ -465,15 +510,26 @@
const int spaceDim = quadrature->spaceDim();
const int numBasis = quadrature->numBasis();
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
+ // Get cells associated with material
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
#endif
// Create arrays for querying
@@ -481,21 +537,15 @@
scalar_array quadPtsGlobal(numQuadPts*spaceDim);
scalar_array strainCell(numQuadPts*tensorSize);
- // Get cells associated with material
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
-
// Create field to hold initial strain state.
const int fiberDim = numQuadPts * tensorSize;
assert(fiberDim > 0);
- initialStrain.newSection(cells, fiberDim);
+ int_array cellsTmp(cells, numCells);
+ initialStrain.newSection(cellsTmp, fiberDim);
initialStrain.allocate();
initialStrain.zero();
- const ALE::Obj<RealSection>& initialStrainSection = initialStrain.section();
+ PetscSection initialStrainSection = initialStrain.petscSection();
+ Vec initialStrainVec = initialStrain.localVector();
// Setup databases for querying
_dbInitialStrain->open();
@@ -537,16 +587,18 @@
const PylithScalar lengthScale = _normalizer->lengthScale();
const PylithScalar pressureScale = _normalizer->pressureScale();
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
quadrature->retrieveGeometry(*c_iter);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Dimensionalize coordinates for querying
@@ -572,7 +624,7 @@
} // if
} // for
- initialStrainSection->updatePoint(*c_iter, &strainCell[0]);
+ err = DMComplexVecSetClosure(dmMesh, initialStrainSection, initialStrainVec, cell, &strainCell[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // for
// Close databases
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -125,13 +125,16 @@
const int spaceDim = quadrature->spaceDim();
// Get cells associated with material
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", _id);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _id, &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
assert(0 != cs);
@@ -140,18 +143,21 @@
_properties->label("properties");
assert(0 != _properties);
int fiberDim = numQuadPts * _numPropsQuadPt;
- _properties->newSection(cells, fiberDim);
+ int_array cellsTmp(cells, numCells);
+
+ _properties->newSection(cellsTmp, fiberDim);
_properties->allocate();
_properties->zero();
- const ALE::Obj<RealSection>& propertiesSection = _properties->section();
- assert(!propertiesSection.isNull());
+ PetscSection propertiesSection = _properties->petscSection();
+ Vec propertiesVec = _properties->localVector();
#if !defined(PRECOMPUTE_GEOMETRY)
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateVec(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
#endif
// Create arrays for querying.
@@ -169,6 +175,8 @@
// Create field to hold state variables. We create the field even
// if there is no initial state, because this we will use this field
// to hold the state variables.
+ PetscSection stateVarsSection = PETSC_NULL;
+ Vec stateVarsVec = PETSC_NULL;
delete _stateVars; _stateVars = new topology::Field<topology::Mesh>(mesh);
_stateVars->label("state variables");
fiberDim = numQuadPts * _numVarsQuadPt;
@@ -178,16 +186,17 @@
_stateVars->newSection(*_properties, fiberDim);
_stateVars->allocate();
_stateVars->zero();
+ stateVarsSection = _stateVars->petscSection();
+ stateVarsVec = _stateVars->localVector();
+ assert(stateVarsSection);assert(stateVarsVec);
} // if
- const ALE::Obj<RealSection>& stateVarsSection =
- (fiberDim > 0) ? _stateVars->section() : 0;
// Create arrays for querying
const int numDBStateVars = _metadata.numDBStateVars();
scalar_array stateVarsQuery;
scalar_array stateVarsCell;
if (0 != _dbInitialState) {
- assert(!stateVarsSection.isNull());
+ assert(stateVarsSection);
assert(numDBStateVars > 0);
assert(_numVarsQuadPt > 0);
stateVarsQuery.resize(numDBStateVars);
@@ -201,17 +210,22 @@
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
-
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ PetscScalar *propertiesArray, *stateVarsArray;
+
+ err = VecGetArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
+ if (stateVarsVec) {err = VecGetArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);}
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- quadrature->retrieveGeometry(*c_iter);
+ quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ quadrature->computeGeometry(coordinatesCell, cell);
+ err = DMComplexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
const scalar_array& quadPtsNonDim = quadrature->quadPts();
@@ -261,10 +275,27 @@
} // for
// Insert cell contribution into fields
- propertiesSection->updatePoint(*c_iter, &propertiesCell[0]);
- if (0 != _dbInitialState)
- stateVarsSection->updatePoint(*c_iter, &stateVarsCell[0]);
+ PetscInt dof, off, d;
+
+ err = PetscSectionGetDof(propertiesSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(propertiesSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(dof == numQuadPts*_numPropsQuadPt);
+ for(PetscInt d = 0; d < dof; ++d) {
+ propertiesArray[off+d] = propertiesCell[d];
+ }
+ if (0 != _dbInitialState) {
+ err = PetscSectionGetDof(stateVarsSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(stateVarsSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(dof == numQuadPts*numDBStateVars);
+ for(PetscInt d = 0; d < dof; ++d) {
+ stateVarsArray[off+d] = stateVarsCell[d];
+ }
+ }
} // for
+ err = VecRestoreArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
+ if (stateVarsVec) {err = VecRestoreArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);}
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
// Close databases
_dbProperties->close();
@@ -336,14 +367,17 @@
} // else
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = field->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", _id);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = field->mesh().dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", _id, &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
topology::FieldBase::VectorFieldEnum fieldType = topology::FieldBase::OTHER;
if (propertyIndex >= 0) { // If field is a property
@@ -354,14 +388,15 @@
const int fiberDim = _metadata.getProperty(propertyIndex).fiberDim;
// Get properties section
- const ALE::Obj<RealSection>& propertiesSection = _properties->section();
- assert(!propertiesSection.isNull());
- const int totalPropsFiberDimLocal = (cells->size() > 0) ?
- propertiesSection->getFiberDimension(*cells->begin()) : 0;
- int totalPropsFiberDim = 0;
+ PetscSection propertiesSection = _properties->petscSection();
+ Vec propertiesVec = _properties->localVector();
+ assert(propertiesSection);
+ PetscInt totalPropsFiberDimLocal = 0;
+ PetscInt totalPropsFiberDim = 0;
+ if (numCells > 0) {err = PetscSectionGetDof(propertiesSection, cells[0], &totalPropsFiberDimLocal);CHECK_PETSC_ERROR(err);}
MPI_Allreduce((void *) &totalPropsFiberDimLocal,
- (void *) &totalPropsFiberDim, 1,
- MPI_INT, MPI_MAX, field->mesh().comm());
+ (void *) &totalPropsFiberDim, 1,
+ MPIU_INT, MPI_MAX, field->mesh().comm());
assert(totalPropsFiberDim > 0);
const int numPropsQuadPt = _numPropsQuadPt;
const int numQuadPts = totalPropsFiberDim / numPropsQuadPt;
@@ -369,27 +404,30 @@
const int totalFiberDim = numQuadPts * fiberDim;
// Allocate buffer for property field if necessary.
- const ALE::Obj<RealSection>& fieldSection = field->section();
- bool useCurrentField = !fieldSection.isNull();
- if (!fieldSection.isNull()) {
+ PetscSection fieldSection = field->petscSection();
+ Vec fieldVec = field->localVector();
+ bool useCurrentField = fieldSection != PETSC_NULL;
+ if (fieldSection) {
// check fiber dimension
- const int totalFiberDimCurrentLocal = (cells->size() > 0) ?
- fieldSection->getFiberDimension(*cells->begin()) : 0;
- int totalFiberDimCurrent = 0;
+ PetscInt totalFiberDimCurrentLocal = 0;
+ PetscInt totalFiberDimCurrent = 0;
+ if (numCells > 0) {err = PetscSectionGetDof(fieldSection, cells[0], &totalFiberDimCurrentLocal);CHECK_PETSC_ERROR(err);}
MPI_Allreduce((void *) &totalFiberDimCurrentLocal,
- (void *) &totalFiberDimCurrent, 1,
- MPI_INT, MPI_MAX, field->mesh().comm());
+ (void *) &totalFiberDimCurrent, 1,
+ MPIU_INT, MPI_MAX, field->mesh().comm());
assert(totalFiberDimCurrent > 0);
useCurrentField = totalFiberDim == totalFiberDimCurrent;
} // if
if (!useCurrentField) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("OutputFields");
- field->newSection(cells, totalFiberDim);
+ int_array cellsTmp(cells, numCells);
+
+ field->newSection(cellsTmp, totalFiberDim);
field->allocate();
logger.stagePop();
} // if
- assert(!fieldSection.isNull());
+ assert(fieldSection);
field->label(name);
field->scale(1.0);
fieldType = _metadata.getProperty(propertyIndex).fieldType;
@@ -399,21 +437,24 @@
scalar_array propertiesCell(numQuadPts*numPropsQuadPt);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- propertiesSection->restrictPoint(*c_iter,
- &propertiesCell[0], propertiesCell.size());
+ PetscScalar *fieldArray, *propertiesArray;
+
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
+ PetscInt off, poff;
+ err = PetscSectionGetOffset(fieldSection, cell, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(propertiesSection, cell, &poff);CHECK_PETSC_ERROR(err);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- _dimProperties(&propertiesCell[iQuad*numPropsQuadPt],
- numPropsQuadPt);
- for (int i=0; i < fiberDim; ++i)
- fieldCell[iQuad*fiberDim+i] =
- propertiesCell[iQuad*numPropsQuadPt+propOffset+i];
+ _dimProperties(&propertiesCell[iQuad*numPropsQuadPt], numPropsQuadPt);
+ for (int i=0; i < fiberDim; ++i)
+ fieldArray[iQuad*fiberDim + off+i] = propertiesArray[iQuad*numPropsQuadPt+propOffset + poff+i];
} // for
- fieldSection->updatePoint(*c_iter, &fieldCell[0]);
} // for
+ err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
} else { // field is a state variable
assert(stateVarIndex >= 0);
@@ -424,14 +465,15 @@
const int fiberDim = _metadata.getStateVar(stateVarIndex).fiberDim;
// Get state variables section
- const ALE::Obj<RealSection>& stateVarsSection = _stateVars->section();
- assert(!stateVarsSection.isNull());
- const int totalVarsFiberDimLocal = (cells->size() > 0) ?
- stateVarsSection->getFiberDimension(*cells->begin()) : 0;
- int totalVarsFiberDim = 0;
+ PetscSection stateVarsSection = _stateVars->petscSection();
+ Vec stateVarsVec = _stateVars->localVector();
+ assert(stateVarsSection);assert(stateVarsVec);
+ PetscInt totalVarsFiberDimLocal = 0;
+ PetscInt totalVarsFiberDim = 0;
+ if (numCells > 0) {err = PetscSectionGetDof(stateVarsSection, cells[0], &totalVarsFiberDimLocal);CHECK_PETSC_ERROR(err);}
MPI_Allreduce((void *) &totalVarsFiberDimLocal,
- (void *) &totalVarsFiberDim, 1,
- MPI_INT, MPI_MAX, field->mesh().comm());
+ (void *) &totalVarsFiberDim, 1,
+ MPIU_INT, MPI_MAX, field->mesh().comm());
assert(totalVarsFiberDim > 0);
const int numVarsQuadPt = _numVarsQuadPt;
const int numQuadPts = totalVarsFiberDim / numVarsQuadPt;
@@ -439,27 +481,30 @@
const int totalFiberDim = numQuadPts * fiberDim;
// Allocate buffer for state variable field if necessary.
- const ALE::Obj<RealSection>& fieldSection = field->section();
- bool useCurrentField = !fieldSection.isNull();
- if (!fieldSection.isNull()) {
+ PetscSection fieldSection = field->petscSection();
+ Vec fieldVec = field->localVector();
+ bool useCurrentField = fieldSection != PETSC_NULL;
+ if (fieldSection) {
// check fiber dimension
- const int totalFiberDimCurrentLocal = (cells->size() > 0) ?
- fieldSection->getFiberDimension(*cells->begin()) : 0;
- int totalFiberDimCurrent = 0;
+ PetscInt totalFiberDimCurrentLocal = 0;
+ PetscInt totalFiberDimCurrent = 0;
+ if (numCells > 0) {err = PetscSectionGetDof(fieldSection, cells[0], &totalFiberDimCurrentLocal);CHECK_PETSC_ERROR(err);}
MPI_Allreduce((void *) &totalFiberDimCurrentLocal,
- (void *) &totalFiberDimCurrent, 1,
- MPI_INT, MPI_MAX, field->mesh().comm());
+ (void *) &totalFiberDimCurrent, 1,
+ MPIU_INT, MPI_MAX, field->mesh().comm());
assert(totalFiberDimCurrent > 0);
useCurrentField = totalFiberDim == totalFiberDimCurrent;
} // if
if (!useCurrentField) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("OutputFields");
- field->newSection(cells, totalFiberDim);
+ int_array cellsTmp(cells, numCells);
+
+ field->newSection(cellsTmp, totalFiberDim);
field->allocate();
logger.stagePop();
} // if
- assert(!fieldSection.isNull());
+ assert(fieldSection);
fieldType = _metadata.getStateVar(stateVarIndex).fieldType;
field->label(name);
field->scale(1.0);
@@ -469,22 +514,27 @@
scalar_array stateVarsCell(numQuadPts*numVarsQuadPt);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- stateVarsSection->restrictPoint(*c_iter,
- &stateVarsCell[0], stateVarsCell.size());
+ PetscScalar *fieldArray, *stateVarsArray;
+
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
+ PetscInt off, soff;
+ err = PetscSectionGetOffset(fieldSection, cell, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(stateVarsSection, cell, &soff);CHECK_PETSC_ERROR(err);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- _dimStateVars(&stateVarsCell[iQuad*numVarsQuadPt],
- numVarsQuadPt);
- for (int i=0; i < fiberDim; ++i)
- fieldCell[iQuad*fiberDim+i] =
- stateVarsCell[iQuad*numVarsQuadPt+varOffset+i];
+ _dimStateVars(&stateVarsCell[iQuad*numVarsQuadPt], numVarsQuadPt);
+ for (int i=0; i < fiberDim; ++i)
+ fieldArray[iQuad*fiberDim + off+i] = stateVarsCell[iQuad*numVarsQuadPt+varOffset + soff+i];
} // for
- fieldSection->updatePoint(*c_iter, &fieldCell[0]);
} // for
+ err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);
} // if/else
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
topology::FieldBase::VectorFieldEnum multiType =
topology::FieldBase::MULTI_OTHER;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -365,7 +365,7 @@
throw std::runtime_error(msg.str());
} // if
- const ALE::Obj<section_type>& srcSection = src.section();
+ const ALE::Obj<section_type>& srcSection = src._section;
if (!srcSection.isNull()) {
_section->setChart(srcSection->getChart());
const chart_type& chart = _section->getChart();
@@ -393,7 +393,7 @@
// Clear memory
clear();
- const ALE::Obj<section_type>& srcSection = src.section();
+ const ALE::Obj<section_type>& srcSection = src._section;
if (!srcSection.isNull() && _section.isNull()) {
newSection();
} // if
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -27,6 +27,14 @@
inline
const ALE::Obj<section_type>&
pylith::topology::Field<mesh_type, section_type>::section(void) const {
+ std::ostringstream msg;
+
+ msg << "Sieve Sections are no longer supported: " << const_cast<Field*>(this)->_metadata["default"].label << "\n"
+ << " Destination section:\n"
+ << " space dim: " << spaceDim() << "\n"
+ << " vector field type: " << const_cast<Field*>(this)->_metadata["default"].vectorFieldType << "\n"
+ << " scale: " << const_cast<Field*>(this)->_metadata["default"].scale;
+ throw std::runtime_error(msg.str());
return _section;
}
@@ -64,6 +72,14 @@
const
mesh_type&
pylith::topology::Field<mesh_type, section_type>::mesh(void) const {
+ std::ostringstream msg;
+
+ msg << "Sieve Section Meshes are no longer supported: " << const_cast<Field*>(this)->_metadata["default"].label << "\n"
+ << " Destination section:\n"
+ << " space dim: " << spaceDim() << "\n"
+ << " vector field type: " << const_cast<Field*>(this)->_metadata["default"].vectorFieldType << "\n"
+ << " scale: " << const_cast<Field*>(this)->_metadata["default"].scale;
+ throw std::runtime_error(msg.str());
return _mesh;
}
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Jacobian.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -34,8 +34,7 @@
_matrix(0),
_valuesChanged(true)
{ // constructor
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = field.mesh().sieveMesh();
- const ALE::Obj<Mesh::RealSection>& fieldSection = field.section();
+ DM dmMesh = field.dmMesh();
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Jacobian");
@@ -43,8 +42,7 @@
// dimension, otherwise use a block size of 1.
const int blockFlag = (blockOkay) ? -1 : 1;
- PetscErrorCode err = DMMeshCreateMatrix(sieveMesh, fieldSection,
- matrixType, &_matrix, blockFlag);
+ PetscErrorCode err = DMCreateMatrix(dmMesh, matrixType, &_matrix);
CHECK_PETSC_ERROR_MSG(err, "Could not create PETSc sparse matrix "
"associated with system Jacobian.");
logger.stagePop();
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -196,10 +196,15 @@
const PylithScalar* valsE = _data->valsResidual;
const int sizeE = _data->spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
#if 0
@@ -215,6 +220,7 @@
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
// ----------------------------------------------------------------------
@@ -236,10 +242,15 @@
const PylithScalar* valsE = _data->valsResidualLumped;
const int sizeE = _data->spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
#if 0
@@ -255,6 +266,7 @@
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);
} // testIntegrateResidualLumped
// ----------------------------------------------------------------------
@@ -337,6 +349,7 @@
jacobian.complete();
const PylithScalar* valsE = _data->valsJacobianLumped;
+ const int sizeE = _data->numVertices * _data->spaceDim;
#if 0 // DEBUGGING
// TEMPORARY
@@ -350,11 +363,15 @@
} // for
#endif // DEBUGGING
- const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
- CPPUNIT_ASSERT(!jacobianSection.isNull());
- const PylithScalar* vals = jacobianSection->restrictSpace();
- const int size = jacobianSection->sizeWithBC();
- const int sizeE = _data->numVertices * _data->spaceDim;
+ PetscSection jacobianSection = jacobian.petscSection();
+ Vec jacobianVec = jacobian.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(jacobianSection);CPPUNIT_ASSERT(jacobianVec);
+ err = VecGetArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(jacobianSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
const PylithScalar tolerance = 1.0e-06;
@@ -363,6 +380,7 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ err = VecRestoreArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateJacobianLumped
// ----------------------------------------------------------------------
@@ -395,6 +413,9 @@
CPPUNIT_ASSERT_EQUAL(pylith::PYLITH_MAXSCALAR, stableTimeStep);
} // testStableTimeStep
+extern PetscErrorCode DMComplexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]);
+extern PetscErrorCode DMComplexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]);
+
// ----------------------------------------------------------------------
// Initialize elasticity integrator.
void
@@ -420,6 +441,10 @@
new SieveMesh::sieve_type(mesh->comm());
CPPUNIT_ASSERT(!sieve.isNull());
+ mesh->createDMMesh(_data->cellDim);
+ DM dmMesh = mesh->dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+
// Cells and vertices
const bool interpolate = false;
ALE::Obj<SieveFlexMesh::sieve_type> s =
@@ -436,7 +461,11 @@
sieveMesh->stratify();
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
_data->vertices);
+ PetscErrorCode err;
+ err = DMComplexBuildFromCellList_Private(dmMesh, _data->numCells, _data->numVertices, _data->numBasis, _data->cells);CHECK_PETSC_ERROR(err);
+ err = DMComplexBuildCoordinates_Private(dmMesh, _data->spaceDim, _data->numCells, _data->numVertices, _data->vertices);CHECK_PETSC_ERROR(err);
+
// Material ids
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->heightStratum(0);
@@ -449,7 +478,13 @@
e_iter != cells->end();
++e_iter)
sieveMesh->setValue(labelMaterials, *e_iter, _data->matId);
+ PetscInt cStart, cEnd, c;
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ err = DMComplexSetLabelValue(dmMesh, "material-id", c, _data->matId);CHECK_PETSC_ERROR(err);
+ }
+
// Setup quadrature
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
_data->basisDerivRef, _data->numQuadPts,
@@ -499,46 +534,42 @@
fields->copyLayout("residual");
const int fieldSize = spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& dispIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- const ALE::Obj<RealSection>& dispTmdtSection =
- fields->get("disp(t-dt)").section();
- const ALE::Obj<RealSection>& velSection =
- fields->get("velocity(t)").section();
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
- CPPUNIT_ASSERT(!dispTSection.isNull());
- CPPUNIT_ASSERT(!dispTmdtSection.isNull());
- CPPUNIT_ASSERT(!velSection.isNull());
- CPPUNIT_ASSERT(!accSection.isNull());
scalar_array velVertex(spaceDim);
scalar_array accVertex(spaceDim);
-
- const int offset = _data->numCells;
- for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
- dispIncrSection->updatePoint(iVertex+offset,
- &_data->fieldTIncr[iVertex*spaceDim]);
- dispTSection->updatePoint(iVertex+offset,
- &_data->fieldT[iVertex*spaceDim]);
- dispTmdtSection->updatePoint(iVertex+offset,
- &_data->fieldTmdt[iVertex*spaceDim]);
-
- for (int iDim=0; iDim < spaceDim; ++iDim) {
+ PetscSection dispTSectionP = fields->get("disp(t)").petscSection();
+ Vec dispTVec = fields->get("disp(t)").localVector();
+ CPPUNIT_ASSERT(dispTSectionP);CPPUNIT_ASSERT(dispTVec);
+ PetscSection dispTIncrSectionP = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ CPPUNIT_ASSERT(dispTIncrSectionP);CPPUNIT_ASSERT(dispTIncrVec);
+ PetscSection dispTmdtSectionP = fields->get("disp(t-dt)").petscSection();
+ Vec dispTmdtVec = fields->get("disp(t-dt)").localVector();
+ CPPUNIT_ASSERT(dispTmdtSectionP);CPPUNIT_ASSERT(dispTmdtVec);
+ PetscSection velSectionP = fields->get("velocity(t)").petscSection();
+ Vec velVec = fields->get("velocity(t)").localVector();
+ CPPUNIT_ASSERT(velSectionP);CPPUNIT_ASSERT(velVec);
+ PetscSection accSectionP = fields->get("acceleration(t)").petscSection();
+ Vec accVec = fields->get("acceleration(t)").localVector();
+ CPPUNIT_ASSERT(accSectionP);CPPUNIT_ASSERT(accVec);
+ PetscInt offset = _data->numCells;
+ for(int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
+ for(int iDim=0; iDim < spaceDim; ++iDim) {
velVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] +
- _data->fieldT[iVertex*spaceDim+iDim] -
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2.0*dt);
+ _data->fieldT[iVertex*spaceDim+iDim] -
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2.0*dt);
accVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] -
- _data->fieldT[iVertex*spaceDim+iDim] +
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
+ _data->fieldT[iVertex*spaceDim+iDim] +
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
} // for
- velSection->updatePoint(iVertex+offset, &velVertex[0]);
- accSection->updatePoint(iVertex+offset, &accVertex[0]);
+ err = DMComplexVecSetClosure(dmMesh, dispTSectionP, dispTVec, iVertex+offset, &_data->fieldT[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTIncrSectionP, dispTIncrVec, iVertex+offset, &_data->fieldTIncr[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTmdtSectionP, dispTmdtVec, iVertex+offset, &_data->fieldTmdt[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, velSectionP, velVec, iVertex+offset, &velVertex[0], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, accSectionP, accVec, iVertex+offset, &accVertex[0], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
} // for
+
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitLgDeform.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitLgDeform.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -115,10 +115,15 @@
std::cout << " " << valsE[i] << "\n";
#endif
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-04;
@@ -127,6 +132,7 @@
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
// ----------------------------------------------------------------------
@@ -148,10 +154,15 @@
const PylithScalar* valsE = _data->valsResidualLumped;
const int sizeE = _data->spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
#if 0
@@ -167,6 +178,7 @@
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);
} // testIntegrateResidualLumped
// ----------------------------------------------------------------------
@@ -249,6 +261,7 @@
jacobian.complete();
const PylithScalar* valsE = _data->valsJacobianLumped;
+ const int sizeE = _data->numVertices * _data->spaceDim;
#if 0 // DEBUGGING
// TEMPORARY
@@ -262,11 +275,15 @@
} // for
#endif // DEBUGGING
- const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
- CPPUNIT_ASSERT(!jacobianSection.isNull());
- const PylithScalar* vals = jacobianSection->restrictSpace();
- const int size = jacobianSection->sizeWithBC();
- const int sizeE = _data->numVertices * _data->spaceDim;
+ PetscSection jacobianSection = jacobian.petscSection();
+ Vec jacobianVec = jacobian.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(jacobianSection);CPPUNIT_ASSERT(jacobianVec);
+ err = VecGetArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(jacobianSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
const PylithScalar tolerance = 1.0e-06;
@@ -275,6 +292,7 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ err = VecRestoreArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateJacobianLumped
// ----------------------------------------------------------------------
@@ -293,6 +311,9 @@
integrator.updateStateVars(t, &fields);
} // testUpdateStateVars
+extern PetscErrorCode DMComplexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]);
+extern PetscErrorCode DMComplexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]);
+
// ----------------------------------------------------------------------
// Initialize elasticity integrator.
void
@@ -322,6 +343,10 @@
new SieveMesh::sieve_type(mesh->comm());
CPPUNIT_ASSERT(!sieve.isNull());
+ mesh->createDMMesh(_data->cellDim);
+ DM dmMesh = mesh->dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+
// Cells and vertices
const bool interpolate = false;
ALE::Obj<SieveFlexMesh::sieve_type> s =
@@ -338,7 +363,11 @@
sieveMesh->stratify();
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
_data->vertices);
+ PetscErrorCode err;
+ err = DMComplexBuildFromCellList_Private(dmMesh, _data->numCells, _data->numVertices, _data->numBasis, _data->cells);CHECK_PETSC_ERROR(err);
+ err = DMComplexBuildCoordinates_Private(dmMesh, _data->spaceDim, _data->numCells, _data->numVertices, _data->vertices);CHECK_PETSC_ERROR(err);
+
// Material ids
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->heightStratum(0);
@@ -351,7 +380,13 @@
e_iter != cells->end();
++e_iter)
sieveMesh->setValue(labelMaterials, *e_iter, _data->matId);
+ PetscInt cStart, cEnd, c;
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ err = DMComplexSetLabelValue(dmMesh, "material-id", c, _data->matId);CHECK_PETSC_ERROR(err);
+ }
+
// Setup quadrature
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
_data->basisDerivRef, _data->numQuadPts,
@@ -395,45 +430,39 @@
residual.zero();
fields->copyLayout("residual");
- const int fieldSize = spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& dispIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- const ALE::Obj<RealSection>& dispTmdtSection =
- fields->get("disp(t-dt)").section();
- const ALE::Obj<RealSection>& velSection =
- fields->get("velocity(t)").section();
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
- CPPUNIT_ASSERT(!dispTSection.isNull());
- CPPUNIT_ASSERT(!dispTmdtSection.isNull());
- CPPUNIT_ASSERT(!velSection.isNull());
- CPPUNIT_ASSERT(!accSection.isNull());
-
scalar_array velVertex(spaceDim);
scalar_array accVertex(spaceDim);
-
const int offset = _data->numCells;
- for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
- dispIncrSection->updatePoint(iVertex+offset,
- &_data->fieldTIncr[iVertex*spaceDim]);
- dispTSection->updatePoint(iVertex+offset,
- &_data->fieldT[iVertex*spaceDim]);
- dispTmdtSection->updatePoint(iVertex+offset,
- &_data->fieldTmdt[iVertex*spaceDim]);
- for (int iDim=0; iDim < spaceDim; ++iDim) {
+ PetscSection dispTSectionP = fields->get("disp(t)").petscSection();
+ Vec dispTVec = fields->get("disp(t)").localVector();
+ CPPUNIT_ASSERT(dispTSectionP);CPPUNIT_ASSERT(dispTVec);
+ PetscSection dispTIncrSectionP = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ CPPUNIT_ASSERT(dispTIncrSectionP);CPPUNIT_ASSERT(dispTIncrVec);
+ PetscSection dispTmdtSectionP = fields->get("disp(t-dt)").petscSection();
+ Vec dispTmdtVec = fields->get("disp(t-dt)").localVector();
+ CPPUNIT_ASSERT(dispTmdtSectionP);CPPUNIT_ASSERT(dispTmdtVec);
+ PetscSection velSectionP = fields->get("velocity(t)").petscSection();
+ Vec velVec = fields->get("velocity(t)").localVector();
+ CPPUNIT_ASSERT(velSectionP);CPPUNIT_ASSERT(velVec);
+ PetscSection accSectionP = fields->get("acceleration(t)").petscSection();
+ Vec accVec = fields->get("acceleration(t)").localVector();
+ CPPUNIT_ASSERT(accSectionP);CPPUNIT_ASSERT(accVec);
+ for(int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
+ for(int iDim=0; iDim < spaceDim; ++iDim) {
velVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] +
- _data->fieldT[iVertex*spaceDim+iDim] -
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2.0*dt);
+ _data->fieldT[iVertex*spaceDim+iDim] -
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2.0*dt);
accVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] -
- _data->fieldT[iVertex*spaceDim+iDim] +
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
+ _data->fieldT[iVertex*spaceDim+iDim] +
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
} // for
- velSection->updatePoint(iVertex+offset, &velVertex[0]);
- accSection->updatePoint(iVertex+offset, &accVertex[0]);
+ err = DMComplexVecSetClosure(dmMesh, dispTSectionP, dispTVec, iVertex+offset, &_data->fieldT[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTIncrSectionP, dispTIncrVec, iVertex+offset, &_data->fieldTIncr[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTmdtSectionP, dispTmdtVec, iVertex+offset, &_data->fieldTmdt[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, velSectionP, velVec, iVertex+offset, &velVertex[0], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, accSectionP, accVec, iVertex+offset, &accVertex[0], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
} // for
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTet4.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -191,10 +191,15 @@
const PylithScalar* valsE = _data->valsResidual;
const int sizeE = _data->spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
#if 0 // DEBUGGING
@@ -210,6 +215,7 @@
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
// ----------------------------------------------------------------------
@@ -231,10 +237,15 @@
const PylithScalar* valsE = _data->valsResidualLumped;
const int sizeE = _data->spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
#if 0 // DEBUGGING
@@ -250,6 +261,7 @@
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);
} // testIntegrateResidualLumped
// ----------------------------------------------------------------------
@@ -353,10 +365,15 @@
std::cout << " " << valsE[i] << "\n";
#endif // DEBUGGING
- const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
- CPPUNIT_ASSERT(!jacobianSection.isNull());
- const PylithScalar* vals = jacobianSection->restrictSpace();
- const int size = jacobianSection->sizeWithBC();
+ PetscSection jacobianSection = jacobian.petscSection();
+ Vec jacobianVec = jacobian.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(jacobianSection);CPPUNIT_ASSERT(jacobianVec);
+ err = VecGetArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(jacobianSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
const PylithScalar tolerance = 1.0e-06;
@@ -365,6 +382,7 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ err = VecRestoreArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateJacobianLumped
// ----------------------------------------------------------------------
@@ -397,6 +415,9 @@
CPPUNIT_ASSERT_EQUAL(pylith::PYLITH_MAXSCALAR, stableTimeStep);
} // testStableTimeStep
+extern PetscErrorCode DMComplexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]);
+extern PetscErrorCode DMComplexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]);
+
// ----------------------------------------------------------------------
// Initialize elasticity integrator.
void
@@ -422,6 +443,10 @@
new SieveMesh::sieve_type(mesh->comm());
CPPUNIT_ASSERT(!sieve.isNull());
+ mesh->createDMMesh(_data->cellDim);
+ DM dmMesh = mesh->dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+
// Cells and vertices
const bool interpolate = false;
ALE::Obj<SieveFlexMesh::sieve_type> s =
@@ -438,7 +463,11 @@
sieveMesh->stratify();
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
_data->vertices);
+ PetscErrorCode err;
+ err = DMComplexBuildFromCellList_Private(dmMesh, _data->numCells, _data->numVertices, _data->numBasis, _data->cells);CHECK_PETSC_ERROR(err);
+ err = DMComplexBuildCoordinates_Private(dmMesh, _data->spaceDim, _data->numCells, _data->numVertices, _data->vertices);CHECK_PETSC_ERROR(err);
+
// Material ids
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->heightStratum(0);
@@ -451,7 +480,13 @@
e_iter != cells->end();
++e_iter)
sieveMesh->setValue(labelMaterials, *e_iter, _data->matId);
+ PetscInt cStart, cEnd, c;
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ err = DMComplexSetLabelValue(dmMesh, "material-id", c, _data->matId);CHECK_PETSC_ERROR(err);
+ }
+
// Setup quadrature
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
_data->basisDerivRef, _data->numQuadPts,
@@ -500,48 +535,39 @@
residual.zero();
fields->copyLayout("residual");
- const int fieldSize = spaceDim * _data->numVertices;
- topology::Field<topology::Mesh>& dispIncr = fields->get("dispIncr(t->t+dt)");
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- topology::Field<topology::Mesh>& dispTmdt = fields->get("disp(t-dt)");
- const ALE::Obj<RealSection>& dispIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- const ALE::Obj<RealSection>& dispTmdtSection =
- fields->get("disp(t-dt)").section();
- const ALE::Obj<RealSection>& velSection =
- fields->get("velocity(t)").section();
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
- CPPUNIT_ASSERT(!dispTSection.isNull());
- CPPUNIT_ASSERT(!dispTmdtSection.isNull());
- CPPUNIT_ASSERT(!velSection.isNull());
- CPPUNIT_ASSERT(!accSection.isNull());
-
scalar_array velVertex(spaceDim);
scalar_array accVertex(spaceDim);
-
const int offset = _data->numCells;
- for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
- dispIncrSection->updatePoint(iVertex+offset,
- &_data->fieldTIncr[iVertex*spaceDim]);
- dispTSection->updatePoint(iVertex+offset,
- &_data->fieldT[iVertex*spaceDim]);
- dispTmdtSection->updatePoint(iVertex+offset,
- &_data->fieldTmdt[iVertex*spaceDim]);
- for (int iDim=0; iDim < spaceDim; ++iDim) {
+ PetscSection dispTSectionP = fields->get("disp(t)").petscSection();
+ Vec dispTVec = fields->get("disp(t)").localVector();
+ CPPUNIT_ASSERT(dispTSectionP);CPPUNIT_ASSERT(dispTVec);
+ PetscSection dispTIncrSectionP = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ CPPUNIT_ASSERT(dispTIncrSectionP);CPPUNIT_ASSERT(dispTIncrVec);
+ PetscSection dispTmdtSectionP = fields->get("disp(t-dt)").petscSection();
+ Vec dispTmdtVec = fields->get("disp(t-dt)").localVector();
+ CPPUNIT_ASSERT(dispTmdtSectionP);CPPUNIT_ASSERT(dispTmdtVec);
+ PetscSection velSectionP = fields->get("velocity(t)").petscSection();
+ Vec velVec = fields->get("velocity(t)").localVector();
+ CPPUNIT_ASSERT(velSectionP);CPPUNIT_ASSERT(velVec);
+ PetscSection accSectionP = fields->get("acceleration(t)").petscSection();
+ Vec accVec = fields->get("acceleration(t)").localVector();
+ CPPUNIT_ASSERT(accSectionP);CPPUNIT_ASSERT(accVec);
+ for(int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
+ for(int iDim=0; iDim < spaceDim; ++iDim) {
velVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] +
- _data->fieldT[iVertex*spaceDim+iDim] -
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2.0*dt);
+ _data->fieldT[iVertex*spaceDim+iDim] -
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2.0*dt);
accVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] -
- _data->fieldT[iVertex*spaceDim+iDim] +
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
+ _data->fieldT[iVertex*spaceDim+iDim] +
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
} // for
- velSection->updatePoint(iVertex+offset, &velVertex[0]);
- accSection->updatePoint(iVertex+offset, &accVertex[0]);
+ err = DMComplexVecSetClosure(dmMesh, dispTSectionP, dispTVec, iVertex+offset, &_data->fieldT[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTIncrSectionP, dispTIncrVec, iVertex+offset, &_data->fieldTIncr[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTmdtSectionP, dispTmdtVec, iVertex+offset, &_data->fieldTmdt[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, velSectionP, velVec, iVertex+offset, &velVertex[0], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, accSectionP, accVec, iVertex+offset, &accVertex[0], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
} // for
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicitTri3.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -191,10 +191,15 @@
const PylithScalar* valsE = _data->valsResidual;
const int sizeE = _data->spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
#if 0 // DEBUGGING
@@ -210,6 +215,7 @@
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
// ----------------------------------------------------------------------
@@ -231,10 +237,15 @@
const PylithScalar* valsE = _data->valsResidualLumped;
const int sizeE = _data->spaceDim * _data->numVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
#if 0 // DEBUGGING
@@ -250,6 +261,7 @@
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);
} // testIntegrateResidualLumped
// ----------------------------------------------------------------------
@@ -353,10 +365,15 @@
std::cout << " " << valsE[i] << "\n";
#endif // DEBUGGING
- const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
- CPPUNIT_ASSERT(!jacobianSection.isNull());
- const PylithScalar* vals = jacobianSection->restrictSpace();
- const int size = jacobianSection->sizeWithBC();
+ PetscSection jacobianSection = jacobian.petscSection();
+ Vec jacobianVec = jacobian.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(jacobianSection);CPPUNIT_ASSERT(jacobianVec);
+ err = VecGetArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(jacobianSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
const PylithScalar tolerance = 1.0e-06;
@@ -365,6 +382,7 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ err = VecRestoreArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateJacobianLumped
// ----------------------------------------------------------------------
@@ -397,6 +415,9 @@
CPPUNIT_ASSERT_EQUAL(pylith::PYLITH_MAXSCALAR, stableTimeStep);
} // testStableTimeStep
+extern PetscErrorCode DMComplexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]);
+extern PetscErrorCode DMComplexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]);
+
// ----------------------------------------------------------------------
// Initialize elasticity integrator.
void
@@ -422,6 +443,10 @@
new SieveMesh::sieve_type(mesh->comm());
CPPUNIT_ASSERT(!sieve.isNull());
+ mesh->createDMMesh(_data->cellDim);
+ DM dmMesh = mesh->dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+
// Cells and vertices
const bool interpolate = false;
ALE::Obj<SieveFlexMesh::sieve_type> s =
@@ -438,7 +463,11 @@
sieveMesh->stratify();
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
_data->vertices);
+ PetscErrorCode err;
+ err = DMComplexBuildFromCellList_Private(dmMesh, _data->numCells, _data->numVertices, _data->numBasis, _data->cells);CHECK_PETSC_ERROR(err);
+ err = DMComplexBuildCoordinates_Private(dmMesh, _data->spaceDim, _data->numCells, _data->numVertices, _data->vertices);CHECK_PETSC_ERROR(err);
+
// Material ids
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->heightStratum(0);
@@ -451,7 +480,13 @@
e_iter != cells->end();
++e_iter)
sieveMesh->setValue(labelMaterials, *e_iter, _data->matId);
+ PetscInt cStart, cEnd, c;
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ err = DMComplexSetLabelValue(dmMesh, "material-id", c, _data->matId);CHECK_PETSC_ERROR(err);
+ }
+
// Setup quadrature
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
_data->basisDerivRef, _data->numQuadPts,
@@ -500,48 +535,39 @@
residual.zero();
fields->copyLayout("residual");
- const int fieldSize = spaceDim * _data->numVertices;
- topology::Field<topology::Mesh>& dispIncr = fields->get("dispIncr(t->t+dt)");
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- topology::Field<topology::Mesh>& dispTmdt = fields->get("disp(t-dt)");
- const ALE::Obj<RealSection>& dispIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- const ALE::Obj<RealSection>& dispTmdtSection =
- fields->get("disp(t-dt)").section();
- const ALE::Obj<RealSection>& velSection =
- fields->get("velocity(t)").section();
- const ALE::Obj<RealSection>& accSection =
- fields->get("acceleration(t)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
- CPPUNIT_ASSERT(!dispTSection.isNull());
- CPPUNIT_ASSERT(!dispTmdtSection.isNull());
- CPPUNIT_ASSERT(!velSection.isNull());
- CPPUNIT_ASSERT(!accSection.isNull());
-
scalar_array velVertex(spaceDim);
scalar_array accVertex(spaceDim);
-
const int offset = _data->numCells;
- for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
- dispIncrSection->updatePoint(iVertex+offset,
- &_data->fieldTIncr[iVertex*spaceDim]);
- dispTSection->updatePoint(iVertex+offset,
- &_data->fieldT[iVertex*spaceDim]);
- dispTmdtSection->updatePoint(iVertex+offset,
- &_data->fieldTmdt[iVertex*spaceDim]);
- for (int iDim=0; iDim < spaceDim; ++iDim) {
+ PetscSection dispTSectionP = fields->get("disp(t)").petscSection();
+ Vec dispTVec = fields->get("disp(t)").localVector();
+ CPPUNIT_ASSERT(dispTSectionP);CPPUNIT_ASSERT(dispTVec);
+ PetscSection dispTIncrSectionP = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ CPPUNIT_ASSERT(dispTIncrSectionP);CPPUNIT_ASSERT(dispTIncrVec);
+ PetscSection dispTmdtSectionP = fields->get("disp(t-dt)").petscSection();
+ Vec dispTmdtVec = fields->get("disp(t-dt)").localVector();
+ CPPUNIT_ASSERT(dispTmdtSectionP);CPPUNIT_ASSERT(dispTmdtVec);
+ PetscSection velSectionP = fields->get("velocity(t)").petscSection();
+ Vec velVec = fields->get("velocity(t)").localVector();
+ CPPUNIT_ASSERT(velSectionP);CPPUNIT_ASSERT(velVec);
+ PetscSection accSectionP = fields->get("acceleration(t)").petscSection();
+ Vec accVec = fields->get("acceleration(t)").localVector();
+ CPPUNIT_ASSERT(accSectionP);CPPUNIT_ASSERT(accVec);
+ for(int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
+ for(int iDim=0; iDim < spaceDim; ++iDim) {
velVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] +
- _data->fieldT[iVertex*spaceDim+iDim] -
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2.0*dt);
+ _data->fieldT[iVertex*spaceDim+iDim] -
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2.0*dt);
accVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] -
- _data->fieldT[iVertex*spaceDim+iDim] +
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
+ _data->fieldT[iVertex*spaceDim+iDim] +
+ _data->fieldTmdt[iVertex*spaceDim+iDim]) / (dt*dt);
} // for
- velSection->updatePoint(iVertex+offset, &velVertex[0]);
- accSection->updatePoint(iVertex+offset, &accVertex[0]);
+ err = DMComplexVecSetClosure(dmMesh, dispTSectionP, dispTVec, iVertex+offset, &_data->fieldT[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTIncrSectionP, dispTIncrVec, iVertex+offset, &_data->fieldTIncr[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTmdtSectionP, dispTmdtVec, iVertex+offset, &_data->fieldTmdt[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, velSectionP, velVec, iVertex+offset, &velVertex[0], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, accSectionP, accVec, iVertex+offset, &accVertex[0], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
} // for
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -392,29 +392,14 @@
residual.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
residual.allocate();
residual.zero();
- std::cout << residual.label() << " section " << residual.petscSection() << std::endl;
fields->copyLayout("residual");
- const int fieldSize = _data->spaceDim * _data->numVertices;
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- const ALE::Obj<RealSection>& dispTSection = dispT.section();
- CPPUNIT_ASSERT(!dispTSection.isNull());
- topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
- const ALE::Obj<RealSection>& dispTIncrSection = dispTIncr.section();
- CPPUNIT_ASSERT(!dispTIncrSection.isNull());
const int offset = _data->numCells;
- for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
- dispTSection->updatePoint(iVertex+offset,
- &_data->fieldT[iVertex*_data->spaceDim]);
- dispTIncrSection->updatePoint(iVertex+offset,
- &_data->fieldTIncr[iVertex*_data->spaceDim]);
- } // for
-
- PetscSection dispTSectionP = dispT.petscSection();
- Vec dispTVec = dispT.localVector();
+ PetscSection dispTSectionP = fields->get("disp(t)").petscSection();
+ Vec dispTVec = fields->get("disp(t)").localVector();
CPPUNIT_ASSERT(dispTSectionP);CPPUNIT_ASSERT(dispTVec);
- PetscSection dispTIncrSectionP = dispTIncr.petscSection();
- Vec dispTIncrVec = dispTIncr.localVector();
+ PetscSection dispTIncrSectionP = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
CPPUNIT_ASSERT(dispTIncrSectionP);CPPUNIT_ASSERT(dispTIncrVec);
for(int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
err = DMComplexVecSetClosure(dmMesh, dispTSectionP, dispTVec, iVertex+offset, &_data->fieldT[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicitLgDeform.cc 2012-09-06 15:09:35 UTC (rev 20693)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicitLgDeform.cc 2012-09-06 15:12:20 UTC (rev 20694)
@@ -115,10 +115,15 @@
std::cout << " " << valsE[i] << "\n";
#endif
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(residualSection);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-04;
@@ -127,6 +132,7 @@
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
// ----------------------------------------------------------------------
@@ -203,6 +209,9 @@
integrator.updateStateVars(t, &fields);
} // testUpdateStateVars
+extern PetscErrorCode DMComplexBuildFromCellList_Private(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const int cells[]);
+extern PetscErrorCode DMComplexBuildCoordinates_Private(DM dm, PetscInt spaceDim, PetscInt numCells, PetscInt numVertices, const double vertexCoords[]);
+
// ----------------------------------------------------------------------
// Initialize elasticity integrator.
void
@@ -229,6 +238,10 @@
new SieveMesh::sieve_type(mesh->comm());
CPPUNIT_ASSERT(!sieve.isNull());
+ mesh->createDMMesh(_data->cellDim);
+ DM dmMesh = mesh->dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+
// Cells and vertices
const bool interpolate = false;
ALE::Obj<SieveFlexMesh::sieve_type> s =
@@ -245,7 +258,11 @@
sieveMesh->stratify();
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, _data->spaceDim,
_data->vertices);
+ PetscErrorCode err;
+ err = DMComplexBuildFromCellList_Private(dmMesh, _data->numCells, _data->numVertices, _data->numBasis, _data->cells);CHECK_PETSC_ERROR(err);
+ err = DMComplexBuildCoordinates_Private(dmMesh, _data->spaceDim, _data->numCells, _data->numVertices, _data->vertices);CHECK_PETSC_ERROR(err);
+
// Material ids
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->heightStratum(0);
@@ -258,7 +275,13 @@
e_iter != cells->end();
++e_iter)
sieveMesh->setValue(labelMaterials, *e_iter, _data->matId);
+ PetscInt cStart, cEnd;
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ err = DMComplexSetLabelValue(dmMesh, "material-id", c, _data->matId);CHECK_PETSC_ERROR(err);
+ }
+
// Setup quadrature
_quadrature->initialize(_data->basis, _data->numQuadPts, _data->numBasis,
_data->basisDerivRef, _data->numQuadPts,
@@ -299,19 +322,16 @@
residual.zero();
fields->copyLayout("residual");
- const int fieldSize = _data->spaceDim * _data->numVertices;
- topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
- const ALE::Obj<RealSection>& dispTSection = dispT.section();
- CPPUNIT_ASSERT(!dispTSection.isNull());
- topology::Field<topology::Mesh>& dispTIncr = fields->get("dispIncr(t->t+dt)");
- const ALE::Obj<RealSection>& dispTIncrSection = dispTIncr.section();
- CPPUNIT_ASSERT(!dispTIncrSection.isNull());
const int offset = _data->numCells;
- for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
- dispTSection->updatePoint(iVertex+offset,
- &_data->fieldT[iVertex*_data->spaceDim]);
- dispTIncrSection->updatePoint(iVertex+offset,
- &_data->fieldTIncr[iVertex*_data->spaceDim]);
+ PetscSection dispTSectionP = fields->get("disp(t)").petscSection();
+ Vec dispTVec = fields->get("disp(t)").localVector();
+ CPPUNIT_ASSERT(dispTSectionP);CPPUNIT_ASSERT(dispTVec);
+ PetscSection dispTIncrSectionP = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ CPPUNIT_ASSERT(dispTIncrSectionP);CPPUNIT_ASSERT(dispTIncrVec);
+ for(int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
+ err = DMComplexVecSetClosure(dmMesh, dispTSectionP, dispTVec, iVertex+offset, &_data->fieldT[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
+ err = DMComplexVecSetClosure(dmMesh, dispTIncrSectionP, dispTIncrVec, iVertex+offset, &_data->fieldTIncr[iVertex*_data->spaceDim], INSERT_ALL_VALUES);CHECK_PETSC_ERROR(err);
} // for
} // _initialize
More information about the CIG-COMMITS
mailing list