[cig-commits] r19128 - in short/3D/PyLith/trunk/libsrc/pylith: bc faults feassemble meshio
brad at geodynamics.org
brad at geodynamics.org
Fri Oct 28 17:36:17 PDT 2011
Author: brad
Date: 2011-10-28 17:36:16 -0700 (Fri, 28 Oct 2011)
New Revision: 19128
Modified:
short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.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/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnSubset.cc
Log:
Merge from stable.
Modified: short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/bc/AbsorbingDampers.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -590,9 +590,13 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
solutionSection);
assert(!globalOrder.isNull());
+
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const int closureSize =
+ int(pow(sieve->getMaxConeSize(), 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/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -1097,7 +1097,7 @@
areaVertex * dLagrangeTpdtVertexGlobal[iDim] / jacobianVertexP[iDim];
// Set increment in relative displacement.
- dispRelVertex[iDim] = -areaVertex * dLagrangeTpdtVertexGlobal[iDim] /
+ dispRelVertex[iDim] = -areaVertex * 2.0*dLagrangeTpdtVertexGlobal[iDim] /
(jacobianVertexN[iDim] + jacobianVertexP[iDim]);
// Update increment in Lagrange multiplier.
@@ -1406,8 +1406,6 @@
// Close properties database
_dbInitialTract->close();
- initialTractions.complete(); // Assemble contributions
-
//initialTractions.view("INITIAL TRACTIONS"); // DEBUGGING
} // _setupInitialTractions
@@ -1662,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 =
+ int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
+ assert(closureSize >= 0);
+ ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type>
+ ncV(*sieve, closureSize);
int_array indicesGlobal(subnrows);
// Get fault Sieve mesh
@@ -1684,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/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -1087,10 +1087,12 @@
const ALE::Obj<SieveMesh::sieve_type>& sieve = mesh.sieveMesh()->getSieve();
assert(!sieve.isNull());
+ const int closureSize =
+ int(pow(sieve->getMaxConeSize(), faultSieveMesh->depth()));
+ assert(closureSize >= 0);
+ 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) {
@@ -1204,11 +1206,10 @@
const scalar_array& quadWts = _quadrature->quadWts();
scalar_array jacobian(jacobianSize);
PylithScalar jacobianDet = 0;
- scalar_array orientationVertex(orientationSize);
- scalar_array coordinatesCell(numBasis * spaceDim);
scalar_array refCoordsVertex(cohesiveDim);
// Allocate orientation field.
+ scalar_array orientationVertex(orientationSize);
_fields->add("orientation", "orientation");
topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
const topology::Field<topology::SubMesh>& dispRel =
@@ -1234,6 +1235,7 @@
// Compute orientation of fault at constraint vertices
// Get section containing coordinates of vertices
+ scalar_array coordinatesCell(numBasis * spaceDim);
const ALE::Obj<RealSection>& coordinatesSection =
faultSieveMesh->getRealSection("coordinates");
assert(!coordinatesSection.isNull());
@@ -1250,10 +1252,11 @@
const ALE::Obj<SieveSubMesh::sieve_type>& sieve = faultSieveMesh->getSieve();
assert(!sieve.isNull());
-
+ const int closureSize =
+ int(pow(sieve->getMaxConeSize(), faultSieveMesh->depth()));
+ assert(closureSize >= 0);
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) {
@@ -1265,9 +1268,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(PylithScalar));
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -725,9 +725,13 @@
sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solnSection);
assert(!globalOrder.isNull());
// We would need to request unique points here if we had an interpolated mesh
- IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ const ALE::Obj<SieveMesh::sieve_type>& sieve = sieveMesh->getSieve();
+ assert(!sieve.isNull());
+ const int closureSize =
+ int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
+ assert(closureSize >= 0);
+ IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
+ closureSize*spaceDim);
scalar_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -649,9 +649,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 =
+ int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
+ assert(closureSize >= 0);
IndicesVisitor jacobianVisitor(*solnSection, *globalOrder,
- (int) pow(sieveMesh->getSieve()->getMaxConeSize(),
- sieveMesh->depth())*spaceDim);
+ closureSize*spaceDim);
scalar_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicit.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -411,10 +411,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 =
+ int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
+ assert(closureSize >= 0);
+ IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ closureSize*spaceDim);
scalar_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityImplicitLgDeform.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -406,10 +406,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 =
+ int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
+ assert(closureSize >= 0);
+ IndicesVisitor jacobianVisitor(*dispTSection, *globalOrder,
+ closureSize*spaceDim);
scalar_array coordinatesCell(numBasis*spaceDim);
const ALE::Obj<RealSection>& coordinates =
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -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,9 +171,10 @@
const Obj<sieve_type>& sieve = sieveMesh->getSieve();
assert(!sieve.isNull());
- ALE::ISieveVisitor::NConeRetriever<sieve_type>
- ncV(*sieve, (size_t) pow((PylithScalar) sieve->getMaxConeSize(),
- std::max(0, sieveMesh->depth())));
+ const int closureSize =
+ int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
+ assert(closureSize >= 0);
+ ALE::ISieveVisitor::NConeRetriever<sieve_type> ncV(*sieve, closureSize);
int k = 0;
const typename label_sequence::const_iterator cellsEnd = cells->end();
@@ -238,6 +240,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 +471,7 @@
const PylithScalar t,
const int rank)
{ // _writeTimeStamp
+ assert(_tstamp);
PetscErrorCode err = 0;
if (0 == rank) {
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -189,10 +189,11 @@
const Obj<sieve_type>& sieve = sieveMesh->getSieve();
assert(!sieve.isNull());
- ALE::ISieveVisitor::NConeRetriever<sieve_type>
- ncV(*sieve, (size_t) pow((PylithScalar) sieve->getMaxConeSize(),
- std::max(0, sieveMesh->depth())));
-
+ const int closureSize =
+ int(pow(sieve->getMaxConeSize(), sieveMesh->depth()));
+ assert(closureSize >= 0);
+ 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();
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputManager.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -118,8 +118,17 @@
const char* label,
const int labelId)
{ // open
- assert(0 != _writer);
+ if (!_writer) {
+ std::ostringstream msg;
+ if (label) {
+ msg << "Writer for output manager for " << label << " not set.";
+ throw std::runtime_error(msg.str());
+ } else {
+ throw std::runtime_error("Writer for output manager not set.");
+ } // if/else
+ } // if
+ assert(_writer);
_writer->open(mesh, numTimeSteps, label, labelId);
} // open
@@ -129,7 +138,7 @@
void
pylith::meshio::OutputManager<mesh_type, field_type>::close(void)
{ // close
- assert(0 != _writer);
+ assert(_writer);
_writer->close();
} // close
@@ -143,7 +152,7 @@
const char* label,
const int labelId)
{ // openTimeStep
- assert(0 != _writer);
+ assert(_writer);
_writer->openTimeStep(t, mesh, label, labelId);
} // openTimeStep
@@ -153,7 +162,7 @@
void
pylith::meshio::OutputManager<mesh_type, field_type>::closeTimeStep(void)
{ // closeTimeStep
- assert(0 != _writer);
+ assert(_writer);
_writer->closeTimeStep();
} // closeTimeStep
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnSubset.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnSubset.cc 2011-10-29 00:26:58 UTC (rev 19127)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/OutputSolnSubset.cc 2011-10-29 00:36:16 UTC (rev 19128)
@@ -77,7 +77,7 @@
pylith::meshio::OutputSolnSubset::subdomainMesh(const topology::Mesh& mesh)
{ // subdomainMesh
delete _submesh; _submesh = new topology::SubMesh(mesh, _label.c_str());
- assert(0 != _submesh);
+ assert(_submesh);
return *_submesh;
} // subdomainMesh
More information about the CIG-COMMITS
mailing list