[cig-commits] r20655 - in short/3D/PyLith/trunk: libsrc/pylith/bc libsrc/pylith/topology unittests/libtests/bc
knepley at geodynamics.org
knepley at geodynamics.org
Fri Aug 31 18:11:25 PDT 2012
Author: knepley
Date: 2012-08-31 18:11:25 -0700 (Fri, 31 Aug 2012)
New Revision: 20655
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
Log:
Added in some DirichletBC support
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2012-08-31 23:43:50 UTC (rev 20654)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/DirichletBC.cc 2012-09-01 01:11:25 UTC (rev 20655)
@@ -85,25 +85,46 @@
const ALE::Obj<RealSection>& section = field.section();
assert(!section.isNull());
+ PetscSection sectionP = field.petscSection();
+ PetscErrorCode err;
+ assert(sectionP);
+
// Set constraints in field
const int numPoints = _points.size();
_offsetLocal.resize(numPoints);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const int fiberDim = section->getFiberDimension(_points[iPoint]);
- const int curNumConstraints =
- section->getConstraintDimension(_points[iPoint]);
+ const PetscInt point = _points[iPoint];
+ const int fiberDim = section->getFiberDimension(point);
+ const int curNumConstraints = section->getConstraintDimension(point);
if (curNumConstraints + numFixedDOF > fiberDim) {
std::ostringstream msg;
msg
<< "Found overly constrained point while setting up constraints for\n"
<< "DirichletBC boundary condition '" << _label << "'.\n"
- << "Number of DOF at point " << _points[iPoint] << " is " << fiberDim
+ << "Number of DOF at point " << point << " is " << fiberDim
<< "\nand number of attempted constraints is "
<< curNumConstraints+numFixedDOF << ".";
throw std::runtime_error(msg.str());
} // if
_offsetLocal[iPoint] = curNumConstraints;
- section->addConstraintDimension(_points[iPoint], numFixedDOF);
+ section->addConstraintDimension(point, numFixedDOF);
+ PetscInt dof, cdof;
+
+ err = PetscSectionGetDof(sectionP, point, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintDof(sectionP, point, &cdof);CHECK_PETSC_ERROR(err);
+ if (cdof + numFixedDOF > dof) {
+ std::ostringstream msg;
+ msg
+ << "Found overly constrained point while setting up constraints for\n"
+ << "DirichletBC boundary condition '" << _label << "'.\n"
+ << "Number of DOF at point " << point << " is " << dof
+ << "\nand number of attempted constraints is " << cdof+numFixedDOF << ".";
+ throw std::runtime_error(msg.str());
+ } // if
+ _offsetLocal[iPoint] = cdof;
+ err = PetscSectionAddConstraintDof(sectionP, point, numFixedDOF);CHECK_PETSC_ERROR(err);
+ // We should be specifying what field the BC is for
+ //err = PetscSectionAddConstraintFieldDof(sectionP, point, field, numFixedDOF);CHECK_PETSC_ERROR(err);
} // for
// We only worry about the conventional DOF in a split field associated with
@@ -129,6 +150,10 @@
const ALE::Obj<RealSection>& section = field.section();
assert(!section.isNull());
+ PetscSection sectionP = field.petscSection();
+ PetscErrorCode err;
+ assert(sectionP);
+
const int numPoints = _points.size();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
@@ -136,28 +161,35 @@
// Get list of currently constrained DOF
const int* curFixedDOF = section->getConstraintDof(point);
const int numTotalConstrained = section->getConstraintDimension(point);
+ PetscInt cdof, *cInd;
+ err = PetscSectionGetConstraintDof(sectionP, point, &cdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetConstraintIndices(sectionP, point, &cInd);CHECK_PETSC_ERROR(err);
+
// Create array holding all constrained DOF
int_array allFixedDOF(curFixedDOF, numTotalConstrained);
+ int_array allCInd(cInd, cdof);
// Verify other BC has not already constrained DOF
const int numPrevious = _offsetLocal[iPoint];
for (int iDOF=0; iDOF < numPrevious; ++iDOF)
for (int jDOF=0; jDOF < numFixedDOF; ++jDOF)
- if (allFixedDOF[iDOF] == _bcDOF[jDOF]) {
- std::ostringstream msg;
- msg << "Found multiple constraints on degrees of freedom at\n"
- << "point while setting up constraints for DirichletBC\n"
- << "boundary condition '" << _label << "'.\n"
+ if ((allFixedDOF[iDOF] == _bcDOF[jDOF]) || (allCInd[iDOF] == _bcDOF[jDOF])) {
+ std::ostringstream msg;
+ msg << "Found multiple constraints on degrees of freedom at\n"
+ << "point while setting up constraints for DirichletBC\n"
+ << "boundary condition '" << _label << "'.\n"
<< "Degree of freedom " << _bcDOF[jDOF]
- << " is already constrained by another Dirichlet BC.";
- throw std::runtime_error(msg.str());
- } // if
+ << " is already constrained by another Dirichlet BC.";
+ throw std::runtime_error(msg.str());
+ } // if
// Add in the ones for this DirichletBC BC
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF) {
assert(_offsetLocal[iPoint]+iDOF < numTotalConstrained);
+ assert(_offsetLocal[iPoint]+iDOF < cdof);
allFixedDOF[_offsetLocal[iPoint]+iDOF] = _bcDOF[iDOF];
+ allCInd[_offsetLocal[iPoint]+iDOF] = _bcDOF[iDOF];
} // for
// Fill in rest of values not yet set (will be set by
@@ -168,13 +200,19 @@
assert(iDOF < numTotalConstrained);
allFixedDOF[iDOF] = 999;
} // for
+ for (int iDOF=_offsetLocal[iPoint]+numFixedDOF; iDOF < cdof; ++iDOF) {
+ allCInd[iDOF] = 999;
+ } // for
// Sort list of constrained DOF
// I need these sorted for my update algorithms to work properly
std::sort(&allFixedDOF[0], &allFixedDOF[numTotalConstrained]);
+ std::sort(&allCInd[0], &allCInd[cdof]);
// Update list of constrained DOF
section->setConstraintDof(point, &allFixedDOF[0]);
+ err = PetscSectionSetConstraintIndices(sectionP, point, &allCInd[0]);CHECK_PETSC_ERROR(err);
+ //err = PetscSectionSetConstraintFieldIndices(sectionP, point, field, &fieldCInd[0]);CHECK_PETSC_ERROR(err);
} // for
// We only worry about the conventional DOF in a split field associated with
@@ -220,19 +258,34 @@
const int fiberDimension =
(numPoints > 0) ? fieldSection->getFiberDimension(_points[0]) : 0;
scalar_array fieldVertex(fiberDimension);
+
+ PetscSection fieldSectionP = field.petscSection();
+ Vec localVec = field.localVector();
+ PetscErrorCode err;
+ assert(fieldSectionP);
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type p_bc = _points[iPoint];
+ assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
fieldSection->restrictPoint(p_bc, &fieldVertex[0], fieldVertex.size());
+ PetscInt dof, off;
- assert(parametersFiberDim == parametersSection->getFiberDimension(p_bc));
+ err = PetscSectionGetDof(fieldSectionP, p_bc, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSectionP, p_bc, &off);CHECK_PETSC_ERROR(err);
+
const PylithScalar* parametersVertex = parametersSection->restrictPoint(p_bc);
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
fieldVertex[_bcDOF[iDOF]] = parametersVertex[valueIndex+iDOF];
fieldSection->updatePointAll(p_bc, &fieldVertex[0]);
+ PetscScalar *array;
+
+ err = VecGetArray(localVec, &array);CHECK_PETSC_ERROR(err);
+ for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
+ array[_bcDOF[iDOF]] = parametersVertex[valueIndex+iDOF];
+ err = VecRestoreArray(localVec, &array);CHECK_PETSC_ERROR(err);
} // for
} // setField
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-08-31 23:43:50 UTC (rev 20654)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-09-01 01:11:25 UTC (rev 20655)
@@ -49,6 +49,7 @@
_dm = PETSC_NULL;
}
_globalVec = PETSC_NULL;
+ _localVec = PETSC_NULL;
} // constructor
// ----------------------------------------------------------------------
@@ -73,6 +74,7 @@
_dm = PETSC_NULL;
}
_globalVec = PETSC_NULL;
+ _localVec = PETSC_NULL;
} // constructor
// ----------------------------------------------------------------------
@@ -95,6 +97,7 @@
}
err = DMCreateSubDM(_mesh.dmMesh(), numFields, fields, PETSC_NULL, &_dm);CHECK_PETSC_ERROR(err);
_globalVec = PETSC_NULL;
+ _localVec = PETSC_NULL;
} // constructor
// ----------------------------------------------------------------------
@@ -103,6 +106,7 @@
pylith::topology::Field<mesh_type, section_type>::~Field(void)
{ // destructor
deallocate();
+ PetscErrorCode err = DMDestroy(&_dm);CHECK_PETSC_ERROR(err);
} // destructor
// ----------------------------------------------------------------------
@@ -124,7 +128,7 @@
} // for
_scatters.clear();
err = VecDestroy(&_globalVec);CHECK_PETSC_ERROR(err);
- err = DMDestroy(&_dm);CHECK_PETSC_ERROR(err);
+ err = VecDestroy(&_localVec);CHECK_PETSC_ERROR(err);
} // deallocate
// ----------------------------------------------------------------------
@@ -207,6 +211,7 @@
const int fiberDim)
{ // newSection
typedef typename mesh_type::SieveMesh::point_type point_type;
+ PetscErrorCode err;
// Clear memory
clear();
@@ -232,8 +237,24 @@
*std::max_element(points->begin(), points->end());
_section->setChart(chart_type(pointMin, pointMax+1));
_section->setFiberDimension(points, fiberDim);
- } else // Create empty chart
+
+ if (_dm) {
+ PetscSection s;
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, pointMin, pointMax+1);CHECK_PETSC_ERROR(err);
+
+ for(typename label_sequence::const_iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
+ err = PetscSectionSetDof(s, *p_iter, fiberDim);CHECK_PETSC_ERROR(err);
+ }
+ }
+ } else {// Create empty chart
_section->setChart(chart_type(0, 0));
+ if (_dm) {
+ PetscSection s;
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, 0, 0);CHECK_PETSC_ERROR(err);
+ }
+ }
logger.stagePop();
} // newSection
@@ -247,6 +268,7 @@
const int fiberDim)
{ // newSection
typedef typename mesh_type::SieveMesh::point_type point_type;
+ PetscErrorCode err;
// Clear memory
clear();
@@ -271,8 +293,24 @@
_section->setChart(chart_type(pointMin, pointMax+1));
for (int i=0; i < npts; ++i)
_section->setFiberDimension(points[i], fiberDim);
- } else // create empty chart
+
+ if (_dm) {
+ PetscSection s;
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, pointMin, pointMax+1);CHECK_PETSC_ERROR(err);
+
+ for (int i=0; i < npts; ++i) {
+ err = PetscSectionSetDof(s, points[i], fiberDim);CHECK_PETSC_ERROR(err);
+ }
+ }
+ } else { // create empty chart
_section->setChart(chart_type(0, 0));
+ if (_dm) {
+ PetscSection s;
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, 0, 0);CHECK_PETSC_ERROR(err);
+ }
+ }
logger.stagePop();
} // newSection
@@ -611,6 +649,7 @@
PetscErrorCode err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
err = PetscSectionSetUp(s);CHECK_PETSC_ERROR(err);
err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+ err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
}
logger.stagePop();
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-08-31 23:43:50 UTC (rev 20654)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-09-01 01:11:25 UTC (rev 20655)
@@ -105,7 +105,25 @@
* @returns DMComplex section.
*/
void dmSection(PetscSection *s, Vec *v) const;
+
+ /** Get PetscSection.
+ *
+ * @returns PetscSection.
+ */
+ PetscSection petscSection() const;
+ /** Get local Vec.
+ *
+ * @returns local Vec.
+ */
+ Vec localVector() const;
+
+ /** Get global Vec.
+ *
+ * @returns global Vec.
+ */
+ Vec globalVector() const;
+
/** Get mesh associated with field.
*
* @returns Finite-element mesh.
@@ -461,6 +479,7 @@
/* New construction */
DM _dm; /* Holds the PetscSection */
Vec _globalVec;
+ Vec _localVec;
std::map<std::string, int> _tmpFields;
// NOT IMPLEMENTED //////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-08-31 23:43:50 UTC (rev 20654)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-09-01 01:11:25 UTC (rev 20655)
@@ -30,6 +30,34 @@
return _section;
}
+// Get PetscSection.
+template<typename mesh_type, typename section_type>
+inline
+PetscSection
+pylith::topology::Field<mesh_type, section_type>::petscSection(void) const {
+ PetscSection s = PETSC_NULL;
+ if (_dm) {
+ PetscErrorCode err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ }
+ return s;
+}
+
+// Get local vector.
+template<typename mesh_type, typename section_type>
+inline
+Vec
+pylith::topology::Field<mesh_type, section_type>::localVector(void) const {
+ return _localVec;
+}
+
+// Get global vector.
+template<typename mesh_type, typename section_type>
+inline
+Vec
+pylith::topology::Field<mesh_type, section_type>::globalVector(void) const {
+ return _globalVec;
+}
+
// Get mesh associated with field.
template<typename mesh_type, typename section_type>
inline
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2012-08-31 23:43:50 UTC (rev 20654)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestDirichletBC.cc 2012-09-01 01:11:25 UTC (rev 20655)
@@ -170,6 +170,7 @@
topology::Field<topology::Mesh> field(mesh);
field.newSection(vertices, fiberDim);
field.splitDefault();
+
const ALE::Obj<RealSection>& fieldSection = field.section();
CPPUNIT_ASSERT(!fieldSection.isNull());
More information about the CIG-COMMITS
mailing list