[cig-commits] r20619 - short/3D/PyLith/trunk/libsrc/pylith/topology
knepley at geodynamics.org
knepley at geodynamics.org
Tue Aug 21 18:09:35 PDT 2012
Author: knepley
Date: 2012-08-21 18:09:35 -0700 (Tue, 21 Aug 2012)
New Revision: 20619
Modified:
short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.icc
Log:
Added DM version to Fields
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.hh 2012-08-21 23:31:18 UTC (rev 20618)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.hh 2012-08-22 01:09:35 UTC (rev 20619)
@@ -85,6 +85,16 @@
const pylith::topology::FieldBase::DomainEnum domain,
const int fiberDim);
+ /** Add field.
+ *
+ * @param name Name of field.
+ * @param numFields The number of fields to select
+ * @param fields Set of field numbers for this vector
+ */
+ void add(const char* name,
+ const int numFields,
+ const int fields[]);
+
/** Delete field.
*
* @param name Name of field.
@@ -109,6 +119,12 @@
*/
field_type& get(const char* name);
+ /** Get DM field.
+ *
+ * @param name Name of field.
+ */
+ DM getDM(const char* name);
+
/** Copy layout to other fields.
*
* @param name Name of field to use as template for layout.
@@ -133,11 +149,13 @@
protected :
typedef std::map< std::string, field_type* > map_type;
+ typedef std::map< std::string, DM> dm_map_type;
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
map_type _fields;
+ dm_map_type _dmFields; /* These DMs will share the underlying mesh with _mesh */
const typename field_type::Mesh& _mesh;
// NOT IMPLEMENTED //////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.icc 2012-08-21 23:31:18 UTC (rev 20618)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.icc 2012-08-22 01:09:35 UTC (rev 20619)
@@ -22,6 +22,7 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
+#include <pylith/utils/petscerror.h> // USES CHECK_PETSC_ERROR
// ----------------------------------------------------------------------
// Default constructor.
@@ -50,6 +51,11 @@
for (typename map_type::iterator iter=begin; iter != end; ++iter) {
delete iter->second; iter->second = 0;
} // for
+ const typename dm_map_type::iterator dmBegin = _dmFields.begin();
+ const typename dm_map_type::iterator dmEnd = _dmFields.end();
+ for (typename dm_map_type::iterator iter=dmBegin; iter != dmEnd; ++iter) {
+ PetscErrorCode err = DMDestroy(&iter->second);CHECK_PETSC_ERROR(err);
+ } // for
} // deallocate
// ----------------------------------------------------------------------
@@ -58,7 +64,9 @@
bool
pylith::topology::Fields<field_type>::hasField(const char* name) const
{ // hasField
+ typename dm_map_type::const_iterator dmIter = _dmFields.find(name);
typename map_type::const_iterator iter = _fields.find(name);
+ assert((dmIter != _dmFields.end()) == (iter != _fields.end()));
return iter != _fields.end();
} // hasField
@@ -78,6 +86,10 @@
_fields[name] = new field_type(_mesh);
_fields[name]->label(label);
+ if (_mesh.dmMesh()) {
+ PetscErrorCode err = PetscObjectReference((PetscObject) _mesh.dmMesh());CHECK_PETSC_ERROR(err);
+ _dmFields[name] = _mesh.dmMesh();
+ }
} // add
// ----------------------------------------------------------------------
@@ -100,9 +112,66 @@
_fields[name] = new field_type(_mesh);
_fields[name]->label(label);
_fields[name]->newSection(domain, fiberDim);
+ if (_mesh.dmMesh()) {
+ DM dm = _mesh.dmMesh(), dmNew;
+ PetscSection sectionNew;
+ PetscInt pStart, pEnd, p;
+ PetscErrorCode err;
+
+ err = DMComplexClone(dm, &dmNew);CHECK_PETSC_ERROR(err);
+ err = PetscSectionCreate(_mesh.comm(), §ionNew);CHECK_PETSC_ERROR(err);
+ switch(domain) {
+ case pylith::topology::FieldBase::VERTICES_FIELD:
+ err = DMComplexGetDepthStratum(dmNew, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
+ case pylith::topology::FieldBase::CELLS_FIELD:
+ err = DMComplexGetHeightStratum(dmNew, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
+ default:
+ throw std::runtime_error("Invalid domain type");
+ }
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ PetscInt value;
+
+ err = DMComplexGetLabelValue(dmNew, label, p, &value);CHECK_PETSC_ERROR(err);
+ if (value) {
+ err = PetscSectionSetDof(sectionNew, p, fiberDim);CHECK_PETSC_ERROR(err);
+ }
+ }
+ err = PetscSectionSetUp(sectionNew);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultSection(dmNew, sectionNew);CHECK_PETSC_ERROR(err);
+ _dmFields[name] = dmNew;
+ }
} // add
// ----------------------------------------------------------------------
+// Add field.
+template<typename field_type>
+void
+pylith::topology::Fields<field_type>::add(
+ const char* name,
+ const int numFields,
+ const int fields[])
+{ // add
+ if (hasField(name)) {
+ std::ostringstream msg;
+ msg << "Could not add field '" << name
+ << "' to fields manager, because it already exists.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ if (_mesh.dmMesh()) {
+ DM dm = _mesh.dmMesh(), dmNew;
+ PetscSection sectionNew;
+ PetscInt pStart, pEnd, p;
+ PetscErrorCode err;
+
+ err = DMCreateSubDM(dm, numFields, fields, PETSC_NULL, &dmNew);CHECK_PETSC_ERROR(err);
+ _dmFields[name] = dmNew;
+ }
+} // add
+
+// ----------------------------------------------------------------------
// Delete field.
template<typename field_type>
void
@@ -117,6 +186,7 @@
} // if
delete iter->second; iter->second = 0;
_fields.erase(name);
+ _dmFields.erase(name);
} // del
// ----------------------------------------------------------------------
@@ -160,6 +230,22 @@
} // if
return *iter->second;
} // get
+
+// ----------------------------------------------------------------------
+// Get field.
+template<typename field_type>
+DM
+pylith::topology::Fields<field_type>::getDM(const char* name)
+{ // get
+ typename dm_map_type::iterator iter = _dmFields.find(name);
+ if (iter == _dmFields.end()) {
+ std::ostringstream msg;
+ msg << "Could not find field '" << name
+ << "' in fields manager for retrieval.";
+ throw std::runtime_error(msg.str());
+ } // if
+ return *iter->second;
+} // get
// ----------------------------------------------------------------------
// Copy layout to other fields.
More information about the CIG-COMMITS
mailing list