[cig-commits] [commit] knepley/upgrade-petsc-interface: Added C++ unit test for Field::copySubfield(). Fix bugs exposed. (e07734d)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 14 09:43:21 PST 2013
Repository : ssh://geoshell/pylith
On branch : knepley/upgrade-petsc-interface
Link : https://github.com/geodynamics/pylith/compare/cc6d236d1d32991191bebf4738dac67befb5fead...e07734d32415380de05c5d97fd96009562f76035
>---------------------------------------------------------------
commit e07734d32415380de05c5d97fd96009562f76035
Author: Brad Aagaard <baagaard at usgs.gov>
Date: Thu Nov 14 09:45:55 2013 -0800
Added C++ unit test for Field::copySubfield(). Fix bugs exposed.
Also added subfield arg to Field::setupSolnDof() to specify which
subfield to use.
>---------------------------------------------------------------
e07734d32415380de05c5d97fd96009562f76035
libsrc/pylith/topology/Field.cc | 61 ++++++++--------
libsrc/pylith/topology/Field.hh | 4 +-
modulesrc/topology/Field.i | 4 +-
unittests/libtests/topology/TestFieldMesh.cc | 100 +++++++++++++++++++++++++++
unittests/libtests/topology/TestFieldMesh.hh | 4 ++
5 files changed, 140 insertions(+), 33 deletions(-)
diff --git a/libsrc/pylith/topology/Field.cc b/libsrc/pylith/topology/Field.cc
index aa400a1..1694056 100644
--- a/libsrc/pylith/topology/Field.cc
+++ b/libsrc/pylith/topology/Field.cc
@@ -287,7 +287,8 @@ pylith::topology::Field::setupSolnChart(void)
// ----------------------------------------------------------------------
// Setup default DOF for solution.
void
-pylith::topology::Field::setupSolnDof(const int fiberDim)
+ pylith::topology::Field::setupSolnDof(const int fiberDim,
+ const char* subfieldName)
{ // setupSolnDof
PYLITH_METHOD_BEGIN;
@@ -297,7 +298,7 @@ pylith::topology::Field::setupSolnDof(const int fiberDim)
// :TODO: Update this to use discretization information after removing FIAT.
// :KLUDGE: Assume solution has DOF over vertices.
- const int indexDisp = subfieldMetadata("displacement").index;
+ const int subfieldIndex = subfieldMetadata(subfieldName).index;
PetscErrorCode err;
// Get range of vertices.
@@ -310,8 +311,8 @@ pylith::topology::Field::setupSolnDof(const int fiberDim)
for (PetscInt p = pStart; p < pEnd; ++p) {
err = PetscSectionSetDof(s, p, fiberDim);PYLITH_CHECK_ERROR(err);
- // Set DOF in displacement subfield
- err = PetscSectionSetFieldDof(s, p, indexDisp, fiberDim);PYLITH_CHECK_ERROR(err);
+ // Set DOF in subfield
+ err = PetscSectionSetFieldDof(s, p, subfieldIndex, fiberDim);PYLITH_CHECK_ERROR(err);
} // for
PYLITH_METHOD_END;
@@ -1350,17 +1351,17 @@ pylith::topology::Field::copySubfield(const Field& field,
_metadata["default"] = subfieldMetadata;
label(subfieldMetadata.label.c_str()); // Use method to insure propagation to subsidiary objects
- const PetscSection& fieldSection = field.petscSection();
- const PetscSection& subfieldSection = this->petscSection();
-
// Check compatibility of sections
const int srcSize = field.chartSize();
const int dstSize = chartSize();
- if (srcSize < dstSize) {
+ if (dstSize < srcSize) {
_extractSubfield(field, name);
} // if
assert(_localVec && field._localVec);
+ const PetscSection& fieldSection = field.petscSection();
+ const PetscSection& subfieldSection = this->petscSection();
+
PetscInt pStart, pEnd;
err = PetscSectionGetChart(subfieldSection, &pStart, &pEnd);PYLITH_CHECK_ERROR(err);
@@ -1404,12 +1405,12 @@ pylith::topology::Field::_extractSubfield(const Field& field,
const int subfieldIndex = subfieldMetadata.index;assert(subfieldIndex >= 0);
PetscErrorCode err;
- PetscDM subfieldDM = NULL;
PetscIS subfieldIS = NULL;
const int numSubfields = 1;
int indicesSubfield[1];
indicesSubfield[0] = subfieldIndex;
- err = DMCreateSubDM(_dm, numSubfields, indicesSubfield, &subfieldIS, &subfieldDM);PYLITH_CHECK_ERROR(err);assert(subfieldDM);assert(subfieldIS);
+ err = DMCreateSubDM(field.dmMesh(), numSubfields, indicesSubfield, &subfieldIS, &_dm);PYLITH_CHECK_ERROR(err);assert(_dm);
+ err = ISDestroy(&subfieldIS);PYLITH_CHECK_ERROR(err);
err = DMCreateLocalVector(_dm, &_localVec);PYLITH_CHECK_ERROR(err);
err = DMCreateGlobalVector(_dm, &_globalVec);PYLITH_CHECK_ERROR(err);
@@ -1417,35 +1418,33 @@ pylith::topology::Field::_extractSubfield(const Field& field,
// Setup section
const PetscSection& fieldSection = field.petscSection();
PetscSection subfieldSection = NULL;
- err = DMGetDefaultSection(this->_dm, &subfieldSection);PYLITH_CHECK_ERROR(err);
- err = DMSetDefaultGlobalSection(this->_dm, NULL);PYLITH_CHECK_ERROR(err);
-
- const PetscInt* points = NULL;
- PetscInt numPoints = 0;
- err = ISGetSize(subfieldIS, &numPoints);PYLITH_CHECK_ERROR(err);
- err = ISGetIndices(subfieldIS, &points);PYLITH_CHECK_ERROR(err);
- const PetscInt pStart = (numPoints > 0) ? points[0] : 0;
- const PetscInt pEnd = (numPoints > 0) ? points[numPoints-1] : 0;
+ err = DMGetDefaultSection(_dm, &subfieldSection);PYLITH_CHECK_ERROR(err);
+ err = DMSetDefaultGlobalSection(_dm, NULL);PYLITH_CHECK_ERROR(err);
+
+ PetscInt pStart = -1, pEnd = -1;
+ err = PetscSectionGetChart(fieldSection, &pStart, &pEnd);PYLITH_CHECK_ERROR(err);
err = PetscSectionSetChart(subfieldSection, pStart, pEnd);PYLITH_CHECK_ERROR(err);
- for (PetscInt i=0; i < numPoints; ++i) {
- const PetscInt point = points[i];
+ for (PetscInt p=pStart; p < pEnd; ++p) {
PetscInt dof;
- err = PetscSectionGetFieldDof(fieldSection, point, subfieldIndex, &dof);PYLITH_CHECK_ERROR(err);
- err = PetscSectionSetDof(subfieldSection, point, dof);PYLITH_CHECK_ERROR(err);
+ err = PetscSectionGetFieldDof(fieldSection, p, subfieldIndex, &dof);PYLITH_CHECK_ERROR(err);
+ if (dof > 0) {
+ err = PetscSectionSetDof(subfieldSection, p, dof);PYLITH_CHECK_ERROR(err);
- err = PetscSectionGetFieldConstraintDof(fieldSection, point, subfieldIndex, &dof);PYLITH_CHECK_ERROR(err);
- err = PetscSectionSetConstraintDof(subfieldSection, point, dof);PYLITH_CHECK_ERROR(err);
+ err = PetscSectionGetFieldConstraintDof(fieldSection, p, subfieldIndex, &dof);PYLITH_CHECK_ERROR(err);
+ err = PetscSectionSetConstraintDof(subfieldSection, p, dof);PYLITH_CHECK_ERROR(err);
+ } // if
} // for
- this->allocate();
+ allocate();
- for (PetscInt i=0; i < numPoints; ++i) {
- const PetscInt point = points[i];
+ for (PetscInt p=pStart; p < pEnd; ++p) {
PetscInt dof;
const PetscInt* indices = NULL;
- err = PetscSectionGetConstraintDof(subfieldSection, point, &dof);PYLITH_CHECK_ERROR(err);
- err = PetscSectionGetFieldConstraintIndices(fieldSection, point, subfieldIndex, &indices);PYLITH_CHECK_ERROR(err);
- err = PetscSectionSetConstraintIndices(subfieldSection, point, indices);PYLITH_CHECK_ERROR(err);
+ err = PetscSectionGetConstraintDof(subfieldSection, p, &dof);PYLITH_CHECK_ERROR(err);
+ if (dof > 0) {
+ err = PetscSectionGetFieldConstraintIndices(fieldSection, p, subfieldIndex, &indices);PYLITH_CHECK_ERROR(err);
+ err = PetscSectionSetConstraintIndices(subfieldSection, p, indices);PYLITH_CHECK_ERROR(err);
+ } // if
} // for
PYLITH_METHOD_END;
diff --git a/libsrc/pylith/topology/Field.hh b/libsrc/pylith/topology/Field.hh
index d5ae449..99a7126 100644
--- a/libsrc/pylith/topology/Field.hh
+++ b/libsrc/pylith/topology/Field.hh
@@ -222,8 +222,10 @@ public :
/** Set default DOF for solution.
*
* @param fiberDim Total number of components in solution.
+ * @param subfieldName Name of subfield for DOF.
*/
- void setupSolnDof(const int fiberDim);
+ void setupSolnDof(const int fiberDim,
+ const char* subfieldName ="displacement");
/** Create PETSc section and set chart and fiber dimesion for a list
* of points.
diff --git a/modulesrc/topology/Field.i b/modulesrc/topology/Field.i
index e55728e..d4b9f4f 100644
--- a/modulesrc/topology/Field.i
+++ b/modulesrc/topology/Field.i
@@ -127,8 +127,10 @@ namespace pylith {
/** Set default DOF for solution.
*
* @param fiberDim Total number of components in solution.
+ * @param subfieldName Name of subfield for DOF.
*/
- void setupSolnDof(const int fiberDim);
+ void setupSolnDof(const int fiberDim,
+ const char* subfieldName ="displacement");
/** Create PETSc section and set chart and fiber dimesion.
*
diff --git a/unittests/libtests/topology/TestFieldMesh.cc b/unittests/libtests/topology/TestFieldMesh.cc
index ebf3021..3e78cc3 100644
--- a/unittests/libtests/topology/TestFieldMesh.cc
+++ b/unittests/libtests/topology/TestFieldMesh.cc
@@ -825,6 +825,106 @@ pylith::topology::TestFieldMesh::testCopy(void)
} // testCopy
// ----------------------------------------------------------------------
+// Test copySubfield().
+void
+pylith::topology::TestFieldMesh::testCopySubfield(void)
+{ // testCopy
+ PYLITH_METHOD_BEGIN;
+
+ const int fiberDimA = 1;
+ const char* fieldA = "one";
+ const int fiberDimB = 2;
+ const char* fieldB = "two";
+ const int fiberDim = fiberDimA + fiberDimB;
+
+ const PylithScalar scale = 2.0;
+ const int npoints = 4;
+ const PylithScalar valuesNondim[npoints*fiberDim] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+ const PylithScalar valuesANondim[npoints*fiberDimA] = {
+ 1.1,
+ 1.2,
+ 1.3,
+ 1.4,
+ };
+ const PylithScalar valuesBNondim[npoints*fiberDimB] = {
+ 2.2, 3.3,
+ 2.3, 3.4,
+ 2.4, 3.5,
+ 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+
+ PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
+ Stratum depthStratum(dmMesh, Stratum::DEPTH, 0);
+ const PetscInt vStart = depthStratum.begin();
+ const PetscInt vEnd = depthStratum.end();
+
+ Field fieldSrc(mesh);
+ { // Setup source field
+ fieldSrc.label("solution");
+ fieldSrc.subfieldAdd("one", fiberDimA);
+ fieldSrc.subfieldAdd("two", fiberDimB);
+ fieldSrc.subfieldsSetup();
+ fieldSrc.newSection(Field::VERTICES_FIELD, fiberDim);
+ fieldSrc.subfieldSetDof("one", Field::VERTICES_FIELD, fiberDimA);
+ fieldSrc.subfieldSetDof("two", Field::VERTICES_FIELD, fiberDimB);
+ fieldSrc.allocate();
+
+ VecVisitorMesh fieldVisitor(fieldSrc);
+ PetscScalar* fieldArray = fieldVisitor.localArray();
+ for (PetscInt v = vStart, i = 0; v < vEnd; ++v) {
+ const PetscInt off = fieldVisitor.sectionOffset(v);
+ CPPUNIT_ASSERT_EQUAL(fiberDim, fieldVisitor.sectionDof(v));
+ for (PetscInt d = 0; d < fiberDim; ++d)
+ fieldArray[off+d] = valuesNondim[i++];
+ } // for
+ } // Setup source field
+
+ { // Test with preallocated field
+ Field field(mesh);
+ field.newSection(Field::VERTICES_FIELD, fiberDimB);
+ field.allocate();
+ field.copySubfield(fieldSrc, "two");
+
+ VecVisitorMesh fieldVisitor(field);
+ const PetscScalar* fieldArray = fieldVisitor.localArray();
+ const PylithScalar tolerance = 1.0e-6;
+ for (PetscInt v = vStart, i = 0; v < vEnd; ++v) {
+ const PetscInt off = fieldVisitor.sectionOffset(v);
+ CPPUNIT_ASSERT_EQUAL(fiberDimB, fieldVisitor.sectionDof(v));
+ for (PetscInt d = 0; d < fiberDimB; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesBNondim[i++], fieldArray[off+d], tolerance);
+ } // for
+ } // for
+ } // Test with preallocated field
+
+ { // Test with unallocated field
+ Field field(mesh);
+ field.copySubfield(fieldSrc, "one");
+
+ VecVisitorMesh fieldVisitor(field);
+ const PetscScalar* fieldArray = fieldVisitor.localArray();
+ const PylithScalar tolerance = 1.0e-6;
+ for (PetscInt v = vStart, i = 0; v < vEnd; ++v) {
+ const PetscInt off = fieldVisitor.sectionOffset(v);
+ CPPUNIT_ASSERT_EQUAL(fiberDimA, fieldVisitor.sectionDof(v));
+ for (PetscInt d = 0; d < fiberDimA; ++d) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesANondim[i++], fieldArray[off+d], tolerance);
+ } // for
+ } // for
+ } // Test with unallocated field
+
+ PYLITH_METHOD_END;
+} // testCopySubfield
+
+// ----------------------------------------------------------------------
// Test operator+=().
void
pylith::topology::TestFieldMesh::testOperatorAdd(void)
diff --git a/unittests/libtests/topology/TestFieldMesh.hh b/unittests/libtests/topology/TestFieldMesh.hh
index 1610231..436dcc7 100644
--- a/unittests/libtests/topology/TestFieldMesh.hh
+++ b/unittests/libtests/topology/TestFieldMesh.hh
@@ -67,6 +67,7 @@ class pylith::topology::TestFieldMesh : public CppUnit::TestFixture
CPPUNIT_TEST( testZeroAll );
CPPUNIT_TEST( testComplete );
CPPUNIT_TEST( testCopy );
+ CPPUNIT_TEST( testCopySubfield );
CPPUNIT_TEST( testOperatorAdd );
CPPUNIT_TEST( testDimensionalize );
CPPUNIT_TEST( testView );
@@ -149,6 +150,9 @@ public :
/// Test copy().
void testCopy(void);
+ /// Test copySubfield().
+ void testCopySubfield(void);
+
/// Test operator+=().
void testOperatorAdd(void);
More information about the CIG-COMMITS
mailing list