[cig-commits] r19124 - in short/3D/PyLith/branches/v1.6-stable/libsrc/pylith: bc faults feassemble meshio
brad at geodynamics.org
brad at geodynamics.org
Fri Oct 28 15:22:40 PDT 2011
Author: brad
Date: 2011-10-28 15:22:39 -0700 (Fri, 28 Oct 2011)
New Revision: 19124
Modified:
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/bc/AbsorbingDampers.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityExplicit.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityImplicit.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
Log:
Fixes to ensure indices and closure visitors have positive size (getMaxConeSize could return -1 when mesh is empty). Prevents bad_alloc exceptions being thrown.
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/bc/AbsorbingDampers.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/bc/AbsorbingDampers.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -587,9 +587,14 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
solutionSection);
assert(!globalOrder.isNull());
+
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, sieveMesh->depth()))));
IndicesVisitor jacobianVisitor(*solutionSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ closureSize*spaceDim);
// Get sparse matrix
const PetscMat jacobianMat = jacobian->matrix();
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -1660,8 +1660,11 @@
assert(!globalOrderDomain.isNull());
const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
assert(!sieve.isNull());
- ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> ncV(*sieve,
- (size_t) pow(sieve->getMaxConeSize(), std::max(0, sieveMesh->depth())));
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, sieveMesh->depth()))));
+ ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type>
+ ncV(*sieve, closureSize);
int_array indicesGlobal(subnrows);
// Get fault Sieve mesh
@@ -1682,9 +1685,7 @@
assert(!globalOrderFault.isNull());
// We would need to request unique points here if we had an interpolated mesh
IndicesVisitor jacobianFaultVisitor(*solutionFaultSection,
- *globalOrderFault,
- (int) pow(faultSieveMesh->getSieve()->getMaxConeSize(),
- faultSieveMesh->depth())*spaceDim);
+ *globalOrderFault, closureSize*spaceDim);
const int iCone = (negativeSide) ? 0 : 1;
for (SieveMesh::label_sequence::iterator c_iter=cellsCohesiveBegin;
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -1087,10 +1087,12 @@
const ALE::Obj<SieveMesh::sieve_type>& sieve = mesh.sieveMesh()->getSieve();
assert(!sieve.isNull());
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, faultSieveMesh->depth()))));
+ ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type>
+ ncV(*sieve, closureSize);
- ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> ncV(*sieve,
- (size_t) pow(sieve->getMaxConeSize(), std::max(0, sieveMesh->depth())));
-
for (SieveMesh::label_sequence::iterator c_iter=cellsBegin;
c_iter != cellsEnd;
++c_iter, ++f_iter) {
@@ -1202,11 +1204,10 @@
const double_array& quadWts = _quadrature->quadWts();
double_array jacobian(jacobianSize);
double jacobianDet = 0;
- double_array orientationVertex(orientationSize);
- double_array coordinatesCell(numBasis * spaceDim);
double_array refCoordsVertex(cohesiveDim);
// Allocate orientation field.
+ double_array orientationVertex(orientationSize);
_fields->add("orientation", "orientation");
topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
const topology::Field<topology::SubMesh>& dispRel =
@@ -1232,6 +1233,7 @@
// Compute orientation of fault at constraint vertices
// Get section containing coordinates of vertices
+ double_array coordinatesCell(numBasis * spaceDim);
const ALE::Obj<RealSection>& coordinatesSection =
faultSieveMesh->getRealSection("coordinates");
assert(!coordinatesSection.isNull());
@@ -1248,10 +1250,11 @@
const ALE::Obj<SieveSubMesh::sieve_type>& sieve = faultSieveMesh->getSieve();
assert(!sieve.isNull());
-
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, faultSieveMesh->depth()))));
ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type>
- ncV(*sieve, (size_t) pow(sieve->getMaxConeSize(), std::max(0,
- faultSieveMesh->depth())));
+ ncV(*sieve, closureSize);
for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin; c_iter
!= cellsEnd; ++c_iter) {
@@ -1263,9 +1266,9 @@
ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve,
*c_iter, ncV);
const int coneSize = ncV.getSize();
- const SieveSubMesh::point_type *cone = ncV.getPoints();
+ const SieveSubMesh::point_type* cone = ncV.getPoints();
- for (int v = 0; v < coneSize; ++v) {
+ for (int v=0; v < coneSize; ++v) {
// Compute Jacobian and determinant of Jacobian at vertex
memcpy(&refCoordsVertex[0], &verticesRef[v * cohesiveDim], cohesiveDim
* sizeof(double));
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityExplicit.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityExplicit.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -723,9 +723,13 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solnSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, sieveMesh->depth()))));
IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ closureSize*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -647,9 +647,13 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solnSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, sieveMesh->depth()))));
IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ closureSize*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityImplicit.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityImplicit.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -410,10 +410,13 @@
dispTSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- IndicesVisitor jacobianVisitor(*dispTSection,
- *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, sieveMesh->depth()))));
+ IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ closureSize*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -405,10 +405,13 @@
dispTSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- IndicesVisitor jacobianVisitor(*dispTSection,
- *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, sieveMesh->depth()))));
+ IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ closureSize*spaceDim);
double_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -98,8 +98,9 @@
_tstampIndex = 0;
const int rank = sieveMesh->commRank();
const int localSize = (!rank) ? 1 : 0;
- err = VecCreateMPI(mesh.comm(), localSize, PETSC_DETERMINE, &_tstamp);
+ err = VecCreateMPI(mesh.comm(), localSize, 1, &_tstamp);
CHECK_PETSC_ERROR(err);
+ assert(_tstamp);
err = VecSetBlockSize(_tstamp, 1); CHECK_PETSC_ERROR(err);
err = PetscObjectSetName((PetscObject) _tstamp, "time");
@@ -170,10 +171,12 @@
const Obj<sieve_type>& sieve = sieveMesh->getSieve();
assert(!sieve.isNull());
- ALE::ISieveVisitor::NConeRetriever<sieve_type>
- ncV(*sieve, (size_t) pow((double) sieve->getMaxConeSize(),
- std::max(0, sieveMesh->depth())));
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, sieveMesh->depth()))));
+ ALE::ISieveVisitor::NConeRetriever<sieve_type> ncV(*sieve, closureSize);
+
int k = 0;
const typename label_sequence::const_iterator cellsEnd = cells->end();
for (typename label_sequence::iterator c_iter=cells->begin();
@@ -238,6 +241,7 @@
PetscErrorCode err = 0;
err = PetscViewerDestroy(&_viewer); CHECK_PETSC_ERROR(err);
err = VecDestroy(&_tstamp); CHECK_PETSC_ERROR(err);
+ assert(!_tstamp);
_timesteps.clear();
_tstampIndex = 0;
@@ -468,6 +472,7 @@
const double t,
const int rank)
{ // _writeTimeStamp
+ assert(_tstamp);
PetscErrorCode err = 0;
if (0 == rank) {
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2011-10-28 14:44:10 UTC (rev 19123)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2011-10-28 22:22:39 UTC (rev 19124)
@@ -186,10 +186,11 @@
const Obj<sieve_type>& sieve = sieveMesh->getSieve();
assert(!sieve.isNull());
- ALE::ISieveVisitor::NConeRetriever<sieve_type>
- ncV(*sieve, (size_t) pow((double) sieve->getMaxConeSize(),
- std::max(0, sieveMesh->depth())));
-
+ const int closureSize =
+ std::max(0, int(pow(sieve->getMaxConeSize(),
+ std::max(0, sieveMesh->depth()))));
+ ALE::ISieveVisitor::NConeRetriever<sieve_type> ncV(*sieve, closureSize);
+
int k = 0;
const typename label_sequence::const_iterator cellsEnd = cells->end();
for (typename label_sequence::iterator c_iter=cells->begin();
More information about the CIG-COMMITS
mailing list