[cig-commits] r20816 - in short/3D/PyLith/trunk: libsrc/pylith/faults libsrc/pylith/topology unittests/libtests/topology
knepley at geodynamics.org
knepley at geodynamics.org
Tue Oct 9 19:15:08 PDT 2012
Author: knepley
Date: 2012-10-09 19:15:08 -0700 (Tue, 09 Oct 2012)
New Revision: 20816
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldSubMesh.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.cc
Log:
More topology tests passing, fixed more sections on the fault
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-10-09 23:53:54 UTC (rev 20815)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-10-10 02:15:08 UTC (rev 20816)
@@ -1288,12 +1288,14 @@
up[i] = upDir[i];
// Get vertices in fault mesh.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices = faultSieveMesh->depthStratum(0);
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM faultDMMesh = _faultMesh->dmMesh();
+ PetscInt vStart, vEnd, cStart, cEnd;
+ PetscErrorCode err;
+ assert(faultDMMesh);
+ err = DMComplexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetHeightStratum(faultDMMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
// Containers for orientation information.
const int cohesiveDim = _faultMesh->dimension();
const int numBasis = _quadrature->numBasis();
@@ -1320,33 +1322,22 @@
orientation.setupFields();
orientation.newSection(dispRel, orientationSize);
// Create components for along-strike, up-dip, and normal directions
- orientation.updateDof("orientation", pylith::topology::FieldBase::POINTS_FIELD, spaceDim);
+ orientation.updateDof("orientation", pylith::topology::FieldBase::VERTICES_FIELD, spaceDim);
orientation.allocate();
orientation.zero();
PetscSection orientationSection = orientation.petscSection();
Vec orientationVec = orientation.localVector();
PetscScalar *orientationArray;
- PetscErrorCode err;
logger.stagePop();
- // Get fault cells.
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
-
// Compute orientation of fault at constraint vertices
// Get section containing coordinates of vertices
- scalar_array coordinatesCell(numBasis * spaceDim);
- const ALE::Obj<RealSection>& coordinatesSection =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordinatesSection.isNull());
- RestrictVisitor coordinatesVisitor(*coordinatesSection,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
// Set orientation function
assert(cohesiveDim == _quadrature->cellDim());
@@ -1355,46 +1346,46 @@
// Loop over cohesive cells, computing orientation weighted by
// jacobian at constraint vertices
- const ALE::Obj<SieveSubMesh::sieve_type>& sieve = faultSieveMesh->getSieve();
- assert(!sieve.isNull());
- const int closureSize =
- int(pow(sieve->getMaxConeSize(), faultSieveMesh->depth()));
- assert(closureSize >= 0);
- ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type>
- ncV(*sieve, closureSize);
-
err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
- for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin; c_iter
- != cellsEnd; ++c_iter) {
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ PetscInt *closure = PETSC_NULL;
+ PetscInt closureSize, q = 0;
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
+ scalar_array coordinatesCell(numBasis * spaceDim);
+
// Get orientations at fault cell's vertices.
- coordinatesVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordinatesVisitor);
+ err = DMComplexVecGetClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetTransitiveClosure(faultDMMesh, c, PETSC_TRUE, &closureSize, &closure);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
- ncV.clear();
- ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve,
- *c_iter, ncV);
- const int coneSize = ncV.getSize();
- const SieveSubMesh::point_type* cone = ncV.getPoints();
-
- for (int v=0; v < coneSize; ++v) {
+ // Filter out non-vertices
+ for(PetscInt p = 0; p < closureSize*2; p += 2) {
+ if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
+ closure[q*2] = closure[p];
+ closure[q*2+1] = closure[p+1];
+ ++q;
+ }
+ }
+ closureSize = q;
+ for(PetscInt v = 0; v < closureSize; ++v) {
// Compute Jacobian and determinant of Jacobian at vertex
- memcpy(&refCoordsVertex[0], &verticesRef[v * cohesiveDim], cohesiveDim
- * sizeof(PylithScalar));
- cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell,
- refCoordsVertex);
+ memcpy(&refCoordsVertex[0], &verticesRef[v * cohesiveDim], cohesiveDim * sizeof(PylithScalar));
+ cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell, refCoordsVertex);
// Compute orientation
- cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet,
- up);
+ cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet, up);
// Update orientation
PetscInt off;
- err = PetscSectionGetOffset(orientationSection, cone[v], &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, closure[v*2], &off);CHECK_PETSC_ERROR(err);
for(PetscInt d = 0; d < orientationSize; ++d) {
orientationArray[off+d] += orientationVertex[d];
}
} // for
+ err = DMComplexVecRestoreClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ err = DMComplexRestoreTransitiveClosure(faultDMMesh, c, PETSC_TRUE, &closureSize, &closure);CHECK_PETSC_ERROR(err);
} // for
err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
@@ -1407,11 +1398,10 @@
scalar_array vertexDir(orientationSize);
int count = 0;
err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
- for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
- != verticesEnd; ++v_iter, ++count) {
+ for(PetscInt v = vStart; v < vEnd; ++v, ++count) {
PetscInt off;
- err = PetscSectionGetOffset(orientationSection, *v_iter, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, v, &off);CHECK_PETSC_ERROR(err);
for(PetscInt d = 0; d < orientationSize; ++d) {
orientationVertex[d] = orientationArray[off+d];
}
@@ -1432,7 +1422,7 @@
PetscLogFlops(count * orientationSize * 4);
err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
- if (1 == cohesiveDim && vertices->size() > 0) {
+ if (1 == cohesiveDim && vEnd > vStart) {
// Default sense of positive slip is left-lateral and
// fault-opening.
//
@@ -1452,10 +1442,10 @@
// opposite of what we would want, but we cannot flip the fault
// normal direction because it is tied to how the cohesive cells
// are created.
- assert(vertices->size() > 0);
+ assert(vEnd > vStart);
PetscInt off;
- err = PetscSectionGetOffset(orientationSection, *vertices->begin(), &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, vStart, &off);CHECK_PETSC_ERROR(err);
for(PetscInt d = 0; d < orientationSize; ++d) {
orientationVertex[d] = orientationArray[off+d];
}
@@ -1472,13 +1462,11 @@
const int inormal = 2;
if (normalDirDot * shearDirDot < 0.0) {
// Flip shear direction
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
+ for(PetscInt v = vStart; v < vEnd; ++v) {
PetscInt dof;
- err = PetscSectionGetDof(orientationSection, *v_iter, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(orientationSection, *v_iter, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(orientationSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, v, &off);CHECK_PETSC_ERROR(err);
for(PetscInt d = 0; d < orientationSize; ++d) {
orientationVertex[d] = orientationArray[off+d];
}
@@ -1494,7 +1482,7 @@
PetscLogFlops(3 + count * 2);
} // if
- } else if (2 == cohesiveDim && vertices->size() > 0) {
+ } else if (2 == cohesiveDim && vEnd > vStart) {
// Check orientation of first vertex, if dot product of fault
// normal with preferred normal is negative, flip up/down dip
// direction.
@@ -1512,10 +1500,10 @@
// are used are the opposite of what we would want, but we cannot
// flip the fault normal direction because it is tied to how the
// cohesive cells are created.
- assert(vertices->size() > 0);
+ assert(vEnd > vStart);
PetscInt off;
- err = PetscSectionGetOffset(orientationSection, *vertices->begin(), &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, vStart, &off);CHECK_PETSC_ERROR(err);
for(PetscInt d = 0; d < orientationSize; ++d) {
orientationVertex[d] = orientationArray[off+d];
}
@@ -1542,12 +1530,11 @@
// if we have (0,0,-1).
// Flip dip direction
- for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
- != verticesEnd; ++v_iter) {
+ for(PetscInt v = vStart; v < vEnd; ++v) {
PetscInt dof;
- err = PetscSectionGetDof(orientationSection, *v_iter, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(orientationSection, *v_iter, &off);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(orientationSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, v, &off);CHECK_PETSC_ERROR(err);
for(PetscInt d = 0; d < orientationSize; ++d) {
orientationVertex[d] = orientationArray[off+d];
}
@@ -1588,14 +1575,14 @@
scalar_array areaCell(numBasis);
// Get vertices in fault mesh.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveSubMesh::label_sequence::iterator verticesBegin =
- vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM faultDMMesh = _faultMesh->dmMesh();
+ PetscInt vStart, vEnd, cStart, cEnd;
+ PetscErrorCode err;
+ assert(faultDMMesh);
+ err = DMComplexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetHeightStratum(faultDMMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("FaultFields");
@@ -1610,37 +1597,28 @@
const PylithScalar lengthScale = _normalizer->lengthScale();
area.scale(pow(lengthScale, (spaceDim-1)));
area.zero();
- const ALE::Obj<RealSection>& areaSection = area.section();
- assert(!areaSection.isNull());
- UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
+ PetscSection areaSection = area.petscSection();
+ Vec areaVec = area.localVector();
+ PetscScalar *areaArray;
logger.stagePop();
- scalar_array coordinatesCell(numBasis * spaceDim);
- const ALE::Obj<RealSection>& coordinates = faultSieveMesh->getRealSection(
- "coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMComplexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
-
// Loop over cells in fault mesh, compute area
- for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin; c_iter
- != cellsEnd; ++c_iter) {
+ for(PetscInt c = cStart; c < cEnd; ++c) {
areaCell = 0.0;
// Compute geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(c);
#else
- coordsVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
+ err = DMComplexVecGetClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Get cell geometry information that depends on cell
@@ -1655,9 +1633,11 @@
areaCell[iBasis] += dArea;
} // for
} // for
- areaVisitor.clear();
- faultSieveMesh->updateClosure(*c_iter, areaVisitor);
+ err = DMComplexVecSetClosure(faultDMMesh, areaSection, areaVec, c, &areaCell[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
+#if not defined(PRECOMPUTE_GEOMETRY)
+ err = DMComplexVecRestoreClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+#endif
PetscLogFlops( numQuadPts*(1+numBasis*2) );
} // for
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-10-09 23:53:54 UTC (rev 20815)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-10-10 02:15:08 UTC (rev 20816)
@@ -41,8 +41,25 @@
_metadata["default"].scale = 1.0;
_metadata["default"].dimsOkay = false;
if (mesh.dmMesh()) {
- PetscSection s;
- PetscErrorCode err = DMComplexClone(mesh.dmMesh(), &_dm);CHECK_PETSC_ERROR(err);
+ DM dm = mesh.dmMesh();
+ Vec coordVec;
+ PetscSection s;
+ PetscErrorCode err;
+
+ err = DMComplexClone(dm, &_dm);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dm, &coordVec);CHECK_PETSC_ERROR(err);
+ if (coordVec) {
+ DM coordDM, newCoordDM;
+ PetscSection coordSection, newCoordSection;
+
+ err = DMGetCoordinateDM(dm, &coordDM);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinateDM(_dm, &newCoordDM);CHECK_PETSC_ERROR(err);
+ err = DMGetDefaultSection(coordDM, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionClone(coordSection, &newCoordSection);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultSection(newCoordDM, newCoordSection);CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) coordVec);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(_dm, coordVec);CHECK_PETSC_ERROR(err);
+ }
err = PetscSectionCreate(mesh.comm(), &s);CHECK_PETSC_ERROR(err);
err = DMSetDefaultSection(_dm, s);CHECK_PETSC_ERROR(err);
} else {
@@ -64,8 +81,23 @@
assert(!section.isNull());
_metadata["default"] = metadata;
if (mesh.dmMesh()) {
- PetscSection s;
- PetscErrorCode err = DMComplexClone(mesh.dmMesh(), &_dm);CHECK_PETSC_ERROR(err);
+ DM dm = mesh.dmMesh(), coordDM, newCoordDM;
+ PetscSection coordSection, newCoordSection;
+ Vec coordVec;
+ PetscSection s;
+ PetscErrorCode err;
+
+ err = DMComplexClone(dm, &_dm);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dm, &coordVec);CHECK_PETSC_ERROR(err);
+ if (coordVec) {
+ err = DMGetCoordinateDM(dm, &coordDM);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinateDM(_dm, &newCoordDM);CHECK_PETSC_ERROR(err);
+ err = DMGetDefaultSection(coordDM, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionClone(coordSection, &newCoordSection);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultSection(newCoordDM, newCoordSection);CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) coordVec);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(_dm, coordVec);CHECK_PETSC_ERROR(err);
+ }
this->dmSection(&s, &_localVec);
err = DMSetDefaultSection(_dm, s);CHECK_PETSC_ERROR(err);
err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
@@ -84,6 +116,9 @@
_mesh(src._mesh),
_section(PETSC_NULL)
{ // constructor
+ DM dm = mesh.dmMesh(), coordDM, newCoordDM;
+ PetscSection coordSection, newCoordSection;
+ Vec coordVec;
PetscSection s;
PetscErrorCode err;
@@ -94,7 +129,17 @@
err = PetscSectionGetFieldName(s, fields[f], &name);CHECK_PETSC_ERROR(err);
_metadata[name] = src._metadata[name];
}
- err = DMCreateSubDM(_mesh.dmMesh(), numFields, fields, PETSC_NULL, &_dm);CHECK_PETSC_ERROR(err);
+ err = DMCreateSubDM(dm, numFields, fields, PETSC_NULL, &_dm);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dm, &coordVec);CHECK_PETSC_ERROR(err);
+ if (coordVec) {
+ err = DMGetCoordinateDM(dm, &coordDM);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinateDM(_dm, &newCoordDM);CHECK_PETSC_ERROR(err);
+ err = DMGetDefaultSection(coordDM, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionClone(coordSection, &newCoordSection);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultSection(newCoordDM, newCoordSection);CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) coordVec);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(_dm, coordVec);CHECK_PETSC_ERROR(err);
+ }
_globalVec = PETSC_NULL;
_localVec = PETSC_NULL;
} // constructor
@@ -124,7 +169,7 @@
err = DMDestroy(&s_iter->second.dm);CHECK_PETSC_ERROR(err);
err = VecDestroy(&s_iter->second.vector);CHECK_PETSC_ERROR(err);
- err = VecScatterDestroy(&s_iter->second.scatter);CHECK_PETSC_ERROR(err);
+ if (s_iter->second.scatter) {err = VecScatterDestroy(&s_iter->second.scatter);CHECK_PETSC_ERROR(err);}
err = VecDestroy(&s_iter->second.scatterVec);CHECK_PETSC_ERROR(err);
} // for
_scatters.clear();
@@ -509,6 +554,7 @@
++s_iter) {
ScatterInfo sinfo;
sinfo.vector = 0;
+ sinfo.scatter = 0;
sinfo.scatterVec = 0;
// Copy DM
@@ -517,26 +563,29 @@
CHECK_PETSC_ERROR(err);
// Copy scatter
- sinfo.scatter = s_iter->second.scatter;
- err = PetscObjectReference((PetscObject) sinfo.scatter);
- CHECK_PETSC_ERROR(err);
+ if (s_iter->second.scatter) {
+ sinfo.scatter = s_iter->second.scatter;
+ err = PetscObjectReference((PetscObject) sinfo.scatter);
+ CHECK_PETSC_ERROR(err);
+ }
// Create scatter Vec
- assert(s_iter->second.scatterVec);
- int blockSize = 1;
- VecGetBlockSize(s_iter->second.scatterVec, &blockSize);
- if (_section->sizeWithBC() > 0) {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, _section->getStorageSize(),
- _section->restrictSpace(),
- &sinfo.scatterVec);
- CHECK_PETSC_ERROR(err);
- } else {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, 0, PETSC_NULL,
- &sinfo.scatterVec);
- CHECK_PETSC_ERROR(err);
- } // else
+ if (s_iter->second.scatterVec) {
+ int blockSize = 1;
+ err = VecGetBlockSize(s_iter->second.scatterVec, &blockSize);CHECK_PETSC_ERROR(err);
+ if (_section->sizeWithBC() > 0) {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF,
+ blockSize, _section->getStorageSize(),
+ _section->restrictSpace(),
+ &sinfo.scatterVec);
+ CHECK_PETSC_ERROR(err);
+ } else {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF,
+ blockSize, 0, PETSC_NULL,
+ &sinfo.scatterVec);
+ CHECK_PETSC_ERROR(err);
+ } // else
+ }
// Create vector using sizes from source section
int vecLocalSize = 0;
@@ -546,6 +595,8 @@
err = VecGetSize(_globalVec, &vecGlobalSize2);CHECK_PETSC_ERROR(err);
if (vecGlobalSize != vecGlobalSize2) {
+ int blockSize = 1;
+ err = VecGetBlockSize(s_iter->second.vector, &blockSize);CHECK_PETSC_ERROR(err);
err = VecCreate(_mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
err = VecSetSizes(sinfo.vector, vecLocalSize, vecGlobalSize);CHECK_PETSC_ERROR(err);
err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
@@ -882,6 +933,10 @@
} // for
} // if
+ if (_localVec && field._localVec) {
+ PetscErrorCode err = VecCopy(field._localVec, _localVec);CHECK_PETSC_ERROR(err);
+ }
+
label(const_cast<Field&>(field)._metadata["default"].label.c_str()); // Update label
_metadata["default"].scale = const_cast<Field&>(field)._metadata["default"].scale;
} // copy
@@ -984,6 +1039,10 @@
} // for
} // if
+ if (_localVec && field._localVec) {
+ PetscErrorCode err = VecAXPY(_localVec, 1.0, field._localVec);CHECK_PETSC_ERROR(err);
+ }
+
return *this;
} // operator+=
@@ -1001,6 +1060,8 @@
throw std::runtime_error(msg.str());
} // if
+ spatialdata::units::Nondimensional normalizer;
+
if (!_section.isNull()) {
const chart_type& chart = _section->getChart();
const typename chart_type::const_iterator chartEnd = chart.end();
@@ -1010,8 +1071,6 @@
_section->getFiberDimension(*chart.begin()) : 0;
scalar_array values(fiberDim);
- spatialdata::units::Nondimensional normalizer;
-
for (typename chart_type::const_iterator c_iter = chart.begin();
c_iter != chartEnd;
++c_iter)
@@ -1023,6 +1082,27 @@
_section->updatePointAll(*c_iter, &values[0]);
} // if
} // if
+
+ if (_localVec) {
+ PetscSection section;
+ PetscScalar *array;
+ PetscInt pStart, pEnd;
+ PetscErrorCode err;
+
+ err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(_localVec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(section, p, &off);CHECK_PETSC_ERROR(err);
+ if (dof) {
+ normalizer.dimensionalize(&array[off], dof, const_cast<Field*>(this)->_metadata["default"].scale);
+ }
+ }
+ err = VecRestoreArray(_localVec, &array);CHECK_PETSC_ERROR(err);
+ }
} // dimensionalize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc 2012-10-09 23:53:54 UTC (rev 20815)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldMesh.cc 2012-10-10 02:15:08 UTC (rev 20816)
@@ -67,8 +67,8 @@
Field<Mesh> field(mesh);
mesh.createSieveMesh();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(section.isNull());
+ PetscSection section = field.petscSection();
+ CPPUNIT_ASSERT(!section);
} // testSection
// ----------------------------------------------------------------------
@@ -139,16 +139,17 @@
{ // testNewSection
Mesh mesh;
_buildMesh(&mesh);
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
Field<Mesh> field(mesh);
const std::string& label = "field A";
field.label(label.c_str());
field.newSection();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
-
- CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+ PetscSection section = field.petscSection();
+ CPPUNIT_ASSERT(section);
} // testNewSection
// ----------------------------------------------------------------------
@@ -160,27 +161,30 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
Field<Mesh> field(mesh);
const std::string& label = "field A";
field.label(label.c_str());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- field.newSection(vertices, fiberDim);
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ field.newSection(topology::FieldBase::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- CPPUNIT_ASSERT(!vertices.isNull());
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
-
- CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof;
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ }
+ const char *name;
+ err = PetscObjectGetName((PetscObject) vec, &name);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(label, std::string(name));
} // testNewSectionPoints
// ----------------------------------------------------------------------
@@ -192,29 +196,27 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
Field<Mesh> field(mesh);
const std::string& label = "field A";
field.label(label.c_str());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const int npts = vertices->size() / 2;
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ const int npts = (vEnd-vStart) / 2;
int_array pointsIn(npts);
- int_array pointsOut(vertices->size() - npts);
+ int_array pointsOut(vEnd-vStart - npts);
int count = 0;
size_t iIn = 0;
size_t iOut = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
+ for(PetscInt v = vStart; v < vEnd; ++v) {
if (count % 2 == 0)
- pointsIn[iIn++] = *v_iter;
+ pointsIn[iIn++] = v;
else
- pointsOut[iOut++] = *v_iter;
+ pointsOut[iOut++] = v;
++count;
} // for
CPPUNIT_ASSERT_EQUAL(iIn, pointsIn.size());
@@ -222,18 +224,30 @@
field.newSection(pointsIn, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
// Points in array should have a fiber dimension of fiberDim.
- for (int i=0; i < pointsIn.size(); ++i)
- CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(pointsIn[i]));
-
+ for(int i = 0; i < pointsIn.size(); ++i) {
+ PetscInt dof;
+ err = PetscSectionGetDof(section, pointsIn[i], &dof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ }
+
// Points not int array should have a fiber dimension of zero.
- for (int i=0; i < pointsOut.size(); ++i)
- CPPUNIT_ASSERT_EQUAL(0, section->getFiberDimension(pointsOut[i]));
-
- CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+ PetscInt pStart, pEnd;
+ err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ for(int i = 0; i < pointsOut.size(); ++i) {
+ if (pointsOut[i] >= pStart && pointsOut[i] < pEnd) {
+ PetscInt dof;
+ err = PetscSectionGetDof(section, pointsOut[i], &dof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(0, dof);
+ }
+ }
+ const char *name;
+ err = PetscObjectGetName((PetscObject) vec, &name);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(label, std::string(name));
} // testNewSectionPointsArray
// ----------------------------------------------------------------------
@@ -245,25 +259,29 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
Field<Mesh> field(mesh);
const std::string& label = "field A";
field.label(label.c_str());
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
-
- CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof;
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ }
+ const char *name;
+ err = PetscObjectGetName((PetscObject) vec, &name);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(label, std::string(name));
} // testNewSectionDomain
// ----------------------------------------------------------------------
@@ -275,32 +293,34 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
// Create field with atlas to use to create new field
Field<Mesh> fieldSrc(mesh);
fieldSrc.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
- const ALE::Obj<Mesh::RealSection>& sectionSrc = fieldSrc.section();
- CPPUNIT_ASSERT(!sectionSrc.isNull());
- const Mesh::RealSection::chart_type& chart = sectionSrc->getChart();
+ fieldSrc.allocate();
const int fiberDim2 = 5;
Field<Mesh> field(mesh);
const std::string& label = "field A";
field.label(label.c_str());
field.newSection(fieldSrc, fiberDim2);
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- CPPUNIT_ASSERT_EQUAL(fiberDim2, section->getFiberDimension(*v_iter));
-
- CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+ field.allocate();
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof;
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim2, dof);
+ }
+ const char *name;
+ err = PetscObjectGetName((PetscObject) vec, &name);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(label, std::string(name));
} // testNewSectionField
// ----------------------------------------------------------------------
@@ -308,9 +328,9 @@
void
pylith::topology::TestFieldMesh::testCloneSection(void)
{ // testCloneSection
- const int fiberDim = 3;
- const int nconstraints[] = { 0, 2, 1, 3 };
- const int constraints[] = {
+ const PetscInt fiberDim = 3;
+ const PetscInt nconstraints[] = { 0, 2, 1, 3 };
+ const PetscInt constraints[] = {
// 0
0, 2, // 1
2, // 2
@@ -319,32 +339,30 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
// Create field with atlas to use to create new field
Field<Mesh> fieldSrc(mesh);
{ // Setup source field
fieldSrc.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
- const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = fieldSrc.petscSection();
+ CPPUNIT_ASSERT(section);
int iV=0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- section->addConstraintDimension(*v_iter, nconstraints[iV++]);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionAddConstraintDof(section, v, nconstraints[iV++]);CHECK_PETSC_ERROR(err);
+ }
fieldSrc.allocate();
int index = 0;
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, index += nconstraints[i++])
- section->setConstraintDof(*v_iter, &constraints[index]);
+ iV = 0;
+ for(PetscInt v = vStart; v < vEnd; ++v, index += nconstraints[iV++]) {
+ err = PetscSectionSetConstraintIndices(section, v, (PetscInt *) &constraints[index]);CHECK_PETSC_ERROR(err);
+ }
fieldSrc.zero();
fieldSrc.createScatter(mesh);
fieldSrc.createScatter(mesh, "A");
@@ -354,24 +372,24 @@
const std::string& label = "field A";
field.label(label.c_str());
field.cloneSection(fieldSrc);
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
int iV = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
- CPPUNIT_ASSERT_EQUAL(nconstraints[iV++],
- section->getConstraintDimension(*v_iter));
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, cdof;
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(section, v, &cdof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ CPPUNIT_ASSERT_EQUAL(nconstraints[iV++], cdof);
} // for
// Verify vector scatters were also copied.
- CPPUNIT_ASSERT_EQUAL(fieldSrc._scatters[""].scatter,
- field._scatters[""].scatter);
- CPPUNIT_ASSERT_EQUAL(fieldSrc._scatters["A"].scatter,
- field._scatters["A"].scatter);
-
- CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+ CPPUNIT_ASSERT_EQUAL(fieldSrc._scatters[""].dm, field._scatters[""].dm);
+ CPPUNIT_ASSERT_EQUAL(fieldSrc._scatters["A"].dm, field._scatters["A"].dm);
+ const char *name;
+ err = PetscObjectGetName((PetscObject) vec, &name);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(label, std::string(name));
} // testCloneSection
// ----------------------------------------------------------------------
@@ -409,38 +427,44 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
const PylithScalar tolerance = 1.0e-6;
i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testAllocate
// ----------------------------------------------------------------------
@@ -459,39 +483,46 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.zero();
const PylithScalar tolerance = 1.0e-6;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[iDim], tolerance);
+ i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testZero
// ----------------------------------------------------------------------
@@ -499,7 +530,7 @@
void
pylith::topology::TestFieldMesh::testZeroAll(void)
{ // testZeroAll
- const int fiberDim = 3;
+ const PetscInt fiberDim = 3;
const PylithScalar scale = 2.0;
const PylithScalar valuesNondim[] = {
1.1, 2.2, 3.3,
@@ -507,8 +538,8 @@
1.3, 2.4, 3.5,
1.4, 2.5, 3.6,
};
- const int nconstraints[] = { 0, 2, 1, 3 };
- const int constraints[] = {
+ const PetscInt nconstraints[] = { 0, 2, 1, 3 };
+ const PetscInt constraints[] = {
// 0
0, 2, // 1
2, // 2
@@ -517,52 +548,58 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Create field and set constraint sizes
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- int iV=0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- section->addConstraintDimension(*v_iter, nconstraints[iV++]);
+ PetscSection section = field.petscSection();
+ PetscInt iV = 0;
+ PetscInt index = 0;
+ CPPUNIT_ASSERT(section);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionAddConstraintDof(section, v, nconstraints[iV++]);CHECK_PETSC_ERROR(err);
+ }
field.allocate();
- int index = 0;
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, index += nconstraints[i++])
- section->setConstraintDof(*v_iter, &constraints[index]);
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(vec);
+ iV = 0;
+ for(PetscInt v = vStart; v < vEnd; ++v, index += nconstraints[iV++]) {
+ err = PetscSectionSetConstraintIndices(section, v, (PetscInt *) &constraints[index]);CHECK_PETSC_ERROR(err);
+ }
field.zero();
- scalar_array values(fiberDim);
- i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePointAll(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.zeroAll();
const PylithScalar tolerance = 1.0e-6;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[iDim], tolerance);
+ i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testZeroAll
// ----------------------------------------------------------------------
@@ -581,41 +618,47 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.complete();
// Expect no change for this serial test
+ const PylithScalar tolerance = 1.0e-6;
i = 0;
- const PylithScalar tolerance = 1.0e-6;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testComplete
// ----------------------------------------------------------------------
@@ -634,49 +677,56 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> fieldSrc(mesh);
{ // Setup source field
fieldSrc.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
fieldSrc.allocate();
- const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = fieldSrc.petscSection();
+ Vec vec = fieldSrc.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // Setup source field
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
field.copy(fieldSrc);
- int i = 0;
- scalar_array values(fiberDim);
const PylithScalar tolerance = 1.0e-6;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testCopy
// ----------------------------------------------------------------------
@@ -701,63 +751,71 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> fieldSrc(mesh);
{ // Setup source field
fieldSrc.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
fieldSrc.allocate();
- const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = fieldSrc.petscSection();
+ Vec vec = fieldSrc.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesA[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesA[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // Setup source field
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
{ // Setup destination field
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesB[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesB[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // Setup destination field
field += fieldSrc;
- int i = 0;
- scalar_array values(fiberDim);
const PylithScalar tolerance = 1.0e-6;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar valueE = valuesA[i] + valuesB[i];
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesA[i] + valuesB[i], array[off+d], tolerance);
++i;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[iDim], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testOperateAdd
// ----------------------------------------------------------------------
@@ -776,43 +834,48 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.scale(scale);
field.addDimensionOkay(true);
field.dimensionalize();
+ const PylithScalar tolerance = 1.0e-6;
i = 0;
- const PylithScalar tolerance = 1.0e-6;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar valueE = valuesNondim[i++]*scale;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[iDim], tolerance);
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++]*scale, array[off+d], tolerance);
} // for
} // for
-
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testDimensionalize
// ----------------------------------------------------------------------
@@ -831,26 +894,31 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.view("Testing view");
} // testView
@@ -864,25 +932,26 @@
Mesh mesh;
_buildMesh(&mesh);
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
const std::string& label = "field A";
field.label(label.c_str());
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const int sizeE = vertices->size() * fiberDim;
+ const int sizeE = (vEnd-vStart) * fiberDim;
CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
field.createScatter(mesh);
CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
const Field<Mesh>::ScatterInfo& sinfo = field._getScatter("");
- CPPUNIT_ASSERT(sinfo.scatter);
- CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.dm);
CPPUNIT_ASSERT(sinfo.vector);
int size = 0;
@@ -902,8 +971,7 @@
field.createScatter(mesh, "B");
CPPUNIT_ASSERT_EQUAL(size_t(2), field._scatters.size());
const Field<Mesh>::ScatterInfo& sinfoB = field._getScatter("B");
- CPPUNIT_ASSERT(sinfoB.scatter);
- CPPUNIT_ASSERT(sinfoB.scatterVec);
+ CPPUNIT_ASSERT(sinfoB.dm);
CPPUNIT_ASSERT(sinfoB.vector);
Field<Mesh> field2(mesh);
@@ -911,13 +979,11 @@
CPPUNIT_ASSERT_EQUAL(size_t(2), field2._scatters.size());
const Field<Mesh>::ScatterInfo& sinfo2 = field2._getScatter("");
- CPPUNIT_ASSERT(sinfo2.scatter);
- CPPUNIT_ASSERT(sinfo2.scatterVec);
+ CPPUNIT_ASSERT(sinfo2.dm);
CPPUNIT_ASSERT(sinfo2.vector);
const Field<Mesh>::ScatterInfo& sinfo2B = field2._getScatter("B");
- CPPUNIT_ASSERT(sinfo2B.scatter);
- CPPUNIT_ASSERT(sinfo2B.scatterVec);
+ CPPUNIT_ASSERT(sinfo2B.dm);
CPPUNIT_ASSERT(sinfo2B.vector);
} // testCreateScatter
@@ -931,25 +997,26 @@
Mesh mesh;
_buildMesh(&mesh);
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
const std::string& label = "field A";
field.label(label.c_str());
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const int sizeE = vertices->size() * fiberDim;
+ const int sizeE = (vEnd-vStart) * fiberDim;
CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
field.createScatterWithBC(mesh);
CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
const Field<Mesh>::ScatterInfo& sinfo = field._getScatter("");
- CPPUNIT_ASSERT(sinfo.scatter);
- CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.dm);
CPPUNIT_ASSERT(sinfo.vector);
int size = 0;
@@ -969,8 +1036,7 @@
field.createScatterWithBC(mesh, "B");
CPPUNIT_ASSERT_EQUAL(size_t(2), field._scatters.size());
const Field<Mesh>::ScatterInfo& sinfoB = field._getScatter("B");
- CPPUNIT_ASSERT(sinfoB.scatter);
- CPPUNIT_ASSERT(sinfoB.scatterVec);
+ CPPUNIT_ASSERT(sinfoB.dm);
CPPUNIT_ASSERT(sinfoB.vector);
Field<Mesh> field2(mesh);
@@ -978,13 +1044,11 @@
CPPUNIT_ASSERT_EQUAL(size_t(2), field2._scatters.size());
const Field<Mesh>::ScatterInfo& sinfo2 = field2._getScatter("");
- CPPUNIT_ASSERT(sinfo2.scatter);
- CPPUNIT_ASSERT(sinfo2.scatterVec);
+ CPPUNIT_ASSERT(sinfo2.dm);
CPPUNIT_ASSERT(sinfo2.vector);
const Field<Mesh>::ScatterInfo& sinfo2B = field2._getScatter("B");
- CPPUNIT_ASSERT(sinfo2B.scatter);
- CPPUNIT_ASSERT(sinfo2B.scatterVec);
+ CPPUNIT_ASSERT(sinfo2B.dm);
CPPUNIT_ASSERT(sinfo2B.vector);
} // testCreateScatterWithBC
@@ -998,6 +1062,13 @@
Mesh mesh;
_buildMesh(&mesh);
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
@@ -1006,21 +1077,14 @@
field.createScatter(mesh);
CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
const Field<Mesh>::ScatterInfo& sinfo = field._getScatter("");
- CPPUNIT_ASSERT(sinfo.scatter);
- CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.dm);
CPPUNIT_ASSERT(sinfo.vector);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
-
const PetscVec vec = field.vector();
CPPUNIT_ASSERT_EQUAL(sinfo.vector, vec);
int size = 0;
- VecGetSize(vec, &size);
- const int sizeE = vertices->size() * fiberDim;
+ err = VecGetSize(vec, &size);CHECK_PETSC_ERROR(err);
+ const int sizeE = (vEnd-vStart) * fiberDim;
CPPUNIT_ASSERT_EQUAL(sizeE, size);
} // testVector
@@ -1040,42 +1104,47 @@
Mesh mesh;
_buildMesh(&mesh);
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscSection section = field.petscSection();
+ Vec secvec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(secvec);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesE[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(secvec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesE[i++];
} // for
+ err = VecRestoreArray(secvec, &array);CHECK_PETSC_ERROR(err);
field.createScatter(mesh, context);
field.scatterSectionToVector(context);
const PetscVec vec = field.vector(context);
CPPUNIT_ASSERT(0 != vec);
- int size = 0;
- VecGetSize(vec, &size);
+ PetscInt size = 0;
+ err = VecGetSize(vec, &size);CHECK_PETSC_ERROR(err);
PylithScalar* valuesVec = 0;
- VecGetArray(vec, &valuesVec);
+ err = VecGetArray(vec, &valuesVec);CHECK_PETSC_ERROR(err);
const PylithScalar tolerance = 1.0e-06;
- const int sizeE = vertices->size() * fiberDim;
+ const int sizeE = (vEnd-vStart) * fiberDim;
CPPUNIT_ASSERT_EQUAL(sizeE, size);
for (int i=0; i < sizeE; ++i)
CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i], valuesVec[i], tolerance);
- VecRestoreArray(vec, &valuesVec);
+ err = VecRestoreArray(vec, &valuesVec);CHECK_PETSC_ERROR(err);
} // testScatterSectionToVector
// ----------------------------------------------------------------------
@@ -1094,41 +1163,48 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<Mesh> field(mesh);
field.newSection(Field<Mesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<Mesh::RealSection>& section = field.section();
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
+ PetscSection section = field.petscSection();
+ CPPUNIT_ASSERT(section);
field.createScatter(mesh, context);
const PetscVec vec = field.vector(context);
CPPUNIT_ASSERT(0 != vec);
- int size = 0;
- VecGetSize(vec, &size);
+ PetscInt size = 0;
+ err = VecGetSize(vec, &size);CHECK_PETSC_ERROR(err);
PylithScalar* valuesVec = 0;
- VecGetArray(vec, &valuesVec);
+ err = VecGetArray(vec, &valuesVec);CHECK_PETSC_ERROR(err);
const PylithScalar tolerance = 1.0e-06;
- const int sizeE = vertices->size() * fiberDim;
+ const int sizeE = (vEnd-vStart) * fiberDim;
CPPUNIT_ASSERT_EQUAL(sizeE, size);
for (int i=0; i < sizeE; ++i)
valuesVec[i] = valuesE[i];
- VecRestoreArray(vec, &valuesVec);
+ err = VecRestoreArray(vec, &valuesVec);CHECK_PETSC_ERROR(err);
field.createScatter(mesh, context);
field.scatterVectorToSection(context);
- scalar_array values(fiberDim);
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], fiberDim);
- for (int iDim=0; iDim < fiberDim; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i++], values[iDim], tolerance);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i++], array[off+d], tolerance);
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testScatterVectorToSection
// ----------------------------------------------------------------------
@@ -1148,15 +1224,16 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
// Create field with atlas to use to create new field
Field<Mesh> fieldSrc(mesh);
+ CPPUNIT_ASSERT(0);
{ // Setup source field
fieldSrc.newSection(Field<Mesh>::VERTICES_FIELD, spaceDim);
fieldSrc.splitDefault();
@@ -1164,14 +1241,12 @@
CPPUNIT_ASSERT(!section.isNull());
int iV=0;
int iC=0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, ++iV) {
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iV) {
const int nconstraintsVertex = nconstraints[iV];
- section->addConstraintDimension(*v_iter, nconstraintsVertex);
+ section->addConstraintDimension(v, nconstraintsVertex);
for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint) {
const int fibration = constraints[iC++];
- section->addConstraintDimension(*v_iter, 1, fibration);
+ section->addConstraintDimension(v, 1, fibration);
} // for
} // for
fieldSrc.allocate();
@@ -1179,15 +1254,13 @@
iC = 0;
iV = 0;
int zero = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, ++iV) {
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iV) {
const int nconstraintsVertex = nconstraints[iV];
if (nconstraintsVertex > 0)
- section->setConstraintDof(*v_iter, &constraints[iC]);
+ section->setConstraintDof(v, &constraints[iC]);
for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint) {
const int fibration = constraints[iC++];
- section->setConstraintDof(*v_iter, &zero, fibration);
+ section->setConstraintDof(v, &zero, fibration);
} // for
} // for
} // Setup source field
@@ -1199,13 +1272,10 @@
for (int fibration=0; fibration < spaceDim; ++fibration) {
const ALE::Obj<Mesh::RealSection>& sectionSplit = section->getFibration(fibration);
CPPUNIT_ASSERT(!sectionSplit.isNull());
- CPPUNIT_ASSERT(!vertices.isNull());
int iV = 0;
int iC = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, ++iV) {
- CPPUNIT_ASSERT_EQUAL(1, section->getFiberDimension(*v_iter, fibration));
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iV) {
+ CPPUNIT_ASSERT_EQUAL(1, section->getFiberDimension(v, fibration));
bool isConstrained = false;
const int nconstraintsVertex = nconstraints[iV];
for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint)
@@ -1213,7 +1283,7 @@
isConstrained = true;
const int constraintDimE = (!isConstrained) ? 0 : 1;
CPPUNIT_ASSERT_EQUAL(constraintDimE,
- section->getConstraintDimension(*v_iter, fibration));
+ section->getConstraintDimension(v, fibration));
} // for
} // for
} // testSplitDefault
@@ -1235,15 +1305,16 @@
Mesh mesh;
_buildMesh(&mesh);
- const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
- const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
// Create field with atlas to use to create new field
Field<Mesh> fieldSrc(mesh);
+ CPPUNIT_ASSERT(0);
{ // Setup source field
fieldSrc.newSection(Field<Mesh>::VERTICES_FIELD, spaceDim);
fieldSrc.splitDefault();
@@ -1251,14 +1322,12 @@
CPPUNIT_ASSERT(!section.isNull());
int iV=0;
int iC=0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, ++iV) {
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iV) {
const int nconstraintsVertex = nconstraints[iV];
- section->addConstraintDimension(*v_iter, nconstraintsVertex);
+ section->addConstraintDimension(v, nconstraintsVertex);
for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint) {
const int fibration = constraints[iC++];
- section->addConstraintDimension(*v_iter, 1, fibration);
+ section->addConstraintDimension(v, 1, fibration);
} // for
} // for
fieldSrc.allocate();
@@ -1266,15 +1335,13 @@
iC = 0;
iV = 0;
int zero = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, ++iV) {
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iV) {
const int nconstraintsVertex = nconstraints[iV];
if (nconstraintsVertex > 0)
- section->setConstraintDof(*v_iter, &constraints[iC]);
+ section->setConstraintDof(v, &constraints[iC]);
for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint) {
const int fibration = constraints[iC++];
- section->setConstraintDof(*v_iter, &zero, fibration);
+ section->setConstraintDof(v, &zero, fibration);
} // for
} // for
} // Setup source field
@@ -1289,13 +1356,10 @@
for (int fibration=0; fibration < spaceDim; ++fibration) {
const ALE::Obj<Mesh::RealSection>& sectionSplit = section->getFibration(fibration);
CPPUNIT_ASSERT(!sectionSplit.isNull());
- CPPUNIT_ASSERT(!vertices.isNull());
int iV = 0;
int iC = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, ++iV) {
- CPPUNIT_ASSERT_EQUAL(1, section->getFiberDimension(*v_iter, fibration));
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iV) {
+ CPPUNIT_ASSERT_EQUAL(1, section->getFiberDimension(v, fibration));
bool isConstrained = false;
const int nconstraintsVertex = nconstraints[iV];
for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint)
@@ -1303,7 +1367,7 @@
isConstrained = true;
const int constraintDimE = (!isConstrained) ? 0 : 1;
CPPUNIT_ASSERT_EQUAL(constraintDimE,
- section->getConstraintDimension(*v_iter, fibration));
+ section->getConstraintDimension(v, fibration));
} // for
} // for
} // testCloneSectionSplit
@@ -1323,6 +1387,10 @@
ALE::Obj<SieveFlexMesh::sieve_type> s =
new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
+
+ mesh->createDMMesh(_TestFieldMesh::cellDim);
+ DM dmMesh = mesh->dmMesh();
+ PetscErrorCode err;
const int cellDim = _TestFieldMesh::cellDim;
const int ncells = _TestFieldMesh::ncells;
@@ -1342,6 +1410,48 @@
ALE::SieveBuilder<Mesh::SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
coordinates);
+ err = DMComplexSetChart(dmMesh, 0, ncells+nvertices);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = 0; c < ncells; ++c) {
+ err = DMComplexSetConeSize(dmMesh, c, ncorners);CHECK_PETSC_ERROR(err);
+ }
+ err = DMSetUp(dmMesh);CHECK_PETSC_ERROR(err);
+ PetscInt *cone = new PetscInt[ncorners];
+ for(PetscInt c = 0; c < ncells; ++c) {
+ for(PetscInt v = 0; v < ncorners; ++v) {
+ cone[v] = cells[c*ncorners+v]+ncells;
+ }
+ err = DMComplexSetCone(dmMesh, c, cone);CHECK_PETSC_ERROR(err);
+ } // for
+ delete [] cone; cone = 0;
+ err = DMComplexSymmetrize(dmMesh);CHECK_PETSC_ERROR(err);
+ err = DMComplexStratify(dmMesh);CHECK_PETSC_ERROR(err);
+ PetscSection coordSection;
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, ncells, ncells+nvertices);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = ncells; v < ncells+nvertices; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(sieveMesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = 0; v < nvertices; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v+ncells, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = coordinates[v*spaceDim+d];
+ }
+ }
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(dmMesh, coordVec);CHECK_PETSC_ERROR(err);
+
spatialdata::geocoords::CSCart cs;
cs.setSpaceDim(spaceDim);
cs.initialize();
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldSubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldSubMesh.cc 2012-10-09 23:53:54 UTC (rev 20815)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestFieldSubMesh.cc 2012-10-10 02:15:08 UTC (rev 20816)
@@ -84,8 +84,8 @@
_buildMesh(&mesh, &submesh);
Field<SubMesh> field(submesh);
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(section.isNull());
+ PetscSection section = field.petscSection();
+ CPPUNIT_ASSERT(section);
} // testSection
// ----------------------------------------------------------------------
@@ -126,8 +126,8 @@
Field<SubMesh> field(submesh);
field.newSection();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ CPPUNIT_ASSERT(section);
} // testNewSection
// ----------------------------------------------------------------------
@@ -140,21 +140,23 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
Field<SubMesh> field(submesh);
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- field.newSection(vertices, fiberDim);
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ field.newSection(topology::FieldBase::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ PetscSection section = field.petscSection();
+ CPPUNIT_ASSERT(section);
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof;
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ }
} // testNewSectionPoints
// ----------------------------------------------------------------------
@@ -167,21 +169,24 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
Field<SubMesh> field(submesh);
field.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof;
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ }
} // testNewSectionDomain
// ----------------------------------------------------------------------
@@ -194,27 +199,29 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
// Create field with atlas to use to create new field
Field<SubMesh> fieldSrc(submesh);
fieldSrc.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
- const ALE::Obj<SubMesh::RealSection>& sectionSrc = fieldSrc.section();
- CPPUNIT_ASSERT(!sectionSrc.isNull());
+ fieldSrc.allocate();
const int fiberDim2 = 4;
Field<SubMesh> field(submesh);
field.newSection(fieldSrc, fiberDim2);
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- CPPUNIT_ASSERT_EQUAL(fiberDim2, section->getFiberDimension(*v_iter));
+ field.allocate();
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof;
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim2, dof);
+ }
} // testNewSectionChart
// ----------------------------------------------------------------------
@@ -234,49 +241,46 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
// Create field with atlas to use to create new field
Field<SubMesh> fieldSrc(submesh);
{ // Setup source field
fieldSrc.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
- const ALE::Obj<SubMesh::RealSection>& section = fieldSrc.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = fieldSrc.petscSection();
+ CPPUNIT_ASSERT(section);
int iV=0;
-
- CPPUNIT_ASSERT(!vertices.isNull());
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter)
- section->addConstraintDimension(*v_iter, nconstraints[iV++]);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionAddConstraintDof(section, v, nconstraints[iV++]);CHECK_PETSC_ERROR(err);
+ }
fieldSrc.allocate();
int index = 0;
- int i = 0;
- for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter, index += nconstraints[i++])
- section->setConstraintDof(*v_iter, &constraints[index]);
+ iV = 0;
+ for(PetscInt v = vStart; v < vEnd; ++v, index += nconstraints[iV++]) {
+ err = PetscSectionSetConstraintIndices(section, v, (PetscInt *) &constraints[index]);CHECK_PETSC_ERROR(err);
+ }
fieldSrc.zero();
} // Setup source field
Field<SubMesh> field(submesh);
field.cloneSection(fieldSrc);
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
int iV = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
- CPPUNIT_ASSERT_EQUAL(nconstraints[iV++],
- section->getConstraintDimension(*v_iter));
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, cdof;
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(section, v, &cdof);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, dof);
+ CPPUNIT_ASSERT_EQUAL(nconstraints[iV++], cdof);
} // for
} // testCloneSection
@@ -317,38 +321,44 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
const PylithScalar tolerance = 1.0e-6;
i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testAllocate
// ----------------------------------------------------------------------
@@ -367,39 +377,46 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.zero();
const PylithScalar tolerance = 1.0e-6;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[iDim], tolerance);
+ i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testZero
// ----------------------------------------------------------------------
@@ -418,41 +435,47 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.complete();
// Expect no change for this serial test
+ const PylithScalar tolerance = 1.0e-6;
i = 0;
- const PylithScalar tolerance = 1.0e-6;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testComplete
// ----------------------------------------------------------------------
@@ -471,49 +494,56 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> fieldSrc(submesh);
{ // Setup source field
fieldSrc.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
fieldSrc.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = fieldSrc.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = fieldSrc.petscSection();
+ Vec vec = fieldSrc.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // Setup source field
Field<SubMesh> field(submesh);
field.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
field.copy(fieldSrc);
- int i = 0;
- scalar_array values(fiberDim);
const PylithScalar tolerance = 1.0e-6;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testCopy
// ----------------------------------------------------------------------
@@ -537,63 +567,71 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> fieldSrc(submesh);
{ // Setup source field
fieldSrc.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
fieldSrc.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = fieldSrc.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = fieldSrc.petscSection();
+ Vec vec = fieldSrc.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesA[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesA[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // Setup source field
Field<SubMesh> field(submesh);
field.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
{ // Setup destination field
- scalar_array values(fiberDim);
- int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesB[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesB[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // Setup destination field
field += fieldSrc;
- int i = 0;
- scalar_array values(fiberDim);
const PylithScalar tolerance = 1.0e-6;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar valueE = valuesA[i] + valuesB[i];
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesA[i] + valuesB[i], array[off+d], tolerance);
++i;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[iDim], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testOperateAdd
// ----------------------------------------------------------------------
@@ -612,44 +650,48 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.scale(scale);
field.addDimensionOkay(true);
field.dimensionalize();
+ const PylithScalar tolerance = 1.0e-6;
i = 0;
- const PylithScalar tolerance = 1.0e-6;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], values.size());
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar valueE = valuesNondim[i++]*scale;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[iDim], tolerance);
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++]*scale, array[off+d], tolerance);
} // for
} // for
-
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testDimensionalize
// ----------------------------------------------------------------------
@@ -668,26 +710,31 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
- const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(Field<SubMesh>::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- scalar_array values(fiberDim);
- int i = 0;
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesNondim[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesNondim[i++];
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
field.view("Testing view");
} // testView
@@ -702,23 +749,24 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const int sizeE = vertices->size() * fiberDim;
-
+ const int sizeE = (vEnd-vStart) * fiberDim;
+
CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
field.createScatter(submesh);
CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
const Field<SubMesh>::ScatterInfo& sinfo = field._getScatter("");
- CPPUNIT_ASSERT(sinfo.scatter);
- CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.dm);
CPPUNIT_ASSERT(sinfo.vector);
int size = 0;
@@ -733,8 +781,7 @@
field.createScatter(submesh, "B");
CPPUNIT_ASSERT_EQUAL(size_t(2), field._scatters.size());
const Field<SubMesh>::ScatterInfo& sinfoB = field._getScatter("B");
- CPPUNIT_ASSERT(sinfoB.scatter);
- CPPUNIT_ASSERT(sinfoB.scatterVec);
+ CPPUNIT_ASSERT(sinfoB.dm);
CPPUNIT_ASSERT(sinfoB.vector);
Field<SubMesh> field2(submesh);
@@ -742,13 +789,11 @@
CPPUNIT_ASSERT_EQUAL(size_t(2), field2._scatters.size());
const Field<SubMesh>::ScatterInfo& sinfo2 = field2._getScatter("");
- CPPUNIT_ASSERT(sinfo2.scatter);
- CPPUNIT_ASSERT(sinfo2.scatterVec);
+ CPPUNIT_ASSERT(sinfo2.dm);
CPPUNIT_ASSERT(sinfo2.vector);
const Field<SubMesh>::ScatterInfo& sinfo2B = field2._getScatter("B");
- CPPUNIT_ASSERT(sinfo2B.scatter);
- CPPUNIT_ASSERT(sinfo2B.scatterVec);
+ CPPUNIT_ASSERT(sinfo2B.dm);
CPPUNIT_ASSERT(sinfo2B.vector);
} // testCreateScatter
@@ -762,23 +807,24 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const int sizeE = vertices->size() * fiberDim;
+ const int sizeE = (vEnd-vStart) * fiberDim;
CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
field.createScatterWithBC(submesh);
CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
const Field<SubMesh>::ScatterInfo& sinfo = field._getScatter("");
- CPPUNIT_ASSERT(sinfo.scatter);
- CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.dm);
CPPUNIT_ASSERT(sinfo.vector);
int size = 0;
@@ -793,8 +839,7 @@
field.createScatterWithBC(submesh, "B");
CPPUNIT_ASSERT_EQUAL(size_t(2), field._scatters.size());
const Field<SubMesh>::ScatterInfo& sinfoB = field._getScatter("B");
- CPPUNIT_ASSERT(sinfoB.scatter);
- CPPUNIT_ASSERT(sinfoB.scatterVec);
+ CPPUNIT_ASSERT(sinfoB.dm);
CPPUNIT_ASSERT(sinfoB.vector);
Field<SubMesh> field2(submesh);
@@ -802,13 +847,11 @@
CPPUNIT_ASSERT_EQUAL(size_t(2), field2._scatters.size());
const Field<SubMesh>::ScatterInfo& sinfo2 = field2._getScatter("");
- CPPUNIT_ASSERT(sinfo2.scatter);
- CPPUNIT_ASSERT(sinfo2.scatterVec);
+ CPPUNIT_ASSERT(sinfo2.dm);
CPPUNIT_ASSERT(sinfo2.vector);
const Field<SubMesh>::ScatterInfo& sinfo2B = field2._getScatter("B");
- CPPUNIT_ASSERT(sinfo2B.scatter);
- CPPUNIT_ASSERT(sinfo2B.scatterVec);
+ CPPUNIT_ASSERT(sinfo2B.dm);
CPPUNIT_ASSERT(sinfo2B.vector);
} // testCreateScatterWithBC
@@ -822,7 +865,13 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
field.allocate();
@@ -831,21 +880,14 @@
field.createScatter(submesh);
CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
const Field<SubMesh>::ScatterInfo& sinfo = field._getScatter("");
- CPPUNIT_ASSERT(sinfo.scatter);
- CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.dm);
CPPUNIT_ASSERT(sinfo.vector);
- const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
-
const PetscVec vec = field.vector();
CPPUNIT_ASSERT_EQUAL(sinfo.vector, vec);
int size = 0;
- VecGetSize(vec, &size);
- const int sizeE = vertices->size() * fiberDim;
+ err = VecGetSize(vec, &size);CHECK_PETSC_ERROR(err);
+ const int sizeE = (vEnd-vStart) * fiberDim;
CPPUNIT_ASSERT_EQUAL(sizeE, size);
} // testVector
@@ -865,39 +907,47 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
- const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
field.allocate();
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
+ PetscSection section = field.petscSection();
+ Vec secvec = field.localVector();
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(secvec);
- scalar_array values(fiberDim);
- int i = 0;
- for (SubMesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- for (int iDim=0; iDim < fiberDim; ++iDim)
- values[iDim] = valuesE[i++];
- section->updatePoint(*v_iter, &values[0]);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(secvec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ array[off+d] = valuesE[i++];
} // for
+ err = VecRestoreArray(secvec, &array);CHECK_PETSC_ERROR(err);
field.createScatter(submesh, context);
field.scatterSectionToVector(context);
const PetscVec vec = field.vector(context);
CPPUNIT_ASSERT(0 != vec);
int size = 0;
- VecGetSize(vec, &size);
+ err = VecGetSize(vec, &size);CHECK_PETSC_ERROR(err);
PylithScalar* valuesVec = 0;
- VecGetArray(vec, &valuesVec);
+ err = VecGetArray(vec, &valuesVec);CHECK_PETSC_ERROR(err);
const PylithScalar tolerance = 1.0e-06;
- const int sizeE = vertices->size() * fiberDim;
+ const int sizeE = (vEnd-vStart) * fiberDim;
CPPUNIT_ASSERT_EQUAL(sizeE, size);
for (int i=0; i < sizeE; ++i)
CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i], valuesVec[i], tolerance);
- VecRestoreArray(vec, &valuesVec);
+ err = VecRestoreArray(vec, &valuesVec);CHECK_PETSC_ERROR(err);
} // testScatterSectionToVector
// ----------------------------------------------------------------------
@@ -916,44 +966,47 @@
Mesh mesh;
SubMesh submesh;
_buildMesh(&mesh, &submesh);
+ DM dmMesh = submesh.dmMesh();
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+
+ PetscInt vStart, vEnd;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
Field<SubMesh> field(submesh);
field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
field.allocate();
+ PetscSection section = field.petscSection();
+ CPPUNIT_ASSERT(section);
field.createScatter(submesh, context);
- const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
-
const PetscVec vec = field.vector(context);
CPPUNIT_ASSERT(0 != vec);
int size = 0;
- VecGetSize(vec, &size);
- const int sizeE = vertices->size() * fiberDim;
- CPPUNIT_ASSERT_EQUAL(sizeE, size);
+ err = VecGetSize(vec, &size);CHECK_PETSC_ERROR(err);
+ PylithScalar* valuesVec = 0;
+ err = VecGetArray(vec, &valuesVec);CHECK_PETSC_ERROR(err);
const PylithScalar tolerance = 1.0e-06;
- PylithScalar* valuesVec = 0;
- VecGetArray(vec, &valuesVec);
+ const int sizeE = (vEnd-vStart) * fiberDim;
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
for (int i=0; i < sizeE; ++i)
valuesVec[i] = valuesE[i];
- VecRestoreArray(vec, &valuesVec);
+ err = VecRestoreArray(vec, &valuesVec);CHECK_PETSC_ERROR(err);
field.scatterVectorToSection(context);
- scalar_array values(fiberDim);
- int i = 0;
- const ALE::Obj<SubMesh::RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
- for (SubMesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != vertices->end();
- ++v_iter) {
- section->restrictPoint(*v_iter, &values[0], fiberDim);
- for (int iDim=0; iDim < fiberDim; ++iDim)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i++], values[iDim], tolerance);
+ PetscScalar *array;
+ PetscInt i = 0;
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < fiberDim; ++d)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i++], array[off+d], tolerance);
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testScatterVectorToSection
// ----------------------------------------------------------------------
@@ -975,6 +1028,10 @@
ALE::Obj<SieveFlexMesh::sieve_type> s =
new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
CPPUNIT_ASSERT(!s.isNull());
+
+ mesh->createDMMesh(_TestFieldSubMesh::cellDim);
+ DM dmMesh = mesh->dmMesh();
+ PetscErrorCode err;
const int cellDim = _TestFieldSubMesh::cellDim;
const int ncells = _TestFieldSubMesh::ncells;
@@ -993,7 +1050,49 @@
sieveMesh->stratify();
ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
coordinates);
+
+ err = DMComplexSetChart(dmMesh, 0, ncells+nvertices);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = 0; c < ncells; ++c) {
+ err = DMComplexSetConeSize(dmMesh, c, ncorners);CHECK_PETSC_ERROR(err);
+ }
+ err = DMSetUp(dmMesh);CHECK_PETSC_ERROR(err);
+ PetscInt *cone = new PetscInt[ncorners];
+ for(PetscInt c = 0; c < ncells; ++c) {
+ for(PetscInt v = 0; v < ncorners; ++v) {
+ cone[v] = cells[c*ncorners+v]+ncells;
+ }
+ err = DMComplexSetCone(dmMesh, c, cone);CHECK_PETSC_ERROR(err);
+ } // for
+ delete [] cone; cone = 0;
+ err = DMComplexSymmetrize(dmMesh);CHECK_PETSC_ERROR(err);
+ err = DMComplexStratify(dmMesh);CHECK_PETSC_ERROR(err);
+ PetscSection coordSection;
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, ncells, ncells+nvertices);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = ncells; v < ncells+nvertices; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(sieveMesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = 0; v < nvertices; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v+ncells, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = coordinates[v*spaceDim+d];
+ }
+ }
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(dmMesh, coordVec);CHECK_PETSC_ERROR(err);
+
typedef Mesh::SieveMesh::int_section_type::chart_type chart_type;
const ALE::Obj<SieveMesh::int_section_type>& groupField =
sieveMesh->getIntSection(_TestFieldSubMesh::label);
@@ -1008,6 +1107,10 @@
1);
sieveMesh->allocate(groupField);
+ for(PetscInt i = 0; i < numPoints; ++i) {
+ err = DMComplexSetLabelValue(dmMesh, _TestFieldSubMesh::label, ncells+_TestFieldSubMesh::groupVertices[i], 1);CHECK_PETSC_ERROR(err);
+ }
+
spatialdata::geocoords::CSCart cs;
cs.setSpaceDim(spaceDim);
cs.initialize();
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.cc 2012-10-09 23:53:54 UTC (rev 20815)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestSubMesh.cc 2012-10-10 02:15:08 UTC (rev 20816)
@@ -217,6 +217,10 @@
ALE::Obj<SieveFlexMesh::sieve_type> s =
new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
+
+ mesh->createDMMesh(_TestSubMesh::cellDim);
+ DM dmMesh = mesh->dmMesh();
+ PetscErrorCode err;
const int cellDim = _TestSubMesh::cellDim;
const int ncells = _TestSubMesh::ncells;
@@ -236,6 +240,48 @@
ALE::SieveBuilder<Mesh::SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
coordinates);
+ err = DMComplexSetChart(dmMesh, 0, ncells+nvertices);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = 0; c < ncells; ++c) {
+ err = DMComplexSetConeSize(dmMesh, c, ncorners);CHECK_PETSC_ERROR(err);
+ }
+ err = DMSetUp(dmMesh);CHECK_PETSC_ERROR(err);
+ PetscInt *cone = new PetscInt[ncorners];
+ for(PetscInt c = 0; c < ncells; ++c) {
+ for(PetscInt v = 0; v < ncorners; ++v) {
+ cone[v] = cells[c*ncorners+v]+ncells;
+ }
+ err = DMComplexSetCone(dmMesh, c, cone);CHECK_PETSC_ERROR(err);
+ } // for
+ delete [] cone; cone = 0;
+ err = DMComplexSymmetrize(dmMesh);CHECK_PETSC_ERROR(err);
+ err = DMComplexStratify(dmMesh);CHECK_PETSC_ERROR(err);
+ PetscSection coordSection;
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, ncells, ncells+nvertices);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = ncells; v < ncells+nvertices; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(sieveMesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = 0; v < nvertices; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v+ncells, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = coordinates[v*spaceDim+d];
+ }
+ }
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(dmMesh, coordVec);CHECK_PETSC_ERROR(err);
+
spatialdata::geocoords::CSCart cs;
cs.setSpaceDim(spaceDim);
cs.initialize();
@@ -253,7 +299,10 @@
for(int i=0; i < numPoints; ++i)
groupField->setFiberDimension(numCells+_TestSubMesh::groupVertices[i], 1);
sieveMesh->allocate(groupField);
-
+
+ for(PetscInt i = 0; i < numPoints; ++i) {
+ err = DMComplexSetLabelValue(dmMesh, _TestSubMesh::label, ncells+_TestSubMesh::groupVertices[i], 1);CHECK_PETSC_ERROR(err);
+ }
} // _buildMesh
More information about the CIG-COMMITS
mailing list