[cig-commits] r13960 - in short/3D/PyLith/branches/pylith-swig/libsrc: . bc topology utils
brad at geodynamics.org
brad at geodynamics.org
Mon Jan 26 21:22:03 PST 2009
Author: brad
Date: 2009-01-26 21:22:03 -0800 (Mon, 26 Jan 2009)
New Revision: 13960
Added:
short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.hh
short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.icc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.hh
short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.icc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh
short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.icc
Removed:
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
short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldUniform.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldUniform.hh
Modified:
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBC.cc
short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBC.hh
short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.cc
short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.hh
short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.icc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.icc
short/3D/PyLith/branches/pylith-swig/libsrc/utils/sievetypes.hh
Log:
Added SubMesh and FieldSubMesh. Cleaned up Field.
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-01-27 05:22:03 UTC (rev 13960)
@@ -93,14 +93,16 @@
meshio/VertexFilterVecNorm.cc \
meshio/UCDFaultFile.cc \
topology/Distributor.cc \
+ topology/FieldBase.cc \
topology/Field.cc \
- topology/FieldUniform.cc \
+ topology/FieldSubMesh.cc \
topology/FieldsManager.cc \
topology/FieldOps.cc \
topology/Mesh.cc \
topology/MeshOps.cc \
topology/MeshRefiner.cc \
topology/RefineUniform.cc \
+ topology/SubMesh.cc \
utils/EventLogger.cc
# bc/Neumann.cc
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBC.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBC.cc 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBC.cc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -79,7 +79,7 @@
if (0 == numFixedDOF)
return;
- const ALE::Obj<SieveRealSection>& section = field.section();
+ const ALE::Obj<MeshRealSection>& section = field.section();
assert(!section.isNull());
const int numPoints = _points.size();
@@ -113,7 +113,7 @@
if (0 == numFixedDOF)
return;
- const ALE::Obj<SieveRealSection>& section = field.section();
+ const ALE::Obj<MeshRealSection>& section = field.section();
assert(!section.isNull());
const int numPoints = _points.size();
@@ -176,7 +176,7 @@
if (0 == numFixedDOF)
return;
- const ALE::Obj<SieveRealSection>& section = field.section();
+ const ALE::Obj<MeshRealSection>& section = field.section();
assert(!section.isNull());
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
assert(!sieveMesh.isNull());
@@ -271,7 +271,7 @@
const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
assert(!sieveMesh.isNull());
- const ALE::Obj<SieveRealSection>& coordinates =
+ const ALE::Obj<MeshRealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBC.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBC.hh 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBC.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -38,7 +38,7 @@
// DirichletBC ------------------------------------------------------
class pylith::bc::DirichletBC : public BoundaryCondition,
- public feassemble::Constraint
+ public feassemble::Constraint
{ // class DirichletBC
friend class TestDirichletBC; // unit testing
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.cc 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.cc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -14,13 +14,13 @@
#include "DirichletBoundary.hh" // implementation of object methods
-#include "pylith/topology/FieldUniform.hh" // USES FieldUniform
#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/FieldSubMesh.hh" // USES FieldSubMesh
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
-#include <Selection.hh> // USES submesh algorithms
-
#include <strings.h> // USES strcasecmp()
#include <cassert> // USES assert()
#include <stdexcept> // USES std::runtime_error
@@ -29,6 +29,7 @@
// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::DirichletBoundary::DirichletBoundary(void) :
+ _boundaryMesh(0),
_tmpField(0)
{ // constructor
} // constructor
@@ -37,6 +38,7 @@
// Destructor.
pylith::bc::DirichletBoundary::~DirichletBoundary(void)
{ // destructor
+ delete _boundaryMesh; _boundaryMesh = 0;
delete _tmpField; _tmpField = 0;
} // destructor
@@ -50,7 +52,7 @@
if (0 == numFixedDOF)
return;
- _createBoundaryMesh(mesh);
+ _boundaryMesh = new topology::SubMesh(mesh, _label.c_str());
_getPoints(mesh);
_setupQueryDatabases();
_queryDatabases(mesh);
@@ -58,45 +60,52 @@
// ----------------------------------------------------------------------
// Get vertex field of BC initial or rate of change of values.
-const pylith::topology::Field&
+const pylith::topology::FieldSubMesh&
pylith::bc::DirichletBoundary::vertexField(const char* name,
const topology::Mesh& mesh,
const topology::SolutionFields& fields)
{ // getVertexField
assert(0 != name);
- assert(!_boundaryMesh.isNull());
+ assert(0 != _boundaryMesh);
assert(0 != _normalizer);
+ const ALE::Obj<SieveSubMesh>& sieveMesh = _boundaryMesh->sieveMesh();
+ assert(!sieveMesh.isNull());
+
const ALE::Obj<SieveMesh::label_sequence>& vertices =
- _boundaryMesh->depthStratum(0);
+ sieveMesh->depthStratum(0);
assert(!vertices.isNull());
const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
assert(0 != cs);
- const int spaceDim = cs->spaceDim();
- double_array values(spaceDim);
+ const int fiberDim = cs->spaceDim();
+ double_array values(fiberDim);
const int numPoints = _points.size();
const int numFixedDOF = _fixedDOF.size();
if (0 == _tmpField) {
- _tmpField = new topology::FieldUniform(_boundaryMesh, spaceDim);
+ _tmpField = new topology::FieldSubMesh(*_boundaryMesh);
assert(0 != _tmpField);
- _tmpField->createSection(vertices);
+ _tmpField->newSection(vertices, fiberDim);
+ _tmpField->allocate();
} // if
+ // ERROR: NEED TO TRANSLATE LABELS FROM MESH INTO SUBMESH
+ assert(0);
+
if (0 == strcasecmp(name, "initial")) {
_tmpField->name("displacement");
_tmpField->vectorFieldType(topology::Field::VECTOR);
_tmpField->scale(_normalizer->lengthScale());
_tmpField->addDimensionOkay(true);
_tmpField->zero();
- const ALE::Obj<SieveRealSection>& section = _tmpField->section();
+ const ALE::Obj<SubMeshRealSection>& section = _tmpField->section();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
- const SieveMesh::point_type point = _points[iPoint];
- assert(spaceDim == section->getFiberDimension(point));
+ const SieveSubMesh::point_type point = _points[iPoint];
+ assert(fiberDim == section->getFiberDimension(point));
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
values[_fixedDOF[iDOF]] = _valuesInitial[iPoint*numFixedDOF+iDOF];
section->updatePointAll(_points[iPoint], &values[0]);
@@ -107,11 +116,11 @@
_tmpField->scale(_normalizer->lengthScale());
_tmpField->addDimensionOkay(true);
_tmpField->zero();
- const ALE::Obj<SieveRealSection>& section = _tmpField->section();
+ const ALE::Obj<SubMeshRealSection>& section = _tmpField->section();
for (int iPoint=0; iPoint < numPoints; ++iPoint) {
const SieveMesh::point_type point = _points[iPoint];
- assert(spaceDim == section->getFiberDimension(point));
+ assert(fiberDim == section->getFiberDimension(point));
for (int iDOF=0; iDOF < numFixedDOF; ++iDOF)
values[_fixedDOF[iDOF]] = _valuesRate[iPoint*numFixedDOF+iDOF];
section->updatePointAll(_points[iPoint], &values[0]);
@@ -127,45 +136,5 @@
return *_tmpField;
} // getVertexField
-// ----------------------------------------------------------------------
-// Extract submesh associated with boundary.
-void
-pylith::bc::DirichletBoundary::_createBoundaryMesh(const topology::Mesh& mesh)
-{ // _createBoundaryMesh
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::int_section_type>& groupField =
- sieveMesh->getIntSection(_label);
- if (groupField.isNull()) {
- std::ostringstream msg;
- msg << "Could not find group of points '" << _label << "' in mesh.";
- throw std::runtime_error(msg.str());
- } // if
- _boundaryMesh =
- ALE::Selection<SieveMesh>::submeshV<SieveSubMesh>(sieveMesh, groupField);
- if (_boundaryMesh.isNull()) {
- std::ostringstream msg;
- msg << "Could not construct boundary mesh for Dirichlet boundary "
- << "condition '" << _label << "'.";
- throw std::runtime_error(msg.str());
- } // if
- _boundaryMesh->setRealSection("coordinates",
- sieveMesh->getRealSection("coordinates"));
- // Create the parallel overlap
- ALE::Obj<SieveSubMesh::send_overlap_type> sendParallelMeshOverlap =
- _boundaryMesh->getSendOverlap();
- ALE::Obj<SieveSubMesh::recv_overlap_type> recvParallelMeshOverlap =
- _boundaryMesh->getRecvOverlap();
- SieveMesh::renumbering_type& renumbering = sieveMesh->getRenumbering();
- // Can I figure this out in a nicer way?
- ALE::SetFromMap<std::map<SieveMesh::point_type,SieveMesh::point_type> > globalPoints(renumbering);
-
- ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering,
- sendParallelMeshOverlap,
- recvParallelMeshOverlap);
- _boundaryMesh->setCalculatedOverlap(true);
-} // _createBoundaryMesh
-
-
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.hh 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -26,11 +26,13 @@
namespace pylith {
namespace bc {
class DirichletBoundary;
+
class TestDirichletBoundary; // unit testing
} // bc
namespace topology {
- class FieldUniform; // USES FieldUniform
+ class SubMesh; // HASA SubMesh
+ class FieldSubMesh; // HASA FieldSubMesh
class SolutionFields; // USES SolutionFields
} // topology
} // pylith
@@ -61,7 +63,7 @@
*
* @return Boundary mesh.
*/
- const ALE::Obj<SieveSubMesh>& boundaryMesh(void) const;
+ const topology::SubMesh& boundaryMesh(void) const;
/** Get vertex field with BC information.
*
@@ -72,20 +74,11 @@
*
* @returns Field over vertices.
*/
- const topology::Field&
+ const topology::FieldSubMesh&
vertexField(const char* name,
const topology::Mesh& mesh,
const topology::SolutionFields& fields);
- // PRIVATE METHODS ////////////////////////////////////////////////////
-private :
-
- /** Extract submesh associated with boundary.
- *
- * @param mesh Finite-element mesh.
- */
- void _createBoundaryMesh(const topology::Mesh& mesh);
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
@@ -98,8 +91,8 @@
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
- ALE::Obj<SieveSubMesh> _boundaryMesh; ///< Boundary mesh.
- topology::FieldUniform* _tmpField; ///< Temporary field for output.
+ topology::SubMesh* _boundaryMesh; ///< Boundary mesh.
+ topology::FieldSubMesh* _tmpField; ///< Temporary field for output.
}; // class DirichletBoundary
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.icc 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/bc/DirichletBoundary.icc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -16,9 +16,9 @@
// Get boundary mesh.
inline
-const ALE::Obj<pylith::SieveSubMesh>&
+const pylith::topology::SubMesh&
pylith::bc::DirichletBoundary::boundaryMesh(void) const {
- return _boundaryMesh;
+ return *_boundaryMesh;
} // boundaryMesh
Deleted: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -1,346 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "Field.hh" // implementation of class methods
-
-#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 ALE::Obj<SieveMesh>& mesh) :
- _mesh(mesh),
- _scale(1.0),
- _name("unknown"),
- _vecFieldType(OTHER),
- _dimensionsOkay(false)
-{ // constructor
- assert(!mesh.isNull());
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor.
-pylith::topology::Field::~Field(void)
-{ // destructor
-} // destructor
-
-// ----------------------------------------------------------------------
-// Get spatial dimension of domain.
-int
-pylith::topology::Field::spaceDim(void) const
-{ // spaceDim
- assert(!_mesh.isNull());
- return _mesh->getDimension();
-} // spaceDim
-
-// ----------------------------------------------------------------------
-// Create seive section.
-void
-pylith::topology::Field::newSection(void)
-{ // newSection
- assert(!_mesh.isNull());
- _section = new SieveRealSection(_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
-
- assert(!_mesh.isNull());
- _section = new SieveRealSection(_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(SieveRealSection::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
- ALE::Obj<SieveMesh::label_sequence> points;
- if (VERTICES_FIELD == domain)
- points = _mesh->depthStratum(0);
- else if (CELLS_FIELD == domain)
- points = _mesh->heightStratum(1);
- else {
- std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
- assert(0);
- } // else
-
- newSection(points, fiberDim);
-} // newSection
-
-// ----------------------------------------------------------------------
-// Create section given atlas.
-void
-pylith::topology::Field::copyLayout(const Field& src)
-{ // createSection
- _vecFieldType = src._vecFieldType;
-
- const ALE::Obj<SieveRealSection>& srcSection = src.section();
- if (!srcSection.isNull() && _section.isNull())
- newSection();
-
- if (!_section.isNull()) {
- _section->setAtlas(srcSection->getAtlas());
- _section->allocateStorage();
- _section->setBC(srcSection->getBC());
- } // if
-} // createSection
-
-// ----------------------------------------------------------------------
-// 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());
-
- _mesh->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
- if (!_section.isNull())
- ALE::Completion::completeSectionAdd(_mesh->getSendOverlap(),
- _mesh->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 SieveRealSection::chart_type& chart = _section->getChart();
- const SieveRealSection::chart_type::const_iterator chartEnd = chart.end();
-
- for (SieveRealSection::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 SieveRealSection::chart_type& chart = _section->getChart();
- const SieveRealSection::chart_type::const_iterator chartEnd = chart.end();
-
- // Assume fiber dimension is uniform
- const int fiberDim = _section->getFiberDimension(*chart.begin());
- double_array values(fiberDim);
-
- for (SieveRealSection::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 SieveRealSection::chart_type& chart = _section->getChart();
- const SieveRealSection::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 (SieveRealSection::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
Deleted: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -1,227 +0,0 @@
-// -*- 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.
- */
-
-#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
-
-#include <string> // HASA std::string
-
-// Forward declarations -------------------------------------------------
-namespace pylith {
- namespace topology {
- class Field;
- class TestField;
- } // topology
-} // pylith
-
-// Field ----------------------------------------------------------------
-class pylith::topology::Field
-{ // Field
- friend class TestField; // unit testing
-
-// PUBLIC ENUMS /////////////////////////////////////////////////////////
-public :
-
- enum VectorFieldEnum {
- SCALAR=0, ///< Scalar.
- VECTOR=1, ///< Vector.
- TENSOR=2, ///< Tensor.
- OTHER=3, ///< Not a scalar, vector, or tensor.
- MULTI_SCALAR=4, ///< Scalar at multiple points.
- MULTI_VECTOR=5, ///< Vector at multiple points.
- MULTI_TENSOR=6, ///< Tensor at multiple points.
- MULTI_OTHER=7, ///< Not a scalar, vector, or tensor at multiple points.
- }; // VectorFieldEnum
-
- enum DomainEnum {
- VERTICES_FIELD=0, ///< Field over vertices.
- CELLS_FIELD=1, ///< Field over cells.
- }; // omainEnum
-
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
-public :
-
- /** Default constructor.
- *
- * @param mesh Sieve mesh.
- */
- Field(const ALE::Obj<SieveMesh>& mesh);
-
- /// Destructor.
- ~Field(void);
-
- /** Get Sieve section.
- *
- * @returns Sieve section.
- */
- const ALE::Obj<SieveRealSection>& section(void) const;
-
- /** Set name of field.
- *
- * @param value Name of field.
- */
- void name(const char* value);
-
- /** Get name of field.
- *
- * @returns Name of field.
- */
- const char* name(void) const;
-
- /** Set vector field type
- *
- * @param value Type of vector field.
- */
- void vectorFieldType(const VectorFieldEnum value);
-
- /** Get vector field type
- *
- * @returns Type of vector field.
- */
- VectorFieldEnum vectorFieldType(void) const;
-
- /** Get spatial dimension of domain.
- *
- * @returns Spatial dimension of domain.
- */
- int spaceDim(void) const;
-
- /** Set scale for dimensionalizing field.
- *
- * @param value Scale associated with field.
- */
- void scale(const double value);
-
- /** Get scale for dimensionalizing field.
- *
- * @returns Scale associated with field.
- */
- double scale(void) const;
-
- /** Set flag indicating whether it is okay to dimensionalize field.
- *
- * @param value True if it is okay to dimensionalize field.
- */
- void addDimensionOkay(const bool value);
-
- /** Set flag indicating whether it is okay to dimensionalize field.
- *
- * @param value True if it is okay to dimensionalize field.
- */
- bool addDimensionOkay(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 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 copyLayout(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 ALE::Obj<SieveMesh>& _mesh; ///< Mesh associated with section
- ALE::Obj<SieveRealSection> _section; ///< Real section with data
-
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
- double _scale; ///< Dimensional scale associated with field
- std::string _name; ///< Name of field
- VectorFieldEnum _vecFieldType; ///< Type of vector field
- bool _dimensionsOkay; ///< Flag indicating it is okay to dimensionalize
-
-// 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
Deleted: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.icc 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.icc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -1,83 +0,0 @@
-// -*- 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::SieveRealSection>&
-pylith::topology::Field::section(void) const {
- return _section;
-}
-
-// Set name of field.
-inline
-void
-pylith::topology::Field::name(const char* value) {
- _name = value;
-}
-
-// Get name of field.
-inline
-const char*
-pylith::topology::Field::name(void) const {
- return _name.c_str();
-}
-
-// Set vector field type
-inline
-void
-pylith::topology::Field::vectorFieldType(const VectorFieldEnum value) {
- _vecFieldType = value;
-}
-
-// Get vector field type
-inline
-pylith::topology::Field::VectorFieldEnum
-pylith::topology::Field::vectorFieldType(void) const {
- return _vecFieldType;
-}
-
-// Set scale for dimensionalizing field.
-inline
-void
-pylith::topology::Field::scale(const double value) {
- _scale = value;
-}
-
-// Get scale for dimensionalizing field.
-inline
-double
-pylith::topology::Field::scale(void) const {
- return _scale;
-}
-
-// Set flag indicating whether it is okay to dimensionalize field.
-inline
-void
-pylith::topology::Field::addDimensionOkay(const bool value) {
- _dimensionsOkay = value;
-}
-
-// Set flag indicating whether it is okay to dimensionalize field.
-inline
-bool
-pylith::topology::Field::addDimensionOkay(void) const {
- return _dimensionsOkay;
-}
-
-#endif
-
-
-// End of file
Copied: short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.cc (from rev 13957, short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc)
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.cc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,34 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "FieldBase.hh" // implementation of class methods
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::topology::FieldBase::FieldBase(void) :
+ _scale(1.0),
+ _name("unknown"),
+ _vecFieldType(OTHER),
+ _dimensionsOkay(false)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::topology::FieldBase::~FieldBase(void)
+{ // destructor
+} // destructor
+
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.hh (from rev 13957, short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh)
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,149 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/topology/FieldBase.hh
+ *
+ * @brief Vector field over the vertices or cells of a finite-element
+ * mesh or subset of a finite-element mesh.
+ *
+ * Extends Sieve real general section by adding metadata.
+ *
+ * We could replace FieldBase, Field, and FieldSubMesh (which
+ * implement a field over a finite-element and subset of the
+ * finite-element mesh) by templating Field over the mesh, but this
+ * would not permit as much insulation against the underlying Sieve
+ * implementation. Separation into Field and FieldSubMesh results in
+ * complete insulation of objects using Field and FieldSubMesh from
+ * the Sieve implementation.
+ */
+
+#if !defined(pylith_topology_fieldbase_hh)
+#define pylith_topology_fieldbase_hh
+
+// Include directives ---------------------------------------------------
+#define NEWPYLITHMESH 1
+#include "pylith/utils/sievetypes.hh" // HASA PETSc real_section_type
+
+#include <string> // HASA std::string
+
+// Forward declarations -------------------------------------------------
+namespace pylith {
+ namespace topology {
+ class FieldBase;
+ class TestFieldBase;
+ } // topology
+} // pylith
+
+// FieldBase ----------------------------------------------------------------
+class pylith::topology::FieldBase
+{ // FieldBase
+ friend class TestFieldBase; // unit testing
+
+// PUBLIC ENUMS /////////////////////////////////////////////////////////
+public :
+
+ enum VectorFieldEnum {
+ SCALAR=0, ///< Scalar.
+ VECTOR=1, ///< Vector.
+ TENSOR=2, ///< Tensor.
+ OTHER=3, ///< Not a scalar, vector, or tensor.
+ MULTI_SCALAR=4, ///< Scalar at multiple points.
+ MULTI_VECTOR=5, ///< Vector at multiple points.
+ MULTI_TENSOR=6, ///< Tensor at multiple points.
+ MULTI_OTHER=7, ///< Not a scalar, vector, or tensor at multiple points.
+ }; // VectorFieldEnum
+
+ enum DomainEnum {
+ VERTICES_FIELD=0, ///< FieldBase over vertices.
+ CELLS_FIELD=1, ///< FieldBase over cells.
+ }; // omainEnum
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ FieldBase(void);
+
+ /// Destructor.
+ ~FieldBase(void);
+
+ /** Set name of field.
+ *
+ * @param value Name of field.
+ */
+ void name(const char* value);
+
+ /** Get name of field.
+ *
+ * @returns Name of field.
+ */
+ const char* name(void) const;
+
+ /** Set vector field type
+ *
+ * @param value Type of vector field.
+ */
+ void vectorFieldType(const VectorFieldEnum value);
+
+ /** Get vector field type
+ *
+ * @returns Type of vector field.
+ */
+ VectorFieldEnum vectorFieldType(void) const;
+
+ /** Set scale for dimensionalizing field.
+ *
+ * @param value Scale associated with field.
+ */
+ void scale(const double value);
+
+ /** Get scale for dimensionalizing field.
+ *
+ * @returns Scale associated with field.
+ */
+ double scale(void) const;
+
+ /** Set flag indicating whether it is okay to dimensionalize field.
+ *
+ * @param value True if it is okay to dimensionalize field.
+ */
+ void addDimensionOkay(const bool value);
+
+ /** Set flag indicating whether it is okay to dimensionalize field.
+ *
+ * @param value True if it is okay to dimensionalize field.
+ */
+ bool addDimensionOkay(void) const;
+
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ double _scale; ///< Dimensional scale associated with field
+ std::string _name; ///< Name of field
+ VectorFieldEnum _vecFieldType; ///< Type of vector field
+ bool _dimensionsOkay; ///< Flag indicating it is okay to dimensionalize
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ FieldBase(const FieldBase&); ///< Not implemented
+ const FieldBase& operator=(const FieldBase&); ///< Not implemented
+
+}; // FieldBase
+
+#include "FieldBase.icc"
+
+#endif // pylith_topology_fieldbase_hh
+
+
+// End of file
Copied: short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.icc (from rev 13957, short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.icc)
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.icc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldBase.icc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,76 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_topology_fieldbase_hh)
+#error "FieldBase.icc must be included only from FieldBase.hh"
+#else
+
+// Set name of field.
+inline
+void
+pylith::topology::FieldBase::name(const char* value) {
+ _name = value;
+}
+
+// Get name of field.
+inline
+const char*
+pylith::topology::FieldBase::name(void) const {
+ return _name.c_str();
+}
+
+// Set vector field type
+inline
+void
+pylith::topology::FieldBase::vectorFieldType(const VectorFieldEnum value) {
+ _vecFieldType = value;
+}
+
+// Get vector field type
+inline
+pylith::topology::FieldBase::VectorFieldEnum
+pylith::topology::FieldBase::vectorFieldType(void) const {
+ return _vecFieldType;
+}
+
+// Set scale for dimensionalizing field.
+inline
+void
+pylith::topology::FieldBase::scale(const double value) {
+ _scale = value;
+}
+
+// Get scale for dimensionalizing field.
+inline
+double
+pylith::topology::FieldBase::scale(void) const {
+ return _scale;
+}
+
+// Set flag indicating whether it is okay to dimensionalize field.
+inline
+void
+pylith::topology::FieldBase::addDimensionOkay(const bool value) {
+ _dimensionsOkay = value;
+}
+
+// Set flag indicating whether it is okay to dimensionalize field.
+inline
+bool
+pylith::topology::FieldBase::addDimensionOkay(void) const {
+ return _dimensionsOkay;
+}
+
+#endif
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.cc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,368 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "FieldSubMesh.hh" // implementation of class methods
+
+#include "SubMesh.hh" // HASA SubMesh
+#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::FieldSubMesh::FieldSubMesh(const SubMesh& mesh) :
+ _mesh(mesh)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::topology::FieldSubMesh::~FieldSubMesh(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Get spatial dimension of domain.
+int
+pylith::topology::FieldSubMesh::spaceDim(void) const
+{ // spaceDim
+ const spatialdata::geocoords::CoordSys* cs = _mesh.coordsys();
+ return (0 != cs) ? cs->spaceDim() : 0;
+} // spaceDim
+
+// ----------------------------------------------------------------------
+// Create seive section.
+void
+pylith::topology::FieldSubMesh::newSection(void)
+{ // newSection
+ _section = new SubMeshRealSection(_mesh.comm(), _mesh.debug());
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion.
+void
+pylith::topology::FieldSubMesh::newSection(
+ const ALE::Obj<SieveSubMesh::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 SubMeshRealSection(_mesh.comm(), _mesh.debug());
+
+ const SieveSubMesh::point_type pointMin =
+ *std::min_element(points->begin(), points->end());
+ const SieveSubMesh::point_type pointMax =
+ *std::max_element(points->begin(), points->end());
+ _section->setChart(SubMeshRealSection::chart_type(pointMin, pointMax+1));
+ _section->setFiberDimension(points, fiberDim);
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion.
+void
+pylith::topology::FieldSubMesh::newSection(const DomainEnum domain,
+ const int fiberDim)
+{ // newSection
+ const ALE::Obj<SieveSubMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ ALE::Obj<SieveSubMesh::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::FieldSubMesh::newSection(
+ const SubMeshRealSection::chart_type& chart,
+ const int fiberDim)
+{ // newSection
+ if (_section.isNull())
+ FieldSubMesh::newSection();
+
+ _section->setChart(chart);
+
+ const SubMeshRealSection::chart_type::const_iterator chartEnd = chart.end();
+ for (SubMeshRealSection::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::FieldSubMesh::newSection(const FieldSubMesh& src)
+{ // newSection
+ _vecFieldType = src._vecFieldType;
+
+ const ALE::Obj<SubMeshRealSection>& 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::FieldSubMesh::clear(void)
+{ // clear
+ if (!_section.isNull())
+ _section->clear();
+
+ _scale = 1.0;
+ _vecFieldType = OTHER;
+ _dimensionsOkay = false;
+} // clear
+
+// ----------------------------------------------------------------------
+// Allocate Sieve section.
+void
+pylith::topology::FieldSubMesh::allocate(void)
+{ // allocate
+ assert(!_section.isNull());
+
+ const ALE::Obj<SieveSubMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ sieveMesh->allocate(_section);
+} // allocate
+
+// ----------------------------------------------------------------------
+// Zero section values.
+void
+pylith::topology::FieldSubMesh::zero(void)
+{ // zero
+ if (!_section.isNull())
+ _section->zero();
+} // zero
+
+// ----------------------------------------------------------------------
+// Complete section by assembling across processors.
+void
+pylith::topology::FieldSubMesh::complete(void)
+{ // complete
+ const ALE::Obj<SieveSubMesh>& 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::FieldSubMesh::copy(const FieldSubMesh& 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 SubMeshRealSection::chart_type& chart = _section->getChart();
+ const SubMeshRealSection::chart_type::const_iterator chartEnd = chart.end();
+
+ for (SubMeshRealSection::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::FieldSubMesh::operator+=(const FieldSubMesh& 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 SubMeshRealSection::chart_type& chart = _section->getChart();
+ const SubMeshRealSection::chart_type::const_iterator chartEnd = chart.end();
+
+ // Assume fiber dimension is uniform
+ const int fiberDim = _section->getFiberDimension(*chart.begin());
+ double_array values(fiberDim);
+
+ for (SubMeshRealSection::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::FieldSubMesh::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 SubMeshRealSection::chart_type& chart = _section->getChart();
+ const SubMeshRealSection::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 (SubMeshRealSection::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::FieldSubMesh::view(const char* label)
+{ // view
+ std::string vecFieldSubMeshString;
+ switch(_vecFieldType)
+ { // switch
+ case SCALAR:
+ vecFieldSubMeshString = "scalar";
+ break;
+ case VECTOR:
+ vecFieldSubMeshString = "vector";
+ break;
+ case TENSOR:
+ vecFieldSubMeshString = "tensor";
+ break;
+ case OTHER:
+ vecFieldSubMeshString = "other";
+ break;
+ case MULTI_SCALAR:
+ vecFieldSubMeshString = "multiple scalars";
+ break;
+ case MULTI_VECTOR:
+ vecFieldSubMeshString = "multiple vectors";
+ break;
+ case MULTI_TENSOR:
+ vecFieldSubMeshString = "multiple tensors";
+ break;
+ case MULTI_OTHER:
+ vecFieldSubMeshString = "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: " << vecFieldSubMeshString << "\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/FieldSubMesh.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,161 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/topology/FieldSubMesh.hh
+ *
+ * @brief Vector field over the vertices or cells of a lower dimension
+ * portion of a finite-element mesh.
+ *
+ * Extends Sieve real general section by adding metadata.
+ */
+
+#include "FieldBase.hh" // ISA FieldBase
+
+#if !defined(pylith_topology_fieldsubmesh_hh)
+#define pylith_topology_fieldsubmesh_hh
+
+// Include directives ---------------------------------------------------
+#define NEWPYLITHMESH 1
+#include "pylith/utils/sievetypes.hh" // HASA PETSc real_section_type
+
+// Forward declarations -------------------------------------------------
+namespace pylith {
+ namespace topology {
+ class FieldSubMesh;
+ class TestFieldSubMesh;
+
+ class SubMesh; // HASA SubMesh
+ } // topology
+} // pylith
+
+// FieldSubMesh ----------------------------------------------------------------
+class pylith::topology::FieldSubMesh : public FieldBase
+{ // FieldSubMesh
+ friend class TestFieldSubMesh; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Default constructor.
+ *
+ * @param mesh Lower dimension finite-element mesh.
+ */
+ FieldSubMesh(const SubMesh& mesh);
+
+ /// Destructor.
+ ~FieldSubMesh(void);
+
+ /** Get Sieve section.
+ *
+ * @returns Sieve section.
+ */
+ const ALE::Obj<SubMeshRealSection>& 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<SieveSubMesh::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 SubMeshRealSection::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 FieldSubMesh& 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 FieldSubMesh to copy.
+ */
+ void copy(const FieldSubMesh& field);
+
+ /** Add two fields, storing the result in one of the fields.
+ *
+ * @param field FieldSubMesh to add.
+ */
+ void operator+=(const FieldSubMesh& 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 SubMesh& _mesh; ///< Mesh associated with section
+ ALE::Obj<SubMeshRealSection> _section; ///< Real section with data
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ FieldSubMesh(const FieldSubMesh&); ///< Not implemented
+ const FieldSubMesh& operator=(const FieldSubMesh&); ///< Not implemented
+
+}; // FieldSubMesh
+
+#include "FieldSubMesh.icc"
+
+#endif // pylith_topology_fieldsubmesh_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.icc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldSubMesh.icc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,27 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_topology_fieldsubmesh_hh)
+#error "FieldSubMesh.icc must be included only from FieldSubMesh.hh"
+#else
+
+// Get Sieve section.
+inline
+const ALE::Obj<pylith::SubMeshRealSection>&
+pylith::topology::FieldSubMesh::section(void) const {
+ return _section;
+}
+
+#endif
+
+
+// End of file
Deleted: short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldUniform.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldUniform.cc 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldUniform.cc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -1,67 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-#include <portinfo>
-
-#include "FieldUniform.hh" // implementation of class methods
-
-#include <stdexcept> // USES std::runtime_error
-#include <sstream> // USES std::ostringstream
-#include <cassert> // USES assert()
-
-// ----------------------------------------------------------------------
-// Default constructor.
-pylith::topology::FieldUniform::FieldUniform(const ALE::Obj<SieveMesh>& mesh,
- const int fiberDim) :
- Field(mesh),
- _fiberDim(fiberDim)
-{ // constructor
- assert(fiberDim >= 0);
-} // constructor
-
-// ----------------------------------------------------------------------
-// Destructor.
-pylith::topology::FieldUniform::~FieldUniform(void)
-{ // destructor
-} // destructor
-
-// ----------------------------------------------------------------------
-// Create section given points.
-void
-pylith::topology::FieldUniform::createSection(
- const ALE::Obj<SieveMesh::label_sequence>& points)
-{ // createSection
- newSection(points, _fiberDim);
- _mesh->allocate(_section);
-} // createSection
-
-// ----------------------------------------------------------------------
-// Create section given chart.
-void
-pylith::topology::FieldUniform::createSection(
- const SieveRealSection::chart_type& chart)
-{ // createSection
- if (_section.isNull())
- newSection();
-
- _section->setChart(chart);
-
- const SieveRealSection::chart_type::const_iterator chartEnd = chart.end();
- for (SieveRealSection::chart_type::const_iterator c_iter = chart.begin();
- c_iter != chartEnd;
- ++c_iter)
- _section->setFiberDimension(*c_iter, _fiberDim);
- _mesh->allocate(_section);
-} // createSection
-
-
-// End of file
Deleted: short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldUniform.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldUniform.hh 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/FieldUniform.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -1,83 +0,0 @@
-// -*- C++ -*-
-//
-// ======================================================================
-//
-// Brad T. Aagaard
-// U.S. Geological Survey
-//
-// {LicenseText}
-//
-// ======================================================================
-//
-
-/**
- * @file pylith/topology/FieldUniform.hh
- *
- * @brief Vector field with uniform fiber dimension (no constraints)
- * over the vertices or cells of a finite-element mesh.
- *
- * Special case of Field with a uniform fiber dimension and no constraints.
- */
-
-#if !defined(pylith_topology_fielduniform_hh)
-#define pylith_topology_fielduniform_hh
-
-// Include directives ---------------------------------------------------
-#include "Field.hh" // ISA Field
-
-// Forward declarations -------------------------------------------------
-namespace pylith {
- namespace topology {
- class FieldUniform;
- class TestFieldUniform;
- } // topology
-} // pylith
-
-// FieldUniform ---------------------------------------------------------
-class pylith::topology::FieldUniform : public Field
-{ // FieldUniform
- friend class TestFieldUniform; // unit testing
-
-// PUBLIC MEMBERS ///////////////////////////////////////////////////////
-public :
-
- /** Default constructor.
- *
- * @param mesh Sieve mesh.
- */
- FieldUniform(const ALE::Obj<SieveMesh>& mesh,
- const int fiberDim);
-
- /// Destructor.
- ~FieldUniform(void);
-
- /** Create section given points.
- *
- * @param points Mesh points over which to define section.
- */
- void createSection(const ALE::Obj<SieveMesh::label_sequence>& points);
-
- /** 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.
- */
- void createSection(const SieveRealSection::chart_type& chart);
-
-// PRIVATE MEMBERS //////////////////////////////////////////////////////
-private :
-
- const int _fiberDim; ///< Fiber dimension
-
-// NOT IMPLEMENTED //////////////////////////////////////////////////////
-private :
-
- FieldUniform(const FieldUniform&); ///< Not implemented
- const FieldUniform& operator=(const FieldUniform&); ///< Not implemented
-
-}; // FieldUniform
-
-#endif // pylith_topology_fielduniform_hh
-
-
-// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Makefile.am 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Makefile.am 2009-01-27 05:22:03 UTC (rev 13960)
@@ -15,16 +15,21 @@
subpkginclude_HEADERS = \
Distributor.hh \
+ FieldBase.hh \
+ FieldBase.icc \
Field.hh \
Field.icc \
- FieldUniform.hh \
+ FieldSubMesh.hh \
+ FieldSubMesh.icc \
FieldsManager.hh \
FieldOps.hh \
Mesh.hh \
Mesh.icc \
MeshOps.hh \
MeshRefiner.hh \
- RefineUniform.hh
+ RefineUniform.hh \
+ SubMesh.hh \
+ SubMesh.icc
noinst_HEADERS =
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.cc 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.cc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -70,51 +70,5 @@
_coordsys->initialize();
} // initialize
-// ----------------------------------------------------------------------
-// Create boundary mesh from mesh and label.
-void
-pylith::topology::Mesh::boundaryMeshFromMesh(const char* label)
-{ // boundaryMeshFromMesh
-#if 0
- assert(!_mesh.isNull());
- const ALE::Obj<SieveMesh::int_section_type>& groupField =
- _mesh->getIntSection(_label);
- if (groupField.isNull()) {
- std::ostringstream msg;
- msg << "Could not find group of points '" << label << "' in mesh.";
- throw std::runtime_error(msg.str());
- } // if
- ALE::Obj<SubSieveMesh> boundarySieveMesh =
- ALE::Selection<SieveMesh>::submeshV<SieveSubMesh>(_mesh, groupField);
- if (boundarySieveMesh.isNull()) {
- std::ostringstream msg;
- msg << "Could not construct boundary mesh for boundary '"
- << label << "'.";
- throw std::runtime_error(msg.str());
- } // if
- boundarySieveMesh->setRealSection("coordinates",
- _mesh->getRealSection("coordinates"));
-
- // Create the parallel overlap
- ALE::Obj<SieveSubMesh::send_overlap_type> sendParallelMeshOverlap =
- boundarySieveMesh->getSendOverlap();
- ALE::Obj<SieveSubMesh::recv_overlap_type> recvParallelMeshOverlap =
- boundarySieveMesh->getRecvOverlap();
- SieveMesh::renumbering_type& renumbering = _mesh->getRenumbering();
- // Can I figure this out in a nicer way?
- ALE::SetFromMap<std::map<SieveMesh::point_type,SieveMesh::point_type> > globalPoints(renumbering);
-
- ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering,
- sendParallelMeshOverlap,
- recvParallelMeshOverlap);
- boundarySieveMesh->setCalculatedOverlap(true);
-
- Mesh boundaryMesh;
-
- return boundaryMesh
-#endif
-} // createBoundaryMesh
-
-
// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -131,12 +131,6 @@
*/
void view(const char* label);
- /** Create boundary mesh from mesh and label.
- *
- * @param label Label for group of vertices forming boundary.
- */
- Mesh createBoundaryMesh(const char* label);
-
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.icc 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Mesh.icc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -55,7 +55,7 @@
inline
int
pylith::topology::Mesh::dimension(void) const {
- return (_mesh.isNull()) ? 0 : _mesh->getDimension();
+ return (!_mesh.isNull()) ? _mesh->getDimension() : 0;
}
// Set MPI communicator associated with mesh.
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.cc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,109 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "SubMesh.hh" // implementation of class methods
+
+#include "Mesh.hh" // USES Mesh
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
+
+#include <Selection.hh> // USES ALE::Selection
+
+// ----------------------------------------------------------------------
+// Default constructor
+pylith::topology::SubMesh::SubMesh(void) :
+ _coordsys(0),
+ _debug(false)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Constructor with mesh and label for vertices marking boundary.
+pylith::topology::SubMesh::SubMesh(const Mesh& mesh,
+ const char* label) :
+ _coordsys(0),
+ _debug(false)
+{ // constructor
+ createSubMesh(mesh, label);
+} // constructor
+
+// ----------------------------------------------------------------------
+// Default destructor
+pylith::topology::SubMesh::~SubMesh(void)
+{ // destructor
+ delete _coordsys; _coordsys = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Create Sieve mesh.
+void
+pylith::topology::SubMesh::createSubMesh(const Mesh& mesh,
+ const char* label)
+{ // createSieveMesh
+ _mesh.destroy();
+
+ const ALE::Obj<SieveMesh>& meshSieveMesh = mesh.sieveMesh();
+ assert(!meshSieveMesh.isNull());
+
+ const ALE::Obj<SieveMesh::int_section_type>& groupField =
+ meshSieveMesh->getIntSection(label);
+ if (groupField.isNull()) {
+ std::ostringstream msg;
+ msg << "Could not find group of points '" << label << "' in mesh.";
+ throw std::runtime_error(msg.str());
+ } // if
+ _mesh =
+ ALE::Selection<SieveMesh>::submeshV<SieveSubMesh>(meshSieveMesh,
+ groupField);
+ if (_mesh.isNull()) {
+ std::ostringstream msg;
+ msg << "Could not construct boundary mesh for boundary '"
+ << label << "'.";
+ throw std::runtime_error(msg.str());
+ } // if
+ _mesh->setRealSection("coordinates",
+ meshSieveMesh->getRealSection("coordinates"));
+
+ // Create the parallel overlap
+ ALE::Obj<SieveSubMesh::send_overlap_type> sendParallelMeshOverlap =
+ _mesh->getSendOverlap();
+ ALE::Obj<SieveSubMesh::recv_overlap_type> recvParallelMeshOverlap =
+ _mesh->getRecvOverlap();
+ SieveMesh::renumbering_type& renumbering = meshSieveMesh->getRenumbering();
+ // Can I figure this out in a nicer way?
+ ALE::SetFromMap<std::map<SieveMesh::point_type,SieveMesh::point_type> > globalPoints(renumbering);
+
+ ALE::OverlapBuilder<>::constructOverlap(globalPoints, renumbering,
+ sendParallelMeshOverlap,
+ recvParallelMeshOverlap);
+ _mesh->setCalculatedOverlap(true);
+
+ // Set data from mesh.
+ _mesh->setDebug(mesh.debug());
+ delete _coordsys; _coordsys = 0;
+ const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
+ if (0 != cs)
+ _coordsys = cs->clone();
+} // createSubMesh
+
+// ----------------------------------------------------------------------
+// Initialize the finite-element mesh.
+void
+pylith::topology::SubMesh::initialize(void)
+{ // initialize
+ if (0 != _coordsys)
+ _coordsys->initialize();
+} // initialize
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,146 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/SubMesh.hh
+ *
+ * @brief C++ PyLith finite-element mesh.
+ *
+ * Extends Sieve mesh over subset of domain to include coordinate
+ * system associated with domain.
+ */
+
+#if !defined(pylith_topology_submesh_hh)
+#define pylith_topology_submesh_hh
+
+// Include directives ---------------------------------------------------
+#define NEWPYLITHMESH 1
+#include "pylith/utils/sievetypes.hh"
+
+// Forward declarations -------------------------------------------------
+namespace pylith {
+ namespace topology {
+ class SubMesh;
+ class TestSubMesh; // unit testing
+
+ class Mesh; // USES Mesh
+ } // topology
+} // pylith
+
+namespace spatialdata {
+ namespace geocoords {
+ class CoordSys; // USES Mesh
+ } // geocoords
+} // spatialdata
+
+// SubMesh -----------------------------------------------------------------
+class pylith::topology::SubMesh
+{ // SubMesh
+ friend class TestSubMesh; // unit testing
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public :
+
+ /// Default constructor.
+ SubMesh(void);
+
+ /** Constructor with mesh and label for vertices marking boundary.
+ *
+ * @param mesh Finite-element mesh over domain.
+ * @param label Label for vertices marking boundary.
+ */
+ SubMesh(const Mesh& mesh,
+ const char* label);
+
+ /// Default destructor
+ ~SubMesh(void);
+
+ /** Create Sieve mesh.
+ *
+ * @param mesh Finite-element mesh over domain.
+ * @param label Label for vertices marking boundary.
+ */
+ void createSubMesh(const Mesh& mesh,
+ const char* label);
+
+ /** Get Sieve mesh.
+ *
+ * @returns Sieve mesh.
+ */
+ const ALE::Obj<SieveSubMesh>& sieveMesh(void) const;
+
+ /** Get Sieve mesh.
+ *
+ * @returns Sieve mesh.
+ */
+ ALE::Obj<SieveSubMesh>& sieveMesh(void);
+
+ /** Get coordinate system.
+ *
+ * @returns Coordinate system.
+ */
+ const spatialdata::geocoords::CoordSys* coordsys(void) const;
+
+ /** Set debug flag.
+ *
+ * @param value Turn on debugging if true.
+ */
+ void debug(const bool value);
+
+ /** Get debug flag.
+ *
+ * @param Get debugging flag.
+ */
+ bool debug(void) const;
+
+ /** Get dimension of mesh.
+ *
+ * @returns Dimension of mesh.
+ */
+ int dimension(void) const;
+
+ /** Get MPI communicator associated with mesh.
+ *
+ * @returns MPI communicator.
+ */
+ const MPI_Comm comm(void) const;
+
+ /// Initialize the finite-element mesh.
+ void initialize(void);
+
+ /** Print mesh to stdout.
+ *
+ * @param label Label for mesh.
+ */
+ void view(const char* label);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ ALE::Obj<SieveSubMesh> _mesh; ///< Sieve mesh.
+ spatialdata::geocoords::CoordSys* _coordsys; ///< Coordinate system.
+ bool _debug; ///< Debugging flag for mesh.
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ SubMesh(const SubMesh&); ///< Not implemented
+ const SubMesh& operator=(const SubMesh&); ///< Not implemented
+
+}; // SubMesh
+
+#include "SubMesh.icc"
+
+#endif // pylith_topology_submesh_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.icc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.icc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/SubMesh.icc 2009-01-27 05:22:03 UTC (rev 13960)
@@ -0,0 +1,79 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_topology_submesh_hh)
+#error "SubMesh.icc must be included only from SubMesh.hh"
+#else
+
+// Get Sieve mesh.
+inline
+const ALE::Obj<pylith::SieveSubMesh>&
+pylith::topology::SubMesh::sieveMesh(void) const {
+ return _mesh;
+}
+
+// Get Sieve mesh.
+inline
+ALE::Obj<pylith::SieveSubMesh>&
+pylith::topology::SubMesh::sieveMesh(void) {
+ return _mesh;
+}
+
+// Get coordinate system.
+inline
+const spatialdata::geocoords::CoordSys*
+pylith::topology::SubMesh::coordsys(void) const {
+ return _coordsys;
+}
+
+// Set debug flag.
+inline
+void
+pylith::topology::SubMesh::debug(const bool value) {
+ _debug = value;
+ if (!_mesh.isNull())
+ _mesh->setDebug(value);
+}
+
+// Get debug flag.
+inline
+bool
+pylith::topology::SubMesh::debug(void) const {
+ return _debug;
+}
+
+// Get dimension of mesh.
+inline
+int
+pylith::topology::SubMesh::dimension(void) const {
+ return (!_mesh.isNull()) ? _mesh->getDimension() : 0;
+}
+
+// Get MPI communicator associated with mesh.
+inline
+const MPI_Comm
+pylith::topology::SubMesh::comm(void) const {
+ return (!_mesh.isNull()) ? _mesh->comm() : 0;
+}
+
+// Print mesh to stdout.
+inline
+void
+pylith::topology::SubMesh::view(const char* label) {
+ _mesh->view(label);
+}
+
+
+#endif
+
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/utils/sievetypes.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/utils/sievetypes.hh 2009-01-27 03:41:13 UTC (rev 13959)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/utils/sievetypes.hh 2009-01-27 05:22:03 UTC (rev 13960)
@@ -25,9 +25,12 @@
#if NEWPYLITHMESH // For use with pylith::topology::Mesh
typedef ALE::IMesh<> SieveMesh;
+ typedef SieveMesh::real_section_type MeshRealSection;
+ typedef SieveMesh::int_section_type MeshIntSection;
+
typedef ALE::IMesh<ALE::LabelSifter<int, SieveMesh::point_type> > SieveSubMesh;
- typedef SieveMesh::real_section_type SieveRealSection;
- typedef SieveMesh::int_section_type SieveIntSection;
+ typedef SieveSubMesh::real_section_type SubMeshRealSection;
+ typedef SieveSubMesh::int_section_type SubMeshIntSection;
#else
typedef ALE::IMesh<> Mesh;
More information about the CIG-COMMITS
mailing list