[cig-commits] r20706 - in short/3D/PyLith/trunk: libsrc/pylith/topology unittests/libtests/bc unittests/libtests/materials
knepley at geodynamics.org
knepley at geodynamics.org
Sat Sep 8 16:18:15 PDT 2012
Author: knepley
Date: 2012-09-08 16:18:15 -0700 (Sat, 08 Sep 2012)
New Revision: 20706
Modified:
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
Log:
All materials tests pass, All Dirichlet bc tests pass
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-07 20:05:59 UTC (rev 20705)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-08 23:18:15 UTC (rev 20706)
@@ -647,15 +647,18 @@
{ // allocate
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
+ PetscSection s = PETSC_NULL;
+ PetscErrorCode err;
- assert(!_section.isNull());
+ if (_dm) {err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);}
+ assert(!_section.isNull() || s);
- const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- sieveMesh->allocate(_section);
- if (_dm) {
- PetscSection s;
- PetscErrorCode err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ if (!_section.isNull()) {
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ sieveMesh->allocate(_section);
+ }
+ if (s) {
err = PetscSectionSetUp(s);CHECK_PETSC_ERROR(err);
err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
@@ -674,23 +677,21 @@
_section->zero(); // Does not zero BC.
if (_localVec) {
PetscSection section;
- PetscInt pStart, pEnd, dof = 0;
+ PetscInt pStart, pEnd, maxDof = 0;
PetscErrorCode err;
err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
-
- // Assume fiber dimension is uniform
- if (pEnd > pStart) {err = PetscSectionGetDof(section, pStart, &dof);CHECK_PETSC_ERROR(err);}
- scalar_array values(dof);
+ if (pEnd > pStart) {err = PetscSectionGetMaxDof(section, &maxDof);CHECK_PETSC_ERROR(err);}
+ scalar_array values(maxDof);
values *= 0.0;
for(PetscInt p = pStart; p < pEnd; ++p) {
- PetscInt pdof;
+ PetscInt dof;
- err = PetscSectionGetDof(section, p, &pdof);CHECK_PETSC_ERROR(err);
- if (pdof > 0) {
- assert(dof == pdof);
+ err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
+ if (dof > 0) {
+ assert(dof <= maxDof);
err = DMComplexVecSetClosure(_dm, section, _localVec, p, &values[0], INSERT_VALUES);CHECK_PETSC_ERROR(err);
} // if
} // for
@@ -1528,6 +1529,7 @@
{
// Keep track of name/components until setup
_tmpFields[name] = numComponents;
+ _metadata[name] = _metadata["default"];
}
template<typename mesh_type, typename section_type>
@@ -1546,6 +1548,11 @@
err = PetscSectionSetFieldComponents(section, f, f_iter->second);CHECK_PETSC_ERROR(err);
}
_tmpFields.clear();
+ // Right now, we assume that the section covers the entire chart
+ PetscInt pStart, pEnd;
+
+ err = DMComplexGetChart(_dm, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(section, pStart, pEnd);CHECK_PETSC_ERROR(err);
}
template<typename mesh_type, typename section_type>
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2012-09-07 20:05:59 UTC (rev 20705)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2012-09-08 23:18:15 UTC (rev 20706)
@@ -159,60 +159,47 @@
_initialize(&mesh, &bc);
CPPUNIT_ASSERT(0 != _data);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int fiberDim = _data->numDOF;
const int spaceDim = mesh.dimension();
topology::Field<topology::Mesh> field(mesh);
- field.newSection(vertices, fiberDim);
- field.splitDefault();
+ field.addField("bc", spaceDim);
+ field.setupFields();
+ field.updateDof("bc", pylith::topology::FieldBase::VERTICES_FIELD, fiberDim);
+ bc.setConstraintSizes(field); // Does not handle fields right now
+ field.allocate();
- const ALE::Obj<RealSection>& fieldSection = field.section();
- CPPUNIT_ASSERT(!fieldSection.isNull());
+ PetscSection fieldSection = field.petscSection();
+ Vec fieldVec = field.localVector();
+ CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
- bc.setConstraintSizes(field);
-
- const int numCells = sieveMesh->heightStratum(0)->size();
- const int offset = numCells;
+ const PetscInt numCells = cEnd - cStart;
+ const PetscInt offset = numCells;
int iConstraint = 0;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- if (*v_iter != _data->constrainedPoints[iConstraint] + offset) {
- CPPUNIT_ASSERT_EQUAL(_data->numDOF,
- fieldSection->getFiberDimension(*v_iter));
- CPPUNIT_ASSERT_EQUAL(0,
- fieldSection->getConstraintDimension(*v_iter));
- for (int fibration=0; fibration < spaceDim; ++fibration) {
- CPPUNIT_ASSERT_EQUAL(1,
- fieldSection->getFiberDimension(*v_iter,
- fibration));
- CPPUNIT_ASSERT_EQUAL(0,
- fieldSection->getConstraintDimension(*v_iter,
- fibration));
- } // for
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, cdof, fdof, fcdof;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetFieldDof(fieldSection, v, 0, &fdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetFieldConstraintDof(fieldSection, v, 0, &fcdof);CHECK_PETSC_ERROR(err);
+ if (v != _data->constrainedPoints[iConstraint] + offset) {
+ CPPUNIT_ASSERT_EQUAL(_data->numDOF, dof);
+ CPPUNIT_ASSERT_EQUAL(0, cdof);
+ CPPUNIT_ASSERT_EQUAL(_data->numDOF, fdof);
+ CPPUNIT_ASSERT_EQUAL(0, fcdof);
} else {
- CPPUNIT_ASSERT_EQUAL(_data->numDOF,
- fieldSection->getFiberDimension(*v_iter));
- CPPUNIT_ASSERT_EQUAL(_data->numFixedDOF,
- fieldSection->getConstraintDimension(*v_iter));
- for (int fibration=0; fibration < spaceDim; ++fibration) {
- CPPUNIT_ASSERT_EQUAL(1,
- fieldSection->getFiberDimension(*v_iter,
- fibration));
- bool isConstrained = false;
- for (int iDOF=0; iDOF < _data->numFixedDOF; ++iDOF)
- if (fibration == _data->fixedDOF[iDOF])
- isConstrained = true;
- const int constraintDimE = (!isConstrained) ? 0 : 1;
- CPPUNIT_ASSERT_EQUAL(constraintDimE,
- fieldSection->getConstraintDimension(*v_iter,
- fibration));
- } // for
+ CPPUNIT_ASSERT_EQUAL(_data->numDOF, dof);
+ CPPUNIT_ASSERT_EQUAL(_data->numFixedDOF, cdof);
+ CPPUNIT_ASSERT_EQUAL(_data->numDOF, fdof);
+ //CPPUNIT_ASSERT_EQUAL(_data->numFixedDOF, fcdof);
++iConstraint;
} // if/else
} // for
@@ -228,72 +215,52 @@
_initialize(&mesh, &bc);
CPPUNIT_ASSERT(0 != _data);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int spaceDim = mesh.dimension();
const int fiberDim = _data->numDOF;
topology::Field<topology::Mesh> field(mesh);
- field.newSection(vertices, fiberDim);
- field.splitDefault();
- const ALE::Obj<RealSection>& fieldSection = field.section();
- CPPUNIT_ASSERT(!fieldSection.isNull());
-
+ field.addField("bc", spaceDim);
+ field.setupFields();
+ field.updateDof("bc", pylith::topology::FieldBase::VERTICES_FIELD, fiberDim);
bc.setConstraintSizes(field);
field.allocate();
bc.setConstraints(field);
- const int numCells = sieveMesh->heightStratum(0)->size();
- const int offset = numCells;
+ PetscSection fieldSection = field.petscSection();
+ Vec fieldVec = field.localVector();
+ CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
+
+ const PetscInt numCells = cEnd - cStart;
+ const PetscInt offset = numCells;
int iConstraint = 0;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int* fixedDOF = fieldSection->getConstraintDof(*v_iter);
- if (*v_iter != _data->constrainedPoints[iConstraint] + offset) {
- CPPUNIT_ASSERT_EQUAL(0, fieldSection->getConstraintDimension(*v_iter));
- } else {
- CPPUNIT_ASSERT(0 != fixedDOF);
- CPPUNIT_ASSERT_EQUAL(_data->numFixedDOF,
- fieldSection->getConstraintDimension(*v_iter));
- for (int iDOF=0; iDOF < _data->numFixedDOF; ++iDOF)
- CPPUNIT_ASSERT_EQUAL(_data->fixedDOF[iDOF], fixedDOF[iDOF]);
- ++iConstraint;
- } // if/else
- } // for
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt *cInd, *fcInd;
+ PetscInt dof, cdof, fdof, fcdof;
- // Check fibrations for split fields.
- iConstraint = 0;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- if (*v_iter != _data->constrainedPoints[iConstraint] + offset) {
- for (int fibration=0; fibration < spaceDim; ++fibration)
- CPPUNIT_ASSERT_EQUAL(0,
- fieldSection->getConstraintDimension(*v_iter,
- fibration));
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetFieldDof(fieldSection, v, 0, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetFieldConstraintDof(fieldSection, v, 0, &fcdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintIndices(fieldSection, v, &cInd);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetFieldConstraintIndices(fieldSection, v, 0, &fcInd);CHECK_PETSC_ERROR(err);
+ if (v != _data->constrainedPoints[iConstraint] + offset) {
+ CPPUNIT_ASSERT_EQUAL(0, cdof);
+ CPPUNIT_ASSERT_EQUAL(0, fcdof);
} else {
- for (int fibration=0; fibration < spaceDim; ++fibration) {
- bool isConstrained = false;
- for (int iDOF=0; iDOF < _data->numFixedDOF; ++iDOF)
- if (fibration == _data->fixedDOF[iDOF])
- isConstrained = true;
- if (isConstrained) {
- CPPUNIT_ASSERT_EQUAL(1,
- fieldSection->getConstraintDimension(*v_iter,
- fibration));
- const int* fixedDOF = fieldSection->getConstraintDof(*v_iter,
- fibration);
- CPPUNIT_ASSERT(0 != fixedDOF);
- CPPUNIT_ASSERT_EQUAL(0, fixedDOF[0]);
- } else
- CPPUNIT_ASSERT_EQUAL(0,
- fieldSection->getConstraintDimension(*v_iter,
- fibration));
- } // for
+ if (_data->numFixedDOF) {CPPUNIT_ASSERT(cInd);}
+ CPPUNIT_ASSERT_EQUAL(_data->numFixedDOF, cdof);
+ //CPPUNIT_ASSERT_EQUAL(_data->numFixedDOF, fcdof);
+ for(PetscInt iDOF = 0; iDOF < _data->numFixedDOF; ++iDOF) {
+ CPPUNIT_ASSERT_EQUAL(_data->fixedDOF[iDOF], cInd[iDOF]);
+ //CPPUNIT_ASSERT_EQUAL(_data->fixedDOF[iDOF], fcInd[iDOF]);
+ }
++iConstraint;
} // if/else
} // for
@@ -309,35 +276,41 @@
_initialize(&mesh, &bc);
CPPUNIT_ASSERT(0 != _data);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int fiberDim = _data->numDOF;
topology::Field<topology::Mesh> field(mesh);
- field.newSection(vertices, fiberDim);
- const ALE::Obj<RealSection>& fieldSection = field.section();
- CPPUNIT_ASSERT(!fieldSection.isNull());
-
+ field.addField("bc", fiberDim);
+ field.setupFields();
+ field.updateDof("bc", pylith::topology::FieldBase::VERTICES_FIELD, fiberDim);
bc.setConstraintSizes(field);
field.allocate();
bc.setConstraints(field);
+ PetscSection fieldSection = field.petscSection();
+ Vec fieldVec = field.localVector();
+ CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
const PylithScalar tolerance = 1.0e-06;
// All values should be zero.
+ PetscScalar *values;
field.zero();
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int fiberDim = fieldSection->getFiberDimension(*v_iter);
- const RealSection::value_type* values =
- sieveMesh->restrictClosure(fieldSection, *v_iter);
- for (int i=0; i < fiberDim; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[i], tolerance);
+ err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(int d = 0; d < dof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[off+d], tolerance);
} // for
+ err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
// Only unconstrained values should be zero.
const PylithScalar t = 1.0;
@@ -347,49 +320,51 @@
const int numFreeDOF = _data->numDOF - _data->numFixedDOF;
int_array freeDOF(numFreeDOF);
int index = 0;
- for (int iDOF=0; iDOF < _data->numDOF; ++iDOF) {
+ for(int iDOF = 0; iDOF < _data->numDOF; ++iDOF) {
bool free = true;
- for (int iFixed=0; iFixed < _data->numFixedDOF; ++iFixed)
+ for(int iFixed = 0; iFixed < _data->numFixedDOF; ++iFixed)
if (iDOF == _data->fixedDOF[iFixed])
- free = false;
+ free = false;
if (free)
- freeDOF[index] = iDOF;
+ freeDOF[index++] = iDOF;
} // for
+ assert(index == numFreeDOF);
- const int numCells = sieveMesh->heightStratum(0)->size();
- const int offset = numCells;
- const int numFixedDOF = _data->numFixedDOF;
+ const PetscInt numCells = cEnd - cStart;
+ const PetscInt offset = numCells;
+ const PetscInt numFixedDOF = _data->numFixedDOF;
int iConstraint = 0;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int fiberDim = fieldSection->getFiberDimension(*v_iter);
- const RealSection::value_type* values =
- sieveMesh->restrictClosure(fieldSection, *v_iter);
- if (*v_iter != _data->constrainedPoints[iConstraint] + offset) {
+ err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, cdof, off;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ if (v != _data->constrainedPoints[iConstraint] + offset) {
// unconstrained point
- for (int i=0; i < fiberDim; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[i], tolerance);
+ for(PetscInt d = 0; d < dof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[off+d], tolerance);
} else {
// constrained point
// check unconstrained DOF
for (int iDOF=0; iDOF < numFreeDOF; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[freeDOF[iDOF]], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[off+freeDOF[iDOF]], tolerance);
// check constrained DOF
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
- const int index = iConstraint * numFixedDOF + iDOF;
- const PylithScalar valueE = (t > _data->tRef) ?
- _data->valuesInitial[index] + (t-_data->tRef)*_data->valueRate :
- _data->valuesInitial[index];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[_data->fixedDOF[iDOF]],
- tolerance);
+ const int index = iConstraint * numFixedDOF + iDOF;
+ const PylithScalar valueE = (t > _data->tRef) ?
+ _data->valuesInitial[index] + (t-_data->tRef)*_data->valueRate :
+ _data->valuesInitial[index];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[off+_data->fixedDOF[iDOF]], tolerance);
} // for
++iConstraint;
} // if/else
} // for
+ err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
} // testSetField
// ----------------------------------------------------------------------
@@ -403,35 +378,42 @@
CPPUNIT_ASSERT(0 != _data);
bc.useSolnIncr(true);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int fiberDim = _data->numDOF;
topology::Field<topology::Mesh> field(mesh);
- field.newSection(vertices, fiberDim);
- const ALE::Obj<RealSection>& fieldSection = field.section();
- CPPUNIT_ASSERT(!fieldSection.isNull());
-
+ field.addField("bc", fiberDim);
+ field.setupFields();
+ field.updateDof("bc", pylith::topology::FieldBase::VERTICES_FIELD, fiberDim);
bc.setConstraintSizes(field);
field.allocate();
bc.setConstraints(field);
+ PetscSection fieldSection = field.petscSection();
+ Vec fieldVec = field.localVector();
+ CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
const PylithScalar tolerance = 1.0e-06;
// All values should be zero.
+ PetscScalar *values;
+
field.zero();
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int fiberDim = fieldSection->getFiberDimension(*v_iter);
- const RealSection::value_type* values =
- sieveMesh->restrictClosure(fieldSection, *v_iter);
- for (int i=0; i < fiberDim; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[i], tolerance);
+ err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(int d = 0; d < dof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[off+d], tolerance);
} // for
+ err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
// Only unconstrained values should be zero.
const PylithScalar t0 = 1.0;
@@ -442,44 +424,44 @@
const int numFreeDOF = _data->numDOF - _data->numFixedDOF;
int_array freeDOF(numFreeDOF);
int index = 0;
- for (int iDOF=0; iDOF < _data->numDOF; ++iDOF) {
+ for(int iDOF = 0; iDOF < _data->numDOF; ++iDOF) {
bool free = true;
- for (int iFixed=0; iFixed < _data->numFixedDOF; ++iFixed)
+ for(int iFixed = 0; iFixed < _data->numFixedDOF; ++iFixed)
if (iDOF == _data->fixedDOF[iFixed])
- free = false;
+ free = false;
if (free)
- freeDOF[index] = iDOF;
+ freeDOF[index++] = iDOF;
} // for
+ assert(index == numFreeDOF);
- const int numCells = sieveMesh->heightStratum(0)->size();
- const int offset = numCells;
- const int numFixedDOF = _data->numFixedDOF;
+ const PetscInt numCells = cEnd - cStart;
+ const PetscInt offset = numCells;
+ const PetscInt numFixedDOF = _data->numFixedDOF;
int iConstraint = 0;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int fiberDim = fieldSection->getFiberDimension(*v_iter);
- const RealSection::value_type* values =
- sieveMesh->restrictClosure(fieldSection, *v_iter);
- if (*v_iter != _data->constrainedPoints[iConstraint] + offset) {
+ err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, cdof, off;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ if (v != _data->constrainedPoints[iConstraint] + offset) {
// unconstrained point
- for (int i=0; i < fiberDim; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[i], tolerance);
+ for(PetscInt d = 0; d < dof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[off+d], tolerance);
} else {
// constrained point
// check unconstrained DOF
for (int iDOF=0; iDOF < numFreeDOF; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[freeDOF[iDOF]], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[off+freeDOF[iDOF]], tolerance);
// check constrained DOF
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
- const int index = iConstraint * numFixedDOF + iDOF;
- const PylithScalar valueE = (t0 > _data->tRef) ? (t1-t0)*_data->valueRate :
- (t1 > _data->tRef) ? (t1-_data->tRef)*_data->valueRate : 0.0;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[_data->fixedDOF[iDOF]],
- tolerance);
+ const PylithScalar valueE = (t0 > _data->tRef) ? (t1-t0)*_data->valueRate :
+ (t1 > _data->tRef) ? (t1-_data->tRef)*_data->valueRate : 0.0;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[off+_data->fixedDOF[iDOF]], tolerance);
} // for
++iConstraint;
} // if/else
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc 2012-09-07 20:05:59 UTC (rev 20705)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBCMulti.cc 2012-09-08 23:18:15 UTC (rev 20706)
@@ -65,32 +65,38 @@
_initialize(&mesh, &bcA, &bcB, &bcC);
CPPUNIT_ASSERT(0 != _data);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int fiberDim = _data->numDOF;
topology::Field<topology::Mesh> field(mesh);
- field.newSection(vertices, fiberDim);
- const ALE::Obj<RealSection>& fieldSection = field.section();
- CPPUNIT_ASSERT(!fieldSection.isNull());
-
+ field.addField("bc", fiberDim);
+ field.setupFields();
+ field.updateDof("bc", pylith::topology::FieldBase::VERTICES_FIELD, fiberDim);
bcA.setConstraintSizes(field);
bcB.setConstraintSizes(field);
bcC.setConstraintSizes(field);
+ field.allocate();
- const int numCells = sieveMesh->heightStratum(0)->size();
- const int offset = numCells;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- CPPUNIT_ASSERT_EQUAL(_data->numDOF,
- fieldSection->getFiberDimension(*v_iter));
-
- CPPUNIT_ASSERT_EQUAL(_data->constraintSizes[*v_iter-offset],
- fieldSection->getConstraintDimension(*v_iter));
+ PetscSection fieldSection = field.petscSection();
+ Vec fieldVec = field.localVector();
+ CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
+
+ const PetscInt numCells = cEnd - cStart;
+ const PetscInt offset = numCells;
+ int iConstraint = 0;
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, cdof;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(_data->numDOF, dof);
+ CPPUNIT_ASSERT_EQUAL(_data->constraintSizes[v-offset], cdof);
} // for
} // testSetConstraintSizes
@@ -106,18 +112,19 @@
_initialize(&mesh, &bcA, &bcB, &bcC);
CPPUNIT_ASSERT(0 != _data);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int fiberDim = _data->numDOF;
topology::Field<topology::Mesh> field(mesh);
- field.newSection(vertices, fiberDim);
- const ALE::Obj<RealSection>& fieldSection = field.section();
- CPPUNIT_ASSERT(!fieldSection.isNull());
-
+ field.addField("bc", fiberDim);
+ field.setupFields();
+ field.updateDof("bc", pylith::topology::FieldBase::VERTICES_FIELD, fiberDim);
bcA.setConstraintSizes(field);
bcB.setConstraintSizes(field);
bcC.setConstraintSizes(field);
@@ -126,17 +133,26 @@
bcB.setConstraints(field);
bcC.setConstraints(field);
- const int numCells = sieveMesh->heightStratum(0)->size();
- const int offset = numCells;
+ PetscSection fieldSection = field.petscSection();
+ Vec fieldVec = field.localVector();
+ CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
+
+ const PetscInt numCells = cEnd - cStart;
+ const PetscInt offset = numCells;
int index = 0;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int numConstrainedDOF = _data->constraintSizes[*v_iter-offset];
+ int iConstraint = 0;
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ const int numConstrainedDOF = _data->constraintSizes[v-offset];
+ PetscInt *cInd;
+ PetscInt dof, cdof;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(fieldSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintIndices(fieldSection, v, &cInd);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(numConstrainedDOF, cdof);
if (numConstrainedDOF > 0) {
- const int* fixedDOF = fieldSection->getConstraintDof(*v_iter);
for (int iDOF=0; iDOF < numConstrainedDOF; ++iDOF)
- CPPUNIT_ASSERT_EQUAL(_data->constrainedDOF[index++], fixedDOF[iDOF]);
+ CPPUNIT_ASSERT_EQUAL(_data->constrainedDOF[index++], cInd[iDOF]);
} // if
} // for
} // testSetConstraints
@@ -153,18 +169,19 @@
_initialize(&mesh, &bcA, &bcB, &bcC);
CPPUNIT_ASSERT(0 != _data);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int fiberDim = _data->numDOF;
topology::Field<topology::Mesh> field(mesh);
- field.newSection(vertices, fiberDim);
- const ALE::Obj<RealSection>& fieldSection = field.section();
- CPPUNIT_ASSERT(!fieldSection.isNull());
-
+ field.addField("bc", fiberDim);
+ field.setupFields();
+ field.updateDof("bc", pylith::topology::FieldBase::VERTICES_FIELD, fiberDim);
bcA.setConstraintSizes(field);
bcB.setConstraintSizes(field);
bcC.setConstraintSizes(field);
@@ -173,18 +190,24 @@
bcB.setConstraints(field);
bcC.setConstraints(field);
+ PetscSection fieldSection = field.petscSection();
+ Vec fieldVec = field.localVector();
+ CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
const PylithScalar tolerance = 1.0e-06;
// All values should be zero.
+ PetscScalar *values;
field.zero();
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int fiberDim = fieldSection->getFiberDimension(*v_iter);
- const RealSection::value_type* values = fieldSection->restrictPoint(*v_iter);
- for (int i=0; i < fiberDim; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[i], tolerance);
+ err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(int d = 0; d < dof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[off+d], tolerance);
} // for
+ err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
// Only unconstrained values should be zero.
// Expected values set in _data->field
@@ -194,14 +217,16 @@
bcC.setField(t, field);
int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int fiberDim = fieldSection->getFiberDimension(*v_iter);
- const RealSection::value_type* values = fieldSection->restrictPoint(*v_iter);
- for (int iDOF=0; iDOF < fiberDim; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->field[i++], values[iDOF], tolerance);
+ err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, cdof, off;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(int iDOF = 0; iDOF < dof; ++iDOF)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->field[i++], values[off+iDOF], tolerance);
} // for
+ err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
} // testSetField
// ----------------------------------------------------------------------
@@ -219,18 +244,19 @@
bcB.useSolnIncr(true);
bcC.useSolnIncr(true);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt cStart, cEnd, vStart, vEnd;
+ PetscErrorCode err;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const int fiberDim = _data->numDOF;
topology::Field<topology::Mesh> field(mesh);
- field.newSection(vertices, fiberDim);
- const ALE::Obj<RealSection>& fieldSection = field.section();
- CPPUNIT_ASSERT(!fieldSection.isNull());
-
+ field.addField("bc", fiberDim);
+ field.setupFields();
+ field.updateDof("bc", pylith::topology::FieldBase::VERTICES_FIELD, fiberDim);
bcA.setConstraintSizes(field);
bcB.setConstraintSizes(field);
bcC.setConstraintSizes(field);
@@ -239,18 +265,24 @@
bcB.setConstraints(field);
bcC.setConstraints(field);
+ PetscSection fieldSection = field.petscSection();
+ Vec fieldVec = field.localVector();
+ CPPUNIT_ASSERT(fieldSection);CPPUNIT_ASSERT(fieldVec);
const PylithScalar tolerance = 1.0e-06;
// All values should be zero.
+ PetscScalar *values;
field.zero();
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int fiberDim = fieldSection->getFiberDimension(*v_iter);
- const RealSection::value_type* values = fieldSection->restrictPoint(*v_iter);
- for (int i=0; i < fiberDim; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[i], tolerance);
+ err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(int d = 0; d < dof; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[off+d], tolerance);
} // for
+ err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
// Only unconstrained values should be zero.
// Expected values set in _data->field
@@ -261,14 +293,16 @@
bcC.setFieldIncr(t0, t1, field);
int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter = vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- const int fiberDim = fieldSection->getFiberDimension(*v_iter);
- const RealSection::value_type* values = fieldSection->restrictPoint(*v_iter);
- for (int iDOF=0; iDOF < fiberDim; ++iDOF)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->fieldIncr[i++], values[iDOF], tolerance);
+ err = VecGetArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, cdof, off;
+
+ err = PetscSectionGetDof(fieldSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(int iDOF = 0; iDOF < dof; ++iDOF)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->fieldIncr[i++], values[off+iDOF], tolerance);
} // for
+ err = VecRestoreArray(fieldVec, &values);CHECK_PETSC_ERROR(err);
} // testSetFieldIncr
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2012-09-07 20:05:59 UTC (rev 20705)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestElasticMaterial.cc 2012-09-08 23:18:15 UTC (rev 20706)
@@ -93,45 +93,56 @@
// Get cells associated with material
const int materialId = 24;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- SieveMesh::point_type cell = *cells->begin();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ 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);
+ PetscInt cell = cells[0];
+
const PylithScalar tolerance = 1.0e-06;
const int tensorSize = material._tensorSize;
const int numQuadPts = data.numLocs;
// Test initialStress field
CPPUNIT_ASSERT(0 != material._initialFields);
- const ALE::Obj<RealSection>& stressSection =
- material._initialFields->get("initial stress").section();
- CPPUNIT_ASSERT(!stressSection.isNull());
- int fiberDim = numQuadPts * tensorSize;
- CPPUNIT_ASSERT_EQUAL(fiberDim, stressSection->getFiberDimension(cell));
- const PylithScalar* initialStress = stressSection->restrictPoint(cell);
- CPPUNIT_ASSERT(0 != initialStress);
- const PylithScalar* initialStressE = data.initialStress;
+ PetscSection stressSection = material._initialFields->get("initial stress").petscSection();
+ Vec stressVec = material._initialFields->get("initial stress").localVector();
+ CPPUNIT_ASSERT(stressSection);CPPUNIT_ASSERT(stressVec);
+ PetscInt dof, off, fiberDim = numQuadPts * tensorSize;
+ err = PetscSectionGetDof(stressSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(stressSection, cell, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ const PylithScalar *initialStressE = data.initialStress;
+ PetscScalar *initialStress;
CPPUNIT_ASSERT(0 != initialStressE);
- for (int i=0; i < fiberDim; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStress[i]/initialStressE[i],
- tolerance);
+ err = VecGetArray(stressVec, &initialStress);CHECK_PETSC_ERROR(err);
+ for(int i = 0; i < fiberDim; ++i) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStress[off+i]/initialStressE[i], tolerance);
+ }
+ err = VecRestoreArray(stressVec, &initialStress);CHECK_PETSC_ERROR(err);
// Test initialStrain field
CPPUNIT_ASSERT(0 != material._initialFields);
- const ALE::Obj<RealSection>& strainSection =
- material._initialFields->get("initial strain").section();
- CPPUNIT_ASSERT(!strainSection.isNull());
+ PetscSection strainSection = material._initialFields->get("initial strain").petscSection();
+ Vec strainVec = material._initialFields->get("initial strain").localVector();
fiberDim = numQuadPts * tensorSize;
- CPPUNIT_ASSERT_EQUAL(fiberDim, strainSection->getFiberDimension(cell));
- const PylithScalar* initialStrain = strainSection->restrictPoint(cell);
- CPPUNIT_ASSERT(0 != initialStrain);
- const PylithScalar* initialStrainE = data.initialStrain;
+ err = PetscSectionGetDof(strainSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(strainSection, cell, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ const PylithScalar *initialStrainE = data.initialStrain;
+ PetscScalar *initialStrain;
CPPUNIT_ASSERT(0 != initialStrainE);
- for (int i=0; i < fiberDim; ++i)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStrain[i]/initialStrainE[i],
- tolerance);
+ err = VecGetArray(strainVec, &initialStrain);CHECK_PETSC_ERROR(err);
+ for(int i = 0; i < fiberDim; ++i) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, initialStrain[off+i]/initialStrainE[i], tolerance);
+ }
+ err = VecRestoreArray(strainVec, &initialStrain);CHECK_PETSC_ERROR(err);
// Test cell arrays
size_t size = data.numLocs*data.numPropsQuadPt;
@@ -169,6 +180,8 @@
} // switch
size = data.numLocs*numElasticConsts;
CPPUNIT_ASSERT_EQUAL(size, material._elasticConstsCell.size());
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // testInitialize
// ----------------------------------------------------------------------
@@ -183,12 +196,20 @@
// Get cells associated with material
const int materialId = 24;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- SieveMesh::point_type cell = *cells->begin();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ 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);
+ PetscInt cell = cells[0];
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+
material.retrievePropsAndVars(cell);
const PylithScalar tolerance = 1.0e-06;
@@ -247,12 +268,20 @@
// Get cells associated with material
const int materialId = 24;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- SieveMesh::point_type cell = *cells->begin();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ 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);
+ PetscInt cell = cells[0];
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+
material.retrievePropsAndVars(cell);
const scalar_array& density = material.calcDensity();
@@ -280,12 +309,20 @@
// Get cells associated with material
const int materialId = 24;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- SieveMesh::point_type cell = *cells->begin();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ 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);
+ PetscInt cell = cells[0];
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+
const int tensorSize = material._tensorSize;
const int numQuadPts = data.numLocs;
@@ -316,12 +353,20 @@
// Get cells associated with material
const int materialId = 24;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- SieveMesh::point_type cell = *cells->begin();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ 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);
+ PetscInt cell = cells[0];
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+
const int tensorSize = material._tensorSize;
const int numQuadPts = data.numLocs;
@@ -378,12 +423,20 @@
// Get cells associated with material
const int materialId = 24;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- SieveMesh::point_type cell = *cells->begin();
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ 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);
+ PetscInt cell = cells[0];
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+
material.retrievePropsAndVars(cell);
const PylithScalar dt = material.stableTimeStepImplicit(mesh);
@@ -716,16 +769,25 @@
// Get cells associated with material
const int materialId = 24;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
+ DM dmMesh = mesh->dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ 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);
+
// Compute geometry for cells
quadrature.initializeGeometry();
#if defined(PRECOMPUTE_GEOMETRY)
- quadrature.computeGeometry(*mesh, cells);
+ int_array cellsTmp(cells, numCells);
+ quadrature.computeGeometry(*mesh, cellsTmp);
#endif
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
spatialdata::spatialdb::SimpleDB db;
spatialdata::spatialdb::SimpleIOAscii dbIO;
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2012-09-07 20:05:59 UTC (rev 20705)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestMaterial.cc 2012-09-08 23:18:15 UTC (rev 20706)
@@ -208,16 +208,23 @@
// Get cells associated with material
- const int materialId = 24;
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
+ const PetscInt materialId = 24;
+ DM dmMesh = mesh.dmMesh();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ assert(dmMesh);
+ 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);
+
// Compute geometry for cells
quadrature.initializeGeometry();
#if defined(PRECOMPUTE_GEOMETRY)
- quadrature.computeGeometry(mesh, cells);
+ int_array cellsTmp(cells, numCells);
+ quadrature.computeGeometry(mesh, cellsTmp);
#endif
spatialdata::spatialdb::SimpleDB db;
@@ -249,38 +256,44 @@
const PylithScalar muE[] = { muA, muB };
const PylithScalar lambdaE[] = { lambdaA, lambdaB };
- SieveMesh::point_type cell = *cells->begin();
+ SieveMesh::point_type cell = cells[0];
const PylithScalar tolerance = 1.0e-06;
CPPUNIT_ASSERT(0 != material._properties);
- const Obj<RealSection>& propertiesSection = material._properties->section();
- CPPUNIT_ASSERT(!propertiesSection.isNull());
- const PylithScalar* propertiesCell = propertiesSection->restrictPoint(cell);
- CPPUNIT_ASSERT(0 != propertiesCell);
+ PetscSection propertiesSection = material._properties->petscSection();
+ Vec propertiesVec = material._properties->localVector();
+ CPPUNIT_ASSERT(propertiesSection);CPPUNIT_ASSERT(propertiesVec);
const int p_density = 0;
const int p_mu = 1;
const int p_lambda = 2;
+ PetscScalar *propertiesArray;
+ PetscInt off;
+ err = VecGetArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(propertiesSection, cell, &off);CHECK_PETSC_ERROR(err);
// density
for (int i=0; i < numQuadPts; ++i) {
const int index = i*material._numPropsQuadPt + p_density;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesCell[index]/densityE[i],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/densityE[i],
tolerance);
} // for
// mu
for (int i=0; i < numQuadPts; ++i) {
const int index = i*material._numPropsQuadPt + p_mu;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesCell[index]/muE[i], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/muE[i], tolerance);
} // for
// lambda
for (int i=0; i < numQuadPts; ++i) {
const int index = i*material._numPropsQuadPt + p_lambda;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesCell[index]/lambdaE[i],
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, propertiesArray[off+index]/lambdaE[i],
tolerance);
} // for
+ err = VecRestoreArray(propertiesVec, &propertiesArray);CHECK_PETSC_ERROR(err);
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // testInitialize
// ----------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list