[cig-commits] r18415 - in short/3D/PyLith/branches/multifields: libsrc/pylith/topology modulesrc/topology unittests/libtests/topology
brad at geodynamics.org
brad at geodynamics.org
Sat May 21 12:47:17 PDT 2011
Author: brad
Date: 2011-05-21 12:47:16 -0700 (Sat, 21 May 2011)
New Revision: 18415
Added:
short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.cc
short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.hh
short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.icc
short/3D/PyLith/branches/multifields/modulesrc/topology/MultiField.i
short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldMesh.cc
short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldMesh.hh
short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldSubMesh.cc
short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldSubMesh.hh
Modified:
short/3D/PyLith/branches/multifields/libsrc/pylith/topology/Makefile.am
short/3D/PyLith/branches/multifields/libsrc/pylith/topology/topologyfwd.hh
short/3D/PyLith/branches/multifields/modulesrc/topology/Makefile.am
short/3D/PyLith/branches/multifields/modulesrc/topology/topology.i
short/3D/PyLith/branches/multifields/unittests/libtests/topology/Makefile.am
Log:
Created MultiField for having multiple fields in a section.
Modified: short/3D/PyLith/branches/multifields/libsrc/pylith/topology/Makefile.am
===================================================================
--- short/3D/PyLith/branches/multifields/libsrc/pylith/topology/Makefile.am 2011-05-21 19:28:25 UTC (rev 18414)
+++ short/3D/PyLith/branches/multifields/libsrc/pylith/topology/Makefile.am 2011-05-21 19:47:16 UTC (rev 18415)
@@ -25,6 +25,9 @@
Field.hh \
Field.icc \
Field.cc \
+ MultiField.hh \
+ MultiField.cc \
+ MultiField.icc \
Fields.hh \
Fields.icc \
PackedFields.hh \
Added: short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.cc
===================================================================
--- short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.cc (rev 0)
+++ short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.cc 2011-05-21 19:47:16 UTC (rev 18415)
@@ -0,0 +1,1181 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "MultiField.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 "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
+#include <stdexcept> // USES std::runtime_error
+#include <sstream> // USES std::ostringstream
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+// Default constructor.
+template<typename mesh_type, typename section_type>
+pylith::topology::MultiField<mesh_type, section_type>::MultiField(const mesh_type& mesh) :
+ _mesh(mesh)
+{ // constructor
+ _metadata.label = "unknown";
+ _metadata.vectorFieldType = OTHER;
+ _metadata.scale = 1.0;
+ _metadata.dimsOkay = false;
+} // constructor
+
+// ----------------------------------------------------------------------
+// Constructor with mesh, section, and metadata.
+template<typename mesh_type, typename section_type>
+pylith::topology::MultiField<mesh_type, section_type>::MultiField(const mesh_type& mesh,
+ const ALE::Obj<section_type>& section,
+ const Metadata& metadata) :
+ _metadata(metadata),
+ _mesh(mesh),
+ _section(section)
+{ // constructor
+ assert(!section.isNull());
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+template<typename mesh_type, typename section_type>
+pylith::topology::MultiField<mesh_type, section_type>::~MultiField(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate PETSc and local data structures.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::deallocate(void)
+{ // deallocate
+ PetscErrorCode err = 0;
+
+ const typename scatter_map_type::const_iterator scattersEnd = _scatters.end();
+ for (typename scatter_map_type::iterator s_iter=_scatters.begin();
+ s_iter != scattersEnd;
+ ++s_iter) {
+ if (s_iter->second.vector) {
+ err = VecDestroy(&s_iter->second.vector);CHECK_PETSC_ERROR(err);
+ } // if
+
+ if (s_iter->second.scatter) {
+ err = VecScatterDestroy(&s_iter->second.scatter);CHECK_PETSC_ERROR(err);
+ } // if
+
+ if (s_iter->second.scatterVec) {
+ err = VecDestroy(&s_iter->second.scatterVec);CHECK_PETSC_ERROR(err);
+ } // if
+ } // for
+} // deallocate
+
+// ----------------------------------------------------------------------
+// Set label for field.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::label(const char* value)
+{ // label
+ _metadata.label = value;
+ if (!_section.isNull()) {
+ _section->setName(value);
+ } // if
+
+ const typename scatter_map_type::const_iterator scattersEnd = _scatters.end();
+ for (typename scatter_map_type::const_iterator s_iter=_scatters.begin();
+ s_iter != scattersEnd;
+ ++s_iter) {
+ if (s_iter->second.vector) {
+ PetscErrorCode err =
+ PetscObjectSetName((PetscObject)s_iter->second.vector, value);
+ CHECK_PETSC_ERROR(err);
+ } // if
+ } // for
+} // label
+
+// ----------------------------------------------------------------------
+// Get spatial dimension of domain.
+template<typename mesh_type, typename section_type>
+int
+pylith::topology::MultiField<mesh_type, section_type>::spaceDim(void) const
+{ // spaceDim
+ const spatialdata::geocoords::CoordSys* cs = _mesh.coordsys();
+ return (cs) ? cs->spaceDim() : 0;
+} // spaceDim
+
+// ----------------------------------------------------------------------
+// Get the chart size.
+template<typename mesh_type, typename section_type>
+int
+pylith::topology::MultiField<mesh_type, section_type>::chartSize(void) const
+{ // chartSize
+ return _section.isNull() ? 0 : _section->getChart().size();
+} // chartSize
+
+// ----------------------------------------------------------------------
+// Get the number of degrees of freedom.
+template<typename mesh_type, typename section_type>
+int
+pylith::topology::MultiField<mesh_type, section_type>::sectionSize(void) const
+{ // sectionSize
+ return _section.isNull() ? 0 : _section->size();
+} // sectionSize
+
+// ----------------------------------------------------------------------
+// Create seive section.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::newSection(void)
+{ // newSection
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+
+ _section = new section_type(_mesh.comm(), _mesh.debug());
+ assert(!_section.isNull());
+ _section->setName(_metadata.label);
+
+ logger.stagePop();
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion for a
+// sequence of points.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::newSection(
+ const ALE::Obj<label_sequence>& points,
+ const int fiberDim)
+{ // newSection
+ typedef typename mesh_type::SieveMesh::point_type point_type;
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+ if (fiberDim < 0) {
+ std::ostringstream msg;
+ msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata.label
+ << "' must be nonnegative.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ _section = new section_type(_mesh.comm(), _mesh.debug());
+ assert(!_section.isNull());
+ _section->setName(_metadata.label);
+
+ if (points->size() > 0) {
+ const point_type pointMin =
+ *std::min_element(points->begin(), points->end());
+ const point_type pointMax =
+ *std::max_element(points->begin(), points->end());
+ _section->setChart(chart_type(pointMin, pointMax+1));
+ _section->setFiberDimension(points, fiberDim);
+ } else // Create empty chart
+ _section->setChart(chart_type(0, 0));
+
+ logger.stagePop();
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion for a list of
+// points.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::newSection(const int_array& points,
+ const int fiberDim)
+{ // newSection
+ typedef typename mesh_type::SieveMesh::point_type point_type;
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+ if (fiberDim < 0) {
+ std::ostringstream msg;
+ msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata.label
+ << "' must be nonnegative.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ _section = new section_type(_mesh.comm(), _mesh.debug());
+ assert(!_section.isNull());
+ _section->setName(_metadata.label);
+
+ const int npts = points.size();
+ if (npts > 0) {
+ const point_type pointMin = points.min();
+ const point_type pointMax = points.max();
+ _section->setChart(chart_type(pointMin, pointMax+1));
+ for (int i=0; i < npts; ++i)
+ _section->setFiberDimension(points[i], fiberDim);
+ } else // create empty chart
+ _section->setChart(chart_type(0, 0));
+
+ logger.stagePop();
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::newSection(const DomainEnum domain,
+ const int fiberDim,
+ const int stratum)
+{ // newSection
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ ALE::Obj<label_sequence> points;
+ if (VERTICES_FIELD == domain)
+ points = sieveMesh->depthStratum(stratum);
+ else if (CELLS_FIELD == domain)
+ points = sieveMesh->heightStratum(stratum);
+ else {
+ std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
+ assert(0);
+ throw std::logic_error("Bad domain enum in MultiField.");
+ } // else
+
+ newSection(points, fiberDim);
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create section given chart.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::newSection(const MultiField& src,
+ const int fiberDim)
+{ // newSection
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+
+ if (_section.isNull()) {
+ logger.stagePop();
+ newSection();
+ logger.stagePush("Field");
+ } // if
+ if (fiberDim < 0) {
+ std::ostringstream msg;
+ msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata.label
+ << "' must be nonnegative.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ const ALE::Obj<section_type>& srcSection = src.section();
+ if (!srcSection.isNull()) {
+ _section->setChart(srcSection->getChart());
+ const chart_type& chart = _section->getChart();
+ const typename chart_type::const_iterator chartBegin = chart.begin();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+
+ for (typename chart_type::const_iterator c_iter = chartBegin;
+ c_iter != chartEnd;
+ ++c_iter)
+ if (srcSection->getFiberDimension(*c_iter) > 0)
+ _section->setFiberDimension(*c_iter, fiberDim);
+ } // if
+
+ logger.stagePop();
+} // newSection
+
+// ----------------------------------------------------------------------
+// Create section with same layout as another section.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::cloneSection(const MultiField& src)
+{ // cloneSection
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+
+ deallocate();
+ std::string origLabel = _metadata.label;
+ _metadata = src._metadata;
+ label(origLabel.c_str());
+
+ const ALE::Obj<section_type>& srcSection = src.section();
+ if (!srcSection.isNull() && _section.isNull()) {
+ logger.stagePop();
+ newSection();
+ logger.stagePush("Field");
+ }
+
+ if (!_section.isNull()) {
+ if (!srcSection->sharedStorage()) {
+ _section->setAtlas(srcSection->getAtlas());
+ _section->allocateStorage();
+ _section->setBC(srcSection->getBC());
+ _section->copySpaces(srcSection);
+ } else {
+ _section->setChart(srcSection->getChart());
+ const chart_type& chart = _section->getChart();
+ const typename chart_type::const_iterator chartBegin = chart.begin();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+ for (typename chart_type::const_iterator c_iter = chartBegin;
+ c_iter != chartEnd;
+ ++c_iter) {
+ const int fiberDim = srcSection->getFiberDimension(*c_iter);
+ if (fiberDim > 0)
+ _section->setFiberDimension(*c_iter, fiberDim);
+ } // for
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ sieveMesh->allocate(_section);
+ _section->setBC(srcSection->getBC());
+ _section->copySpaces(srcSection); // :BUG: NEED TO REBUILD SPACES
+ } // if/else
+
+ // Reuse scatters in clone
+ PetscErrorCode err = 0;
+ if (src._scatters.size() > 0) {
+ const typename scatter_map_type::const_iterator scattersEnd = src._scatters.end();
+ for (typename scatter_map_type::const_iterator s_iter=src._scatters.begin();
+ s_iter != scattersEnd;
+ ++s_iter) {
+ ScatterInfo sinfo;
+ sinfo.vector = 0;
+ sinfo.scatterVec = 0;
+
+ // Copy scatter
+ sinfo.scatter = s_iter->second.scatter;
+ err = PetscObjectReference((PetscObject) sinfo.scatter);
+ CHECK_PETSC_ERROR(err);
+
+ // Create scatter Vec
+ if (_section->sizeWithBC() > 0) {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF,
+ _section->getStorageSize(),
+ _section->restrictSpace(),
+ &sinfo.scatterVec);
+ CHECK_PETSC_ERROR(err);
+ } else {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, 0, 0,
+ &sinfo.scatterVec);
+ CHECK_PETSC_ERROR(err);
+ } // else
+
+ // Create vector using sizes from source section
+ int vecLocalSize = 0;
+ int vecGlobalSize = 0;
+ err = VecGetLocalSize(s_iter->second.vector, &vecLocalSize);
+ CHECK_PETSC_ERROR(err);
+ err = VecGetSize(s_iter->second.vector, &vecGlobalSize); CHECK_PETSC_ERROR(err);
+
+ err = VecCreate(_mesh.comm(), &sinfo.vector); CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata.label.c_str());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(sinfo.vector, vecLocalSize, vecGlobalSize);
+ CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
+
+ _scatters[s_iter->first] = sinfo;
+ } // for
+ } // if
+ } // if
+ logger.stagePop();
+} // cloneSection
+
+// ----------------------------------------------------------------------
+// Clear variables associated with section.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::clear(void)
+{ // clear
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+
+ deallocate();
+ if (!_section.isNull())
+ _section->clear();
+
+ _metadata.scale = 1.0;
+ _metadata.vectorFieldType = OTHER;
+ _metadata.dimsOkay = false;
+
+ logger.stagePop();
+} // clear
+
+// ----------------------------------------------------------------------
+// Allocate Sieve section.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::allocate(void)
+{ // allocate
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+
+ assert(!_section.isNull());
+
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ sieveMesh->allocate(_section);
+
+ logger.stagePop();
+} // allocate
+
+// ----------------------------------------------------------------------
+// Zero section values (excluding constrained DOF).
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::zero(void)
+{ // zero
+ if (!_section.isNull())
+ _section->zero(); // Does not zero BC.
+} // zero
+
+// ----------------------------------------------------------------------
+// Zero section values (including constrained DOF).
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::zeroAll(void)
+{ // zeroAll
+ if (!_section.isNull()) {
+ const chart_type& chart = _section->getChart();
+ const typename chart_type::const_iterator chartBegin = chart.begin();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+
+ // Assume fiber dimension is uniform
+ const int fiberDim = (chart.size() > 0) ?
+ _section->getFiberDimension(*chartBegin) : 0;
+ double_array values(fiberDim);
+ values *= 0.0;
+
+ for (typename chart_type::const_iterator c_iter = chartBegin;
+ c_iter != chartEnd;
+ ++c_iter) {
+ if (_section->getFiberDimension(*c_iter) > 0) {
+ assert(fiberDim == _section->getFiberDimension(*c_iter));
+ _section->updatePointAll(*c_iter, &values[0]);
+ } // if
+ } // for
+ } // if
+} // zeroAll
+
+// ----------------------------------------------------------------------
+// Complete section by assembling across processors.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::complete(void)
+{ // complete
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Completion");
+
+ const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+
+ if (!_section.isNull())
+ ALE::Completion::completeSectionAdd(sieveMesh->getSendOverlap(),
+ sieveMesh->getRecvOverlap(),
+ _section, _section);
+
+ logger.stagePop();
+} // complete
+
+// ----------------------------------------------------------------------
+// Copy field values and metadata.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::copy(const MultiField& field)
+{ // copy
+ // Check compatibility of sections
+ const int srcSize = field.chartSize();
+ const int dstSize = chartSize();
+ if (field.spaceDim() != spaceDim() ||
+ field._metadata.vectorFieldType != _metadata.vectorFieldType ||
+ srcSize != dstSize) {
+ std::ostringstream msg;
+
+ msg << "Cannot copy values from section '" << field._metadata.label
+ << "' to section '" << _metadata.label
+ << "'. Sections are incompatible.\n"
+ << " Source section:\n"
+ << " space dim: " << field.spaceDim() << "\n"
+ << " vector field type: " << field._metadata.vectorFieldType << "\n"
+ << " scale: " << field._metadata.scale << "\n"
+ << " size: " << srcSize << "\n"
+ << " Destination section:\n"
+ << " space dim: " << spaceDim() << "\n"
+ << " vector field type: " << _metadata.vectorFieldType << "\n"
+ << " scale: " << _metadata.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 chart_type& chart = _section->getChart();
+ const typename chart_type::const_iterator chartBegin = chart.begin();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+
+ for (typename chart_type::const_iterator c_iter = chartBegin;
+ c_iter != chartEnd;
+ ++c_iter) {
+ assert(field._section->getFiberDimension(*c_iter) ==
+ _section->getFiberDimension(*c_iter));
+ if (_section->getFiberDimension(*c_iter) > 0)
+ _section->updatePointAll(*c_iter,
+ field._section->restrictPoint(*c_iter));
+ } // for
+ } // if
+
+ label(field._metadata.label.c_str()); // Update label
+ _metadata.scale = field._metadata.scale;
+} // copy
+
+// ----------------------------------------------------------------------
+// Copy field values.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::copy(const ALE::Obj<section_type>& osection)
+{ // copy
+ // Check compatibility of sections
+ const int srcSize = osection->getChart().size();
+ const int dstSize = chartSize();
+ if (srcSize != dstSize) {
+ std::ostringstream msg;
+
+ msg << "Cannot copy values from Sieve section "
+ << _metadata.label << "'. Sections are incompatible.\n"
+ << " Source section:\n"
+ << " size: " << srcSize << "\n"
+ << " Destination section:\n"
+ << " space dim: " << spaceDim() << "\n"
+ << " vector field type: " << _metadata.vectorFieldType << "\n"
+ << " scale: " << _metadata.scale << "\n"
+ << " size: " << dstSize;
+ throw std::runtime_error(msg.str());
+ } // if
+ assert( (_section.isNull() && osection.isNull()) ||
+ (!_section.isNull() && !osection.isNull()) );
+
+ if (!_section.isNull()) {
+ // Copy values from field
+ const chart_type& chart = _section->getChart();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+
+ for (typename chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter) {
+ assert(osection->getFiberDimension(*c_iter) ==
+ _section->getFiberDimension(*c_iter));
+ if (osection->getFiberDimension(*c_iter))
+ _section->updatePoint(*c_iter, osection->restrictPoint(*c_iter));
+ } // for
+ } // if
+} // copy
+
+// ----------------------------------------------------------------------
+// Add two fields, storing the result in one of the fields.
+template<typename mesh_type, typename section_type>
+pylith::topology::MultiField<mesh_type, section_type>&
+pylith::topology::MultiField<mesh_type, section_type>::operator+=(const MultiField& field)
+{ // operator+=
+ // Check compatibility of sections
+ const int srcSize = field.chartSize();
+ const int dstSize = chartSize();
+ if (field.spaceDim() != spaceDim() ||
+ field._metadata.vectorFieldType != _metadata.vectorFieldType ||
+ field._metadata.scale != _metadata.scale ||
+ srcSize != dstSize) {
+ std::ostringstream msg;
+
+ msg << "Cannot add values from section '" << field._metadata.label
+ << "' to section '" << _metadata.label
+ << "'. Sections are incompatible.\n"
+ << " Source section:\n"
+ << " space dim: " << field.spaceDim() << "\n"
+ << " vector field type: " << field._metadata.vectorFieldType << "\n"
+ << " scale: " << field._metadata.scale << "\n"
+ << " size: " << srcSize << "\n"
+ << " Destination section:\n"
+ << " space dim: " << spaceDim() << "\n"
+ << " vector field type: " << _metadata.vectorFieldType << "\n"
+ << " scale: " << _metadata.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 chart_type& chart = _section->getChart();
+ const typename chart_type::const_iterator chartBegin = chart.begin();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+
+ // Assume fiber dimension is uniform
+ const int fiberDim = (chart.size() > 0) ?
+ _section->getFiberDimension(*chartBegin) : 0;
+ double_array values(fiberDim);
+
+ for (typename chart_type::const_iterator c_iter = chartBegin;
+ c_iter != chartEnd;
+ ++c_iter) {
+ if (field._section->getFiberDimension(*c_iter) > 0) {
+ assert(fiberDim == field._section->getFiberDimension(*c_iter));
+ assert(fiberDim == _section->getFiberDimension(*c_iter));
+ field._section->restrictPoint(*c_iter, &values[0], values.size());
+ _section->updatePointAllAdd(*c_iter, &values[0]);
+ } // if
+ } // for
+ } // if
+
+ return *this;
+} // operator+=
+
+// ----------------------------------------------------------------------
+// Dimensionalize field.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::dimensionalize(void) const
+{ // dimensionalize
+ if (!_metadata.dimsOkay) {
+ std::ostringstream msg;
+ msg << "Cannot dimensionalize field '" << _metadata.label
+ << "' because the flag "
+ << "has been set to keep field nondimensional.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ if (!_section.isNull()) {
+ const chart_type& chart = _section->getChart();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+
+ // Assume fiber dimension is uniform
+ const int fiberDim = (chart.size() > 0) ?
+ _section->getFiberDimension(*chart.begin()) : 0;
+ double_array values(fiberDim);
+
+ spatialdata::units::Nondimensional normalizer;
+
+ for (typename chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter)
+ if (_section->getFiberDimension(*c_iter) > 0) {
+ assert(fiberDim == _section->getFiberDimension(*c_iter));
+
+ _section->restrictPoint(*c_iter, &values[0], values.size());
+ normalizer.dimensionalize(&values[0], values.size(), _metadata.scale);
+ _section->updatePointAll(*c_iter, &values[0]);
+ } // if
+ } // if
+} // dimensionalize
+
+// ----------------------------------------------------------------------
+// Print field to standard out.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::view(const char* label) const
+{ // view
+ std::string vecFieldString;
+ switch(_metadata.vectorFieldType)
+ { // 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 '" << _metadata.vectorFieldType
+ << "'." << std::endl;
+ assert(0);
+ throw std::logic_error("Bad vector field type in MultiField.");
+ } // switch
+
+ std::cout << "Viewing field '" << _metadata.label << "' "<< label << ".\n"
+ << " vector field type: " << vecFieldString << "\n"
+ << " scale: " << _metadata.scale << "\n"
+ << " dimensionalize flag: " << _metadata.dimsOkay << std::endl;
+ if (!_section.isNull())
+ _section->view(label);
+} // view
+
+// ----------------------------------------------------------------------
+// Create PETSc vector scatter for field. This is used to transfer
+// information from the "global" PETSc vector view to the "local"
+// Sieve section view.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::createScatter(const char* context)
+{ // createScatter
+ assert(context);
+ assert(!_section.isNull());
+ assert(!_mesh.sieveMesh().isNull());
+
+ PetscErrorCode err = 0;
+ const bool createScatterOk = true;
+ ScatterInfo& sinfo = _getScatter(context, createScatterOk);
+ if (sinfo.scatter) {
+ assert(sinfo.scatterVec);
+ assert(sinfo.vector);
+ return;
+ } // if
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("GlobalOrder");
+
+ // Get global order (create if necessary).
+ const std::string& orderLabel = _section->getName();
+ const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
+ _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, orderLabel,
+ _section);
+ assert(!order.isNull());
+
+ // Create scatter
+ err = DMMeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order, false, &sinfo.scatter);
+ CHECK_PETSC_ERROR(err);
+
+ // Create scatterVec
+ if (_section->sizeWithBC() > 0) {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, _section->getStorageSize(),
+ _section->restrictSpace(),
+ &sinfo.scatterVec); CHECK_PETSC_ERROR(err);
+ } else {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, 0, 0,
+ &sinfo.scatterVec); CHECK_PETSC_ERROR(err);
+ } // if/else
+
+ // Create vector
+ err = VecCreate(_mesh.comm(), &sinfo.vector);
+ CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)sinfo.vector,
+ _metadata.label.c_str());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(sinfo.vector,
+ order->getLocalSize(), order->getGlobalSize());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
+
+ logger.stagePop();
+} // createScatter
+
+// ----------------------------------------------------------------------
+// Create PETSc vector scatter for field. This is used to transfer
+// information from the "global" PETSc vector view to the "local"
+// Sieve section view. The PETSc vector does not contain constrained
+// DOF. Use createScatterWithBC() to include the constrained DOF in
+// the PETSc vector.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::createScatter(const typename ALE::Obj<typename SieveMesh::numbering_type> numbering,
+ const char* context)
+{ // createScatter
+ assert(!numbering.isNull());
+ assert(context);
+ assert(!_section.isNull());
+ assert(!_mesh.sieveMesh().isNull());
+
+ PetscErrorCode err = 0;
+ const bool createScatterOk = true;
+ ScatterInfo& sinfo = _getScatter(context, createScatterOk);
+
+ // Only create if scatter and scatterVec do not alreay exist.
+ if (sinfo.scatter) {
+ assert(sinfo.scatterVec);
+ assert(sinfo.vector);
+ return;
+ } // if
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("GlobalOrder");
+
+ // Get global order (create if necessary).
+ const std::string& orderLabel =
+ (strlen(context) > 0) ?
+ _section->getName() + std::string("_") + std::string(context) :
+ _section->getName();
+ const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
+ _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order =
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, orderLabel,
+ numbering->getChart().begin(),
+ numbering->getChart().end(),
+ _section);
+ assert(!order.isNull());
+
+ // Create scatter
+ err = DMMeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order, false, &sinfo.scatter);
+ CHECK_PETSC_ERROR(err);
+
+ // Create scatterVec
+ if (_section->sizeWithBC() > 0) {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, _section->getStorageSize(),
+ _section->restrictSpace(),
+ &sinfo.scatterVec); CHECK_PETSC_ERROR(err);
+ } else {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, 0, 0,
+ &sinfo.scatterVec); CHECK_PETSC_ERROR(err);
+ } // if/else
+
+ // Create vector
+ err = VecCreate(_mesh.comm(), &sinfo.vector);
+ CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)sinfo.vector,
+ _metadata.label.c_str());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(sinfo.vector,
+ order->getLocalSize(), order->getGlobalSize());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
+
+#if 0
+ std::cout << "CONTEXT: " << context
+ << ", orderLabel: " << orderLabel
+ << ", section size w/BC: " << _section->sizeWithBC()
+ << ", section size: " << _section->size()
+ << ", global numbering size: " << numbering->getGlobalSize()
+ << ", global size: " << order->getGlobalSize()
+ << std::endl;
+#endif
+
+ logger.stagePop();
+} // createScatter
+
+// ----------------------------------------------------------------------
+// Create PETSc vector scatter for field. This is used to transfer
+// information from the "global" PETSc vector view to the "local"
+// Sieve section view. The PETSc vector does not contain constrained
+// DOF. Use createScatterWithBC() to include the constrained DOF in
+// the PETSc vector.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::createScatterWithBC(const char* context)
+{ // createScatterWithBC
+ assert(context);
+ assert(!_section.isNull());
+ assert(!_mesh.sieveMesh().isNull());
+
+
+ PetscErrorCode err = 0;
+ const bool createScatterOk = true;
+ ScatterInfo& sinfo = _getScatter(context, createScatterOk);
+ if (sinfo.scatter) {
+ assert(sinfo.scatterVec);
+ assert(sinfo.vector);
+ return;
+ } // if
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("GlobalOrder");
+
+ // Get global order (create if necessary).
+ const std::string& orderLabel = _section->getName();
+ const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
+ _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order =
+ sieveMesh->getFactory()->getGlobalOrderWithBC(sieveMesh, orderLabel,
+ _section);
+ assert(!order.isNull());
+
+ // Create scatter
+ err = DMMeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order, true, &sinfo.scatter);
+ CHECK_PETSC_ERROR(err);
+
+ // Create scatterVec
+ if (_section->sizeWithBC() > 0) {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, _section->getStorageSize(),
+ _section->restrictSpace(),
+ &sinfo.scatterVec); CHECK_PETSC_ERROR(err);
+ } else {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, 0, 0,
+ &sinfo.scatterVec); CHECK_PETSC_ERROR(err);
+ } // if/else
+
+ // Create vector
+ err = VecCreate(_mesh.comm(), &sinfo.vector);
+ CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)sinfo.vector,
+ _metadata.label.c_str());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(sinfo.vector,
+ order->getLocalSize(), order->getGlobalSize());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
+
+ logger.stagePop();
+} // createScatterWithBC
+
+// ----------------------------------------------------------------------
+// Create PETSc vector scatter for field. This is used to transfer
+// information from the "global" PETSc vector view to the "local"
+// Sieve section view. The PETSc vector includes constrained DOF. Use
+// createScatter() if constrained DOF should be omitted from the PETSc
+// vector.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::createScatterWithBC(const typename ALE::Obj<typename SieveMesh::numbering_type> numbering,
+ const char* context)
+{ // createScatterWithBC
+ assert(!numbering.isNull());
+ assert(context);
+ assert(!_section.isNull());
+
+ PetscErrorCode err = 0;
+ const bool createScatterOk = true;
+ ScatterInfo& sinfo = _getScatter(context, createScatterOk);
+
+ // Only create if scatter and scatterVec do not alreay exist.
+ if (sinfo.scatter) {
+ assert(sinfo.scatterVec);
+ assert(sinfo.vector);
+ return;
+ } // if
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("GlobalOrder");
+
+ // Get global order (create if necessary).
+ const std::string& orderLabel =
+ (strlen(context) > 0) ?
+ _section->getName() + std::string("_") + std::string(context) :
+ _section->getName();
+ const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
+ _mesh.sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order =
+ sieveMesh->getFactory()->getGlobalOrderWithBC(sieveMesh, orderLabel,
+ numbering->getChart().begin(),
+ numbering->getChart().end(),
+ _section);
+ assert(!order.isNull());
+ //order->view("GLOBAL ORDER"); // DEBUG
+
+ // Create scatter
+ err = DMMeshCreateGlobalScatter(_mesh.sieveMesh(), _section, order, true, &sinfo.scatter);
+ CHECK_PETSC_ERROR(err);
+
+ // Create scatterVec
+ if (_section->sizeWithBC() > 0) {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, _section->getStorageSize(),
+ _section->restrictSpace(),
+ &sinfo.scatterVec); CHECK_PETSC_ERROR(err);
+ } else {
+ err = VecCreateSeqWithArray(PETSC_COMM_SELF, 0, 0,
+ &sinfo.scatterVec); CHECK_PETSC_ERROR(err);
+ } // if/else
+
+ // Create vector
+ err = VecCreate(_mesh.comm(), &sinfo.vector);
+ CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)sinfo.vector,
+ _metadata.label.c_str());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(sinfo.vector,
+ order->getLocalSize(), order->getGlobalSize());
+ CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
+
+#if 0
+ std::cout << "CONTEXT: " << context
+ << ", orderLabel: " << orderLabel
+ << ", section size w/BC: " << _section->sizeWithBC()
+ << ", section size: " << _section->size()
+ << ", section storage size: " << _section->getStorageSize()
+ << ", global numbering size: " << numbering->getGlobalSize()
+ << ", global size: " << order->getGlobalSize()
+ << ", scatter from size: " << sinfo.scatter->from_n
+ << std::endl;
+#endif
+
+ logger.stagePop();
+} // createScatterWithBC
+
+// ----------------------------------------------------------------------
+// Get PETSc vector associated with field.
+template<typename mesh_type, typename section_type>
+PetscVec
+pylith::topology::MultiField<mesh_type, section_type>::vector(const char* context)
+{ // vector
+ ScatterInfo& sinfo = _getScatter(context);
+ return sinfo.vector;
+} // vector
+
+// ----------------------------------------------------------------------
+// Get PETSc vector associated with field.
+template<typename mesh_type, typename section_type>
+const PetscVec
+pylith::topology::MultiField<mesh_type, section_type>::vector(const char* context) const
+{ // vector
+ const ScatterInfo& sinfo = _getScatter(context);
+ return sinfo.vector;
+} // vector
+
+// ----------------------------------------------------------------------
+// Scatter section information across processors to update the
+// PETSc vector view of the field.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::scatterSectionToVector(const char* context) const
+{ // scatterSectionToVector
+ assert(context);
+
+ const ScatterInfo& sinfo = _getScatter(context);
+ scatterSectionToVector(sinfo.vector, context);
+} // scatterSectionToVector
+
+// ----------------------------------------------------------------------
+// Scatter section information across processors to update the
+// PETSc vector view of the field.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::scatterSectionToVector(const PetscVec vector,
+ const char* context) const
+{ // scatterSectionToVector
+ assert(vector);
+ assert(context);
+ assert(!_section.isNull());
+
+ PetscErrorCode err = 0;
+ const ScatterInfo& sinfo = _getScatter(context);
+ err = VecScatterBegin(sinfo.scatter, sinfo.scatterVec, vector,
+ INSERT_VALUES, SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
+ err = VecScatterEnd(sinfo.scatter, sinfo.scatterVec, vector,
+ INSERT_VALUES, SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
+} // scatterSectionToVector
+
+// ----------------------------------------------------------------------
+// Scatter PETSc vector information across processors to update the
+// section view of the field.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::scatterVectorToSection(const char* context) const
+{ // scatterVectorToSection
+ assert(context);
+
+ const ScatterInfo& sinfo = _getScatter(context);
+ scatterVectorToSection(sinfo.vector, context);
+} // scatterVectorToSection
+
+// ----------------------------------------------------------------------
+// Scatter PETSc vector information across processors to update the
+// section view of the field.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::scatterVectorToSection(const PetscVec vector,
+ const char* context) const
+{ // scatterVectorToSection
+ assert(vector);
+ assert(context);
+ assert(!_section.isNull());
+
+ PetscErrorCode err = 0;
+ const ScatterInfo& sinfo = _getScatter(context);
+ err = VecScatterBegin(sinfo.scatter, vector, sinfo.scatterVec,
+ INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
+ err = VecScatterEnd(sinfo.scatter, vector, sinfo.scatterVec,
+ INSERT_VALUES, SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
+} // scatterVectorToSection
+
+// ----------------------------------------------------------------------
+// Setup split field with all one space per spatial dimension.
+template<typename mesh_type, typename section_type>
+void
+pylith::topology::MultiField<mesh_type, section_type>::splitDefault(void)
+{ // splitDefault
+ assert(!_section.isNull());
+ const int spaceDim = _mesh.dimension();
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _section->addSpace(); // displacements
+
+ const chart_type& chart = _section->getChart();
+
+ const typename chart_type::const_iterator chartBegin = chart.begin();
+ const typename chart_type::const_iterator chartEnd = chart.end();
+ for (int fibration=0; fibration < spaceDim; ++fibration)
+ for (typename chart_type::const_iterator c_iter = chart.begin();
+ c_iter != chartEnd;
+ ++c_iter) {
+ assert(spaceDim == _section->getFiberDimension(*c_iter));
+ _section->setFiberDimension(*c_iter, 1, fibration);
+ } // for
+} // splitDefault
+
+// ----------------------------------------------------------------------
+// Get scatter for given context.
+template<typename mesh_type, typename section_type>
+typename pylith::topology::MultiField<mesh_type, section_type>::ScatterInfo&
+pylith::topology::MultiField<mesh_type, section_type>::_getScatter(const char* context,
+ const bool createOk)
+{ // _getScatter
+ assert(context);
+
+ const bool isNewScatter = _scatters.find(context) == _scatters.end();
+
+ if (isNewScatter && !createOk) {
+ std::ostringstream msg;
+ msg << "Scatter for context '" << context << "' does not exist.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ ScatterInfo& sinfo = _scatters[context];
+ if (isNewScatter) {
+ sinfo.vector = 0;
+ sinfo.scatter = 0;
+ sinfo.scatterVec = 0;
+ } // if
+ assert(_scatters.find(context) != _scatters.end());
+
+ return sinfo;
+} // _getScatter
+
+// ----------------------------------------------------------------------
+// Get scatter for given context.
+template<typename mesh_type, typename section_type>
+const typename pylith::topology::MultiField<mesh_type, section_type>::ScatterInfo&
+pylith::topology::MultiField<mesh_type, section_type>::_getScatter(const char* context) const
+{ // _getScatter
+ assert(context);
+
+ const typename scatter_map_type::const_iterator s_iter =
+ _scatters.find(context);
+ if (s_iter == _scatters.end()) {
+ std::ostringstream msg;
+ msg << "Scatter for context '" << context << "' does not exist.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ return s_iter->second;
+} // _getScatter
+
+
+// End of file
Added: short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.hh
===================================================================
--- short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.hh (rev 0)
+++ short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.hh 2011-05-21 19:47:16 UTC (rev 18415)
@@ -0,0 +1,409 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/MultiField.hh
+ *
+ * @brief Vector field over the vertices or cells of a finite-element
+ * mesh.
+ */
+
+#if !defined(pylith_topology_multifield_hh)
+#define pylith_topology_multifield_hh
+
+// Include directives ---------------------------------------------------
+#include "FieldBase.hh" // ISA FieldBase
+
+#include "pylith/utils/arrayfwd.hh" // HASA int_array
+#include "pylith/utils/petscfwd.h" // HASA PetscVec
+
+#include <petscdmmesh.hh>
+
+#include <map> // USES std::map
+#include <string> // USES std::string
+
+// MultiField ----------------------------------------------------------------
+/** @brief Vector field over the vertices or cells of a finite-element
+ * mesh.
+ *
+ * Extends Sieve real general section by adding metadata.
+ */
+template<typename mesh_type,
+ typename section_type>
+class pylith::topology::MultiField : public FieldBase
+{ // MultiField
+ friend class TestMultiFieldMesh; // unit testing
+ friend class TestMultiFieldSubMesh; // unit testing
+
+// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
+public:
+
+ // Convenience typedefs
+ typedef mesh_type Mesh;
+
+ typedef ALE::ISieveVisitor::RestrictVisitor<section_type> RestrictVisitor;
+ typedef ALE::ISieveVisitor::UpdateAddVisitor<section_type> UpdateAddVisitor;
+ typedef ALE::ISieveVisitor::UpdateAllVisitor<section_type> UpdateAllVisitor;
+
+// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
+private:
+
+ // Convenience typedefs
+ typedef typename mesh_type::SieveMesh SieveMesh;
+ typedef typename SieveMesh::label_sequence label_sequence;
+ typedef typename section_type::chart_type chart_type;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Default constructor.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ MultiField(const mesh_type& mesh);
+
+ /** Constructor with mesh, section, and metadata.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ MultiField(const mesh_type& mesh,
+ const ALE::Obj<section_type>& section,
+ const Metadata& metadata);
+
+ /// Destructor.
+ ~MultiField(void);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(void);
+
+ /** Get Sieve section.
+ *
+ * @returns Sieve section.
+ */
+ const ALE::Obj<section_type>& section(void) const;
+
+ /** Get mesh associated with field.
+ *
+ * @returns Finite-element mesh.
+ */
+ const mesh_type& mesh(void) const;
+
+ /** Set label for field.
+ *
+ * @param value Label for field.
+ */
+ void label(const char* value);
+
+ /** Get label for field.
+ *
+ * @returns Label for field.
+ */
+ const char* label(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;
+
+ /** Get spatial dimension of domain.
+ *
+ * @returns Spatial dimension of domain.
+ */
+ int spaceDim(void) const;
+
+ /** Get the number of sieve points in the chart.
+ *
+ * @returns the chart size.
+ */
+ int chartSize(void) const;
+
+ /** Get the number of degrees of freedom.
+ *
+ * @returns the number of degrees of freedom.
+ */
+ int sectionSize(void) const;
+
+ /// Create sieve section.
+ void newSection(void);
+
+ /** Create sieve section and set chart and fiber dimesion for
+ * sequence of points.
+ *
+ * @param points Points over which to define section.
+ * @param dim Fiber dimension for section.
+ */
+ void newSection(const ALE::Obj<label_sequence>& points,
+ const int fiberDim);
+
+ /** Create sieve section and set chart and fiber dimesion for a list
+ * of points.
+ *
+ * @param points Points over which to define section.
+ * @param dim Fiber dimension for section.
+ */
+ void newSection(const int_array& 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.
+ * @param stratum Stratum depth (for vertices) and height (for cells).
+ */
+ void newSection(const DomainEnum domain,
+ const int fiberDim,
+ const int stratum =0);
+
+ /** Create section using src field as template with given fiber dimension.
+ *
+ * @param sec MultiField defining layout.
+ * @param fiberDim Fiber dimension.
+ */
+ void newSection(const MultiField& src,
+ 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 MultiField defining layout.
+ */
+ void cloneSection(const MultiField& src);
+
+ /// Clear variables associated with section.
+ void clear(void);
+
+ /// Allocate field.
+ void allocate(void);
+
+ /// Zero section values (does not zero constrained values).
+ void zero(void);
+
+ /// Zero section values (including constrained values).
+ void zeroAll(void);
+
+ /// Complete section by assembling across processors.
+ void complete(void);
+
+ /** Copy field values and metadata.
+ *
+ * @param field MultiField to copy.
+ */
+ void copy(const MultiField& field);
+
+ /** Copy field values.
+ *
+ * @param field MultiField to copy.
+ */
+ void copy(const ALE::Obj<section_type>& field);
+
+ /** Add two fields, storing the result in one of the fields.
+ *
+ * @param field MultiField to add.
+ */
+ MultiField& operator+=(const MultiField& field);
+
+ /** Dimensionalize field. Throws runtime_error if field is not
+ * allowed to be dimensionalized.
+ */
+ void dimensionalize(void) const;
+
+ /** Print field to standard out.
+ *
+ * @param label Label for output.
+ */
+ void view(const char* label) const;
+
+ /** Create PETSc vector scatter for field. This is used to transfer
+ * information from the "global" PETSc vector view to the "local"
+ * Sieve section view. The PETSc vector does not contain constrained
+ * DOF. Use createScatterWithBC() to include the constrained DOF in
+ * the PETSc vector.
+ *
+ * @param context Label for context associated with vector.
+ */
+ void createScatter(const char* context ="");
+
+
+ /** Create PETSc vector scatter for field. This is used to transfer
+ * information from the "global" PETSc vector view to the "local"
+ * Sieve section view. The PETSc vector does not contain constrained
+ * DOF. Use createScatterWithBC() to include the constrained DOF in
+ * the PETSc vector.
+ *
+ * @param numbering Numbering used to select points in section.
+ * @param context Label for context associated with vector.
+ */
+ void createScatter(const typename ALE::Obj<typename SieveMesh::numbering_type> numbering,
+ const char* context ="");
+
+ /** Create PETSc vector scatter for field. This is used to transfer
+ * information from the "global" PETSc vector view to the "local"
+ * Sieve section view. The PETSc vector includes constrained
+ * DOF. Use createScatter() if constrained DOF should be omitted
+ * from the PETSc vector.
+ *
+ * @param context Label for context associated with vector.
+ */
+ void createScatterWithBC(const char* context ="");
+
+
+ /** Create PETSc vector scatter for field. This is used to transfer
+ * information from the "global" PETSc vector view to the "local"
+ * Sieve section view. The PETSc vector includes constrained
+ * DOF. Use createScatter() if constrained DOF should be omitted
+ * from the PETSc vector.
+ *
+ * @param numbering Numbering used to select points in section.
+ * @param context Label for context associated with vector.
+ */
+ void createScatterWithBC(const typename ALE::Obj<typename SieveMesh::numbering_type> numbering,
+ const char* context ="");
+
+ /** Get PETSc vector associated with field.
+ *
+ * @param context Label for context associated with vector.
+ * @returns PETSc vector.
+ */
+ PetscVec vector(const char* context ="");
+
+ /** Get PETSc vector associated with field.
+ *
+ * @param context Label for context associated with vector.
+ * @returns PETSc vector.
+ */
+ const PetscVec vector(const char* context ="") const;
+
+ /// Scatter section information across processors to update the
+ /// PETSc vector view of the field.
+ void scatterSectionToVector(const char* context ="") const;
+
+ /** Scatter section information across processors to update the
+ * PETSc vector view of the field.
+ *
+ * @param vector PETSc vector to update.
+ */
+ void scatterSectionToVector(const PetscVec vector,
+ const char* context ="") const;
+
+ /// Scatter PETSc vector information across processors to update the
+ /// Sieve section view of the field.
+ void scatterVectorToSection(const char* context ="") const;
+
+ /** Scatter section information across processors to update the
+ * PETSc vector view of the field.
+ *
+ * @param vector PETSc vector used in update.
+ */
+ void scatterVectorToSection(const PetscVec vector,
+ const char* context ="") const;
+
+ /// Setup split field with all entries set to a default space of 0.
+ void splitDefault(void);
+
+// PRIVATE STRUCTS //////////////////////////////////////////////////////
+private :
+
+ /// Data structures used in scattering to/from PETSc Vecs.
+ struct ScatterInfo {
+ PetscVec vector; ///< PETSc vector associated with field.
+ PetscVecScatter scatter; ///< PETSc scatter associated with field.
+ PetscVec scatterVec; ///< PETSC vector used in scattering.
+ }; // ScatterInfo
+
+// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
+private :
+
+ typedef std::map<std::string, ScatterInfo> scatter_map_type;
+
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /** Get scatter for given context.
+ *
+ * @param context Context for scatter.
+ * @param createOk If true, okay to create new scatter for
+ * context, if false will throw std::runtime_error if scatter for
+ * context doesn't already exist.
+ */
+ ScatterInfo& _getScatter(const char* context,
+ const bool createOk =false);
+
+ /** Get scatter for given context.
+ *
+ * @param context Context for scatter.
+ */
+ const ScatterInfo& _getScatter(const char* context) const;
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ Metadata _metadata;
+ const mesh_type& _mesh; ///< Mesh associated with section.
+ ALE::Obj<section_type> _section; ///< Real section with data.
+ scatter_map_type _scatters; ///< Collection of scatters.
+
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ MultiField(const MultiField&); ///< Not implemented
+ const MultiField& operator=(const MultiField&); ///< Not implemented
+
+}; // MultiField
+
+#include "MultiField.icc"
+#include "MultiField.cc"
+
+#endif // pylith_topology_multifield_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.icc
===================================================================
--- short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.icc (rev 0)
+++ short/3D/PyLith/branches/multifields/libsrc/pylith/topology/MultiField.icc 2011-05-21 19:47:16 UTC (rev 18415)
@@ -0,0 +1,101 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+#if !defined(pylith_topology_multifield_hh)
+#error "MultiField.icc must be included only from MultiField.hh"
+#else
+
+#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
+
+// Get Sieve section.
+template<typename mesh_type, typename section_type>
+inline
+const ALE::Obj<section_type>&
+pylith::topology::MultiField<mesh_type, section_type>::section(void) const {
+ return _section;
+}
+
+// Get mesh associated with field.
+template<typename mesh_type, typename section_type>
+inline
+const
+mesh_type&
+pylith::topology::MultiField<mesh_type, section_type>::mesh(void) const {
+ return _mesh;
+}
+
+// Get label for field.
+template<typename mesh_type, typename section_type>
+inline
+const char*
+pylith::topology::MultiField<mesh_type, section_type>::label(void) const {
+ return _metadata.label.c_str();
+}
+
+// Set vector field type
+template<typename mesh_type, typename section_type>
+inline
+void
+pylith::topology::MultiField<mesh_type, section_type>::vectorFieldType(const VectorFieldEnum value) {
+ _metadata.vectorFieldType = value;
+}
+
+// Get vector field type
+template<typename mesh_type, typename section_type>
+inline
+typename pylith::topology::MultiField<mesh_type, section_type>::VectorFieldEnum
+pylith::topology::MultiField<mesh_type, section_type>::vectorFieldType(void) const {
+ return _metadata.vectorFieldType;
+}
+
+// Set scale for dimensionalizing field.
+template<typename mesh_type, typename section_type>
+inline
+void
+pylith::topology::MultiField<mesh_type, section_type>::scale(const double value) {
+ _metadata.scale = value;
+}
+
+// Get scale for dimensionalizing field.
+template<typename mesh_type, typename section_type>
+inline
+double
+pylith::topology::MultiField<mesh_type, section_type>::scale(void) const {
+ return _metadata.scale;
+}
+
+// Set flag indicating whether it is okay to dimensionalize field.
+template<typename mesh_type, typename section_type>
+inline
+void
+pylith::topology::MultiField<mesh_type, section_type>::addDimensionOkay(const bool value) {
+ _metadata.dimsOkay = value;
+}
+
+// Set flag indicating whether it is okay to dimensionalize field.
+template<typename mesh_type, typename section_type>
+inline
+bool
+pylith::topology::MultiField<mesh_type, section_type>::addDimensionOkay(void) const {
+ return _metadata.dimsOkay;
+}
+
+#endif
+
+
+// End of file
Modified: short/3D/PyLith/branches/multifields/libsrc/pylith/topology/topologyfwd.hh
===================================================================
--- short/3D/PyLith/branches/multifields/libsrc/pylith/topology/topologyfwd.hh 2011-05-21 19:28:25 UTC (rev 18414)
+++ short/3D/PyLith/branches/multifields/libsrc/pylith/topology/topologyfwd.hh 2011-05-21 19:47:16 UTC (rev 18415)
@@ -58,6 +58,8 @@
class FieldBase;
template<typename mesh_type,
typename section_type =ALE::IGeneralSection<pylith::SieveMesh::point_type, double> > class Field;
+ template<typename mesh_type,
+ typename section_type =ALE::IGeneralSection<pylith::SieveMesh::point_type, double> > class MultiField;
template<typename field_type> class Fields;
template<typename mesh_type> class PackedFields;
Modified: short/3D/PyLith/branches/multifields/modulesrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/branches/multifields/modulesrc/topology/Makefile.am 2011-05-21 19:28:25 UTC (rev 18414)
+++ short/3D/PyLith/branches/multifields/modulesrc/topology/Makefile.am 2011-05-21 19:47:16 UTC (rev 18415)
@@ -30,6 +30,7 @@
MeshOps.i \
FieldBase.i \
Field.i \
+ MultiField.i \
Fields.i \
PackedFields.i \
SolutionFields.i \
Added: short/3D/PyLith/branches/multifields/modulesrc/topology/MultiField.i
===================================================================
--- short/3D/PyLith/branches/multifields/modulesrc/topology/MultiField.i (rev 0)
+++ short/3D/PyLith/branches/multifields/modulesrc/topology/MultiField.i 2011-05-21 19:47:16 UTC (rev 18415)
@@ -0,0 +1,257 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file modulesrc/topology/MultiField.hh
+ *
+ * @brief Python interface to C++ MultiField object.
+ */
+
+// Typemaps for returning reference in operator+=. The default
+// behavior is that the Python object will gain ownership. We want the
+// C++ object to retain ownership. This could be alleviated by using a
+// shared pointer.
+%typemap(out) pylith::topology::MultiField<pylith::topology::Mesh>& {
+ $result = SWIG_NewPointerObj(SWIG_as_voidptr($1), $descriptor(pylith::topology::MultiField<pylith::topology::Mesh>*), 0 | 0);
+ return $result;
+ }
+%typemap(out) pylith::topology::MultiField<pylith::topology::SubMesh>& {
+ $result = SWIG_NewPointerObj(SWIG_as_voidptr($1), $descriptor(pylith::topology::MultiField<pylith::topology::SubMesh>*), 0 | 0);
+ return $result;
+ }
+
+
+namespace pylith {
+ namespace topology {
+
+ template<typename mesh_type>
+ class MultiField : public FieldBase
+ { // MultiField
+
+ // PRIVATE TYPEDEFS ///////////////////////////////////////////////
+ private:
+
+ // Convenience typedefs
+ typedef typename mesh_type::RealSection RealSection;
+ typedef typename mesh_type::SieveMesh SieveMesh;
+ typedef typename SieveMesh::label_sequence label_sequence;
+ typedef typename RealSection::chart_type chart_type;
+
+ // PUBLIC MEMBERS /////////////////////////////////////////////////
+ public :
+
+ /** Default constructor.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ MultiField(const mesh_type& mesh);
+
+ /// Destructor.
+ ~MultiField(void);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(void);
+
+ /** Get mesh associated with field.
+ *
+ * @returns Finite-element mesh.
+ */
+ const mesh_type& mesh(void) const;
+
+ /** Set label for field.
+ *
+ * @param value Label for field.
+ */
+ void label(const char* value);
+
+ /** Get label for field.
+ *
+ * @returns Label for field.
+ */
+ const char* label(void) const;
+
+ /** Set vector field type
+ *
+ * @param value Type of vector field.
+ */
+ void vectorFieldType(const pylith::topology::FieldBase::VectorFieldEnum value);
+
+ /** Get vector field type
+ *
+ * @returns Type of vector field.
+ */
+ pylith::topology::FieldBase::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;
+
+ /** Get spatial dimension of domain.
+ *
+ * @returns Spatial dimension of domain.
+ */
+ int spaceDim(void) const;
+
+ /** Get the number of sieve points in the chart.
+ *
+ * @returns the chart size.
+ */
+ int chartSize(void) const;
+
+ /** Get the number of degrees of freedom.
+ *
+ * @returns the number of degrees of freedom.
+ */
+ int sectionSize(void) const;
+
+ /// Create sieve section.
+ void newSection(void);
+
+ /** 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.
+ * @param stratum Stratum depth (for vertices) and height (for cells).
+ */
+ void newSection(const pylith::topology::FieldBase::DomainEnum domain,
+ const int fiberDim,
+ const int stratum =0);
+
+ /** 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 cloneSection(const MultiField& src);
+
+ /// Clear variables associated with section.
+ void clear(void);
+
+ /// Allocate field.
+ void allocate(void);
+
+ /// Zero section values (excluding constrained DOF).
+ void zero(void);
+
+ /// Zero section values (including constrained DOF).
+ void zeroAll(void);
+
+ /// Complete section by assembling across processors.
+ void complete(void);
+
+ /** Copy field values and metadata.
+ *
+ * @param field MultiField to copy.
+ */
+ void copy(const MultiField& field);
+
+ /** Add two fields, storing the result in one of the fields.
+ *
+ * @param field MultiField to add.
+ */
+ MultiField& operator+=(const MultiField& 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);
+
+ /** Create PETSc vector scatter for field. This is used to transfer
+ * information from the "global" PETSc vector view to the "local"
+ * Sieve section view.
+ *
+ * @param context Label for context associated with vector.
+ */
+ void createScatter(const char* context ="");
+
+ /** Get PETSc vector associated with field.
+ *
+ * @param context Label for context associated with vector.
+ * @returns PETSc vector.
+ */
+ PetscVec vector(const char* context ="");
+
+ /** Get PETSc vector associated with field.
+ *
+ * @param context Label for context associated with vector.
+ * @returns PETSc vector.
+ */
+ const PetscVec vector(const char* context ="") const;
+
+ /// Scatter section information across processors to update the
+ /// PETSc vector view of the field.
+ void scatterSectionToVector(const char* context ="") const;
+
+ /** Scatter section information across processors to update the
+ * PETSc vector view of the field.
+ *
+ * @param vector PETSc vector to update.
+ */
+ void scatterSectionToVector(const PetscVec vector,
+ const char* context ="") const;
+
+ /// Scatter PETSc vector information across processors to update the
+ /// Sieve section view of the field.
+ void scatterVectorToSection(const char* context ="") const;
+
+ /** Scatter section information across processors to update the
+ * PETSc vector view of the field.
+ *
+ * @param vector PETSc vector used in update.
+ */
+ void scatterVectorToSection(const PetscVec vector,
+ const char* context ="") const;
+
+ /// Setup split field with all entries set to a default space of 0.
+ void splitDefault(void);
+
+ }; // MultiField
+
+ } // topology
+} // pylith
+
+
+// End of file
Modified: short/3D/PyLith/branches/multifields/modulesrc/topology/topology.i
===================================================================
--- short/3D/PyLith/branches/multifields/modulesrc/topology/topology.i 2011-05-21 19:28:25 UTC (rev 18414)
+++ short/3D/PyLith/branches/multifields/modulesrc/topology/topology.i 2011-05-21 19:47:16 UTC (rev 18415)
@@ -27,6 +27,7 @@
#include "pylith/topology/FieldBase.hh"
#include "pylith/topology/Field.hh"
#include "pylith/topology/Fields.hh"
+#include "pylith/topology/MultiField.hh"
#include "pylith/topology/PackedFields.hh"
#include "pylith/topology/SolutionFields.hh"
#include "pylith/topology/Jacobian.hh"
@@ -65,6 +66,7 @@
%include "MeshOps.i"
%include "FieldBase.i"
%include "Field.i"
+%include "MultiField.i"
%include "Fields.i"
%include "PackedFields.i"
%include "SolutionFields.i"
@@ -74,11 +76,20 @@
%include "ReverseCuthillMcKee.i"
// Template instatiation
+
+ // Field
%template(MeshField) pylith::topology::Field<pylith::topology::Mesh>;
%template(SubMeshField) pylith::topology::Field<pylith::topology::SubMesh>;
+
+// MultiField
+%template(MeshMultiField) pylith::topology::MultiField<pylith::topology::Mesh>;
+%template(SubMeshMultiField) pylith::topology::MultiField<pylith::topology::SubMesh>;
+
+// Fields
%template(MeshFields) pylith::topology::Fields<pylith::topology::Field<pylith::topology::Mesh> >;
%template(SubMeshFields) pylith::topology::Fields<pylith::topology::Field<pylith::topology::SubMesh> >;
+// PackedFields
%template(MeshPackedFields) pylith::topology::PackedFields<pylith::topology::Mesh>;
%template(SubMeshPackedFields) pylith::topology::PackedFields<pylith::topology::SubMesh>;
Modified: short/3D/PyLith/branches/multifields/unittests/libtests/topology/Makefile.am
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/topology/Makefile.am 2011-05-21 19:28:25 UTC (rev 18414)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/topology/Makefile.am 2011-05-21 19:47:16 UTC (rev 18415)
@@ -34,6 +34,8 @@
TestFieldSubMesh.cc \
TestFieldsMesh.cc \
TestFieldsSubMesh.cc \
+ TestMultiFieldMesh.cc \
+ TestMultiFieldSubMesh.cc \
TestPackedFieldsMesh.cc \
TestSolutionFields.cc \
TestJacobian.cc \
@@ -49,6 +51,8 @@
TestFieldSubMesh.hh \
TestFieldsMesh.hh \
TestFieldsSubMesh.hh \
+ TestMultiFieldMesh.hh \
+ TestMultiFieldSubMesh.hh \
TestPackedFieldsMesh.hh \
TestSolutionFields.hh \
TestJacobian.hh \
Added: short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldMesh.cc
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldMesh.cc (rev 0)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldMesh.cc 2011-05-21 19:47:16 UTC (rev 18415)
@@ -0,0 +1,1352 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestMultiFieldMesh.hh" // Implementation of class methods
+
+#include "pylith/topology/MultiField.hh" // USES MultiField
+#include "pylith/topology/Mesh.hh" // USES Mesh
+
+#include "pylith/utils/array.hh" // USES double_array
+
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::topology::TestMultiFieldMesh );
+
+// ----------------------------------------------------------------------
+namespace pylith {
+ namespace topology {
+ namespace _TestMultiFieldMesh {
+ const int cellDim = 2;
+ const int nvertices = 4;
+ const int ncells = 1;
+ const int ncorners = 4;
+ const int cells[] = { 0, 1, 2, 3 };
+ const double coordinates[] = {
+ 0.0, 0.0,
+ 1.0, 0.0,
+ 0.0, 1.0,
+ 1.0, 1.0,
+ };
+ } // _TestMultiFieldMesh
+ } // topology
+} // pylith
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::topology::TestMultiFieldMesh::testConstructor(void)
+{ // testConstructor
+ Mesh mesh;
+ MultiField<Mesh> field(mesh);
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test section().
+void
+pylith::topology::TestMultiFieldMesh::testSection(void)
+{ // testSection
+ Mesh mesh;
+ MultiField<Mesh> field(mesh);
+
+ mesh.createSieveMesh();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(section.isNull());
+} // testSection
+
+// ----------------------------------------------------------------------
+// Test mesh().
+void
+pylith::topology::TestMultiFieldMesh::testMesh(void)
+{ // testMesh
+ Mesh mesh;
+ _buildMesh(&mesh);
+ MultiField<Mesh> field(mesh);
+
+ const Mesh& mesh2 = field.mesh();
+ CPPUNIT_ASSERT_EQUAL(_TestMultiFieldMesh::cellDim, mesh2.dimension());
+} // testMesh
+
+// ----------------------------------------------------------------------
+// Test label().
+void
+pylith::topology::TestMultiFieldMesh::testLabel(void)
+{ // testLabel
+ const std::string label = "velocity";
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ MultiField<Mesh> field(mesh);
+
+ field.label(label.c_str());
+ CPPUNIT_ASSERT_EQUAL(label, std::string(field.label()));
+} // testLabel
+
+// ----------------------------------------------------------------------
+// Test vectorFieldType().
+void
+pylith::topology::TestMultiFieldMesh::testVectorFieldType(void)
+{ // testVectorFieldType
+} // testVectorFieldType
+
+// ----------------------------------------------------------------------
+// Test scale().
+void
+pylith::topology::TestMultiFieldMesh::testScale(void)
+{ // testScale
+} // testScale
+
+// ----------------------------------------------------------------------
+// Test addDimensionsOkay().
+void
+pylith::topology::TestMultiFieldMesh::testAddDimensionsOkay(void)
+{ // testAddDimensionsOkay
+} // testAddDimensionsOkay
+
+// ----------------------------------------------------------------------
+// Test spaceDim().
+void
+pylith::topology::TestMultiFieldMesh::testSpaceDim(void)
+{ // testSpaceDim
+ Mesh mesh;
+ _buildMesh(&mesh);
+ MultiField<Mesh> field(mesh);
+
+ CPPUNIT_ASSERT_EQUAL(_TestMultiFieldMesh::cellDim, field.spaceDim());
+} // testSpaceDim
+
+// ----------------------------------------------------------------------
+// Test newSection().
+void
+pylith::topology::TestMultiFieldMesh::testNewSection(void)
+{ // testNewSection
+ Mesh mesh;
+ _buildMesh(&mesh);
+
+ MultiField<Mesh> field(mesh);
+ const std::string& label = "field A";
+ field.label(label.c_str());
+
+ field.newSection();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+} // testNewSection
+
+// ----------------------------------------------------------------------
+// Test newSection(points).
+void
+pylith::topology::TestMultiFieldMesh::testNewSectionPoints(void)
+{ // testNewSectionPoints
+ const int fiberDim = 2;
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ MultiField<Mesh> field(mesh);
+ const std::string& label = "field A";
+ field.label(label.c_str());
+
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ field.newSection(vertices, fiberDim);
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ CPPUNIT_ASSERT(!vertices.isNull());
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
+
+ CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+} // testNewSectionPoints
+
+// ----------------------------------------------------------------------
+// Test newSection(int_array).
+void
+pylith::topology::TestMultiFieldMesh::testNewSectionPointsArray(void)
+{ // testNewSectionPointsArray
+ const int fiberDim = 2;
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ MultiField<Mesh> field(mesh);
+ const std::string& label = "field A";
+ field.label(label.c_str());
+
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const int npts = vertices->size() / 2;
+ int_array pointsIn(npts);
+ int_array pointsOut(vertices->size() - npts);
+ int count = 0;
+ size_t iIn = 0;
+ size_t iOut = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ if (count % 2 == 0)
+ pointsIn[iIn++] = *v_iter;
+ else
+ pointsOut[iOut++] = *v_iter;
+ ++count;
+ } // for
+ CPPUNIT_ASSERT_EQUAL(iIn, pointsIn.size());
+ CPPUNIT_ASSERT_EQUAL(iOut, pointsOut.size());
+
+ field.newSection(pointsIn, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ // Points in array should have a fiber dimension of fiberDim.
+ for (int i=0; i < pointsIn.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(pointsIn[i]));
+
+ // Points not int array should have a fiber dimension of zero.
+ for (int i=0; i < pointsOut.size(); ++i)
+ CPPUNIT_ASSERT_EQUAL(0, section->getFiberDimension(pointsOut[i]));
+
+ CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+} // testNewSectionPointsArray
+
+// ----------------------------------------------------------------------
+// Test newSection(domain).
+void
+pylith::topology::TestMultiFieldMesh::testNewSectionDomain(void)
+{ // testNewSectionDomain
+ const int fiberDim = 2;
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ MultiField<Mesh> field(mesh);
+ const std::string& label = "field A";
+ field.label(label.c_str());
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
+
+ CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+} // testNewSectionDomain
+
+// ----------------------------------------------------------------------
+// Test newSection(field).
+void
+pylith::topology::TestMultiFieldMesh::testNewSectionField(void)
+{ // testNewSectionField
+ const int fiberDim = 3;
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ // Create field with atlas to use to create new field
+ MultiField<Mesh> fieldSrc(mesh);
+ fieldSrc.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ const ALE::Obj<Mesh::RealSection>& sectionSrc = fieldSrc.section();
+ CPPUNIT_ASSERT(!sectionSrc.isNull());
+ const Mesh::RealSection::chart_type& chart = sectionSrc->getChart();
+
+ const int fiberDim2 = 5;
+ MultiField<Mesh> field(mesh);
+ const std::string& label = "field A";
+ field.label(label.c_str());
+ field.newSection(fieldSrc, fiberDim2);
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ CPPUNIT_ASSERT_EQUAL(fiberDim2, section->getFiberDimension(*v_iter));
+
+ CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+} // testNewSectionField
+
+// ----------------------------------------------------------------------
+// Test cloneSection().
+void
+pylith::topology::TestMultiFieldMesh::testCloneSection(void)
+{ // testCloneSection
+ const int fiberDim = 3;
+ const int nconstraints[] = { 0, 2, 1, 3 };
+ const int constraints[] = {
+ // 0
+ 0, 2, // 1
+ 2, // 2
+ 0, 1, 2, // 3
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ // Create field with atlas to use to create new field
+ MultiField<Mesh> fieldSrc(mesh);
+ { // Setup source field
+ fieldSrc.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ int iV=0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ section->addConstraintDimension(*v_iter, nconstraints[iV++]);
+ fieldSrc.allocate();
+
+ int index = 0;
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, index += nconstraints[i++])
+ section->setConstraintDof(*v_iter, &constraints[index]);
+ fieldSrc.zero();
+ fieldSrc.createScatter();
+ fieldSrc.createScatter("A");
+ } // Setup source field
+
+ MultiField<Mesh> field(mesh);
+ const std::string& label = "field A";
+ field.label(label.c_str());
+ field.cloneSection(fieldSrc);
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ int iV = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
+ CPPUNIT_ASSERT_EQUAL(nconstraints[iV++],
+ section->getConstraintDimension(*v_iter));
+ } // for
+
+ // Verify vector scatters were also copied.
+ CPPUNIT_ASSERT_EQUAL(fieldSrc._scatters[""].scatter,
+ field._scatters[""].scatter);
+ CPPUNIT_ASSERT_EQUAL(fieldSrc._scatters["A"].scatter,
+ field._scatters["A"].scatter);
+
+ CPPUNIT_ASSERT_EQUAL(label, std::string(section->getName()));
+} // testCloneSection
+
+// ----------------------------------------------------------------------
+// Test clear().
+void
+pylith::topology::TestMultiFieldMesh::testClear(void)
+{ // testClear
+ Mesh mesh(_TestMultiFieldMesh::cellDim);
+ MultiField<Mesh> field(mesh);
+
+ field.scale(2.0);
+ field.vectorFieldType(MultiField<Mesh>::TENSOR);
+ field.addDimensionOkay(true);
+
+ field.clear();
+
+ CPPUNIT_ASSERT_EQUAL(1.0, field._metadata.scale);
+ CPPUNIT_ASSERT_EQUAL(MultiField<Mesh>::OTHER, field._metadata.vectorFieldType);
+ CPPUNIT_ASSERT_EQUAL(false, field._metadata.dimsOkay);
+} // testClear
+
+// ----------------------------------------------------------------------
+// Test allocate().
+void
+pylith::topology::TestMultiFieldMesh::testAllocate(void)
+{ // testAllocate
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ const double tolerance = 1.0e-6;
+ i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ } // for
+ } // for
+} // testAllocate
+
+// ----------------------------------------------------------------------
+// Test zero().
+void
+pylith::topology::TestMultiFieldMesh::testZero(void)
+{ // testZero
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.zero();
+
+ const double tolerance = 1.0e-6;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[iDim], tolerance);
+ } // for
+ } // for
+} // testZero
+
+// ----------------------------------------------------------------------
+// Test zero().
+void
+pylith::topology::TestMultiFieldMesh::testZeroAll(void)
+{ // testZeroAll
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+ const int nconstraints[] = { 0, 2, 1, 3 };
+ const int constraints[] = {
+ // 0
+ 0, 2, // 1
+ 2, // 2
+ 0, 1, 2, // 3
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ // Create field and set constraint sizes
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ int iV=0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ section->addConstraintDimension(*v_iter, nconstraints[iV++]);
+ field.allocate();
+ int index = 0;
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, index += nconstraints[i++])
+ section->setConstraintDof(*v_iter, &constraints[index]);
+ field.zero();
+
+ double_array values(fiberDim);
+ i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePointAll(*v_iter, &values[0]);
+ } // for
+
+ field.zeroAll();
+
+ const double tolerance = 1.0e-6;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[iDim], tolerance);
+ } // for
+ } // for
+} // testZeroAll
+
+// ----------------------------------------------------------------------
+// Test complete().
+void
+pylith::topology::TestMultiFieldMesh::testComplete(void)
+{ // testComplete
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.complete();
+
+ // Expect no change for this serial test
+ i = 0;
+ const double tolerance = 1.0e-6;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ } // for
+ } // for
+} // testComplete
+
+// ----------------------------------------------------------------------
+// Test copy().
+void
+pylith::topology::TestMultiFieldMesh::testCopy(void)
+{ // testCopy
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<Mesh> fieldSrc(mesh);
+ { // Setup source field
+ fieldSrc.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ fieldSrc.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+ } // Setup source field
+
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ field.copy(fieldSrc);
+
+ int i = 0;
+ double_array values(fiberDim);
+ const double tolerance = 1.0e-6;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ } // for
+ } // for
+} // testCopy
+
+// ----------------------------------------------------------------------
+// Test operator+=().
+void
+pylith::topology::TestMultiFieldMesh::testOperatorAdd(void)
+{ // testOperateAdd
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesA[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+ const double valuesB[] = {
+ 10.1, 20.2, 30.3,
+ 10.2, 20.3, 30.4,
+ 10.3, 20.4, 30.5,
+ 10.4, 20.5, 30.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<Mesh> fieldSrc(mesh);
+ { // Setup source field
+ fieldSrc.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ fieldSrc.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesA[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+ } // Setup source field
+
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ { // Setup destination field
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesB[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+ } // Setup destination field
+
+ field += fieldSrc;
+
+ int i = 0;
+ double_array values(fiberDim);
+ const double tolerance = 1.0e-6;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ const double valueE = valuesA[i] + valuesB[i];
+ ++i;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[iDim], tolerance);
+ } // for
+ } // for
+} // testOperateAdd
+
+// ----------------------------------------------------------------------
+// Test dimensionalize().
+void
+pylith::topology::TestMultiFieldMesh::testDimensionalize(void)
+{ // testDimensionalize
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.scale(scale);
+ field.addDimensionOkay(true);
+ field.dimensionalize();
+
+ i = 0;
+ const double tolerance = 1.0e-6;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ const double valueE = valuesNondim[i++]*scale;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[iDim], tolerance);
+ } // for
+ } // for
+
+} // testDimensionalize
+
+// ----------------------------------------------------------------------
+// Test view().
+void
+pylith::topology::TestMultiFieldMesh::testView(void)
+{ // testView
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.view("Testing view");
+} // testView
+
+// ----------------------------------------------------------------------
+// Test createScatter().
+void
+pylith::topology::TestMultiFieldMesh::testCreateScatter(void)
+{ // testCreateScatter
+ const int fiberDim = 3;
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ MultiField<Mesh> field(mesh);
+ const std::string& label = "field A";
+ field.label(label.c_str());
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const int sizeE = vertices->size() * fiberDim;
+
+ CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
+ field.createScatter();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+ const MultiField<Mesh>::ScatterInfo& sinfo = field._getScatter("");
+ CPPUNIT_ASSERT(sinfo.scatter);
+ CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.vector);
+
+ int size = 0;
+ VecGetSize(sinfo.vector, &size);
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+ // Check vector name
+ const char* vecname = 0;
+ PetscObjectGetName((PetscObject)sinfo.vector, &vecname);
+ CPPUNIT_ASSERT_EQUAL(label, std::string(vecname));
+
+ // Make sure we can do multiple calls to createScatter().
+ field.createScatter();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+
+ // Create another scatter.
+ field.createScatter("B");
+ CPPUNIT_ASSERT_EQUAL(size_t(2), field._scatters.size());
+ const MultiField<Mesh>::ScatterInfo& sinfoB = field._getScatter("B");
+ CPPUNIT_ASSERT(sinfoB.scatter);
+ CPPUNIT_ASSERT(sinfoB.scatterVec);
+ CPPUNIT_ASSERT(sinfoB.vector);
+
+ MultiField<Mesh> field2(mesh);
+ field2.cloneSection(field);
+ CPPUNIT_ASSERT_EQUAL(size_t(2), field2._scatters.size());
+
+ const MultiField<Mesh>::ScatterInfo& sinfo2 = field2._getScatter("");
+ CPPUNIT_ASSERT(sinfo2.scatter);
+ CPPUNIT_ASSERT(sinfo2.scatterVec);
+ CPPUNIT_ASSERT(sinfo2.vector);
+
+ const MultiField<Mesh>::ScatterInfo& sinfo2B = field2._getScatter("B");
+ CPPUNIT_ASSERT(sinfo2B.scatter);
+ CPPUNIT_ASSERT(sinfo2B.scatterVec);
+ CPPUNIT_ASSERT(sinfo2B.vector);
+
+} // testCreateScatter
+
+// ----------------------------------------------------------------------
+// Test createScatterWithBC().
+void
+pylith::topology::TestMultiFieldMesh::testCreateScatterWithBC(void)
+{ // testCreateScatterWithBC
+ const int fiberDim = 3;
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ MultiField<Mesh> field(mesh);
+ const std::string& label = "field A";
+ field.label(label.c_str());
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const int sizeE = vertices->size() * fiberDim;
+
+ CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
+ field.createScatterWithBC();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+ const MultiField<Mesh>::ScatterInfo& sinfo = field._getScatter("");
+ CPPUNIT_ASSERT(sinfo.scatter);
+ CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.vector);
+
+ int size = 0;
+ VecGetSize(sinfo.vector, &size);
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+ // Check vector name
+ const char* vecname = 0;
+ PetscObjectGetName((PetscObject)sinfo.vector, &vecname);
+ CPPUNIT_ASSERT_EQUAL(label, std::string(vecname));
+
+ // Make sure we can do multiple calls to createScatterWithBC().
+ field.createScatterWithBC();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+
+ // Create another scatter.
+ field.createScatterWithBC("B");
+ CPPUNIT_ASSERT_EQUAL(size_t(2), field._scatters.size());
+ const MultiField<Mesh>::ScatterInfo& sinfoB = field._getScatter("B");
+ CPPUNIT_ASSERT(sinfoB.scatter);
+ CPPUNIT_ASSERT(sinfoB.scatterVec);
+ CPPUNIT_ASSERT(sinfoB.vector);
+
+ MultiField<Mesh> field2(mesh);
+ field2.cloneSection(field);
+ CPPUNIT_ASSERT_EQUAL(size_t(2), field2._scatters.size());
+
+ const MultiField<Mesh>::ScatterInfo& sinfo2 = field2._getScatter("");
+ CPPUNIT_ASSERT(sinfo2.scatter);
+ CPPUNIT_ASSERT(sinfo2.scatterVec);
+ CPPUNIT_ASSERT(sinfo2.vector);
+
+ const MultiField<Mesh>::ScatterInfo& sinfo2B = field2._getScatter("B");
+ CPPUNIT_ASSERT(sinfo2B.scatter);
+ CPPUNIT_ASSERT(sinfo2B.scatterVec);
+ CPPUNIT_ASSERT(sinfo2B.vector);
+
+} // testCreateScatterWithBC
+
+// ----------------------------------------------------------------------
+// Test vector().
+void
+pylith::topology::TestMultiFieldMesh::testVector(void)
+{ // testVector
+ const int fiberDim = 3;
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+
+ CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
+ field.createScatter();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+ const MultiField<Mesh>::ScatterInfo& sinfo = field._getScatter("");
+ CPPUNIT_ASSERT(sinfo.scatter);
+ CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.vector);
+
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ const PetscVec vec = field.vector();
+ CPPUNIT_ASSERT_EQUAL(sinfo.vector, vec);
+ int size = 0;
+ VecGetSize(vec, &size);
+ const int sizeE = vertices->size() * fiberDim;
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+} // testVector
+
+// ----------------------------------------------------------------------
+// Test scatterSectionToVector().
+void
+pylith::topology::TestMultiFieldMesh::testScatterSectionToVector(void)
+{ // testScatterSectionToVector
+ const char* context = "abc";
+ const int fiberDim = 3;
+ const double valuesE[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesE[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.createScatter(context);
+ field.scatterSectionToVector(context);
+ const PetscVec vec = field.vector(context);
+ CPPUNIT_ASSERT(0 != vec);
+ int size = 0;
+ VecGetSize(vec, &size);
+ double* valuesVec = 0;
+ VecGetArray(vec, &valuesVec);
+
+ const double tolerance = 1.0e-06;
+ const int sizeE = vertices->size() * fiberDim;
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+ for (int i=0; i < sizeE; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i], valuesVec[i], tolerance);
+ VecRestoreArray(vec, &valuesVec);
+} // testScatterSectionToVector
+
+// ----------------------------------------------------------------------
+// Test scatterVectorToSection().
+void
+pylith::topology::TestMultiFieldMesh::testScatterVectorToSection(void)
+{ // testScatterVectorToSection
+ const char* context = "abcd";
+ const int fiberDim = 3;
+ const double valuesE[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ 1.4, 2.5, 3.6,
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ MultiField<Mesh> field(mesh);
+ field.newSection(MultiField<Mesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<Mesh::RealSection>& section = field.section();
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+
+ field.createScatter(context);
+ const PetscVec vec = field.vector(context);
+ CPPUNIT_ASSERT(0 != vec);
+ int size = 0;
+ VecGetSize(vec, &size);
+ double* valuesVec = 0;
+ VecGetArray(vec, &valuesVec);
+
+ const double tolerance = 1.0e-06;
+ const int sizeE = vertices->size() * fiberDim;
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+ for (int i=0; i < sizeE; ++i)
+ valuesVec[i] = valuesE[i];
+ VecRestoreArray(vec, &valuesVec);
+
+ field.createScatter(context);
+ field.scatterVectorToSection(context);
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], fiberDim);
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i++], values[iDim], tolerance);
+ } // for
+} // testScatterVectorToSection
+
+// ----------------------------------------------------------------------
+// Test splitDefault().
+void
+pylith::topology::TestMultiFieldMesh::testSplitDefault(void)
+{ // testSplitDefault
+ const int spaceDim = _TestMultiFieldMesh::cellDim;
+ const int numFibrations = spaceDim;
+ const int nconstraints[4] = { 1, 2, 0, 1 };
+ const int constraints[4] = {
+ 1, // 0
+ 0, 1, // 1
+ // 2
+ 0, // 3
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ // Create field with atlas to use to create new field
+ MultiField<Mesh> fieldSrc(mesh);
+ { // Setup source field
+ fieldSrc.newSection(MultiField<Mesh>::VERTICES_FIELD, spaceDim);
+ fieldSrc.splitDefault();
+ const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ int iV=0;
+ int iC=0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, ++iV) {
+ const int nconstraintsVertex = nconstraints[iV];
+ section->addConstraintDimension(*v_iter, nconstraintsVertex);
+ for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint) {
+ const int fibration = constraints[iC++];
+ section->addConstraintDimension(*v_iter, 1, fibration);
+ } // for
+ } // for
+ fieldSrc.allocate();
+
+ iC = 0;
+ iV = 0;
+ int zero = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, ++iV) {
+ const int nconstraintsVertex = nconstraints[iV];
+ if (nconstraintsVertex > 0)
+ section->setConstraintDof(*v_iter, &constraints[iC]);
+ for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint) {
+ const int fibration = constraints[iC++];
+ section->setConstraintDof(*v_iter, &zero, fibration);
+ } // for
+ } // for
+ } // Setup source field
+
+ const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ CPPUNIT_ASSERT_EQUAL(numFibrations, section->getNumSpaces());
+
+ for (int fibration=0; fibration < spaceDim; ++fibration) {
+ const ALE::Obj<Mesh::RealSection>& sectionSplit = section->getFibration(fibration);
+ CPPUNIT_ASSERT(!sectionSplit.isNull());
+ CPPUNIT_ASSERT(!vertices.isNull());
+ int iV = 0;
+ int iC = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, ++iV) {
+ CPPUNIT_ASSERT_EQUAL(1, section->getFiberDimension(*v_iter, fibration));
+ bool isConstrained = false;
+ const int nconstraintsVertex = nconstraints[iV];
+ for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint)
+ if (constraints[iC++] == fibration)
+ isConstrained = true;
+ const int constraintDimE = (!isConstrained) ? 0 : 1;
+ CPPUNIT_ASSERT_EQUAL(constraintDimE,
+ section->getConstraintDimension(*v_iter, fibration));
+ } // for
+ } // for
+} // testSplitDefault
+
+// ----------------------------------------------------------------------
+// Test cloneSection() with split field.
+void
+pylith::topology::TestMultiFieldMesh::testCloneSectionSplit(void)
+{ // testCloneSectionSplit
+ const int spaceDim = _TestMultiFieldMesh::cellDim;
+ const int numFibrations = spaceDim;
+ const int nconstraints[4] = { 1, 2, 0, 1 };
+ const int constraints[4] = {
+ 1, // 0
+ 0, 1, // 1
+ // 2
+ 0, // 3
+ };
+
+ Mesh mesh;
+ _buildMesh(&mesh);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ const ALE::Obj<Mesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ // Create field with atlas to use to create new field
+ MultiField<Mesh> fieldSrc(mesh);
+ { // Setup source field
+ fieldSrc.newSection(MultiField<Mesh>::VERTICES_FIELD, spaceDim);
+ fieldSrc.splitDefault();
+ const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ int iV=0;
+ int iC=0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, ++iV) {
+ const int nconstraintsVertex = nconstraints[iV];
+ section->addConstraintDimension(*v_iter, nconstraintsVertex);
+ for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint) {
+ const int fibration = constraints[iC++];
+ section->addConstraintDimension(*v_iter, 1, fibration);
+ } // for
+ } // for
+ fieldSrc.allocate();
+
+ iC = 0;
+ iV = 0;
+ int zero = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, ++iV) {
+ const int nconstraintsVertex = nconstraints[iV];
+ if (nconstraintsVertex > 0)
+ section->setConstraintDof(*v_iter, &constraints[iC]);
+ for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint) {
+ const int fibration = constraints[iC++];
+ section->setConstraintDof(*v_iter, &zero, fibration);
+ } // for
+ } // for
+ } // Setup source field
+
+ MultiField<Mesh> field(mesh);
+ field.cloneSection(fieldSrc);
+
+ const ALE::Obj<Mesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ CPPUNIT_ASSERT_EQUAL(numFibrations, section->getNumSpaces());
+
+ for (int fibration=0; fibration < spaceDim; ++fibration) {
+ const ALE::Obj<Mesh::RealSection>& sectionSplit = section->getFibration(fibration);
+ CPPUNIT_ASSERT(!sectionSplit.isNull());
+ CPPUNIT_ASSERT(!vertices.isNull());
+ int iV = 0;
+ int iC = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, ++iV) {
+ CPPUNIT_ASSERT_EQUAL(1, section->getFiberDimension(*v_iter, fibration));
+ bool isConstrained = false;
+ const int nconstraintsVertex = nconstraints[iV];
+ for (int iConstraint=0; iConstraint < nconstraintsVertex; ++iConstraint)
+ if (constraints[iC++] == fibration)
+ isConstrained = true;
+ const int constraintDimE = (!isConstrained) ? 0 : 1;
+ CPPUNIT_ASSERT_EQUAL(constraintDimE,
+ section->getConstraintDimension(*v_iter, fibration));
+ } // for
+ } // for
+} // testCloneSectionSplit
+
+// ----------------------------------------------------------------------
+void
+pylith::topology::TestMultiFieldMesh::_buildMesh(Mesh* mesh)
+{ // _buildMesh
+ assert(0 != mesh);
+
+ mesh->createSieveMesh(_TestMultiFieldMesh::cellDim);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh->sieveMesh();
+
+ ALE::Obj<Mesh::SieveMesh::sieve_type> sieve =
+ new Mesh::SieveMesh::sieve_type(sieveMesh->comm());
+ CPPUNIT_ASSERT(!sieve.isNull());
+
+ ALE::Obj<SieveFlexMesh::sieve_type> s =
+ new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
+
+ const int cellDim = _TestMultiFieldMesh::cellDim;
+ const int ncells = _TestMultiFieldMesh::ncells;
+ const int* cells = _TestMultiFieldMesh::cells;
+ const int nvertices = _TestMultiFieldMesh::nvertices;
+ const int ncorners = _TestMultiFieldMesh::ncorners;
+ const int spaceDim = _TestMultiFieldMesh::cellDim;
+ const double* coordinates = _TestMultiFieldMesh::coordinates;
+ const bool interpolate = false;
+ ALE::SieveBuilder<SieveFlexMesh>::buildTopology(s, cellDim, ncells, (int*) cells,
+ nvertices, interpolate,
+ ncorners);
+ std::map<Mesh::SieveMesh::point_type,Mesh::SieveMesh::point_type> renumbering;
+ ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
+ sieveMesh->setSieve(sieve);
+ sieveMesh->stratify();
+ ALE::SieveBuilder<Mesh::SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
+ coordinates);
+
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(spaceDim);
+ cs.initialize();
+ mesh->coordsys(&cs);
+} // _buildMesh
+
+
+// End of file
Added: short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldMesh.hh
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldMesh.hh (rev 0)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldMesh.hh 2011-05-21 19:47:16 UTC (rev 18415)
@@ -0,0 +1,189 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/topology/TestMultiFieldMesh.hh
+ *
+ * @brief C++ unit testing for Field.
+ */
+
+#if !defined(pylith_topology_testmultifieldmesh_hh)
+#define pylith_topology_testmultifieldmesh_hh
+
+// Include directives ---------------------------------------------------
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "pylith/topology/topologyfwd.hh" // forward declarations
+
+// Forward declarations -------------------------------------------------
+/// Namespace for pylith package
+namespace pylith {
+ namespace topology {
+ class TestMultiFieldMesh;
+ } // topology
+} // pylith
+
+// TestMultiFieldMesh -------------------------------------------------------------
+/// C++ unit testing for Field.
+class pylith::topology::TestMultiFieldMesh : public CppUnit::TestFixture
+{ // class TestMultiFieldMesh
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestMultiFieldMesh );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testMesh );
+ CPPUNIT_TEST( testLabel );
+ CPPUNIT_TEST( testVectorFieldType );
+ CPPUNIT_TEST( testScale );
+ CPPUNIT_TEST( testAddDimensionsOkay );
+ CPPUNIT_TEST( testSpaceDim );
+ CPPUNIT_TEST( testNewSection );
+ CPPUNIT_TEST( testNewSectionPoints );
+ CPPUNIT_TEST( testNewSectionPointsArray );
+ CPPUNIT_TEST( testNewSectionDomain );
+ CPPUNIT_TEST( testNewSectionField );
+ CPPUNIT_TEST( testCloneSection );
+ CPPUNIT_TEST( testClear );
+ CPPUNIT_TEST( testAllocate );
+ CPPUNIT_TEST( testZero );
+ CPPUNIT_TEST( testZeroAll );
+ CPPUNIT_TEST( testComplete );
+ CPPUNIT_TEST( testCopy );
+ CPPUNIT_TEST( testOperatorAdd );
+ CPPUNIT_TEST( testDimensionalize );
+ CPPUNIT_TEST( testView );
+ CPPUNIT_TEST( testCreateScatter );
+ CPPUNIT_TEST( testCreateScatterWithBC );
+ CPPUNIT_TEST( testVector );
+ CPPUNIT_TEST( testScatterSectionToVector );
+ CPPUNIT_TEST( testScatterVectorToSection );
+ CPPUNIT_TEST( testSplitDefault );
+ CPPUNIT_TEST( testCloneSectionSplit );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test section().
+ void testSection(void);
+
+ /// Test mesh().
+ void testMesh(void);
+
+ /// Test label().
+ void testLabel(void);
+
+ /// Test vectorFieldType().
+ void testVectorFieldType(void);
+
+ /// Test scale().
+ void testScale(void);
+
+ /// Test addDimensionsOkay().
+ void testAddDimensionsOkay(void);
+
+ /// Test spaceDim().
+ void testSpaceDim(void);
+
+ /// Test newSection().
+ void testNewSection(void);
+
+ /// Test newSection(points).
+ void testNewSectionPoints(void);
+
+ /// Test newSection(int_array).
+ void testNewSectionPointsArray(void);
+
+ /// Test newSection(domain).
+ void testNewSectionDomain(void);
+
+ /// Test newSection(field).
+ void testNewSectionField(void);
+
+ /// Test cloneSection().
+ void testCloneSection(void);
+
+ /// Test clear().
+ void testClear(void);
+
+ /// Test allocate().
+ void testAllocate(void);
+
+ /// Test zero().
+ void testZero(void);
+
+ /// Test zeroAll().
+ void testZeroAll(void);
+
+ /// Test complete().
+ void testComplete(void);
+
+ /// Test copy().
+ void testCopy(void);
+
+ /// Test operator+=().
+ void testOperatorAdd(void);
+
+ /// Test dimensionalize().
+ void testDimensionalize(void);
+
+ /// Test view().
+ void testView(void);
+
+ /// Test createScatter().
+ void testCreateScatter(void);
+
+ /// Test createScatterWithBC().
+ void testCreateScatterWithBC(void);
+
+ /// Test vector().
+ void testVector(void);
+
+ /// Test scatterSectionToVector().
+ void testScatterSectionToVector(void);
+
+ /// Test scatterVectorToSection().
+ void testScatterVectorToSection(void);
+
+ /// Test splitDefault().
+ void testSplitDefault(void);
+
+ /// Test cloneSection() with split field.
+ void testCloneSectionSplit(void);
+
+// PRIVATE METHODS /////////////////////////////////////////////////////
+private :
+
+ /** Build mesh.
+ *
+ * @param mesh Finite-element mesh.
+ */
+ static
+ void _buildMesh(Mesh* mesh);
+
+}; // class TestMultiFieldMesh
+
+#endif // pylith_topology_testmultifieldmesh_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldSubMesh.cc
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldSubMesh.cc (rev 0)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldSubMesh.cc 2011-05-21 19:47:16 UTC (rev 18415)
@@ -0,0 +1,1020 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestMultiFieldSubMesh.hh" // Implementation of class methods
+
+#include "pylith/topology/MultiField.hh" // USES MultiField
+#include "pylith/topology/SubMesh.hh" // USES SubMesh
+
+#include "pylith/utils/array.hh" // USES double_array
+
+#include "spatialdata/geocoords/CSCart.hh" // USES CSCart
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::topology::TestMultiFieldSubMesh );
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
+
+// ----------------------------------------------------------------------
+namespace pylith {
+ namespace topology {
+ namespace _TestMultiFieldSubMesh {
+ const int cellDim = 2;
+ const int nvertices = 4;
+ const int ncells = 2;
+ const int ncorners = 3;
+ const int cells[] = {
+ 0, 1, 3,
+ 0, 3, 2,
+ };
+ const double coordinates[] = {
+ 0.0, 0.0,
+ 1.0, 0.0,
+ 0.0, 1.0,
+ 1.0, 1.0,
+ };
+ const char* label = "bc";
+ const int groupSize = 3;
+ const int groupVertices[] = {
+ 1, 2, 3
+ };
+ const int submeshVertices[] = {
+ 3, 4, 5
+ };
+ } // _TestMultiFieldSubMesh
+ } // topology
+} // pylith
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::topology::TestMultiFieldSubMesh::testConstructor(void)
+{ // testConstructor
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test newSection().
+void
+pylith::topology::TestMultiFieldSubMesh::testSection(void)
+{ // testSection
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(section.isNull());
+} // testSection
+
+// ----------------------------------------------------------------------
+// Test mesh().
+void
+pylith::topology::TestMultiFieldSubMesh::testMesh(void)
+{ // testMesh
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+
+ const SubMesh& mesh2 = field.mesh();
+ CPPUNIT_ASSERT_EQUAL(_TestMultiFieldSubMesh::cellDim-1, mesh2.dimension());
+} // testMesh
+
+// ----------------------------------------------------------------------
+// Test spaceDim().
+void
+pylith::topology::TestMultiFieldSubMesh::testSpaceDim(void)
+{ // testSpaceDim
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+
+ CPPUNIT_ASSERT_EQUAL(_TestMultiFieldSubMesh::cellDim, field.spaceDim());
+} // testSpaceDim
+
+// ----------------------------------------------------------------------
+// Test newSection().
+void
+pylith::topology::TestMultiFieldSubMesh::testNewSection(void)
+{ // testNewSection
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+
+ field.newSection();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+} // testNewSection
+
+// ----------------------------------------------------------------------
+// Test newSection(points).
+void
+pylith::topology::TestMultiFieldSubMesh::testNewSectionPoints(void)
+{ // testNewSectionPoints
+ const int fiberDim = 2;
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ MultiField<SubMesh> field(submesh);
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ field.newSection(vertices, fiberDim);
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
+} // testNewSectionPoints
+
+// ----------------------------------------------------------------------
+// Test newSection(domain).
+void
+pylith::topology::TestMultiFieldSubMesh::testNewSectionDomain(void)
+{ // testNewSectionDomain
+ const int fiberDim = 2;
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ MultiField<SubMesh> field(submesh);
+ field.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
+} // testNewSectionDomain
+
+// ----------------------------------------------------------------------
+// Test newSection(field).
+void
+pylith::topology::TestMultiFieldSubMesh::testNewSectionField(void)
+{ // testNewSectionField
+ const int fiberDim = 3;
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ // Create field with atlas to use to create new field
+ MultiField<SubMesh> fieldSrc(submesh);
+ fieldSrc.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ const ALE::Obj<SubMesh::RealSection>& sectionSrc = fieldSrc.section();
+ CPPUNIT_ASSERT(!sectionSrc.isNull());
+
+ const int fiberDim2 = 4;
+ MultiField<SubMesh> field(submesh);
+ field.newSection(fieldSrc, fiberDim2);
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ CPPUNIT_ASSERT_EQUAL(fiberDim2, section->getFiberDimension(*v_iter));
+} // testNewSectionChart
+
+// ----------------------------------------------------------------------
+// Test cloneSection().
+void
+pylith::topology::TestMultiFieldSubMesh::testCloneSection(void)
+{ // testCloneSection
+ const int fiberDim = 3;
+ const int nconstraints[] = { 0, 2, 1, 3 };
+ const int constraints[] = {
+ // 0
+ 0, 2, // 1
+ 2, // 2
+ 0, 1, 2, // 3
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ // Create field with atlas to use to create new field
+ MultiField<SubMesh> fieldSrc(submesh);
+ { // Setup source field
+ fieldSrc.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ const ALE::Obj<SubMesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ int iV=0;
+
+ CPPUNIT_ASSERT(!vertices.isNull());
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter)
+ section->addConstraintDimension(*v_iter, nconstraints[iV++]);
+ fieldSrc.allocate();
+
+ int index = 0;
+ int i = 0;
+ for (Mesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter, index += nconstraints[i++])
+ section->setConstraintDof(*v_iter, &constraints[index]);
+ fieldSrc.zero();
+ } // Setup source field
+
+
+ MultiField<SubMesh> field(submesh);
+ field.cloneSection(fieldSrc);
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ int iV = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ CPPUNIT_ASSERT_EQUAL(fiberDim, section->getFiberDimension(*v_iter));
+ CPPUNIT_ASSERT_EQUAL(nconstraints[iV++],
+ section->getConstraintDimension(*v_iter));
+ } // for
+} // testCloneSection
+
+// ----------------------------------------------------------------------
+// Test clear().
+void
+pylith::topology::TestMultiFieldSubMesh::testClear(void)
+{ // testClear
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+
+ field.scale(2.0);
+ field.vectorFieldType(MultiField<SubMesh>::TENSOR);
+ field.addDimensionOkay(true);
+
+ field.clear();
+
+ CPPUNIT_ASSERT_EQUAL(1.0, field._metadata.scale);
+ CPPUNIT_ASSERT_EQUAL(MultiField<SubMesh>::OTHER, field._metadata.vectorFieldType);
+ CPPUNIT_ASSERT_EQUAL(false, field._metadata.dimsOkay);
+} // testClear
+
+// ----------------------------------------------------------------------
+// Test allocate().
+void
+pylith::topology::TestMultiFieldSubMesh::testAllocate(void)
+{ // testAllocate
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<SubMesh> field(submesh);
+ field.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ const double tolerance = 1.0e-6;
+ i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ } // for
+ } // for
+} // testAllocate
+
+// ----------------------------------------------------------------------
+// Test zero().
+void
+pylith::topology::TestMultiFieldSubMesh::testZero(void)
+{ // testZero
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<SubMesh> field(submesh);
+ field.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.zero();
+
+ const double tolerance = 1.0e-6;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, values[iDim], tolerance);
+ } // for
+ } // for
+} // testZero
+
+// ----------------------------------------------------------------------
+// Test complete().
+void
+pylith::topology::TestMultiFieldSubMesh::testComplete(void)
+{ // testComplete
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<SubMesh> field(submesh);
+ field.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.complete();
+
+ // Expect no change for this serial test
+ i = 0;
+ const double tolerance = 1.0e-6;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ } // for
+ } // for
+} // testComplete
+
+// ----------------------------------------------------------------------
+// Test copy().
+void
+pylith::topology::TestMultiFieldSubMesh::testCopy(void)
+{ // testCopy
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<SubMesh> fieldSrc(submesh);
+ { // Setup source field
+ fieldSrc.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ fieldSrc.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+ } // Setup source field
+
+ MultiField<SubMesh> field(submesh);
+ field.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ field.copy(fieldSrc);
+
+ int i = 0;
+ double_array values(fiberDim);
+ const double tolerance = 1.0e-6;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesNondim[i++], values[iDim], tolerance);
+ } // for
+ } // for
+} // testCopy
+
+// ----------------------------------------------------------------------
+// Test operator+=().
+void
+pylith::topology::TestMultiFieldSubMesh::testOperatorAdd(void)
+{ // testOperateAdd
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesA[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+ const double valuesB[] = {
+ 10.1, 20.2, 30.3,
+ 10.2, 20.3, 30.4,
+ 10.3, 20.4, 30.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ MultiField<SubMesh> fieldSrc(submesh);
+ { // Setup source field
+ fieldSrc.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ fieldSrc.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = fieldSrc.section();
+ CPPUNIT_ASSERT(!section.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesA[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+ } // Setup source field
+
+ MultiField<SubMesh> field(submesh);
+ field.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ { // Setup destination field
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesB[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+ } // Setup destination field
+
+ field += fieldSrc;
+
+ int i = 0;
+ double_array values(fiberDim);
+ const double tolerance = 1.0e-6;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ const double valueE = valuesA[i] + valuesB[i];
+ ++i;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[iDim], tolerance);
+ } // for
+ } // for
+} // testOperateAdd
+
+// ----------------------------------------------------------------------
+// Test dimensionalize().
+void
+pylith::topology::TestMultiFieldSubMesh::testDimensionalize(void)
+{ // testDimensionalize
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ MultiField<SubMesh> field(submesh);
+ field.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.scale(scale);
+ field.addDimensionOkay(true);
+ field.dimensionalize();
+
+ i = 0;
+ const double tolerance = 1.0e-6;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], values.size());
+ for (int iDim=0; iDim < fiberDim; ++iDim) {
+ const double valueE = valuesNondim[i++]*scale;
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valueE, values[iDim], tolerance);
+ } // for
+ } // for
+
+} // testDimensionalize
+
+// ----------------------------------------------------------------------
+// Test view().
+void
+pylith::topology::TestMultiFieldSubMesh::testView(void)
+{ // testView
+ const int fiberDim = 3;
+ const double scale = 2.0;
+ const double valuesNondim[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ const ALE::Obj<SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ MultiField<SubMesh> field(submesh);
+ field.newSection(MultiField<SubMesh>::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesNondim[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.view("Testing view");
+} // testView
+
+// ----------------------------------------------------------------------
+// Test createScatter().
+void
+pylith::topology::TestMultiFieldSubMesh::testCreateScatter(void)
+{ // testCreateScatter
+ const int fiberDim = 3;
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+ field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
+ field.allocate();
+
+ const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const int sizeE = vertices->size() * fiberDim;
+
+ CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
+ field.createScatter();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+ const MultiField<SubMesh>::ScatterInfo& sinfo = field._getScatter("");
+ CPPUNIT_ASSERT(sinfo.scatter);
+ CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.vector);
+
+ int size = 0;
+ VecGetSize(sinfo.vector, &size);
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+ // Make sure we can do multiple calls to createScatter().
+ field.createScatter();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+
+ // Create another scatter.
+ field.createScatter("B");
+ CPPUNIT_ASSERT_EQUAL(size_t(2), field._scatters.size());
+ const MultiField<SubMesh>::ScatterInfo& sinfoB = field._getScatter("B");
+ CPPUNIT_ASSERT(sinfoB.scatter);
+ CPPUNIT_ASSERT(sinfoB.scatterVec);
+ CPPUNIT_ASSERT(sinfoB.vector);
+
+ MultiField<SubMesh> field2(submesh);
+ field2.cloneSection(field);
+ CPPUNIT_ASSERT_EQUAL(size_t(2), field2._scatters.size());
+
+ const MultiField<SubMesh>::ScatterInfo& sinfo2 = field2._getScatter("");
+ CPPUNIT_ASSERT(sinfo2.scatter);
+ CPPUNIT_ASSERT(sinfo2.scatterVec);
+ CPPUNIT_ASSERT(sinfo2.vector);
+
+ const MultiField<SubMesh>::ScatterInfo& sinfo2B = field2._getScatter("B");
+ CPPUNIT_ASSERT(sinfo2B.scatter);
+ CPPUNIT_ASSERT(sinfo2B.scatterVec);
+ CPPUNIT_ASSERT(sinfo2B.vector);
+} // testCreateScatter
+
+// ----------------------------------------------------------------------
+// Test createScatterWithBC().
+void
+pylith::topology::TestMultiFieldSubMesh::testCreateScatterWithBC(void)
+{ // testCreateScatterWithBC
+ const int fiberDim = 3;
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+ field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
+ field.allocate();
+
+ const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const int sizeE = vertices->size() * fiberDim;
+
+ CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
+ field.createScatterWithBC();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+ const MultiField<SubMesh>::ScatterInfo& sinfo = field._getScatter("");
+ CPPUNIT_ASSERT(sinfo.scatter);
+ CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.vector);
+
+ int size = 0;
+ VecGetSize(sinfo.vector, &size);
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+ // Make sure we can do multiple calls to createScatterWithBC().
+ field.createScatterWithBC();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+
+ // Create another scatter.
+ field.createScatterWithBC("B");
+ CPPUNIT_ASSERT_EQUAL(size_t(2), field._scatters.size());
+ const MultiField<SubMesh>::ScatterInfo& sinfoB = field._getScatter("B");
+ CPPUNIT_ASSERT(sinfoB.scatter);
+ CPPUNIT_ASSERT(sinfoB.scatterVec);
+ CPPUNIT_ASSERT(sinfoB.vector);
+
+ MultiField<SubMesh> field2(submesh);
+ field2.cloneSection(field);
+ CPPUNIT_ASSERT_EQUAL(size_t(2), field2._scatters.size());
+
+ const MultiField<SubMesh>::ScatterInfo& sinfo2 = field2._getScatter("");
+ CPPUNIT_ASSERT(sinfo2.scatter);
+ CPPUNIT_ASSERT(sinfo2.scatterVec);
+ CPPUNIT_ASSERT(sinfo2.vector);
+
+ const MultiField<SubMesh>::ScatterInfo& sinfo2B = field2._getScatter("B");
+ CPPUNIT_ASSERT(sinfo2B.scatter);
+ CPPUNIT_ASSERT(sinfo2B.scatterVec);
+ CPPUNIT_ASSERT(sinfo2B.vector);
+} // testCreateScatterWithBC
+
+// ----------------------------------------------------------------------
+// Test vector().
+void
+pylith::topology::TestMultiFieldSubMesh::testVector(void)
+{ // testVector
+ const int fiberDim = 3;
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+
+ MultiField<SubMesh> field(submesh);
+ field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
+ field.allocate();
+
+ CPPUNIT_ASSERT_EQUAL(size_t(0), field._scatters.size());
+ field.createScatter();
+ CPPUNIT_ASSERT_EQUAL(size_t(1), field._scatters.size());
+ const MultiField<SubMesh>::ScatterInfo& sinfo = field._getScatter("");
+ CPPUNIT_ASSERT(sinfo.scatter);
+ CPPUNIT_ASSERT(sinfo.scatterVec);
+ CPPUNIT_ASSERT(sinfo.vector);
+
+ const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ const PetscVec vec = field.vector();
+ CPPUNIT_ASSERT_EQUAL(sinfo.vector, vec);
+ int size = 0;
+ VecGetSize(vec, &size);
+ const int sizeE = vertices->size() * fiberDim;
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+} // testVector
+
+// ----------------------------------------------------------------------
+// Test scatterSectionToVector().
+void
+pylith::topology::TestMultiFieldSubMesh::testScatterSectionToVector(void)
+{ // testScatterSectionToVector
+ const char* context = "abc";
+ const int fiberDim = 3;
+ const double valuesE[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+ const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
+ field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+
+ double_array values(fiberDim);
+ int i = 0;
+ for (SubMesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ values[iDim] = valuesE[i++];
+ section->updatePoint(*v_iter, &values[0]);
+ } // for
+
+ field.createScatter(context);
+ field.scatterSectionToVector(context);
+ const PetscVec vec = field.vector(context);
+ CPPUNIT_ASSERT(0 != vec);
+ int size = 0;
+ VecGetSize(vec, &size);
+ double* valuesVec = 0;
+ VecGetArray(vec, &valuesVec);
+
+ const double tolerance = 1.0e-06;
+ const int sizeE = vertices->size() * fiberDim;
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+ for (int i=0; i < sizeE; ++i)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i], valuesVec[i], tolerance);
+ VecRestoreArray(vec, &valuesVec);
+} // testScatterSectionToVector
+
+// ----------------------------------------------------------------------
+// Test scatterVectorToSection().
+void
+pylith::topology::TestMultiFieldSubMesh::testScatterVectorToSection(void)
+{ // testScatterVectorToSection
+ const char* context = "abcd";
+ const int fiberDim = 3;
+ const double valuesE[] = {
+ 1.1, 2.2, 3.3,
+ 1.2, 2.3, 3.4,
+ 1.3, 2.4, 3.5,
+ };
+
+ Mesh mesh;
+ SubMesh submesh;
+ _buildMesh(&mesh, &submesh);
+ MultiField<SubMesh> field(submesh);
+ field.newSection(FieldBase::VERTICES_FIELD, fiberDim);
+ field.allocate();
+ field.createScatter(context);
+
+ const ALE::Obj<SubMesh::SieveMesh>& sieveMesh = submesh.sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+ const ALE::Obj<SubMesh::SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+
+ const PetscVec vec = field.vector(context);
+ CPPUNIT_ASSERT(0 != vec);
+ int size = 0;
+ VecGetSize(vec, &size);
+ const int sizeE = vertices->size() * fiberDim;
+ CPPUNIT_ASSERT_EQUAL(sizeE, size);
+
+ const double tolerance = 1.0e-06;
+ double* valuesVec = 0;
+ VecGetArray(vec, &valuesVec);
+ for (int i=0; i < sizeE; ++i)
+ valuesVec[i] = valuesE[i];
+ VecRestoreArray(vec, &valuesVec);
+
+ field.scatterVectorToSection(context);
+
+ double_array values(fiberDim);
+ int i = 0;
+ const ALE::Obj<SubMesh::RealSection>& section = field.section();
+ CPPUNIT_ASSERT(!section.isNull());
+ for (SubMesh::SieveMesh::label_sequence::iterator v_iter=vertices->begin();
+ v_iter != vertices->end();
+ ++v_iter) {
+ section->restrictPoint(*v_iter, &values[0], fiberDim);
+ for (int iDim=0; iDim < fiberDim; ++iDim)
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valuesE[i++], values[iDim], tolerance);
+ } // for
+} // testScatterVectorToSection
+
+// ----------------------------------------------------------------------
+void
+pylith::topology::TestMultiFieldSubMesh::_buildMesh(Mesh* mesh,
+ SubMesh* submesh)
+{ // _buildMesh
+ assert(0 != mesh);
+ assert(0 != submesh);
+
+ mesh->createSieveMesh(_TestMultiFieldSubMesh::cellDim);
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = mesh->sieveMesh();
+ CPPUNIT_ASSERT(!sieveMesh.isNull());
+
+ ALE::Obj<Mesh::SieveMesh::sieve_type> sieve =
+ new Mesh::SieveMesh::sieve_type(sieveMesh->comm());
+ CPPUNIT_ASSERT(!sieve.isNull());
+
+ ALE::Obj<SieveFlexMesh::sieve_type> s =
+ new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
+ CPPUNIT_ASSERT(!s.isNull());
+
+ const int cellDim = _TestMultiFieldSubMesh::cellDim;
+ const int ncells = _TestMultiFieldSubMesh::ncells;
+ const int* cells = _TestMultiFieldSubMesh::cells;
+ const int nvertices = _TestMultiFieldSubMesh::nvertices;
+ const int ncorners = _TestMultiFieldSubMesh::ncorners;
+ const int spaceDim = _TestMultiFieldSubMesh::cellDim;
+ const double* coordinates = _TestMultiFieldSubMesh::coordinates;
+ const bool interpolate = false;
+ ALE::SieveBuilder<SieveFlexMesh>::buildTopology(s, cellDim, ncells, (int*) cells,
+ nvertices, interpolate,
+ ncorners);
+ std::map<Mesh::SieveMesh::point_type,Mesh::SieveMesh::point_type> renumbering;
+ ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
+ sieveMesh->setSieve(sieve);
+ sieveMesh->stratify();
+ ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
+ coordinates);
+
+ typedef Mesh::SieveMesh::int_section_type::chart_type chart_type;
+ const ALE::Obj<SieveMesh::int_section_type>& groupField =
+ sieveMesh->getIntSection(_TestMultiFieldSubMesh::label);
+ assert(!groupField.isNull());
+
+ const int numPoints = _TestMultiFieldSubMesh::groupSize;
+ const int numVertices = sieveMesh->depthStratum(0)->size();
+ const int numCells = sieveMesh->heightStratum(0)->size();
+ groupField->setChart(chart_type(numCells, numCells+numVertices));
+ for(int i=0; i < numPoints; ++i)
+ groupField->setFiberDimension(numCells+_TestMultiFieldSubMesh::groupVertices[i],
+ 1);
+ sieveMesh->allocate(groupField);
+
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim(spaceDim);
+ cs.initialize();
+ mesh->coordsys(&cs);
+
+ submesh->createSubMesh(*mesh, _TestMultiFieldSubMesh::label);
+} // _buildMesh
+
+
+// End of file
Added: short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldSubMesh.hh
===================================================================
--- short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldSubMesh.hh (rev 0)
+++ short/3D/PyLith/branches/multifields/unittests/libtests/topology/TestMultiFieldSubMesh.hh 2011-05-21 19:47:16 UTC (rev 18415)
@@ -0,0 +1,160 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/topology/TestMultiFieldSubMesh.hh
+ *
+ * @brief C++ unit testing for Field.
+ */
+
+#if !defined(pylith_topology_testmultifieldsubmesh_hh)
+#define pylith_topology_testmultifieldsubmesh_hh
+
+// Include directives ---------------------------------------------------
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "pylith/topology/topologyfwd.hh" // forward declarations
+
+// Forward declarations -------------------------------------------------
+/// Namespace for pylith package
+namespace pylith {
+ namespace topology {
+ class TestMultiFieldSubMesh;
+ } // topology
+} // pylith
+
+// TestMultiFieldSubMesh -----------------------------------------------------
+/// C++ unit testing for Field.
+class pylith::topology::TestMultiFieldSubMesh : public CppUnit::TestFixture
+{ // class TestMultiFieldSubMesh
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestMultiFieldSubMesh );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testSection );
+ CPPUNIT_TEST( testMesh );
+ CPPUNIT_TEST( testSpaceDim );
+ CPPUNIT_TEST( testNewSection );
+ CPPUNIT_TEST( testNewSectionPoints );
+ CPPUNIT_TEST( testNewSectionDomain );
+ CPPUNIT_TEST( testNewSectionField );
+ CPPUNIT_TEST( testCloneSection );
+ CPPUNIT_TEST( testClear );
+ CPPUNIT_TEST( testAllocate );
+ CPPUNIT_TEST( testZero );
+ CPPUNIT_TEST( testComplete );
+ CPPUNIT_TEST( testCopy );
+ CPPUNIT_TEST( testOperatorAdd );
+ CPPUNIT_TEST( testDimensionalize );
+ CPPUNIT_TEST( testView );
+ CPPUNIT_TEST( testCreateScatter );
+ CPPUNIT_TEST( testCreateScatterWithBC );
+ CPPUNIT_TEST( testVector );
+ CPPUNIT_TEST( testScatterSectionToVector );
+ CPPUNIT_TEST( testScatterVectorToSection );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test section().
+ void testSection(void);
+
+ /// Test mesh().
+ void testMesh(void);
+
+ /// Test spaceDim().
+ void testSpaceDim(void);
+
+ /// Test newSection().
+ void testNewSection(void);
+
+ /// Test newSection(points).
+ void testNewSectionPoints(void);
+
+ /// Test newSection(domain).
+ void testNewSectionDomain(void);
+
+ /// Test newSection(field).
+ void testNewSectionField(void);
+
+ /// Test cloneSection().
+ void testCloneSection(void);
+
+ /// Test clear().
+ void testClear(void);
+
+ /// Test allocate().
+ void testAllocate(void);
+
+ /// Test zero().
+ void testZero(void);
+
+ /// Test complete().
+ void testComplete(void);
+
+ /// Test copy().
+ void testCopy(void);
+
+ /// Test operator+=().
+ void testOperatorAdd(void);
+
+ /// Test dimensionalize().
+ void testDimensionalize(void);
+
+ /// Test view().
+ void testView(void);
+
+ /// Test createScatter().
+ void testCreateScatter(void);
+
+ /// Test createScatterWithBC().
+ void testCreateScatterWithBC(void);
+
+ /// Test vector().
+ void testVector(void);
+
+ /// Test scatterSectionToVector().
+ void testScatterSectionToVector(void);
+
+ /// Test scatterVectorToSection().
+ void testScatterVectorToSection(void);
+
+// PRIVATE METHODS /////////////////////////////////////////////////////
+private :
+
+ /** Build mesh.
+ *
+ * @param mesh Finite-element mesh.
+ * @param submesh Boundary mesh.
+ */
+ static
+ void _buildMesh(Mesh* mesh,
+ SubMesh* submesh);
+
+}; // class TestMultiFieldSubMesh
+
+#endif // pylith_topology_testmultifieldsubmesh_hh
+
+
+// End of file
More information about the CIG-COMMITS
mailing list