[cig-commits] r13970 - short/3D/PyLith/branches/pylith-swig/libsrc/topology
brad at geodynamics.org
brad at geodynamics.org
Tue Jan 27 17:35:54 PST 2009
Author: brad
Date: 2009-01-27 17:35:54 -0800 (Tue, 27 Jan 2009)
New Revision: 13970
Added:
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.icc
Log:
Added missing Field files.
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc 2009-01-28 01:35:54 UTC (rev 13970)
@@ -0,0 +1,367 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "Field.hh" // implementation of class methods
+
+#include "Mesh.hh" // HASA Mesh
+#include "pylith/utils/array.hh" // USES double_array
+
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::topology::Field::Field(const Mesh& mesh) :
+ _mesh(mesh)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::topology::Field::~Field(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Get spatial dimension of domain.
+int
+pylith::topology::Field::spaceDim(void) const
+{ // spaceDim
+ const spatialdata::geocoords::CoordSys* cs = _mesh.coordsys();
+ return (0 != cs) ? cs->spaceDim() : 0;
+} // spaceDim
+
+// ----------------------------------------------------------------------
+// Create seive section.
+void
+pylith::topology::Field::newSection(void)
+{ // newSection
+ _section = new MeshRealSection(_mesh.comm(), _mesh.debug());
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion.
+void
+pylith::topology::Field::newSection(
+ const ALE::Obj<SieveMesh::label_sequence>& points,
+ const int fiberDim)
+{ // newSection
+ if (fiberDim < 0) {
+ std::ostringstream msg;
+ msg
+ << "Fiber dimension (" << fiberDim << ") for field '" << _name
+ << "' must be nonnegative.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ _section = new MeshRealSection(_mesh.comm(), _mesh.debug());
+
+ const SieveMesh::point_type pointMin =
+ *std::min_element(points->begin(), points->end());
+ const SieveMesh::point_type pointMax =
+ *std::max_element(points->begin(), points->end());
+ _section->setChart(MeshRealSection::chart_type(pointMin, pointMax+1));
+ _section->setFiberDimension(points, fiberDim);
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion.
+void
+pylith::topology::Field::newSection(const DomainEnum domain,
+ const int fiberDim)
+{ // newSection
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ ALE::Obj<SieveMesh::label_sequence> points;
+ if (VERTICES_FIELD == domain)
+ points = sieveMesh->depthStratum(0);
+ else if (CELLS_FIELD == domain)
+ points = sieveMesh->heightStratum(1);
+ else {
+ std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
+ assert(0);
+ } // else
+
+ newSection(points, fiberDim);
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create section given chart.
+void
+pylith::topology::Field::newSection(const MeshRealSection::chart_type& chart,
+ const int fiberDim)
+{ // newSection
+ if (_section.isNull())
+ newSection();
+
+ _section->setChart(chart);
+
+ const MeshRealSection::chart_type::const_iterator chartEnd = chart.end();
+ for (MeshRealSection::chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter)
+ _section->setFiberDimension(*c_iter, fiberDim);
+ allocate();
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create section with same layout as another section.
+void
+pylith::topology::Field::newSection(const Field& src)
+{ // newSection
+ _vecFieldType = src._vecFieldType;
+
+ const ALE::Obj<MeshRealSection>& srcSection = src.section();
+ if (!srcSection.isNull() && _section.isNull())
+ newSection();
+
+ if (!_section.isNull()) {
+ _section->setAtlas(srcSection->getAtlas());
+ _section->allocateStorage();
+ _section->setBC(srcSection->getBC());
+ } // if
+} // newSection
+
+// ----------------------------------------------------------------------
+// Clear variables associated with section.
+void
+pylith::topology::Field::clear(void)
+{ // clear
+ if (!_section.isNull())
+ _section->clear();
+
+ _scale = 1.0;
+ _vecFieldType = OTHER;
+ _dimensionsOkay = false;
+} // clear
+
+// ----------------------------------------------------------------------
+// Allocate Sieve section.
+void
+pylith::topology::Field::allocate(void)
+{ // allocate
+ assert(!_section.isNull());
+
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ sieveMesh->allocate(_section);
+} // allocate
+
+// ----------------------------------------------------------------------
+// Zero section values.
+void
+pylith::topology::Field::zero(void)
+{ // zero
+ if (!_section.isNull())
+ _section->zero();
+} // zero
+
+// ----------------------------------------------------------------------
+// Complete section by assembling across processors.
+void
+pylith::topology::Field::complete(void)
+{ // complete
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ if (!_section.isNull())
+ ALE::Completion::completeSectionAdd(sieveMesh->getSendOverlap(),
+ sieveMesh->getRecvOverlap(),
+ _section, _section);
+} // complete
+
+// ----------------------------------------------------------------------
+// Copy field values and metadata.
+void
+pylith::topology::Field::copy(const Field& field)
+{ // copy
+ // Check compatibility of sections
+ const int srcSize = (!field._section.isNull()) ? field._section->size() : 0;
+ const int dstSize = (!_section.isNull()) ? _section->size() : 0;
+ if (field.spaceDim() != spaceDim() ||
+ field._vecFieldType != _vecFieldType ||
+ field._scale != _scale ||
+ srcSize != dstSize) {
+ std::ostringstream msg;
+
+ msg << "Cannot copy values from section '" << field._name
+ << "' to section '" << _name << "'. Sections are incompatible.\n"
+ << " Source section:\n"
+ << " space dim: " << field.spaceDim() << "\n"
+ << " vector field type: " << field._vecFieldType << "\n"
+ << " scale: " << field._scale << "\n"
+ << " size: " << srcSize
+ << " Destination section:\n"
+ << " space dim: " << spaceDim() << "\n"
+ << " vector field type: " << _vecFieldType << "\n"
+ << " scale: " << _scale << "\n"
+ << " size: " << dstSize;
+ throw std::runtime_error(msg.str());
+ } // if
+ assert( (_section.isNull() && field._section.isNull()) ||
+ (!_section.isNull() && !field._section.isNull()) );
+
+ if (!_section.isNull()) {
+ // Copy values from field
+ const MeshRealSection::chart_type& chart = _section->getChart();
+ const MeshRealSection::chart_type::const_iterator chartEnd = chart.end();
+
+ for (MeshRealSection::chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter) {
+ assert(field._section->getFiberDimension(*c_iter) ==
+ _section->getFiberDimension(*c_iter));
+ _section->updatePoint(*c_iter, field._section->restrictPoint(*c_iter));
+ } // for
+ } // if
+} // copy
+
+// ----------------------------------------------------------------------
+// Add two fields, storing the result in one of the fields.
+void
+pylith::topology::Field::operator+=(const Field& field)
+{ // operator+=
+ // Check compatibility of sections
+ const int srcSize = (!field._section.isNull()) ? field._section->size() : 0;
+ const int dstSize = (!_section.isNull()) ? _section->size() : 0;
+ if (field.spaceDim() != spaceDim() ||
+ field._vecFieldType != _vecFieldType ||
+ field._scale != _scale ||
+ srcSize != dstSize) {
+ std::ostringstream msg;
+
+ msg << "Cannot add values from section '" << field._name
+ << "' to section '" << _name << "'. Sections are incompatible.\n"
+ << " Source section:\n"
+ << " space dim: " << field.spaceDim() << "\n"
+ << " vector field type: " << field._vecFieldType << "\n"
+ << " scale: " << field._scale << "\n"
+ << " size: " << srcSize
+ << " Destination section:\n"
+ << " space dim: " << spaceDim() << "\n"
+ << " vector field type: " << _vecFieldType << "\n"
+ << " scale: " << _scale << "\n"
+ << " size: " << dstSize;
+ throw std::runtime_error(msg.str());
+ } // if
+ assert( (_section.isNull() && field._section.isNull()) ||
+ (!_section.isNull() && !field._section.isNull()) );
+
+ if (!_section.isNull()) {
+ // Add values from field
+ const MeshRealSection::chart_type& chart = _section->getChart();
+ const MeshRealSection::chart_type::const_iterator chartEnd = chart.end();
+
+ // Assume fiber dimension is uniform
+ const int fiberDim = _section->getFiberDimension(*chart.begin());
+ double_array values(fiberDim);
+
+ for (MeshRealSection::chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter) {
+ assert(fiberDim == field._section->getFiberDimension(*c_iter));
+ assert(fiberDim == _section->getFiberDimension(*c_iter));
+ field._section->restrictPoint(*c_iter, &values[0], values.size());
+ _section->updateAddPoint(*c_iter, &values[0]);
+ } // for
+ } // if
+} // operator+=
+
+// ----------------------------------------------------------------------
+// Dimensionalize field.
+void
+pylith::topology::Field::dimensionalize(void)
+{ // dimensionalize
+ if (!_dimensionsOkay) {
+ std::ostringstream msg;
+ msg << "Cannot dimensionalize field '" << _name << "' because the flag "
+ << "has been set to keep field nondimensional.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ if (!_section.isNull()) {
+ const MeshRealSection::chart_type& chart = _section->getChart();
+ const MeshRealSection::chart_type::const_iterator chartEnd = chart.end();
+
+ // Assume fiber dimension is uniform
+ const int fiberDim = _section->getFiberDimension(*chart.begin());
+ double_array values(fiberDim);
+
+ spatialdata::units::Nondimensional normalizer;
+
+ for (MeshRealSection::chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter) {
+ assert(fiberDim == _section->getFiberDimension(*c_iter));
+
+ _section->restrictPoint(*c_iter, &values[0], values.size());
+ normalizer.dimensionalize(&values[0], values.size(), _scale);
+ _section->updatePoint(*c_iter, &values[0]);
+ } // for
+ } // if
+} // dimensionalize
+
+// ----------------------------------------------------------------------
+// Print field to standard out.
+void
+pylith::topology::Field::view(const char* label)
+{ // view
+ std::string vecFieldString;
+ switch(_vecFieldType)
+ { // switch
+ case SCALAR:
+ vecFieldString = "scalar";
+ break;
+ case VECTOR:
+ vecFieldString = "vector";
+ break;
+ case TENSOR:
+ vecFieldString = "tensor";
+ break;
+ case OTHER:
+ vecFieldString = "other";
+ break;
+ case MULTI_SCALAR:
+ vecFieldString = "multiple scalars";
+ break;
+ case MULTI_VECTOR:
+ vecFieldString = "multiple vectors";
+ break;
+ case MULTI_TENSOR:
+ vecFieldString = "multiple tensors";
+ break;
+ case MULTI_OTHER:
+ vecFieldString = "multiple other values";
+ break;
+ default :
+ std::cerr << "Unknown vector field value '" << _vecFieldType
+ << "'." << std::endl;
+ assert(0);
+ } // switch
+
+ std::cout << "Viewing field '" << _name << "' "<< label << ".\n"
+ << " vector field type: " << vecFieldString << "\n"
+ << " scale: " << _scale << "\n"
+ << " dimensionalize flag: " << _dimensionsOkay << std::endl;
+ if (!_section.isNull())
+ _section->view(label);
+} // view
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh 2009-01-28 01:35:54 UTC (rev 13970)
@@ -0,0 +1,161 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/topology/Field.hh
+ *
+ * @brief Vector field over the vertices or cells of a finite-element
+ * mesh.
+ *
+ * Extends Sieve real general section by adding metadata.
+ */
+
+#include "FieldBase.hh" // ISA FieldBase
+
+#if !defined(pylith_topology_field_hh)
+#define pylith_topology_field_hh
+
+// Include directives ---------------------------------------------------
+#define NEWPYLITHMESH 1
+#include "pylith/utils/sievetypes.hh" // HASA PETSc real_section_type
+
+// Forward declarations -------------------------------------------------
+namespace pylith {
+ namespace topology {
+ class Field;
+ class TestField;
+
+ class Mesh; // HASA Mesh
+ } // topology
+} // pylith
+
+// Field ----------------------------------------------------------------
+class pylith::topology::Field : public FieldBase
+{ // Field
+ friend class TestField; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Default constructor.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ Field(const Mesh& mesh);
+
+ /// Destructor.
+ ~Field(void);
+
+ /** Get Sieve section.
+ *
+ * @returns Sieve section.
+ */
+ const ALE::Obj<MeshRealSection>& section(void) const;
+
+ /** Get spatial dimension of domain.
+ *
+ * @returns Spatial dimension of domain.
+ */
+ int spaceDim(void) const;
+
+ /// Create sieve section.
+ void newSection(void);
+
+ /** Create sieve section and set chart and fiber dimesion.
+ *
+ * @param points Points over which to define section.
+ * @param dim Fiber dimension for section.
+ */
+ void newSection(const ALE::Obj<SieveMesh::label_sequence>& points,
+ const int fiberDim);
+
+ /** Create sieve section and set chart and fiber dimesion.
+ *
+ * @param domain Type of points over which to define section.
+ * @param dim Fiber dimension for section.
+ */
+ void newSection(const DomainEnum domain,
+ const int fiberDim);
+
+ /** Create section given chart. This allows a chart to be reused
+ * across multiple fields, reducing memory usage.
+ *
+ * @param chart Chart defining points over which section is defined.
+ * @param fiberDim Fiber dimension.
+ */
+ void newSection(const MeshRealSection::chart_type& chart,
+ const int fiberDim);
+
+ /** Create section with same layout (fiber dimension and
+ * constraints) as another section. This allows the layout data
+ * structures to be reused across multiple fields, reducing memory
+ * usage.
+ *
+ * @param sec Section defining layout.
+ */
+ void newSection(const Field& src);
+
+ /// Clear variables associated with section.
+ void clear(void);
+
+ /// Allocate field.
+ void allocate(void);
+
+ /// Zero section values.
+ void zero(void);
+
+ /// Complete section by assembling across processors.
+ void complete(void);
+
+ /** Copy field values and metadata.
+ *
+ * @param field Field to copy.
+ */
+ void copy(const Field& field);
+
+ /** Add two fields, storing the result in one of the fields.
+ *
+ * @param field Field to add.
+ */
+ void operator+=(const Field& field);
+
+ /** Dimensionalize field. Throws runtime_error if field is not
+ * allowed to be dimensionalized.
+ */
+ void dimensionalize(void);
+
+ /** Print field to standard out.
+ *
+ * @param label Label for output.
+ */
+ void view(const char* label);
+
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ const Mesh& _mesh; ///< Mesh associated with section
+ ALE::Obj<MeshRealSection> _section; ///< Real section with data
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ Field(const Field&); ///< Not implemented
+ const Field& operator=(const Field&); ///< Not implemented
+
+}; // Field
+
+#include "Field.icc"
+
+#endif // pylith_topology_field_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.icc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.icc 2009-01-28 01:35:54 UTC (rev 13970)
@@ -0,0 +1,27 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_topology_field_hh)
+#error "Field.icc must be included only from Field.hh"
+#else
+
+// Get Sieve section.
+inline
+const ALE::Obj<pylith::MeshRealSection>&
+pylith::topology::Field::section(void) const {
+ return _section;
+}
+
+#endif
+
+
+// End of file
More information about the CIG-COMMITS
mailing list