[cig-commits] r21384 - in short/3D/PyLith/trunk: libsrc/pylith/bc libsrc/pylith/faults libsrc/pylith/feassemble libsrc/pylith/friction libsrc/pylith/materials libsrc/pylith/meshio libsrc/pylith/problems libsrc/pylith/topology modulesrc/faults modulesrc/friction modulesrc/topology unittests/libtests/bc unittests/libtests/faults unittests/libtests/feassemble unittests/libtests/friction unittests/libtests/topology
knepley at geodynamics.org
knepley at geodynamics.org
Fri Feb 22 11:56:51 PST 2013
Author: knepley
Date: 2013-02-22 11:56:50 -0800 (Fri, 22 Feb 2013)
New Revision: 21384
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc
short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.hh
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh
short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Explicit.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLumped.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
short/3D/PyLith/trunk/libsrc/pylith/topology/FieldsNew.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/topologyfwd.hh
short/3D/PyLith/trunk/modulesrc/faults/TractPerturbation.i
short/3D/PyLith/trunk/modulesrc/friction/FrictionModel.i
short/3D/PyLith/trunk/modulesrc/topology/Field.i
short/3D/PyLith/trunk/modulesrc/topology/topology.i
short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc
short/3D/PyLith/trunk/unittests/libtests/topology/Makefile.am
Log:
Removed ALE::Section and FieldsNew
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -41,12 +41,7 @@
// ----------------------------------------------------------------------
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::RealUniformSection SubRealUniformSection;
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-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,PylithInt> IndicesVisitor;
// ----------------------------------------------------------------------
// Default constructor.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/BoundaryConditionPoints.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -70,29 +70,30 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("BoundaryConditions");
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
+ DM dmMesh = mesh.dmMesh();
+ DMLabel label;
+ IS pointIS;
+ const PetscInt *points;
+ PetscInt numPoints;
+ PetscBool has;
+ PetscErrorCode err;
- if (!sieveMesh->hasIntSection(_label)) {
+ assert(dmMesh);
+ err = DMPlexHasLabel(dmMesh, _label.c_str(), &has);CHECK_PETSC_ERROR(err);
+ if (!has) {
std::ostringstream msg;
msg << "Could not find group of points '" << _label << "' in mesh.";
throw std::runtime_error(msg.str());
} // if
- const ALE::Obj<SieveMesh::int_section_type>& groupField =
- sieveMesh->getIntSection(_label);
- assert(!groupField.isNull());
- const chart_type& chart = groupField->getChart();
- const chart_type::const_iterator& chartEnd = chart.end();
- const int numPoints = groupField->size();
+ err = DMPlexGetLabel(dmMesh, _label.c_str(), &label);CHECK_PETSC_ERROR(err);
+ err = DMLabelGetStratumIS(label, 1, &pointIS);CHECK_PETSC_ERROR(err);
+ err = ISGetLocalSize(pointIS, &numPoints);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(pointIS, &points);CHECK_PETSC_ERROR(err);
_points.resize(numPoints);
- int i = 0;
- for(chart_type::const_iterator c_iter = chart.begin();
- c_iter != chartEnd;
- ++c_iter)
- if (groupField->getFiberDimension(*c_iter))
- _points[i++] = *c_iter;
-
+ for(PetscInt p = 0; p < numPoints; ++p) {_points[p] = points[p];}
+ err = ISRestoreIndices(pointIS, &points);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&pointIS);CHECK_PETSC_ERROR(err);
logger.stagePop();
} // _getPoints
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/Neumann.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -41,9 +41,6 @@
typedef pylith::topology::SubMesh::RealUniformSection SubRealUniformSection;
typedef pylith::topology::Mesh::RealSection RealSection;
-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)
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -53,14 +53,8 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::Mesh::RealUniformSection RealUniformSection;
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-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,PylithInt> IndicesVisitor;
-
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(void) :
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -50,11 +50,7 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-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,PylithInt> IndicesVisitor;
// ----------------------------------------------------------------------
// Default constructor.
@@ -108,20 +104,17 @@
logger.stagePush("FaultFields");
// Allocate dispRel field
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- assert(!vertices.isNull());
_fields->add("relative disp", "relative_disp");
topology::Field<topology::SubMesh>& dispRel = _fields->get("relative disp");
- dispRel.newSection(vertices, cs->spaceDim());
+ dispRel.newSection(topology::FieldBase::VERTICES_FIELD, cs->spaceDim());
dispRel.allocate();
dispRel.vectorFieldType(topology::FieldBase::VECTOR);
dispRel.scale(_normalizer->lengthScale());
logger.stagePop();
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -28,13 +28,7 @@
// ----------------------------------------------------------------------
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::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
-typedef pylith::topology::Field<pylith::topology::SubMesh>::UpdateAddVisitor UpdateAddVisitor;
-typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveSubMesh::order_type,PylithInt> IndicesVisitor;
-
// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveTract::FaultCohesiveTract(void) :
@@ -206,14 +200,13 @@
up[i] = upDir[i];
// Get 'fault' cells.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM faultDMMesh = _faultMesh->dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(faultDMMesh);
+ err = DMPlexGetHeightStratum(faultDMMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
// Quadrature related values.
const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
@@ -234,29 +227,32 @@
// Get sections.
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMPlexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
// :TODO: Use spaces to create subsections like in FaultCohesiveKin.
_fields->add("orientation", "orientation",
topology::FieldBase::CELLS_FIELD, fiberDim);
topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
orientation.allocate();
- const ALE::Obj<RealSection>& orientationSection = orientation.section();
- assert(!orientationSection.isNull());
+ PetscSection orientationSection = orientation.petscSection();
+ Vec orientationVec = orientation.localVector();
+ PetscScalar *orientationArray;
+ assert(orientationSection);assert(orientationVec);
// Loop over cells in fault mesh and compute orientations.
- for(SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- // Compute geometry information for current cell
- coordsVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ err = VecGetArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt c = cStart; c < cEnd; ++c) {
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
+ err = DMPlexVecGetClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ err = DMPlexVecRestoreClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ _quadrature->computeGeometry(coordinatesCell, c);
+
// Reset orientation to zero.
orientationCell = 0.0;
@@ -279,9 +275,15 @@
memcpy(&orientationCell[iQuad*orientationSize],
&orientationQuadPt[0], orientationSize*sizeof(PylithScalar));
} // for
+ PetscInt dof, off;
- orientationSection->updatePoint(*c_iter, &orientationCell[0]);
+ err = PetscSectionGetDof(orientationSection, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(orientationSection, c, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ orientationArray[off+d] = orientationCell[d];
+ }
} // for
+ err = VecRestoreArray(orientationVec, &orientationArray);CHECK_PETSC_ERROR(err);
// debugging
orientation.view("FAULT ORIENTATION");
@@ -308,8 +310,10 @@
traction.allocate();
traction.scale(pressureScale);
traction.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
- const ALE::Obj<RealSection>& tractionSection = traction.section();
- assert(!tractionSection.isNull());
+ PetscSection tractionSection = traction.petscSection();
+ Vec tractionVec = traction.localVector();
+ PetscScalar *tractionArray;
+ assert(tractionSection);assert(tractionVec);
_dbInitial->open();
switch (spaceDim)
@@ -338,14 +342,13 @@
} // switch
// Get 'fault' cells.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM faultDMMesh = _faultMesh->dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+ assert(faultDMMesh);
+ err = DMPlexGetHeightStratum(faultDMMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+
const int cellDim = _quadrature->cellDim() > 0 ? _quadrature->cellDim() : 1;
const int numBasis = _quadrature->numBasis();
const int numQuadPts = _quadrature->numQuadPts();
@@ -358,24 +361,25 @@
// Get sections.
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMPlexGetCoordinateSection(faultDMMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(faultDMMesh, &coordVec);CHECK_PETSC_ERROR(err);
const spatialdata::geocoords::CoordSys* cs = _faultMesh->coordsys();
// Compute quadrature information
// Loop over cells in boundary mesh and perform queries.
- for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- // Compute geometry information for current cell
- coordsVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ err = VecGetArray(tractionVec, &tractionArray);CHECK_PETSC_ERROR(err);
+ for (PetscInt c = cStart; c < cEnd; ++c) {
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
+
+ err = DMPlexVecGetClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ err = DMPlexVecRestoreClosure(faultDMMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ _quadrature->computeGeometry(coordinatesCell, c);
const scalar_array& quadPtsNondim = _quadrature->quadPts();
quadPtsGlobal = quadPtsNondim;
@@ -383,29 +387,33 @@
lengthScale);
tractionCell = 0.0;
- for (int iQuad=0, iSpace=0;
- iQuad < numQuadPts;
- ++iQuad, iSpace+=spaceDim) {
- const int err = _dbInitial->query(&tractionCell[iQuad*spaceDim], spaceDim,
- &quadPtsGlobal[iSpace], spaceDim, cs);
- if (err) {
- std::ostringstream msg;
- msg << "Could not find initial tractions at (";
- for (int i=0; i < spaceDim; ++i)
- msg << " " << quadPtsGlobal[i+iSpace];
- msg << ") for dynamic fault interface " << label() << "\n"
- << "using spatial database " << _dbInitial->label() << ".";
- throw std::runtime_error(msg.str());
- } // if
-
+ for (int iQuad=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iSpace+=spaceDim) {
+ const int err = _dbInitial->query(&tractionCell[iQuad*spaceDim], spaceDim,
+ &quadPtsGlobal[iSpace], spaceDim, cs);
+ if (err) {
+ std::ostringstream msg;
+ msg << "Could not find initial tractions at (";
+ for (int i=0; i < spaceDim; ++i)
+ msg << " " << quadPtsGlobal[i+iSpace];
+ msg << ") for dynamic fault interface " << label() << "\n"
+ << "using spatial database " << _dbInitial->label() << ".";
+ throw std::runtime_error(msg.str());
+ } // if
} // for
_normalizer->nondimensionalize(&tractionCell[0], tractionCell.size(),
- pressureScale);
+ pressureScale);
// Update section
- assert(tractionCell.size() == tractionSection->getFiberDimension(*c_iter));
- tractionSection->updatePoint(*c_iter, &tractionCell[0]);
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(tractionSection, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(tractionSection, c, &off);CHECK_PETSC_ERROR(err);
+ assert(tractionCell.size() == dof);
+ for(PetscInt d = 0; d < dof; ++d) {
+ tractionArray[off+d] = tractionCell[d];
+ }
} // for
+ err = VecRestoreArray(tractionVec, &tractionArray);CHECK_PETSC_ERROR(err);
_dbInitial->close();
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.hh 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/TractPerturbation.hh 2013-02-22 19:56:50 UTC (rev 21384)
@@ -29,8 +29,6 @@
#include "faultsfwd.hh" // forward declarations
#include "pylith/bc/TimeDependent.hh" // ISA TimeDependent
-#include "pylith/topology/topologyfwd.hh" // HOLSA FieldsNew
-
#include "spatialdata/units/unitsfwd.hh" // USES Nondimensional
// TractPerturbation -----------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -46,12 +46,7 @@
// ----------------------------------------------------------------------
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,PylithInt> IndicesVisitor;
-
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::ElasticityImplicit::ElasticityImplicit(void) :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.hh 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.hh 2013-02-22 19:56:50 UTC (rev 21384)
@@ -241,12 +241,7 @@
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,PylithInt> IndicesVisitor;
-
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticityLgDeform.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -106,58 +106,69 @@
strainCell = 0.0;
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
+ DM dmMesh = fields->mesh().dmMesh();
+ assert(!dmMesh);
const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetLocalSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get fields
const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
- const ALE::Obj<RealSection>& dispSection = disp.section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ PetscSection dispSection = disp.petscSection();
+ Vec dispVec = disp.localVector();
+ assert(dispSection);assert(dispVec);
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for (PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Retrieve geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
+ _quadrature->retrieveGeometry(cell);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
+
+ err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ _quadrature->computeGeometry(coordinatesCell, cell);
#endif
// Restrict input fields to cell
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
+ const PetscScalar *disp = PETSC_NULL;
+ PetscInt dispSize;
+ err = DMPlexVecGetClosure(dmMesh, dispSection, dispVec, cell, &dispSize, &disp);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < dispSize; ++i) {dispCell[i] = disp[i];}
+ err = DMPlexVecRestoreClosure(dmMesh, dispSection, dispVec, cell, &dispSize, &disp);CHECK_PETSC_ERROR(err);
+
// Get cell geometry information that depends on cell
const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Compute deformation tensor.
_calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCell,
- numBasis, numQuadPts, spaceDim);
+ numBasis, numQuadPts, spaceDim);
// Compute strains
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
// Update material state
- _material->updateStateVars(strainCell, *c_iter);
+ _material->updateStateVars(strainCell, cell);
} // for
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // updateStateVars
// ----------------------------------------------------------------------
@@ -202,64 +213,82 @@
scalar_array stressCell(numQuadPts*tensorSize);
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = field->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
+ DM dmMesh = fields->mesh().dmMesh();
+ assert(!dmMesh);
const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetLocalSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get field
const topology::Field<topology::Mesh>& disp = fields->get("disp(t)");
- const ALE::Obj<RealSection>& dispSection = disp.section();
- assert(!dispSection.isNull());
- RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
+ PetscSection dispSection = disp.petscSection();
+ Vec dispVec = disp.localVector();
+ assert(dispSection);assert(dispVec);
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
- const ALE::Obj<RealSection>& fieldSection = field->section();
- assert(!fieldSection.isNull());
+ PetscSection fieldSection = field->petscSection();
+ Vec fieldVec = field->localVector();
+ PetscScalar *fieldArray;
+ assert(fieldSection);assert(fieldVec);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
+ for (PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
// Retrieve geometry information for current cell
#if defined(PRECOMPUTE_GEOMETRY)
_quadrature->retrieveGeometry(*c_iter);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
+
+ err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ _quadrature->computeGeometry(coordinatesCell, cell);
#endif
// Restrict input fields to cell
- dispVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, dispVisitor);
+ const PetscScalar *disp = PETSC_NULL;
+ PetscInt dispSize;
+ err = DMPlexVecGetClosure(dmMesh, dispSection, dispVec, cell, &dispSize, &disp);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < dispSize; ++i) {dispCell[i] = disp[i];}
+ err = DMPlexVecRestoreClosure(dmMesh, dispSection, dispVec, cell, &dispSize, &disp);CHECK_PETSC_ERROR(err);
+
// Get cell geometry information that depends on cell
const scalar_array& basisDeriv = _quadrature->basisDeriv();
// Compute deformation tensor.
_calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCell,
- numBasis, numQuadPts, spaceDim);
+ numBasis, numQuadPts, spaceDim);
// Compute strains
calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
+ PetscInt dof, off;
- if (!calcStress)
- fieldSection->updatePoint(*c_iter, &strainCell[0]);
- else {
- _material->retrievePropsAndVars(*c_iter);
+ err = PetscSectionGetDof(fieldSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, cell, &off);CHECK_PETSC_ERROR(err);
+ if (!calcStress) {
+ for (PetscInt d = 0; d < dof; ++d) {
+ fieldArray[off+d] = strainCell[d];
+ }
+ } else {
+ _material->retrievePropsAndVars(cell);
stressCell = _material->calcStress(strainCell);
- fieldSection->updatePoint(*c_iter, &stressCell[0]);
+ for (PetscInt d = 0; d < dof; ++d) {
+ fieldArray[off+d] = stressCell[d];
+ }
} // else
} // for
} // _calcStrainStressField
@@ -286,28 +315,45 @@
stressCell = 0.0;
// Get cell information
- const ALE::Obj<SieveMesh>& sieveMesh = field->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
+ DM dmMesh = field->mesh().dmMesh();
+ assert(!dmMesh);
const int materialId = _material->id();
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", materialId);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ err = DMPlexGetStratumIS(dmMesh, "material-id", materialId, &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetLocalSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Get field
- const ALE::Obj<RealSection>& fieldSection = field->section();
- assert(!fieldSection.isNull());
+ PetscSection fieldSection = field->petscSection();
+ Vec fieldVec = field->localVector();
+ PetscScalar *fieldArray;
+ assert(fieldSection);assert(fieldVec);
// Loop over cells
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- fieldSection->restrictPoint(*c_iter, &strainCell[0], strainCell.size());
- _material->retrievePropsAndVars(*c_iter);
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ for (PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
+ PetscInt dof, off;
+
+ _material->retrievePropsAndVars(cell);
+ err = PetscSectionGetDof(fieldSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(strainCell.size() == dof);
+ for (PetscInt d = 0; d < dof; ++d) {
+ strainCell[d] = fieldArray[off+d];
+ }
stressCell = _material->calcStress(strainCell);
- fieldSection->updatePoint(*c_iter, &stressCell[0]);
+ for (PetscInt d = 0; d < dof; ++d) {
+ fieldArray[off+d] = stressCell[d];
+ }
} // for
+ err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // _calcStressFromStrain
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -263,8 +263,151 @@
} // computeGeometry
// ----------------------------------------------------------------------
+// Compute geometric quantities for each cell.
template<typename mesh_type>
void
+pylith::feassemble::Quadrature<mesh_type>::computeGeometry(
+ const mesh_type& mesh,
+ PetscInt cStart, PetscInt cEnd)
+{ // computeGeometry
+ assert(0 != _engine);
+
+ const char* loggingStage = "Quadrature";
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush(loggingStage);
+
+ delete _geometryFields;
+ _geometryFields = new topology::Fields<topology::Field<mesh_type> >;
+
+ // Allocate field and cell buffer for quadrature points
+ _geometryFields->add("quadrature points", "quadrature_points");
+ topology::Field<mesh_type>& quadPtsField =
+ _geometryFields->get("quadrature points");
+ int fiberDim = _numQuadPts * _spaceDim;
+ quadPtsField.newSection(cStart, cEnd, fiberDim);
+ quadPtsField.allocate();
+
+ // Allocate field and cell buffer for Jacobian at quadrature points
+ std::cout << "Jacobian: cell dim: " << _cellDim << std::endl;
+ _geometryFields->add("jacobian", "jacobian");
+ topology::Field<mesh_type>& jacobianField =
+ _geometryFields->get("jacobian");
+ fiberDim = (_cellDim > 0) ?
+ _numQuadPts * _cellDim * _spaceDim :
+ _numQuadPts * 1 * _spaceDim;
+ jacobianField->newSection(cStart, cEnd, fiberDim);
+ jacobianField->allocate();
+
+ // Allocate field and cell buffer for determinant of Jacobian at quad pts
+ std::cout << "Jacobian det:" << std::endl;
+ _geometryFields->add("determinant(jacobian)", "determinant_jacobian");
+ topology::Field<mesh_type>& jacobianDetField =
+ _geometryFields->get("determinant(jacobian)");
+ fiberDim = _numQuadPts;
+ jacobianDetField.newSection(cStart, cEnd, fiberDim);
+ jacobianDetField.allocate();
+
+ // Allocate field for derivatives of basis functions at quad pts
+ std::cout << "Basis derivatives: num basis: " << _numBasis << std::endl;
+ _geometryFields->add("derivative basis functions",
+ "derivative_basis_functions");
+ topology::Field<mesh_type>& basisDerivField =
+ _geometryFields->get("jacobian");
+ fiberDim = _numQuadPts * _numBasis * _spaceDim;
+ basisDerivField.newSection(cStart, cEnd, fiberDim);
+ basisDerivField.allocate();
+
+ logger.stagePop();
+
+#if defined(ALE_MEM_LOGGING)
+ std::cout
+ << loggingStage << ": "
+ << logger.getNumAllocations(loggingStage)
+ << " allocations " << logger.getAllocationTotal(loggingStage)
+ << " bytes"
+ << std::endl
+ << loggingStage << ": "
+ << logger.getNumDeallocations(loggingStage)
+ << " deallocations " << logger.getDeallocationTotal(loggingStage)
+ << " bytes"
+ << std::endl;
+#endif
+
+ assert(0 != _geometry);
+ const int numBasis = _numBasis;
+ PetscErrorCode err;
+
+ DM dm = mesh.dmMesh();
+ scalar_array coordinatesCell(numBasis*_spaceDim);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMPlexGetCoordinateSection(dm, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dm, &coordVec);CHECK_PETSC_ERROR(err);
+
+ PetscSection quadPtsSection = quadPtsField->petscSection();
+ Vec quadPtsVec = quadPtsField->localVector();
+ PetscScalar *quadPtsArray;
+ PetscSection jacobianSection = jacobianField->petscSection();
+ Vec jacobianVec = jacobianField->localVector();
+ PetscScalar *jacobianArray;
+ PetscSection jacobianDetSection = jacobianDetField->petscSection();
+ Vec jacobianDetVec = jacobianDetField->localVector();
+ PetscScalar *jacobianDetArray;
+ PetscSection basisDerivSection = basisDerivField->petscSection();
+ Vec basisDerivVec = basisDerivField->localVector();
+ PetscScalar *basisDerivArray;
+
+ const scalar_array& quadPts = _engine->quadPts();
+ const scalar_array& jacobian = _engine->jacobian();
+ const scalar_array& jacobianDet = _engine->jacobianDet();
+ const scalar_array& basisDeriv = _engine->basisDeriv();
+
+ err = VecGetArray(quadPtsVec, &quadPtsArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(jacobianVec, &jacobianArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(jacobianDetVec, &jacobianDetArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(basisDerivVec, &basisDerivArray);CHECK_PETSC_ERROR(err);
+ for (PetscInt c = cStart; c < cEnd; ++c) {
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
+
+ err = DMPlexVecGetClosure(dm, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ err = DMPlexVecRestoreClosure(dm, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+
+ _engine->computeGeometry(coordinatesCell, c);
+ // Update fields with cell data
+ PetscInt dof, off;
+
+ err = PetscSectionGetDof(quadPtsSection, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(quadPtsSection, c, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ quadPtsArray[off+d] = quadPts[d];
+ }
+ err = PetscSectionGetDof(jacobianSection, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(jacobianSection, c, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ jacobianArray[off+d] = jacobian[d];
+ }
+ err = PetscSectionGetDof(jacobianDetSection, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(jacobianDetSection, c, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ jacobianDetArray[off+d] = jacobianDet[d];
+ }
+ err = PetscSectionGetDof(basisDerivSection, c, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(basisDerivSection, c, &off);CHECK_PETSC_ERROR(err);
+ for(PetscInt d = 0; d < dof; ++d) {
+ basisDerivArray[off+d] = basisDeriv[d];
+ }
+ } // for
+ err = VecRestoreArray(quadPtsVec, &quadPtsArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(jacobianVec, &jacobianArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(jacobianDetVec, &jacobianDetArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(basisDerivVec, &basisDerivArray);CHECK_PETSC_ERROR(err);
+} // computeGeometry
+
+// ----------------------------------------------------------------------
+template<typename mesh_type>
+void
pylith::feassemble::Quadrature<mesh_type>::retrieveGeometry(const typename mesh_type::SieveMesh::point_type& cell)
{ // retrieveGeometry
assert(0 != _geometryFields);
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/Quadrature.hh 2013-02-22 19:56:50 UTC (rev 21384)
@@ -125,6 +125,13 @@
void computeGeometry(const mesh_type& mesh,
const ALE::Obj<typename mesh_type::SieveMesh::label_sequence>& cells);
+ /** Compute geometric quantities for each cell.
+ *
+ * @param mesh Finite-element mesh
+ * @param cells Finite-element cells for geometry.
+ */
+ void computeGeometry(const mesh_type& mesh, PetscInt cStart, PetscInt cEnd);
+
/** Retrieve precomputed geometric quantities for a cell.
*
* @param mesh Finite-element mesh
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -38,12 +38,7 @@
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveSubMesh SieveSubMesh;
-typedef pylith::topology::Mesh::RealSection RealSection;
-typedef pylith::topology::SubMesh::RealUniformSection SubRealUniformSection;
-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) :
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -38,11 +38,7 @@
// ----------------------------------------------------------------------
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,
@@ -258,9 +254,20 @@
&_initialStressCell[iQuad*_tensorSize], _tensorSize,
&_initialStrainCell[iQuad*_tensorSize], _tensorSize);
- const ALE::Obj<RealSection>& stateVarsSection = _stateVars->section();
- assert(!stateVarsSection.isNull());
- stateVarsSection->updatePoint(cell, &_stateVarsCell[0]);
+ PetscSection stateVarsSection = _stateVars->petscSection();
+ Vec stateVarsVec = _stateVars->localVector();
+ PetscScalar *stateVarsArray;
+ PetscInt dof, off;
+ PetscErrorCode err;
+ assert(stateVarsSection);assert(stateVarsVec);
+
+ err = PetscSectionGetDof(stateVarsSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(stateVarsSection, cell, &off);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);
+ for (PetscInt d = 0; d < dof; ++d) {
+ stateVarsArray[off+d] = _stateVarsCell[d];
+ }
+ err = VecRestoreArray(stateVarsVec, &stateVarsArray);CHECK_PETSC_ERROR(err);
} // updateStateVars
// ----------------------------------------------------------------------
@@ -385,39 +392,45 @@
assert(_initialStrainCell.size() == numQuadPts*_tensorSize);
// Get cells associated with material
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->getLabelStratum("material-id", id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = mesh.dmMesh();
+ assert(dmMesh);
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+ err = DMPlexGetStratumIS(dmMesh, "material-id", id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetLocalSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+
// Setup field if necessary.
- ALE::Obj<RealSection> fieldSection;
+ PetscSection fieldSection;
+ Vec fieldVec;
+ PetscScalar *fieldArray;
if (field) {
const int fiberDim = 1*numQuadPts;
- fieldSection = field->section();
+ fieldSection = field->petscSection();
+ fieldVec = field->localVector();
bool useCurrentField = false;
- if (!fieldSection.isNull()) {
+ if (fieldSection) {
// check fiber dimension
- const int fiberDimCurrentLocal = (cells->size() > 0) ? fieldSection->getFiberDimension(*cells->begin()) : 0;
- int fiberDimCurrent = 0;
- MPI_Allreduce((void *) &fiberDimCurrentLocal,
- (void *) &fiberDimCurrent, 1,
- MPI_INT, MPI_MAX, field->mesh().comm());
+ PetscInt fiberDimCurrentLocal = 0;
+ if (numCells > 0) {err = PetscSectionGetDof(fieldSection, cells[0], &fiberDimCurrentLocal);CHECK_PETSC_ERROR(err);}
+ PetscInt fiberDimCurrent = 0;
+ MPI_Allreduce(&fiberDimCurrentLocal, &fiberDimCurrent, 1, MPIU_INT, MPI_MAX, field->mesh().comm());
assert(fiberDimCurrent > 0);
useCurrentField = fiberDim == fiberDimCurrent;
} // if
if (!useCurrentField) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("OutputFields");
- field->newSection(cells, fiberDim);
+ field->newSection(cells, numCells, fiberDim);
field->allocate();
- fieldSection = field->section();
+ fieldSection = field->petscSection();
+ fieldVec = field->localVector();
logger.stagePop();
} // if
- assert(!fieldSection.isNull());
+ assert(fieldSection);assert(fieldVec);
field->label("stable_dt_explicit");
assert(_normalizer);
field->scale(_normalizer->timeScale());
@@ -427,20 +440,23 @@
const int spaceDim = quadrature->spaceDim();
const int numBasis = quadrature->numBasis();
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
PylithScalar dtStable = pylith::PYLITH_MAXSCALAR;
scalar_array dtStableCell(numQuadPts);
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- retrievePropsAndVars(*c_iter);
+ if (field) {err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);}
+ for (PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
+ retrievePropsAndVars(cell);
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, cell, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
const double minCellWidth = quadrature->minCellWidth(coordinatesCell);
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
@@ -452,16 +468,25 @@
minCellWidth);
dtStableCell[iQuad] = dt;
if (dt < dtStable) {
- dtStable = dt;
+ dtStable = dt;
} // if
} // for
if (field) {
- assert(!fieldSection.isNull());
- assert(numQuadPts == fieldSection->getFiberDimension(*c_iter));
- fieldSection->updatePoint(*c_iter, &dtStableCell[0]);
+ PetscInt dof, off;
+
+ assert(fieldSection);
+ err = PetscSectionGetDof(fieldSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(numQuadPts == dof);
+ for (PetscInt d = 0; d < dof; ++d) {
+ fieldArray[off+d] = dtStableCell[d];
+ }
} // if
} // for
-
+ if (field) {err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);}
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
+
return dtStable;
} // stableTimeStepExplicit
@@ -476,35 +501,42 @@
if (field) {
const int numQuadPts = _numQuadPts;
- const ALE::Obj<RealSection>& fieldSection = field->section();
+ PetscSection fieldSection = field->petscSection();
+ Vec fieldVec = field->localVector();
+ PetscScalar *fieldArray;
// Get cells associated with material
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->getLabelStratum("material-id", id());
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
+ DM dmMesh = mesh.dmMesh();
+ assert(dmMesh);
+ IS cellIS;
+ const PetscInt *cells;
+ PetscInt numCells;
+ PetscErrorCode err;
+
+ err = DMPlexGetStratumIS(dmMesh, "material-id", id(), &cellIS);CHECK_PETSC_ERROR(err);
+ err = ISGetLocalSize(cellIS, &numCells);CHECK_PETSC_ERROR(err);
+ err = ISGetIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
const int fiberDim = 1*numQuadPts;
bool useCurrentField = false;
- if (!fieldSection.isNull()) {
+ if (fieldSection) {
// check fiber dimension
- const int fiberDimCurrentLocal = (cells->size() > 0) ? fieldSection->getFiberDimension(*cells->begin()) : 0;
- int fiberDimCurrent = 0;
- MPI_Allreduce((void *) &fiberDimCurrentLocal,
- (void *) &fiberDimCurrent, 1,
- MPI_INT, MPI_MAX, field->mesh().comm());
+ PetscInt fiberDimCurrentLocal = 0;
+ if (numCells > 0) {err = PetscSectionGetDof(fieldSection, cells[0], &fiberDimCurrentLocal);CHECK_PETSC_ERROR(err);}
+ PetscInt fiberDimCurrent = 0;
+ MPI_Allreduce(&fiberDimCurrentLocal, &fiberDimCurrent, 1, MPIU_INT, MPI_MAX, field->mesh().comm());
assert(fiberDimCurrent > 0);
useCurrentField = fiberDim == fiberDimCurrent;
} // if
if (!useCurrentField) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("OutputFields");
- field->newSection(cells, fiberDim);
+ field->newSection(cells, numCells, fiberDim);
field->allocate();
+ fieldSection = field->petscSection();
+ fieldVec = field->localVector();
logger.stagePop();
} // if
- assert(!fieldSection.isNull());
+ assert(fieldSection);
field->label("stable_dt_implicit");
assert(_normalizer);
field->scale(_normalizer->timeScale());
@@ -512,12 +544,22 @@
scalar_array dtStableCell(numQuadPts);
dtStableCell = PYLITH_MAXSCALAR;
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- assert(numQuadPts == fieldSection->getFiberDimension(*c_iter));
- fieldSection->updatePoint(*c_iter, &dtStableCell[0]);
+ err = VecGetArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ for (PetscInt c = 0; c < numCells; ++c) {
+ const PetscInt cell = cells[c];
+ PetscInt dof, off;
+
+ assert(fieldSection);
+ err = PetscSectionGetDof(fieldSection, cell, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(fieldSection, cell, &off);CHECK_PETSC_ERROR(err);
+ assert(numQuadPts == dof);
+ for (PetscInt d = 0; d < dof; ++d) {
+ fieldArray[off+d] = dtStableCell[d];
+ }
} // for
+ err = VecRestoreArray(fieldVec, &fieldArray);CHECK_PETSC_ERROR(err);
+ err = ISRestoreIndices(cellIS, &cells);CHECK_PETSC_ERROR(err);
+ err = ISDestroy(&cellIS);CHECK_PETSC_ERROR(err);
} // if
return dtStable;
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -38,11 +38,7 @@
// ----------------------------------------------------------------------
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::Material::Material(const int dimension,
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -25,6 +25,8 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
+#include <petscviewerhdf5.h>
+
// ----------------------------------------------------------------------
// Constructor
template<typename mesh_type, typename field_type>
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -219,60 +219,19 @@
typedef typename field_type::Mesh::RealSection RealSection;
try {
- DM complexMesh = mesh.dmMesh();
+ DM dmMesh = mesh.dmMesh();
+ assert(dmMesh);
- if (complexMesh) {
- /* DMPlex */
- Vec v = field.localVector(); /* Could check the field.petscSection() matches the default section from VecGetDM() */
- assert(v);
+ /* DMPlex */
+ Vec v = field.localVector(); /* Could check the field.petscSection() matches the default section from VecGetDM() */
+ assert(v);
- /* Will change to just VecView() once I setup the vectors correctly (use VecSetOperation() to change the view method) */
- PetscViewerVTKFieldType ft = field.vectorFieldType() != topology::FieldBase::VECTOR ? PETSC_VTK_POINT_FIELD : PETSC_VTK_POINT_VECTOR_FIELD;
- PetscErrorCode err = PetscViewerVTKAddField(_viewer, (PetscObject) complexMesh, DMPlexVTKWriteAll, ft, (PetscObject) v);CHECK_PETSC_ERROR(err);
- err = PetscObjectReference((PetscObject) v);CHECK_PETSC_ERROR(err); /* Needed because viewer destroys the Vec */
+ /* Will change to just VecView() once I setup the vectors correctly (use VecSetOperation() to change the view method) */
+ PetscViewerVTKFieldType ft = field.vectorFieldType() != topology::FieldBase::VECTOR ? PETSC_VTK_POINT_FIELD : PETSC_VTK_POINT_VECTOR_FIELD;
+ PetscErrorCode err = PetscViewerVTKAddField(_viewer, (PetscObject) dmMesh, DMPlexVTKWriteAll, ft, (PetscObject) v);CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) v);CHECK_PETSC_ERROR(err); /* Needed because viewer destroys the Vec */
- _wroteVertexHeader = true;
- } else {
- int rank = 0;
- MPI_Comm_rank(field.mesh().comm(), &rank);
-
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
-
- const std::string labelName =
- (sieveMesh->hasLabel("censored depth")) ? "censored depth" : "depth";
- const ALE::Obj<typename SieveMesh::numbering_type>& numbering =
- sieveMesh->getFactory()->getNumbering(sieveMesh, labelName, 0);
- assert(!numbering.isNull());
-
- PetscSection sectionP = field.petscSection();
- PetscInt dof = 0;
-
- assert(sectionP);
- assert(!sieveMesh->getLabelStratum(labelName, 0).isNull());
- if (sieveMesh->getLabelStratum(labelName, 0)->size() > 0) {
- PetscErrorCode err = PetscSectionGetDof(sectionP, *sieveMesh->getLabelStratum(labelName, 0)->begin(), &dof);CHECK_PETSC_ERROR(err);
- }
- int fiberDimLocal = dof;
- int fiberDim = 0;
- MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
- field.mesh().comm());
- assert(fiberDim > 0);
- const int enforceDim =
- (field.vectorFieldType() != topology::FieldBase::VECTOR) ? fiberDim : 3;
-
- PetscErrorCode err = 0;
- if (!_wroteVertexHeader) {
- err = PetscViewerASCIIPrintf(_viewer, "POINT_DATA %d\n",
- numbering->getGlobalSize());
- CHECK_PETSC_ERROR(err);
- _wroteVertexHeader = true;
- } // if
-
- err = VTKViewer::writeField(field.section(), field.label(), fiberDim, numbering,
- _viewer, enforceDim, _precision);
- CHECK_PETSC_ERROR(err);
- }
+ _wroteVertexHeader = true;
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
@@ -300,69 +259,20 @@
typedef typename field_type::Mesh::RealSection RealSection;
try {
- DM complexMesh = field.mesh().dmMesh();
+ DM dmMesh = field.mesh().dmMesh();
+ assert(dmMesh);
- if (complexMesh) {
- /* DMPlex */
- PetscContainer c;
- PetscSection s = field.petscSection();
- Vec v = field.localVector();
- assert(s);assert(v);
+ /* DMPlex */
+ PetscSection s = field.petscSection();
+ Vec v = field.localVector();
+ assert(s);assert(v);
- /* Will change to just VecView() once I setup the vectors correctly (use VecSetOperation() to change the view) */
- PetscViewerVTKFieldType ft = field.vectorFieldType() != topology::FieldBase::VECTOR ? PETSC_VTK_CELL_FIELD : PETSC_VTK_CELL_VECTOR_FIELD;
- PetscErrorCode err = PetscViewerVTKAddField(_viewer, (PetscObject) complexMesh, DMPlexVTKWriteAll, ft, (PetscObject) v); CHECK_PETSC_ERROR(err);
- err = PetscObjectReference((PetscObject) v);CHECK_PETSC_ERROR(err); /* Needed because viewer destroys the Vec */
- err = PetscContainerCreate(((PetscObject) v)->comm, &c);CHECK_PETSC_ERROR(err);
- err = PetscContainerSetPointer(c, s);CHECK_PETSC_ERROR(err);
- err = PetscObjectCompose((PetscObject) v, "section", (PetscObject) c);CHECK_PETSC_ERROR(err);
- err = PetscContainerDestroy(&c);CHECK_PETSC_ERROR(err);
+ /* Will change to just VecView() once I setup the vectors correctly (use VecSetOperation() to change the view) */
+ PetscViewerVTKFieldType ft = field.vectorFieldType() != topology::FieldBase::VECTOR ? PETSC_VTK_CELL_FIELD : PETSC_VTK_CELL_VECTOR_FIELD;
+ PetscErrorCode err = PetscViewerVTKAddField(_viewer, (PetscObject) dmMesh, DMPlexVTKWriteAll, ft, (PetscObject) v); CHECK_PETSC_ERROR(err);
+ err = PetscObjectReference((PetscObject) v);CHECK_PETSC_ERROR(err); /* Needed because viewer destroys the Vec */
- _wroteCellHeader = true;
- } else {
- PetscErrorCode err = 0;
-
- // Correctly handle boundary and fault meshes
- // Cannot just use mesh->depth() because boundaries report the wrong thing
- const ALE::Obj<SieveMesh>& sieveMesh = field.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
- int cellDepth = 0;
- err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
- sieveMesh->comm());CHECK_PETSC_ERROR(err);
- const int depth = (!label) ? cellDepth : labelId;
- const std::string labelName = (!label) ?
- ((sieveMesh->hasLabel("censored depth")) ?
- "censored depth" : "depth") : label;
- assert(!sieveMesh->getFactory().isNull());
- const ALE::Obj<typename SieveMesh::numbering_type>& numbering =
- sieveMesh->getFactory()->getNumbering(sieveMesh, labelName, depth);
- assert(!numbering.isNull());
- const ALE::Obj<RealSection>& section = field.section();
- assert(!section.isNull());
-
- assert(!sieveMesh->getLabelStratum(labelName, depth).isNull());
- int fiberDimLocal =
- (sieveMesh->getLabelStratum(labelName, depth)->size() > 0) ?
- section->getFiberDimension(*sieveMesh->getLabelStratum(labelName, depth)->begin()) : 0;
- int fiberDim = 0;
- MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
- field.mesh().comm());
- assert(fiberDim > 0);
-
- const int enforceDim =
- (field.vectorFieldType() != topology::FieldBase::VECTOR) ? fiberDim : 3;
-
- if (!_wroteCellHeader) {
- err = PetscViewerASCIIPrintf(_viewer, "CELL_DATA %d\n",
- numbering->getGlobalSize());
- CHECK_PETSC_ERROR(err);
- _wroteCellHeader = true;
- } // if
-
- VTKViewer::writeField(section, field.label(), fiberDim, numbering,
- _viewer, enforceDim, _precision);
- }
+ _wroteCellHeader = true;
} catch (const std::exception& err) {
std::ostringstream msg;
msg << "Error while writing field '" << field.label() << "' at time "
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Explicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Explicit.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Explicit.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -61,62 +61,85 @@
const int spaceDim = cs->spaceDim();
// Get sections.
- const ALE::Obj<RealSection>& dispIncrSection = dispIncr.section();
- assert(!dispIncrSection.isNull());
+ PetscSection dispIncrSection = dispIncr.petscSection();
+ Vec dispIncrVec = dispIncr.localVector();
+ PetscScalar *dispIncrArray;
+ assert(dispIncrSection);assert(dispIncrVec);
- const ALE::Obj<RealSection>& dispTSection = _fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
+ PetscSection dispTSection = _fields->get("disp(t)").petscSection();
+ Vec dispTVec = _fields->get("disp(t)").localVector();
+ PetscScalar *dispTArray;
+ assert(dispTSection);assert(dispTVec);
- const ALE::Obj<RealSection>& dispTmdtSection =
- _fields->get("disp(t-dt)").section();
- assert(!dispTmdtSection.isNull());
+ PetscSection dispTmdtSection = _fields->get("disp(t-dt)").petscSection();
+ Vec dispTmdtVec = _fields->get("disp(t-dt)").localVector();
+ PetscScalar *dispTmdtArray;
+ assert(dispTmdtSection);assert(dispTmdtVec);
- scalar_array velVertex(spaceDim);
- const ALE::Obj<RealSection>& velSection =
- _fields->get("velocity(t)").section();
- assert(!velSection.isNull());
+ PetscSection velSection = _fields->get("velocity(t)").petscSection();
+ Vec velVec = _fields->get("velocity(t)").localVector();
+ PetscScalar *velArray;
+ assert(velSection);assert(velVec);
- scalar_array accVertex(spaceDim);
- const ALE::Obj<RealSection>& accSection =
- _fields->get("acceleration(t)").section();
- assert(!accSection.isNull());
+ PetscSection accSection = _fields->get("acceleration(t)").petscSection();
+ Vec accVec = _fields->get("acceleration(t)").localVector();
+ PetscScalar *accArray;
+ assert(accSection);assert(accVec);
// Get mesh vertices.
- const ALE::Obj<SieveMesh>& sieveMesh = dispIncr.mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- assert(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
-
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- assert(spaceDim == dispIncrSection->getFiberDimension(*v_iter));
- const PylithScalar* dispIncrVertex = dispIncrSection->restrictPoint(*v_iter);
+ DM dmMesh = dispIncr.mesh().dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
- assert(spaceDim == dispTSection->getFiberDimension(*v_iter));
- const PylithScalar* dispTVertex = dispTSection->restrictPoint(*v_iter);
+ assert(dmMesh);
+ err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
- assert(spaceDim == dispTmdtSection->getFiberDimension(*v_iter));
- const PylithScalar* dispTmdtVertex = dispTmdtSection->restrictPoint(*v_iter);
+ err = VecGetArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(dispTmdtVec, &dispTmdtArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(accVec, &accArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt didof, dioff;
- for (int i=0; i < spaceDim; ++i) {
- velVertex[i] =
- (dispIncrVertex[i] + dispTVertex[i] - dispTmdtVertex[i]) / twodt;
- accVertex[i] =
- (dispIncrVertex[i] - dispTVertex[i] + dispTmdtVertex[i]) / dt2;
+ err = PetscSectionGetDof(dispIncrSection, v, &didof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispIncrSection, v, &dioff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == didof);
+ PetscInt dtdof, dtoff;
+
+ err = PetscSectionGetDof(dispTSection, v, &dtdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTSection, v, &dtoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dtdof);
+ PetscInt dmdof, dmoff;
+
+ err = PetscSectionGetDof(dispTmdtSection, v, &dmdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(dispTmdtSection, v, &dmoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == dmdof);
+ PetscInt vdof, voff;
+
+ err = PetscSectionGetDof(velSection, v, &vdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(velSection, v, &voff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == vdof);
+ PetscInt adof, aoff;
+
+ err = PetscSectionGetDof(accSection, v, &adof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(accSection, v, &aoff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == adof);
+
+ // TODO: I am not sure why these were updateAll() before, but if BCs need to be changed, then
+ // the global update will probably need to be modified
+ for (PetscInt i = 0; i < spaceDim; ++i) {
+ velArray[voff+i] = (dispIncrArray[dioff+i] + dispTArray[dtoff+i] - dispTmdtArray[dmoff+i]) / twodt;
+ accArray[aoff+i] = (dispIncrArray[dioff+i] - dispTArray[dtoff+i] + dispTmdtArray[dmoff+i]) / dt2;
} // for
-
- assert(velSection->getFiberDimension(*v_iter) == spaceDim);
- velSection->updatePointAll(*v_iter, &velVertex[0]);
-
- assert(accSection->getFiberDimension(*v_iter) == spaceDim);
- accSection->updatePointAll(*v_iter, &accVertex[0]);
} // for
+ err = VecRestoreArray(dispIncrVec, &dispIncrArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTVec, &dispTArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(dispTmdtVec, &dispTmdtArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(velVec, &velArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(accVec, &accArray);CHECK_PETSC_ERROR(err);
- PetscLogFlops(vertices->size() * 6*spaceDim);
+ PetscLogFlops((vEnd - vStart) * 6*spaceDim);
} // calcRateFields
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLumped.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLumped.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/SolverLumped.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -89,47 +89,61 @@
const int spaceDim = cs->spaceDim();
// Get mesh vertices.
- const ALE::Obj<SieveMesh>& sieveMesh = solution->mesh().sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<SieveMesh::label_sequence>& vertices =
- sieveMesh->depthStratum(0);
- assert(!vertices.isNull());
- const SieveMesh::label_sequence::iterator verticesBegin =
- vertices->begin();
- const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ DM dmMesh = solution->mesh().dmMesh();
+ PetscInt vStart, vEnd;
+ PetscErrorCode err;
+
+ assert(dmMesh);
+ err = DMPlexGetDepthStratum(dmMesh, 0, &vStart, &vEnd);CHECK_PETSC_ERROR(err);
// Get sections.
- scalar_array solutionVertex(spaceDim);
- const ALE::Obj<RealSection>& solutionSection = solution->section();
- assert(!solutionSection.isNull());
+ PetscSection solutionSection = solution->petscSection();
+ Vec solutionVec = solution->localVector();
+ PetscScalar *solutionArray;
+ assert(solutionSection);assert(solutionVec);
- const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
- assert(!jacobianSection.isNull());
+ PetscSection jacobianSection = jacobian.petscSection();
+ Vec jacobianVec = jacobian.localVector();
+ PetscScalar *jacobianArray;
+ assert(jacobianSection);assert(jacobianVec);
- const ALE::Obj<RealSection>& residualSection = residual.section();
- assert(!residualSection.isNull());
+ PetscSection residualSection = residual.petscSection();
+ Vec residualVec = residual.localVector();
+ PetscScalar *residualArray;
+ assert(residualSection);assert(residualVec);
_logger->eventEnd(setupEvent);
_logger->eventBegin(solveEvent);
- for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- assert(spaceDim == jacobianSection->getFiberDimension(*v_iter));
- const PylithScalar* jacobianVertex = jacobianSection->restrictPoint(*v_iter);
+ err = VecGetArray(jacobianVec, &jacobianArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(solutionVec, &solutionArray);CHECK_PETSC_ERROR(err);
+ for(PetscInt v = vStart; v < vEnd; ++v) {
+ PetscInt jdof, joff;
- assert(spaceDim == residualSection->getFiberDimension(*v_iter));
- const PylithScalar* residualVertex = residualSection->restrictPoint(*v_iter);
+ err = PetscSectionGetDof(jacobianSection, v, &jdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(jacobianSection, v, &joff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == jdof);
+ PetscInt rdof, roff;
+ err = PetscSectionGetDof(residualSection, v, &rdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(residualSection, v, &roff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == rdof);
+ PetscInt sdof, soff;
+
+ err = PetscSectionGetDof(solutionSection, v, &sdof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(solutionSection, v, &soff);CHECK_PETSC_ERROR(err);
+ assert(spaceDim == sdof);
+
for (int i=0; i < spaceDim; ++i) {
- assert(jacobianVertex[i] != 0.0);
- solutionVertex[i] = residualVertex[i] / jacobianVertex[i];
+ assert(jacobianArray[joff+i] != 0.0);
+ solutionArray[soff+i] = residualArray[roff+i] / jacobianArray[joff+i];
} // for
-
- assert(solutionSection->getFiberDimension(*v_iter) == spaceDim);
- solutionSection->updatePoint(*v_iter, &solutionVertex[0]);
} // for
- PetscLogFlops(vertices->size() * spaceDim);
+ PetscLogFlops((vEnd - vStart) * spaceDim);
+ err = VecRestoreArray(jacobianVec, &jacobianArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(residualVec, &residualArray);CHECK_PETSC_ERROR(err);
+ err = VecRestoreArray(solutionVec, &solutionArray);CHECK_PETSC_ERROR(err);
_logger->eventEnd(solveEvent);
_logger->eventBegin(adjustEvent);
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -123,23 +123,25 @@
partition.scale(1.0);
partition.label("partition");
partition.vectorFieldType(topology::FieldBase::SCALAR);
- const ALE::Obj<RealSection>& partitionSection = partition.section();
- assert(!partitionSection.isNull());
+ PetscSection partitionSection = partition.petscSection();
+ Vec partitionVec = partition.localVector();
+ PetscScalar *partitionArray;
+ assert(partitionSection);assert(partitionVec);
- const ALE::Obj<SieveMesh> sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- PylithScalar rankReal = PylithScalar(commRank);
- assert(sieveMesh->height() > 0);
- const ALE::Obj<SieveMesh::label_sequence>& cells =
- sieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- partitionSection->updatePoint(*c_iter, &rankReal);
+ PylithScalar rankReal = PylithScalar(commRank);
+ DM dmMesh = mesh.dmMesh();
+ PetscInt cStart, cEnd;
+ PetscErrorCode err;
+
+ assert(dmMesh);
+ err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(partitionVec, &partitionArray);CHECK_PETSC_ERROR(err);
+ for (PetscInt c = cStart; c < cEnd; ++c) {
+ PetscInt off;
+ err = PetscSectionGetOffset(partitionSection, c, &off);CHECK_PETSC_ERROR(err);
+ partitionArray[off] = rankReal;
} // for
+ err = VecRestoreArray(partitionVec, &partitionArray);CHECK_PETSC_ERROR(err);
//partition->view("PARTITION");
const PylithScalar t = 0.0;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -32,8 +32,8 @@
// ----------------------------------------------------------------------
// Default constructor.
-template<typename mesh_type, typename section_type>
-pylith::topology::Field<mesh_type, section_type>::Field(const mesh_type& mesh) :
+template<typename mesh_type>
+pylith::topology::Field<mesh_type>::Field(const mesh_type& mesh) :
_mesh(mesh)
{ // constructor
_metadata["default"].label = "unknown";
@@ -69,49 +69,9 @@
} // constructor
// ----------------------------------------------------------------------
-// Constructor with mesh, section, and metadata.
-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) :
- _mesh(mesh),
- _section(section)
-{ // constructor
- assert(!section.isNull());
- _metadata["default"] = metadata;
- if (mesh.dmMesh()) {
- DM dm = mesh.dmMesh(), coordDM, newCoordDM;
- PetscSection coordSection, newCoordSection;
- Vec coordVec;
- PetscSection s;
- PetscErrorCode err;
-
- err = DMPlexClone(dm, &_dm);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinatesLocal(dm, &coordVec);CHECK_PETSC_ERROR(err);
- if (coordVec) {
- err = DMGetCoordinateDM(dm, &coordDM);CHECK_PETSC_ERROR(err);
- err = DMGetCoordinateDM(_dm, &newCoordDM);CHECK_PETSC_ERROR(err);
- err = DMGetDefaultSection(coordDM, &coordSection);CHECK_PETSC_ERROR(err);
- err = PetscSectionClone(coordSection, &newCoordSection);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultSection(newCoordDM, newCoordSection);CHECK_PETSC_ERROR(err);
- err = DMSetCoordinatesLocal(_dm, coordVec);CHECK_PETSC_ERROR(err);
- }
- this->dmSection(&s, &_localVec);
- err = DMSetDefaultSection(_dm, s);CHECK_PETSC_ERROR(err);
- err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) _localVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- } else {
- _dm = PETSC_NULL;
- _globalVec = PETSC_NULL;
- _localVec = PETSC_NULL;
- }
-} // constructor
-
-// ----------------------------------------------------------------------
// Constructor with mesh, DM, and metadata
-template<typename mesh_type, typename section_type>
-pylith::topology::Field<mesh_type, section_type>::Field(const mesh_type& mesh, DM dm, const Metadata& metadata) :
+template<typename mesh_type>
+pylith::topology::Field<mesh_type>::Field(const mesh_type& mesh, DM dm, const Metadata& metadata) :
_mesh(mesh),
_dm(dm)
{ // constructor
@@ -126,11 +86,9 @@
// ----------------------------------------------------------------------
// Constructor with field and subfields
-template<typename mesh_type, typename section_type>
-pylith::topology::Field<mesh_type, section_type>::Field(const Field& src,
- const int fields[], int numFields) :
- _mesh(src._mesh),
- _section(PETSC_NULL)
+template<typename mesh_type>
+pylith::topology::Field<mesh_type>::Field(const Field& src, const int fields[], int numFields) :
+ _mesh(src._mesh)
{ // constructor
DM dm = mesh.dmMesh(), coordDM, newCoordDM;
PetscSection coordSection, newCoordSection;
@@ -162,8 +120,8 @@
// ----------------------------------------------------------------------
// Destructor.
-template<typename mesh_type, typename section_type>
-pylith::topology::Field<mesh_type, section_type>::~Field(void)
+template<typename mesh_type>
+pylith::topology::Field<mesh_type>::~Field(void)
{ // destructor
deallocate();
PetscErrorCode err = DMDestroy(&_dm);CHECK_PETSC_ERROR(err);
@@ -171,9 +129,9 @@
// ----------------------------------------------------------------------
// Deallocate PETSc and local data structures.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::deallocate(void)
+pylith::topology::Field<mesh_type>::deallocate(void)
{ // deallocate
PetscErrorCode err = 0;
@@ -195,38 +153,31 @@
// ----------------------------------------------------------------------
// Set label for field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::label(const char* value)
+pylith::topology::Field<mesh_type>::label(const char* value)
{ // label
+ PetscErrorCode err;
+
_metadata["default"].label = value;
- if (!_section.isNull()) {
- _section->setName(value);
- } // if
- if (_localVec) {
- PetscErrorCode err = PetscObjectSetName((PetscObject) _localVec, value);CHECK_PETSC_ERROR(err);
- }
- if (_globalVec) {
- PetscErrorCode err = PetscObjectSetName((PetscObject) _globalVec, value);CHECK_PETSC_ERROR(err);
- }
+ if (_localVec) {err = PetscObjectSetName((PetscObject) _localVec, value);CHECK_PETSC_ERROR(err);}
+ if (_globalVec) {err = PetscObjectSetName((PetscObject) _globalVec, value);CHECK_PETSC_ERROR(err);}
const typename scatter_map_type::const_iterator scattersEnd = _scatters.end();
for (typename scatter_map_type::const_iterator s_iter=_scatters.begin();
s_iter != scattersEnd;
++s_iter) {
if (s_iter->second.vector) {
- PetscErrorCode err =
- PetscObjectSetName((PetscObject)s_iter->second.vector, value);
- CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject)s_iter->second.vector, value);CHECK_PETSC_ERROR(err);
} // if
} // for
} // label
// ----------------------------------------------------------------------
// Get spatial dimension of domain.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
int
-pylith::topology::Field<mesh_type, section_type>::spaceDim(void) const
+pylith::topology::Field<mesh_type>::spaceDim(void) const
{ // spaceDim
const spatialdata::geocoords::CoordSys* cs = _mesh.coordsys();
return (cs) ? cs->spaceDim() : 0;
@@ -234,33 +185,31 @@
// ----------------------------------------------------------------------
// Get the chart size.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
int
-pylith::topology::Field<mesh_type, section_type>::chartSize(void) const
+pylith::topology::Field<mesh_type>::chartSize(void) const
{ // chartSize
- if (_dm) {
- PetscSection s;
- PetscInt pStart, pEnd;
- PetscErrorCode err;
+ assert(_dm);
+ PetscSection s;
+ PetscInt pStart, pEnd;
+ PetscErrorCode err;
- err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetChart(s, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- return pEnd-pStart;
- }
- return _section.isNull() ? 0 : _section->getChart().size();
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetChart(s, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ return pEnd-pStart;
} // chartSize
// ----------------------------------------------------------------------
// Get the number of degrees of freedom.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
int
-pylith::topology::Field<mesh_type, section_type>::sectionSize(void) const
+pylith::topology::Field<mesh_type>::sectionSize(void) const
{ // sectionSize
- PetscInt size = 0;
- PetscErrorCode err;
+ PetscInt size = 0;
if (_dm) {
- PetscSection s;
+ PetscSection s;
+ PetscErrorCode err;
err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
err = PetscSectionGetStorageSize(s, &size);CHECK_PETSC_ERROR(err);
@@ -270,78 +219,58 @@
// ----------------------------------------------------------------------
// Create seive section.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::newSection(void)
+pylith::topology::Field<mesh_type>::newSection(void)
{ // newSection
// Clear memory
clear();
-
- ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Field");
-
- _section = new section_type(_mesh.comm(), _mesh.debug());
- assert(!_section.isNull());
- _section->setName(_metadata["default"].label);
-
- logger.stagePop();
} // newSection
// ----------------------------------------------------------------------
-// Create sieve section and set chart and fiber dimesion for a
-// sequence of points.
-template<typename mesh_type, typename section_type>
+// Create sieve section and set chart and fiber dimesion for a list of
+// points.
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::newSection(
- const ALE::Obj<label_sequence>& points,
- const int fiberDim)
+pylith::topology::Field<mesh_type>::newSection(const int_array& points,
+ const int fiberDim)
{ // newSection
typedef typename mesh_type::SieveMesh::point_type point_type;
PetscErrorCode err;
// Clear memory
clear();
-
- ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Field");
-
+ assert(_dm);
if (fiberDim < 0) {
std::ostringstream msg;
msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata["default"].label
<< "' must be nonnegative.";
throw std::runtime_error(msg.str());
} // if
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+ const PetscInt npts = points.size();
+ if (npts > 0) {
+ PetscSection s;
+ PetscInt pointMin = 0, pointMax = 0;
- _section = new section_type(_mesh.comm(), _mesh.debug());
- assert(!_section.isNull());
- _section->setName(_metadata["default"].label);
-
- if (points->size() > 0) {
- const point_type pointMin =
- *std::min_element(points->begin(), points->end());
- const point_type pointMax =
- *std::max_element(points->begin(), points->end());
- _section->setChart(chart_type(pointMin, pointMax+1));
- _section->setFiberDimension(points, fiberDim);
-
- if (_dm) {
- PetscSection s;
- err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetChart(s, pointMin, pointMax+1);CHECK_PETSC_ERROR(err);
-
- for(typename label_sequence::const_iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
- err = PetscSectionSetDof(s, *p_iter, fiberDim);CHECK_PETSC_ERROR(err);
- }
+ for (PetscInt i = 0; i < npts; ++i) {
+ pointMin = std::min(pointMin, points[i]);
+ pointMax = std::max(pointMax, points[i]);
}
- } else {// Create empty chart
- _section->setChart(chart_type(0, 0));
- if (_dm) {
- PetscSection s;
- err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetChart(s, 0, 0);CHECK_PETSC_ERROR(err);
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, pointMin, pointMax+1);CHECK_PETSC_ERROR(err);
+ for (PetscInt i = 0; i < npts; ++i) {
+ err = PetscSectionSetDof(s, points[i], fiberDim);CHECK_PETSC_ERROR(err);
}
+ } else { // create empty chart
+ PetscSection s;
+
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, 0, 0);CHECK_PETSC_ERROR(err);
}
logger.stagePop();
@@ -350,19 +279,16 @@
// ----------------------------------------------------------------------
// Create sieve section and set chart and fiber dimesion for a list of
// points.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::newSection(const int_array& points,
- const int fiberDim)
+pylith::topology::Field<mesh_type>::newSection(const PetscInt *points, const PetscInt num,
+ const int fiberDim)
{ // newSection
- typedef typename mesh_type::SieveMesh::point_type point_type;
PetscErrorCode err;
// Clear memory
clear();
-
- ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Field");
+ assert(_dm);
if (fiberDim < 0) {
std::ostringstream msg;
msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata["default"].label
@@ -370,36 +296,28 @@
throw std::runtime_error(msg.str());
} // if
- _section = new section_type(_mesh.comm(), _mesh.debug());
- assert(!_section.isNull());
- _section->setName(_metadata["default"].label);
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("Field");
+ if (num > 0) {
+ PetscSection s;
+ PetscInt pointMin = 0, pointMax = 0;
- const int npts = points.size();
- if (npts > 0) {
- const point_type pointMin = points.min();
- const point_type pointMax = points.max();
- _section->setChart(chart_type(pointMin, pointMax+1));
- for (int i=0; i < npts; ++i)
- _section->setFiberDimension(points[i], fiberDim);
-
- if (_dm) {
- PetscSection s;
- err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetChart(s, pointMin, pointMax+1);CHECK_PETSC_ERROR(err);
-
- for (int i=0; i < npts; ++i) {
- err = PetscSectionSetDof(s, points[i], fiberDim);CHECK_PETSC_ERROR(err);
- }
+ for (PetscInt i = 0; i < num; ++i) {
+ pointMin = std::min(pointMin, points[i]);
+ pointMax = std::max(pointMax, points[i]);
}
- } else { // create empty chart
- _section->setChart(chart_type(0, 0));
- if (_dm) {
- PetscSection s;
- err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetChart(s, 0, 0);CHECK_PETSC_ERROR(err);
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, pointMin, pointMax+1);CHECK_PETSC_ERROR(err);
+ for (PetscInt i = 0; i < num; ++i) {
+ err = PetscSectionSetDof(s, points[i], fiberDim);CHECK_PETSC_ERROR(err);
}
+ } else { // create empty chart
+ PetscSection s;
+
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, 0, 0);CHECK_PETSC_ERROR(err);
}
logger.stagePop();
@@ -407,81 +325,73 @@
// ----------------------------------------------------------------------
// Create sieve section and set chart and fiber dimesion.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::newSection(const DomainEnum domain,
- const int fiberDim,
- const int stratum)
+pylith::topology::Field<mesh_type>::newSection(const DomainEnum domain,
+ const int fiberDim,
+ const int stratum)
{ // newSection
// Changing this because cells/vertices are numbered differently in the new scheme
- if (_dm) {
- PetscSection s;
- PetscInt pStart, pEnd;
- PetscErrorCode err;
+ assert(_dm);
+ PetscSection s;
+ PetscInt pStart, pEnd;
+ PetscErrorCode err;
- switch(domain) {
- case VERTICES_FIELD:
- err = DMPlexGetDepthStratum(_dm, stratum, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- break;
- case CELLS_FIELD:
- err = DMPlexGetHeightStratum(_dm, stratum, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- break;
- case FACES_FIELD:
- err = DMPlexGetHeightStratum(_dm, stratum+1, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- break;
- case POINTS_FIELD:
- err = DMPlexGetChart(_dm, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- break;
- default:
- std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
- throw std::logic_error("Bad domain enum in Field.");
- }
- err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetChart(s, pStart, pEnd);CHECK_PETSC_ERROR(err);
-
- for(PetscInt p = pStart; p < pEnd; ++p) {
- err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
- }
- } else {
- const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
- assert(!sieveMesh.isNull());
-
- ALE::Obj<label_sequence> points;
- if (VERTICES_FIELD == domain)
- points = sieveMesh->depthStratum(stratum);
- else if (CELLS_FIELD == domain)
- points = sieveMesh->heightStratum(stratum);
- else if (FACES_FIELD == domain)
- points = sieveMesh->heightStratum(stratum+1);
- else {
+ switch(domain) {
+ case VERTICES_FIELD:
+ err = DMPlexGetDepthStratum(_dm, stratum, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
+ case CELLS_FIELD:
+ err = DMPlexGetHeightStratum(_dm, stratum, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
+ case FACES_FIELD:
+ err = DMPlexGetHeightStratum(_dm, stratum+1, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
+ case POINTS_FIELD:
+ err = DMPlexGetChart(_dm, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ break;
+ default:
std::cerr << "Unknown value for DomainEnum: " << domain << std::endl;
- assert(0);
throw std::logic_error("Bad domain enum in Field.");
- } // else
+ }
+ newSection(pStart, pEnd, fiberDim);
+} // newSection
- newSection(points, fiberDim);
+// ----------------------------------------------------------------------
+// Create sieve section and set chart and fiber dimesion.
+template<typename mesh_type>
+void
+pylith::topology::Field<mesh_type>::newSection(const PetscInt pStart, const PetscInt pEnd,
+ const int fiberDim)
+{ // newSection
+ // Changing this because cells/vertices are numbered differently in the new scheme
+ assert(_dm);
+ PetscSection s;
+ PetscErrorCode err;
+
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, pStart, pEnd);CHECK_PETSC_ERROR(err);
+
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
}
} // newSection
// ----------------------------------------------------------------------
// Create section given chart.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::newSection(const Field& src,
- const int fiberDim)
+pylith::topology::Field<mesh_type>::newSection(const Field& src,
+ const int fiberDim)
{ // newSection
// Clear memory
clear();
+ assert(_dm);assert(src._dm);
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
- if (_section.isNull()) {
- logger.stagePop();
- newSection();
- logger.stagePush("Field");
- } // if
if (fiberDim < 0) {
std::ostringstream msg;
msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata["default"].label
@@ -489,322 +399,115 @@
throw std::runtime_error(msg.str());
} // if
- const ALE::Obj<section_type>& srcSection = src._section;
- if (!srcSection.isNull()) {
- _section->setChart(srcSection->getChart());
- const chart_type& chart = _section->getChart();
- const typename chart_type::const_iterator chartBegin = chart.begin();
- const typename chart_type::const_iterator chartEnd = chart.end();
+ PetscSection srcs, s;
+ PetscInt pStart, pEnd;
+ PetscErrorCode err;
- for (typename chart_type::const_iterator c_iter = chartBegin;
- c_iter != chartEnd;
- ++c_iter)
- if (srcSection->getFiberDimension(*c_iter) > 0)
- _section->setFiberDimension(*c_iter, fiberDim);
+ err = DMGetDefaultSection(src._dm, &srcs);CHECK_PETSC_ERROR(err);
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetChart(srcs, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ err = PetscSectionSetChart(s, pStart, pEnd);CHECK_PETSC_ERROR(err);
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
+ }
- if (_dm) {
- PetscSection s;
- PetscInt pStart = srcSection->getChart().min(), pEnd = srcSection->getChart().max();
- PetscErrorCode err;
-
- err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetChart(s, pStart, pEnd);CHECK_PETSC_ERROR(err);
- for(PetscInt p = pStart; p < pEnd; ++p) {
- err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
- }
- }
- } else if (src._dm) {
- PetscSection srcs, s;
- PetscInt pStart, pEnd;
- PetscErrorCode err;
-
- err = DMGetDefaultSection(src._dm, &srcs);CHECK_PETSC_ERROR(err);
- err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultGlobalSection(_dm, PETSC_NULL);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetChart(srcs, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetChart(s, pStart, pEnd);CHECK_PETSC_ERROR(err);
- for(PetscInt p = pStart; p < pEnd; ++p) {
- err = PetscSectionSetDof(s, p, fiberDim);CHECK_PETSC_ERROR(err);
- }
- } // if
-
logger.stagePop();
} // newSection
// ----------------------------------------------------------------------
// Create section with same layout as another section.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::cloneSection(const Field& src)
+pylith::topology::Field<mesh_type>::cloneSection(const Field& src)
{ // cloneSection
std::string origLabel = _metadata["default"].label;
// Clear memory
clear();
- const ALE::Obj<section_type>& srcSection = src._section;
- if (!srcSection.isNull() && _section.isNull()) {
- newSection();
- } // if
-
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
_metadata["default"] = const_cast<Field&>(src)._metadata["default"];
label(origLabel.c_str());
- if (!_section.isNull()) {
- if (!srcSection->sharedStorage()) {
- _section->setAtlas(srcSection->getAtlas());
- _section->allocateStorage();
- _section->setBC(srcSection->getBC());
- _section->copySpaces(srcSection);
- } else {
- _section->setChart(srcSection->getChart());
- const chart_type& chart = _section->getChart();
- const typename chart_type::const_iterator chartBegin = chart.begin();
- const typename chart_type::const_iterator chartEnd = chart.end();
- for (typename chart_type::const_iterator c_iter = chartBegin;
- c_iter != chartEnd;
- ++c_iter) {
- const int fiberDim = srcSection->getFiberDimension(*c_iter);
- if (fiberDim > 0)
- _section->setFiberDimension(*c_iter, fiberDim);
- } // for
- const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- sieveMesh->allocate(_section);
- _section->setBC(srcSection->getBC());
- _section->copySpaces(srcSection); // :BUG: NEED TO REBUILD SPACES
- } // if/else
- } // if
PetscSection section = src.petscSection();
PetscSection newSection;
PetscErrorCode err;
- if (_dm) {
- err = PetscSectionClone(section, &newSection);CHECK_PETSC_ERROR(err);
- err = DMSetDefaultSection(_dm, newSection);CHECK_PETSC_ERROR(err);
- err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
- err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) _localVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- }
+ assert(_dm);
+ err = PetscSectionClone(section, &newSection);CHECK_PETSC_ERROR(err);
+ err = DMSetDefaultSection(_dm, newSection);CHECK_PETSC_ERROR(err);
+ err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+ err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) _localVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- // Reuse scatters in clone
- if (src._scatters.size() > 0) {
- const typename scatter_map_type::const_iterator scattersEnd = src._scatters.end();
- for (typename scatter_map_type::const_iterator s_iter=src._scatters.begin();
- s_iter != scattersEnd;
- ++s_iter) {
- ScatterInfo sinfo;
- sinfo.vector = 0;
- sinfo.scatter = 0;
- sinfo.scatterVec = 0;
+ // Reuse scatters in clone
+ if (src._scatters.size() > 0) {
+ const typename scatter_map_type::const_iterator scattersEnd = src._scatters.end();
+ for (typename scatter_map_type::const_iterator s_iter=src._scatters.begin();
+ s_iter != scattersEnd;
+ ++s_iter) {
+ ScatterInfo sinfo;
+ sinfo.vector = 0;
+ sinfo.scatter = 0;
+ sinfo.scatterVec = 0;
- // Copy DM
- sinfo.dm = s_iter->second.dm;
- err = PetscObjectReference((PetscObject) sinfo.dm);
- CHECK_PETSC_ERROR(err);
-
- // Copy scatter
- if (s_iter->second.scatter) {
- sinfo.scatter = s_iter->second.scatter;
- err = PetscObjectReference((PetscObject) sinfo.scatter);
+ // Copy DM
+ sinfo.dm = s_iter->second.dm;
+ err = PetscObjectReference((PetscObject) sinfo.dm);
CHECK_PETSC_ERROR(err);
- }
-
- // Create scatter Vec
- if (s_iter->second.scatterVec) {
- int blockSize = 1;
- err = VecGetBlockSize(s_iter->second.scatterVec, &blockSize);CHECK_PETSC_ERROR(err);
- if (_section->sizeWithBC() > 0) {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, _section->getStorageSize(),
- _section->restrictSpace(),
- &sinfo.scatterVec);
+
+ // Copy scatter
+ if (s_iter->second.scatter) {
+ sinfo.scatter = s_iter->second.scatter;
+ err = PetscObjectReference((PetscObject) sinfo.scatter);
CHECK_PETSC_ERROR(err);
+ }
+
+ // Create scatter Vec
+ sinfo.scatterVec = _localVec;
+ err = PetscObjectReference((PetscObject) sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+
+ // Create vector using sizes from source section
+ PetscInt vecLocalSize = 0;
+ PetscInt vecGlobalSize = 0, vecGlobalSize2 = 0;
+ err = VecGetLocalSize(s_iter->second.vector, &vecLocalSize);CHECK_PETSC_ERROR(err);
+ err = VecGetSize(s_iter->second.vector, &vecGlobalSize);CHECK_PETSC_ERROR(err);
+ err = VecGetSize(_globalVec, &vecGlobalSize2);CHECK_PETSC_ERROR(err);
+
+ if (vecGlobalSize != vecGlobalSize2) {
+ int blockSize = 1;
+ err = VecGetBlockSize(s_iter->second.vector, &blockSize);CHECK_PETSC_ERROR(err);
+ err = VecCreate(_mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(sinfo.vector, vecLocalSize, vecGlobalSize);CHECK_PETSC_ERROR(err);
+ err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
+ err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
} else {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, 0, PETSC_NULL,
- &sinfo.scatterVec);
+ sinfo.vector = _globalVec;
+ err = PetscObjectReference((PetscObject) sinfo.vector);
CHECK_PETSC_ERROR(err);
- } // else
- }
-
- // Create vector using sizes from source section
- int vecLocalSize = 0;
- int vecGlobalSize = 0, vecGlobalSize2 = 0;
- err = VecGetLocalSize(s_iter->second.vector, &vecLocalSize);CHECK_PETSC_ERROR(err);
- err = VecGetSize(s_iter->second.vector, &vecGlobalSize);CHECK_PETSC_ERROR(err);
- err = VecGetSize(_globalVec, &vecGlobalSize2);CHECK_PETSC_ERROR(err);
-
- if (vecGlobalSize != vecGlobalSize2) {
- int blockSize = 1;
- err = VecGetBlockSize(s_iter->second.vector, &blockSize);CHECK_PETSC_ERROR(err);
- err = VecCreate(_mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
- err = VecSetSizes(sinfo.vector, vecLocalSize, vecGlobalSize);CHECK_PETSC_ERROR(err);
- err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
- err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
- } else {
- sinfo.vector = _globalVec;
- err = PetscObjectReference((PetscObject) sinfo.vector);
- CHECK_PETSC_ERROR(err);
- }
- err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+ }
+ err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- _scatters[s_iter->first] = sinfo;
- } // for
- } // if
+ _scatters[s_iter->first] = sinfo;
+ } // for
+ } // if
logger.stagePop();
} // cloneSection
// ----------------------------------------------------------------------
-// Get DMPlex section.
-template<typename mesh_type, typename section_type>
-void
-pylith::topology::Field<mesh_type, section_type>::dmSection(PetscSection *s, Vec *v) const {
- PetscSection section;
- PetscInt size;
- PetscInt numNormalCells, numCohesiveCells, numNormalVertices, numShadowVertices, numLagrangeVertices;
- PetscErrorCode err;
-
- err = DMMeshConvertSection(_mesh.sieveMesh(), _section, §ion);CHECK_PETSC_ERROR(err);
- _mesh.getPointTypeSizes(&numNormalCells, &numCohesiveCells, &numNormalVertices, &numShadowVertices, &numLagrangeVertices);
- if (numNormalCells+numCohesiveCells+numNormalVertices+numShadowVertices+numLagrangeVertices > 0) {
- PetscInt numFields, numComp, pMax, pStart, pEnd, qStart, qEnd;
-
- err = DMPlexGetChart(_mesh.dmMesh(), PETSC_NULL, &pMax);CHECK_PETSC_ERROR(err);
- err = PetscSectionCreate(_mesh.sieveMesh()->comm(), s);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetNumFields(section, &numFields);CHECK_PETSC_ERROR(err);
- if (numFields) {
- err = PetscSectionSetNumFields(*s, numFields);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldComponents(section, f, &numComp);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldComponents(*s, f, numComp);CHECK_PETSC_ERROR(err);
- }
- }
- err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- if (pStart > 0) {
- qStart = pStart + numCohesiveCells;
- } else {
- qStart = pStart;
- }
- qEnd = PetscMin(pEnd + numCohesiveCells, pMax);
- err = PetscSectionSetChart(*s, qStart, qEnd);CHECK_PETSC_ERROR(err);
-#if 0
- if (pEnd - pStart != (numNormalCells + numCohesiveCells + numNormalVertices + numShadowVertices + numLagrangeVertices)) {
- PetscPrintf(PETSC_COMM_SELF, "numCells %d != totCells %d\n", pEnd - pStart, numNormalCells + numCohesiveCells + numNormalVertices + numShadowVertices + numLagrangeVertices);
- CHECK_PETSC_ERROR(PETSC_ERR_ARG_SIZ);
- }
-#endif
- /* The old-style point numbering: [normalCells, normalVertices, shadowVertices, lagrangeVertices, cohesiveCells]
- The new-style point numbering: [normalCells, cohesiveCells, normalVertices, shadowVertices, lagrangeVertices] */
- for(PetscInt p = pStart; p < numNormalCells; ++p) {
- PetscInt dof, fdof, cdof, cfdof, q = p;
-
- err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetDof(*s, q, dof);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldDof(section, p, f, &fdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldDof(*s, q, f, fdof);CHECK_PETSC_ERROR(err);
- }
- err = PetscSectionGetConstraintDof(section, p, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetConstraintDof(*s, q, cdof);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldConstraintDof(section, p, f, &cfdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldConstraintDof(*s, q, f, cfdof);CHECK_PETSC_ERROR(err);
- }
- }
- for(PetscInt p = PetscMax(pStart, numNormalCells); p < PetscMin(pEnd, numNormalCells+numNormalVertices+numShadowVertices+numLagrangeVertices); ++p) {
- PetscInt dof, fdof, cdof, cfdof, q = p + numCohesiveCells;
-
- err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetDof(*s, q, dof);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldDof(section, p, f, &fdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldDof(*s, q, f, fdof);CHECK_PETSC_ERROR(err);
- }
- err = PetscSectionGetConstraintDof(section, p, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetConstraintDof(*s, q, cdof);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldConstraintDof(section, p, f, &cfdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldConstraintDof(*s, q, f, cfdof);CHECK_PETSC_ERROR(err);
- }
- }
- for(PetscInt p = PetscMax(pStart, numNormalCells+numNormalVertices+numShadowVertices+numLagrangeVertices); p < pEnd; ++p) {
- PetscInt dof, fdof, cdof, cfdof, q = p - (numNormalVertices+numShadowVertices+numLagrangeVertices);
-
- err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetDof(*s, q, dof);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldDof(section, p, f, &fdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldDof(*s, q, f, fdof);CHECK_PETSC_ERROR(err);
- }
- err = PetscSectionGetConstraintDof(section, p, &cdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetConstraintDof(*s, q, cdof);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldConstraintDof(section, p, f, &cfdof);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldConstraintDof(*s, q, f, cfdof);CHECK_PETSC_ERROR(err);
- }
- }
- err = PetscSectionSetUp(*s);CHECK_PETSC_ERROR(err);
- for(PetscInt p = pStart; p < pStart+numNormalCells; ++p) {
- const PetscInt *cind, *cfind;
- PetscInt q = p;
-
- err = PetscSectionGetConstraintIndices(section, p, &cind);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetConstraintIndices(*s, q, cind);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldConstraintIndices(section, p, f, &cfind);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldConstraintIndices(*s, q, f, cfind);CHECK_PETSC_ERROR(err);
- }
- }
- for(PetscInt p = pStart+numNormalCells; p < pStart+numNormalCells+numNormalVertices+numShadowVertices+numLagrangeVertices; ++p) {
- const PetscInt *cind, *cfind;
- PetscInt q = p + numCohesiveCells;
-
- err = PetscSectionGetConstraintIndices(section, p, &cind);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetConstraintIndices(*s, q, cind);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldConstraintIndices(section, p, f, &cfind);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldConstraintIndices(*s, q, f, cfind);CHECK_PETSC_ERROR(err);
- }
- }
- for(PetscInt p = pStart+numNormalCells+numNormalVertices+numShadowVertices+numLagrangeVertices; p < pEnd; ++p) {
- const PetscInt *cind, *cfind;
- PetscInt q = p - (numNormalVertices+numShadowVertices+numLagrangeVertices);
-
- err = PetscSectionGetConstraintIndices(section, p, &cind);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetConstraintIndices(*s, q, cind);CHECK_PETSC_ERROR(err);
- for(PetscInt f = 0; f < numFields; ++f) {
- err = PetscSectionGetFieldConstraintIndices(section, p, f, &cfind);CHECK_PETSC_ERROR(err);
- err = PetscSectionSetFieldConstraintIndices(*s, q, f, cfind);CHECK_PETSC_ERROR(err);
- }
- }
- err = PetscSectionDestroy(§ion);CHECK_PETSC_ERROR(err);
- } else {
- *s = section;
- }
- err = PetscSectionGetStorageSize(*s, &size);CHECK_PETSC_ERROR(err);
- err = VecCreateMPIWithArray(_section->comm(), 1, size, PETSC_DETERMINE, _section->restrictSpace(), v);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) *v, _section->getName().c_str());CHECK_PETSC_ERROR(err);
-}
-
-// ----------------------------------------------------------------------
// Clear variables associated with section.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::clear(void)
+pylith::topology::Field<mesh_type>::clear(void)
{ // clear
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
deallocate();
- if (!_section.isNull())
- _section->clear();
-
_metadata["default"].scale = 1.0;
_metadata["default"].vectorFieldType = OTHER;
_metadata["default"].dimsOkay = false;
@@ -814,9 +517,9 @@
// ----------------------------------------------------------------------
// Allocate Sieve section.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::allocate(void)
+pylith::topology::Field<mesh_type>::allocate(void)
{ // allocate
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Field");
@@ -824,121 +527,80 @@
PetscErrorCode err;
if (_dm) {err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);}
- assert(!_section.isNull() || s);
+ assert(s);
- if (!_section.isNull()) {
- const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- sieveMesh->allocate(_section);
- }
- if (s) {
- err = PetscSectionSetUp(s);CHECK_PETSC_ERROR(err);
- err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
- err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject) _localVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- }
+ err = PetscSectionSetUp(s);CHECK_PETSC_ERROR(err);
+ err = DMCreateGlobalVector(_dm, &_globalVec);CHECK_PETSC_ERROR(err);
+ err = DMCreateLocalVector(_dm, &_localVec);CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) _globalVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
+ err = PetscObjectSetName((PetscObject) _localVec, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
logger.stagePop();
} // allocate
// ----------------------------------------------------------------------
// Zero section values (excluding constrained DOF).
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::zero(void)
+pylith::topology::Field<mesh_type>::zero(void)
{ // zero
- if (!_section.isNull())
- _section->zero(); // Does not zero BC.
- if (_localVec) {
- PetscSection section;
- PetscInt pStart, pEnd, maxDof = 0;
- PetscErrorCode err;
+ assert(_localVec);
+ PetscSection section;
+ PetscInt pStart, pEnd, maxDof = 0;
+ PetscErrorCode err;
- err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- if (pEnd > pStart) {err = PetscSectionGetMaxDof(section, &maxDof);CHECK_PETSC_ERROR(err);}
- scalar_array values(maxDof);
- values *= 0.0;
+ err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ if (pEnd > pStart) {err = PetscSectionGetMaxDof(section, &maxDof);CHECK_PETSC_ERROR(err);}
+ scalar_array values(maxDof);
+ values *= 0.0;
- for(PetscInt p = pStart; p < pEnd; ++p) {
- PetscInt dof;
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ PetscInt dof;
- err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
- if (dof > 0) {
- assert(dof <= maxDof);
- err = DMPlexVecSetClosure(_dm, section, _localVec, p, &values[0], INSERT_VALUES);CHECK_PETSC_ERROR(err);
- } // if
- } // for
- }
+ err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
+ if (dof > 0) {
+ assert(dof <= maxDof);
+ err = DMPlexVecSetClosure(_dm, section, _localVec, p, &values[0], INSERT_VALUES);CHECK_PETSC_ERROR(err);
+ } // if
+ } // for
} // zero
// ----------------------------------------------------------------------
// Zero section values (including constrained DOF).
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::zeroAll(void)
+pylith::topology::Field<mesh_type>::zeroAll(void)
{ // zeroAll
- if (!_section.isNull()) {
- const chart_type& chart = _section->getChart();
- const typename chart_type::const_iterator chartBegin = chart.begin();
- const typename chart_type::const_iterator chartEnd = chart.end();
-
- // Assume fiber dimension is uniform
- const int fiberDim = (chart.size() > 0) ?
- _section->getFiberDimension(*chartBegin) : 0;
- scalar_array values(fiberDim);
- values *= 0.0;
-
- for (typename chart_type::const_iterator c_iter = chartBegin;
- c_iter != chartEnd;
- ++c_iter) {
- if (_section->getFiberDimension(*c_iter) > 0) {
- assert(fiberDim == _section->getFiberDimension(*c_iter));
- _section->updatePointAll(*c_iter, &values[0]);
- } // if
- } // for
- } // if
- if (_localVec) {
- PetscErrorCode err = VecSet(_localVec, 0.0);CHECK_PETSC_ERROR(err);
- }
+ assert(_localVec);
+ PetscErrorCode err = VecSet(_localVec, 0.0);CHECK_PETSC_ERROR(err);
} // zeroAll
// ----------------------------------------------------------------------
// Complete section by assembling across processors.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::complete(void)
+pylith::topology::Field<mesh_type>::complete(void)
{ // complete
+ assert(_dm);
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Completion");
+ // Not sure if DMLocalToLocal() would work
+ PetscErrorCode err;
- const ALE::Obj<SieveMesh>& sieveMesh = _mesh.sieveMesh();
- assert(!sieveMesh.isNull());
-
- if (!_section.isNull())
- ALE::Completion::completeSectionAdd(sieveMesh->getSendOverlap(),
- sieveMesh->getRecvOverlap(),
- _section, _section);
-
- if (_dm) {
- // Not sure if DMLocalToLocal() would work
- PetscErrorCode err;
-
- err = VecSet(_globalVec, 0.0);CHECK_PETSC_ERROR(err);
- err = DMLocalToGlobalBegin(_dm, _localVec, ADD_VALUES, _globalVec);CHECK_PETSC_ERROR(err);
- err = DMLocalToGlobalEnd(_dm, _localVec, ADD_VALUES, _globalVec);CHECK_PETSC_ERROR(err);
- err = DMGlobalToLocalBegin(_dm, _globalVec, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
- err = DMGlobalToLocalEnd(_dm, _globalVec, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
- }
+ err = VecSet(_globalVec, 0.0);CHECK_PETSC_ERROR(err);
+ err = DMLocalToGlobalBegin(_dm, _localVec, ADD_VALUES, _globalVec);CHECK_PETSC_ERROR(err);
+ err = DMLocalToGlobalEnd(_dm, _localVec, ADD_VALUES, _globalVec);CHECK_PETSC_ERROR(err);
+ err = DMGlobalToLocalBegin(_dm, _globalVec, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
+ err = DMGlobalToLocalEnd(_dm, _globalVec, INSERT_VALUES, _localVec);CHECK_PETSC_ERROR(err);
logger.stagePop();
} // complete
// ----------------------------------------------------------------------
// Copy field values and metadata.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::copy(const Field& field)
+pylith::topology::Field<mesh_type>::copy(const Field& field)
{ // copy
// Check compatibility of sections
const int srcSize = field.chartSize();
@@ -965,77 +627,16 @@
} // if
assert(_localVec && field._localVec);
- if (!_section.isNull() && !field._section.isNull()) {
- // Copy values from field
- const chart_type& chart = _section->getChart();
- const typename chart_type::const_iterator chartBegin = chart.begin();
- const typename chart_type::const_iterator chartEnd = chart.end();
+ PetscErrorCode err = VecCopy(field._localVec, _localVec);CHECK_PETSC_ERROR(err);
- for (typename chart_type::const_iterator c_iter = chartBegin;
- c_iter != chartEnd;
- ++c_iter) {
- assert(field._section->getFiberDimension(*c_iter) ==
- _section->getFiberDimension(*c_iter));
- if (_section->getFiberDimension(*c_iter) > 0)
- _section->updatePointAll(*c_iter,
- field._section->restrictPoint(*c_iter));
- } // for
- } // if
-
- if (_localVec && field._localVec) {
- PetscErrorCode err = VecCopy(field._localVec, _localVec);CHECK_PETSC_ERROR(err);
- }
-
label(const_cast<Field&>(field)._metadata["default"].label.c_str()); // Update label
_metadata["default"].scale = const_cast<Field&>(field)._metadata["default"].scale;
} // copy
-// ----------------------------------------------------------------------
-// Copy field values.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::copy(const ALE::Obj<section_type>& osection)
+pylith::topology::Field<mesh_type>::copy(PetscSection osection, PetscInt field, PetscInt component, Vec ovec)
{ // copy
- // Check compatibility of sections
- const int srcSize = osection->getChart().size();
- const int dstSize = chartSize();
- if (srcSize != dstSize) {
- std::ostringstream msg;
-
- msg << "Cannot copy values from Sieve section "
- << _metadata["default"].label << "'. Sections are incompatible.\n"
- << " Source section:\n"
- << " size: " << srcSize << "\n"
- << " Destination section:\n"
- << " space dim: " << spaceDim() << "\n"
- << " vector field type: " << _metadata["default"].vectorFieldType << "\n"
- << " scale: " << _metadata["default"].scale << "\n"
- << " size: " << dstSize;
- throw std::runtime_error(msg.str());
- } // if
- assert( (_section.isNull() && osection.isNull()) ||
- (!_section.isNull() && !osection.isNull()) );
-
- if (!_section.isNull()) {
- // Copy values from field
- const chart_type& chart = _section->getChart();
- const typename chart_type::const_iterator chartEnd = chart.end();
-
- for (typename chart_type::const_iterator c_iter = chart.begin();
- c_iter != chartEnd;
- ++c_iter) {
- assert(osection->getFiberDimension(*c_iter) ==
- _section->getFiberDimension(*c_iter));
- if (osection->getFiberDimension(*c_iter))
- _section->updatePoint(*c_iter, osection->restrictPoint(*c_iter));
- } // for
- } // if
-} // copy
-
-template<typename mesh_type, typename section_type>
-void
-pylith::topology::Field<mesh_type, section_type>::copy(PetscSection osection, PetscInt field, PetscInt component, Vec ovec)
-{ // copy
PetscSection section;
PetscScalar *array, *oarray;
PetscInt numFields, numComp, pStart, pEnd, qStart, qEnd;
@@ -1107,9 +708,9 @@
// ----------------------------------------------------------------------
// Add two fields, storing the result in one of the fields.
-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)
+template<typename mesh_type>
+pylith::topology::Field<mesh_type>&
+pylith::topology::Field<mesh_type>::operator+=(const Field& field)
{ // operator+=
// Check compatibility of sections
const int srcSize = field.chartSize();
@@ -1135,44 +736,16 @@
<< " size: " << dstSize;
throw std::runtime_error(msg.str());
} // if
- assert( (_section.isNull() && field._section.isNull()) ||
- (!_section.isNull() && !field._section.isNull()) );
-
- if (!_section.isNull()) {
- // Add values from field
- const chart_type& chart = _section->getChart();
- const typename chart_type::const_iterator chartBegin = chart.begin();
- const typename chart_type::const_iterator chartEnd = chart.end();
-
- // Assume fiber dimension is uniform
- const int fiberDim = (chart.size() > 0) ?
- _section->getFiberDimension(*chartBegin) : 0;
- scalar_array values(fiberDim);
-
- for (typename chart_type::const_iterator c_iter = chartBegin;
- c_iter != chartEnd;
- ++c_iter) {
- if (field._section->getFiberDimension(*c_iter) > 0) {
- assert(fiberDim == field._section->getFiberDimension(*c_iter));
- assert(fiberDim == _section->getFiberDimension(*c_iter));
- field._section->restrictPoint(*c_iter, &values[0], values.size());
- _section->updatePointAllAdd(*c_iter, &values[0]);
- } // if
- } // for
- } // if
-
- if (_localVec && field._localVec) {
- PetscErrorCode err = VecAXPY(_localVec, 1.0, field._localVec);CHECK_PETSC_ERROR(err);
- }
-
+ assert(_localVec && field._localVec);
+ PetscErrorCode err = VecAXPY(_localVec, 1.0, field._localVec);CHECK_PETSC_ERROR(err);
return *this;
} // operator+=
// ----------------------------------------------------------------------
// Dimensionalize field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::dimensionalize(void) const
+pylith::topology::Field<mesh_type>::dimensionalize(void) const
{ // dimensionalize
if (!const_cast<Field*>(this)->_metadata["default"].dimsOkay) {
std::ostringstream msg;
@@ -1183,55 +756,32 @@
} // if
spatialdata::units::Nondimensional normalizer;
+ assert(_localVec);
+ PetscSection section;
+ PetscScalar *array;
+ PetscInt pStart, pEnd;
+ PetscErrorCode err;
- if (!_section.isNull()) {
- const chart_type& chart = _section->getChart();
- const typename chart_type::const_iterator chartEnd = chart.end();
+ err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ err = VecGetArray(_localVec, &array);CHECK_PETSC_ERROR(err);
+ for(PetscInt p = pStart; p < pEnd; ++p) {
+ PetscInt dof, off;
- // Assume fiber dimension is uniform
- const int fiberDim = (chart.size() > 0) ?
- _section->getFiberDimension(*chart.begin()) : 0;
- scalar_array values(fiberDim);
-
- for (typename chart_type::const_iterator c_iter = chart.begin();
- c_iter != chartEnd;
- ++c_iter)
- if (_section->getFiberDimension(*c_iter) > 0) {
- assert(fiberDim == _section->getFiberDimension(*c_iter));
-
- _section->restrictPoint(*c_iter, &values[0], values.size());
- normalizer.dimensionalize(&values[0], values.size(), const_cast<Field*>(this)->_metadata["default"].scale);
- _section->updatePointAll(*c_iter, &values[0]);
- } // if
- } // if
-
- if (_localVec) {
- PetscSection section;
- PetscScalar *array;
- PetscInt pStart, pEnd;
- PetscErrorCode err;
-
- err = DMGetDefaultSection(_dm, §ion);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetChart(section, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
- err = VecGetArray(_localVec, &array);CHECK_PETSC_ERROR(err);
- for(PetscInt p = pStart; p < pEnd; ++p) {
- PetscInt dof, off;
-
- err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
- err = PetscSectionGetOffset(section, p, &off);CHECK_PETSC_ERROR(err);
- if (dof) {
- normalizer.dimensionalize(&array[off], dof, const_cast<Field*>(this)->_metadata["default"].scale);
- }
+ err = PetscSectionGetDof(section, p, &dof);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetOffset(section, p, &off);CHECK_PETSC_ERROR(err);
+ if (dof) {
+ normalizer.dimensionalize(&array[off], dof, const_cast<Field*>(this)->_metadata["default"].scale);
}
- err = VecRestoreArray(_localVec, &array);CHECK_PETSC_ERROR(err);
}
+ err = VecRestoreArray(_localVec, &array);CHECK_PETSC_ERROR(err);
} // dimensionalize
// ----------------------------------------------------------------------
// Print field to standard out.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::view(const char* label) const
+pylith::topology::Field<mesh_type>::view(const char* label) const
{ // view
std::string vecFieldString;
switch(const_cast<Field*>(this)->_metadata["default"].vectorFieldType)
@@ -1271,8 +821,6 @@
<< " vector field type: " << vecFieldString << "\n"
<< " scale: " << const_cast<Field*>(this)->_metadata["default"].scale << "\n"
<< " dimensionalize flag: " << const_cast<Field*>(this)->_metadata["default"].dimsOkay << std::endl;
- if (!_section.isNull())
- _section->view(label);
if (_dm) {
PetscSection section;
PetscErrorCode err;
@@ -1288,10 +836,10 @@
// Create PETSc vector scatter for field. This is used to transfer
// information from the "global" PETSc vector view to the "local"
// Sieve section view.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
template<typename scatter_mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::createScatter(const scatter_mesh_type& mesh,
+pylith::topology::Field<mesh_type>::createScatter(const scatter_mesh_type& mesh,
const char* context)
{ // createScatter
assert(context);
@@ -1307,43 +855,6 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("GlobalOrder");
-
- if (!_section.isNull()) {
- assert(!mesh.sieveMesh().isNull());
- // Get global order (create if necessary).
- const std::string& orderLabel = _section->getName();
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, orderLabel, _section);
- assert(!order.isNull());
-
- // Create scatter
- err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, false, &sinfo.scatter);
- CHECK_PETSC_ERROR(err);
-
- // Create scatterVec
- const int blockSize = 1;
- if (_section->sizeWithBC() > 0) {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, _section->getStorageSize(),
- _section->restrictSpace(),
- &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
- } else {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, 0, PETSC_NULL,
- &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
- } // else
-
-#if 0
- // Create vector
- err = VecCreate(_mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- err = VecSetSizes(sinfo.vector, order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
- err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
- err = VecSetFromOptions(sinfo.vector);CHECK_PETSC_ERROR(err);
-#endif
- }
PetscInt localSize, globalSize;
err = PetscObjectReference((PetscObject) _dm);CHECK_PETSC_ERROR(err);
@@ -1365,10 +876,10 @@
// Sieve section view. The PETSc vector does not contain constrained
// DOF. Use createScatterWithBC() to include the constrained DOF in
// the PETSc vector.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
template<typename scatter_mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::createScatter(
+pylith::topology::Field<mesh_type>::createScatter(
const scatter_mesh_type& mesh,
const typename ALE::Obj<typename SieveMesh::numbering_type> numbering,
const char* context)
@@ -1389,48 +900,6 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("GlobalOrder");
-
- if (!_section.isNull()); {
- assert(!mesh.sieveMesh().isNull());
- // Get global order (create if necessary).
- const std::string& orderLabel =
- (strlen(context) > 0) ?
- _section->getName() + std::string("_") + std::string(context) :
- _section->getName();
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, orderLabel,
- numbering->getChart().begin(),
- numbering->getChart().end(),
- _section);
- assert(!order.isNull());
-
- // Create scatter
- err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, false, &sinfo.scatter);CHECK_PETSC_ERROR(err);
-
- // Create scatterVec
- const int blockSize = 1;
- if (_section->sizeWithBC() > 0) {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, _section->getStorageSize(),
- _section->restrictSpace(),
- &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
- } else {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, 0, PETSC_NULL,
- &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
- } // else
-
-#if 0
- // Create vector
- err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- err = VecSetSizes(sinfo.vector, order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
- err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
- err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);
-#endif
- }
PetscInt localSize, globalSize;
err = PetscObjectReference((PetscObject) _dm);CHECK_PETSC_ERROR(err);
@@ -1462,10 +931,10 @@
// Sieve section view. The PETSc vector does not contain constrained
// DOF. Use createScatterWithBC() to include the constrained DOF in
// the PETSc vector.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
template<typename scatter_mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::createScatterWithBC(
+pylith::topology::Field<mesh_type>::createScatterWithBC(
const scatter_mesh_type& mesh,
const char* context)
{ // createScatterWithBC
@@ -1483,43 +952,6 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("GlobalOrder");
- if (!_section.isNull()) {
- assert(!mesh.sieveMesh().isNull());
-
- // Get global order (create if necessary).
- const std::string& orderLabel = _section->getName();
- const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
- assert(!sieveMesh.isNull());
- const ALE::Obj<typename mesh_type::SieveMesh::order_type>& order =
- sieveMesh->getFactory()->getGlobalOrderWithBC(sieveMesh, orderLabel, _section);
- assert(!order.isNull());
-
- // Create scatter
- err = DMMeshCreateGlobalScatter(sieveMesh, _section, order, true, &sinfo.scatter);
- CHECK_PETSC_ERROR(err);
-
- // Create scatterVec
- const int blockSize = _getFiberDim();
- if (_section->sizeWithBC() > 0) {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, _section->getStorageSize(),
- _section->restrictSpace(),
- &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
- } else {
- err = VecCreateSeqWithArray(PETSC_COMM_SELF,
- blockSize, 0, PETSC_NULL,
- &sinfo.scatterVec);CHECK_PETSC_ERROR(err);
- } // else
-#if 0
- // Create vector
- err = VecCreate(mesh.comm(), &sinfo.vector);CHECK_PETSC_ERROR(err);
- err = PetscObjectSetName((PetscObject)sinfo.vector, _metadata["default"].label.c_str());CHECK_PETSC_ERROR(err);
- err = VecSetSizes(sinfo.vector, order->getLocalSize(), order->getGlobalSize());CHECK_PETSC_ERROR(err);
- err = VecSetBlockSize(sinfo.vector, blockSize);CHECK_PETSC_ERROR(err);
- err = VecSetFromOptions(sinfo.vector);CHECK_PETSC_ERROR(err);
-#endif
- }
-
PetscSection section, newSection, gsection;
PetscSF sf;
@@ -1548,10 +980,10 @@
// Sieve section view. The PETSc vector includes constrained DOF. Use
// createScatter() if constrained DOF should be omitted from the PETSc
// vector.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
template<typename scatter_mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::createScatterWithBC(
+pylith::topology::Field<mesh_type>::createScatterWithBC(
const scatter_mesh_type& mesh,
const std::string& labelName,
PetscInt labelValue,
@@ -1670,9 +1102,9 @@
// ----------------------------------------------------------------------
// Get PETSc vector associated with field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
PetscVec
-pylith::topology::Field<mesh_type, section_type>::vector(const char* context)
+pylith::topology::Field<mesh_type>::vector(const char* context)
{ // vector
std::ostringstream msg;
@@ -1682,9 +1114,9 @@
// ----------------------------------------------------------------------
// Get PETSc vector associated with field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
const PetscVec
-pylith::topology::Field<mesh_type, section_type>::vector(const char* context) const
+pylith::topology::Field<mesh_type>::vector(const char* context) const
{ // vector
std::ostringstream msg;
@@ -1695,9 +1127,9 @@
// ----------------------------------------------------------------------
// Scatter section information across processors to update the
// PETSc vector view of the field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::scatterSectionToVector(const char* context) const
+pylith::topology::Field<mesh_type>::scatterSectionToVector(const char* context) const
{ // scatterSectionToVector
assert(context);
@@ -1708,9 +1140,9 @@
// ----------------------------------------------------------------------
// Scatter section information across processors to update the
// PETSc vector view of the field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::scatterSectionToVector(const PetscVec vector,
+pylith::topology::Field<mesh_type>::scatterSectionToVector(const PetscVec vector,
const char* context) const
{ // scatterSectionToVector
assert(vector);
@@ -1734,9 +1166,9 @@
// ----------------------------------------------------------------------
// Scatter PETSc vector information across processors to update the
// section view of the field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::scatterVectorToSection(const char* context) const
+pylith::topology::Field<mesh_type>::scatterVectorToSection(const char* context) const
{ // scatterVectorToSection
assert(context);
@@ -1747,9 +1179,9 @@
// ----------------------------------------------------------------------
// Scatter PETSc vector information across processors to update the
// section view of the field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::scatterVectorToSection(const PetscVec vector,
+pylith::topology::Field<mesh_type>::scatterVectorToSection(const PetscVec vector,
const char* context) const
{ // scatterVectorToSection
assert(vector);
@@ -1772,52 +1204,30 @@
} // scatterVectorToSection
// ----------------------------------------------------------------------
-// Setup split field with all one space per spatial dimension.
-template<typename mesh_type, typename section_type>
-void
-pylith::topology::Field<mesh_type, section_type>::splitDefault(void)
-{ // splitDefault
- assert(!_section.isNull());
- const int spaceDim = _mesh.dimension();
- for (int iDim=0; iDim < spaceDim; ++iDim)
- _section->addSpace(); // displacements
-
- const chart_type& chart = _section->getChart();
-
- const typename chart_type::const_iterator chartBegin = chart.begin();
- const typename chart_type::const_iterator chartEnd = chart.end();
- for (int fibration=0; fibration < spaceDim; ++fibration)
- for (typename chart_type::const_iterator c_iter = chart.begin();
- c_iter != chartEnd;
- ++c_iter) {
- if (_section->getFiberDimension(*c_iter) > 0) {
- assert(spaceDim == _section->getFiberDimension(*c_iter));
- _section->setFiberDimension(*c_iter, 1, fibration);
- } // if
- } // for
-} // splitDefault
-
-// ----------------------------------------------------------------------
// Get fiber dimension associated with section (only works if fiber
// dimension is uniform).
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
int
-pylith::topology::Field<mesh_type, section_type>::_getFiberDim(void)
+pylith::topology::Field<mesh_type>::_getFiberDim(void)
{ // _getFiberDim
-
- int fiberDimLocal = (this->chartSize() > 0) ? _section->getFiberDimension(*_section->getChart().begin()) : 0;
- int fiberDim = 0;
- MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX,
- _mesh.comm());
+ assert(_dm);
+ PetscSection s;
+ PetscInt pStart, pEnd;
+ int fiberDimLocal, fiberDim = 0;
+ PetscErrorCode err;
+ err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
+ err = PetscSectionGetChart(s, &pStart, &pEnd);CHECK_PETSC_ERROR(err);
+ if (pEnd > pStart) {err = PetscSectionGetDof(s, pStart, &fiberDimLocal);CHECK_PETSC_ERROR(err);}
+ MPI_Allreduce(&fiberDimLocal, &fiberDim, 1, MPI_INT, MPI_MAX, _mesh.comm());
return fiberDim;
} // _getFiberDim
// ----------------------------------------------------------------------
// Get scatter for given context.
-template<typename mesh_type, typename section_type>
-typename pylith::topology::Field<mesh_type, section_type>::ScatterInfo&
-pylith::topology::Field<mesh_type, section_type>::_getScatter(const char* context,
+template<typename mesh_type>
+typename pylith::topology::Field<mesh_type>::ScatterInfo&
+pylith::topology::Field<mesh_type>::_getScatter(const char* context,
const bool createOk)
{ // _getScatter
assert(context);
@@ -1869,9 +1279,9 @@
// ----------------------------------------------------------------------
// Get scatter for given context.
-template<typename mesh_type, typename section_type>
-const typename pylith::topology::Field<mesh_type, section_type>::ScatterInfo&
-pylith::topology::Field<mesh_type, section_type>::_getScatter(const char* context) const
+template<typename mesh_type>
+const typename pylith::topology::Field<mesh_type>::ScatterInfo&
+pylith::topology::Field<mesh_type>::_getScatter(const char* context) const
{ // _getScatter
assert(context);
@@ -1887,18 +1297,18 @@
} // _getScatter
// Experimental
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::addField(const char *name, int numComponents)
+pylith::topology::Field<mesh_type>::addField(const char *name, int numComponents)
{
// Keep track of name/components until setup
_tmpFields[name] = numComponents;
_metadata[name] = _metadata["default"];
}
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::setupFields()
+pylith::topology::Field<mesh_type>::setupFields()
{
assert(_dm);
// Keep track of name/components until setup
@@ -1921,9 +1331,9 @@
#endif
}
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
void
-pylith::topology::Field<mesh_type, section_type>::updateDof(const char *name, const DomainEnum domain, int fiberDim)
+pylith::topology::Field<mesh_type>::updateDof(const char *name, const DomainEnum domain, int fiberDim)
{
PetscSection section;
PetscInt pStart, pEnd, f = 0;
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.hh 2013-02-22 19:56:50 UTC (rev 21384)
@@ -43,8 +43,7 @@
*
* Extends Sieve real general section by adding metadata.
*/
-template<typename mesh_type,
- typename section_type>
+template<typename mesh_type>
class pylith::topology::Field : public FieldBase
{ // Field
friend class TestFieldMesh; // unit testing
@@ -56,18 +55,9 @@
// Convenience typedefs
typedef mesh_type Mesh;
- typedef ALE::ISieveVisitor::RestrictVisitor<section_type> RestrictVisitor;
- typedef ALE::ISieveVisitor::UpdateAddVisitor<section_type> UpdateAddVisitor;
- typedef ALE::ISieveVisitor::UpdateAllVisitor<section_type> UpdateAllVisitor;
-
// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
private:
- // Convenience typedefs
- typedef typename mesh_type::SieveMesh SieveMesh;
- typedef typename SieveMesh::label_sequence label_sequence;
- typedef typename section_type::chart_type chart_type;
-
// PUBLIC MEMBERS ///////////////////////////////////////////////////////
public :
@@ -87,9 +77,7 @@
*
* @param mesh Finite-element mesh.
*/
- Field(const mesh_type& mesh,
- const ALE::Obj<section_type>& section,
- const Metadata& metadata);
+ Field(const mesh_type& mesh, const Metadata& metadata);
Field(const Field& src, const int fields[], int numFields);
@@ -100,18 +88,6 @@
/// Deallocate PETSc and local data structures.
void deallocate(void);
- /** Get Sieve section.
- *
- * @returns Sieve section.
- */
- const ALE::Obj<section_type>& section(void) const;
-
- /** Get DMPlex section.
- *
- * @returns DMPlex section.
- */
- void dmSection(PetscSection *s, Vec *v) const;
-
/** Get PetscSection.
*
* @returns PetscSection.
@@ -228,16 +204,16 @@
*/
void newSection(void);
- /** Create sieve section and set chart and fiber dimesion for
- * sequence of points.
+ /** Create sieve section and set chart and fiber dimesion for a list
+ * of points.
*
- * @param points Points over which to define section.
+ * @param pStart First point
+ * @param pEnd Upper bound for points
* @param dim Fiber dimension for section.
*
* @note Don't forget to call label(), especially if reusing a field.
*/
- void newSection(const ALE::Obj<label_sequence>& points,
- const int fiberDim);
+ void newSection(const PetscInt pStart, const PetscInt pEnd, const int fiberDim);
/** Create sieve section and set chart and fiber dimesion for a list
* of points.
@@ -248,8 +224,21 @@
* @note Don't forget to call label(), especially if reusing a field.
*/
void newSection(const int_array& points,
- const int fiberDim);
+ const int fiberDim);
+ /** Create sieve section and set chart and fiber dimesion for a list
+ * of points.
+ *
+ * @param points Points over which to define section.
+ * @param num The number of points
+ * @param dim Fiber dimension for section.
+ *
+ * @note Don't forget to call label(), especially if reusing a field.
+ */
+ void newSection(const PetscInt *points,
+ const PetscInt num,
+ const int fiberDim);
+
/** Create sieve section and set chart and fiber dimesion.
*
* @param domain Type of points over which to define section.
@@ -312,12 +301,6 @@
/** Copy field values.
*
- * @param field Field to copy.
- */
- void copy(const ALE::Obj<section_type>& field);
-
- /** Copy field values.
- *
* @param osection Field to copy.
* @param field Section field or -1
* @param component Section field component or -1
@@ -440,9 +423,6 @@
void scatterVectorToSection(const PetscVec vector,
const char* context ="") const;
- /// Setup split field with all entries set to a default space of 0.
- void splitDefault(void);
-
// PRIVATE STRUCTS //////////////////////////////////////////////////////
private :
@@ -499,7 +479,6 @@
map_type _metadata;
/* Old construction */
const mesh_type& _mesh; ///< Mesh associated with section.
- ALE::Obj<section_type> _section; ///< Real section with data.
scatter_map_type _scatters; ///< Collection of scatters.
/* New construction */
DM _dm; /* Holds the PetscSection */
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.icc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -22,29 +22,11 @@
#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
-// Get Sieve section.
-template<typename mesh_type, typename section_type>
-inline
-const ALE::Obj<section_type>&
-pylith::topology::Field<mesh_type, section_type>::section(void) const {
- std::ostringstream msg;
-
- char name[1024];
- PetscGetProgramName(name, 1024);
- msg << "Sieve Sections are no longer supported: " << const_cast<Field*>(this)->_metadata["default"].label << "\n"
- << " Destination section:\n"
- << " space dim: " << spaceDim() << "\n"
- << " vector field type: " << const_cast<Field*>(this)->_metadata["default"].vectorFieldType << "\n"
- << " scale: " << const_cast<Field*>(this)->_metadata["default"].scale;
- throw std::runtime_error(msg.str());
- return _section;
-}
-
// Get PetscSection.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
PetscSection
-pylith::topology::Field<mesh_type, section_type>::petscSection(void) const {
+pylith::topology::Field<mesh_type>::petscSection(void) const {
PetscSection s = PETSC_NULL;
if (_dm) {
PetscErrorCode err = DMGetDefaultSection(_dm, &s);CHECK_PETSC_ERROR(err);
@@ -53,106 +35,106 @@
}
// Get local vector.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
Vec
-pylith::topology::Field<mesh_type, section_type>::localVector(void) const {
+pylith::topology::Field<mesh_type>::localVector(void) const {
return _localVec;
}
// Get global vector.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
Vec
-pylith::topology::Field<mesh_type, section_type>::globalVector(void) const {
+pylith::topology::Field<mesh_type>::globalVector(void) const {
return _globalVec;
}
// Get mesh associated with field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
const
mesh_type&
-pylith::topology::Field<mesh_type, section_type>::mesh(void) const {
+pylith::topology::Field<mesh_type>::mesh(void) const {
return _mesh;
}
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
DM
-pylith::topology::Field<mesh_type, section_type>::dmMesh(void) const {
+pylith::topology::Field<mesh_type>::dmMesh(void) const {
return _dm;
}
// Get label for field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
const char*
-pylith::topology::Field<mesh_type, section_type>::label(void) const {
+pylith::topology::Field<mesh_type>::label(void) const {
return const_cast<Field*>(this)->_metadata["default"].label.c_str();
}
// Set vector field type
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
void
-pylith::topology::Field<mesh_type, section_type>::vectorFieldType(const VectorFieldEnum value) {
+pylith::topology::Field<mesh_type>::vectorFieldType(const VectorFieldEnum value) {
_metadata["default"].vectorFieldType = value;
}
// Set vector field type
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
void
-pylith::topology::Field<mesh_type, section_type>::vectorFieldType(const std::string& name, const VectorFieldEnum value) {
+pylith::topology::Field<mesh_type>::vectorFieldType(const std::string& name, const VectorFieldEnum value) {
const_cast<Field *>(this)->_metadata[name].vectorFieldType = value;
}
// Get vector field type
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
-typename pylith::topology::Field<mesh_type, section_type>::VectorFieldEnum
-pylith::topology::Field<mesh_type, section_type>::vectorFieldType(void) const {
+typename pylith::topology::Field<mesh_type>::VectorFieldEnum
+pylith::topology::Field<mesh_type>::vectorFieldType(void) const {
return const_cast<Field *>(this)->_metadata["default"].vectorFieldType;
}
// Set scale for dimensionalizing field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
void
-pylith::topology::Field<mesh_type, section_type>::scale(const PylithScalar value) {
+pylith::topology::Field<mesh_type>::scale(const PylithScalar value) {
_metadata["default"].scale = value;
}
// Set scale for dimensionalizing field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
void
-pylith::topology::Field<mesh_type, section_type>::scale(const std::string& name, const PylithScalar value) {
+pylith::topology::Field<mesh_type>::scale(const std::string& name, const PylithScalar value) {
_metadata[name].scale = value;
}
// Get scale for dimensionalizing field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
PylithScalar
-pylith::topology::Field<mesh_type, section_type>::scale(void) const {
+pylith::topology::Field<mesh_type>::scale(void) const {
return const_cast<Field *>(this)->_metadata["default"].scale;
}
// Set flag indicating whether it is okay to dimensionalize field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
void
-pylith::topology::Field<mesh_type, section_type>::addDimensionOkay(const bool value) {
+pylith::topology::Field<mesh_type>::addDimensionOkay(const bool value) {
_metadata["default"].dimsOkay = value;
}
// Set flag indicating whether it is okay to dimensionalize field.
-template<typename mesh_type, typename section_type>
+template<typename mesh_type>
inline
bool
-pylith::topology::Field<mesh_type, section_type>::addDimensionOkay(void) const {
+pylith::topology::Field<mesh_type>::addDimensionOkay(void) const {
return const_cast<Field*>(this)->_metadata["default"].dimsOkay;
}
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/FieldsNew.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/FieldsNew.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/FieldsNew.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -238,13 +238,10 @@
if (!f_iter->second.field) {
delete f_iter->second.field; f_iter->second.field = 0;
- assert(!_section.isNull());
- f_iter->second.field =
- new Field<mesh_type>(_mesh, _section->getFibration(fibration),
- f_iter->second.metadata);
+ f_iter->second.field = new Field<mesh_type>(_mesh, f_iter->second.metadata);
+ //new Field<mesh_type>(_mesh, _section->getFibration(fibration), f_iter->second.metadata);
assert(0 != f_iter->second.field);
} // if
-
return *f_iter->second.field;
} // get
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/topologyfwd.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/topologyfwd.hh 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/topologyfwd.hh 2013-02-22 19:56:50 UTC (rev 21384)
@@ -57,10 +57,8 @@
class MeshOps;
class FieldBase;
- template<typename mesh_type,
- typename section_type =ALE::IGeneralSection<pylith::SieveMesh::point_type, PylithScalar> > class Field;
+ template<typename mesh_type> class Field;
template<typename field_type> class Fields;
- template<typename mesh_type> class FieldsNew;
class SolutionFields;
Modified: short/3D/PyLith/trunk/modulesrc/faults/TractPerturbation.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/TractPerturbation.i 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/modulesrc/faults/TractPerturbation.i 2013-02-22 19:56:50 UTC (rev 21384)
@@ -51,7 +51,7 @@
*
* @returns Parameter fields.
*/
- const pylith::topology::FieldsNew<pylith::topology::SubMesh>* parameterFields(void) const;
+ const pylith::topology::Fields<topology::Field<pylith::topology::SubMesh> >* parameterFields(void) const;
/** Initialize slip time function.
*
Modified: short/3D/PyLith/trunk/modulesrc/friction/FrictionModel.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/friction/FrictionModel.i 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/modulesrc/friction/FrictionModel.i 2013-02-22 19:56:50 UTC (rev 21384)
@@ -122,7 +122,7 @@
*
* @returns Properties field.
*/
- const pylith::topology::FieldsNew<pylith::topology::SubMesh>& fieldsPropsStateVars() const;
+ const pylith::topology::Fields<pylith::topology::Field<pylith::topology::SubMesh> >& fieldsPropsStateVars() const;
/** Retrieve parameters for physical properties and state variables
* for vertex.
Modified: short/3D/PyLith/trunk/modulesrc/topology/Field.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/Field.i 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/modulesrc/topology/Field.i 2013-02-22 19:56:50 UTC (rev 21384)
@@ -254,9 +254,6 @@
void scatterVectorToSection(const PetscVec vector,
const char* context ="") const;
- /// Setup split field with all entries set to a default space of 0.
- void splitDefault(void);
-
}; // Field
} // topology
Modified: short/3D/PyLith/trunk/modulesrc/topology/topology.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/topology.i 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/modulesrc/topology/topology.i 2013-02-22 19:56:50 UTC (rev 21384)
@@ -27,7 +27,6 @@
#include "pylith/topology/FieldBase.hh"
#include "pylith/topology/Field.hh"
#include "pylith/topology/Fields.hh"
-#include "pylith/topology/FieldsNew.hh"
#include "pylith/topology/SolutionFields.hh"
#include "pylith/topology/Jacobian.hh"
#include "pylith/topology/Distributor.hh"
@@ -67,7 +66,6 @@
%include "FieldBase.i"
%include "Field.i"
%include "Fields.i"
-%include "FieldsNew.i"
%include "SolutionFields.i"
%include "Jacobian.i"
%include "Distributor.i"
@@ -81,9 +79,6 @@
%template(MeshFields) pylith::topology::Fields<pylith::topology::Field<pylith::topology::Mesh> >;
%template(SubMeshFields) pylith::topology::Fields<pylith::topology::Field<pylith::topology::SubMesh> >;
-%template(MeshFieldsNew) pylith::topology::FieldsNew<pylith::topology::Mesh>;
-%template(SubMeshFieldsNew) pylith::topology::FieldsNew<pylith::topology::SubMesh>;
-
%extend pylith::topology::Field<pylith::topology::Mesh> {
%template(createScatterMesh) createScatter<pylith::topology::Mesh>;
}
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestNeumann.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -46,12 +46,8 @@
// ----------------------------------------------------------------------
typedef pylith::topology::SubMesh::SieveMesh SieveMesh;
-typedef pylith::topology::SubMesh::RealSection RealSection;
typedef pylith::topology::SubMesh::SieveMesh SieveSubMesh;
-typedef pylith::topology::SubMesh::RealUniformSection SubRealUniformSection;
-typedef pylith::topology::Field<pylith::topology::SubMesh>::RestrictVisitor RestrictVisitor;
-
// ----------------------------------------------------------------------
namespace pylith {
namespace bc {
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -776,7 +776,6 @@
splitField.addField("multipliers", spaceDim);
splitField.setupFields();
splitField.newSection(disp, spaceDim);
- splitField.splitDefault();
fault.splitField(&splitField);
splitField.allocate();
splitField.zero();
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -33,8 +33,6 @@
// ----------------------------------------------------------------------
CPPUNIT_TEST_SUITE_REGISTRATION( pylith::feassemble::TestQuadrature );
-typedef pylith::topology::Field<pylith::topology::Mesh>::RestrictVisitor RestrictVisitor;
-
// ----------------------------------------------------------------------
// Test copy constuctor.
void
@@ -197,28 +195,16 @@
// Create mesh with test cell
topology::Mesh mesh(data.cellDim);
- const ALE::Obj<SieveMesh>& sieveMesh = mesh.sieveMesh();
- CPPUNIT_ASSERT(!sieveMesh.isNull());
- ALE::Obj<SieveMesh::sieve_type> sieve =
- new SieveMesh::sieve_type(mesh.comm());
- CPPUNIT_ASSERT(!sieve.isNull());
+ DM dmMesh;
+ PetscErrorCode err;
// Cells and vertices
- const bool interpolate = false;
- ALE::Obj<SieveFlexMesh::sieve_type> s =
- new SieveFlexMesh::sieve_type(sieve->comm(), sieve->debug());
-
- ALE::SieveBuilder<SieveFlexMesh>::buildTopology(s, cellDim, numCells,
- const_cast<int*>(data.cells),
- data.numVertices,
- interpolate, numBasis);
- std::map<SieveFlexMesh::point_type,SieveFlexMesh::point_type> renumbering;
- ALE::ISieveConverter::convertSieve(*s, *sieve, renumbering);
- sieveMesh->setSieve(sieve);
- sieveMesh->stratify();
- ALE::SieveBuilder<SieveMesh>::buildCoordinates(sieveMesh, spaceDim,
- data.vertices);
+ PetscBool interpolate = PETSC_FALSE;
+ err = DMPlexCreateFromCellList(mesh.comm(), cellDim, numCells, data.numVertices, numBasis, interpolate, const_cast<int*>(data.cells), data.vertices, &dmMesh);CHECK_PETSC_ERROR(err);
+ CPPUNIT_ASSERT(dmMesh);
+ mesh.setDMMesh(dmMesh);
+
// Setup quadrature and compute geometry
GeometryTri2D geometry;
Quadrature<topology::Mesh> quadrature;
@@ -230,34 +216,35 @@
data.quadWts, numQuadPts,
spaceDim);
- const ALE::Obj<SieveMesh::label_sequence>& cells = sieveMesh->heightStratum(0);
- CPPUNIT_ASSERT(!cells.isNull());
+ PetscInt cStart, cEnd;
+
+ err = DMPlexGetHeightStratum(dmMesh, 0, &cStart, &cEnd);CHECK_PETSC_ERROR(err);
quadrature.initializeGeometry();
#if defined(PRECOMPUTE_GEOMETRY)
quadrature.computeGeometry(mesh, cells);
#else
scalar_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<topology::Mesh::RealSection>& coordinates =
- sieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(), &coordinatesCell[0]);
+ PetscSection coordSection;
+ Vec coordVec;
+ err = DMPlexGetCoordinateSection(dmMesh, &coordSection);CHECK_PETSC_ERROR(err);
+ err = DMGetCoordinatesLocal(dmMesh, &coordVec);CHECK_PETSC_ERROR(err);
#endif
size_t size = 0;
// Check values from computeGeometry()
const PylithScalar tolerance = 1.0e-06;
- const SieveMesh::label_sequence::iterator cellsEnd = cells->end();
- for (SieveMesh::label_sequence::iterator c_iter=cells->begin();
- c_iter != cellsEnd;
- ++c_iter) {
+ for (PetscInt c = cStart; c < cEnd; ++c) {
#if defined(PRECOMPUTE_GEOMETRY)
- quadrature.retrieveGeometry(*c_iter);
+ quadrature.retrieveGeometry(c);
#else
- coordsVisitor.clear();
- sieveMesh->restrictClosure(*c_iter, coordsVisitor);
- quadrature.computeGeometry(coordinatesCell, *c_iter);
+ const PetscScalar *coords = PETSC_NULL;
+ PetscInt coordsSize;
+
+ err = DMPlexVecGetClosure(dmMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ for(PetscInt i = 0; i < coordsSize; ++i) {coordinatesCell[i] = coords[i];}
+ err = DMPlexVecRestoreClosure(dmMesh, coordSection, coordVec, c, &coordsSize, &coords);CHECK_PETSC_ERROR(err);
+ quadrature.computeGeometry(coordinatesCell, c);
#endif
const scalar_array& quadPts = quadrature.quadPts();
Modified: short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/unittests/libtests/friction/TestFrictionModel.cc 2013-02-22 19:56:50 UTC (rev 21384)
@@ -28,7 +28,6 @@
#include "pylith/topology/SubMesh.hh" // USES SubMesh
#include "pylith/topology/Field.hh" // USES Field
-#include "pylith/topology/FieldsNew.hh" // USES FieldsNew
#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
#include "pylith/faults/FaultCohesiveDyn.hh" // USES FaultCohesiveDyn
#include "pylith/feassemble/Quadrature.hh" // USES Quadrature
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/Makefile.am 2013-02-21 22:14:40 UTC (rev 21383)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/Makefile.am 2013-02-22 19:56:50 UTC (rev 21384)
@@ -35,7 +35,6 @@
TestFieldSubMesh.cc \
TestFieldsMesh.cc \
TestFieldsSubMesh.cc \
- TestFieldsNewMesh.cc \
TestSolutionFields.cc \
TestJacobian.cc \
TestRefineUniform.cc \
More information about the CIG-COMMITS
mailing list