[cig-commits] r20686 - in short/3D/PyLith/trunk/libsrc/pylith: bc feassemble
knepley at geodynamics.org
knepley at geodynamics.org
Tue Sep 4 21:00:24 PDT 2012
Author: knepley
Date: 2012-09-04 21:00:23 -0700 (Tue, 04 Sep 2012)
New Revision: 20686
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
Log:
More work on integrators, finished DirichletBC
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2012-09-05 00:37:39 UTC (rev 20685)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2012-09-05 04:00:23 UTC (rev 20686)
@@ -82,9 +82,6 @@
if (0 == numFixedDOF)
return;
- const ALE::Obj<RealSection>& section = field.section();
- assert(!section.isNull());
-
PetscSection sectionP = field.petscSection();
PetscErrorCode err;
assert(sectionP);
@@ -93,22 +90,8 @@
const int numPoints = _points.size();
_offsetLocal.resize(numPoints);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const PetscInt point = _points[iPoint];
- const int fiberDim = section->getFiberDimension(point);
- const int curNumConstraints = section->getConstraintDimension(point);
- if (curNumConstraints + numFixedDOF > fiberDim) {
- std::ostringstream msg;
- msg
- << "Found overly constrained point while setting up constraints for\n"
- << "DirichletBC boundary condition '" << _label << "'.\n"
- << "Number of DOF at point " << point << " is " << fiberDim
- << "\nand number of attempted constraints is "
- << curNumConstraints+numFixedDOF << ".";
- throw std::runtime_error(msg.str());
- } // if
- _offsetLocal[iPoint] = curNumConstraints;
- section->addConstraintDimension(point, numFixedDOF);
- PetscInt dof, cdof;
+ const PetscInt point = _points[iPoint];
+ PetscInt dof, cdof;
err = PetscSectionGetDof(sectionP, point, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetConstraintDof(sectionP, point, &cdof);CHECK_PETSC_ERROR(err);
@@ -126,16 +109,6 @@
// We should be specifying what field the BC is for
//err = PetscSectionAddConstraintFieldDof(sectionP, point, field, numFixedDOF);CHECK_PETSC_ERROR(err);
} // for
-
- // We only worry about the conventional DOF in a split field associated with
- // fibrations.
- const int numFibrations = section->getNumSpaces(); // > 1 for split field
- if (numFibrations > 0)
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
- const int fibration = _bcDOF[iDOF];
- for (int iPoint=0; iPoint < numPoints; ++iPoint)
- section->addConstraintDimension(_points[iPoint], 1, fibration);
- } // for
} // setConstraintSizes
// ----------------------------------------------------------------------
@@ -147,9 +120,6 @@
if (0 == numFixedDOF)
return;
- const ALE::Obj<RealSection>& section = field.section();
- assert(!section.isNull());
-
PetscSection sectionP = field.petscSection();
PetscErrorCode err;
assert(sectionP);
@@ -159,22 +129,19 @@
const SieveMesh::point_type point = _points[iPoint];
// Get list of currently constrained DOF
- const int* curFixedDOF = section->getConstraintDof(point);
- const int numTotalConstrained = section->getConstraintDimension(point);
PetscInt cdof, *cInd;
err = PetscSectionGetConstraintDof(sectionP, point, &cdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetConstraintIndices(sectionP, point, &cInd);CHECK_PETSC_ERROR(err);
// Create array holding all constrained DOF
- int_array allFixedDOF(curFixedDOF, numTotalConstrained);
int_array allCInd(cInd, cdof);
// Verify other BC has not already constrained DOF
const int numPrevious = _offsetLocal[iPoint];
for (int iDOF=0; iDOF < numPrevious; ++iDOF)
for (int jDOF=0; jDOF < numFixedDOF; ++jDOF)
- if ((allFixedDOF[iDOF] == _bcDOF[jDOF]) || (allCInd[iDOF] == _bcDOF[jDOF])) {
+ if (allCInd[iDOF] == _bcDOF[jDOF]) {
std::ostringstream msg;
msg << "Found multiple constraints on degrees of freedom at\n"
<< "point while setting up constraints for DirichletBC\n"
@@ -186,46 +153,23 @@
// Add in the ones for this DirichletBC BC
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
- assert(_offsetLocal[iPoint]+iDOF < numTotalConstrained);
assert(_offsetLocal[iPoint]+iDOF < cdof);
- allFixedDOF[_offsetLocal[iPoint]+iDOF] = _bcDOF[iDOF];
allCInd[_offsetLocal[iPoint]+iDOF] = _bcDOF[iDOF];
} // for
- // Fill in rest of values not yet set (will be set by
- // another DirichletBC BC)
- for (int iDOF=_offsetLocal[iPoint]+numFixedDOF;
- iDOF < numTotalConstrained;
- ++iDOF) {
- assert(iDOF < numTotalConstrained);
- allFixedDOF[iDOF] = 999;
- } // for
+ // Fill in rest of values not yet set (will be set by another DirichletBC BC)
for (int iDOF=_offsetLocal[iPoint]+numFixedDOF; iDOF < cdof; ++iDOF) {
allCInd[iDOF] = 999;
} // for
// Sort list of constrained DOF
// I need these sorted for my update algorithms to work properly
- std::sort(&allFixedDOF[0], &allFixedDOF[numTotalConstrained]);
std::sort(&allCInd[0], &allCInd[cdof]);
// Update list of constrained DOF
- section->setConstraintDof(point, &allFixedDOF[0]);
err = PetscSectionSetConstraintIndices(sectionP, point, &allCInd[0]);CHECK_PETSC_ERROR(err);
//err = PetscSectionSetConstraintFieldIndices(sectionP, point, field, &fieldCInd[0]);CHECK_PETSC_ERROR(err);
} // for
-
- // We only worry about the conventional DOF in a split field associated with
- // fibrations.
- const int numFibrations = section->getNumSpaces(); // > 1 for split field
- if (numFibrations > 0) {
- int zero = 0;
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
- const int fibration = _bcDOF[iDOF];
- for (int iPoint=0; iPoint < numPoints; ++iPoint)
- section->setConstraintDof(_points[iPoint], &zero, fibration);
- } // for
- } // if
} // setConstraints
// ----------------------------------------------------------------------
@@ -253,37 +197,23 @@
assert(valueFiberDim == numFixedDOF);
assert(valueIndex+valueFiberDim <= parametersFiberDim);
- const ALE::Obj<RealSection>& fieldSection = field.section();
- assert(!fieldSection.isNull());
- const int fiberDimension =
- (numPoints > 0) ? fieldSection->getFiberDimension(_points[0]) : 0;
- scalar_array fieldVertex(fiberDimension);
-
- PetscSection fieldSectionP = field.petscSection();
- Vec localVec = field.localVector();
+ PetscSection fieldSectionP = field.petscSection();
+ Vec localVec = field.localVector();
+ PetscScalar *array;
PetscErrorCode err;
assert(fieldSectionP);
err = VecGetArray(localVec, &array);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const SieveMesh::point_type p_bc = _points[iPoint];
-
- assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
- fieldSection->restrictPoint(p_bc, &fieldVertex[0], fieldVertex.size());
+ PetscInt p_bc = _points[iPoint];
PetscInt dof, off;
+ assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
err = PetscSectionGetDof(fieldSectionP, p_bc, &dof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(fieldSectionP, p_bc, &off);CHECK_PETSC_ERROR(err);
+ const PylithScalar *parametersVertex = parametersSection->restrictPoint(p_bc);
- const PylithScalar* parametersVertex = parametersSection->restrictPoint(p_bc);
-
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- fieldVertex[_bcDOF[iDOF]] = parametersVertex[valueIndex+iDOF];
-
- fieldSection->updatePointAll(p_bc, &fieldVertex[0]);
- PetscScalar *array;
-
- for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
array[_bcDOF[iDOF]+off] = parametersVertex[valueIndex+iDOF];
} // for
err = VecRestoreArray(localVec, &array);CHECK_PETSC_ERROR(err);
@@ -317,26 +247,26 @@
assert(valueFiberDim == numFixedDOF);
assert(valueIndex+valueFiberDim <= parametersFiberDim);
- const ALE::Obj<RealSection>& fieldSection = field.section();
- assert(!fieldSection.isNull());
- const int fiberDimension =
- (numPoints > 0) ? fieldSection->getFiberDimension(_points[0]) : 0;
- scalar_array fieldVertex(fiberDimension);
+ PetscSection fieldSectionP = field.petscSection();
+ Vec localVec = field.localVector();
+ PetscScalar *array;
+ PetscErrorCode err;
+ assert(fieldSectionP);
+ err = VecGetArray(localVec, &array);CHECK_PETSC_ERROR(err);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const SieveMesh::point_type p_bc = _points[iPoint];
+ PetscInt p_bc = _points[iPoint];
+ PetscInt dof, off;
- fieldSection->restrictPoint(p_bc, &fieldVertex[0], fieldVertex.size());
-
assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
+ err = PetscSectionGetDof(fieldSectionP, p_bc, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSectionP, p_bc, &off);CHECK_PETSC_ERROR(err);
const PylithScalar* parametersVertex = parametersSection->restrictPoint(p_bc);
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
- fieldVertex[_bcDOF[iDOF]] = parametersVertex[valueIndex+iDOF];
-
- fieldSection->updatePointAll(p_bc, &fieldVertex[0]);
+ array[_bcDOF[iDOF]+off] = parametersVertex[valueIndex+iDOF];
} // for
-
+ err = VecRestoreArray(localVec, &array);CHECK_PETSC_ERROR(err);
} // setFieldIncr
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2012-09-05 00:37:39 UTC (rev 20685)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2012-09-05 04:00:23 UTC (rev 20686)
@@ -108,207 +108,7 @@
_material->useElasticBehavior(!_useSolnIncr);
} // useSolnIncr
-// ----------------------------------------------------------------------
-// Integrate constributions to residual term (r) for operator.
void
-pylith::feassemble::ElasticityImplicit::integrateResidualOld(
- const topology::Field<topology::Mesh>& residual,
- const PylithScalar t,
- topology::SolutionFields* const fields)
-{ // integrateResidualOld
- /// Member prototype for _elasticityResidualXD()
- typedef void (pylith::feassemble::ElasticityImplicit::*elasticityResidual_fn_type)
- (const scalar_array&);
-
- assert(0 != _quadrature);
- assert(0 != _material);
- assert(0 != _logger);
- assert(0 != fields);
-
- const int setupEvent = _logger->eventId("ElIR setup");
- const int geometryEvent = _logger->eventId("ElIR geometry");
- const int computeEvent = _logger->eventId("ElIR compute");
- const int restrictEvent = _logger->eventId("ElIR restrict");
- const int stateVarsEvent = _logger->eventId("ElIR stateVars");
- const int stressEvent = _logger->eventId("ElIR stress");
- const int updateEvent = _logger->eventId("ElIR update");
-
- _logger->eventBegin(setupEvent);
-
- // Get cell geometry information that doesn't depend on cell
- const int numQuadPts = _quadrature->numQuadPts();
- const scalar_array& quadWts = _quadrature->quadWts();
- assert(quadWts.size() == numQuadPts);
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
- const int cellDim = _quadrature->cellDim();
- const int tensorSize = _material->tensorSize();
- if (cellDim != spaceDim)
- throw std::logic_error("Integration for cells with spatial dimensions "
- "different than the spatial dimension of the "
- "domain not implemented yet.");
-
- // Set variables dependent on dimension of cell
- totalStrain_fn_type calcTotalStrainFn;
- elasticityResidual_fn_type elasticityResidualFn;
- if (1 == cellDim) {
- elasticityResidualFn =
- &pylith::feassemble::ElasticityImplicit::_elasticityResidual1D;
- calcTotalStrainFn =
- &pylith::feassemble::IntegratorElasticity::_calcTotalStrain1D;
- } else if (2 == cellDim) {
- elasticityResidualFn =
- &pylith::feassemble::ElasticityImplicit::_elasticityResidual2D;
- calcTotalStrainFn =
- &pylith::feassemble::IntegratorElasticity::_calcTotalStrain2D;
- } else if (3 == cellDim) {
- elasticityResidualFn =
- &pylith::feassemble::ElasticityImplicit::_elasticityResidual3D;
- calcTotalStrainFn =
- &pylith::feassemble::IntegratorElasticity::_calcTotalStrain3D;
- } else
- assert(0);
-
- // Allocate vectors for cell values.
- scalar_array dispTpdtCell(numBasis*spaceDim);
- scalar_array strainCell(numQuadPts*tensorSize);
- strainCell = 0.0;
- scalar_array gravVec(spaceDim);
- 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();
-
- // Get sections
- scalar_array dispTCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
- RestrictVisitor dispTVisitor(*dispTSection, dispTCell.size(), &dispTCell[0]);
-
- scalar_array dispTIncrCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dispTIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- assert(!dispTIncrSection.isNull());
- RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
- dispTIncrCell.size(), &dispTIncrCell[0]);
-
- 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]);
-
- assert(0 != _normalizer);
- const PylithScalar lengthScale = _normalizer->lengthScale();
- const PylithScalar gravityScale =
- _normalizer->pressureScale() / (_normalizer->lengthScale() *
- _normalizer->densityScale());
-
- _logger->eventEnd(setupEvent);
- _logger->eventBegin(computeEvent);
-
- // Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- // 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);
-#endif
-
- // Get state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
-
- // 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);
-
- // Get cell geometry information that depends on cell
- const scalar_array& basis = _quadrature->basis();
- const scalar_array& basisDeriv = _quadrature->basisDeriv();
- const scalar_array& jacobianDet = _quadrature->jacobianDet();
- const scalar_array& quadPtsNondim = _quadrature->quadPts();
-
- // Compute current estimate of displacement at time t+dt using
- // solution increment.
- dispTpdtCell = dispTCell + dispTIncrCell;
-
- // Compute body force vector if gravity is being used.
- if (0 != _gravityField) {
- const spatialdata::geocoords::CoordSys* cs = fields->mesh().coordsys();
- assert(0 != cs);
-
- // Get density at quadrature points for this cell
- const scalar_array& density = _material->calcDensity();
-
- quadPtsGlobal = quadPtsNondim;
- _normalizer->dimensionalize(&quadPtsGlobal[0], quadPtsGlobal.size(),
- lengthScale);
-
- // Compute action for element body forces
- spatialdata::spatialdb::SpatialDB* db = _gravityField;
- for (int iQuad = 0; iQuad < numQuadPts; ++iQuad) {
- const int err = db->query(&gravVec[0], gravVec.size(),
- &quadPtsGlobal[0], spaceDim, cs);
- if (err)
- throw std::runtime_error("Unable to get gravity vector for point.");
- _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];
- for (int iDim = 0; iDim < spaceDim; ++iDim) {
- _cellVector[iBasis * spaceDim + iDim] += valI * gravVec[iDim];
- } // for
- } // for
- } // for
- PetscLogFlops(numQuadPts * (2 + numBasis * (1 + 2 * spaceDim)));
- } // if
-
- // residualSection->view("After gravity contribution");
- // Compute B(transpose) * sigma, first computing strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispTpdtCell,
- numBasis, numQuadPts);
- const scalar_array& stressCell = _material->calcStress(strainCell, true);
-
- CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell);
-
-#if 0 // DEBUGGING
- std::cout << "Updating residual for cell " << *c_iter << std::endl;
- for(int i = 0; i < _quadrature->spaceDim() * _quadrature->numBasis(); ++i) {
- std::cout << " v["<<i<<"]: " << _cellVector[i] << std::endl;
- }
-#endif
- // Assemble cell contribution into field
- residualVisitor.clear();
- sieveMesh->updateClosure(*c_iter, residualVisitor);
- } // for
-
- _logger->eventEnd(computeEvent);
-} // integrateResidualOld
-
-void
pylith::feassemble::ElasticityImplicit::integrateResidual(
const topology::Field<topology::Mesh>& residual,
const PylithScalar t,
@@ -381,8 +181,7 @@
const PetscInt *cells;
PetscInt numCells;
assert(dmMesh);
- const int materialId = _material->id();
- err = DMComplexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
+ 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);
@@ -400,12 +199,14 @@
PetscSection residualSection = residual.petscSection();
Vec residualVec = residual.localVector();
+#if !defined(PRECOMPUTE_GEOMETRY)
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);
+#endif
assert(0 != _normalizer);
const PylithScalar lengthScale = _normalizer->lengthScale();
@@ -499,8 +300,8 @@
// Assemble cell contribution into field
err = DMComplexVecSetClosure(dmMesh, residualSection, residualVec, cell, &_cellVector[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
}
-
err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
_logger->eventEnd(computeEvent);
} // integrateResidual
@@ -693,6 +494,8 @@
err = DMComplexMatSetClosure(dmMesh, dispTSection, 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();
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc 2012-09-05 00:37:39 UTC (rev 20685)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc 2012-09-05 04:00:23 UTC (rev 20686)
@@ -108,16 +108,15 @@
_initializeLogger();
+ // Compute geometry for quadrature operations.
+ _quadrature->initializeGeometry();
+#if defined(PRECOMPUTE_GEOMETRY)
// Get cell information
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
assert(!sieveMesh.isNull());
const int materialId = _material->id();
const ALE::Obj<SieveMesh::label_sequence>& cells =
sieveMesh->getLabelStratum("material-id", materialId);
-
- // Compute geometry for quadrature operations.
- _quadrature->initializeGeometry();
-#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->computeGeometry(mesh, cells);
#endif
@@ -193,61 +192,69 @@
strainCell = 0.0;
// 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 fields
- scalar_array dispCell(numBasis*spaceDim);
- const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
- const ALE::Obj<RealSection>& dispSection = disp.section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ scalar_array dispTCell(numBasis*spaceDim);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
#if !defined(PRECOMPUTE_GEOMETRY)
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);
#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];
// Retrieve 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 physical properties and state variables for cell.
- _material->retrievePropsAndVars(*c_iter);
+ _material->retrievePropsAndVars(cell);
// Restrict input fields to cell
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
+ const PetscScalar *dispTArray;
+ PetscInt dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispTCell[i] = dispTArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Get cell geometry information that depends on cell
const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
- numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTCell, numBasis, numQuadPts);
// Update material state
- _material->updateStateVars(strainCell, *c_iter);
+ _material->updateStateVars(strainCell, cell);
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // updateStateVars
// ----------------------------------------------------------------------
@@ -284,27 +291,34 @@
} // if
const int numCorners = _quadrature->refGeometry().numCorners();
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", _material->id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- const int cellNumCorners = sieveMesh->getNumCellCorners(*c_iter);
+ DM dmMesh = 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);
+
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
+ PetscInt cellNumCorners;
+
+ err = DMComplexGetConeSize(dmMesh, cell, &cellNumCorners);CHECK_PETSC_ERROR(err);
if (numCorners != cellNumCorners) {
std::ostringstream msg;
msg << "Quadrature is incompatible with cell in material '"
<< _material->label() << "'.\n"
- << "Cell " << *c_iter << " has " << cellNumCorners
+ << "Cell " << cell << " has " << cellNumCorners
<< " vertices but quadrature reference cell has "
<< numCorners << " vertices.";
throw std::runtime_error(msg.str());
} // if
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // verifyConfiguration
// ----------------------------------------------------------------------
@@ -329,7 +343,7 @@
assert(0 != fields);
_allocateTensorField(mesh);
topology::Field<topology::Mesh>& buffer =
- _outputFields->get("buffer (tensor)");
+ _outputFields->get("buffer (tensor)");
_material->getField(&buffer, "total_strain");
buffer.addDimensionOkay(true);
return buffer;
@@ -338,7 +352,7 @@
assert(0 != fields);
_allocateTensorField(mesh);
topology::Field<topology::Mesh>& buffer =
- _outputFields->get("buffer (tensor)");
+ _outputFields->get("buffer (tensor)");
buffer.label("total_strain");
buffer.scale(1.0);
buffer.addDimensionOkay(true);
@@ -353,7 +367,7 @@
assert(0 != fields);
_allocateTensorField(mesh);
topology::Field<topology::Mesh>& buffer =
- _outputFields->get("buffer (tensor)");
+ _outputFields->get("buffer (tensor)");
_material->getField(&buffer, "stress");
buffer.addDimensionOkay(true);
return buffer;
@@ -450,13 +464,18 @@
assert(0 != _material);
assert(0 != _outputFields);
- const ALE::Obj<SieveMesh>& sieveMesh = 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());
+ DM dmMesh = 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);
+ int_array cellsTmp(cells, numCells);
+
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
@@ -468,12 +487,14 @@
_outputFields->add("buffer (tensor)", "buffer");
topology::Field<topology::Mesh>& buffer =
_outputFields->get("buffer (tensor)");
- buffer.newSection(cells, numQuadPts*tensorSize);
+ buffer.newSection(cellsTmp, numQuadPts*tensorSize);
buffer.allocate();
buffer.vectorFieldType(topology::FieldBase::MULTI_TENSOR);
buffer.addDimensionOkay(true);
logger.stagePop();
} // if
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // _allocateTensorField
// ----------------------------------------------------------------------
@@ -519,65 +540,74 @@
stressCell = 0.0;
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = field->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 field
- const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
- const ALE::Obj<RealSection>& dispSection = disp.section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ scalar_array dispTCell(numBasis*spaceDim);
+ topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
+ PetscSection dispTSection = dispT.petscSection();
+ Vec dispTVec = dispT.localVector();
+ assert(dispTSection);assert(dispTVec);
#if !defined(PRECOMPUTE_GEOMETRY)
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);
#endif
- const ALE::Obj<RealSection>& fieldSection = field->section();
- assert(!fieldSection.isNull());
+ PetscSection fieldSection = field->petscSection();
+ assert(fieldSection);
// 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];
// Retrieve 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
// Restrict input fields to cell
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
+ const PetscScalar *dispTArray;
+ PetscInt dispTSize;
+ err = DMComplexVecGetClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < dispTSize; ++i) {dispTCell[i] = dispTArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, dispTSection, dispTVec, cell, &dispTSize, &dispTArray);CHECK_PETSC_ERROR(err);
// Get cell geometry information that depends on cell
const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Compute strains
- calcTotalStrainFn(&strainCell, basisDeriv, dispCell,
- numBasis, numQuadPts);
+ calcTotalStrainFn(&strainCell, basisDeriv, dispTCell, numBasis, numQuadPts);
- if (!calcStress)
- fieldSection->updatePoint(*c_iter, &strainCell[0]);
- else {
- _material->retrievePropsAndVars(*c_iter);
+ if (!calcStress) {
+ err = DMComplexVecSetClosure(dmMesh, fieldSection, PETSC_NULL, cell, &strainCell[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+ } else {
+ _material->retrievePropsAndVars(cell);
stressCell = _material->calcStress(strainCell);
- fieldSection->updatePoint(*c_iter, &stressCell[0]);
+ err = DMComplexVecSetClosure(dmMesh, fieldSection, PETSC_NULL, cell, &stressCell[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // else
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // _calcStrainStressField
// ----------------------------------------------------------------------
@@ -602,28 +632,38 @@
stressCell = 0.0;
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = field->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 = field->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 field
- const ALE::Obj<RealSection>& fieldSection = field->section();
- assert(!fieldSection.isNull());
+ PetscSection fieldSection = field->petscSection();
+ Vec fieldVec = field->localVector();
+ assert(fieldSection);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- fieldSection->restrictPoint(*c_iter, &strainCell[0], strainCell.size());
- _material->retrievePropsAndVars(*c_iter);
+ for(PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
+ PetscInt fieldSize;
+ const PetscScalar *fieldArray;
+
+ err = DMComplexVecGetClosure(dmMesh, fieldSection, fieldVec, cell, &fieldSize, &fieldArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < fieldSize; ++i) {strainCell[i] = fieldArray[i];}
+ err = DMComplexVecRestoreClosure(dmMesh, fieldSection, fieldVec, cell, &fieldSize, &fieldArray);CHECK_PETSC_ERROR(err);
+
+ _material->retrievePropsAndVars(cell);
stressCell = _material->calcStress(strainCell);
- fieldSection->updatePoint(*c_iter, &stressCell[0]);
+ err = DMComplexVecSetClosure(dmMesh, fieldSection, PETSC_NULL, cell, &stressCell[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // _calcStressFromStrain
// ----------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list