[cig-commits] r20632 - in short/3D/PyLith/trunk: libsrc/pylith/topology pylith/problems
knepley at geodynamics.org
knepley at geodynamics.org
Sat Aug 25 12:54:25 PDT 2012
Author: knepley
Date: 2012-08-25 12:54:25 -0700 (Sat, 25 Aug 2012)
New Revision: 20632
Modified:
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/libsrc/pylith/topology/Fields.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.icc
short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
Changed to a model where Field holds a DM, put some prelim stuff in Formulation
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-08-25 15:24:08 UTC (rev 20631)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2012-08-25 19:54:25 UTC (rev 20632)
@@ -36,10 +36,19 @@
pylith::topology::Field<mesh_type, section_type>::Field(const mesh_type& mesh) :
_mesh(mesh)
{ // constructor
- _metadata.label = "unknown";
- _metadata.vectorFieldType = OTHER;
- _metadata.scale = 1.0;
- _metadata.dimsOkay = false;
+ _metadata["default"].label = "unknown";
+ _metadata["default"].vectorFieldType = OTHER;
+ _metadata["default"].scale = 1.0;
+ _metadata["default"].dimsOkay = false;
+ if (mesh.dmMesh()) {
+ PetscSection s;
+ PetscErrorCode err = DMComplexClone(mesh.dmMesh(), &_dm);CHECK_PETSC_ERROR(err);
+ err = PetscSectionCreate(mesh.comm(), &s);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultSection(_dm, s);CHECK_PETSC_ERROR(err);
+ } else {
+ _dm = PETSC_NULL;
+ }
+ _globalVec = PETSC_NULL;
} // constructor
// ----------------------------------------------------------------------
@@ -48,14 +57,48 @@
pylith::topology::Field<mesh_type, section_type>::Field(const mesh_type& mesh,
const ALE::Obj<section_type>& section,
const Metadata& metadata) :
- _metadata(metadata),
_mesh(mesh),
_section(section)
{ // constructor
assert(!section.isNull());
+ _metadata["default"] = metadata;
+ if (mesh.dmMesh()) {
+ Vec lv;
+ PetscSection s;
+ PetscErrorCode err = DMComplexClone(mesh.dmMesh(), &_dm);CHECK_PETSC_ERROR(err);
+ this->dmSection(&s, &lv);
+ err = VecDestroy(&lv);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultSection(_dm, s);CHECK_PETSC_ERROR(err);
+ } else {
+ _dm = PETSC_NULL;
+ }
+ _globalVec = PETSC_NULL;
} // constructor
// ----------------------------------------------------------------------
+// Constructor with field and subfields
+template<typename mesh_type, typename section_type>
+pylith::topology::Field<mesh_type, section_type>::Field(const Field& src,
+ int numFields,
+ const int fields[]) :
+ _mesh(src._mesh),
+ _section(PETSC_NULL)
+{ // constructor
+ PetscSection s;
+ PetscErrorCode err;
+
+ err = DMGetDefaultSection(src._dm, &s);CHECK_PETSC_ERROR(err);
+ for(PetscInt f = 0; f < numFields; ++f) {
+ const char *name;
+
+ err = PetscSectionGetFieldName(s, fields[f], &name);CHECK_PETSC_ERROR(err);
+ _metadata[name] = src._metadata[name];
+ }
+ err = DMCreateSubDM(_mesh.dmMesh(), numFields, fields, PETSC_NULL, &_dm);CHECK_PETSC_ERROR(err);
+ _globalVec = PETSC_NULL;
+} // constructor
+
+// ----------------------------------------------------------------------
// Destructor.
template<typename mesh_type, typename section_type>
pylith::topology::Field<mesh_type, section_type>::~Field(void)
@@ -81,6 +124,8 @@
err = VecDestroy(&s_iter->second.scatterVec);CHECK_PETSC_ERROR(err);
} // for
_scatters.clear();
+ err = VecDestroy(&_globalVec);CHECK_PETSC_ERROR(err);
+ err = DMDestroy(&_dm);CHECK_PETSC_ERROR(err);
} // deallocate
// ----------------------------------------------------------------------
@@ -89,7 +134,7 @@
void
pylith::topology::Field<mesh_type, section_type>::label(const char* value)
{ // label
- _metadata.label = value;
+ _metadata["default"].label = value;
if (!_section.isNull()) {
_section->setName(value);
} // if
@@ -148,7 +193,7 @@
_section = new section_type(_mesh.comm(), _mesh.debug());
assert(!_section.isNull());
- _section->setName(_metadata.label);
+ _section->setName(_metadata["default"].label);
logger.stagePop();
} // newSection
@@ -172,14 +217,14 @@
if (fiberDim < 0) {
std::ostringstream msg;
- msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata.label
+ msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata["default"].label
<< "' must be nonnegative.";
throw std::runtime_error(msg.str());
} // if
_section = new section_type(_mesh.comm(), _mesh.debug());
assert(!_section.isNull());
- _section->setName(_metadata.label);
+ _section->setName(_metadata["default"].label);
if (points->size() > 0) {
const point_type pointMin =
@@ -211,14 +256,14 @@
logger.stagePush("Field");
if (fiberDim < 0) {
std::ostringstream msg;
- msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata.label
+ msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata["default"].label
<< "' must be nonnegative.";
throw std::runtime_error(msg.str());
} // if
_section = new section_type(_mesh.comm(), _mesh.debug());
assert(!_section.isNull());
- _section->setName(_metadata.label);
+ _section->setName(_metadata["default"].label);
const int npts = points.size();
if (npts > 0) {
@@ -278,7 +323,7 @@
} // if
if (fiberDim < 0) {
std::ostringstream msg;
- msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata.label
+ msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata["default"].label
<< "' must be nonnegative.";
throw std::runtime_error(msg.str());
} // if
@@ -306,7 +351,7 @@
void
pylith::topology::Field<mesh_type, section_type>::cloneSection(const Field& src)
{ // cloneSection
- std::string origLabel = _metadata.label;
+ std::string origLabel = _metadata["default"].label;
// Clear memory
clear();
@@ -319,7 +364,7 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
- _metadata = src._metadata;
+ _metadata["default"] = const_cast<Field&>(src)._metadata["default"];
label(origLabel.c_str());
if (!_section.isNull()) {
@@ -387,7 +432,7 @@
err = VecGetSize(s_iter->second.vector, &vecGlobalSize);CHECK_PETSC_ERROR(err);
err = VecCreate(_mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata.label.c_str());CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
err = VecSetSizes(sinfo.vector, vecLocalSize, vecGlobalSize);CHECK_PETSC_ERROR(err);
err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
@@ -541,9 +586,9 @@
if (!_section.isNull())
_section->clear();
- _metadata.scale = 1.0;
- _metadata.vectorFieldType = OTHER;
- _metadata.dimsOkay = false;
+ _metadata["default"].scale = 1.0;
+ _metadata["default"].vectorFieldType = OTHER;
+ _metadata["default"].dimsOkay = false;
logger.stagePop();
} // clear
@@ -562,6 +607,12 @@
const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
assert(!sieveMesh.isNull());
sieveMesh->allocate(_section);
+ if (_dm) {
+ PetscSection s;
+ PetscErrorCode err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetUp(s);CHECK_PETSC_ERROR(err);
+ err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+ }
logger.stagePop();
} // allocate
@@ -634,22 +685,22 @@
const int srcSize = field.chartSize();
const int dstSize = chartSize();
if (field.spaceDim() != spaceDim() ||
- field._metadata.vectorFieldType != _metadata.vectorFieldType ||
+ const_cast<Field&>(field)._metadata["default"].vectorFieldType != _metadata["default"].vectorFieldType ||
srcSize != dstSize) {
std::ostringstream msg;
- msg << "Cannot copy values from section '" << field._metadata.label
- << "' to section '" << _metadata.label
+ msg << "Cannot copy values from section '" << const_cast<Field&>(field)._metadata["default"].label
+ << "' to section '" << _metadata["default"].label
<< "'. Sections are incompatible.\n"
<< " Source section:\n"
<< " space dim: " << field.spaceDim() << "\n"
- << " vector field type: " << field._metadata.vectorFieldType << "\n"
- << " scale: " << field._metadata.scale << "\n"
+ << " vector field type: " << const_cast<Field&>(field)._metadata["default"].vectorFieldType << "\n"
+ << " scale: " << const_cast<Field&>(field)._metadata["default"].scale << "\n"
<< " size: " << srcSize << "\n"
<< " Destination section:\n"
<< " space dim: " << spaceDim() << "\n"
- << " vector field type: " << _metadata.vectorFieldType << "\n"
- << " scale: " << _metadata.scale << "\n"
+ << " vector field type: " << _metadata["default"].vectorFieldType << "\n"
+ << " scale: " << _metadata["default"].scale << "\n"
<< " size: " << dstSize;
throw std::runtime_error(msg.str());
} // if
@@ -673,8 +724,8 @@
} // for
} // if
- label(field._metadata.label.c_str()); // Update label
- _metadata.scale = field._metadata.scale;
+ label(const_cast<Field&>(field)._metadata["default"].label.c_str()); // Update label
+ _metadata["default"].scale = const_cast<Field&>(field)._metadata["default"].scale;
} // copy
// ----------------------------------------------------------------------
@@ -690,13 +741,13 @@
std::ostringstream msg;
msg << "Cannot copy values from Sieve section "
- << _metadata.label << "'. Sections are incompatible.\n"
+ << _metadata["default"].label << "'. Sections are incompatible.\n"
<< " Source section:\n"
<< " size: " << srcSize << "\n"
<< " Destination section:\n"
<< " space dim: " << spaceDim() << "\n"
- << " vector field type: " << _metadata.vectorFieldType << "\n"
- << " scale: " << _metadata.scale << "\n"
+ << " vector field type: " << _metadata["default"].vectorFieldType << "\n"
+ << " scale: " << _metadata["default"].scale << "\n"
<< " size: " << dstSize;
throw std::runtime_error(msg.str());
} // if
@@ -729,23 +780,23 @@
const int srcSize = field.chartSize();
const int dstSize = chartSize();
if (field.spaceDim() != spaceDim() ||
- field._metadata.vectorFieldType != _metadata.vectorFieldType ||
- field._metadata.scale != _metadata.scale ||
+ const_cast<Field&>(field)._metadata["default"].vectorFieldType != _metadata["default"].vectorFieldType ||
+ const_cast<Field&>(field)._metadata["default"].scale != _metadata["default"].scale ||
srcSize != dstSize) {
std::ostringstream msg;
- msg << "Cannot add values from section '" << field._metadata.label
- << "' to section '" << _metadata.label
+ msg << "Cannot add values from section '" << const_cast<Field&>(field)._metadata["default"].label
+ << "' to section '" << _metadata["default"].label
<< "'. Sections are incompatible.\n"
<< " Source section:\n"
<< " space dim: " << field.spaceDim() << "\n"
- << " vector field type: " << field._metadata.vectorFieldType << "\n"
- << " scale: " << field._metadata.scale << "\n"
+ << " vector field type: " << const_cast<Field&>(field)._metadata["default"].vectorFieldType << "\n"
+ << " scale: " << const_cast<Field&>(field)._metadata["default"].scale << "\n"
<< " size: " << srcSize << "\n"
<< " Destination section:\n"
<< " space dim: " << spaceDim() << "\n"
- << " vector field type: " << _metadata.vectorFieldType << "\n"
- << " scale: " << _metadata.scale << "\n"
+ << " vector field type: " << _metadata["default"].vectorFieldType << "\n"
+ << " scale: " << _metadata["default"].scale << "\n"
<< " size: " << dstSize;
throw std::runtime_error(msg.str());
} // if
@@ -784,9 +835,9 @@
void
pylith::topology::Field<mesh_type, section_type>::dimensionalize(void) const
{ // dimensionalize
- if (!_metadata.dimsOkay) {
+ if (!const_cast<Field*>(this)->_metadata["default"].dimsOkay) {
std::ostringstream msg;
- msg << "Cannot dimensionalize field '" << _metadata.label
+ msg << "Cannot dimensionalize field '" << const_cast<Field*>(this)->_metadata["default"].label
<< "' because the flag "
<< "has been set to keep field nondimensional.";
throw std::runtime_error(msg.str());
@@ -810,7 +861,7 @@
assert(fiberDim == _section->getFiberDimension(*c_iter));
_section->restrictPoint(*c_iter, &values[0], values.size());
- normalizer.dimensionalize(&values[0], values.size(), _metadata.scale);
+ normalizer.dimensionalize(&values[0], values.size(), const_cast<Field*>(this)->_metadata["default"].scale);
_section->updatePointAll(*c_iter, &values[0]);
} // if
} // if
@@ -823,7 +874,7 @@
pylith::topology::Field<mesh_type, section_type>::view(const char* label) const
{ // view
std::string vecFieldString;
- switch(_metadata.vectorFieldType)
+ switch(const_cast<Field*>(this)->_metadata["default"].vectorFieldType)
{ // switch
case SCALAR:
vecFieldString = "scalar";
@@ -850,16 +901,16 @@
vecFieldString = "multiple other values";
break;
default :
- std::cerr << "Unknown vector field value '" << _metadata.vectorFieldType
+ std::cerr << "Unknown vector field value '" << const_cast<Field*>(this)->_metadata["default"].vectorFieldType
<< "'." << std::endl;
assert(0);
throw std::logic_error("Bad vector field type in Field.");
} // switch
- std::cout << "Viewing field '" << _metadata.label << "' "<< label << ".\n"
+ std::cout << "Viewing field '" << const_cast<Field*>(this)->_metadata["default"].label << "' "<< label << ".\n"
<< " vector field type: " << vecFieldString << "\n"
- << " scale: " << _metadata.scale << "\n"
- << " dimensionalize flag: " << _metadata.dimsOkay << std::endl;
+ << " scale: " << const_cast<Field*>(this)->_metadata["default"].scale << "\n"
+ << " dimensionalize flag: " << const_cast<Field*>(this)->_metadata["default"].dimsOkay << std::endl;
if (!_section.isNull())
_section->view(label);
} // view
@@ -922,7 +973,7 @@
err = VecCreate(_mesh.comm(), &sinfo.vector);
CHECK_PETSC_ERROR(err);
err = PetscObjectSetName((PetscObject)sinfo.vector,
- _metadata.label.c_str());CHECK_PETSC_ERROR(err);
+ _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
err = VecSetSizes(sinfo.vector,
order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
@@ -1000,7 +1051,7 @@
// Create vector
err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
err = PetscObjectSetName((PetscObject)sinfo.vector,
- _metadata.label.c_str());CHECK_PETSC_ERROR(err);
+ _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
err = VecSetSizes(sinfo.vector,
order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
@@ -1079,7 +1130,7 @@
// Create vector
err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata.label.c_str());CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
err = VecSetSizes(sinfo.vector, order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
@@ -1155,7 +1206,7 @@
// Create vector
err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata.label.c_str());CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
err = VecSetSizes(sinfo.vector,order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
@@ -1378,5 +1429,63 @@
return s_iter->second;
} // _getScatter
+// Experimental
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::Field<mesh_type, section_type>::addField(const std::string& name, int numComponents)
+{
+ // Keep track of name/components until setup
+ _tmpFields[name] = numComponents;
+}
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::Field<mesh_type, section_type>::setupFields()
+{
+ assert(_dm);
+ // Keep track of name/components until setup
+ PetscSection section;
+ PetscInt f = 0;
+ PetscErrorCode err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
+ assert(section);
+ err = PetscSectionSetNumFields(section, _tmpFields.size());CHECK_PETSC_ERROR(err);
+ for(std::map<std::string, int>::const_iterator f_iter = _tmpFields.begin(); f_iter != _tmpFields.end(); ++f_iter, ++f) {
+ err = PetscSectionSetFieldName(section, f, f_iter->first.c_str());CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetFieldComponents(section, f, f_iter->second);CHECK_PETSC_ERROR(err);
+ }
+ _tmpFields.clear();
+}
+
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::Field<mesh_type, section_type>::updateDof(const std::string& name, const DomainEnum domain, int fiberDim)
+{
+ PetscSection section;
+ PetscInt pStart, pEnd, f = 0;
+ PetscErrorCode err;
+
+ assert(_dm);
+ switch(domain) {
+ case VERTICES_FIELD:
+ err = DMComplexGetDepthStratum(_dm, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
+ case CELLS_FIELD:
+ err = DMComplexGetHeightStratum(_dm, 0, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
+ default:
+ std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
+ throw std::logic_error("Bad domain enum in Field.");
+ }
+ err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
+ assert(section);
+ for(map_type::const_iterator f_iter = _metadata.begin(); f_iter != _metadata.end(); ++f_iter, ++f) {
+ if (f_iter->first == name) break;
+ }
+ assert(f < _metadata.size());
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ err = PetscSectionAddDof(section, p, fiberDim);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetFieldDof(section, p, f, fiberDim);CHECK_PETSC_ERROR(err);
+ }
+}
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-08-25 15:24:08 UTC (rev 20631)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2012-08-25 19:54:25 UTC (rev 20632)
@@ -85,6 +85,9 @@
const ALE::Obj<section_type>& section,
const Metadata& metadata);
+
+ Field(const Field& src, int numFields, const int fields[]);
+
/// Destructor.
~Field(void);
@@ -127,6 +130,13 @@
*/
void vectorFieldType(const VectorFieldEnum value);
+ /** Set vector field type
+ *
+ * @param name Field name
+ * @param value Type of vector field.
+ */
+ void vectorFieldType(const std::string& name, const VectorFieldEnum value);
+
/** Get vector field type
*
* @returns Type of vector field.
@@ -139,6 +149,13 @@
*/
void scale(const PylithScalar value);
+ /** Set scale for dimensionalizing field.
+ *
+ * @param name Field name
+ * @param value Scale associated with field.
+ */
+ void scale(const std::string& name, const PylithScalar value);
+
/** Get scale for dimensionalizing field.
*
* @returns Scale associated with field.
@@ -236,6 +253,12 @@
*/
void cloneSection(const Field& src);
+ void addField(const std::string& name, int numComponents);
+
+ void setupFields();
+
+ void updateDof(const std::string& name, const DomainEnum domain, const int fiberDim);
+
/// Clear variables associated with section.
void clear(void);
@@ -422,17 +445,24 @@
*/
const ScatterInfo& _getScatter(const char* context) const;
+// PROTECTED TYPEDEFS ///////////////////////////////////////////////////
+protected :
+
+ typedef std::map< std::string, Metadata > map_type;
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
- Metadata _metadata;
+ map_type _metadata;
+ /* Old construction */
const mesh_type& _mesh; ///< Mesh associated with section.
ALE::Obj<section_type> _section; ///< Real section with data.
- PetscSection _newSection;
- Vec _newLocalVec;
scatter_map_type _scatters; ///< Collection of scatters.
+ /* New construction */
+ DM _dm; /* Holds the PetscSection */
+ Vec _globalVec;
+ std::map<std::string, int> _tmpFields;
-
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-08-25 15:24:08 UTC (rev 20631)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2012-08-25 19:54:25 UTC (rev 20632)
@@ -44,7 +44,7 @@
inline
const char*
pylith::topology::Field<mesh_type, section_type>::label(void) const {
- return _metadata.label.c_str();
+ return const_cast<Field*>(this)->_metadata["default"].label.c_str();
}
// Set vector field type
@@ -52,15 +52,23 @@
inline
void
pylith::topology::Field<mesh_type, section_type>::vectorFieldType(const VectorFieldEnum value) {
- _metadata.vectorFieldType = value;
+ _metadata["default"].vectorFieldType = value;
}
+// Set vector field type
+template<typename mesh_type, typename section_type>
+inline
+void
+pylith::topology::Field<mesh_type, section_type>::vectorFieldType(const std::string& name, const VectorFieldEnum value) {
+ const_cast<Field *>(this)->_metadata[name].vectorFieldType = value;
+}
+
// Get vector field type
template<typename mesh_type, typename section_type>
inline
typename pylith::topology::Field<mesh_type, section_type>::VectorFieldEnum
pylith::topology::Field<mesh_type, section_type>::vectorFieldType(void) const {
- return _metadata.vectorFieldType;
+ return const_cast<Field *>(this)->_metadata["default"].vectorFieldType;
}
// Set scale for dimensionalizing field.
@@ -68,15 +76,23 @@
inline
void
pylith::topology::Field<mesh_type, section_type>::scale(const PylithScalar value) {
- _metadata.scale = value;
+ _metadata["default"].scale = value;
}
+// Set scale for dimensionalizing field.
+template<typename mesh_type, typename section_type>
+inline
+void
+pylith::topology::Field<mesh_type, section_type>::scale(const std::string& name, const PylithScalar value) {
+ _metadata[name].scale = value;
+}
+
// Get scale for dimensionalizing field.
template<typename mesh_type, typename section_type>
inline
PylithScalar
pylith::topology::Field<mesh_type, section_type>::scale(void) const {
- return _metadata.scale;
+ return const_cast<Field *>(this)->_metadata["default"].scale;
}
// Set flag indicating whether it is okay to dimensionalize field.
@@ -84,7 +100,7 @@
inline
void
pylith::topology::Field<mesh_type, section_type>::addDimensionOkay(const bool value) {
- _metadata.dimsOkay = value;
+ _metadata["default"].dimsOkay = value;
}
// Set flag indicating whether it is okay to dimensionalize field.
@@ -92,7 +108,7 @@
inline
bool
pylith::topology::Field<mesh_type, section_type>::addDimensionOkay(void) const {
- return _metadata.dimsOkay;
+ return const_cast<Field*>(this)->_metadata["default"].dimsOkay;
}
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.hh 2012-08-25 15:24:08 UTC (rev 20631)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.hh 2012-08-25 19:54:25 UTC (rev 20632)
@@ -92,6 +92,7 @@
* @param fields Set of field numbers for this vector
*/
void add(const char* name,
+ const char* oldName,
const int numFields,
const int fields[]);
@@ -119,12 +120,6 @@
*/
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.
@@ -149,13 +144,11 @@
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-25 15:24:08 UTC (rev 20631)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Fields.icc 2012-08-25 19:54:25 UTC (rev 20632)
@@ -51,11 +51,6 @@
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
// ----------------------------------------------------------------------
@@ -64,9 +59,7 @@
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
@@ -86,10 +79,6 @@
_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
// ----------------------------------------------------------------------
@@ -112,36 +101,6 @@
_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
// ----------------------------------------------------------------------
@@ -150,6 +109,7 @@
void
pylith::topology::Fields<field_type>::add(
const char* name,
+ const char* oldName,
const int numFields,
const int fields[])
{ // add
@@ -160,15 +120,7 @@
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;
- }
+ _fields[name] = new field_type(_fields[oldName], numFields, fields);
} // add
// ----------------------------------------------------------------------
@@ -186,7 +138,6 @@
} // if
delete iter->second; iter->second = 0;
_fields.erase(name);
- _dmFields.erase(name);
} // del
// ----------------------------------------------------------------------
@@ -230,22 +181,6 @@
} // 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.
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2012-08-25 15:24:08 UTC (rev 20631)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2012-08-25 19:54:25 UTC (rev 20632)
@@ -532,22 +532,41 @@
self.fields.solutionName("dispIncr(t->t+dt)")
lengthScale = normalizer.lengthScale()
- solution = self.fields.get("dispIncr(t->t+dt)")
- solution.newSection(solution.VERTICES_FIELD, dimension)
- solution.vectorFieldType(solution.VECTOR)
- solution.scale(lengthScale.value)
- if self.splitFields():
- solution.splitDefault()
+ if 1:
+ solution = self.fields.get("dispIncr(t->t+dt)")
+ solution.newSection(solution.VERTICES_FIELD, dimension)
+ solution.vectorFieldType(solution.VECTOR)
+ solution.scale(lengthScale.value)
+ if self.splitFields():
+ solution.splitDefault()
+ for integrator in self.integratorsMesh + self.integratorsSubMesh:
+ integrator.splitField(solution)
+ for constraint in self.constraints:
+ constraint.setConstraintSizes(solution)
+ solution.allocate()
+ for constraint in self.constraints:
+ constraint.setConstraints(solution)
for integrator in self.integratorsMesh + self.integratorsSubMesh:
- integrator.splitField(solution)
- for constraint in self.constraints:
- constraint.setConstraintSizes(solution)
- solution.allocate()
- for constraint in self.constraints:
- constraint.setConstraints(solution)
- for integrator in self.integratorsMesh + self.integratorsSubMesh:
- integrator.checkConstraints(solution)
+ integrator.checkConstraints(solution)
+ else:
+ solution.addField("displacement", dimension)
+ solution.addField("pressure", 1)
+ solution.addField("temperature", 1)
+ solution.setupFields()
+ solution.updateDof("displacement", solution.VERTICES_FIELD, dimension)
+ solution.updateDof("pressure", solution.CELLS_FIELD, 1)
+ solution.updateDof("temperature", solution.VERTICES_FIELD, 1)
+ solution.vectorFieldType("displacement", solution.VECTOR)
+ solution.scale("displacement", lengthScale.value)
+ for constraint in self.constraints:
+ constraint.setConstraintSizes(solution)
+ solution.allocate()
+ for constraint in self.constraints:
+ constraint.setConstraints(solution)
+ for integrator in self.integratorsMesh + self.integratorsSubMesh:
+ integrator.checkConstraints(solution)
+
memoryLogger.stagePop()
# This also creates a global order.
More information about the CIG-COMMITS
mailing list