[cig-commits] r21009 - in short/3D/PyLith/trunk: libsrc/pylith/faults unittests/libtests/faults unittests/libtests/faults/data
knepley at geodynamics.org
knepley at geodynamics.org
Sat Nov 10 06:41:39 PST 2012
Author: knepley
Date: 2012-11-10 06:41:39 -0800 (Sat, 10 Nov 2012)
New Revision: 21009
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.hh
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
Log:
Now CohesiveKinLine2 tests work (except splitting)
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-11-10 14:41:39 UTC (rev 21009)
@@ -29,7 +29,6 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
#include "pylith/topology/Fields.hh" // USES Fields
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
#include "pylith/topology/Jacobian.hh" // USES Jacobian
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/friction/FrictionModel.hh" // USES FrictionModel
@@ -210,56 +209,49 @@
const int spaceDim = _quadrature->spaceDim();
// Get sections associated with cohesive cells
- scalar_array residualVertexN(spaceDim);
- scalar_array residualVertexP(spaceDim);
- scalar_array residualVertexL(spaceDim);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- assert(!residualSection.isNull());
+ DM residualDM = residual.dmMesh();
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscSection residualGlobalSection;
+ PetscScalar *residualArray;
+ PetscErrorCode err;
+ assert(residualSection);assert(residualVec);
+ err = DMGetDefaultGlobalSection(residualDM, &residualGlobalSection);CHECK_PETSC_ERROR(err);
- const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
+ PetscSection dispTSection = fields->get("disp(t)").petscSection();
+ Vec dispTVec = fields->get("disp(t)").localVector();
+ PetscScalar *dispTArray;
+ assert(dispTSection);assert(dispTVec);
- const ALE::Obj<RealSection>& dispTIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- assert(!dispTIncrSection.isNull());
+ PetscSection dispTIncrSection = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ PetscScalar *dispTIncrArray;
+ assert(dispTIncrSection);assert(dispTIncrVec);
- scalar_array dispTpdtVertexN(spaceDim);
- scalar_array dispTpdtVertexP(spaceDim);
- scalar_array dispTpdtVertexL(spaceDim);
-
scalar_array tractPerturbVertex(spaceDim);
- ALE::Obj<RealUniformSection> tractPerturbSection;
- int tractPerturbIndex = 0;
- int tractPerturbFiberDim = 0;
+ PetscSection valuesSection;
+ Vec valuesVec;
+ PetscScalar *tractionsArray;
if (_tractPerturbation) {
_tractPerturbation->calculate(t);
- const topology::FieldsNew<topology::SubMesh>* params = _tractPerturbation->parameterFields();
+ const topology::Fields<topology::Field<topology::SubMesh> >* params = _tractPerturbation->parameterFields();
assert(params);
- tractPerturbSection = params->section();
- assert(!tractPerturbSection.isNull());
-
- tractPerturbFiberDim = params->fiberDim();
- tractPerturbIndex = params->sectionIndex("value");
- const int valueFiberDim = params->sectionFiberDim("value");
- assert(valueFiberDim == spaceDim);
+ valuesSection = params->get("value").petscSection();
+ valuesVec = params->get("value").localVector();
+ assert(valuesSection);assert(valuesVec);
} // if
- const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
- assert(!areaSection.isNull());
+ PetscSection areaSection = _fields->get("area").petscSection();
+ Vec areaVec = _fields->get("area").localVector();
+ PetscScalar *areaArray;
+ assert(areaSection);assert(areaVec);
- const ALE::Obj<RealSection>& orientationSection =
- _fields->get("orientation").section();
- assert(!orientationSection.isNull());
+ PetscSection orientationSection = _fields->get("orientation").petscSection();
+ Vec orientationVec = _fields->get("orientation").localVector();
+ PetscScalar *orientationArray;
+ assert(orientationSection);assert(orientationVec);
- // Get fault information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- residualSection);
- assert(!globalOrder.isNull());
-
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(computeEvent);
@@ -267,15 +259,21 @@
// Loop over fault vertices
const int numVertices = _cohesiveVertices.size();
+ err = VecGetArray(areaVec, &areaArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
const int v_fault = _cohesiveVertices[iVertex].fault;
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+ PetscInt goff;
// Compute contribution only if Lagrange constraint is local.
- if (!globalOrder->isLocal(v_lagrange))
- continue;
+ err = PetscSectionGetOffset(residualGlobalSection, v_lagrange, &goff);CHECK_PETSC_ERROR(err);
+ if (goff < 0) continue;
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventBegin(restrictEvent);
@@ -283,122 +281,122 @@
// Get prescribed traction perturbation at fault vertex.
if (_tractPerturbation) {
- assert(tractPerturbFiberDim == tractPerturbSection->getFiberDimension(v_fault));
- const PylithScalar* paramsVertex = tractPerturbSection->restrictPoint(v_fault);
- assert(paramsVertex);
+ PetscInt vdof, voff;
+
+ err = PetscSectionGetDof(valuesSection, v_fault, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valuesSection, v_fault, &voff);CHECK_PETSC_ERROR(err);
+ assert(vdof == spaceDim);
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- tractPerturbVertex[iDim] = paramsVertex[tractPerturbIndex+iDim];
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ tractPerturbVertex[d] = tractionsArray[voff+d];
} // for
} else {
tractPerturbVertex = 0.0;
} // if/else
// Get orientation associated with fault vertex.
- assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
- const PylithScalar* orientationVertex =
- orientationSection->restrictPoint(v_fault);
- assert(orientationVertex);
+ PetscInt odof, ooff;
+ err = PetscSectionGetDof(orientationSection, v_fault, &odof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, v_fault, &ooff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim*spaceDim == odof);
+
// Get area associated with fault vertex.
- assert(1 == areaSection->getFiberDimension(v_fault));
- assert(areaSection->restrictPoint(v_fault));
- const PylithScalar areaVertex = *areaSection->restrictPoint(v_fault);
+ PetscInt adof, aoff;
+ err = PetscSectionGetDof(areaSection, v_fault, &adof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(areaSection, v_fault, &aoff);CHECK_PETSC_ERROR(err);
+ assert(1 == adof);
+
// Get disp(t) at conventional vertices and Lagrange vertex.
- assert(spaceDim == dispTSection->getFiberDimension(v_negative));
- const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
- assert(dispTVertexN);
+ PetscInt dtndof, dtnoff;
- assert(spaceDim == dispTSection->getFiberDimension(v_positive));
- const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
- assert(dispTVertexP);
+ err = PetscSectionGetDof(dispTSection, v_negative, &dtndof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTSection, v_negative, &dtnoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dtndof);
+ PetscInt dtpdof, dtpoff;
- assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
- const PylithScalar* dispTVertexL = dispTSection->restrictPoint(v_lagrange);
- assert(dispTVertexL);
+ err = PetscSectionGetDof(dispTSection, v_positive, &dtpdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTSection, v_positive, &dtpoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dtpdof);
+ PetscInt dtldof, dtloff;
+ err = PetscSectionGetDof(dispTSection, v_lagrange, &dtldof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTSection, v_lagrange, &dtloff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dtldof);
+
// Get dispIncr(t->t+dt) at conventional vertices and Lagrange vertex.
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_negative));
- const PylithScalar* dispTIncrVertexN =
- dispTIncrSection->restrictPoint(v_negative);
- assert(dispTIncrVertexN);
+ PetscInt dindof, dinoff;
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_positive));
- const PylithScalar* dispTIncrVertexP =
- dispTIncrSection->restrictPoint(v_positive);
- assert(dispTIncrVertexP);
+ err = PetscSectionGetDof(dispTIncrSection, v_negative, &dindof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTIncrSection, v_negative, &dinoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dindof);
+ PetscInt dipdof, dipoff;
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
- const PylithScalar* dispTIncrVertexL =
- dispTIncrSection->restrictPoint(v_lagrange);
- assert(dispTIncrVertexL);
+ err = PetscSectionGetDof(dispTIncrSection, v_positive, &dipdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTIncrSection, v_positive, &dipoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dipdof);
+ PetscInt dildof, diloff;
+ err = PetscSectionGetDof(dispTIncrSection, v_lagrange, &dildof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTIncrSection, v_lagrange, &diloff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dildof);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
#endif
- // Compute current estimate of displacement at time t+dt using
- // solution increment.
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- dispTpdtVertexN[iDim] = dispTVertexN[iDim] + dispTIncrVertexN[iDim];
- dispTpdtVertexP[iDim] = dispTVertexP[iDim] + dispTIncrVertexP[iDim];
- dispTpdtVertexL[iDim] = dispTVertexL[iDim] + dispTIncrVertexL[iDim];
- } // for
-
// Compute slip (in fault coordinates system) from displacements.
- PylithScalar slipNormal = 0.0;
- PylithScalar tractionNormal = 0.0;
- const int indexN = spaceDim - 1;
- for (int jDim=0; jDim < spaceDim; ++jDim) {
- slipNormal += orientationVertex[indexN*spaceDim+jDim] *
- (dispTpdtVertexP[jDim] - dispTpdtVertexN[jDim]);
- tractionNormal +=
- orientationVertex[indexN*spaceDim+jDim] * dispTpdtVertexL[jDim];
+ PylithScalar slipNormal = 0.0;
+ PylithScalar tractionNormal = 0.0;
+ const PetscInt indexN = spaceDim - 1;
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ slipNormal += orientationArray[ooff+indexN*spaceDim+d] * (dispTArray[dtpoff+d] + dispTIncrArray[dipoff+d] - dispTArray[dtnoff+d] - dispTIncrArray[dinoff+d]);
+ tractionNormal += orientationArray[ooff+indexN*spaceDim+d] * (dispTArray[dtloff+d] + dispTIncrArray[diloff+d]);
} // for
- residualVertexN = 0.0;
- residualVertexL = 0.0;
+
+#if defined(DETAILED_EVENT_LOGGING)
+ _logger->eventEnd(computeEvent);
+ _logger->eventBegin(updateEvent);
+#endif
if (slipNormal < _zeroTolerance || !_openFreeSurf) {
- // if no opening or flag indicates to still impose initial
- // tractions when fault is open.
- //
- // Initial (external) tractions oppose (internal) tractions
- // associated with Lagrange multiplier.
- residualVertexN = areaVertex * (dispTpdtVertexL - tractPerturbVertex);
+ // if no opening or flag indicates to still impose initial tractions when fault is open.
+ // Assemble contributions into field
+ PetscInt rndof, rnoff, rpdof, rpoff;
+ err = PetscSectionGetDof(residualSection, v_negative, &rndof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(residualSection, v_negative, &rnoff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(residualSection, v_positive, &rpdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(residualSection, v_positive, &rpoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == rndof);assert(spaceDim == rpdof);
+ // Initial (external) tractions oppose (internal) tractions associated with Lagrange multiplier.
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ residualArray[rnoff+d] += areaArray[aoff] * (dispTArray[dtloff+d] + dispTIncrArray[diloff+d] - tractPerturbVertex[d]);
+ residualArray[rpoff+d] += -areaArray[aoff] * (dispTArray[dtloff+d] + dispTIncrArray[diloff+d] - tractPerturbVertex[d]);
+ }
} else { // opening, normal traction should be zero
if (fabs(tractionNormal) > _zeroTolerance) {
- std::cerr << "WARNING! Fault opening with nonzero traction."
- << ", v_fault: " << v_fault
- << ", opening: " << slipNormal
- << ", normal traction: " << tractionNormal
- << std::endl;
+ std::cerr << "WARNING! Fault opening with nonzero traction."
+ << ", v_fault: " << v_fault
+ << ", opening: " << slipNormal
+ << ", normal traction: " << tractionNormal
+ << std::endl;
} // if
assert(fabs(tractionNormal) < _zeroTolerance);
} // if/else
- residualVertexP = -residualVertexN;
#if defined(DETAILED_EVENT_LOGGING)
- _logger->eventEnd(computeEvent);
- _logger->eventBegin(updateEvent);
-#endif
-
- // Assemble contributions into field
- assert(residualVertexN.size() ==
- residualSection->getFiberDimension(v_negative));
- residualSection->updateAddPoint(v_negative, &residualVertexN[0]);
-
- assert(residualVertexP.size() ==
- residualSection->getFiberDimension(v_positive));
- residualSection->updateAddPoint(v_positive, &residualVertexP[0]);
-
-#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
PetscLogFlops(numVertices*spaceDim*8);
+ err = VecRestoreArray(areaVec, &areaArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-11-10 14:41:39 UTC (rev 21009)
@@ -1687,6 +1687,7 @@
logger.stagePop();
+ scalar_array coordinatesCell(numBasis*spaceDim);
PetscSection coordSection;
Vec coordVec;
err = DMComplexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
@@ -1703,6 +1704,9 @@
const PetscScalar *coords = PETSC_NULL;
PetscInt coordsSize;
err = DMComplexVecGetClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ _quadrature->computeGeometry(coordinatesCell, c);
+ err = DMComplexVecRestoreClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
#endif
// Get cell geometry information that depends on cell
@@ -1718,11 +1722,7 @@
} // for
} // for
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) );
+ PetscLogFlops(numQuadPts*(1+numBasis*2));
} // for
// Assemble area information
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc 2012-11-10 14:41:39 UTC (rev 21009)
@@ -23,7 +23,7 @@
#include "FaultCohesiveLagrange.hh" // USES faultToGlobal()
#include "pylith/topology/SubMesh.hh" // USES SubMesh
-#include "pylith/topology/FieldsNew.hh" // HOLDSA FieldsNew
+#include "pylith/topology/Fields.hh" // HOLDSA Fields
#include "pylith/topology/Field.hh" // USES Field
#include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
@@ -76,7 +76,7 @@
// ----------------------------------------------------------------------
// Get parameter fields.
-const pylith::topology::FieldsNew<pylith::topology::SubMesh>*
+const pylith::topology::Fields<pylith::topology::Field<pylith::topology::SubMesh> >*
pylith::faults::TractPerturbation::parameterFields(void) const
{ // parameterFields
return _parameters;
@@ -102,23 +102,39 @@
const int spaceDim = cs->spaceDim();
delete _parameters;
- _parameters = new topology::FieldsNew<topology::SubMesh>(faultMesh);
+ _parameters = new topology::Fields<topology::Field<topology::SubMesh> >(faultMesh);
// Create section to hold time dependent values
- _parameters->add("value", "traction", spaceDim, topology::FieldBase::VECTOR, pressureScale);
- if (_dbInitial)
- _parameters->add("initial", "traction_initial", spaceDim, topology::FieldBase::VECTOR, pressureScale);
+ _parameters->add("value", "traction", topology::FieldBase::VERTICES_FIELD, spaceDim);
+ _parameters->get("value").vectorFieldType(topology::FieldBase::VECTOR);
+ _parameters->get("value").scale(pressureScale);
+ _parameters->get("value").allocate();
+ if (_dbInitial) {
+ _parameters->add("initial", "traction_initial", topology::FieldBase::VERTICES_FIELD, spaceDim);
+ _parameters->get("value").vectorFieldType(topology::FieldBase::VECTOR);
+ _parameters->get("value").scale(pressureScale);
+ _parameters->get("value").allocate();
+ }
if (_dbRate) {
- _parameters->add("rate", "traction_rate", spaceDim, topology::FieldBase::VECTOR, rateScale);
- _parameters->add("rate time", "rate_start_time", 1, topology::FieldBase::SCALAR, timeScale);
+ _parameters->add("rate", "traction_rate", topology::FieldBase::VERTICES_FIELD, spaceDim);
+ _parameters->get("rate").vectorFieldType(topology::FieldBase::VECTOR);
+ _parameters->get("rate").scale(rateScale);
+ _parameters->get("rate").allocate();
+ _parameters->add("rate time", "rate_start_time", topology::FieldBase::VERTICES_FIELD, 1);
+ _parameters->get("rate time").vectorFieldType(topology::FieldBase::SCALAR);
+ _parameters->get("rate time").scale(timeScale);
+ _parameters->get("rate time").allocate();
} // if
if (_dbChange) {
- _parameters->add("change", "traction_change", spaceDim, topology::FieldBase::VECTOR, pressureScale);
- _parameters->add("change time", "change_start_time", 1, topology::FieldBase::SCALAR, timeScale);
+ _parameters->add("change", "traction_change", topology::FieldBase::VERTICES_FIELD, spaceDim);
+ _parameters->get("change").vectorFieldType(topology::FieldBase::VECTOR);
+ _parameters->get("change").scale(pressureScale);
+ _parameters->get("change").allocate();
+ _parameters->add("change time", "change_start_time", topology::FieldBase::FACES_FIELD, 1);
+ _parameters->get("change time").vectorFieldType(topology::FieldBase::SCALAR);
+ _parameters->get("change time").scale(timeScale);
+ _parameters->get("change time").allocate();
} // if
- _parameters->allocate(topology::FieldBase::VERTICES_FIELD, 0);
- const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
- assert(!parametersSection.isNull());
if (_dbInitial) { // Setup initial values, if provided.
_dbInitial->open();
@@ -241,96 +257,130 @@
const PylithScalar timeScale = _timeScale;
// Get vertices.
- const ALE::Obj<SieveSubMesh>& sieveMesh = _parameters->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
- assert(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = _parameters->mesh().dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
const spatialdata::geocoords::CoordSys* cs = _parameters->mesh().coordsys();
assert(cs);
const int spaceDim = cs->spaceDim();
- const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- scalar_array parametersVertex(parametersFiberDim);
-
- const int valueIndex = _parameters->sectionIndex("value");
- const int valueFiberDim = _parameters->sectionFiberDim("value");
- assert(spaceDim == valueFiberDim);
+ PetscSection initialSection, rateSection, rateTimeSection, changeSection, changeTimeSection;
+ Vec initialVec, rateVec, rateTimeVec, changeVec, changeTimeVec;
+ PetscScalar *valuesArray, *initialArray, *rateArray, *rateTimeArray, *changeArray, *changeTimeArray;
- const int initialIndex = (_dbInitial) ? _parameters->sectionIndex("initial") : -1;
- const int initialFiberDim = (_dbInitial) ? _parameters->sectionFiberDim("initial") : 0;
+ PetscSection valuesSection = _parameters->get("value").petscSection();
+ Vec valuesVec = _parameters->get("value").localVector();
+ assert(valuesSection);assert(valuesVec);
+ err = VecGetArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ if (_dbInitial) {
+ initialSection = _parameters->get("initial").petscSection();
+ initialVec = _parameters->get("initial").localVector();
+ assert(initialSection);assert(initialVec);
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbRate) {
+ rateSection = _parameters->get("rate").petscSection();
+ rateTimeSection = _parameters->get("rate time").petscSection();
+ rateVec = _parameters->get("rate").localVector();
+ rateTimeVec = _parameters->get("rate time").localVector();
+ assert(rateSection);assert(rateTimeSection);assert(rateVec);assert(rateTimeVec);
+ err = VecGetArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbChange) {
+ changeSection = _parameters->get("change").petscSection();
+ changeTimeSection = _parameters->get("change time").petscSection();
+ changeVec = _parameters->get("change").localVector();
+ changeTimeVec = _parameters->get("change time").localVector();
+ assert(changeSection);assert(changeTimeSection);assert(changeVec);assert(changeTimeVec);
+ err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
- const int rateIndex = (_dbRate) ? _parameters->sectionIndex("rate") : -1;
- const int rateFiberDim = (_dbRate) ? _parameters->sectionFiberDim("rate") : 0;
- const int rateTimeIndex = (_dbRate) ? _parameters->sectionIndex("rate time") : -1;
- const int rateTimeFiberDim = (_dbRate) ? _parameters->sectionFiberDim("rate time") : 0;
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt vdof, voff;
- const int changeIndex = (_dbChange) ? _parameters->sectionIndex("change") : -1;
- const int changeFiberDim = (_dbChange) ? _parameters->sectionFiberDim("change") : 0;
- const int changeTimeIndex = (_dbChange) ? _parameters->sectionIndex("change time") : -1;
- const int changeTimeFiberDim = (_dbChange) ? _parameters->sectionFiberDim("change time") : 0;
+ err = PetscSectionGetDof(valuesSection, v, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valuesSection, v, &voff);CHECK_PETSC_ERROR(err);
+ assert(vdof == spaceDim);
+ for(PetscInt d = 0; d < vdof; ++d)
+ valuesArray[voff+d] = 0.0;
- for(SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter != verticesEnd; ++v_iter) {
- assert(parametersFiberDim == parametersSection->getFiberDimension(*v_iter));
- parametersSection->restrictPoint(*v_iter, ¶metersVertex[0], parametersVertex.size());
- for (int i=0; i < valueFiberDim; ++i)
- parametersVertex[valueIndex+i] = 0.0;
-
// Contribution from initial value
if (_dbInitial) {
- assert(initialIndex >= 0);
- assert(initialFiberDim == valueFiberDim);
- for (int i=0; i < initialFiberDim; ++i)
- parametersVertex[valueIndex+i] += parametersVertex[initialIndex+i];
+ PetscInt idof, ioff;
+
+ err = PetscSectionGetDof(initialSection, v, &idof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(initialSection, v, &ioff);CHECK_PETSC_ERROR(err);
+ assert(idof == vdof);
+ for(PetscInt d = 0; d < idof; ++d)
+ valuesArray[voff+d] += initialArray[ioff+d];
} // if
// Contribution from rate of change of value
if (_dbRate) {
- assert(rateIndex >= 0);
- assert(rateFiberDim == valueFiberDim);
- assert(rateTimeIndex >= 0);
- assert(rateTimeFiberDim == 1);
-
- const PylithScalar tRel = t - parametersVertex[rateTimeIndex];
+ PetscInt rdof, roff, rtdof, rtoff;
+
+ err = PetscSectionGetDof(rateSection, v, &rdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(rateSection, v, &roff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(rateTimeSection, v, &rtdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(rateTimeSection, v, &rtoff);CHECK_PETSC_ERROR(err);
+ assert(rdof == vdof);
+ assert(rtdof == 1);
+
+ const PylithScalar tRel = t - rateTimeArray[rtoff];
if (tRel > 0.0) // rate of change integrated over time
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- parametersVertex[valueIndex+iDim] += parametersVertex[rateIndex+iDim] * tRel;
- } // for
+ for(PetscInt d = 0; d < spaceDim; ++d)
+ valuesArray[voff+d] += rateArray[roff+d] * tRel;
} // if
// Contribution from change of value
if (_dbChange) {
- assert(changeIndex >= 0);
- assert(changeFiberDim == valueFiberDim);
- assert(changeTimeIndex >= 0);
- assert(changeTimeFiberDim == 1);
+ PetscInt cdof, coff, ctdof, ctoff;
- const PylithScalar tRel = t - parametersVertex[changeTimeIndex];
+ err = PetscSectionGetDof(changeSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(changeSection, v, &coff);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetDof(changeTimeSection, v, &ctdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(changeTimeSection, v, &ctoff);CHECK_PETSC_ERROR(err);
+ assert(cdof == vdof);
+ assert(ctdof == 1);
+
+ const PylithScalar tRel = t - changeTimeArray[ctoff];
if (tRel >= 0) { // change in value over time
- PylithScalar scale = 1.0;
- if (_dbTimeHistory) {
- PylithScalar tDim = tRel*timeScale;
- const int err = _dbTimeHistory->query(&scale, tDim);
- if (err) {
- std::ostringstream msg;
- msg << "Error querying for time '" << tDim
- << "' in time history database "
- << _dbTimeHistory->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
- } // if
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- parametersVertex[valueIndex+iDim] += parametersVertex[changeIndex+iDim] * scale;
- } // for
+ PylithScalar scale = 1.0;
+ if (0 != _dbTimeHistory) {
+ PylithScalar tDim = tRel * timeScale;
+ /* _getNormalizer().dimensionalize(&tDim, 1, timeScale);*/
+ const int err = _dbTimeHistory->query(&scale, tDim);
+ if (0 != err) {
+ std::ostringstream msg;
+ msg << "Error querying for time '" << tDim
+ << "' in time history database "
+ << _dbTimeHistory->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ } // if
+ for(PetscInt d = 0; d < spaceDim; ++d)
+ valuesArray[voff+d] += changeArray[coff+d] * scale;
} // if
} // if
-
- parametersSection->updatePoint(*v_iter, ¶metersVertex[0]);
} // for
+ err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ if (_dbInitial) {
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbRate) {
+ err = VecRestoreArray(rateVec, &rateArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(rateTimeVec, &rateTimeArray);CHECK_PETSC_ERROR(err);
+ }
+ if (_dbChange) {
+ err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ }
} // calculate
@@ -409,13 +459,13 @@
assert(_parameters);
// Get vertices.
- const ALE::Obj<SieveSubMesh>& sieveMesh = _parameters->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
- assert(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = _parameters->mesh().dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ assert(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
const spatialdata::geocoords::CoordSys* cs = _parameters->mesh().coordsys();
assert(cs);
const int spaceDim = cs->spaceDim();
@@ -425,33 +475,42 @@
// Containers for database query results and quadrature coordinates in
// reference geometry.
scalar_array valuesVertex(querySize);
- scalar_array coordsVertexGlobal(spaceDim);
// Get sections.
- const ALE::Obj<RealSection>& coordsSection = sieveMesh->getRealSection("coordinates");
- assert(!coordsSection.isNull());
+ scalar_array coordsVertexGlobal(spaceDim);
+ PetscSection coordSection;
+ Vec coordVec;
+ PetscScalar *coordArray;
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ assert(coordSection);assert(coordVec);
- const ALE::Obj<SubRealUniformSection>& parametersSection = _parameters->section();
- assert(!parametersSection.isNull());
- const int parametersFiberDim = _parameters->fiberDim();
- const int valueIndex = _parameters->sectionIndex(name);
- const int valueFiberDim = _parameters->sectionFiberDim(name);
- assert(valueIndex+valueFiberDim <= parametersFiberDim);
- scalar_array parametersVertex(parametersFiberDim);
+ PetscSection valueSection = _parameters->get("value").petscSection();
+ Vec valueVec = _parameters->get("value").localVector();
+ PetscScalar *tractionsArray;
+ assert(valueSection);assert(valueVec);
// Loop over cells in boundary mesh and perform queries.
- for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin; v_iter != verticesEnd; ++v_iter) {
- assert(spaceDim == coordsSection->getFiberDimension(*v_iter));
- coordsSection->restrictPoint(*v_iter, &coordsVertexGlobal[0], coordsVertexGlobal.size());
+ err = VecGetArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(valueVec, &tractionsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt cdof, coff;
+
+ err = PetscSectionGetDof(coordSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(coordSection, v, &coff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == cdof);
+ for(PetscInt d = 0; d < cdof; ++d) {
+ coordsVertexGlobal[d] = coordArray[coff+d];
+ }
normalizer.dimensionalize(&coordsVertexGlobal[0], coordsVertexGlobal.size(), lengthScale);
valuesVertex = 0.0;
- const int err = db->query(&valuesVertex[0], querySize, &coordsVertexGlobal[0], spaceDim, cs);
+ err = db->query(&valuesVertex[0], querySize, &coordsVertexGlobal[0], spaceDim, cs);
if (err) {
std::ostringstream msg;
msg << "Could not find values at (";
for (int i=0; i < spaceDim; ++i)
- msg << " " << coordsVertexGlobal[i];
+ msg << " " << coordsVertexGlobal[i];
msg << ") for traction boundary condition " << _label << "\n"
<< "using spatial database " << db->label() << ".";
throw std::runtime_error(msg.str());
@@ -460,13 +519,15 @@
normalizer.nondimensionalize(&valuesVertex[0], valuesVertex.size(), scale);
// Update section
- assert(parametersFiberDim == parametersSection->getFiberDimension(*v_iter));
- parametersSection->restrictPoint(*v_iter, ¶metersVertex[0], parametersVertex.size());
- for (int i=0; i < valueFiberDim; ++i)
- parametersVertex[valueIndex+i] = valuesVertex[i];
-
- parametersSection->updatePoint(*v_iter, ¶metersVertex[0]);
+ PetscInt vdof, voff;
+
+ err = PetscSectionGetDof(valueSection, v, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, v, &voff);CHECK_PETSC_ERROR(err);
+ assert(vdof == spaceDim);
+ for(PetscInt d = 0; d < vdof; ++d)
+ tractionsArray[voff+d] = valuesVertex[d];
} // for
+ err = VecRestoreArray(coordVec, &coordArray);CHECK_PETSC_ERROR(err);
} // _queryDB
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.hh 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.hh 2012-11-10 14:41:39 UTC (rev 21009)
@@ -66,7 +66,7 @@
*
* @returns Parameter fields.
*/
- const topology::FieldsNew<topology::SubMesh>* parameterFields(void) const;
+ const topology::Fields<topology::Field<topology::SubMesh> >* parameterFields(void) const;
/** Initialize slip time function.
*
@@ -129,7 +129,7 @@
private :
/// Parameters for perturbations.
- topology::FieldsNew<topology::SubMesh>* _parameters;
+ topology::Fields<topology::Field<topology::SubMesh> >* _parameters;
/// Time scale for current time.
PylithScalar _timeScale;
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-11-10 14:41:39 UTC (rev 21009)
@@ -137,79 +137,74 @@
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
- int iVertex = 0;
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- CPPUNIT_ASSERT(renumbering.find(_data->constraintVertices[iVertex]) !=
- renumbering.end());
- CPPUNIT_ASSERT_EQUAL(renumbering[_data->constraintVertices[iVertex]],
- *v_iter);
+ DM dmMesh = fault._faultMesh->dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscInt vStart, vEnd, numPoints;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(dmMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT(subpointMap);
+ err = ISGetSize(subpointMap, &numPoints);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt faultPoint;
+
+ err = PetscFindInt(_data->constraintVertices[v-vStart], numPoints, points, &faultPoint);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT(faultPoint >= 0);
+ CPPUNIT_ASSERT_EQUAL(faultPoint, v);
} // for
- CPPUNIT_ASSERT_EQUAL(_data->numConstraintVert, iVertex);
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(_data->numConstraintVert, vEnd-vStart);
// Check orientation
//fault._fields->get("orientation").view("ORIENTATION"); // DEBUGGING
-
- const ALE::Obj<RealSection>& orientationSection =
- fault._fields->get("orientation").section();
- CPPUNIT_ASSERT(!orientationSection.isNull());
+ PetscSection orientationSection = fault._fields->get("orientation").petscSection();
+ Vec orientationVec = fault._fields->get("orientation").localVector();
+ PetscScalar *orientationArray;
+ CPPUNIT_ASSERT(orientationSection);CPPUNIT_ASSERT(orientationVec);
const int spaceDim = _data->spaceDim;
const int orientationSize = spaceDim*spaceDim;
- iVertex = 0;
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- const int fiberDim = orientationSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(orientationSize, fiberDim);
- const PylithScalar* orientationVertex =
- orientationSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != orientationVertex);
+ int iVertex = 0;
+ err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
+ err = PetscSectionGetDof(orientationSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(orientationSize, dof);
+
const PylithScalar tolerance = 1.0e-06;
- for (int i=0; i < orientationSize; ++i) {
- const int index = iVertex*orientationSize+i;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->orientation[index],
- orientationVertex[i], tolerance);
+ for(PetscInt d = 0; d < orientationSize; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->orientation[iVertex*orientationSize+d], orientationArray[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
// Prescribed traction perturbation
if (fault._tractPerturbation) {
// :KLUDGE: Only check initial value
- const ALE::Obj<RealSection>& initialTractionsSection =
- fault.vertexField("traction_initial_value").section();
- CPPUNIT_ASSERT(!initialTractionsSection.isNull());
-
- //initialTractionsSection->view("INITIAL TRACTIONS"); // DEBUGGING
-
- const int spaceDim = _data->spaceDim;
+ PetscSection initialTractionsSection = fault.vertexField("traction_initial_value").petscSection();
+ Vec initialTractionsVec = fault.vertexField("traction_initial_value").localVector();
+ PetscScalar *initialTractionsArray;
+ CPPUNIT_ASSERT(initialTractionsSection);CPPUNIT_ASSERT(initialTractionsVec);
iVertex = 0;
- for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- const int fiberDim = initialTractionsSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const PylithScalar* initialTractionsVertex =
- initialTractionsSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(initialTractionsVertex);
+ err = VecGetArray(initialTractionsVec, &initialTractionsArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
+ err = PetscSectionGetDof(initialTractionsSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(initialTractionsSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, dof);
+
const PylithScalar tolerance = 1.0e-06;
- for (int i = 0; i < spaceDim; ++i) {
- const int index = iVertex * spaceDim + i;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->initialTractions[index],
- initialTractionsVertex[i], tolerance);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->initialTractions[iVertex * spaceDim + d], initialTractionsArray[d], tolerance);
} // for
} // for
+ err = VecRestoreArray(initialTractionsVec, &initialTractionsArray);CHECK_PETSC_ERROR(err);
} // if
} // testInitialize
@@ -224,7 +219,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields.solution(), "seqdense");
+ topology::Jacobian jacobian(fields.solution());
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrStick);
const int spaceDim = _data->spaceDim;
@@ -242,97 +237,85 @@
{ // Check solution values
// No change to Lagrange multipliers for stick case.
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = mesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Get section containing solution (disp + Lagrange multipliers)
- const ALE::Obj<RealSection>& dispIncrSection =
- fields.get("dispIncr(t->t+dt)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
+ PetscSection dispTIncrSection = fields.get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields.get("dispIncr(t->t+dt)").localVector();
+ PetscScalar *dispTIncrArray;
+ assert(dispTIncrSection);assert(dispTIncrVec);
- //dispIncrSection->view("DISP INCREMENT"); // DEBUGGING
-
// Get expected values
const PylithScalar* valsE = _data->fieldIncrStick; // No change in dispIncr
int iVertex = 0; // variable to use as index into valsE array
const int fiberDimE = spaceDim; // number of values per point
const PylithScalar tolerance = 1.0e-06;
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) { // loop over all vertices in mesh
- // Check fiber dimension (number of values at point)
- const int fiberDim = dispIncrSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = dispIncrSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ err = VecGetArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
+ err = PetscSectionGetDof(dispTIncrSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTIncrSection, v, &off);CHECK_PETSC_ERROR(err);
+ assert(fiberDimE == dof);
+
// Check values at point
- for (int i = 0; i < fiberDimE; ++i) {
- const int index = iVertex * spaceDim + i;
- const PylithScalar valE = valsE[index];
+ for(PetscInt d = 0; d < fiberDimE; ++d) {
+ const PylithScalar valE = valsE[iVertex*spaceDim+d];
if (fabs(valE) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, dispTIncrArray[off+d]/valE, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, dispTIncrArray[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
} // Check solution values
{ // Check slip values
// Slip should be zero for the stick case.
// Get fault vertex info
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM faultDMMesh = fault._faultMesh->dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(faultDMMesh);
+ err = DMComplexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Get section containing slip
- const ALE::Obj<RealSection>& slipSection =
- fault.vertexField("slip").section();
- CPPUNIT_ASSERT(!slipSection.isNull());
+ PetscSection slipSection = fault.vertexField("slip").petscSection();
+ Vec slipVec = fault.vertexField("slip").localVector();
+ PetscScalar *slipArray;
+ CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
- //slipSection->view("SLIP"); // DEBUGGING
-
// Get expected values
CPPUNIT_ASSERT(_data->slipStickE);
const PylithScalar* valsE = _data->slipStickE;
int iVertex = 0; // variable to use as index into valsE array
const int fiberDimE = spaceDim; // number of values per point
const PylithScalar tolerance = 1.0e-06;
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) { // loop over fault vertices
- // Check fiber dimension (number of values at point)
- const int fiberDim = slipSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(vals);
+ err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
+ err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
+
// Check values at point
- for (int i = 0; i < fiberDimE; ++i) {
- const int index = iVertex * spaceDim + i;
- const PylithScalar valE = valsE[index];
-#if 0 // DEBUGGING
- std::cout << "SLIP valE: " << valE
- << ", val: " << vals[i]
- << ", error: " << fabs(1.0-vals[i]/valE)
- << std::endl;
-#endif // DEBUGGING
+ for(PetscInt d = 0; d < dof; ++d) {
+ const PylithScalar valE = valsE[iVertex*spaceDim+d];
if (fabs(valE) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, slipArray[off+d]/valE, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
- } // for
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, slipArray[off+d], tolerance);
+ }
} // for
+ err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
} // Check slip values
} // testConstrainSolnSpaceStick
@@ -348,7 +331,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields.solution(), "seqdense");
+ topology::Jacobian jacobian(fields.solution());
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrSlip);
const int spaceDim = _data->spaceDim;
@@ -367,49 +350,40 @@
{ // Check solution values
// Lagrange multipliers should be adjusted according to friction
// as reflected in the fieldIncrSlipE data member.
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = mesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Get section containing solution (disp + Lagrange multipliers)
- const ALE::Obj<RealSection>& dispIncrSection =
- fields.get("dispIncr(t->t+dt)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
+ PetscSection dispIncrSection = fields.get("dispIncr(t->t+dt)").petscSection();
+ Vec dispIncrVec = fields.get("dispIncr(t->t+dt)").localVector();
+ PetscScalar *dispIncrArray;
+ CPPUNIT_ASSERT(dispIncrSection);CPPUNIT_ASSERT(dispIncrVec);
// Get expected values
const PylithScalar* valsE = _data->fieldIncrSlipE; // Expected values for dispIncr
int iVertex = 0; // variable to use as index into valsE array
const int fiberDimE = spaceDim; // number of values per point
const PylithScalar tolerance = 1.0e-06;
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) { // loop over all vertices in mesh
- // Check fiber dimension (number of values at point)
- const int fiberDim = dispIncrSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = dispIncrSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ err = VecGetArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
- // Check values at point
- for (int i = 0; i < fiberDimE; ++i) {
- const int index = iVertex * spaceDim + i;
- const PylithScalar valE = valsE[index];
-#if 0 // DEBUGGING
- std::cout << "SOLUTION valE: " << valE
- << ", val: " << vals[i]
- << ", error: " << fabs(1.0-vals[i]/valE)
- << std::endl;
-#endif // DEBUGGING
- if (fabs(valE) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ err = PetscSectionGetDof(dispIncrSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispIncrSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
+ for(PetscInt d = 0; d < dof; ++d) {
+ const PylithScalar valE = valsE[iVertex*spaceDim+d];
+ if (fabs(valE) > tolerance)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, dispIncrArray[off+d]/valE, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
- } // for
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, dispIncrArray[off+d], tolerance);
+ }
} // for
+ err = VecRestoreArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
} // Check solution values
{ // Check slip values
@@ -417,52 +391,40 @@
// Lagrange multipliers as reflected in the slipSlipE data member.
// Get fault vertex info
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM faultDMMesh = fault._faultMesh->dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(faultDMMesh);
+ err = DMComplexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Get section containing slip
- const ALE::Obj<RealSection>& slipSection =
- fault.vertexField("slip").section();
- CPPUNIT_ASSERT(!slipSection.isNull());
+ PetscSection slipSection = fault.vertexField("slip").petscSection();
+ Vec slipVec = fault.vertexField("slip").localVector();
+ PetscScalar *slipArray;
+ CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
- //slipSection->view("SLIP"); // DEBUGGING
-
// Get expected values
const PylithScalar* valsE = _data->slipSlipE;
int iVertex = 0; // variable to use as index into valsE array
const int fiberDimE = spaceDim; // number of values per point
- const PylithScalar tolerance =
- (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-5;
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
- != verticesEnd;
- ++v_iter, ++iVertex) { // loop over fault vertices
- // Check fiber dimension (number of values at point)
- const int fiberDim = slipSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-5;
+ err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
- // Check values at point
- for (int i = 0; i < fiberDimE; ++i) {
- const int index = iVertex * spaceDim + i;
- const PylithScalar valE = valsE[index];
-#if 0 // DEBUGGING
- std::cout << "SLIP valE: " << valE
- << ", val: " << vals[i]
- << ", error: " << fabs(1.0-vals[i]/valE)
- << std::endl;
-#endif // DEBUGGING
+ err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
+ for(PetscInt d = 0; d < dof; ++d) {
+ const PylithScalar valE = valsE[iVertex*spaceDim+d];
if (fabs(valE) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, slipArray[off+d]/valE, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
- } // for
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, slipArray[off+d], tolerance);
+ }
} // for
+ err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
} // Check slip values
} // testConstrainSolnSpaceSlip
@@ -478,7 +440,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields.solution(), "seqdense");
+ topology::Jacobian jacobian(fields.solution());
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrOpen);
const int spaceDim = _data->spaceDim;
@@ -499,48 +461,40 @@
{ // Check solution values
// Lagrange multipliers should be set to zero as reflected in the
// fieldIncrOpenE data member.
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = mesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Get section containing solution (disp + Lagrange multipliers)
- const ALE::Obj<RealSection>& dispIncrSection =
- fields.get("dispIncr(t->t+dt)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
+ PetscSection dispIncrSection = fields.get("dispIncr(t->t+dt)").petscSection();
+ Vec dispIncrVec = fields.get("dispIncr(t->t+dt)").localVector();
+ PetscScalar *dispIncrArray;
+ CPPUNIT_ASSERT(dispIncrSection);CPPUNIT_ASSERT(dispIncrVec);
// Get expected values
const PylithScalar* valsE = _data->fieldIncrOpenE; // Expected values for dispIncr
int iVertex = 0; // variable to use as index into valsE array
const int fiberDimE = spaceDim; // number of values per point
const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-05;
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) { // loop over all vertices in mesh
- // Check fiber dimension (number of values at point)
- const int fiberDim = dispIncrSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = dispIncrSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ err = VecGetArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
- // Check values at point
- for (int i = 0; i < fiberDimE; ++i) {
- const int index = iVertex * spaceDim + i;
- const PylithScalar valE = valsE[index];
-#if 0 // DEBUGGING
- std::cout << "valE: " << valE
- << ", val: " << vals[i]
- << std::endl;
-#endif // DEBUGGING
+ err = PetscSectionGetDof(dispIncrSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispIncrSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
+ for(PetscInt d = 0; d < dof; ++d) {
+ const PylithScalar valE = valsE[iVertex*spaceDim+d];
if (fabs(valE) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, dispIncrArray[off+d]/valE, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
- } // for
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, dispIncrArray[off+d], tolerance);
+ }
} // for
+ err = VecRestoreArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
} // Check solution values
{ // Check slip values
@@ -548,48 +502,40 @@
// Lagrange multipliers as reflected in the slipOpenE data member.
// Get fault vertex info
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM faultDMMesh = fault._faultMesh->dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(faultDMMesh);
+ err = DMComplexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Get section containing slip
- const ALE::Obj<RealSection>& slipSection =
- fault.vertexField("slip").section();
- CPPUNIT_ASSERT(!slipSection.isNull());
+ PetscSection slipSection = fault.vertexField("slip").petscSection();
+ Vec slipVec = fault.vertexField("slip").localVector();
+ PetscScalar *slipArray;
+ CPPUNIT_ASSERT(slipSection);CPPUNIT_ASSERT(slipVec);
// Get expected values
const PylithScalar* valsE = _data->slipOpenE;
int iVertex = 0; // variable to use as index into valsE array
const int fiberDimE = spaceDim; // number of values per point
const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-05;
- for (SieveMesh::label_sequence::iterator v_iter = verticesBegin; v_iter
- != verticesEnd;
- ++v_iter, ++iVertex) { // loop over fault vertices
- // Check fiber dimension (number of values at point)
- const int fiberDim = slipSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = slipSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ err = VecGetArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
- // Check values at point
- for (int i = 0; i < fiberDimE; ++i) {
- const int index = iVertex * spaceDim + i;
- const PylithScalar valE = valsE[index];
-#if 0 // DEBUGGING
- std::cout << "valE: " << valE
- << ", val: " << vals[i]
- << std::endl;
-#endif // DEBUGGING
+ err = PetscSectionGetDof(slipSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(slipSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
+ for(PetscInt d = 0; d < dof; ++d) {
+ const PylithScalar valE = valsE[iVertex*spaceDim+d];
if (fabs(valE) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, slipArray[off+d]/valE, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
- } // for
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, slipArray[off+d], tolerance);
+ }
} // for
+ err = VecRestoreArray(slipVec, &slipArray);CHECK_PETSC_ERROR(err);
} // Check slip values
} // testConstrainSolnSpaceOpen
@@ -605,7 +551,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields.solution(), "seqdense");
+ topology::Jacobian jacobian(fields.solution());
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrSlip);
const int spaceDim = _data->spaceDim;
@@ -630,7 +576,7 @@
FaultCohesiveDyn fault;
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- topology::Jacobian jacobian(fields.solution(), "seqdense");
+ topology::Jacobian jacobian(fields.solution());
_setFieldsJacobian(&mesh, &fault, &fields, &jacobian, _data->fieldIncrStick);
CPPUNIT_ASSERT(0 != fault._faultMesh);
@@ -639,77 +585,64 @@
tractions.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
tractions.allocate();
tractions.zero();
- const ALE::Obj<RealSection>& tractionsSection = tractions.section();
- CPPUNIT_ASSERT(!tractionsSection.isNull());
+ PetscSection tractionsSection = tractions.petscSection();
+ Vec tractionsVec = tractions.localVector();
+ PetscScalar *tractionsArray;
+ CPPUNIT_ASSERT(tractionsSection);CPPUNIT_ASSERT(tractionsVec);
const PylithScalar t = 0;
fault.updateStateVars(t, &fields);
fault._calcTractions(&tractions, fields.get("disp(t)"));
- //tractions.view("TRACTIONS"); // DEBUGGING
+ DM faultDMMesh = fault._faultMesh->dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscInt vStart, vEnd, numPoints;
+ PetscErrorCode err;
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin =
- vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
- const SieveMesh::renumbering_type::const_iterator renumberingBegin =
- renumbering.begin();
- const SieveMesh::renumbering_type::const_iterator renumberingEnd =
- renumbering.end();
+ CPPUNIT_ASSERT(faultDMMesh);
+ err = DMComplexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(faultDMMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT(subpointMap);
+ err = ISGetSize(subpointMap, &numPoints);CHECK_PETSC_ERROR(err);
- const ALE::Obj<RealSection>& dispSection = fields.get("disp(t)").section();
- CPPUNIT_ASSERT(!dispSection.isNull());
+ PetscSection dispSection = fields.get("disp(t)").petscSection();
+ Vec dispVec = fields.get("disp(t)").localVector();
+ PetscScalar *dispArray;
+ CPPUNIT_ASSERT(dispSection);CPPUNIT_ASSERT(dispVec);
int iVertex = 0;
const PylithScalar tolerance = 1.0e-06;
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- SieveMesh::point_type meshVertex = -1;
- bool found = false;
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(tractionsVec, &tractionsArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispVec, &dispArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt tdof, toff, ddof, doff;
- for (SieveMesh::renumbering_type::const_iterator r_iter = renumberingBegin;
- r_iter != renumberingEnd;
- ++r_iter) {
- if (r_iter->second == *v_iter) {
- meshVertex = r_iter->first;
- found = true;
- break;
- } // if
- } // for
- CPPUNIT_ASSERT(found);
- int fiberDim = tractionsSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const PylithScalar* tractionsVertex = tractionsSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(tractionsVertex);
+ err = PetscSectionGetDof(tractionsSection, v, &tdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(tractionsSection, v, &toff);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, tdof);
+ err = PetscSectionGetDof(dispSection, points[v-vStart], &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispSection, points[v-vStart], &doff);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, ddof);
- const PylithScalar* tractionsVertexGlobalE =
- dispSection->restrictPoint(meshVertex);
- CPPUNIT_ASSERT(tractionsVertexGlobalE);
- const PylithScalar* orientationVertex =
- &_data->orientation[iVertex*spaceDim*spaceDim];
+ const PylithScalar *orientationVertex = &_data->orientation[iVertex*spaceDim*spaceDim];
CPPUNIT_ASSERT(orientationVertex);
- for (int iDim=0; iDim < spaceDim; ++iDim) {
+ for(PetscInt d = 0; d < spaceDim; ++d) {
PylithScalar tractionE = 0.0;
- for (int jDim=0; jDim < spaceDim; ++jDim) {
- tractionE +=
- orientationVertex[iDim*spaceDim+jDim]*tractionsVertexGlobalE[jDim];
+ for(PetscInt e = 0; e < spaceDim; ++e) {
+ tractionE += orientationVertex[d*spaceDim+e]*dispArray[doff+e];
} // for
if (tractionE != 0.0)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, tractionsVertex[iDim]/tractionE,
- tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, tractionsArray[toff+d]/tractionE, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionE, tractionsVertex[iDim],
- tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(tractionE, tractionsArray[toff+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(tractionsVec, &tractionsArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispVec, &dispArray);CHECK_PETSC_ERROR(err);
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
} // testCalcTractions
// ----------------------------------------------------------------------
@@ -768,11 +701,15 @@
_friction = friction;
fault->frictionModel(friction);
- int firstFaultVertex = 0;
- int firstLagrangeVertex = mesh->sieveMesh()->getIntSection(_data->label)->size();
- int firstFaultCell = mesh->sieveMesh()->getIntSection(_data->label)->size();
+ PetscInt labelSize;
+ PetscErrorCode err;
+ err = DMComplexGetStratumSize(mesh->dmMesh(), _data->label, 1, &labelSize);CHECK_PETSC_ERROR(err);
+
+ PetscInt firstFaultVertex = 0;
+ PetscInt firstLagrangeVertex = labelSize;
+ PetscInt firstFaultCell = labelSize;
if (fault->useLagrangeConstraints())
- firstFaultCell += mesh->sieveMesh()->getIntSection(_data->label)->size();
+ firstFaultCell += labelSize;
fault->id(_data->id);
fault->label(_data->label);
fault->quadrature(_quadrature);
@@ -827,37 +764,45 @@
const int spaceDim = _data->spaceDim;
// Get vertices in mesh
- const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = mesh->dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
// Set displacement values
- const ALE::Obj<RealSection>& dispSection =
- fields->get("disp(t)").section();
- CPPUNIT_ASSERT(!dispSection.isNull());
+ PetscSection dispSection = fields->get("disp(t)").petscSection();
+ Vec dispVec = fields->get("disp(t)").localVector();
+ PetscScalar *dispArray;
+ CPPUNIT_ASSERT(dispSection);CPPUNIT_ASSERT(dispVec);
int iVertex = 0;
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex)
- dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ err = VecGetArray(dispVec, &dispArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt off;
+ err = PetscSectionGetOffset(dispSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < spaceDim; ++d) dispArray[off+d] = _data->fieldT[iVertex*spaceDim+d];
+ }
+ err = VecRestoreArray(dispVec, &dispArray);CHECK_PETSC_ERROR(err);
// Set increment values
- const ALE::Obj<RealSection>& dispIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- CPPUNIT_ASSERT(!dispIncrSection.isNull());
+ PetscSection dispIncrSection = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ PetscScalar *dispIncrArray;
+ CPPUNIT_ASSERT(dispIncrSection);CPPUNIT_ASSERT(dispIncrVec);
iVertex = 0;
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex)
- dispIncrSection->updatePoint(*v_iter, &fieldIncr[iVertex*spaceDim]);
+ err = VecGetArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt off;
+ err = PetscSectionGetOffset(dispIncrSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < spaceDim; ++d) dispIncrArray[off+d] = fieldIncr[iVertex*spaceDim+d];
+ }
+ err = VecRestoreArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
// Setup Jacobian matrix
- const int nrows = dispIncrSection->sizeWithBC();
- const int ncols = nrows;
+ PetscInt nrows;
+ err = PetscSectionGetStorageSize(dispIncrSection, &nrows);CHECK_PETSC_ERROR(err);
+ PetscInt ncols = nrows;
int nrowsM = 0;
int ncolsM = 0;
PetscMat jacobianMat = jacobian->matrix();
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc 2012-11-10 14:41:39 UTC (rev 21009)
@@ -153,71 +153,74 @@
topology::SolutionFields fields(mesh);
_initialize(&mesh, &fault, &fields);
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = fault._faultMesh->sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- SieveSubMesh::renumbering_type& renumbering =
- faultSieveMesh->getRenumbering();
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
- int iVertex = 0;
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- CPPUNIT_ASSERT(renumbering.find(_data->constraintVertices[iVertex]) !=
- renumbering.end());
- CPPUNIT_ASSERT_EQUAL(renumbering[_data->constraintVertices[iVertex]],
- *v_iter);
+ DM dmMesh = fault._faultMesh->dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscInt vStart, vEnd, numPoints;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(dmMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT(subpointMap);
+ err = ISGetSize(subpointMap, &numPoints);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt faultPoint;
+
+ err = PetscFindInt(_data->constraintVertices[v-vStart], numPoints, points, &faultPoint);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT(faultPoint >= 0);
+ CPPUNIT_ASSERT_EQUAL(faultPoint, v);
} // for
- CPPUNIT_ASSERT_EQUAL(_data->numConstraintVert, iVertex);
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(_data->numConstraintVert, vEnd-vStart);
// Check orientation
//fault._fields->get("orientation").view("ORIENTATION"); // DEBUGGING
- const ALE::Obj<RealSection>& orientationSection =
- fault._fields->get("orientation").section();
- CPPUNIT_ASSERT(!orientationSection.isNull());
+ PetscSection orientationSection = fault._fields->get("orientation").petscSection();
+ Vec orientationVec = fault._fields->get("orientation").localVector();
+ PetscScalar *orientationArray;
+ CPPUNIT_ASSERT(orientationSection);CPPUNIT_ASSERT(orientationVec);
const int spaceDim = _data->spaceDim;
const int orientationSize = spaceDim*spaceDim;
- iVertex = 0;
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- const int fiberDim = orientationSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(orientationSize, fiberDim);
- const PylithScalar* orientationVertex =
- orientationSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != orientationVertex);
+ int iVertex = 0;
+ err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
+ err = PetscSectionGetDof(orientationSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(orientationSize, dof);
+
const PylithScalar tolerance = 1.0e-06;
- for (int i=0; i < orientationSize; ++i) {
- const int index = iVertex*orientationSize+i;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->orientation[index],
- orientationVertex[i], tolerance);
+ for(PetscInt d = 0; d < orientationSize; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->orientation[iVertex*orientationSize+d], orientationArray[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
// Initial tractions
if (fault._dbImpulseAmp) {
//fault._fields->get("initial traction").view("INITIAL TRACTIONS"); // DEBUGGING
- const ALE::Obj<RealSection>& amplitudeSection =
- fault._fields->get("impulse amplitude").section();
- CPPUNIT_ASSERT(!amplitudeSection.isNull());
+ PetscSection amplitudeSection = fault._fields->get("impulse amplitude").petscSection();
+ Vec amplitudeVec = fault._fields->get("impulse amplitude").localVector();
+ PetscScalar *amplitudeArray;
+ CPPUNIT_ASSERT(amplitudeSection);CPPUNIT_ASSERT(amplitudeVec);
const int spaceDim = _data->spaceDim;
iVertex = 0;
- for (SieveSubMesh::label_sequence::iterator v_iter = verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- const int fiberDim = amplitudeSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(1, fiberDim);
- const PylithScalar* amplitudeVertex = amplitudeSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(amplitudeVertex);
+ err = VecGetArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
+ err = PetscSectionGetDof(amplitudeSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(amplitudeSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(1, dof);
+
const PylithScalar tolerance = 1.0e-06;
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->amplitude[iVertex], amplitudeVertex[0], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->amplitude[iVertex], amplitudeArray[off], tolerance);
} // for
+ err = VecRestoreArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
} // if
} // testInitialize
@@ -237,23 +240,30 @@
const int spaceDim = _data->spaceDim;
topology::Field<topology::Mesh>& residual = fields.get("residual");
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *residualArray;
+ CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
- const ALE::Obj<RealSection>& dispSection = fields.get("disp(t)").section();
- CPPUNIT_ASSERT(!dispSection.isNull());
+ PetscSection dispSection = fields.get("disp(t)").petscSection();
+ Vec dispVec = fields.get("disp(t)").localVector();
+ PetscScalar *dispArray;
+ CPPUNIT_ASSERT(dispSection);CPPUNIT_ASSERT(dispVec);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = mesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
int iVertex = 0;
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex)
- dispSection->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ err = VecGetArray(dispVec, &dispArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt off;
+ err = PetscSectionGetOffset(dispSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < spaceDim; ++d) dispArray[off+d] = _data->fieldT[iVertex*spaceDim+d];
+ }
+ err = VecRestoreArray(dispVec, &dispArray);CHECK_PETSC_ERROR(err);
const PylithScalar t = 1.0;
const PylithScalar dt = 1.0;
@@ -270,21 +280,19 @@
iVertex = 0;
const int fiberDimE = spaceDim;
const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-05;
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- const int fiberDim = residualSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = residualSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vals);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(residualSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(residualSection, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
- for (int i = 0; i < fiberDimE; ++i) {
- const int index = iVertex * spaceDim + i;
- const PylithScalar valE = valsE[index];
+ for(PetscInt d = 0; d < fiberDimE; ++d) {
+ const PylithScalar valE = valsE[iVertex*spaceDim+d];
if (fabs(valE) > tolerance)
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valE, tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, residualArray[off+d]/valE, tolerance);
else
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+d], tolerance);
} // for
} // for
} // Integrate residual with disp increment.
@@ -330,11 +338,15 @@
CPPUNIT_ASSERT(_data->impulseDOF);
fault->impulseDOF(_data->impulseDOF, _data->numComponents);
- int firstFaultVertex = 0;
- int firstLagrangeVertex = mesh->sieveMesh()->getIntSection(_data->label)->size();
- int firstFaultCell = mesh->sieveMesh()->getIntSection(_data->label)->size();
+ PetscInt labelSize;
+ PetscErrorCode err;
+ err = DMComplexGetStratumSize(mesh->dmMesh(), _data->label, 1, &labelSize);CHECK_PETSC_ERROR(err);
+
+ PetscInt firstFaultVertex = 0;
+ PetscInt firstLagrangeVertex = labelSize;
+ PetscInt firstFaultCell = labelSize;
if (fault->useLagrangeConstraints())
- firstFaultCell += mesh->sieveMesh()->getIntSection(_data->label)->size();
+ firstFaultCell += labelSize;
fault->id(_data->id);
fault->label(_data->label);
fault->quadrature(_quadrature);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2012-11-10 14:41:39 UTC (rev 21009)
@@ -156,9 +156,7 @@
err = PetscFindInt(_data->verticesLagrange[v-vStart], numPoints, points, &faultPoint);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT(faultPoint >= 0);
-#if 1
CPPUNIT_ASSERT_EQUAL(faultPoint, v);
-#endif
} // for
err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(_data->numFaultVertices, vEnd-vStart);
@@ -206,7 +204,7 @@
CPPUNIT_ASSERT_EQUAL(orientationSize, dof);
const PylithScalar tolerance = 1.0e-06;
- for(int d = 0; d < orientationSize; ++d) {
+ for(PetscInt d = 0; d < orientationSize; ++d) {
CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->orientation[iVertex*orientationSize+d], orientationArray[off+d], tolerance);
} // for
} // for
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc 2012-11-10 14:41:39 UTC (rev 21009)
@@ -26,7 +26,7 @@
#include "pylith/topology/Mesh.hh" // USES Mesh
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
+#include "pylith/topology/Fields.hh" // USES Fields
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
@@ -141,35 +141,38 @@
CPPUNIT_ASSERT(cs);
const int spaceDim = cs->spaceDim();
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
- CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices = faultSieveMesh->depthStratum(0);
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM faultDMMesh = faultMesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(faultDMMesh);
+ err = DMComplexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+
//traction.view("TRACTION"); // DEBUGGING
const PylithScalar tolerance = 1.0e-06;
int iPoint = 0;
CPPUNIT_ASSERT(tract._parameters);
- const ALE::Obj<RealUniformSection>& paramsSection = tract._parameters->section();
- CPPUNIT_ASSERT(!paramsSection.isNull());
- const int fiberDim = tract._parameters->fiberDim();
- const int valueIndex = tract._parameters->sectionIndex("value");
- const int valueFiberDim = tract._parameters->sectionFiberDim("value");
- CPPUNIT_ASSERT_EQUAL(spaceDim, valueFiberDim);
- CPPUNIT_ASSERT(valueIndex + valueFiberDim < fiberDim);
-
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin(); v_iter != verticesEnd; ++v_iter, ++iPoint) {
- CPPUNIT_ASSERT_EQUAL(fiberDim, paramsSection->getFiberDimension(*v_iter));
- const PylithScalar* vals = paramsSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(vals);
+ PetscSection valuesSection = tract._parameters->get("values").petscSection();
+ Vec valuesVec = tract._parameters->get("values").localVector();
+ PetscScalar *valuesArray;
+ assert(valuesSection);assert(valuesVec);
+
+ err = VecGetArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+ PetscInt vdof, voff;
+
+ err = PetscSectionGetDof(valuesSection, v, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valuesSection, v, &voff);CHECK_PETSC_ERROR(err);
+ assert(vdof == spaceDim);
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- const PylithScalar valueE = tractionE[iPoint*spaceDim+iDim];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, vals[valueIndex+iDim], tolerance);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ const PylithScalar valueE = tractionE[iPoint*spaceDim+d];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, valuesArray[voff+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
} // testCalculate
// ----------------------------------------------------------------------
@@ -188,9 +191,9 @@
TractPerturbation tract;
_initialize(&mesh, &faultMesh, &tract);
- const topology::FieldsNew<topology::SubMesh>* parameters = tract.parameterFields();
- const ALE::Obj<RealUniformSection>& parametersSection = parameters->section();
- CPPUNIT_ASSERT(!parametersSection.isNull());
+ const topology::Fields<topology::Field<topology::SubMesh> >* parameters = tract.parameterFields();
+ //const ALE::Obj<RealUniformSection>& parametersSection = parameters->section();
+ //CPPUNIT_ASSERT(!parametersSection.isNull());
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh.sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
@@ -200,14 +203,14 @@
const PylithScalar tolerance = 1.0e-06;
int iPoint = 0;
for (SieveMesh::label_sequence::iterator v_iter=vertices->begin(); v_iter != verticesEnd; ++v_iter, ++iPoint) {
- const int fiberDim = parametersSection->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = parametersSection->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(vals);
+ //const int fiberDim = parametersSection->getFiberDimension(*v_iter);
+ //CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const PylithScalar* vals;// = parametersSection->restrictPoint(*v_iter);
+ //CPPUNIT_ASSERT(vals);
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar valueE = parametersE[iPoint*fiberDim+iDim];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, vals[iDim], tolerance);
+ for(PetscInt d = 0; d < fiberDimE; ++d) {
+ const PylithScalar valueE = parametersE[iPoint*fiberDimE+d];
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, vals[d], tolerance);
} // for
} // for
} // testParameterFields
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2012-11-09 22:56:24 UTC (rev 21008)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2012-11-10 14:41:39 UTC (rev 21009)
@@ -23,7 +23,8 @@
* 2 -------- 3 -------- 4
*
* After adding cohesive elements
- * 2 -------- 3-6-5 -------- 4
+ * 2
+ * 3 -------- 4-7-6 -------- 5
*/
#include "CohesiveKinDataLine2.hh"
@@ -102,13 +103,13 @@
1
};
const int pylith::faults::CohesiveKinDataLine2::_verticesLagrange[] = {
- 6
+ 7
};
const int pylith::faults::CohesiveKinDataLine2::_verticesPositive[] = {
- 5
+ 6
};
const int pylith::faults::CohesiveKinDataLine2::_verticesNegative[] = {
- 3
+ 4
};
const int pylith::faults::CohesiveKinDataLine2::_numCohesiveCells = 1;
@@ -116,7 +117,7 @@
0
};
const int pylith::faults::CohesiveKinDataLine2::_cellMappingCohesive[] = {
- 7
+ 2
};
More information about the CIG-COMMITS
mailing list