[cig-commits] r20734 - in short/3D/PyLith/trunk: libsrc/pylith/bc libsrc/pylith/meshio libsrc/pylith/topology unittests/libtests/bc
knepley at geodynamics.org
knepley at geodynamics.org
Wed Sep 19 09:45:05 PDT 2012
Author: knepley
Date: 2012-09-19 09:45:05 -0700 (Wed, 19 Sep 2012)
New Revision: 20734
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
Log:
Fix for CellFilterAvg with a label, some work on SubMesh
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2012-09-19 16:14:39 UTC (rev 20733)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2012-09-19 16:45:05 UTC (rev 20734)
@@ -296,14 +296,15 @@
assert(!parametersSection.isNull());
// Use _cellVector for cell residual.
- const ALE::Obj<RealSection>& residualSection = residual.section();
- assert(!residualSection.isNull());
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ assert(residualSection);assert(residualVec);
UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
scalar_array velCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& velSection =
- fields->get("velocity(t)").section();
- assert(!velSection.isNull());
+ PetscSection velSection = fields->get("velocity(t)").petscSection();
+ PetscSection velVec = fields->get("velocity(t)").localVector();
+ assert(velSection);assert(velVec);
RestrictVisitor velVisitor(*velSection, velCell.size(), &velCell[0]);
#if !defined(PRECOMPUTE_GEOMETRY)
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc 2012-09-19 16:14:39 UTC (rev 20733)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc 2012-09-19 16:45:05 UTC (rev 20734)
@@ -113,6 +113,9 @@
} else {
err = DMComplexGetStratumIS(dmMesh, label, 1, &cellIS);CHECK_PETSC_ERROR(err);
err = ISGetSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ cStart = cells[0];
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
}
// Only processors with cells for output get the correct fiber dimension.
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc 2012-09-19 16:14:39 UTC (rev 20733)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.cc 2012-09-19 16:45:05 UTC (rev 20734)
@@ -76,12 +76,22 @@
const ALE::Obj<DomainSieveMesh>& meshSieveMesh = mesh.sieveMesh();
assert(!meshSieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ assert(dmMesh);
if (!meshSieveMesh->hasIntSection(label)) {
std::ostringstream msg;
msg << "Could not find group of points '" << label << "' in mesh.";
throw std::runtime_error(msg.str());
} // if
+ PetscBool hasLabel;
+ PetscErrorCode err;
+ err = DMComplexHasLabel(dmMesh, label, &hasLabel);CHECK_PETSC_ERROR(err);
+ if (!hasLabel) {
+ std::ostringstream msg;
+ msg << "Could not find group of points '" << label << "' in DM mesh.";
+ throw std::runtime_error(msg.str());
+ } // if
const ALE::Obj<IntSection>& groupField = meshSieveMesh->getIntSection(label);
assert(!groupField.isNull());
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2012-09-19 16:14:39 UTC (rev 20733)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2012-09-19 16:45:05 UTC (rev 20734)
@@ -180,29 +180,35 @@
const PylithScalar t = 0.0;
bc.integrateResidual(residual, t, &fields);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- CPPUNIT_ASSERT(!sieveMesh->depthStratum(0).isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscInt vStart, vEnd;
+ const PylithScalar* valsE = _data->valsResidual;
+ PetscErrorCode err;
- const PylithScalar* valsE = _data->valsResidual;
- const int totalNumVertices = sieveMesh->depthStratum(0)->size();
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ const int totalNumVertices = vEnd - vStart;
const int sizeE = _data->spaceDim * totalNumVertices;
- const ALE::Obj<RealSection>& residualSection = residual.section();
- CPPUNIT_ASSERT(!residualSection.isNull());
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *vals;
+ PetscInt size;
- const PylithScalar* vals = residualSection->restrictSpace();
- const int size = residualSection->sizeWithBC();
+ CPPUNIT_ASSERT(residualSection);CPPUNIT_ASSERT(residualVec);
+ err = PetscSectionGetStorageSize(residualSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
//residual->view("RESIDUAL");
const PylithScalar tolerance = 1.0e-06;
- for (int i=0; i < size; ++i)
+ err = VecGetArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
+ for(int i = 0; i < size; ++i)
if (fabs(valsE[i]) > 1.0)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
+ err = VecRestoreArray(residualVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateResidual
// ----------------------------------------------------------------------
@@ -217,16 +223,10 @@
topology::SolutionFields fields(mesh);
_initialize(&mesh, &bc, &fields);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
-
const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
topology::Field<topology::Mesh>& solution = fields.solution();
- const ALE::Obj<RealSection>& solutionSection = solution.section();
- CPPUNIT_ASSERT(!solutionSection.isNull());
-
topology::Jacobian jacobian(solution);
const PylithScalar t = 1.0;
@@ -234,22 +234,26 @@
CPPUNIT_ASSERT_EQUAL(false, bc.needNewJacobian());
jacobian.assemble("final_assembly");
- CPPUNIT_ASSERT(!sieveMesh->depthStratum(0).isNull());
+ DM dmMesh = mesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const PylithScalar* valsE = _data->valsJacobian;
- const int totalNumVertices = sieveMesh->depthStratum(0)->size();
+ const int totalNumVertices = vEnd - vStart;
const int nrowsE = totalNumVertices * _data->spaceDim;
const int ncolsE = totalNumVertices * _data->spaceDim;
const PetscMat jacobianMat = jacobian.matrix();
int nrows = 0;
int ncols = 0;
- MatGetSize(jacobianMat, &nrows, &ncols);
+ err = MatGetSize(jacobianMat, &nrows, &ncols);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(nrowsE, nrows);
CPPUNIT_ASSERT_EQUAL(ncolsE, ncols);
PetscMat jDense;
- MatConvert(jacobianMat, MATSEQDENSE, MAT_INITIAL_MATRIX, &jDense);
+ err = MatConvert(jacobianMat, MATSEQDENSE, MAT_INITIAL_MATRIX, &jDense);CHECK_PETSC_ERROR(err);
scalar_array vals(nrows*ncols);
int_array rows(nrows);
@@ -258,7 +262,7 @@
rows[iRow] = iRow;
for (int iCol=0; iCol < ncols; ++iCol)
cols[iCol] = iCol;
- MatGetValues(jDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);
+ err = MatGetValues(jDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);CHECK_PETSC_ERROR(err);
#if 0
std::cout << "JACOBIAN\n";
@@ -276,7 +280,7 @@
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[index], vals[index], tolerance);
} // for
- MatDestroy(&jDense);
+ err = MatDestroy(&jDense);CHECK_PETSC_ERROR(err);
} // testIntegrateJacobian
// ----------------------------------------------------------------------
@@ -297,24 +301,23 @@
jacobian.newSection(topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
jacobian.allocate();
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
-
const topology::SubMesh& boundaryMesh = *bc._boundaryMesh;
const ALE::Obj<SieveSubMesh>& submesh = boundaryMesh.sieveMesh();
CPPUNIT_ASSERT(!submesh.isNull());
- topology::Field<topology::Mesh>& solution = fields.solution();
- const ALE::Obj<RealSection>& solutionSection = solution.section();
- CPPUNIT_ASSERT(!solutionSection.isNull());
-
const PylithScalar t = 1.0;
bc.integrateJacobian(&jacobian, t, &fields);
CPPUNIT_ASSERT_EQUAL(false, bc.needNewJacobian());
jacobian.complete();
+ DM dmMesh = mesh.dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
const PylithScalar* valsMatrixE = _data->valsJacobian;
- const int totalNumVertices = sieveMesh->depthStratum(0)->size();
+ const int totalNumVertices = vEnd - vStart;
const int sizeE = totalNumVertices * _data->spaceDim;
scalar_array valsE(sizeE);
const int spaceDim = _data->spaceDim;
@@ -338,19 +341,22 @@
} // for
#endif // DEBUGGING
- const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
- CPPUNIT_ASSERT(!jacobianSection.isNull());
- const PylithScalar* vals = jacobianSection->restrictSpace();
- const int size = jacobianSection->sizeWithBC();
+ PetscSection jacobianSection = jacobian.petscSection();
+ Vec jacobianVec = jacobian.localVector();
+ PetscScalar *vals;
+ PetscInt size;
+
+ CPPUNIT_ASSERT(jacobianSection);CPPUNIT_ASSERT(jacobianVec);
+ err = PetscSectionGetStorageSize(jacobianSection, &size);CHECK_PETSC_ERROR(err);
CPPUNIT_ASSERT_EQUAL(sizeE, size);
-
const PylithScalar tolerance = 1.0e-06;
- for (int i=0; i < size; ++i)
+ err = VecGetArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
+ for(int i = 0; i < size; ++i)
if (fabs(valsE[i]) > 1.0)
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, vals[i]/valsE[i], tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(valsE[i], vals[i], tolerance);
-
+ err = VecRestoreArray(jacobianVec, &vals);CHECK_PETSC_ERROR(err);
} // testIntegrateJacobianLumped
// ----------------------------------------------------------------------
@@ -413,52 +419,73 @@
fields->add("disp(t-dt)", "displacement");
fields->add("velocity(t)", "velocity");
fields->solutionName("dispIncr(t->t+dt)");
+
+ const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
topology::Field<topology::Mesh>& residual = fields->get("residual");
- const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- CPPUNIT_ASSERT(!vertices.isNull());
- residual.newSection(vertices, _data->spaceDim);
+ DM dmMesh = mesh->dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+
+ CPPUNIT_ASSERT(dmMesh);
+ err = DMComplexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
+ residual.newSection(pylith::topology::FieldBase::VERTICES_FIELD, _data->spaceDim);
residual.allocate();
residual.zero();
fields->copyLayout("residual");
- const int totalNumVertices = sieveMesh->depthStratum(0)->size();
- const int numMeshCells = sieveMesh->heightStratum(0)->size();
- int fieldSize = _data->spaceDim * totalNumVertices;
- const ALE::Obj<RealSection>& dispTIncrSection =
- fields->get("dispIncr(t->t+dt)").section();
- const ALE::Obj<RealSection>& dispTSection =
- fields->get("disp(t)").section();
- const ALE::Obj<RealSection>& dispTmdtSection =
- fields->get("disp(t-dt)").section();
- const ALE::Obj<RealSection>& velSection =
- fields->get("velocity(t)").section();
- CPPUNIT_ASSERT(!dispTIncrSection.isNull());
- CPPUNIT_ASSERT(!dispTSection.isNull());
- CPPUNIT_ASSERT(!dispTmdtSection.isNull());
- CPPUNIT_ASSERT(!velSection.isNull());
- const int offset = numMeshCells;
+ const int totalNumVertices = vEnd - vStart;
+ const int fieldSize = _data->spaceDim * totalNumVertices;
+ PetscSection dispTIncrSection = fields->get("dispIncr(t->t+dt)").petscSection();
+ Vec dispTIncrVec = fields->get("dispIncr(t->t+dt)").localVector();
+ PetscSection dispTSection = fields->get("disp(t)").petscSection();
+ Vec dispTVec = fields->get("disp(t)").localVector();
+ PetscSection dispTmdtSection = fields->get("disp(t-dt)").petscSection();
+ Vec dispTmdtVec = fields->get("disp(t-dt)").localVector();
+ PetscSection velSection = fields->get("velocity(t)").petscSection();
+ Vec velVec = fields->get("velocity(t)").localVector();
+ PetscScalar *dispTIncrArray, *dispTArray, *dispTmdtArray, *velArray;
const int spaceDim = _data->spaceDim;
const PylithScalar dt = _data->dt;
- scalar_array velVertex(spaceDim);
- for (int iVertex=0; iVertex < totalNumVertices; ++iVertex) {
- dispTIncrSection->updatePoint(iVertex+offset,
- &_data->fieldTIncr[iVertex*spaceDim]);
- dispTSection->updatePoint(iVertex+offset,
- &_data->fieldT[iVertex*spaceDim]);
- dispTmdtSection->updatePoint(iVertex+offset,
- &_data->fieldTmdt[iVertex*spaceDim]);
-
- for (int iDim=0; iDim < spaceDim; ++iDim)
- velVertex[iDim] = (_data->fieldTIncr[iVertex*spaceDim+iDim] +
- _data->fieldT[iVertex*spaceDim+iDim] -
- _data->fieldTmdt[iVertex*spaceDim+iDim]) / (2*dt);
- velSection->updatePoint(iVertex+offset, &velVertex[0]);
+ CPPUNIT_ASSERT(dispTIncrSection);CPPUNIT_ASSERT(dispTIncrVec);
+ CPPUNIT_ASSERT(dispTSection);CPPUNIT_ASSERT(dispTVec);
+ CPPUNIT_ASSERT(dispTmdtSection);CPPUNIT_ASSERT(dispTmdtVec);
+ CPPUNIT_ASSERT(velSection);CPPUNIT_ASSERT(velVec);
+ err = VecGetArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTmdtVec, &dispTmdtArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(dispTIncrSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTIncrSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ dispTIncrArray[off+d] = _data->fieldTIncr[(v - vStart)*spaceDim+d];
+ }
+ err = PetscSectionGetDof(dispTSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ dispTIncrArray[off+d] = _data->fieldT[(v - vStart)*spaceDim+d];
+ }
+ err = PetscSectionGetDof(dispTmdtSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTmdtSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ dispTIncrArray[off+d] = _data->fieldTmdt[(v - vStart)*spaceDim+d];
+ }
+ err = PetscSectionGetDof(velSection, v, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(velSection, v, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ velArray[off+d] = (_data->fieldTIncr[(v - vStart)*spaceDim+d] +
+ _data->fieldT[(v - vStart)*spaceDim+d] -
+ _data->fieldTmdt[(v - vStart)*spaceDim+d]) / (2*dt);
+ }
} // for
+ err = VecRestoreArray(dispTIncrVec, &dispTIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTmdtVec, &dispTmdtArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
} catch (const ALE::Exception& err) {
throw std::runtime_error(err.msg());
} // catch
More information about the CIG-COMMITS
mailing list