[cig-commits] r17177 - in short/3D/PyLith/trunk: libsrc/bc libsrc/faults libsrc/feassemble libsrc/friction libsrc/materials libsrc/topology unittests/libtests/bc unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Tue Sep 7 15:35:37 PDT 2010
Author: brad
Date: 2010-09-07 15:35:37 -0700 (Tue, 07 Sep 2010)
New Revision: 17177
Added:
short/3D/PyLith/trunk/libsrc/topology/UniformSectionDS.cc
short/3D/PyLith/trunk/libsrc/topology/UniformSectionDS.hh
Modified:
short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/trunk/libsrc/friction/FrictionModel.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/materials/Material.cc
short/3D/PyLith/trunk/libsrc/topology/Field.cc
short/3D/PyLith/trunk/libsrc/topology/Field.hh
short/3D/PyLith/trunk/libsrc/topology/Field.icc
short/3D/PyLith/trunk/libsrc/topology/FieldsNew.cc
short/3D/PyLith/trunk/libsrc/topology/FieldsNew.hh
short/3D/PyLith/trunk/libsrc/topology/Makefile.am
short/3D/PyLith/trunk/libsrc/topology/Mesh.hh
short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh
short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
Log:
Adjusted Visitor typedefs (moved to Field). Added UniformSection (will eventually move to Sieve).
Modified: short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/bc/AbsorbingDampers.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -43,8 +43,9 @@
typedef pylith::topology::SubMesh::RealSection SubRealSection;
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::UpdateAddVisitor UpdateAddVisitor;
+typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveMesh::order_type,PetscInt> IndicesVisitor;
// ----------------------------------------------------------------------
// Default constructor.
@@ -149,9 +150,8 @@
const ALE::Obj<RealSection>& coordinates =
sieveSubMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
@@ -588,10 +588,9 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
solutionSection);
assert(!globalOrder.isNull());
- topology::Mesh::IndicesVisitor jacobianVisitor(*solutionSection,
- *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianVisitor(*solutionSection, *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
@@ -754,8 +753,7 @@
const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
assert(!jacobianSection.isNull());
- topology::Mesh::UpdateAddVisitor jacobianVisitor(*jacobianSection,
- &_cellVector[0]);
+ UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -40,8 +40,10 @@
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
typedef pylith::topology::SubMesh::RealSection SubRealSection;
typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::UpdateAddVisitor UpdateAddVisitor;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::bc::Neumann::Neumann(void)
@@ -115,8 +117,7 @@
_parameters->get("value").section();
assert(!tractionSection.isNull());
const ALE::Obj<RealSection>& residualSection = residual.section();
- topology::SubMesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
@@ -425,9 +426,8 @@
const ALE::Obj<RealSection>& coordinates =
subSieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
const ALE::Obj<RealSection>& section = field->section();
assert(!section.isNull());
@@ -538,9 +538,8 @@
const ALE::Obj<RealSection>& coordinates =
subSieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
const ALE::Obj<RealSection>& initialSection = (0 != _dbInitial) ?
_parameters->get("initial").section() : 0;
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -54,8 +54,11 @@
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::UpdateAddVisitor UpdateAddVisitor;
+typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveSubMesh::order_type,PetscInt> IndicesVisitor;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(void) :
@@ -1126,8 +1129,8 @@
assert(!forcesInitialSection.isNull());
double_array forcesInitialCell(numBasis*spaceDim);
double_array tractionQuadPt(spaceDim);
- topology::Mesh::UpdateAddVisitor forcesInitialVisitor(*forcesInitialSection,
- &forcesInitialCell[0]);
+ UpdateAddVisitor forcesInitialVisitor(*forcesInitialSection,
+ &forcesInitialCell[0]);
assert(0 != _dbInitialTract);
_dbInitialTract->open();
@@ -1637,10 +1640,10 @@
faultSieveMesh->getFactory()->getGlobalOrder(faultSieveMesh, "default", solutionFaultSection);
assert(!globalOrderFault.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::SubMesh::IndicesVisitor jacobianFaultVisitor(*solutionFaultSection,
- *globalOrderFault,
- (int) pow(faultSieveMesh->getSieve()->getMaxConeSize(),
- faultSieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianFaultVisitor(*solutionFaultSection,
+ *globalOrderFault,
+ (int) pow(faultSieveMesh->getSieve()->getMaxConeSize(),
+ faultSieveMesh->depth())*spaceDim);
const int iCone = (negativeSide) ? 0 : 1;
for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveLagrange.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -52,7 +52,9 @@
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::UpdateAddVisitor UpdateAddVisitor;
+typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveSubMesh::order_type,PetscInt> IndicesVisitor;
// ----------------------------------------------------------------------
// Default constructor.
@@ -857,7 +859,7 @@
_fields->get("orientation").section();
assert(!orientationSection.isNull());
RestrictVisitor orientationVisitor(*orientationSection,
- orientationCell.size(),
+ orientationCell.size(),
&orientationCell[0]);
const int numConstraintVert = numBasis;
@@ -885,9 +887,9 @@
solutionSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*solutionSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianVisitor(*solutionSection, *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
@@ -1618,8 +1620,9 @@
const ALE::Obj<RealSection>& coordinatesSection =
faultSieveMesh->getRealSection("coordinates");
assert(!coordinatesSection.isNull());
- topology::Mesh::RestrictVisitor coordinatesVisitor(*coordinatesSection,
- coordinatesCell.size(), &coordinatesCell[0]);
+ RestrictVisitor coordinatesVisitor(*coordinatesSection,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
// Set orientation function
assert(cohesiveDim == _quadrature->cellDim());
@@ -1840,14 +1843,14 @@
area.zero();
const ALE::Obj<RealSection>& areaSection = area.section();
assert(!areaSection.isNull());
- topology::Mesh::UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
+ UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
double_array coordinatesCell(numBasis * spaceDim);
const ALE::Obj<RealSection>& coordinates = faultSieveMesh->getRealSection(
"coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveTract.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -27,11 +27,14 @@
// ----------------------------------------------------------------------
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
+typedef pylith::topology::SubMesh::DomainSieveMesh SieveMesh;
typedef pylith::topology::SubMesh::RealSection SubRealSection;
typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::UpdateAddVisitor UpdateAddVisitor;
+typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveSubMesh::order_type,PetscInt> IndicesVisitor;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveTract::FaultCohesiveTract(void) :
@@ -232,9 +235,8 @@
const ALE::Obj<RealSection>& coordinates =
faultSieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
// :TODO: Use spaces to create subsections like in FaultCohesiveKin.
_fields->add("orientation", "orientation",
@@ -357,9 +359,8 @@
const ALE::Obj<RealSection>& coordinates =
faultSieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -44,12 +44,6 @@
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityExplicit::ElasticityExplicit(void) :
_dtm1(-1.0)
@@ -675,17 +669,16 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solnSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -828,16 +821,14 @@
// Get sections
const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
assert(!jacobianSection.isNull());
- topology::Mesh::UpdateAddVisitor jacobianVisitor(*jacobianSection,
- &_cellVector[0]);
+ UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitLgDeform.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitLgDeform.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -43,12 +43,6 @@
//#define PRECOMPUTE_GEOMETRY
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityExplicitLgDeform::ElasticityExplicitLgDeform(void) :
_dtm1(-1.0)
@@ -194,16 +188,14 @@
RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
const ALE::Obj<RealSection>& residualSection = residual.section();
- topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
@@ -408,16 +400,14 @@
RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
const ALE::Obj<RealSection>& residualSection = residual.section();
- topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
@@ -592,17 +582,16 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solnSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
@@ -718,16 +707,14 @@
// Get sections
const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
assert(!jacobianSection.isNull());
- topology::Mesh::UpdateAddVisitor jacobianVisitor(*jacobianSection,
- &_cellVector[0]);
+ UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTet4.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -42,12 +42,6 @@
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
-
-// ----------------------------------------------------------------------
const int pylith::feassemble::ElasticityExplicitTet4::_spaceDim = 3;
const int pylith::feassemble::ElasticityExplicitTet4::_cellDim = 3;
const int pylith::feassemble::ElasticityExplicitTet4::_tensorSize = 6;
@@ -752,17 +746,16 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*dispSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianVisitor(*dispSection, *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -893,16 +886,14 @@
// Get sections
const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
assert(!jacobianSection.isNull());
- topology::Mesh::UpdateAddVisitor jacobianVisitor(*jacobianSection,
- &_cellVector[0]);
+ UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicitTri3.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -42,12 +42,6 @@
//#define DETAILED_EVENT_LOGGING
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
-
-// ----------------------------------------------------------------------
const int pylith::feassemble::ElasticityExplicitTri3::_spaceDim = 2;
const int pylith::feassemble::ElasticityExplicitTri3::_cellDim = 2;
const int pylith::feassemble::ElasticityExplicitTri3::_tensorSize = 3;
@@ -684,17 +678,16 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", dispSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*dispSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianVisitor(*dispSection, *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -825,16 +818,14 @@
// Get sections
const ALE::Obj<RealSection>& jacobianSection = jacobian->section();
assert(!jacobianSection.isNull());
- topology::Mesh::UpdateAddVisitor jacobianVisitor(*jacobianSection,
- &_cellVector[0]);
+ UpdateAddVisitor jacobianVisitor(*jacobianSection, &_cellVector[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -47,9 +47,11 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Mesh::UpdateAddVisitor UpdateAddVisitor;
+typedef pylith::topology::Field<pylith::topology::Mesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::Mesh>::UpdateAddVisitor UpdateAddVisitor;
+typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveMesh::order_type,PetscInt> IndicesVisitor;
+
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityImplicit::ElasticityImplicit(void) :
@@ -415,10 +417,10 @@
dispTSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*dispTSection,
- *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianVisitor(*dispTSection,
+ *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicitLgDeform.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicitLgDeform.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -45,10 +45,6 @@
//#define PRECOMPUTE_GEOMETRY
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityImplicitLgDeform::ElasticityImplicitLgDeform(void) :
_dtm1(-1.0)
@@ -188,27 +184,23 @@
const ALE::Obj<RealSection>& dispTSection =
fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim,
- &dispTCell[0]);
+ RestrictVisitor dispTVisitor(*dispTSection,
+ numBasis*spaceDim, &dispTCell[0]);
const ALE::Obj<RealSection>& dispTIncrSection =
fields->get("dispIncr(t->t+dt)").section();
assert(!dispTIncrSection.isNull());
- topology::Mesh::RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
- numBasis*spaceDim,
- &dispTIncrCell[0]);
+ RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ numBasis*spaceDim, &dispTIncrCell[0]);
const ALE::Obj<RealSection>& residualSection = residual.section();
- topology::Mesh::UpdateAddVisitor residualVisitor(*residualSection,
- &_cellVector[0]);
+ UpdateAddVisitor residualVisitor(*residualSection, &_cellVector[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
assert(0 != _normalizer);
const double lengthScale = _normalizer->lengthScale();
@@ -392,15 +384,13 @@
const ALE::Obj<RealSection>& dispTSection =
fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- topology::Mesh::RestrictVisitor dispTVisitor(*dispTSection,
- numBasis*spaceDim,
- &dispTCell[0]);
+ RestrictVisitor dispTVisitor(*dispTSection,
+ numBasis*spaceDim, &dispTCell[0]);
const ALE::Obj<RealSection>& dispTIncrSection =
fields->get("dispIncr(t->t+dt)").section();
assert(!dispTIncrSection.isNull());
- topology::Mesh::RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
- numBasis*spaceDim,
- &dispTIncrCell[0]);
+ RestrictVisitor dispTIncrVisitor(*dispTIncrSection,
+ numBasis*spaceDim, &dispTIncrCell[0]);
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
@@ -415,18 +405,17 @@
dispTSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- topology::Mesh::IndicesVisitor jacobianVisitor(*dispTSection,
- *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ IndicesVisitor jacobianVisitor(*dispTSection,
+ *globalOrder,
+ (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
+ sieveMesh->depth())*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
_logger->eventEnd(setupEvent);
_logger->eventBegin(computeEvent);
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -41,10 +41,6 @@
//#define PRECOMPUTE_GEOMETRY
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::IntegratorElasticity::IntegratorElasticity(void) :
_material(0),
@@ -211,17 +207,15 @@
const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
const ALE::Obj<RealSection>& dispSection = disp.section();
assert(!dispSection.isNull());
- topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
- dispCell.size(), &dispCell[0]);
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
#endif
// Loop over cells
@@ -537,17 +531,15 @@
const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
const ALE::Obj<RealSection>& dispSection = disp.section();
assert(!dispSection.isNull());
- topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
- dispCell.size(), &dispCell[0]);
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
#if !defined(PRECOMPUTE_GEOMETRY)
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
#endif
const ALE::Obj<RealSection>& fieldSection = field->section();
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.hh 2010-09-07 22:35:37 UTC (rev 17177)
@@ -270,6 +270,16 @@
/// Buffers for output.
topology::Fields<topology::Field<topology::Mesh> >* _outputFields;
+// PROTECTED TYPEDEFS ///////////////////////////////////////////////////
+protected :
+
+ typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+ typedef pylith::topology::Mesh::RealSection RealSection;
+
+ typedef pylith::topology::Field<pylith::topology::Mesh>::RestrictVisitor RestrictVisitor;
+ typedef pylith::topology::Field<pylith::topology::Mesh>::UpdateAddVisitor UpdateAddVisitor;
+ typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveMesh::order_type,PetscInt> IndicesVisitor;
+
// NOT IMPLEMENTED //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -41,10 +41,6 @@
//#define PRECOMPUTE_GEOMETRY
// ----------------------------------------------------------------------
-typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-
-// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::IntegratorElasticityLgDeform::IntegratorElasticityLgDeform(void)
{ // constructor
@@ -123,16 +119,14 @@
const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
const ALE::Obj<RealSection>& dispSection = disp.section();
assert(!dispSection.isNull());
- topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
- dispCell.size(), &dispCell[0]);
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
// Loop over cells
for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
@@ -221,16 +215,14 @@
const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
const ALE::Obj<RealSection>& dispSection = disp.section();
assert(!dispSection.isNull());
- topology::Mesh::RestrictVisitor dispVisitor(*dispSection,
- dispCell.size(), &dispCell[0]);
+ RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
const ALE::Obj<RealSection>& fieldSection = field->section();
assert(!fieldSection.isNull());
Modified: short/3D/PyLith/trunk/libsrc/friction/FrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/friction/FrictionModel.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/friction/FrictionModel.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -39,8 +39,10 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::UpdateAddVisitor UpdateAddVisitor;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::friction::FrictionModel::FrictionModel(const materials::Metadata& metadata) :
@@ -171,8 +173,7 @@
const ALE::Obj<RealSection>& propertiesSection = _properties->section();
assert(!propertiesSection.isNull());
double_array propertiesCell(numBasis*numDBProperties);
- topology::Mesh::UpdateAddVisitor propertiesVisitor(*propertiesSection,
- &propertiesCell[0]);
+ UpdateAddVisitor propertiesVisitor(*propertiesSection, &propertiesCell[0]);
// Setup database for querying for physical properties
assert(0 != _dbProperties);
@@ -287,8 +288,7 @@
const ALE::Obj<RealSection>& stateVarsSection =_stateVars->section();
assert(!stateVarsSection.isNull());
double_array stateVarsCell(numBasis*numDBStateVars);
- topology::Mesh::UpdateAddVisitor stateVarsVisitor(*stateVarsSection,
- &stateVarsCell[0]);
+ UpdateAddVisitor stateVarsVisitor(*stateVarsSection, &stateVarsCell[0]);
for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -40,6 +40,9 @@
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
+typedef pylith::topology::Field<pylith::topology::Mesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::Mesh>::UpdateAddVisitor UpdateAddVisitor;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::ElasticMaterial::ElasticMaterial(const int dimension,
@@ -323,9 +326,8 @@
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
#endif
// Create arrays for querying
@@ -470,9 +472,8 @@
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
#endif
// Create arrays for querying
Modified: short/3D/PyLith/trunk/libsrc/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Material.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/materials/Material.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -39,8 +39,10 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::Mesh>::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::Mesh>::UpdateAddVisitor UpdateAddVisitor;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::Material::Material(const int dimension,
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -32,8 +32,8 @@
// ----------------------------------------------------------------------
// Default constructor.
-template<typename mesh_type>
-pylith::topology::Field<mesh_type>::Field(const mesh_type& mesh) :
+template<typename mesh_type, typename section_type>
+pylith::topology::Field<mesh_type, section_type>::Field(const mesh_type& mesh) :
_mesh(mesh),
_vector(0),
_scatter(0),
@@ -47,9 +47,9 @@
// ----------------------------------------------------------------------
// Constructor with mesh, section, and metadata.
-template<typename mesh_type>
-pylith::topology::Field<mesh_type>::Field(const mesh_type& mesh,
- const ALE::Obj<RealSection>& section,
+template<typename mesh_type, typename section_type>
+pylith::topology::Field<mesh_type, section_type>::Field(const mesh_type& mesh,
+ const ALE::Obj<section_type>& section,
const Metadata& metadata) :
_metadata(metadata),
_mesh(mesh),
@@ -63,17 +63,17 @@
// ----------------------------------------------------------------------
// Destructor.
-template<typename mesh_type>
-pylith::topology::Field<mesh_type>::~Field(void)
+template<typename mesh_type, typename section_type>
+pylith::topology::Field<mesh_type, section_type>::~Field(void)
{ // destructor
deallocate();
} // destructor
// ----------------------------------------------------------------------
// Deallocate PETSc and local data structures.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::deallocate(void)
+pylith::topology::Field<mesh_type, section_type>::deallocate(void)
{ // deallocate
PetscErrorCode err = 0;
if (0 != _vector) {
@@ -94,9 +94,9 @@
// ----------------------------------------------------------------------
// Get spatial dimension of domain.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
int
-pylith::topology::Field<mesh_type>::spaceDim(void) const
+pylith::topology::Field<mesh_type, section_type>::spaceDim(void) const
{ // spaceDim
const spatialdata::geocoords::CoordSys* cs = _mesh.coordsys();
return (0 != cs) ? cs->spaceDim() : 0;
@@ -104,32 +104,32 @@
// ----------------------------------------------------------------------
// Get the chart size.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
int
-pylith::topology::Field<mesh_type>::chartSize(void) const
+pylith::topology::Field<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>
+template<typename mesh_type, typename section_type>
int
-pylith::topology::Field<mesh_type>::sectionSize(void) const
+pylith::topology::Field<mesh_type, section_type>::sectionSize(void) const
{ // sectionSize
return _section.isNull() ? 0 : _section->size();
} // sectionSize
// ----------------------------------------------------------------------
// Create seive section.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::newSection(void)
+pylith::topology::Field<mesh_type, section_type>::newSection(void)
{ // newSection
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
- _section = new RealSection(_mesh.comm(), _mesh.debug());
+ _section = new section_type(_mesh.comm(), _mesh.debug());
logger.stagePop();
} // newSection
@@ -137,9 +137,9 @@
// ----------------------------------------------------------------------
// Create sieve section and set chart and fiber dimesion for a
// sequence of points.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::newSection(
+pylith::topology::Field<mesh_type, section_type>::newSection(
const ALE::Obj<label_sequence>& points,
const int fiberDim)
{ // newSection
@@ -154,7 +154,7 @@
throw std::runtime_error(msg.str());
} // if
- _section = new RealSection(_mesh.comm(), _mesh.debug());
+ _section = new section_type(_mesh.comm(), _mesh.debug());
if (points->size() > 0) {
const point_type pointMin =
@@ -172,9 +172,9 @@
// ----------------------------------------------------------------------
// Create sieve section and set chart and fiber dimesion for a list of
// points.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::newSection(const int_array& points,
+pylith::topology::Field<mesh_type, section_type>::newSection(const int_array& points,
const int fiberDim)
{ // newSection
typedef typename mesh_type::SieveMesh::point_type point_type;
@@ -188,7 +188,7 @@
throw std::runtime_error(msg.str());
} // if
- _section = new RealSection(_mesh.comm(), _mesh.debug());
+ _section = new section_type(_mesh.comm(), _mesh.debug());
const int npts = points.size();
if (npts > 0) {
@@ -205,9 +205,9 @@
// ----------------------------------------------------------------------
// Create sieve section and set chart and fiber dimesion.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::newSection(const DomainEnum domain,
+pylith::topology::Field<mesh_type, section_type>::newSection(const DomainEnum domain,
const int fiberDim,
const int stratum)
{ // newSection
@@ -230,9 +230,9 @@
// ----------------------------------------------------------------------
// Create section given chart.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::newSection(const Field& src,
+pylith::topology::Field<mesh_type, section_type>::newSection(const Field& src,
const int fiberDim)
{ // newSection
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
@@ -250,7 +250,7 @@
throw std::runtime_error(msg.str());
} // if
- const ALE::Obj<RealSection>& srcSection = src.section();
+ const ALE::Obj<section_type>& srcSection = src.section();
if (!srcSection.isNull()) {
_section->setChart(srcSection->getChart());
const chart_type& chart = _section->getChart();
@@ -269,9 +269,9 @@
// ----------------------------------------------------------------------
// Create section with same layout as another section.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::cloneSection(const Field& src)
+pylith::topology::Field<mesh_type, section_type>::cloneSection(const Field& src)
{ // cloneSection
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
@@ -281,7 +281,7 @@
_metadata = src._metadata;
_metadata.label = label;
- const ALE::Obj<RealSection>& srcSection = src.section();
+ const ALE::Obj<section_type>& srcSection = src.section();
if (!srcSection.isNull() && _section.isNull()) {
logger.stagePop();
newSection();
@@ -313,9 +313,9 @@
// ----------------------------------------------------------------------
// Clear variables associated with section.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::clear(void)
+pylith::topology::Field<mesh_type, section_type>::clear(void)
{ // clear
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
@@ -333,9 +333,9 @@
// ----------------------------------------------------------------------
// Allocate Sieve section.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::allocate(void)
+pylith::topology::Field<mesh_type, section_type>::allocate(void)
{ // allocate
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
@@ -351,9 +351,9 @@
// ----------------------------------------------------------------------
// Zero section values (excluding constrained DOF).
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::zero(void)
+pylith::topology::Field<mesh_type, section_type>::zero(void)
{ // zero
if (!_section.isNull())
_section->zero(); // Does not zero BC.
@@ -361,9 +361,9 @@
// ----------------------------------------------------------------------
// Zero section values (including constrained DOF).
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::zeroAll(void)
+pylith::topology::Field<mesh_type, section_type>::zeroAll(void)
{ // zeroAll
if (!_section.isNull()) {
const chart_type& chart = _section->getChart();
@@ -389,9 +389,9 @@
// ----------------------------------------------------------------------
// Complete section by assembling across processors.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::complete(void)
+pylith::topology::Field<mesh_type, section_type>::complete(void)
{ // complete
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Completion");
@@ -409,9 +409,9 @@
// ----------------------------------------------------------------------
// Copy field values and metadata.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::copy(const Field& field)
+pylith::topology::Field<mesh_type, section_type>::copy(const Field& field)
{ // copy
// Check compatibility of sections
const int srcSize = field.chartSize();
@@ -460,9 +460,9 @@
// ----------------------------------------------------------------------
// Copy field values.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::copy(const ALE::Obj<typename mesh_type::RealSection>& osection)
+pylith::topology::Field<mesh_type, section_type>::copy(const ALE::Obj<section_type>& osection)
{ // copy
// Check compatibility of sections
const int srcSize = osection->getChart().size();
@@ -501,9 +501,9 @@
// ----------------------------------------------------------------------
// Add two fields, storing the result in one of the fields.
-template<typename mesh_type>
-pylith::topology::Field<mesh_type>&
-pylith::topology::Field<mesh_type>::operator+=(const Field& field)
+template<typename mesh_type, typename section_type>
+pylith::topology::Field<mesh_type, section_type>&
+pylith::topology::Field<mesh_type, section_type>::operator+=(const Field& field)
{ // operator+=
// Check compatibility of sections
const int srcSize = field.chartSize();
@@ -560,9 +560,9 @@
// ----------------------------------------------------------------------
// Dimensionalize field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::dimensionalize(void) const
+pylith::topology::Field<mesh_type, section_type>::dimensionalize(void) const
{ // dimensionalize
if (!_metadata.dimsOkay) {
std::ostringstream msg;
@@ -598,9 +598,9 @@
// ----------------------------------------------------------------------
// Print field to standard out.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::view(const char* label) const
+pylith::topology::Field<mesh_type, section_type>::view(const char* label) const
{ // view
std::string vecFieldString;
switch(_metadata.vectorFieldType)
@@ -646,9 +646,9 @@
// ----------------------------------------------------------------------
// Create PETSc vector for field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::createVector(void)
+pylith::topology::Field<mesh_type, section_type>::createVector(void)
{ // createVector
PetscErrorCode err = 0;
@@ -680,9 +680,9 @@
// 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>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::createScatter(void)
+pylith::topology::Field<mesh_type, section_type>::createScatter(void)
{ // createScatter
assert(!_section.isNull());
assert(!_mesh.sieveMesh().isNull());
@@ -708,9 +708,9 @@
// ----------------------------------------------------------------------
// Scatter section information across processors to update the
// PETSc vector view of the field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::scatterSectionToVector(void) const
+pylith::topology::Field<mesh_type, section_type>::scatterSectionToVector(void) const
{ // scatterSectionToVector
scatterSectionToVector(_vector);
} // scatterSectionToVector
@@ -718,9 +718,9 @@
// ----------------------------------------------------------------------
// Scatter section information across processors to update the
// PETSc vector view of the field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::scatterSectionToVector(const PetscVec vector) const
+pylith::topology::Field<mesh_type, section_type>::scatterSectionToVector(const PetscVec vector) const
{ // scatterSectionToVector
assert(!_section.isNull());
assert(0 != _scatter);
@@ -737,9 +737,9 @@
// ----------------------------------------------------------------------
// Scatter PETSc vector information across processors to update the
// section view of the field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::scatterVectorToSection(void) const
+pylith::topology::Field<mesh_type, section_type>::scatterVectorToSection(void) const
{ // scatterVectorToSection
scatterVectorToSection(_vector);
} // scatterVectorToSection
@@ -747,9 +747,9 @@
// ----------------------------------------------------------------------
// Scatter PETSc vector information across processors to update the
// section view of the field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::scatterVectorToSection(const PetscVec vector) const
+pylith::topology::Field<mesh_type, section_type>::scatterVectorToSection(const PetscVec vector) const
{ // scatterVectorToSection
assert(!_section.isNull());
assert(0 != _scatter);
@@ -765,9 +765,9 @@
// ----------------------------------------------------------------------
// Setup split field with all one space per spatial dimension.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
void
-pylith::topology::Field<mesh_type>::splitDefault(void)
+pylith::topology::Field<mesh_type, section_type>::splitDefault(void)
{ // splitDefault
assert(!_section.isNull());
const int spaceDim = _mesh.dimension();
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.hh 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.hh 2010-09-07 22:35:37 UTC (rev 17177)
@@ -40,7 +40,8 @@
*
* Extends Sieve real general section by adding metadata.
*/
-template<typename mesh_type>
+template<typename mesh_type,
+ typename section_type>
class pylith::topology::Field : public FieldBase
{ // Field
friend class TestFieldMesh; // unit testing
@@ -52,14 +53,17 @@
// 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::RealSection RealSection;
typedef typename mesh_type::SieveMesh SieveMesh;
typedef typename SieveMesh::label_sequence label_sequence;
- typedef typename RealSection::chart_type chart_type;
+ typedef typename section_type::chart_type chart_type;
// PUBLIC MEMBERS ///////////////////////////////////////////////////////
public :
@@ -75,7 +79,7 @@
* @param mesh Finite-element mesh.
*/
Field(const mesh_type& mesh,
- const ALE::Obj<RealSection>& section,
+ const ALE::Obj<section_type>& section,
const Metadata& metadata);
/// Destructor.
@@ -88,7 +92,7 @@
*
* @returns Sieve section.
*/
- const ALE::Obj<RealSection>& section(void) const;
+ const ALE::Obj<section_type>& section(void) const;
/** Get mesh associated with field.
*
@@ -235,7 +239,7 @@
*
* @param field Field to copy.
*/
- void copy(const ALE::Obj<typename mesh_type::RealSection>& field);
+ void copy(const ALE::Obj<section_type>& field);
/** Add two fields, storing the result in one of the fields.
*
@@ -304,7 +308,7 @@
Metadata _metadata;
const mesh_type& _mesh; ///< Mesh associated with section.
- ALE::Obj<RealSection> _section; ///< Real section with data.
+ ALE::Obj<section_type> _section; ///< Real section with data.
PetscVec _vector; ///< PETSc vector associated with field.
PetscVecScatter _scatter; ///< PETSc scatter associated with field.
PetscVec _scatterVec; ///< PETSC vector used in scattering.
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.icc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.icc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -21,99 +21,99 @@
#else
// Get Sieve section.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
-const ALE::Obj<typename mesh_type::RealSection>&
-pylith::topology::Field<mesh_type>::section(void) const {
+const ALE::Obj<section_type>&
+pylith::topology::Field<mesh_type, section_type>::section(void) const {
return _section;
}
// Get mesh associated with field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
const
mesh_type&
-pylith::topology::Field<mesh_type>::mesh(void) const {
+pylith::topology::Field<mesh_type, section_type>::mesh(void) const {
return _mesh;
}
// Set label for field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
void
-pylith::topology::Field<mesh_type>::label(const char* value) {
+pylith::topology::Field<mesh_type, section_type>::label(const char* value) {
_metadata.label = value;
}
// Get label for field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
const char*
-pylith::topology::Field<mesh_type>::label(void) const {
+pylith::topology::Field<mesh_type, section_type>::label(void) const {
return _metadata.label.c_str();
}
// Set vector field type
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
void
-pylith::topology::Field<mesh_type>::vectorFieldType(const VectorFieldEnum value) {
+pylith::topology::Field<mesh_type, section_type>::vectorFieldType(const VectorFieldEnum value) {
_metadata.vectorFieldType = value;
}
// Get vector field type
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
-typename pylith::topology::Field<mesh_type>::VectorFieldEnum
-pylith::topology::Field<mesh_type>::vectorFieldType(void) const {
+typename pylith::topology::Field<mesh_type, section_type>::VectorFieldEnum
+pylith::topology::Field<mesh_type, section_type>::vectorFieldType(void) const {
return _metadata.vectorFieldType;
}
// Set scale for dimensionalizing field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
void
-pylith::topology::Field<mesh_type>::scale(const double value) {
+pylith::topology::Field<mesh_type, section_type>::scale(const double value) {
_metadata.scale = value;
}
// Get scale for dimensionalizing field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
double
-pylith::topology::Field<mesh_type>::scale(void) const {
+pylith::topology::Field<mesh_type, section_type>::scale(void) const {
return _metadata.scale;
}
// Set flag indicating whether it is okay to dimensionalize field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
void
-pylith::topology::Field<mesh_type>::addDimensionOkay(const bool value) {
+pylith::topology::Field<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>
+template<typename mesh_type, typename section_type>
inline
bool
-pylith::topology::Field<mesh_type>::addDimensionOkay(void) const {
+pylith::topology::Field<mesh_type, section_type>::addDimensionOkay(void) const {
return _metadata.dimsOkay;
}
// Get PETSc vector associated with field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
PetscVec
-pylith::topology::Field<mesh_type>::vector(void) {
+pylith::topology::Field<mesh_type, section_type>::vector(void) {
return _vector;
}
// Get PETSc vector associated with field.
-template<typename mesh_type>
+template<typename mesh_type, typename section_type>
inline
const PetscVec
-pylith::topology::Field<mesh_type>::vector(void) const {
+pylith::topology::Field<mesh_type, section_type>::vector(void) const {
return _vector;
}
Modified: short/3D/PyLith/trunk/libsrc/topology/FieldsNew.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/FieldsNew.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/FieldsNew.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -87,7 +87,6 @@
pylith::topology::FieldsNew<mesh_type>::allocate(const ALE::Obj<typename mesh_type::SieveMesh::label_sequence>& points)
{ // allocate
typedef typename mesh_type::SieveMesh::point_type point_type;
- typedef typename mesh_type::RealSection RealSection;
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Fields");
@@ -95,7 +94,7 @@
// Set fiber dimension
const int fiberDim = _fiberDim();
assert(fiberDim > 0);
- _section = new RealSection(_mesh.comm(), _mesh.debug());
+ _section = new section_type(_mesh.comm(), fiberDim, _mesh.debug());
assert(!_section.isNull());
// Set spaces
@@ -110,7 +109,7 @@
*std::min_element(points->begin(), points->end());
const point_type pointMax =
*std::max_element(points->begin(), points->end());
- _section->setChart(typename RealSection::chart_type(pointMin, pointMax+1));
+ _section->setChart(typename section_type::chart_type(pointMin, pointMax+1));
_section->setFiberDimension(points, fiberDim);
int fibration = 0;
@@ -119,7 +118,7 @@
++f_iter, ++fibration)
_section->setFiberDimension(points, f_iter->second.fiberDim, fibration);
} else // Create empty chart
- _section->setChart(typename RealSection::chart_type(0, 0));
+ _section->setChart(typename section_type::chart_type(0, 0));
// Allocate section
const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = _mesh.sieveMesh();
@@ -136,7 +135,6 @@
pylith::topology::FieldsNew<mesh_type>::allocate(const int_array& points)
{ // allocate
typedef typename mesh_type::SieveMesh::point_type point_type;
- typedef typename mesh_type::RealSection RealSection;
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
@@ -144,7 +142,7 @@
// Set fiber dimension
const int fiberDim = _fiberDim();
assert(fiberDim > 0);
- _section = new RealSection(_mesh.comm(), _mesh.debug());
+ _section = new section_type(_mesh.comm(), fiberDim, _mesh.debug());
assert(!_section.isNull());
// Set spaces
@@ -158,7 +156,7 @@
if (npts > 0) {
const point_type pointMin = points.min();
const point_type pointMax = points.max();
- _section->setChart(typename RealSection::chart_type(pointMin, pointMax+1));
+ _section->setChart(typename section_type::chart_type(pointMin, pointMax+1));
for (int i=0; i < npts; ++i)
_section->setFiberDimension(points[i], fiberDim);
@@ -170,7 +168,7 @@
_section->setFiberDimension(points[i],
f_iter->second.fiberDim, fibration);
} else // create empty chart
- _section->setChart(typename RealSection::chart_type(0, 0));
+ _section->setChart(typename section_type::chart_type(0, 0));
// Allocate section
const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = _mesh.sieveMesh();
Modified: short/3D/PyLith/trunk/libsrc/topology/FieldsNew.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/FieldsNew.hh 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/FieldsNew.hh 2010-09-07 22:35:37 UTC (rev 17177)
@@ -2,11 +2,17 @@
//
// ======================================================================
//
-// Brad T. Aagaard
-// U.S. Geological Survey
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
//
-// {LicenseText}
+// 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.
+//
// ======================================================================
//
@@ -134,6 +140,8 @@
typedef std::map< std::string, FieldInfo > map_type;
+ typedef typename mesh_type::RealSection section_type;
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
@@ -147,7 +155,7 @@
protected :
map_type _fields; ///< Fields without constraints over a common set of points.
- ALE::Obj<typename mesh_type::RealSection> _section; ///< Section containing fields.
+ ALE::Obj<section_type> _section; ///< Section containing fields.
const mesh_type& _mesh; ///< Mesh associated with fields.
// NOT IMPLEMENTED //////////////////////////////////////////////////////
Modified: short/3D/PyLith/trunk/libsrc/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/Makefile.am 2010-09-07 22:35:37 UTC (rev 17177)
@@ -40,6 +40,8 @@
SubMesh.hh \
SubMesh.icc \
SubMesh.cc \
+ UniformSectionDS.hh \
+ UniformSectionDS.cc \
topologyfwd.hh
noinst_HEADERS =
Modified: short/3D/PyLith/trunk/libsrc/topology/Mesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Mesh.hh 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/Mesh.hh 2010-09-07 22:35:37 UTC (rev 17177)
@@ -30,6 +30,8 @@
#include "spatialdata/geocoords/geocoordsfwd.hh" // forward declarations
#include "spatialdata/units/unitsfwd.hh" // forward declarations
+#include "UniformSectionDS.hh" // USES IUniformSectionDS
+
#include <petscmesh.hh> // HASA ALE::IMesh
// Mesh -----------------------------------------------------------------
@@ -57,13 +59,11 @@
*/
//@{
typedef ALE::IMesh<> SieveMesh;
+ typedef ALE::IMesh<ALE::LabelSifter<int, SieveMesh::point_type> > SieveSubMesh;
+
+ typedef SieveMesh::int_section_type IntSection;
typedef SieveMesh::real_section_type RealSection;
- typedef SieveMesh::int_section_type IntSection;
- typedef ALE::IMesh<ALE::LabelSifter<int, SieveMesh::point_type> > SieveSubMesh;
- typedef ALE::ISieveVisitor::RestrictVisitor<RealSection> RestrictVisitor;
- typedef ALE::ISieveVisitor::UpdateAddVisitor<RealSection> UpdateAddVisitor;
- typedef ALE::ISieveVisitor::UpdateAllVisitor<RealSection> UpdateAllVisitor;
- typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveMesh::order_type,PetscInt> IndicesVisitor;
+ typedef ALE::IUniformSectionDS<SieveMesh::point_type, double> RealUniformSection;
//@}
Modified: short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/SubMesh.hh 2010-09-07 22:35:37 UTC (rev 17177)
@@ -45,18 +45,17 @@
// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
public:
+ // Sieve mesh for higher level domain (mesh, not submesh)
+ typedef Mesh::SieveMesh DomainSieveMesh;
+
// Typedefs for basic types associated with Sieve mesh.
// SieveMesh, RealSection, and IntSection are used in templated code.
typedef Mesh::SieveSubMesh SieveMesh;
- typedef Mesh::RealSection RealSection;
+
typedef Mesh::IntSection IntSection;
- typedef Mesh::RestrictVisitor RestrictVisitor;
- typedef Mesh::UpdateAddVisitor UpdateAddVisitor;
- typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveMesh::order_type,PetscInt> IndicesVisitor;
+ typedef Mesh::RealSection RealSection;
+ typedef Mesh::RealUniformSection RealUniformSection;
- // Sieve mesh for higher level domain (mesh, not submesh)
- typedef Mesh::SieveMesh DomainSieveMesh;
-
// PUBLIC METHODS ///////////////////////////////////////////////////////
public :
Added: short/3D/PyLith/trunk/libsrc/topology/UniformSectionDS.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/UniformSectionDS.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/UniformSectionDS.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -0,0 +1,737 @@
+// -*- 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.
+//
+// ======================================================================
+//
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::IUniformSectionDS(MPI_Comm comm,
+ const int fiberDim,
+ const int debug) :
+ ParallelObject(comm, debug)
+{ // constructor
+ if (fiberDim <= 0) {
+ std::ostringstream msg;
+ msg << "Fiber dimension '" << fiberDim << "' for section must be positive.";
+ throw ALE::Exception(msg.str().c_str());
+ } // if
+ _fiberDim = fiberDim;
+
+ atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
+ atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(comm, debug));
+ this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
+ this->_array = NULL;
+ this->_emptyValue.v = new value_type[fiberDim];
+ for(int i = 0; i < fiberDim; ++i)
+ this->_emptyValue.v[i] = value_type();
+} // constructor
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::IUniformSectionDS(MPI_Comm comm,
+ const int fiberDim,
+ const point_type& min,
+ const point_type& max,
+ const int debug) :
+ ParallelObject(comm, debug)
+{ // constructor
+ if (fiberDim <= 0) {
+ std::ostringstream msg;
+ msg << "Fiber dimension '" << fiberDim << "' for section must be positive.";
+ throw ALE::Exception(msg.str().c_str());
+ } // if
+ _fiberDim = fiberDim;
+
+ atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
+ atlas_alloc_type(this->_allocator).construct(pAtlas,
+ atlas_type(comm, min, max, fiberDim, debug));
+ this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
+ this->_array = NULL;
+ this->_emptyValue.v = new value_type[fiberDim];
+ for(int i = 0; i < fiberDim; ++i)
+ this->_emptyValue.v[i] = value_type();
+} // constructor
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::IUniformSectionDS(const Obj<atlas_type>& atlas,
+ const int fiberDim) :
+ ParallelObject(atlas->comm(), atlas->debug()), _atlas(atlas)
+{ // constructor
+ if (fiberDim <= 0) {
+ std::ostringstream msg;
+ msg << "Fiber dimension '" << fiberDim << "' for section must be positive.";
+ throw ALE::Exception(msg.str().c_str());
+ } // if
+ _fiberDim = fiberDim;
+
+ this->_atlas->update(*this->_atlas->getChart().begin(), &fiberDim);
+ this->_array = NULL;
+ this->_emptyValue.v = new value_type[fiberDim];
+ for(int i = 0; i < fiberDim; ++i)
+ this->_emptyValue.v[i] = value_type();
+} // constructor
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::~IUniformSectionDS(void)
+{ // destructor
+ delete[] this->_emptyValue.v; this->_emptyValue.v = NULL;
+
+ if (this->_array) {
+ const int chartEnd = this->getChart().max()*_fiberDim;
+ for(int i = this->getChart().min()*_fiberDim;
+ i < chartEnd;
+ ++i)
+ this->_allocator.destroy(this->_array+i);
+ this->_array += this->getChart().min()*_fiberDim;
+ this->_allocator.deallocate(this->_array, this->sizeWithBC());
+ this->_array = NULL;
+ this->_atlas = NULL;
+ } // if
+} // destructor
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+typename ALE::IUniformSectionDS<point_type, value_type, alloc_type>::value_type*
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getRawArray(const int size)
+{ // getRawArray
+ static value_type* array = NULL;
+ static int maxSize = 0;
+
+ if (size > maxSize) {
+ const value_type dummy(0);
+
+ if (array) {
+ for(int i = 0; i < maxSize; ++i)
+ this->_allocator.destroy(array+i);
+ this->_allocator.deallocate(array, maxSize);
+ } // if
+ maxSize = size;
+ array = this->_allocator.allocate(maxSize);
+ for (int i = 0; i < maxSize; ++i)
+ this->_allocator.construct(array+i, dummy);
+ }
+ return array;
+} // getRawArray
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+bool
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::hasPoint(const point_type& point) const
+{ // hasPoint
+ return this->_atlas->hasPoint(point);
+} // hasPoint
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::checkDimension(const int& dim)
+{ // checkDimension
+ if (dim != _fiberDim) {
+ ostringstream msg;
+ msg << "Invalid fiber dimension '" << dim << "' must be '" << _fiberDim
+ << "'." << std::endl;
+ throw ALE::Exception(msg.str().c_str());
+ } // if
+} // checkDimension
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+const typename ALE::IUniformSectionDS<point_type, value_type, alloc_type>::chart_type&
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getChart(void) const
+{ // getChart
+ return this->_atlas->getChart();
+} // getChart
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::setChart(const chart_type& chart)
+{ // setChart
+ this->_atlas->setChart(chart);
+ int dim = _fiberDim;
+ this->_atlas->updatePoint(*this->getChart().begin(), &dim);
+} // setChart
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+bool
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::resizeChart(const chart_type& chart)
+{ // resizeChart
+ if ((chart.min() >= this->getChart().min()) &&
+ (chart.max() <= this->getChart().max()))
+ return false;
+ this->setChart(chart);
+ return true;
+} // resizeChart
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+const ALE::Obj<typename ALE::IUniformSectionDS<point_type, value_type, alloc_type>::atlas_type>&
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getAtlas(void) const
+{ // getAtlas
+ return this->_atlas;
+} // getAtlas
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::setAtlas(const Obj<atlas_type>& atlas)
+{ // setAtlas
+ this->_atlas = atlas;
+} // setAtlas
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::addPoint(const point_type& point)
+{ // addPoint
+ this->setFiberDimension(point, _fiberDim);
+} // addPoint
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+template<typename Points>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::addPoint(const Obj<Points>& points) {
+ const typename Points::const_iterator pointsEnd = points->end();
+ for(typename Points::iterator p_iter=points->begin();
+ p_iter != pointsEnd;
+ ++p_iter)
+ this->setFiberDimension(*p_iter, _fiberDim);
+} // addPoint
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::copy(const Obj<IUniformSectionDS>& section)
+{ // copy
+ this->getAtlas()->copy(section->getAtlas());
+ 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)
+ this->updatePoint(*c_iter, section->restrictPoint(*c_iter));
+} // copy
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+const typename ALE::IUniformSectionDS<point_type, value_type, alloc_type>::value_type*
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getDefault(void) const
+{ // getDefault
+ return this->_emptyValue.v;
+} // getDefault
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::setDefault(const value_type v[])
+{ // setDefault
+ for (int i = 0; i < _fiberDim; ++i)
+ this->_emptyValue.v[i] = v[i];
+} // setDefault
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::clear(void)
+{ // clear
+ this->zero();
+ this->_atlas->clear();
+} // clear
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+int
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getFiberDimension(const point_type& p) const
+{ // getFiberDimension
+ return this->_atlas->restrictPoint(p)[0];
+} // getFiberDimension
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::setFiberDimension(const point_type& p,
+ int dim)
+{ // setFiberDimension
+ this->checkDimension(dim);
+ this->_atlas->addPoint(p);
+ this->_atlas->updatePoint(p, &dim);
+} // setFiberDimension
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+template<typename Sequence>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::setFiberDimension(const Obj<Sequence>& points,
+ int dim)
+{ // setFiberDimension
+ const typename Sequence::const_iterator pointsEnd = points->end();
+ for (typename Sequence::iterator p_iter=points->begin();
+ p_iter != pointsEnd;
+ ++p_iter)
+ this->setFiberDimension(*p_iter, dim);
+} // setFiberDimension
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::setFiberDimension(const std::set<point_type>& points,
+ int dim)
+{ // setFiberDimension
+ const typename std::set<point_type>::const_iterator pointsEnd = points.end();
+ for (typename std::set<point_type>::iterator p_iter=points.begin();
+ p_iter != pointsEnd;
+ ++p_iter)
+ this->setFiberDimension(*p_iter, dim);
+} // setFiberDimension
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::addFiberDimension(const point_type& p,
+ int dim)
+{ // addFiberDimension
+ if (this->hasPoint(p)) {
+ ostringstream msg;
+ msg << "Invalid addition to fiber dimension " << dim
+ << " cannot exceed " << _fiberDim << std::endl;
+ throw ALE::Exception(msg.str().c_str());
+ } else
+ this->setFiberDimension(p, dim);
+} // addFiberDimension
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+int
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::size(void) const
+{ // size
+ return this->_atlas->getChart().size()*_fiberDim;
+} // size
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+int
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::sizeWithBC(void) const
+{ // sizeWithBC
+ return this->size();
+} // sizeWithBC
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::allocatePoint(void)
+{ // allocatePoint
+ this->_array = this->_allocator.allocate(this->sizeWithBC());
+ this->_array -= this->getChart().min()*_fiberDim;
+ const index_type chartEnd = this->getChart().max()*_fiberDim;
+ for(index_type i=this->getChart().min()*_fiberDim;
+ i < chartEnd;
+ ++i)
+ this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);
+} // allocatePoint
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+bool
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::reallocatePoint(const chart_type& chart,
+ values_type* oldData)
+{ // reallocatePoint
+ const chart_type oldChart = this->getChart();
+ const int oldSize = this->sizeWithBC();
+ values_type oldArray = this->_array;
+ if (!this->resizeChart(chart))
+ return false;
+ const int size = this->sizeWithBC();
+
+ this->_array = this->_allocator.allocate(size);
+ this->_array -= this->getChart().min()*_fiberDim;
+ const index_type chartEnd = this->getChart().max()*_fiberDim;
+ for(index_type i=this->getChart().min()*_fiberDim;
+ i < chartEnd;
+ ++i)
+ this->_allocator.construct(this->_array+i, this->_emptyValue.v[0]);
+
+ const index_type oldChartEnd = oldChart.max()*_fiberDim;
+ for(index_type i=oldChart.min()*_fiberDim; i < oldChartEnd; ++i)
+ this->_array[i] = oldArray[i];
+ if (!oldData) {
+ for(index_type i=oldChart.min()*_fiberDim;
+ i < oldChartEnd;
+ ++i)
+ this->_allocator.destroy(oldArray+i);
+ oldArray += this->getChart().min()*_fiberDim;
+ this->_allocator.deallocate(oldArray, oldSize);
+ ///std::cout << "Freed IUniformSection data" << std::endl;
+ } else {
+ ///std::cout << "Did not free IUniformSection data" << std::endl;
+ *oldData = oldArray;
+ } // if/else
+ return true;
+} // reallocatePoint
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+template<typename Iterator, typename Extractor>
+bool
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::reallocatePoint(const Iterator& begin,
+ const Iterator& end,
+ const Extractor& extractor)
+{ // reallocatePoint
+ point_type min = this->getChart().min();
+ point_type max = this->getChart().max()-1;
+
+ for(Iterator p_iter = begin; p_iter != end; ++p_iter) {
+ min = std::min(extractor(*p_iter), min);
+ max = std::max(extractor(*p_iter), max);
+ } // for
+ return reallocatePoint(chart_type(min, max+1));
+} // reallocatePoint
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::zero(void)
+{ // zero
+ memset(this->_array+(this->getChart().min()*_fiberDim), 0,
+ this->sizeWithBC()*sizeof(value_type));
+} // zero
+
+// ----------------------------------------------------------------------
+// Return a pointer to the entire contiguous storage array
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+const typename ALE::IUniformSectionDS<point_type, value_type, alloc_type>::values_type&
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::restrictSpace(void) const
+{ // restrictSpace
+ return this->_array;
+} // restrictSpace
+
+// ----------------------------------------------------------------------
+// Return only the values associated to this point, not its closure
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+const typename ALE::IUniformSectionDS<point_type, value_type, alloc_type>::value_type*
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::restrictPoint(const point_type& p) const
+{ // restrictPoint
+ if (!this->hasPoint(p))
+ return this->_emptyValue.v;
+ return &this->_array[p*_fiberDim];
+} // restrictPoint
+
+// ----------------------------------------------------------------------
+// Update only the values associated to this point, not its closure
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::updatePoint(const point_type& p,
+ const value_type v[])
+{ // updatePoint
+ for(int i = 0, idx = p*_fiberDim; i < _fiberDim; ++i, ++idx)
+ this->_array[idx] = v[i];
+} // updatePoint
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+// Update only the values associated to this point, not its closure
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::updateAddPoint(const point_type& p,
+ const value_type v[])
+{ // updateAddPoint
+ for(int i = 0, idx = p*_fiberDim; i < _fiberDim; ++i, ++idx)
+ this->_array[idx] += v[i];
+} // updateAddPoint
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::updatePointAll(const point_type& p,
+ const value_type v[])
+{ // updatePointAll
+ this->updatePoint(p, v);
+} // updatePointAll
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::view(const std::string& name,
+ MPI_Comm comm) const
+{ // view
+ ostringstream txt;
+ int rank;
+
+ if (comm == MPI_COMM_NULL) {
+ comm = this->comm();
+ rank = this->commRank();
+ } else {
+ MPI_Comm_rank(comm, &rank);
+ }
+ if (name == "") {
+ if(rank == 0) {
+ txt << "viewing an IUniformSection" << std::endl;
+ }
+ } else {
+ if(rank == 0) {
+ txt << "viewing IUniformSection '" << name << "'" << std::endl;
+ }
+ }
+ const typename atlas_type::chart_type& chart = this->_atlas->getChart();
+ values_type array = this->_array;
+
+ for(typename atlas_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
+ const int idx = (*p_iter)*_fiberDim;
+
+ if (_fiberDim != 0) {
+ txt << "[" << this->commRank() << "]: " << *p_iter << " dim " << _fiberDim << " ";
+ for(int i = 0; i < _fiberDim; i++) {
+ txt << " " << array[idx+i];
+ }
+ txt << std::endl;
+ }
+ }
+ if (chart.size() == 0) {
+ txt << "[" << this->commRank() << "]: empty" << std::endl;
+ }
+ PetscSynchronizedPrintf(comm, txt.str().c_str());
+ PetscSynchronizedFlush(comm);
+} // view
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+int
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getNumSpaces(void) const
+{ // getNumSpaces
+ return this->_spaces.size();
+} // getNumSpaces
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+const std::vector<ALE::Obj<typename ALE::IUniformSectionDS<point_type, value_type, alloc_type>::atlas_type> >&
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getSpaces(void)
+{ // getSpaces
+ return this->_spaces;
+} // getSpaces
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::addSpace(void)
+{ // addSpace
+ Obj<atlas_type> space = new atlas_type(this->comm(), this->debug());
+ this->_spaces.push_back(space);
+} // addSpace
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+int
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getFiberDimension(
+ const point_type& p,
+ const int space) const
+{ // getFiberDimension
+ return *this->_spaces[space]->restrictPoint(p);
+} // getFiberDimension
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::setFiberDimension(
+ const point_type& p,
+ int dim,
+ const int space)
+{ // setFiberDimension
+ const index_type idx = -1;
+ this->_spaces[space]->addPoint(p);
+ this->_spaces[space]->updatePoint(p, &idx);
+} // setFiberDimension
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+template<typename Sequence>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::setFiberDimension(
+ const Obj<Sequence>& points,
+ int dim,
+ const int space)
+{ // setFiberDimension
+ const typename Sequence::const_iterator pointsEnd = points->end();
+ for(typename Sequence::iterator p_iter = points->begin();
+ p_iter != pointsEnd;
+ ++p_iter)
+ this->setFiberDimension(*p_iter, dim, space);
+} // setFiberDimension
+
+// ----------------------------------------------------------------------
+// Return the total number of free dofs
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+int
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::size(const int space) const
+{ // size
+ const chart_type& points = this->getChart();
+ int size = 0;
+
+ const typename chart_type::const_iterator pointsEnd = points.end();
+ for (typename chart_type::const_iterator p_iter = points.begin();
+ p_iter != pointsEnd;
+ ++p_iter)
+ size += this->getConstrainedFiberDimension(*p_iter, space);
+
+ return size;
+} // size
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+template<typename OtherSection>
+void
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::copyFibration(const Obj<OtherSection>& section)
+{ // copyFibration
+ const std::vector<Obj<atlas_type> >& spaces = section->getSpaces();
+
+ this->_spaces.clear();
+ const typename std::vector<Obj<atlas_type> >::const_iteraor spacesEnd = spaces.end();
+ for(typename std::vector<Obj<atlas_type> >::const_iterator s_iter=spaces.begin();
+ s_iter != spacesEnd;
+ ++s_iter)
+ this->_spaces.push_back(*s_iter);
+} // copyFibration
+
+// ----------------------------------------------------------------------
+template<typename point_type,
+ typename value_type,
+ typename alloc_type>
+ALE::Obj<ALE::IUniformSectionDS<point_type, value_type, alloc_type> >
+ALE::IUniformSectionDS<point_type, value_type, alloc_type>::getFibration(const int space) const {
+ const chart_type& chart = this->getChart();
+
+ const int localFiberDim = this->getFiberDimension(*chart.begin(), space);
+ int fiberDim = 0;
+ MPI_Allreduce((void *) &localFiberDim, (void *) &fiberDim, 1,
+ MPI_INT, MPI_MAX, this->comm());
+ assert(fiberDim > 0);
+ Obj<IUniformSectionDS> field = new IUniformSectionDS(this->comm(),
+ fiberDim,
+ *chart.begin(),
+ *chart.end(),
+ this->debug());
+
+ // Copy sizes
+ const typename chart_type::const_iterator chartEnd = chart.end();
+ for(typename chart_type::const_iterator c_iter=chart.begin();
+ c_iter != chartEnd;
+ ++c_iter) {
+ const int fDim = this->getFiberDimension(*c_iter, space);
+ if (fDim)
+ field->setFiberDimension(*c_iter, fDim);
+ } // for
+
+ // Copy values
+ const chart_type& newChart = field->getChart();
+ const typename chart_type::const_iterator newChartEnd = newChart.end();
+ for(typename chart_type::const_iterator c_iter = newChart.begin();
+ c_iter != newChartEnd;
+ ++c_iter)
+ field->updatePoint(*c_iter, this->restrictPoint(*c_iter));
+
+ return field;
+} // getFibration
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/UniformSectionDS.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/UniformSectionDS.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/UniformSectionDS.hh 2010-09-07 22:35:37 UTC (rev 17177)
@@ -0,0 +1,197 @@
+// -*- 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/UniformSectionDS.hh
+ *
+ * @brief Sieve section with uniform size set at runtime in contrast
+ * to ALE::IUniformSection, which has a uniform size set at compile
+ * time.
+ */
+
+#if !defined(pylith_topology_uniformsection_hh)
+#define pylith_topology_uniformsection_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+#include "IField.hh"
+
+// UniformSectionDS -----------------------------------------------------
+/// Sieve section with uniform size set at runtime.
+/// All fibers are the same dimension
+/// Note we can use a IConstantSection for this Atlas
+/// Each point may have a different vector
+/// Thus we need storage for values, and hence must implement completion
+template<typename Point_,
+ typename Value_,
+ typename Alloc_ =ALE::malloc_allocator<Value_> >
+class ALE::IUniformSectionDS : public ALE::ParallelObject {
+
+// PUBLIC TYPEDEFS //////////////////////////////////////////////////////
+public:
+ typedef Point_ point_type;
+ typedef Value_ value_type;
+ typedef Alloc_ alloc_type;
+ typedef typename alloc_type::template rebind<point_type>::other point_alloc_type;
+ typedef IConstantSection<point_type, int, point_alloc_type> atlas_type;
+ typedef typename atlas_type::chart_type chart_type;
+ typedef point_type index_type;
+ typedef struct { value_type* v;} fiber_type;
+ typedef value_type* values_type;
+ typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
+ typedef typename atlas_alloc_type::pointer atlas_ptr;
+
+// PUBLIC METHODS ///////////////////////////////////////////////////////
+public:
+
+ IUniformSectionDS(MPI_Comm comm,
+ const int fiberDim,
+ const int debug =0);
+
+ IUniformSectionDS(MPI_Comm comm,
+ const int diberDim,
+ const point_type& min,
+ const point_type& max,
+ const int debug =0);
+
+ IUniformSectionDS(const Obj<atlas_type>& atlas,
+ const int fiberDim);
+
+ /// Destructor.
+ virtual ~IUniformSectionDS(void);
+
+ value_type* getRawArray(const int size);
+
+ bool hasPoint(const point_type& point) const;
+
+ void checkDimension(const int& dim);
+
+ const chart_type& getChart(void) const;
+
+ void setChart(const chart_type& chart);
+
+ bool resizeChart(const chart_type& chart);
+
+ const Obj<atlas_type>& getAtlas(void) const;
+
+ void setAtlas(const Obj<atlas_type>& atlas);
+ void addPoint(const point_type& point);
+
+ template<typename Points>
+ void addPoint(const Obj<Points>& points);
+
+ void copy(const Obj<IUniformSectionDS>& section);
+
+ const value_type* getDefault(void) const;
+
+ void setDefault(const value_type v[]);
+
+ void clear(void);
+
+ int getFiberDimension(const point_type& p) const;
+
+ void setFiberDimension(const point_type& p,
+ int dim);
+
+ template<typename Sequence>
+ void setFiberDimension(const Obj<Sequence>& points,
+ int dim);
+
+ void setFiberDimension(const std::set<point_type>& points,
+ int dim);
+
+ void addFiberDimension(const point_type& p,
+ int dim);
+
+ int size(void) const;
+
+ int sizeWithBC(void) const;
+
+ void allocatePoint(void);
+
+ bool reallocatePoint(const chart_type& chart,
+ values_type *oldData =NULL);
+
+ template<typename Iterator,
+ typename Extractor>
+ bool reallocatePoint(const Iterator& begin,
+ const Iterator& end,
+ const Extractor& extractor);
+
+ void zero(void);
+
+ // Return a pointer to the entire contiguous storage array
+ const values_type& restrictSpace(void) const;
+
+ // Return only the values associated to this point, not its closure
+ const value_type *restrictPoint(const point_type& p) const;
+
+ // Update only the values associated to this point, not its closure
+ void updatePoint(const point_type& p,
+ const value_type v[]);
+
+ // Update only the values associated to this point, not its closure
+ void updateAddPoint(const point_type& p,
+ const value_type v[]);
+
+ void updatePointAll(const point_type& p,
+ const value_type v[]);
+
+ void view(const std::string& name,
+ MPI_Comm comm =MPI_COMM_NULL) const;
+
+ int getNumSpaces(void) const;
+
+ const std::vector<Obj<atlas_type> >& getSpaces(void);
+ void addSpace(void);
+ int getFiberDimension(const point_type& p,
+ const int space) const;
+ void setFiberDimension(const point_type& p,
+ int dim,
+ const int space);
+
+ template<typename Sequence>
+ void setFiberDimension(const Obj<Sequence>& points,
+ int dim,
+ const int space);
+
+ // Return the total number of free dofs
+ int size(const int space) const;
+
+ template<typename OtherSection>
+ void copyFibration(const Obj<OtherSection>& section);
+
+ Obj<IUniformSectionDS> getFibration(const int space) const;
+
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected:
+ int _fiberDim;
+ Obj<atlas_type> _atlas;
+ std::vector<Obj<atlas_type> > _spaces;
+ values_type _array;
+ fiber_type _emptyValue;
+ alloc_type _allocator;
+
+}; // class IUniformSectionDS
+
+#include "UniformSectionDS.cc" // template definitions
+
+#endif // uniformsectionds_hh
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/libsrc/topology/topologyfwd.hh 2010-09-07 22:35:37 UTC (rev 17177)
@@ -27,6 +27,14 @@
#if !defined(pylith_topology_topologyfwd_hh)
#define pylith_topology_topologyfwd_hh
+#include "pylith/utils/sievetypes.hh"
+
+namespace ALE {
+ template<typename point_type,
+ typename value_type,
+ typename allocator> class IUniformSectionDS;
+} // ALE
+
namespace pylith {
namespace topology {
@@ -35,7 +43,8 @@
class MeshOps;
class FieldBase;
- template<typename mesh_type> class Field;
+ template<typename mesh_type,
+ typename section_type =ALE::IGeneralSection<pylith::Mesh::point_type, double> > class Field;
template<typename field_type> class Fields;
template<typename mesh_type> class FieldsNew;
class SolutionFields;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -45,7 +45,6 @@
typedef pylith::topology::SubMesh::RealSection RealSection;
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
typedef pylith::topology::SubMesh::RealSection SubRealSection;
-typedef pylith::topology::SubMesh::RestrictVisitor RestrictVisitor;
// ----------------------------------------------------------------------
// Setup testing data.
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -48,8 +48,9 @@
typedef pylith::topology::SubMesh::RealSection RealSection;
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
typedef pylith::topology::SubMesh::RealSection SubRealSection;
-typedef pylith::topology::SubMesh::RestrictVisitor RestrictVisitor;
+typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
+
// ----------------------------------------------------------------------
namespace pylith {
namespace bc {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2010-09-07 07:05:09 UTC (rev 17176)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2010-09-07 22:35:37 UTC (rev 17177)
@@ -33,6 +33,8 @@
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestQuadrature );
+typedef pylith::topology::Field<pylith::topology::Mesh>::RestrictVisitor RestrictVisitor;
+
// ----------------------------------------------------------------------
// Test copy constuctor.
void
@@ -238,9 +240,8 @@
const ALE::Obj<topology::Mesh::RealSection>& coordinates =
sieveMesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
+ RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(), &coordinatesCell[0]);
#endif
size_t size = 0;
More information about the CIG-COMMITS
mailing list