[cig-commits] r21082 - in short/3D/PyLith/trunk: libsrc/pylith/faults libsrc/pylith/meshio libsrc/pylith/topology unittests/libtests/faults unittests/libtests/faults/data unittests/libtests/meshio unittests/libtests/meshio/data
knepley at geodynamics.org
knepley at geodynamics.org
Tue Nov 27 16:18:23 PST 2012
Author: knepley
Date: 2012-11-27 16:18:22 -0800 (Tue, 27 Nov 2012)
New Revision: 21082
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.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/meshio/DataWriterHDF5.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.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/TestTractPerturbation.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterVTKFaultMeshCases.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshHex8.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshTet4.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshTri3.cc
short/3D/PyLith/trunk/unittests/libtests/meshio/data/MeshDataCubitTet.cc
Log:
Most fault tests passing, VTK all but Points, working on HDF5
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -876,6 +876,7 @@
err = VecRestoreArray(coordinatesVec, &coords);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(newCoordinatesVec, &newCoords);CHECK_PETSC_ERROR(err);
err = DMSetCoordinatesLocal(newMesh, newCoordinatesVec);CHECK_PETSC_ERROR(err);
+ err = VecDestroy(&newCoordinatesVec);CHECK_PETSC_ERROR(err);
if (debug) coordinates->view("Coordinates with shadow vertices");
logger.stagePop();
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -659,6 +659,7 @@
const int numVertices = _cohesiveVertices.size();
err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTIncrAdjVec, &dispTIncrAdjArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(sensitivityVec, &sensitivityArray);CHECK_PETSC_ERROR(err);
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
@@ -845,6 +846,7 @@
} // for
err = VecRestoreArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTIncrAdjVec, &dispTIncrAdjArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(sensitivityVec, &sensitivityArray);CHECK_PETSC_ERROR(err);
@@ -971,6 +973,7 @@
err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTIncrAdjVec, &dispTIncrAdjArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(sensitivityVec, &sensitivityArray);CHECK_PETSC_ERROR(err);
err = VecGetArray(sensitivityRelVec, &sensitivityRelArray);CHECK_PETSC_ERROR(err);
@@ -1205,6 +1208,7 @@
} // for
err = VecRestoreArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTIncrAdjVec, &dispTIncrAdjArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(sensitivityVec, &sensitivityArray);CHECK_PETSC_ERROR(err);
err = VecRestoreArray(sensitivityRelVec, &sensitivityRelArray);CHECK_PETSC_ERROR(err);
@@ -1631,7 +1635,6 @@
buffer.label("slip");
FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
-
} else if (0 == strcasecmp("slip_rate", name)) {
const topology::Field<topology::SubMesh>& velRel = _fields->get("relative velocity");
_allocateBufferVectorField();
@@ -1640,49 +1643,38 @@
buffer.label("slip_rate");
FaultCohesiveLagrange::globalToFault(&buffer, orientation);
return buffer;
-
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection = _fields->get("orientation").section();
- assert(!orientationSection.isNull());
- const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(0);
- assert(!dirSection.isNull());
+ PetscSection orientationSection = _fields->get("orientation").petscSection();
+ Vec orientationVec = _fields->get("orientation").localVector();
+ assert(orientationSection);assert(orientationVec);
_allocateBufferVectorField();
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
+ buffer.copy(orientationSection, 0, 0, orientationVec);
buffer.label("strike_dir");
buffer.scale(1.0);
- buffer.copy(dirSection);
return buffer;
} else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection = _fields->get(
- "orientation").section();
- assert(!orientationSection.isNull());
- const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(
- 1);
+ PetscSection orientationSection = _fields->get("orientation").petscSection();
+ Vec orientationVec = _fields->get("orientation").localVector();
+ assert(orientationSection);assert(orientationVec);
_allocateBufferVectorField();
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
+ buffer.copy(orientationSection, 0, 1, orientationVec);
buffer.label("dip_dir");
buffer.scale(1.0);
- buffer.copy(dirSection);
return buffer;
-
} else if (0 == strcasecmp("normal_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection = _fields->get(
- "orientation").section();
- assert(!orientationSection.isNull());
+ PetscSection orientationSection = _fields->get("orientation").petscSection();
+ Vec orientationVec = _fields->get("orientation").localVector();
+ assert(orientationSection);assert(orientationVec);
const int space = (0 == cohesiveDim) ? 0 : (1 == cohesiveDim) ? 1 : 2;
- const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(
- space);
- assert(!dirSection.isNull());
_allocateBufferVectorField();
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
+ buffer.copy(orientationSection, 0, space, orientationVec);
buffer.label("normal_dir");
buffer.scale(1.0);
- buffer.copy(dirSection);
return buffer;
-
} else if (0 == strcasecmp("traction", name)) {
assert(fields);
const topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
@@ -1690,10 +1682,8 @@
topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
_calcTractions(&buffer, dispT);
return buffer;
-
} else if (_friction->hasPropStateVar(name)) {
return _friction->getField(name);
-
} else if (_tractPerturbation && _tractPerturbation->hasParameter(name)) {
const topology::Field<topology::SubMesh>& param = _tractPerturbation->vertexField(name, fields);
if (param.vectorFieldType() == topology::FieldBase::VECTOR) {
@@ -1820,88 +1810,106 @@
assert(_fields);
const int spaceDim = _quadrature->spaceDim();
+ PetscErrorCode err;
// Get section information
- 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>& dispIncrSection =
- fields.get("dispIncr(t->t+dt)").section();
- 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);
- scalar_array dispRelVertex(spaceDim);
- const ALE::Obj<RealSection>& dispRelSection =
- _fields->get("relative disp").section();
- assert(!dispRelSection.isNull());
+ PetscSection dispRelSection = _fields->get("relative disp").petscSection();
+ Vec dispRelVec = _fields->get("relative disp").localVector();
+ PetscScalar *dispRelArray;
+ assert(dispRelSection);assert(dispRelVec);
- const ALE::Obj<RealSection>& velocitySection =
- fields.get("velocity(t)").section();
- assert(!velocitySection.isNull());
+ PetscSection velocitySection = fields.get("velocity(t)").petscSection();
+ Vec velocityVec = fields.get("velocity(t)").localVector();
+ PetscScalar *velocityArray;
+ assert(velocitySection);assert(velocityVec);
- scalar_array velRelVertex(spaceDim);
- const ALE::Obj<RealSection>& velRelSection =
- _fields->get("relative velocity").section();
- assert(!velRelSection.isNull());
+ PetscSection velRelSection = _fields->get("relative velocity").petscSection();
+ Vec velRelVec = _fields->get("relative velocity").localVector();
+ PetscScalar *velRelArray;
+ assert(velRelSection);assert(velRelVec);
const int numVertices = _cohesiveVertices.size();
+ err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(velRelVec, &velRelArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(velocityVec, &velocityArray);CHECK_PETSC_ERROR(err);
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_fault = _cohesiveVertices[iVertex].fault;
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
// Get displacement values
- 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 == dispIncrSection->getFiberDimension(v_negative));
- const PylithScalar* dispIncrVertexN =
- dispIncrSection->restrictPoint(v_negative);
- assert(dispIncrVertexN);
+ err = PetscSectionGetDof(dispTSection, v_positive, &dtpdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTSection, v_positive, &dtpoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dtpdof);
+ PetscInt dindof, dinoff;
- assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
- const PylithScalar* dispIncrVertexP =
- dispIncrSection->restrictPoint(v_positive);
- assert(dispIncrVertexP);
+ 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;
- // Compute relative displacememt
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- const PylithScalar value =
- dispTVertexP[iDim] + dispIncrVertexP[iDim]
- - dispTVertexN[iDim] - dispIncrVertexN[iDim];
- dispRelVertex[iDim] = fabs(value) > _zeroTolerance ? value : 0.0;
- } // for
+ err = PetscSectionGetDof(dispTIncrSection, v_positive, &dipdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTIncrSection, v_positive, &dipoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dipdof);
// Update relative displacement field.
- assert(dispRelVertex.size() ==
- dispRelSection->getFiberDimension(v_fault));
- dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
+ PetscInt drdof, droff;
+ err = PetscSectionGetDof(dispRelSection, v_fault, &drdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispRelSection, v_fault, &droff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == drdof);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ const PylithScalar value = dispTArray[dtpoff+d] + dispTIncrArray[dipoff+d] - dispTArray[dtnoff+d] - dispTIncrArray[dinoff+d];
+ dispRelArray[droff+d] = fabs(value) > _zeroTolerance ? value : 0.0;
+ } // for
+
// Get velocity values
- assert(spaceDim == velocitySection->getFiberDimension(v_negative));
- const PylithScalar* velocityVertexN = velocitySection->restrictPoint(v_negative);
- assert(velocityVertexN);
+ PetscInt vndof, vnoff;
- assert(spaceDim == velocitySection->getFiberDimension(v_positive));
- const PylithScalar* velocityVertexP = velocitySection->restrictPoint(v_positive);
- assert(velocityVertexP);
+ err = PetscSectionGetDof(velocitySection, v_negative, &vndof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(velocitySection, v_negative, &vnoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == vndof);
+ PetscInt vpdof, vpoff;
- // Compute relative velocity
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- const PylithScalar value = velocityVertexP[iDim] - velocityVertexN[iDim];
- velRelVertex[iDim] = fabs(value) > _zeroTolerance ? value : 0.0;
- } // for
+ err = PetscSectionGetDof(velocitySection, v_positive, &vpdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(velocitySection, v_positive, &vpoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == vpdof);
// Update relative velocity field.
- assert(velRelVertex.size() ==
- velRelSection->getFiberDimension(v_fault));
- velRelSection->updatePoint(v_fault, &velRelVertex[0]);
+ PetscInt vrdof, vroff;
+
+ err = PetscSectionGetDof(velRelSection, v_fault, &vrdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(velRelSection, v_fault, &vroff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == vrdof);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ const PylithScalar value = velocityArray[vpoff+d] - velocityArray[vnoff+d];
+ velRelArray[vroff+d] = fabs(value) > _zeroTolerance ? value : 0.0;
+ } // for
} // for
+ err = VecRestoreArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(velRelVec, &velRelArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(velocityVec, &velocityArray);CHECK_PETSC_ERROR(err);
PetscLogFlops(numVertices*spaceDim*spaceDim*4);
} // _updateRelMotion
@@ -1919,42 +1927,34 @@
// Setup fields involved in sensitivity solve.
if (!_fields->hasField("sensitivity solution")) {
_fields->add("sensitivity solution", "sensitivity_soln");
- topology::Field<topology::SubMesh>& solution =
- _fields->get("sensitivity solution");
- const topology::Field<topology::SubMesh>& dispRel =
- _fields->get("relative disp");
+ topology::Field<topology::SubMesh>& solution = _fields->get("sensitivity solution");
+ const topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
solution.cloneSection(dispRel);
solution.createScatter(solution.mesh());
} // if
- const topology::Field<topology::SubMesh>& solution =
- _fields->get("sensitivity solution");
+ const topology::Field<topology::SubMesh>& solution = _fields->get("sensitivity solution");
if (!_fields->hasField("sensitivity residual")) {
_fields->add("sensitivity residual", "sensitivity_residual");
- topology::Field<topology::SubMesh>& residual =
- _fields->get("sensitivity residual");
+ topology::Field<topology::SubMesh>& residual = _fields->get("sensitivity residual");
residual.cloneSection(solution);
residual.createScatter(solution.mesh());
} // if
if (!_fields->hasField("sensitivity relative disp")) {
_fields->add("sensitivity relative disp", "sensitivity_relative_disp");
- topology::Field<topology::SubMesh>& dispRel =
- _fields->get("sensitivity relative disp");
+ topology::Field<topology::SubMesh>& dispRel = _fields->get("sensitivity relative disp");
dispRel.cloneSection(solution);
} // if
- topology::Field<topology::SubMesh>& dispRel =
- _fields->get("sensitivity relative disp");
+ topology::Field<topology::SubMesh>& dispRel = _fields->get("sensitivity relative disp");
dispRel.zero();
if (!_fields->hasField("sensitivity dLagrange")) {
_fields->add("sensitivity dLagrange", "sensitivity_dlagrange");
- topology::Field<topology::SubMesh>& dLagrange =
- _fields->get("sensitivity dLagrange");
+ topology::Field<topology::SubMesh>& dLagrange = _fields->get("sensitivity dLagrange");
dLagrange.cloneSection(solution);
} // if
- topology::Field<topology::SubMesh>& dLagrange =
- _fields->get("sensitivity dLagrange");
+ topology::Field<topology::SubMesh>& dLagrange = _fields->get("sensitivity dLagrange");
dLagrange.zero();
// Setup Jacobian sparse matrix for sensitivity solve.
@@ -1966,26 +1966,24 @@
// Setup PETSc KSP linear solver.
if (0 == _ksp) {
PetscErrorCode err = 0;
- err = KSPCreate(_faultMesh->comm(), &_ksp); CHECK_PETSC_ERROR(err);
- err = KSPSetInitialGuessNonzero(_ksp, PETSC_FALSE); CHECK_PETSC_ERROR(err);
+ err = KSPCreate(_faultMesh->comm(), &_ksp);CHECK_PETSC_ERROR(err);
+ err = KSPSetInitialGuessNonzero(_ksp, PETSC_FALSE);CHECK_PETSC_ERROR(err);
PylithScalar rtol = 0.0;
PylithScalar atol = 0.0;
PylithScalar dtol = 0.0;
int maxIters = 0;
- err = KSPGetTolerances(_ksp, &rtol, &atol, &dtol, &maxIters);
- CHECK_PETSC_ERROR(err);
+ err = KSPGetTolerances(_ksp, &rtol, &atol, &dtol, &maxIters);CHECK_PETSC_ERROR(err);
rtol = 1.0e-3*_zeroTolerance;
atol = 1.0e-5*_zeroTolerance;
- err = KSPSetTolerances(_ksp, rtol, atol, dtol, maxIters);
- CHECK_PETSC_ERROR(err);
+ err = KSPSetTolerances(_ksp, rtol, atol, dtol, maxIters);CHECK_PETSC_ERROR(err);
PC pc;
- err = KSPGetPC(_ksp, &pc); CHECK_PETSC_ERROR(err);
- err = PCSetType(pc, PCJACOBI); CHECK_PETSC_ERROR(err);
- err = KSPSetType(_ksp, KSPGMRES); CHECK_PETSC_ERROR(err);
+ err = KSPGetPC(_ksp, &pc);CHECK_PETSC_ERROR(err);
+ err = PCSetType(pc, PCJACOBI);CHECK_PETSC_ERROR(err);
+ err = KSPSetType(_ksp, KSPGMRES);CHECK_PETSC_ERROR(err);
- err = KSPAppendOptionsPrefix(_ksp, "friction_");
- err = KSPSetFromOptions(_ksp); CHECK_PETSC_ERROR(err);
+ err = KSPAppendOptionsPrefix(_ksp, "friction_");CHECK_PETSC_ERROR(err);
+ err = KSPSetFromOptions(_ksp);CHECK_PETSC_ERROR(err);
} // if
} // _sensitivitySetup
@@ -2006,81 +2004,84 @@
// Get solution field
const topology::Field<topology::Mesh>& solutionDomain = fields.solution();
- const ALE::Obj<RealSection>& solutionDomainSection = solutionDomain.section();
- assert(!solutionDomainSection.isNull());
+ DM solutionDomainDM = solutionDomain.dmMesh();
+ PetscSection solutionDomainSection = solutionDomain.petscSection();
+ Vec solutionDomainVec = solutionDomain.localVector();
+ PetscSection solutionDomainGlobalSection;
+ PetscScalar *solutionDomainArray;
+ PetscErrorCode err;
+ assert(solutionDomainSection);assert(solutionDomainVec);
+ err = DMGetDefaultGlobalSection(solutionDomainDM, &solutionDomainGlobalSection);CHECK_PETSC_ERROR(err);
// Get cohesive cells
- const ALE::Obj<SieveMesh>& sieveMesh = fields.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cellsCohesive = sieveMesh->getLabelStratum("material-id", id());
- assert(!cellsCohesive.isNull());
- const SieveMesh::label_sequence::iterator cellsCohesiveBegin = cellsCohesive->begin();
- const SieveMesh::label_sequence::iterator cellsCohesiveEnd = cellsCohesive->end();
+ DM dmMesh = fields.mesh().dmMesh();
+ IS cellsCohesiveIS;
+ const PetscInt *cellsCohesive;
+ PetscInt numCohesiveCells, vStart, vEnd;
+ assert(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetStratumIS(dmMesh, "material-id", id(), &cellsCohesiveIS);CHECK_PETSC_ERROR(err);
+ err = ISGetLocalSize(cellsCohesiveIS, &numCohesiveCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellsCohesiveIS, &cellsCohesive);CHECK_PETSC_ERROR(err);
+
// Visitor for Jacobian matrix associated with domain.
scalar_array jacobianSubCell(submatrixSize);
const PetscMat jacobianDomainMatrix = jacobian.matrix();
assert(jacobianDomainMatrix);
- const ALE::Obj<SieveMesh::order_type>& globalOrderDomain =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solutionDomainSection);
- assert(!globalOrderDomain.isNull());
- const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
- assert(!sieve.isNull());
- const int closureSize = int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
- assert(closureSize >= 0);
- ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> ncV(*sieve, closureSize);
// Get fault Sieve mesh
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
+ DM faultDMMesh = _faultMesh->dmMesh();
+ assert(faultDMMesh);
// Get sensitivity solution field
- const ALE::Obj<RealSection>& solutionFaultSection =
- _fields->get("sensitivity solution").section();
- assert(!solutionFaultSection.isNull());
+ DM solutionFaultDM = _fields->get("sensitivity solution").dmMesh();
+ PetscSection solutionFaultSection = _fields->get("sensitivity solution").petscSection();
+ Vec solutionFaultVec = _fields->get("sensitivity solution").localVector();
+ PetscSection solutionFaultGlobalSection;
+ PetscScalar *solutionFaultArray;
+ assert(solutionFaultSection);assert(solutionFaultVec);
+ err = DMGetDefaultGlobalSection(solutionFaultDM, &solutionFaultGlobalSection);CHECK_PETSC_ERROR(err);
// Visitor for Jacobian matrix associated with fault.
assert(_jacobian);
const PetscMat jacobianFaultMatrix = _jacobian->matrix();
assert(jacobianFaultMatrix);
- const ALE::Obj<SieveSubMesh::order_type>& globalOrderFault =
- faultSieveMesh->getFactory()->getGlobalOrder(faultSieveMesh, "default", solutionFaultSection);
- assert(!globalOrderFault.isNull());
- // We would need to request unique points here if we had an interpolated mesh
- IndicesVisitor jacobianFaultVisitor(*solutionFaultSection,
- *globalOrderFault, closureSize*spaceDim);
const int iCone = (negativeSide) ? 0 : 1;
- const int numCohesiveCells = cellsCohesive->size();
- IS* cellsIS = (numCohesiveCells > 0) ? new IS[numCohesiveCells] : 0;
+ IS *cellsIS = (numCohesiveCells > 0) ? new IS[numCohesiveCells] : 0;
int_array indicesGlobal(subnrows);
int_array indicesLocal(numCohesiveCells*subnrows);
int_array indicesPerm(subnrows);
-
- PetscErrorCode err = 0;
- int iCohesiveCell = 0;
- for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter, ++iCohesiveCell) {
+ for(PetscInt c = 0; c < numCohesiveCells; ++c) {
// Get cone for cohesive cell
- ncV.clear();
- ALE::ISieveTraversal<SieveMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
- const int coneSize = ncV.getSize();
- assert(coneSize == 3*numBasis);
- const SieveMesh::point_type* cohesiveCone = ncV.getPoints();
- assert(cohesiveCone);
+ PetscInt *closure = PETSC_NULL;
+ PetscInt closureSize, q = 0;
+ err = DMComplexGetTransitiveClosure(dmMesh, cellsCohesive[c], PETSC_TRUE, &closureSize, &closure);CHECK_PETSC_ERROR(err);
+ // Filter out non-vertices
+ for(PetscInt p = 0; p < closureSize*2; p += 2) {
+ if ((closure[p] >= vStart) && (closure[p] < vEnd)) {
+ closure[q] = closure[p];
+ ++q;
+ }
+ }
+ closureSize = q;
+ assert(closureSize == 3*numBasis);
// Get indices
- for (int iBasis = 0; iBasis < numBasis; ++iBasis) {
+ for(int iBasis = 0; iBasis < numBasis; ++iBasis) {
// negative side of the fault: iCone=0
// positive side of the fault: iCone=1
- const int v_domain = cohesiveCone[iCone*numBasis+iBasis];
+ const int v_domain = closure[iCone*numBasis+iBasis];
+ PetscInt goff;
- for (int iDim=0, iB=iBasis*spaceDim; iDim < spaceDim; ++iDim) {
- indicesGlobal[iB+iDim] = globalOrderDomain->getIndex(v_domain) + iDim;
+ err = PetscSectionGetOffset(solutionDomainGlobalSection, v_domain, &goff);CHECK_PETSC_ERROR(err);
+ for(int iDim = 0, iB = iBasis*spaceDim, gind = goff < 0 ? -(goff+1) : goff; iDim < spaceDim; ++iDim) {
+ indicesGlobal[iB+iDim] = gind + iDim;
} // for
} // for
+ err = DMComplexRestoreTransitiveClosure(dmMesh, cellsCohesive[c], PETSC_TRUE, &closureSize, &closure);CHECK_PETSC_ERROR(err);
for (int i=0; i < subnrows; ++i) {
indicesPerm[i] = i;
@@ -2088,41 +2089,34 @@
err = PetscSortIntWithArray(indicesGlobal.size(), &indicesGlobal[0], &indicesPerm[0]);CHECK_PETSC_ERROR(err);
for (int i=0; i < subnrows; ++i) {
- indicesLocal[iCohesiveCell*subnrows+indicesPerm[i]] = i;
+ indicesLocal[c*subnrows+indicesPerm[i]] = i;
} // for
- cellsIS[iCohesiveCell] = PETSC_NULL;
- err = ISCreateGeneral(PETSC_COMM_SELF, indicesGlobal.size(), &indicesGlobal[0], PETSC_COPY_VALUES, &cellsIS[iCohesiveCell]);CHECK_PETSC_ERROR(err);
-
+ cellsIS[c] = PETSC_NULL;
+ err = ISCreateGeneral(PETSC_COMM_SELF, indicesGlobal.size(), &indicesGlobal[0], PETSC_COPY_VALUES, &cellsIS[c]);CHECK_PETSC_ERROR(err);
} // for
PetscMat* submatrices = 0;
err = MatGetSubMatrices(jacobianDomainMatrix, numCohesiveCells, cellsIS, cellsIS, MAT_INITIAL_MATRIX, &submatrices);CHECK_PETSC_ERROR(err);
- iCohesiveCell = 0;
- for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
- c_iter != cellsCohesiveEnd;
- ++c_iter, ++iCohesiveCell) {
+ for(PetscInt c = 0; c < numCohesiveCells; ++c) {
// Get values for submatrix associated with cohesive cell
jacobianSubCell = 0.0;
- err = MatGetValues(submatrices[iCohesiveCell],
- subnrows, &indicesLocal[iCohesiveCell*subnrows],
- subnrows, &indicesLocal[iCohesiveCell*subnrows],
- &jacobianSubCell[0]);CHECK_PETSC_ERROR_MSG(err, "Restrict from PETSc Mat failed.");
+ err = MatGetValues(submatrices[c], subnrows, &indicesLocal[c*subnrows], subnrows, &indicesLocal[c*subnrows],
+ &jacobianSubCell[0]);CHECK_PETSC_ERROR_MSG(err, "Restrict from PETSc Mat failed.");
// Insert cell contribution into PETSc Matrix
- const SieveMesh::point_type c_fault = _cohesiveToFault[*c_iter];
- jacobianFaultVisitor.clear();
- err = updateOperator(jacobianFaultMatrix, *faultSieveMesh->getSieve(),
- jacobianFaultVisitor, c_fault,
- &jacobianSubCell[0], INSERT_VALUES);
+ PetscInt c_fault = _cohesiveToFault[cellsCohesive[c]];
+ err = DMComplexMatSetClosure(faultDMMesh, solutionFaultSection, solutionFaultGlobalSection, jacobianFaultMatrix, c_fault, &jacobianSubCell[0], INSERT_VALUES);
CHECK_PETSC_ERROR_MSG(err, "Update to PETSc Mat failed.");
// Destory IS for cohesiveCell
- err = ISDestroy(&cellsIS[iCohesiveCell]);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellsIS[c]);CHECK_PETSC_ERROR(err);
} // for
err = MatDestroyMatrices(numCohesiveCells, &submatrices);CHECK_PETSC_ERROR(err);
delete[] cellsIS; cellsIS = 0;
+ err = ISRestoreIndices(cellsCohesiveIS, &cellsCohesive);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellsCohesiveIS);CHECK_PETSC_ERROR(err);
_jacobian->assemble("final_assembly");
@@ -2158,50 +2152,46 @@
scalar_array basisProducts(numBasis*numBasis);
// Get fault cell information
- const ALE::Obj<SieveMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- 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();
- const int numCells = cells->size();
+ DM faultDMMesh = _faultMesh->dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(faultDMMesh);
+ err = DMComplexGetHeightStratum(faultDMMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ const int numCells = cEnd-cStart;
+
// Get sections
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);
- scalar_array dLagrangeCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& dLagrangeSection =
- _fields->get("sensitivity dLagrange").section();
- assert(!dLagrangeSection.isNull());
- RestrictVisitor dLagrangeVisitor(*dLagrangeSection,
- dLagrangeCell.size(), &dLagrangeCell[0]);
+ PetscSection dLagrangeSection = _fields->get("sensitivity dLagrange").petscSection();
+ Vec dLagrangeVec = _fields->get("sensitivity dLagrange").localVector();
+ assert(dLagrangeSection);assert(dLagrangeVec);
scalar_array residualCell(numBasis*spaceDim);
- topology::Field<topology::SubMesh>& residual =
- _fields->get("sensitivity residual");
- const ALE::Obj<RealSection>& residualSection = residual.section();
- UpdateAddVisitor residualVisitor(*residualSection, &residualCell[0]);
-
+ topology::Field<topology::SubMesh>& residual = _fields->get("sensitivity residual");
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ assert(residualSection);assert(residualVec);
residual.zero();
// Loop over cells
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for(PetscInt c = cStart; c < cEnd; ++c) {
// Compute geometry
- 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);
+ 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);
// Restrict input fields to cell
- dLagrangeVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, dLagrangeVisitor);
+ const PetscScalar *dLagrangeArray = PETSC_NULL;
+ PetscInt dLagrangeSize;
+ err = DMComplexVecGetClosure(faultDMMesh, dLagrangeSection, dLagrangeVec, c, &dLagrangeSize, &dLagrangeArray);CHECK_PETSC_ERROR(err);
// Get cell geometry information that depends on cell
const scalar_array& basis = _quadrature->basis();
@@ -2213,31 +2203,30 @@
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const PylithScalar wt = quadWts[iQuad] * jacobianDet[iQuad];
- for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ for(int iBasis = 0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const PylithScalar valI = wt*basis[iQ+iBasis];
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ for(int jBasis = 0; jBasis < numBasis; ++jBasis) {
- basisProducts[iBasis*numBasis+jBasis] += valI*basis[iQ+jBasis];
- } // for
+ basisProducts[iBasis*numBasis+jBasis] += valI*basis[iQ+jBasis];
+ } // for
} // for
} // for
residualCell = 0.0;
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- for (int jBasis=0; jBasis < numBasis; ++jBasis) {
- const PylithScalar l = signFault * basisProducts[iBasis*numBasis+jBasis];
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- residualCell[iBasis*spaceDim+iDim] +=
- l * dLagrangeCell[jBasis*spaceDim+iDim];
- } // for
+ for(int iBasis = 0; iBasis < numBasis; ++iBasis) {
+ for(int jBasis = 0; jBasis < numBasis; ++jBasis) {
+ const PylithScalar l = signFault * basisProducts[iBasis*numBasis+jBasis];
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ residualCell[iBasis*spaceDim+d] += l * dLagrangeArray[jBasis*spaceDim+d];
+ } // for
} // for
} // for
+ err = DMComplexVecRestoreClosure(faultDMMesh, dLagrangeSection, dLagrangeVec, c, &dLagrangeSize, &dLagrangeArray);CHECK_PETSC_ERROR(err);
// Assemble cell contribution into field
- residualVisitor.clear();
- faultSieveMesh->updateClosure(*c_iter, residualVisitor);
+ err = DMComplexVecSetClosure(faultDMMesh, residualSection, residualVec, c, &residualCell[0], ADD_VALUES);CHECK_PETSC_ERROR(err);
} // for
} // _sensitivityReformResidual
@@ -2287,45 +2276,61 @@
assert(_quadrature);
const int spaceDim = _quadrature->spaceDim();
+ PetscErrorCode err;
scalar_array dispVertex(spaceDim);
- const ALE::Obj<RealSection>& solutionSection =
- _fields->get("sensitivity solution").section();
- assert(!solutionSection.isNull());
- const ALE::Obj<RealSection>& dispRelSection =
- _fields->get("sensitivity relative disp").section();
- assert(!dispRelSection.isNull());
- const ALE::Obj<RealSection>& dLagrangeTpdtSection =
- _fields->get("sensitivity dLagrange").section();
- assert(!dLagrangeTpdtSection.isNull());
+ PetscSection solutionSection = _fields->get("sensitivity solution").petscSection();
+ Vec solutionVec = _fields->get("sensitivity solution").localVector();
+ PetscScalar *solutionArray;
+ assert(solutionSection);assert(solutionVec);
+ PetscSection dispRelSection = _fields->get("sensitivity relative disp").petscSection();
+ Vec dispRelVec = _fields->get("sensitivity relative disp").localVector();
+ PetscScalar *dispRelArray;
+ assert(dispRelSection);assert(dispRelVec);
+ PetscSection dLagrangeTpdtSection = _fields->get("sensitivity dLagrange").petscSection();
+ Vec dLagrangeTpdtVec = _fields->get("sensitivity dLagrange").localVector();
+ PetscScalar *dLagrangeTpdtArray;
+ assert(dLagrangeTpdtSection);assert(dLagrangeTpdtVec);
const PylithScalar sign = (negativeSide) ? -1.0 : 1.0;
const int numVertices = _cohesiveVertices.size();
+ err = VecGetArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(solutionVec, &solutionArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dLagrangeTpdtVec, &dLagrangeTpdtArray);CHECK_PETSC_ERROR(err);
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_fault = _cohesiveVertices[iVertex].fault;
+ PetscInt sdof, soff;
- solutionSection->restrictPoint(v_fault, &dispVertex[0], dispVertex.size());
+ err = PetscSectionGetDof(solutionSection, v_fault, &sdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(solutionSection, v_fault, &soff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == sdof);
- // If no change in the Lagrange multiplier computed from friction
- // criterion, there are no updates, so continue.
- assert(spaceDim == dLagrangeTpdtSection->getFiberDimension(v_fault));
- const PylithScalar* dLagrangeTpdtVertex = dLagrangeTpdtSection->restrictPoint(v_fault);
- assert(dLagrangeTpdtVertex);
+ // If no change in the Lagrange multiplier computed from friction criterion, there are no updates, so continue.
+ PetscInt dldof, dloff;
+
+ err = PetscSectionGetDof(dLagrangeTpdtSection, v_fault, &dldof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dLagrangeTpdtSection, v_fault, &dloff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dldof);
PylithScalar dLagrangeTpdtVertexMag = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- dLagrangeTpdtVertexMag += dLagrangeTpdtVertex[iDim]*dLagrangeTpdtVertex[iDim];
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ dLagrangeTpdtVertexMag += dLagrangeTpdtArray[dloff+d]*dLagrangeTpdtArray[dloff+d];
} // for
- if (0.0 == dLagrangeTpdtVertexMag)
- continue;
+ if (0.0 == dLagrangeTpdtVertexMag) continue;
- // Update relative displacements associated with sensitivity solve
- // solution
- dispVertex *= sign;
+ // Update relative displacements associated with sensitivity solve solution
+ PetscInt drdof, droff;
- assert(dispVertex.size() == dispRelSection->getFiberDimension(v_fault));
- dispRelSection->updateAddPoint(v_fault, &dispVertex[0]);
+ err = PetscSectionGetDof(dispRelSection, v_fault, &drdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispRelSection, v_fault, &droff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == drdof);
+ for(PetscInt d = 0; d < drdof; ++d) {
+ dispRelArray[droff+d] = sign*solutionArray[soff+d];
+ }
} // for
+ err = VecRestoreArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(solutionVec, &solutionArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dLagrangeTpdtVec, &dLagrangeTpdtArray);CHECK_PETSC_ERROR(err);
} // _sensitivityUpdateSoln
@@ -2387,57 +2392,71 @@
Vec orientationVec = _fields->get("orientation").localVector();
PetscScalar *orientationArray;
- const ALE::Obj<RealSection>& dLagrangeTpdtSection = _fields->get("sensitivity dLagrange").section();
- assert(!dLagrangeTpdtSection.isNull());
+ PetscSection dLagrangeTpdtSection = _fields->get("sensitivity dLagrange").petscSection();
+ Vec dLagrangeTpdtVec = _fields->get("sensitivity dLagrange").localVector();
+ PetscScalar *dLagrangeTpdtArray;
+ assert(dLagrangeTpdtSection);assert(dLagrangeTpdtVec);
- const ALE::Obj<RealSection>& sensDispRelSection = _fields->get("sensitivity relative disp").section();
- assert(!sensDispRelSection.isNull());
+ PetscSection sensDispRelSection = _fields->get("sensitivity relative disp").petscSection();
+ Vec sensDispRelVec = _fields->get("sensitivity relative disp").localVector();
+ PetscScalar *sensDispRelArray;
- 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>& dispIncrSection = fields->get("dispIncr(t->t+dt)").section();
- assert(!dispIncrSection.isNull());
+ DM dispTIncrDM = fields->get("dispIncr(t->t+dt)").dmMesh();
+ PetscSection dispTIncrSection = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ PetscSection dispTIncrGlobalSection;
+ PetscScalar *dispTIncrArray;
+ assert(dispTIncrSection);assert(dispTIncrVec);
+ err = DMGetDefaultGlobalSection(dispTIncrDM, &dispTIncrGlobalSection);CHECK_PETSC_ERROR(err);
-
- // 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", dispIncrSection);
- assert(!globalOrder.isNull());
-
bool isOpening = false;
PylithScalar norm2 = 0.0;
int numVertices = _cohesiveVertices.size();
+ 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(sensDispRelVec, &sensDispRelArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dLagrangeTpdtVec, &dLagrangeTpdtArray);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(dispTIncrGlobalSection, v_lagrange, &goff);CHECK_PETSC_ERROR(err);
+ if (goff < 0) continue;
// Get displacement values
- 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;
+ err = PetscSectionGetDof(dispTSection, v_positive, &dtpdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTSection, v_positive, &dtpoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dtpdof);
+
// Get displacement increment values.
- assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
- const PylithScalar* dispIncrVertexN = dispIncrSection->restrictPoint(v_negative);
- assert(dispIncrVertexN);
+ PetscInt dindof, dinoff;
- assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
- const PylithScalar* dispIncrVertexP = dispIncrSection->restrictPoint(v_positive);
- assert(dispIncrVertexP);
+ 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;
+ err = PetscSectionGetDof(dispTIncrSection, v_positive, &dipdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTIncrSection, v_positive, &dipoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dipdof);
+
// Get orientation
PetscInt odof, ooff;
@@ -2446,43 +2465,43 @@
assert(spaceDim*spaceDim == odof);
// Get change in relative displacement from sensitivity solve.
- assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
- const PylithScalar* dDispRelVertex = sensDispRelSection->restrictPoint(v_fault);
- assert(dDispRelVertex);
+ PetscInt sdrdof, sdroff;
+ err = PetscSectionGetDof(sensDispRelSection, v_fault, &sdrdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(sensDispRelSection, v_fault, &sdroff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == sdrdof);
+
// Get Lagrange multiplier values from disp(t), and dispIncr(t->t+dt)
- assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
- const PylithScalar* lagrangeTVertex = dispTSection->restrictPoint(v_lagrange);
- assert(lagrangeTVertex);
+ PetscInt dtldof, dtloff;
- assert(spaceDim == dispIncrSection->getFiberDimension(v_lagrange));
- const PylithScalar* lagrangeTIncrVertex = dispIncrSection->restrictPoint(v_lagrange);
- assert(lagrangeTIncrVertex);
+ err = PetscSectionGetDof(dispTSection, v_lagrange, &dtldof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTSection, v_lagrange, &dtloff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dtldof);
+ PetscInt dildof, diloff;
- assert(spaceDim == dLagrangeTpdtSection->getFiberDimension(v_fault));
- const PylithScalar* dLagrangeTpdtVertex = dLagrangeTpdtSection->restrictPoint(v_fault);
- assert(dLagrangeTpdtVertex);
+ err = PetscSectionGetDof(dispTIncrSection, v_lagrange, &dildof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTIncrSection, v_lagrange, &diloff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dildof);
+ PetscInt sdldof, sdloff;
+ err = PetscSectionGetDof(dLagrangeTpdtSection, v_fault, &sdldof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dLagrangeTpdtSection, v_fault, &sdloff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == sdldof);
+
// Compute slip, slip rate, and traction at time t+dt as part of
// line search.
slipTpdtVertex = 0.0;
slipRateVertex = 0.0;
tractionTpdtVertex = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- for (int jDim=0; jDim < spaceDim; ++jDim) {
- slipTpdtVertex[iDim] += orientationArray[ooff+iDim*spaceDim+jDim] *
- (dispTVertexP[jDim] + dispIncrVertexP[jDim]
- - dispTVertexN[jDim] - dispIncrVertexN[jDim] +
- alpha*dDispRelVertex[jDim]);
- slipRateVertex[iDim] += orientationArray[ooff+iDim*spaceDim+jDim] *
- (dispIncrVertexP[jDim] - dispIncrVertexN[jDim] +
- alpha*dDispRelVertex[jDim]) / dt;
- tractionTpdtVertex[iDim] += orientationArray[ooff+iDim*spaceDim+jDim] *
- (lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim] +
- alpha*dLagrangeTpdtVertex[jDim]);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ for(PetscInt e = 0; e < spaceDim; ++e) {
+ slipTpdtVertex[d] += orientationArray[ooff+d*spaceDim+e] *
+ (dispTArray[dtpoff+e] + dispTIncrArray[dipoff+e] - dispTArray[dtnoff+e] - dispTIncrArray[dinoff+e] + alpha*sensDispRelArray[sdroff+e]);
+ slipRateVertex[d] += orientationArray[ooff+d*spaceDim+e] * (dispTIncrArray[dipoff+e] - dispTIncrArray[dinoff+e] + alpha*sensDispRelArray[sdroff+e]) / dt;
+ tractionTpdtVertex[d] += orientationArray[ooff+d*spaceDim+e] * (dispTArray[dtloff+e] + dispTIncrArray[diloff+e] + alpha*dLagrangeTpdtArray[sdloff+e]);
} // for
- if (fabs(slipRateVertex[iDim]) < _zeroTolerance) {
- slipRateVertex[iDim] = 0.0;
+ if (fabs(slipRateVertex[d]) < _zeroTolerance) {
+ slipRateVertex[d] = 0.0;
} // if
} // for
if (fabs(slipTpdtVertex[indexN]) < _zeroTolerance) {
@@ -2501,20 +2520,20 @@
// traction and slip to zero (should be appropriate if problem
// is nondimensionalized correctly).
if (fabs(slipTpdtVertex[indexN]) > fabs(tractionTpdtVertex[indexN])) {
- // fault opening is bigger, so force normal traction back to zero
- tractionTpdtVertex[indexN] = 0.0;
+ // fault opening is bigger, so force normal traction back to zero
+ tractionTpdtVertex[indexN] = 0.0;
} else {
- // traction is bigger, so force fault opening back to zero
- slipTpdtVertex[indexN] = 0.0;
+ // traction is bigger, so force fault opening back to zero
+ slipTpdtVertex[indexN] = 0.0;
} // if/else
} else if (slipTpdtVertex[indexN] > _zeroTolerance) {
- // Step b: Insure fault traction is zero when opening (if
+ // Step b: Ensure fault traction is zero when opening (if
// alpha=1 this should be enforced already, but will not be
// properly enforced when alpha < 1).
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- tractionTpdtVertex[iDim] = 0.0;
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ tractionTpdtVertex[d] = 0.0;
} // for
} else if (slipTpdtVertex[indexN] < 0.0) {
// Step c: Prevent interpenetration.
@@ -2537,11 +2556,9 @@
// friction.
tractionMisfitVertex = 0.0;
const bool iterating = true; // Iterating to get friction
- CALL_MEMBER_FN(*this,
- constrainSolnSpaceFn)(&tractionMisfitVertex, t,
- slipTpdtVertex, slipRateVertex,
- tractionTpdtVertex,
- iterating);
+ CALL_MEMBER_FN(*this, constrainSolnSpaceFn)(&tractionMisfitVertex, t,
+ slipTpdtVertex, slipRateVertex, tractionTpdtVertex,
+ iterating);
#if 0 // DEBUGGING
std::cout << "alpha: " << alpha
@@ -2560,28 +2577,29 @@
} // for
std::cout << ", dDispRel:";
for (int iDim=0; iDim < spaceDim; ++iDim) {
- std::cout << " " << dDispRelVertex[iDim];
+ std::cout << " " << sensDispRelArray[sdroff+iDim];
} // for
std::cout << std::endl;
#endif
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- norm2 += tractionMisfitVertex[iDim]*tractionMisfitVertex[iDim];
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ norm2 += tractionMisfitVertex[d]*tractionMisfitVertex[d];
} // for
} // for
+ 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(sensDispRelVec, &sensDispRelArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dLagrangeTpdtVec, &dLagrangeTpdtArray);CHECK_PETSC_ERROR(err);
if (isOpening && alpha < 1.0) {
norm2 = PYLITH_MAXFLOAT;
} // if
- PylithScalar norm2Total = 0.0;
- int numVerticesTotal = 0;
- if (sizeof(PylithScalar) == 8) {
- MPI_Allreduce(&norm2, &norm2Total, 1, MPI_DOUBLE, MPI_SUM, fields->mesh().comm());
- } else {
- MPI_Allreduce(&norm2, &norm2Total, 1, MPI_FLOAT, MPI_SUM, fields->mesh().comm());
- } // if/else
- MPI_Allreduce(&numVertices, &numVerticesTotal, 1, MPI_INT, MPI_SUM, fields->mesh().comm());
+ PetscScalar norm2Total = 0.0;
+ PetscInt numVerticesTotal = 0;
+ err = MPI_Allreduce(&norm2, &norm2Total, 1, MPIU_SCALAR, MPI_SUM, fields->mesh().comm());
+ err = MPI_Allreduce(&numVertices, &numVerticesTotal, 1, MPIU_INT, MPI_SUM, fields->mesh().comm());
assert(numVerticesTotal > 0);
return sqrt(norm2Total) / numVerticesTotal;
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveImpulses.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -207,53 +207,41 @@
return buffer;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection = _fields->get(
- "orientation").section();
- assert(!orientationSection.isNull());
- const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(
- 0);
- assert(!dirSection.isNull());
+ PetscSection orientationSection = _fields->get("orientation").petscSection();
+ Vec orientationVec = _fields->get("orientation").localVector();
+ assert(orientationSection);assert(orientationVec);
_allocateBufferVectorField();
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
- buffer.copy(dirSection);
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
+ buffer.copy(orientationSection, 0, 0, orientationVec);
buffer.label("strike_dir");
buffer.scale(1.0);
return buffer;
} else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection = _fields->get(
- "orientation").section();
- assert(!orientationSection.isNull());
- const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(
- 1);
+ PetscSection orientationSection = _fields->get("orientation").petscSection();
+ Vec orientationVec = _fields->get("orientation").localVector();
+ assert(orientationSection);assert(orientationVec);
_allocateBufferVectorField();
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
- buffer.copy(dirSection);
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
+ buffer.copy(orientationSection, 0, 1, orientationVec);
buffer.label("dip_dir");
buffer.scale(1.0);
return buffer;
} else if (0 == strcasecmp("normal_dir", name)) {
- const ALE::Obj<RealSection>& orientationSection = _fields->get(
- "orientation").section();
- assert(!orientationSection.isNull());
+ PetscSection orientationSection = _fields->get("orientation").petscSection();
+ Vec orientationVec = _fields->get("orientation").localVector();
+ assert(orientationSection);assert(orientationVec);
const int space = (0 == cohesiveDim) ? 0 : (1 == cohesiveDim) ? 1 : 2;
- const ALE::Obj<RealSection>& dirSection = orientationSection->getFibration(
- space);
- assert(!dirSection.isNull());
_allocateBufferVectorField();
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
- buffer.copy(dirSection);
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
+ buffer.copy(orientationSection, 0, space, orientationVec);
buffer.label("normal_dir");
buffer.scale(1.0);
return buffer;
} else if (0 == strcasecmp("impulse_amplitude", name)) {
- topology::Field<topology::SubMesh>& amplitude =
- _fields->get("impulse amplitude");
+ topology::Field<topology::SubMesh>& amplitude = _fields->get("impulse amplitude");
return amplitude;
} else if (0 == strcasecmp("area", name)) {
@@ -264,8 +252,7 @@
assert(fields);
const topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
_allocateBufferVectorField();
- topology::Field<topology::SubMesh>& buffer =
- _fields->get("buffer (vector)");
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (vector)");
_calcTractionsChange(&buffer, dispT);
return buffer;
@@ -300,31 +287,34 @@
const PylithScalar lengthScale = _normalizer->lengthScale();
const int spaceDim = _quadrature->spaceDim();
+ PetscErrorCode err;
// Create section to hold amplitudes of impulses.
_fields->add("impulse amplitude", "impulse_amplitude");
- topology::Field<topology::SubMesh>& amplitude =
- _fields->get("impulse amplitude");
+ topology::Field<topology::SubMesh>& amplitude = _fields->get("impulse amplitude");
topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
const int fiberDim = 1;
amplitude.newSection(dispRel, fiberDim);
amplitude.allocate();
amplitude.scale(lengthScale);
- PylithScalar amplitudeVertex;
- const ALE::Obj<RealSection>& amplitudeSection = amplitude.section();
- assert(!amplitudeSection.isNull());
+ PetscSection amplitudeSection = amplitude.petscSection();
+ Vec amplitudeVec = amplitude.localVector();
+ PetscScalar *amplitudeArray;
+ assert(amplitudeSection);assert(amplitudeVec);
const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
assert(cs);
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
+ DM faultDMMesh = _faultMesh->dmMesh();
+ assert(faultDMMesh);
scalar_array coordsVertex(spaceDim);
- const ALE::Obj<RealSection>& coordsSection =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordsSection.isNull());
+ PetscSection coordSection;
+ Vec coordVec;
+ PetscScalar *coords;
+ err = DMComplexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
assert(_dbImpulseAmp);
_dbImpulseAmp->open();
@@ -334,39 +324,48 @@
std::map<int, int> pointOrder;
int count = 0;
const int numVertices = _cohesiveVertices.size();
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_fault = _cohesiveVertices[iVertex].fault;
+ PetscInt cdof, coff;
- coordsSection->restrictPoint(v_fault, &coordsVertex[0], coordsVertex.size());
- _normalizer->dimensionalize(&coordsVertex[0], coordsVertex.size(),
- lengthScale);
+ err = PetscSectionGetDof(coordSection, v_fault, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(coordSection, v_fault, &coff);CHECK_PETSC_ERROR(err);
+ assert(cdof == spaceDim);
- amplitudeVertex = 0.0;
- int err = _dbImpulseAmp->query(&litudeVertex, 1,
- &coordsVertex[0], coordsVertex.size(), cs);
+ for(PetscInt d = 0; d < cdof; ++d) coordsVertex[d] = coords[coff+d];
+ _normalizer->dimensionalize(&coordsVertex[0], coordsVertex.size(), lengthScale);
+ PetscInt adof, aoff;
+
+ err = PetscSectionGetDof(amplitudeSection, v_fault, &adof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(amplitudeSection, v_fault, &aoff);CHECK_PETSC_ERROR(err);
+ assert(adof == fiberDim);
+
+ amplitudeArray[aoff] = 0.0;
+ int err = _dbImpulseAmp->query(&litudeArray[aoff], 1, &coordsVertex[0], coordsVertex.size(), cs);
if (err) {
std::ostringstream msg;
msg << "Could not find amplitude for Green's function impulses at \n" << "(";
for (int i = 0; i < spaceDim; ++i)
- msg << " " << coordsVertex[i];
+ msg << " " << coordsVertex[i];
msg << ") in fault " << label() << "\n"
- << "using spatial database '" << _dbImpulseAmp->label() << "'.";
+ << "using spatial database '" << _dbImpulseAmp->label() << "'.";
throw std::runtime_error(msg.str());
} // if
- if (fabs(amplitudeVertex) < _threshold) {
- amplitudeVertex = 0.0;
+ if (fabs(amplitudeArray[aoff]) < _threshold) {
+ amplitudeArray[aoff] = 0.0;
} // if
- _normalizer->nondimensionalize(&litudeVertex, 1, lengthScale);
+ _normalizer->nondimensionalize(&litudeArray[aoff], 1, lengthScale);
- if (fabs(amplitudeVertex) > 0.0) {
+ if (fabs(amplitudeArray[aoff]) > 0.0) {
pointOrder[iVertex] = count;
++count;
} // if
-
- assert(1 == amplitudeSection->getFiberDimension(v_fault));
- amplitudeSection->updatePoint(v_fault, &litudeVertex);
} // for
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
// Close properties database
_dbImpulseAmp->close();
@@ -385,16 +384,18 @@
// Order of impulses is set by processor rank and order of points in
// mesh, using only those points with nonzero amplitudes.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
+ DM faultDMMesh = _faultMesh->dmMesh();
+ PetscErrorCode err;
+ assert(faultDMMesh);
// Gather number of points on each processor.
const int numImpulsesLocal = pointOrder.size();
- const int commSize = faultSieveMesh->commSize();
- const int commRank = faultSieveMesh->commRank();
+ MPI_Comm comm = ((PetscObject) faultDMMesh)->comm;
+ PetscMPIInt commSize, commRank;
+ err = MPI_Comm_size(comm, &commSize);CHECK_PETSC_ERROR(err);
+ err = MPI_Comm_rank(comm, &commRank);CHECK_PETSC_ERROR(err);
int_array numImpulsesAll(commSize);
- MPI_Comm comm = faultSieveMesh->comm();
- MPI_Allgather((void*)&numImpulsesLocal, 1, MPI_INT, (void*)&numImpulsesAll[0], commSize, MPI_INT, comm);
+ err = MPI_Allgather((void *) &numImpulsesLocal, 1, MPI_INT, (void *) &numImpulsesAll[0], commSize, MPI_INT, comm);CHECK_PETSC_ERROR(err);
int localOffset = 0;
for (int i=0; i < commRank; ++i) {
@@ -449,34 +450,45 @@
const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
assert(cs);
const int spaceDim = cs->spaceDim();
+ PetscErrorCode err;
- const ALE::Obj<RealSection>& amplitudeSection = _fields->get("impulse amplitude").section();
- assert(!amplitudeSection.isNull());
+ PetscSection amplitudeSection = _fields->get("impulse amplitude").petscSection();
+ Vec amplitudeVec = _fields->get("impulse amplitude").localVector();
+ PetscScalar *amplitudeArray;
+ assert(amplitudeSection);assert(amplitudeVec);
- const ALE::Obj<RealSection>& dispRelSection = dispRel.section();
- assert(!dispRelSection.isNull());
+ PetscSection dispRelSection = dispRel.petscSection();
+ Vec dispRelVec = dispRel.localVector();
+ PetscScalar *dispRelArray;
+ assert(dispRelSection);assert(dispRelVec);
- scalar_array dispRelVertex(spaceDim);
- dispRelVertex = 0.0;
-
const srcs_type::const_iterator& impulseInfo = _impulsePoints.find(impulse);
+ err = VecGetArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
if (impulseInfo != _impulsePoints.end()) {
const int iVertex = impulseInfo->second.indexCohesive;
const int v_fault = _cohesiveVertices[iVertex].fault;
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
// Get amplitude of slip impulse
- assert(1 == amplitudeSection->getFiberDimension(v_fault));
- const PylithScalar* amplitudeVertex = amplitudeSection->restrictPoint(v_fault);
- assert(amplitudeVertex);
+ PetscInt adof, aoff;
+ err = PetscSectionGetDof(amplitudeSection, v_fault, &adof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(amplitudeSection, v_fault, &aoff);CHECK_PETSC_ERROR(err);
+ assert(adof == 1);
+ PetscInt drdof, droff;
+
+ err = PetscSectionGetDof(dispRelSection, v_fault, &drdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispRelSection, v_fault, &droff);CHECK_PETSC_ERROR(err);
+ assert(drdof == spaceDim);
+ for(PetscInt d = 0; d < drdof; ++d) {dispRelArray[droff+d] = 0.0;}
const int indexDOF = impulseInfo->second.indexDOF;
+
assert(indexDOF >= 0 && indexDOF < spaceDim);
- dispRelVertex[indexDOF] = amplitudeVertex[0];
-
- assert(dispRelVertex.size() == dispRelSection->getFiberDimension(v_fault));
- dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
+ dispRelArray[droff+indexDOF] = amplitudeArray[aoff];
} // if
+ err = VecRestoreArray(amplitudeVec, &litudeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispRelVec, &dispRelArray);CHECK_PETSC_ERROR(err);
#if 0 // DEBUGGING
std::cout << "impulse: " << impulse << std::endl;
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -1687,6 +1687,7 @@
PetscSection areaSection = area.petscSection();
Vec areaVec = area.localVector();
PetscScalar *areaArray;
+ assert(areaSection);assert(areaVec);
logger.stagePop();
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -130,7 +130,7 @@
_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->add("change time", "change_start_time", topology::FieldBase::VERTICES_FIELD, 1);
_parameters->get("change time").vectorFieldType(topology::FieldBase::SCALAR);
_parameters->get("change time").scale(timeScale);
_parameters->get("change time").allocate();
@@ -485,8 +485,8 @@
err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
assert(coordSection);assert(coordVec);
- PetscSection valueSection = _parameters->get("value").petscSection();
- Vec valueVec = _parameters->get("value").localVector();
+ PetscSection valueSection = _parameters->get(name).petscSection();
+ Vec valueVec = _parameters->get(name).localVector();
PetscScalar *tractionsArray;
assert(valueSection);assert(valueVec);
@@ -523,7 +523,7 @@
err = PetscSectionGetDof(valueSection, v, &vdof);CHECK_PETSC_ERROR(err);
err = PetscSectionGetOffset(valueSection, v, &voff);CHECK_PETSC_ERROR(err);
- assert(vdof == spaceDim);
+ assert(vdof == querySize);
for(PetscInt d = 0; d < vdof; ++d)
tractionsArray[voff+d] = valuesVertex[d];
} // for
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -78,11 +78,6 @@
const char* label,
const int labelId)
{ // open
- typedef typename mesh_type::SieveMesh SieveMesh;
- typedef typename mesh_type::SieveMesh::label_sequence label_sequence;
- typedef typename mesh_type::SieveMesh::numbering_type numbering_type;
- typedef typename mesh_type::SieveMesh::sieve_type sieve_type;
-
DataWriter<mesh_type, field_type>::open(mesh, numTimeSteps, label, labelId);
try {
@@ -92,13 +87,13 @@
const std::string& filename = _hdf5Filename();
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- DM dmMesh = mesh.dmMesh();
- assert(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ assert(dmMesh);
_timesteps.clear();
_tstampIndex = 0;
- const int commRank = sieveMesh->commRank();
+ PetscMPIInt commRank;
+ err = MPI_Comm_rank(mesh.comm(), &commRank);CHECK_PETSC_ERROR(err);
const int localSize = (!commRank) ? 1 : 0;
err = VecCreateMPI(mesh.comm(), localSize, 1, &_tstamp);
CHECK_PETSC_ERROR(err);
@@ -106,67 +101,78 @@
err = VecSetBlockSize(_tstamp, 1); CHECK_PETSC_ERROR(err);
err = PetscObjectSetName((PetscObject) _tstamp, "time");
- err = PetscViewerHDF5Open(mesh.comm(), filename.c_str(), FILE_MODE_WRITE,
- &_viewer);
+ err = PetscViewerHDF5Open(mesh.comm(), filename.c_str(), FILE_MODE_WRITE, &_viewer);
CHECK_PETSC_ERROR(err);
- ALE::Obj<numbering_type> vNumbering =
- sieveMesh->hasLabel("censored depth") ?
- sieveMesh->getFactory()->getNumbering(sieveMesh, "censored depth", 0) :
- sieveMesh->getFactory()->getNumbering(sieveMesh, 0);
- assert(!vNumbering.isNull());
- //vNumbering->view("VERTEX NUMBERING");
- if (dmMesh) {
- PetscSection coordSection;
- Vec coordinates;
- PetscReal lengthScale;
- PetscInt vStart, vEnd, vMax, verticesSize, dim, dimLocal = 0;
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ assert(cs);
- const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- assert(cs);
- err = DMComplexGetScale(dmMesh, PETSC_UNIT_LENGTH, &lengthScale);CHECK_PETSC_ERROR(err);
- err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dmMesh, &coordinates);CHECK_PETSC_ERROR(err);
- err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
- err = DMComplexGetVTKBounds(dmMesh, PETSC_NULL, &vMax);CHECK_PETSC_ERROR(err);
- if (vMax >= 0) {vEnd = PetscMin(vEnd, vMax);}
- for(PetscInt vertex = vStart; vertex < vEnd; ++vertex) {
- err = PetscSectionGetDof(coordSection, vertex, &dimLocal);CHECK_PETSC_ERROR(err);
- if (dimLocal) break;
- }
- err = MPI_Allreduce(&dimLocal, &dim, 1, MPIU_INT, MPI_MAX, sieveMesh->comm());CHECK_PETSC_ERROR(err);
- verticesSize = vEnd - vStart;
+#if 1
+ const char *context = DataWriter<mesh_type, field_type>::_context.c_str();
+ DM dmCoord;
+ PetscReal lengthScale;
- PetscVec coordVec;
- PetscScalar *coords, *c;
+ err = DMComplexGetScale(dmMesh, PETSC_UNIT_LENGTH, &lengthScale);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinateDM(dmMesh, &dmCoord);CHECK_PETSC_ERROR(err);
+ topology::Field<mesh_type> coordinates(mesh, dmCoord, topology::FieldBase::Metadata());
+ coordinates.label("vertices");
+ coordinates.createScatterWithBC(mesh, PETSC_NULL, context);
+ coordinates.scatterSectionToVector(context);
+ PetscVec coordVec = coordinates.vector(context);
+ assert(coordVec);
+ err = VecScale(coordVec, lengthScale);CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PushGroup(_viewer, "/geometry");CHECK_PETSC_ERROR(err);
+ err = VecView(coordVec, _viewer);CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
+#else
+ PetscSection coordSection;
+ Vec coordinates;
+ PetscReal lengthScale;
+ PetscInt vStart, vEnd, vMax, verticesSize, dim, dimLocal = 0;
- err = VecCreate(sieveMesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
- err = VecSetSizes(coordVec, verticesSize*dim, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
- err = VecSetBlockSize(coordVec, dim);CHECK_PETSC_ERROR(err);
- err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
- err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
- err = VecGetArray(coordinates, &c);CHECK_PETSC_ERROR(err);
- for(PetscInt v = 0; v < vEnd - vStart; ++v) {
- for(PetscInt d = 0; d < dim; ++d) {
+ /* TODO Get rid of this and use the createScatterWithBC(numbering) code */
+ err = DMComplexGetScale(dmMesh, PETSC_UNIT_LENGTH, &lengthScale);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordinates);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetVTKBounds(dmMesh, PETSC_NULL, &vMax);CHECK_PETSC_ERROR(err);
+ if (vMax >= 0) {vEnd = PetscMin(vEnd, vMax);}
+ for(PetscInt vertex = vStart; vertex < vEnd; ++vertex) {
+ err = PetscSectionGetDof(coordSection, vertex, &dimLocal);CHECK_PETSC_ERROR(err);
+ if (dimLocal) break;
+ }
+ err = MPI_Allreduce(&dimLocal, &dim, 1, MPIU_INT, MPI_MAX, mesh.comm());CHECK_PETSC_ERROR(err);
+ verticesSize = vEnd - vStart;
+
+ PetscVec coordVec;
+ PetscScalar *coords, *c;
+
+ err = VecCreate(mesh.comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, verticesSize*dim, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetBlockSize(coordVec, dim);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordinates, &c);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = 0; v < vEnd - vStart; ++v) {
+ for(PetscInt d = 0; d < dim; ++d) {
coords[v*dim+d] = c[v*dim+d];
- }
}
- err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
- err = VecRestoreArray(coordinates, &c);CHECK_PETSC_ERROR(err);
- err = VecScale(coordVec, lengthScale);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) coordVec, "vertices");CHECK_PETSC_ERROR(err);
- err = PetscViewerHDF5PushGroup(_viewer, "/geometry");CHECK_PETSC_ERROR(err);
- err = VecView(coordVec, _viewer);CHECK_PETSC_ERROR(err);
- err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
- err = VecDestroy(&coordVec);CHECK_PETSC_ERROR(err);
- } else {
+ }
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordinates, &c);CHECK_PETSC_ERROR(err);
+ err = VecScale(coordVec, lengthScale);CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) coordVec, "vertices");CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PushGroup(_viewer, "/geometry");CHECK_PETSC_ERROR(err);
+ err = VecView(coordVec, _viewer);CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
+ err = VecDestroy(&coordVec);CHECK_PETSC_ERROR(err);
+#endif
+#if 0
const ALE::Obj<typename mesh_type::RealSection>& coordinatesSection =
sieveMesh->hasRealSection("coordinates_dimensioned") ?
sieveMesh->getRealSection("coordinates_dimensioned") :
sieveMesh->getRealSection("coordinates");
assert(!coordinatesSection.isNull());
- const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
- assert(cs);
topology::FieldBase::Metadata metadata;
// :KLUDGE: We would like to use field_type for the coordinates
// field. However, the mesh coordinates are Field<mesh_type> and
@@ -183,62 +189,61 @@
err = PetscViewerHDF5PushGroup(_viewer, "/geometry");CHECK_PETSC_ERROR(err);
err = VecView(coordinatesVector, _viewer);CHECK_PETSC_ERROR(err);
err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
- }
+#endif
- if (dmMesh) {
- PetscInt *cones;
- PetscSection coneSection;
- PetscInt vStart, vEnd, cStart, cEnd, cMax, dof, conesSize, numCorners, numCornersLocal = 0;
+ PetscInt *cones;
+ PetscSection coneSection;
+ PetscInt vStart, vEnd, cStart, cEnd, cMax, dof, conesSize, numCorners, numCornersLocal = 0;
- err = DMComplexGetConeSection(dmMesh, &coneSection);CHECK_PETSC_ERROR(err);
- err = DMComplexGetCones(dmMesh, &cones);CHECK_PETSC_ERROR(err);
- err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
- err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
- err = DMComplexGetVTKBounds(dmMesh, &cMax, PETSC_NULL);CHECK_PETSC_ERROR(err);
- if (cMax >= 0) {cEnd = PetscMin(cEnd, cMax);}
+ err = DMComplexGetConeSection(dmMesh, &coneSection);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCones(dmMesh, &cones);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetVTKBounds(dmMesh, &cMax, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ if (cMax >= 0) {cEnd = PetscMin(cEnd, cMax);}
+ for(PetscInt cell = cStart; cell < cEnd; ++cell) {
+ err = DMComplexGetConeSize(dmMesh, cell, &numCornersLocal);CHECK_PETSC_ERROR(err);
+ if (numCornersLocal) break;
+ }
+ err = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, mesh.comm());CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(coneSection, cEnd, &conesSize);CHECK_PETSC_ERROR(err);
+ if (label) {
+ conesSize = 0;
for(PetscInt cell = cStart; cell < cEnd; ++cell) {
- err = DMComplexGetConeSize(dmMesh, cell, &numCornersLocal);CHECK_PETSC_ERROR(err);
- if (numCornersLocal) break;
- }
- err = MPI_Allreduce(&numCornersLocal, &numCorners, 1, MPIU_INT, MPI_MAX, sieveMesh->comm());CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(coneSection, cEnd, &conesSize);CHECK_PETSC_ERROR(err);
- if (label) {
- conesSize = 0;
- for(PetscInt cell = cStart; cell < cEnd; ++cell) {
- PetscInt value;
+ PetscInt value;
- err = DMComplexGetLabelValue(dmMesh, label, cell, &value);CHECK_PETSC_ERROR(err);
- if (value == labelId) ++conesSize;
- }
- conesSize *= numCorners;
+ err = DMComplexGetLabelValue(dmMesh, label, cell, &value);CHECK_PETSC_ERROR(err);
+ if (value == labelId) ++conesSize;
}
+ conesSize *= numCorners;
+ }
- PetscVec cellVec;
- PetscScalar *vertices;
+ PetscVec cellVec;
+ PetscScalar *vertices;
- err = VecCreate(sieveMesh->comm(), &cellVec);CHECK_PETSC_ERROR(err);
- err = VecSetSizes(cellVec, conesSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
- err = VecSetBlockSize(cellVec, numCorners);CHECK_PETSC_ERROR(err);
- err = VecSetFromOptions(cellVec);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) cellVec, "cells");CHECK_PETSC_ERROR(err);
- err = VecGetArray(cellVec, &vertices);CHECK_PETSC_ERROR(err);
- for(PetscInt cell = cStart, v = 0; cell < cEnd; ++cell) {
- if (label) {
- PetscInt value;
+ err = VecCreate(mesh.comm(), &cellVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(cellVec, conesSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetBlockSize(cellVec, numCorners);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(cellVec);CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) cellVec, "cells");CHECK_PETSC_ERROR(err);
+ err = VecGetArray(cellVec, &vertices);CHECK_PETSC_ERROR(err);
+ for(PetscInt cell = cStart, v = 0; cell < cEnd; ++cell) {
+ if (label) {
+ PetscInt value;
- err = DMComplexGetLabelValue(dmMesh, label, cell, &value);CHECK_PETSC_ERROR(err);
- if (value != labelId) continue;
- }
- for(PetscInt corner = 0; corner < numCorners; ++corner, ++v) {
- vertices[v] = cones[cell*numCorners+corner] - vStart;
- }
+ err = DMComplexGetLabelValue(dmMesh, label, cell, &value);CHECK_PETSC_ERROR(err);
+ if (value != labelId) continue;
}
- err = VecRestoreArray(cellVec, &vertices);CHECK_PETSC_ERROR(err);
- err = PetscViewerHDF5PushGroup(_viewer, "/topology");CHECK_PETSC_ERROR(err);
- err = VecView(cellVec, _viewer);CHECK_PETSC_ERROR(err);
- err = PetscViewerHDF5PopGroup(_viewer);CHECK_PETSC_ERROR(err);
- err = VecDestroy(&cellVec);CHECK_PETSC_ERROR(err);
- } else {
+ for(PetscInt corner = 0; corner < numCorners; ++corner, ++v) {
+ vertices[v] = cones[cell*numCorners+corner] - vStart;
+ }
+ }
+ err = VecRestoreArray(cellVec, &vertices);CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PushGroup(_viewer, "/topology");CHECK_PETSC_ERROR(err);
+ err = VecView(cellVec, _viewer);CHECK_PETSC_ERROR(err);
+ err = PetscViewerHDF5PopGroup(_viewer);CHECK_PETSC_ERROR(err);
+ err = VecDestroy(&cellVec);CHECK_PETSC_ERROR(err);
+#if 0
// Account for censored cells
int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
int cellDepth = 0;
@@ -307,7 +312,7 @@
err = PetscViewerHDF5PopGroup(_viewer);CHECK_PETSC_ERROR(err);
err = VecDestroy(&elemVec);CHECK_PETSC_ERROR(err);
delete[] tmpVertices; tmpVertices = 0;
- }
+#endif
hid_t h5 = -1;
err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
@@ -364,29 +369,14 @@
const mesh_type& mesh)
{ // writeVertexField
typedef typename mesh_type::SieveMesh::numbering_type numbering_type;
+ PetscErrorCode err;
assert(_viewer);
try {
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- ALE::Obj<numbering_type> vNumbering =
- sieveMesh->hasLabel("censored depth") ?
- sieveMesh->getFactory()->getNumbering(sieveMesh, "censored depth", 0) :
- sieveMesh->getFactory()->getNumbering(sieveMesh, 0);
- assert(!vNumbering.isNull());
-
-#if 0
- std::cout << "WRITE VERTEX FIELD" << std::endl;
- mesh.view("MESH");
- field.view("FIELD");;
- vNumbering->view("NUMBERING");
- std::cout << std::endl;
-#endif
-
- field.createScatterWithBC(mesh, vNumbering, context);
+ field.createScatterWithBC(mesh, PETSC_NULL, context);
field.scatterSectionToVector(context);
PetscVec vector = field.vector(context);
assert(vector);
@@ -397,15 +387,11 @@
_timesteps[field.label()] += 1;
const int istep = _timesteps[field.label()];
// Add time stamp to "/time" if necessary.
- const int commRank = sieveMesh->commRank();
+ PetscMPIInt commRank;
+ err = MPI_Comm_rank(mesh.comm(), &commRank);CHECK_PETSC_ERROR(err);
if (_tstampIndex == istep)
_writeTimeStamp(t, commRank);
-#if 0 // debugging
- field.view("writeVertexField");
- VecView(vector, PETSC_VIEWER_STDOUT_WORLD);
-#endif
-
PetscErrorCode err = 0;
err = PetscViewerHDF5PushGroup(_viewer, "/vertex_fields");CHECK_PETSC_ERROR(err);
err = PetscViewerHDF5SetTimestep(_viewer, istep);CHECK_PETSC_ERROR(err);
@@ -443,30 +429,13 @@
const char* label,
const int labelId)
{ // writeCellField
- typedef typename mesh_type::SieveMesh::numbering_type numbering_type;
-
assert(_viewer);
try {
const char* context = DataWriter<mesh_type, field_type>::_context.c_str();
PetscErrorCode err = 0;
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
- field.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
- int cellDepth = 0;
- err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
- sieveMesh->comm());CHECK_PETSC_ERROR(err);
- const int depth = (0 == label) ? cellDepth : labelId;
- const std::string labelName = (0 == label) ?
- ((sieveMesh->hasLabel("censored depth")) ?
- "censored depth" : "depth") : label;
- assert(!sieveMesh->getFactory().isNull());
- const ALE::Obj<typename mesh_type::SieveMesh::numbering_type>& numbering =
- sieveMesh->getFactory()->getNumbering(sieveMesh, labelName, depth);
- assert(!numbering.isNull());
- field.createScatterWithBC(field.mesh(), numbering, context);
+ field.createScatterWithBC(field.mesh(), PETSC_NULL, context);
field.scatterSectionToVector(context);
PetscVec vector = field.vector(context);
assert(vector);
@@ -477,7 +446,8 @@
_timesteps[field.label()] += 1;
const int istep = _timesteps[field.label()];
// Add time stamp to "/time" if necessary.
- const int commRank = sieveMesh->commRank();
+ PetscMPIInt commRank;
+ err = MPI_Comm_rank(field.mesh().comm(), &commRank);CHECK_PETSC_ERROR(err);
if (_tstampIndex == istep)
_writeTimeStamp(t, commRank);
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -218,6 +218,7 @@
}
err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
err = DMSetCoordinatesLocal(complexMesh, coordVec);CHECK_PETSC_ERROR(err);
+ err = VecDestroy(&coordVec);CHECK_PETSC_ERROR(err);
logger.stagePop(); // Coordinates
sieveMesh->getFactory()->clear();
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -65,7 +65,7 @@
newMesh->coordsys(origMesh.coordsys());
DM newDM;
- PetscErrorCode err = DMComplexDistribute(origMesh.dmMesh(), partitioner, &newDM);CHECK_PETSC_ERROR(err);
+ PetscErrorCode err = DMComplexDistribute(origMesh.dmMesh(), partitioner, 0, &newDM);CHECK_PETSC_ERROR(err);
newMesh->setDMMesh(newDM);
if (0 == strcasecmp(partitioner, "")) {
if (0 == commRank) {
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -57,7 +57,6 @@
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);
@@ -95,7 +94,6 @@
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);
@@ -105,10 +103,28 @@
err = PetscObjectSetName((PetscObject) _localVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
} else {
_dm = PETSC_NULL;
+ _globalVec = PETSC_NULL;
+ _localVec = PETSC_NULL;
}
} // constructor
// ----------------------------------------------------------------------
+// Constructor with mesh, DM, and metadata
+template<typename mesh_type, typename section_type>
+pylith::topology::Field<mesh_type, section_type>::Field(const mesh_type& mesh, DM dm, const Metadata& metadata) :
+ _mesh(mesh),
+ _dm(dm)
+{ // constructor
+ PetscErrorCode err;
+
+ _metadata["default"] = metadata;
+ err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
+ err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) _localVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+} // constructor
+
+// ----------------------------------------------------------------------
// Constructor with field and subfields
template<typename mesh_type, typename section_type>
pylith::topology::Field<mesh_type, section_type>::Field(const Field& src,
@@ -138,7 +154,6 @@
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;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-11-28 00:18:22 UTC (rev 21082)
@@ -77,6 +77,12 @@
*/
Field(const mesh_type& mesh);
+ /** Constructor with mesh, DM, and metadata
+ *
+ * @param mesh Finite-element mesh.
+ */
+ Field(const mesh_type& mesh, DM dm, const Metadata& metadata);
+
/** Constructor with mesh, section, and metadata.
*
* @param mesh Finite-element mesh.
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -201,7 +201,7 @@
const PylithScalar tolerance = 1.0e-06;
for(PetscInt d = 0; d < spaceDim; ++d) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->initialTractions[iVertex * spaceDim + d], initialTractionsArray[d], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(_data->initialTractions[iVertex * spaceDim + d], initialTractionsArray[off+d], tolerance);
} // for
} // for
err = VecRestoreArray(initialTractionsVec, &initialTractionsArray);CHECK_PETSC_ERROR(err);
@@ -622,8 +622,8 @@
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);
+ err = PetscSectionGetDof(dispSection, points[v], &ddof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispSection, points[v], &doff);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(spaceDim, ddof);
const PylithScalar *orientationVertex = &_data->orientation[iVertex*spaceDim*spaceDim];
@@ -809,6 +809,8 @@
err = MatGetSize(jacobianMat, &nrowsM, &ncolsM);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(nrows, nrowsM);
CPPUNIT_ASSERT_EQUAL(ncols, ncolsM);
+ // We ignore the sparsity patterns in our tests
+ err = MatSetOption(jacobianMat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHECK_PETSC_ERROR(err);
int_array rows(nrows);
int_array cols(ncols);
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveImpulses.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -280,6 +280,7 @@
iVertex = 0;
const int fiberDimE = spaceDim;
const PylithScalar tolerance = (sizeof(double) == sizeof(PylithScalar)) ? 1.0e-06 : 1.0e-05;
+ err = VecGetArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
for(PetscInt v = vStart; v < vEnd; ++v, ++iVertex) {
PetscInt dof, off;
@@ -295,6 +296,7 @@
CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
} // Integrate residual with disp increment.
} // testIntegrateResidual
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestTractPerturbation.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -154,25 +154,25 @@
int iPoint = 0;
CPPUNIT_ASSERT(tract._parameters);
- PetscSection valuesSection = tract._parameters->get("values").petscSection();
- Vec valuesVec = tract._parameters->get("values").localVector();
- PetscScalar *valuesArray;
- assert(valuesSection);assert(valuesVec);
+ PetscSection valueSection = tract._parameters->get("value").petscSection();
+ Vec valueVec = tract._parameters->get("value").localVector();
+ PetscScalar *valueArray;
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
- err = VecGetArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(valueVec, &valueArray);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);
+ err = PetscSectionGetDof(valueSection, v, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, v, &voff);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, vdof);
for(PetscInt d = 0; d < spaceDim; ++d) {
const PylithScalar valueE = tractionE[iPoint*spaceDim+d];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, valuesArray[voff+d], tolerance);
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, valueArray[voff+d], tolerance);
} // for
} // for
- err = VecRestoreArray(valuesVec, &valuesArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
} // testCalculate
// ----------------------------------------------------------------------
@@ -180,6 +180,7 @@
void
pylith::faults::TestTractPerturbation::testParameterFields(void)
{ // testParameterFields
+ const int spaceDim = 2;
const int fiberDimE = 7;
const PylithScalar parametersE[2*fiberDimE] = {
0.0, 0.0, 2.0, -1.0, -1.0, 0.5, 1.5,
@@ -192,27 +193,78 @@
_initialize(&mesh, &faultMesh, &tract);
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());
- 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);
+
+ CPPUNIT_ASSERT(tract._parameters);
+ PetscSection valueSection = tract._parameters->get("value").petscSection();
+ Vec valueVec = tract._parameters->get("value").localVector();
+ PetscScalar *valueArray;
+ CPPUNIT_ASSERT(valueSection);CPPUNIT_ASSERT(valueVec);
+ PetscSection initialSection = tract._parameters->get("initial").petscSection();
+ Vec initialVec = tract._parameters->get("initial").localVector();
+ PetscScalar *initialArray;
+ CPPUNIT_ASSERT(initialSection);CPPUNIT_ASSERT(initialVec);
+ PetscSection changeSection = tract._parameters->get("change").petscSection();
+ Vec changeVec = tract._parameters->get("change").localVector();
+ PetscScalar *changeArray;
+ CPPUNIT_ASSERT(changeSection);CPPUNIT_ASSERT(changeVec);
+ PetscSection changeTimeSection = tract._parameters->get("change time").petscSection();
+ Vec changeTimeVec = tract._parameters->get("change time").localVector();
+ PetscScalar *changeTimeArray;
+ CPPUNIT_ASSERT(changeTimeSection);CPPUNIT_ASSERT(changeTimeVec);
+
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);
+ err = VecGetArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+ PetscInt vdof, voff, e = 0;
- for(PetscInt d = 0; d < fiberDimE; ++d) {
- const PylithScalar valueE = parametersE[iPoint*fiberDimE+d];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, vals[d], tolerance);
+ err = PetscSectionGetDof(valueSection, v, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(valueSection, v, &voff);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, vdof);
+ PetscInt idof, ioff;
+
+ err = PetscSectionGetDof(initialSection, v, &idof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(initialSection, v, &ioff);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, idof);
+ PetscInt cdof, coff;
+
+ err = PetscSectionGetDof(changeSection, v, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(changeSection, v, &coff);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim, cdof);
+ PetscInt ctdof, ctoff;
+
+ err = PetscSectionGetDof(changeTimeSection, v, &ctdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(changeTimeSection, v, &ctoff);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(1, ctdof);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, vdof+idof+cdof+ctdof);
+
+ for(PetscInt d = 0; d < vdof; ++d, ++e) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(parametersE[iPoint*fiberDimE+e], valueArray[voff+d], tolerance);
} // for
+ for(PetscInt d = 0; d < idof; ++d, ++e) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(parametersE[iPoint*fiberDimE+e], initialArray[ioff+d], tolerance);
+ } // for
+ for(PetscInt d = 0; d < cdof; ++d, ++e) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(parametersE[iPoint*fiberDimE+e], changeArray[coff+d], tolerance);
+ } // for
+ for(PetscInt d = 0; d < ctdof; ++d, ++e) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(parametersE[iPoint*fiberDimE+e], changeTimeArray[ctoff+d], tolerance);
+ } // for
} // for
+ err = VecRestoreArray(valueVec, &valueArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(initialVec, &initialArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeVec, &changeArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(changeTimeVec, &changeTimeArray);CHECK_PETSC_ERROR(err);
} // testParameterFields
// ----------------------------------------------------------------------
@@ -233,27 +285,33 @@
_initialize(&mesh, &faultMesh, &tract);
const topology::Field<topology::SubMesh>& field = tract.vertexField(label);
- const ALE::Obj<RealSection>& section = field.section();
- CPPUNIT_ASSERT(!section.isNull());
+ PetscSection section = field.petscSection();
+ Vec vec = field.localVector();
+ PetscScalar *array;
+ CPPUNIT_ASSERT(section);CPPUNIT_ASSERT(vec);
- 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);
+
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 = section->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
- const PylithScalar* vals = section->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(vals);
+ err = VecGetArray(vec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v, ++iPoint) {
+ PetscInt dof, off;
- for (int iDim=0; iDim < fiberDim; ++iDim) {
- const PylithScalar valueE = fieldE[iPoint*fiberDim+iDim];
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, vals[iDim], tolerance);
+ err = PetscSectionGetDof(section, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(section, v, &off);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, dof);
+
+ for(PetscInt d = 0; d < dof; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(fieldE[iPoint*fiberDimE+d], array[off+d], tolerance);
} // for
} // for
+ err = VecRestoreArray(vec, &array);CHECK_PETSC_ERROR(err);
} // testVertexField
// ----------------------------------------------------------------------
@@ -289,13 +347,17 @@
mesh->coordsys(&cs);
const int spaceDim = cs.spaceDim();
+ PetscInt labelSize;
+ PetscErrorCode err;
+ err = DMComplexGetStratumSize(mesh->dmMesh(), faultLabel, 1, &labelSize);CHECK_PETSC_ERROR(err);
+
// Create fault mesh
- int firstFaultVertex = 0;
- int firstLagrangeVertex = mesh->sieveMesh()->getIntSection(faultLabel)->size();
- int firstFaultCell = mesh->sieveMesh()->getIntSection(faultLabel)->size();
+ PetscInt firstFaultVertex = 0;
+ PetscInt firstLagrangeVertex = labelSize;
+ PetscInt firstFaultCell = labelSize;
const bool useLagrangeConstraints = true;
if (useLagrangeConstraints) {
- firstFaultCell += mesh->sieveMesh()->getIntSection(faultLabel)->size();
+ firstFaultCell += labelSize;
} // if
ALE::Obj<SieveFlexMesh> faultBoundary = 0;
const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
@@ -311,9 +373,46 @@
// using create() instead of createParallel().
const ALE::Obj<SieveSubMesh>& faultSieveMesh = faultMesh->sieveMesh();
CPPUNIT_ASSERT(!faultSieveMesh.isNull());
- faultSieveMesh->setRealSection("coordinates",
- sieveMesh->getRealSection("coordinates"));
+ const ALE::Obj<RealSection>& oldCoordSection = sieveMesh->getRealSection("coordinates");
+ faultSieveMesh->setRealSection("coordinates", oldCoordSection);
+ DM faultDMMesh = faultMesh->dmMesh();
+ IS subpointMap;
+ const PetscInt *points;
+ PetscSection coordSection;
+ PetscInt vStart, vEnd;
+ CPPUNIT_ASSERT(faultDMMesh);
+ err = DMComplexGetDepthStratum(faultDMMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetSubpointMap(faultDMMesh, &subpointMap);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(coordSection, vStart, vEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ err = PetscSectionSetDof(coordSection, v, spaceDim);CHECK_PETSC_ERROR(err);
+ }
+ err = PetscSectionSetUp(coordSection);CHECK_PETSC_ERROR(err);
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt coordSize;
+
+ err = PetscSectionGetStorageSize(coordSection, &coordSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(mesh->comm(), &coordVec);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(coordVec, coordSize, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(coordVec);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt off;
+
+ err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
+ const PetscScalar *oldCoords = oldCoordSection->restrictPoint(points[v]);
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ coords[off+d] = oldCoords[d];
+ }
+ }
+ err = ISRestoreIndices(subpointMap, &points);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ err = DMSetCoordinatesLocal(faultDMMesh, coordVec);CHECK_PETSC_ERROR(err);
+
// Setup databases
spatialdata::spatialdb::SimpleDB dbInitial("initial traction");
spatialdata::spatialdb::SimpleIOAscii ioInitial;
@@ -326,17 +425,25 @@
dbChange.ioHandler(&ioChange);
// Setup fault orientation
- const ALE::Obj<SieveMesh::label_sequence>& vertices = faultSieveMesh->depthStratum(0);
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
topology::Field<topology::SubMesh> faultOrientation(*faultMesh);
- faultOrientation.newSection(vertices, spaceDim*spaceDim);
+ faultOrientation.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim*spaceDim);
faultOrientation.allocate();
- const ALE::Obj<RealSection>& orientationSection = faultOrientation.section();
- CPPUNIT_ASSERT(!orientationSection.isNull());
- for (SieveMesh::label_sequence::iterator v_iter=vertices->begin(); v_iter != verticesEnd; ++v_iter) {
- CPPUNIT_ASSERT_EQUAL(spaceDim*spaceDim, orientationSection->getFiberDimension(*v_iter));
- orientationSection->updatePoint(*v_iter, orientationVertex);
+ PetscSection orientationSection = faultOrientation.petscSection();
+ Vec orientationVec = faultOrientation.localVector();
+ PetscScalar *orientationArray;
+ CPPUNIT_ASSERT(orientationSection);CPPUNIT_ASSERT(orientationVec);
+ err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt odof, ooff;
+
+ err = PetscSectionGetDof(orientationSection, v, &odof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, v, &ooff);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(spaceDim*spaceDim, odof);
+ for(PetscInt d = 0; d < odof; ++d) {
+ orientationArray[ooff+d] = orientationVertex[d];
+ }
} // for
+ err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
spatialdata::units::Nondimensional normalizer;
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataHex8.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -27,10 +27,10 @@
*
* After adding cohesive elements
*
- * Cells are 0-1,16 and vertices are 4-15.
+ * Cells are 0-1,2 and vertices are 3-18,19-22.
*
- * 2,3,4,5 -------- 6,7,8,9 -- 14,15,16,17 -------- 10,11,12,13
- * 18,19,20,21
+ * 3,4,5,6 -------- 7,8,9,10 -- 15,16,17,18 -------- 11,12,13,14
+ * 19,20,21,22
* ^^^^^^^^^^^^^^^^^^^^^^ Cohesive element
*
*/
@@ -1350,7 +1350,7 @@
const int pylith::faults::CohesiveDynDataHex8::_numConstraintVert = 4;
const int pylith::faults::CohesiveDynDataHex8::_constraintVertices[] = {
- 18, 19, 20, 21
+ 19, 20, 21, 22
};
// ----------------------------------------------------------------------
// Stick case
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataQuad4.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -34,9 +34,9 @@
*
* After adding cohesive elements
*
- * Cells are 0-1,10 vertices are 2-9.
+ * Cells are 0-1,2 vertices are 3-10,11-12.
*
- * 3 -------- 5 -11-- 9 -------- 7
+ * 4 -------- 6 -12--10 -------- 8
* | | | |
* | | | |
* | | | |
@@ -45,7 +45,7 @@
* | | | |
* | | | |
* | | | |
- * 2 -------- 4 -10-- 8 -------- 6
+ * 3 -------- 5 -11-- 9 -------- 7
*/
#include "CohesiveDynDataQuad4.hh"
@@ -289,7 +289,7 @@
const int pylith::faults::CohesiveDynDataQuad4::_numConstraintVert = 2;
const int pylith::faults::CohesiveDynDataQuad4::_constraintVertices[] = {
- 10, 11
+ 11, 12
};
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTet4.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -27,10 +27,10 @@
*
* After adding cohesive elements
*
- * Cells are 0-1,10, vertices are 2-9.
+ * Cells are 0-1,2, vertices are 3-10,11-13.
*
- * 2 3,4,5 7,8,9 6
- * 10,11,12
+ * 3 4,5,6 8,9,10 7
+ * 11,12,13
* ^^^^^^^^^^^^ Cohesive element in x-y plane.
*/
@@ -494,7 +494,7 @@
const int pylith::faults::CohesiveDynDataTet4::_numConstraintVert = 3;
const int pylith::faults::CohesiveDynDataTet4::_constraintVertices[] = {
- 10, 11, 12
+ 11, 12, 13
};
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveDynDataTri3d.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -41,27 +41,27 @@
*
* After adding cohesive elements
*
- * Cells are 0-3, 13-14, vertices are 4-12.
+ * Cells are 0-3, 4-5, vertices are 6-14,15-17.
*
- * 9
+ * 11
* / \
* / \
* / \
* / \
- * 8--------- 5
- * 15 | 13/|
- * 12--------10 |
+ * 10--------- 7
+ * 17 | 15/|
+ * 14--------12 |
* \ /| |\
* \ / | | \
* \ / | | \
* \ / | | \
- * 4 | | 7
+ * 6 | | 9
* \ | | /
* \ | | /
* \ | | /
* \| |/
- * 11-6
- * 14
+ * 13-8
+ * 16
*/
@@ -433,7 +433,7 @@
const int pylith::faults::CohesiveDynDataTri3d::_numConstraintVert = 3;
const int pylith::faults::CohesiveDynDataTri3d::_constraintVertices[] = {
- 13, 14, 15
+ 15, 16, 17
};
const PylithScalar pylith::faults::CohesiveDynDataTri3d::_initialTractions[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataHex8.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataHex8.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -26,10 +26,10 @@
*
* After adding cohesive elements
*
- * Cells are 0-1,16 and vertices are 4-15.
+ * Cells are 0-1,2 and vertices are 3-18,19-22.
*
- * 2,3,4,5 -------- 6,7,8,9 -- 14,15,16,17 -------- 10,11,12,13
- * 18,19,20,21
+ * 3,4,5,6 -------- 7,8,9,10 -- 15,16,17,18 -------- 11,12,13,14
+ * 19,20,21,22
* ^^^^^^^^^^^^^^^^^^^^^^ Cohesive element
*
*/
@@ -178,7 +178,7 @@
const int pylith::faults::CohesiveImpulsesDataHex8::_numConstraintVert = 4;
const int pylith::faults::CohesiveImpulsesDataHex8::_constraintVertices[4] = {
- 18, 19, 20, 21
+ 19, 20, 21, 22
};
const PylithScalar pylith::faults::CohesiveImpulsesDataHex8::_residualIncr[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataQuad4.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataQuad4.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -33,9 +33,9 @@
*
* After adding cohesive elements
*
- * Cells are 0-1,10 vertices are 2-9.
+ * Cells are 0-1,2 vertices are 3-10,11-12.
*
- * 3 -------- 5 -11-- 9 -------- 7
+ * 4 -------- 6 -12--10 -------- 8
* | | | |
* | | | |
* | | | |
@@ -44,7 +44,7 @@
* | | | |
* | | | |
* | | | |
- * 2 -------- 4 -10-- 8 -------- 6
+ * 3 -------- 5 -11-- 9 -------- 7
*/
#include "CohesiveImpulsesDataQuad4.hh"
@@ -143,7 +143,7 @@
const int pylith::faults::CohesiveImpulsesDataQuad4::_numConstraintVert = 2;
const int pylith::faults::CohesiveImpulsesDataQuad4::_constraintVertices[] = {
- 10, 11
+ 11, 12
};
const PylithScalar pylith::faults::CohesiveImpulsesDataQuad4::_residualIncr[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataTet4.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataTet4.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -26,10 +26,10 @@
*
* After adding cohesive elements
*
- * Cells are 0-1,13, vertices are 2-9.
+ * Cells are 0-1,2, vertices are 3-10,11-13.
*
- * 2 3,4,5 7,8,9 6
- * 10,11,12
+ * 3 4,5,6 8,9,10 7
+ * 11,12,13
* ^^^^^^^^^^^^ Cohesive element in x-y plane.
*/
@@ -146,7 +146,7 @@
const int pylith::faults::CohesiveImpulsesDataTet4::_numConstraintVert = 3;
const int pylith::faults::CohesiveImpulsesDataTet4::_constraintVertices[] = {
- 10, 11, 12,
+ 11, 12, 13,
};
const PylithScalar pylith::faults::CohesiveImpulsesDataTet4::_residualIncr[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataTri3.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveImpulsesDataTri3.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -35,19 +35,19 @@
*
* After adding cohesive elements
*
- * Cells are 0-1, 10, vertices are 2-7.
+ * Cells are 0-1, 2, vertices are 3-8,9-10.
*
- * 6 -8- 3
+ * 7 -9- 4
* /| |\
* / | | \
* / | | \
* / | | \
- * 2 | | 5
+ * 3 | | 6
* \ | | /
* \ | | /
* \ | | /
* \| |/
- * 7 -9- 4
+ * 8-10- 5
*/
#include "CohesiveImpulsesDataTri3.hh"
@@ -142,7 +142,7 @@
const int pylith::faults::CohesiveImpulsesDataTri3::_numConstraintVert = 2;
const int pylith::faults::CohesiveImpulsesDataTri3::_constraintVertices[] = {
- 8, 9
+ 9, 10
};
const PylithScalar pylith::faults::CohesiveImpulsesDataTri3::_residualIncr[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterVTKFaultMeshCases.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterVTKFaultMeshCases.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestDataWriterVTKFaultMeshCases.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -40,7 +40,7 @@
{ // setUp
TestDataWriterVTKFaultMesh::setUp();
_data = new DataWriterVTKDataFaultMeshTri3;
- _flipFault = false;
+ _flipFault = true;
_initialize();
} // setUp
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/TestMeshIO.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -69,47 +69,61 @@
new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
ALE::SieveBuilder<SieveFlexMesh>::buildTopology(s, data.cellDim, data.numCells,
- const_cast<int*>(data.cells),
- data.numVertices,
- interpolate, data.numCorners);
+ const_cast<int*>(data.cells),
+ data.numVertices,
+ interpolate, data.numCorners);
std::map<SieveFlexMesh::point_type,SieveFlexMesh::point_type> renumbering;
ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
sieveMesh->setSieve(sieve);
sieveMesh->stratify();
- ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, data.spaceDim,
- data.vertices);
+ ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, data.spaceDim, data.vertices);
+ DM dmMesh;
+ PetscBool interpolateMesh = interpolate ? PETSC_TRUE : PETSC_FALSE;
+ PetscErrorCode err;
+
+ err = DMComplexCreateFromCellList(mesh->comm(), data.cellDim, data.numCells, data.numVertices, data.numCorners, interpolateMesh, data.cells, data.vertices, &dmMesh);CHECK_PETSC_ERROR(err);
+ mesh->setDMMesh(dmMesh);
+
// Material ids
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->heightStratum(0);
+ const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
CPPUNIT_ASSERT(!cells.isNull());
- const ALE::Obj<SieveMesh::label_type>& labelMaterials =
- sieveMesh->createLabel("material-id");
+ const ALE::Obj<SieveMesh::label_type>& labelMaterials = sieveMesh->createLabel("material-id");
CPPUNIT_ASSERT(!labelMaterials.isNull());
int i = 0;
- for(SieveMesh::label_sequence::iterator e_iter=cells->begin();
- e_iter != cells->end();
- ++e_iter)
+ for(SieveMesh::label_sequence::iterator e_iter=cells->begin(); e_iter != cells->end(); ++e_iter)
sieveMesh->setValue(labelMaterials, *e_iter, data.materialIds[i++]);
+ PetscInt cStart, cEnd;
+
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ err = DMComplexSetLabelValue(dmMesh, "material-id", c, data.materialIds[c-cStart]);CHECK_PETSC_ERROR(err);
+ }
+
// Groups
for (int iGroup=0, index=0; iGroup < data.numGroups; ++iGroup) {
- const ALE::Obj<SieveMesh::int_section_type>& groupField =
- sieveMesh->getIntSection(data.groupNames[iGroup]);
+ const ALE::Obj<SieveMesh::int_section_type>& groupField = sieveMesh->getIntSection(data.groupNames[iGroup]);
CPPUNIT_ASSERT(!groupField.isNull());
groupField->setChart(sieveMesh->getSieve()->getChart());
+ err = DMComplexCreateLabel(dmMesh, data.groupNames[iGroup]);CHECK_PETSC_ERROR(err);
+
MeshIO::GroupPtType type;
const int numPoints = data.groupSizes[iGroup];
if (0 == strcasecmp("cell", data.groupTypes[iGroup])) {
type = MeshIO::CELL;
- for(int i=0; i < numPoints; ++i)
- groupField->setFiberDimension(data.groups[index++], 1);
+ for(int i=0; i < numPoints; ++i, ++index) {
+ groupField->setFiberDimension(data.groups[index], 1);
+ err = DMComplexSetLabelValue(dmMesh, data.groupNames[iGroup], data.groups[index], 1);CHECK_PETSC_ERROR(err);
+ }
} else if (0 == strcasecmp("vertex", data.groupTypes[iGroup])) {
type = MeshIO::VERTEX;
const int numCells = sieveMesh->heightStratum(0)->size();
- for(int i=0; i < numPoints; ++i)
- groupField->setFiberDimension(data.groups[index++]+numCells, 1);
+ for(int i=0; i < numPoints; ++i, ++index) {
+ groupField->setFiberDimension(data.groups[index]+numCells, 1);
+ err = DMComplexSetLabelValue(dmMesh, data.groupNames[iGroup], data.groups[index]+numCells, 1);CHECK_PETSC_ERROR(err);
+ }
} else
throw std::logic_error("Could not parse group type.");
sieveMesh->allocate(groupField);
@@ -140,6 +154,9 @@
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
CPPUNIT_ASSERT(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ CPPUNIT_ASSERT(dmMesh);
+
// Check vertices
const ALE::Obj<SieveMesh::label_sequence>& vertices =
sieveMesh->depthStratum(0);
@@ -169,7 +186,35 @@
tolerance);
}
} // for
+ PetscSection coordSection;
+ Vec coordVec;
+ PetscScalar *coords;
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT(coordSection);CPPUNIT_ASSERT(coordVec);
+ err = VecGetArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(coordSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(coordSection, v, &off);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dof);
+ const PylithScalar tolerance = 1.0e-06;
+ for(PetscInt d = 0; d < spaceDim; ++d) {
+ const PetscInt eoff = (v-vStart)*spaceDim;
+ if (data.vertices[eoff+d] < 1.0) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(data.vertices[eoff+d], coords[off+d], tolerance);
+ } else {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, coords[off+d]/data.vertices[eoff+d], tolerance);
+ }
+ }
+ }
+ err = VecRestoreArray(coordVec, &coords);CHECK_PETSC_ERROR(err);
+
// check cells
const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
@@ -192,10 +237,24 @@
}
pV.clear();
} // for
+ PetscInt cStart, cEnd;
+ err = DMComplexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(data.numCells, cEnd-cStart);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ const PetscInt *cone;
+ PetscInt coneSize;
+
+ err = DMComplexGetConeSize(dmMesh, c, &coneSize);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetCone(dmMesh, c, &cone);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(data.numCorners, coneSize);
+ for(PetscInt p = 0; p < coneSize; ++p) {
+ CPPUNIT_ASSERT_EQUAL(data.cells[(c-cStart)*coneSize+p], cone[p]-offset);
+ }
+ }
+
// check materials
- const ALE::Obj<SieveMesh::label_type>& labelMaterials =
- sieveMesh->getLabel("material-id");
+ const ALE::Obj<SieveMesh::label_type>& labelMaterials = sieveMesh->getLabel("material-id");
const int idDefault = -999;
const int size = numCells;
int_array materialIds(size);
@@ -208,7 +267,15 @@
for (int iCell=0; iCell < numCells; ++iCell)
CPPUNIT_ASSERT_EQUAL(data.materialIds[iCell], materialIds[iCell]);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ PetscInt matId;
+
+ err = DMComplexGetLabelValue(dmMesh, "material-id", c, &matId);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(data.materialIds[c-cStart], matId);
+ }
+
// Check groups
+#if 0
const ALE::Obj<std::set<std::string> >& groupNames =
sieveMesh->getIntSections();
if (data.numGroups > 0) {
@@ -251,6 +318,46 @@
for (int i=0; i < numPoints; ++i)
CPPUNIT_ASSERT_EQUAL(data.groups[index++], points[i]);
} // for
+#endif
+
+ PetscInt numLabels, pStart, pEnd;
+
+ err = DMComplexGetChart(dmMesh, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ err = DMComplexGetNumLabels(dmMesh, &numLabels);CHECK_PETSC_ERROR(err);
+ numLabels -= 2; /* Remove depth and material labels */
+ CPPUNIT_ASSERT_EQUAL(data.numGroups, numLabels);
+ PetscInt iGroup = 0;
+ PetscInt index = 0;
+ for(PetscInt l = numLabels-1; l >= 0; --l, ++iGroup) {
+ const char *name;
+ PetscInt firstPoint;
+
+ err = DMComplexGetLabelName(dmMesh, l, &name);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(std::string(data.groupNames[iGroup]), std::string(name));
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ PetscInt val;
+
+ err = DMComplexGetLabelValue(dmMesh, name, p, &val);CHECK_PETSC_ERROR(err);
+ if (val >= 0) {
+ firstPoint = p;
+ break;
+ } // if
+ } // for
+ std::string groupType = (firstPoint >= cStart && firstPoint < cEnd) ? "cell" : "vertex";
+ CPPUNIT_ASSERT_EQUAL(std::string(data.groupTypes[iGroup]), groupType);
+ PetscInt numPoints;
+ err = DMComplexGetStratumSize(dmMesh, name, 1, &numPoints);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT_EQUAL(data.groupSizes[iGroup], numPoints);
+ IS pointIS;
+ const PetscInt *points;
+ const PetscInt offset = ("vertex" == groupType) ? numCells : 0;
+ err = DMComplexGetStratumIS(dmMesh, name, 1, &pointIS);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(pointIS, &points);CHECK_PETSC_ERROR(err);
+ for(PetscInt p = 0; p < numPoints; ++p) {
+ CPPUNIT_ASSERT_EQUAL(data.groups[index++], points[p]-offset);
+ }
+ err = ISRestoreIndices(pointIS, &points);CHECK_PETSC_ERROR(err);
+ }
} // _checkVals
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshHex8.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshHex8.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -53,7 +53,11 @@
{ "other", topology::FieldBase::OTHER, 2 },
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshHex8::_vertexFieldScalar[12*1] = {
+#if 0
2.1, 3.2, 4.3, 5.4, 6.5, 7.6, 8.7, 9.8, 10.9, 11.8, 12.7, 13.6
+#else
+ 8.7, 9.8, 10.9, 11.8, 12.7, 13.6
+#endif
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshHex8::_vertexFieldVector[12*3] = {
7.8, 8.9, 9.0,
@@ -64,6 +68,7 @@
10.3, 11.4, 12.5,
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshHex8::_vertexFieldTensor[12*6] = {
+#if 0
1.1, 1.2, 1.3, 1.4, 1.5, 1.6,
2.1, 2.2, 2.3, 2.4, 2.5, 2.6,
3.1, 3.2, 3.3, 3.4, 3.5, 3.6,
@@ -76,6 +81,14 @@
10.1, 10.2, 10.3, 10.4, 10.5, 10.6,
11.1, 11.2, 11.3, 11.4, 11.5, 11.6,
12.1, 12.2, 12.3, 12.4, 12.5, 12.6,
+#else
+ 7.1, 7.2, 7.3, 7.4, 7.5, 7.6,
+ 8.1, 8.2, 8.3, 8.4, 8.5, 8.6,
+ 9.1, 9.2, 9.3, 9.4, 9.5, 9.6,
+ 10.1, 10.2, 10.3, 10.4, 10.5, 10.6,
+ 11.1, 11.2, 11.3, 11.4, 11.5, 11.6,
+ 12.1, 12.2, 12.3, 12.4, 12.5, 12.6,
+#endif
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshHex8::_vertexFieldOther[12*2] = {
5.7, 6.8,
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshQuad4.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshQuad4.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -50,7 +50,11 @@
{ "other", topology::FieldBase::OTHER, 2 },
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshQuad4::_vertexFieldScalar[6*1] = {
+#if 0
2.1, 3.2, 4.3, 5.4, 6.5, 7.6,
+#else
+ 2.1, 4.3, 6.5,
+#endif
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshQuad4::_vertexFieldVector[6*2] = {
1.1, 2.2,
@@ -58,12 +62,18 @@
9.9, 10.0,
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshQuad4::_vertexFieldTensor[6*3] = {
+#if 0
1.1, 1.2, 1.3,
2.1, 2.2, 2.3,
3.1, 3.2, 3.3,
4.1, 4.2, 4.3,
5.1, 5.2, 5.3,
6.1, 6.2, 6.3,
+#else
+ 1.1, 1.2, 1.3,
+ 3.1, 3.2, 3.3,
+ 5.1, 5.2, 5.3,
+#endif
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshQuad4::_vertexFieldOther[6*3] = {
1.2, 2.3,
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshTet4.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshTet4.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -53,7 +53,11 @@
{ "other", topology::FieldBase::OTHER, 2 },
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshTet4::_vertexFieldScalar[5*1] = {
+#if 0
2.1, 3.2, 4.3, 5.4, 6.5,
+#else
+ 2.1, 3.2, 5.4, 6.5,
+#endif
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshTet4::_vertexFieldVector[5*3] = {
1.1, 2.2, 3.3,
@@ -62,11 +66,18 @@
13.3, 14.4, 15.5,
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshTet4::_vertexFieldTensor[5*6] = {
+#if 0
1.1, 1.2, 1.3, 1.4, 1.5, 1.6,
2.1, 2.2, 2.3, 2.4, 2.5, 2.6,
3.1, 3.2, 3.3, 3.4, 3.5, 3.6,
4.1, 4.2, 4.3, 4.4, 4.5, 4.6,
5.1, 5.2, 5.3, 5.4, 5.5, 5.6,
+#else
+ 1.1, 1.2, 1.3, 1.4, 1.5, 1.6,
+ 2.1, 2.2, 2.3, 2.4, 2.5, 2.6,
+ 4.1, 4.2, 4.3, 4.4, 4.5, 4.6,
+ 5.1, 5.2, 5.3, 5.4, 5.5, 5.6,
+#endif
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshTet4::_vertexFieldOther[5*2] = {
1.2, 2.3,
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshTri3.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/DataWriterVTKDataSubMeshTri3.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -54,13 +54,18 @@
{ "other", topology::FieldBase::OTHER, 2 },
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshTri3::_vertexFieldScalar[8*1] = {
+#if 0
2.1, 3.2, 4.3, 5.4, 6.5, 7.6, 8.7, 9.8,
+#else
+ 3.2, 5.4,
+#endif
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshTri3::_vertexFieldVector[8*2] = {
3.3, 4.4,
7.7, 8.8,
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshTri3::_vertexFieldTensor[8*3] = {
+#if 0
1.1, 1.2, 1.3,
2.1, 2.2, 2.3,
3.1, 3.2, 3.3,
@@ -69,6 +74,10 @@
6.1, 6.2, 6.3,
7.1, 7.2, 7.3,
8.1, 8.2, 8.3,
+#else
+ 2.1, 2.2, 2.3,
+ 4.1, 4.2, 4.3,
+#endif
};
const PylithScalar pylith::meshio::DataWriterVTKDataSubMeshTri3::_vertexFieldOther[8*2] = {
3.4, 4.5,
Modified: short/3D/PyLith/trunk/unittests/libtests/meshio/data/MeshDataCubitTet.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/meshio/data/MeshDataCubitTet.cc 2012-11-27 23:55:19 UTC (rev 21081)
+++ short/3D/PyLith/trunk/unittests/libtests/meshio/data/MeshDataCubitTet.cc 2012-11-28 00:18:22 UTC (rev 21082)
@@ -47,16 +47,16 @@
const int pylith::meshio::MeshDataCubitTet::_numGroups = 2;
const int pylith::meshio::MeshDataCubitTet::_groupSizes[] =
- { 4, 3 };
+ { 3, 4 };
const int pylith::meshio::MeshDataCubitTet::_groups[] = {
+ 1, 2, 3,
0, 1, 2, 3,
- 1, 2, 3,
};
const char* pylith::meshio::MeshDataCubitTet::_groupNames[] = {
+ "mid_face",
"bottom_face",
- "mid_face",
};
const char* pylith::meshio::MeshDataCubitTet::_groupTypes[] = {
More information about the CIG-COMMITS
mailing list