[cig-commits] r18797 - in short/3D/PyLith/trunk: examples/bar_shearwave/tri3/tmp libsrc/pylith/faults libsrc/pylith/feassemble libsrc/pylith/materials libsrc/pylith/meshio libsrc/pylith/topology pylith/problems

brad at geodynamics.org brad at geodynamics.org
Fri Jul 22 13:35:08 PDT 2011


Author: brad
Date: 2011-07-22 13:35:08 -0700 (Fri, 22 Jul 2011)
New Revision: 18797

Modified:
   short/3D/PyLith/trunk/examples/bar_shearwave/tri3/tmp/pylithapp.cfg
   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/feassemble/IntegratorElasticity.cc
   short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
   short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriter.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/topology/Distributor.cc
   short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
   short/3D/PyLith/trunk/pylith/problems/ExplicitLumped.py
   short/3D/PyLith/trunk/pylith/problems/Formulation.py
Log:
Merge from stable.

Modified: short/3D/PyLith/trunk/examples/bar_shearwave/tri3/tmp/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/bar_shearwave/tri3/tmp/pylithapp.cfg	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/examples/bar_shearwave/tri3/tmp/pylithapp.cfg	2011-07-22 20:35:08 UTC (rev 18797)
@@ -96,6 +96,7 @@
 id = 1
 
 # Spatial database with physical properties for elastic material
+db_properties.label = Elastic properties
 db_properties.iohandler.filename = ../matprops.spatialdb
 
 # Set the basis functions and quadrature:
@@ -182,12 +183,15 @@
 [pylithapp.timedependent.interfaces.fault.eq_srcs.rupture.slip_function]
 
 # Database specifying the final slip.
+slip.label = Final slip
 slip.iohandler.filename = ../shearwave_slip.spatialdb
 
 # Database specifying rise time.
+rise_time.label = Rise time
 rise_time.iohandler.filename = ../shearwave_risetime.spatialdb
 
 # Database specifying time at which slip begins at each point.
+slip_time.label = Slip time
 slip_time.iohandler.filename = ../shearwave_sliptime.spatialdb
 
 # ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -1035,9 +1035,9 @@
     _allocateBufferVectorField();
     topology::Field<topology::SubMesh>& buffer =
         _fields->get("buffer (vector)");
-    buffer.copy(dirSection);
     buffer.label("strike_dir");
     buffer.scale(1.0);
+    buffer.copy(dirSection);
     return buffer;
 
   } else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
@@ -1049,9 +1049,9 @@
     _allocateBufferVectorField();
     topology::Field<topology::SubMesh>& buffer =
         _fields->get("buffer (vector)");
-    buffer.copy(dirSection);
     buffer.label("dip_dir");
     buffer.scale(1.0);
+    buffer.copy(dirSection);
     return buffer;
 
   } else if (0 == strcasecmp("normal_dir", name)) {
@@ -1065,9 +1065,9 @@
     _allocateBufferVectorField();
     topology::Field<topology::SubMesh>& buffer =
         _fields->get("buffer (vector)");
-    buffer.copy(dirSection);
     buffer.label("normal_dir");
     buffer.scale(1.0);
+    buffer.copy(dirSection);
     return buffer;
 
   } else if (0 == strcasecmp("initial_traction", name)) {
@@ -1097,11 +1097,12 @@
     throw std::runtime_error(msg.str());
   } // else
 
-
   // Satisfy return values
   assert(0 != _fields);
   const topology::Field<topology::SubMesh>& buffer = _fields->get(
     "buffer (vector)");
+  throw std::logic_error("Internal error.");
+
   return buffer;
 } // vertexField
 
@@ -1964,13 +1965,18 @@
     if (tractionShearMag > frictionStress || slipRateMag > 0.0) {
       // traction is limited by friction, so have sliding OR
       // friction exceeds traction due to overshoot in slip
-      
-      // Update traction increment based on value required to stick
-      // versus friction
-      const double dlp = -(tractionShearMag - frictionStress) * area *
-	tractionTpdt[0] / tractionShearMag;
-      (*dLagrangeTpdt)[0] = dlp;
-      (*dLagrangeTpdt)[1] = 0.0;
+
+      if (tractionShearMag > 0.0) {
+	// Update traction increment based on value required to stick
+	// versus friction
+	const double dlp = -(tractionShearMag - frictionStress) * area *
+	  tractionTpdt[0] / tractionShearMag;
+	(*dLagrangeTpdt)[0] = dlp;
+	(*dLagrangeTpdt)[1] = 0.0;
+      } else {
+	(*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
+	(*dLagrangeTpdt)[1] = 0.0;
+      } // if/else
     } else {
       // friction exceeds value necessary to stick
       // no changes to solution
@@ -2014,16 +2020,22 @@
       // traction is limited by friction, so have sliding OR
       // friction exceeds traction due to overshoot in slip
       
-      // Update traction increment based on value required to stick
-      // versus friction
-      const double dlp = -(tractionShearMag - frictionStress) * area *
-	tractionTpdt[0] / tractionShearMag;
-      const double dlq = -(tractionShearMag - frictionStress) * area *
-	tractionTpdt[1] / tractionShearMag;
-
-      (*dLagrangeTpdt)[0] = dlp;
-      (*dLagrangeTpdt)[1] = dlq;
-      (*dLagrangeTpdt)[2] = 0.0;
+      if (tractionShearMag > 0.0) {
+	// Update traction increment based on value required to stick
+	// versus friction
+	const double dlp = -(tractionShearMag - frictionStress) * area *
+	  tractionTpdt[0] / tractionShearMag;
+	const double dlq = -(tractionShearMag - frictionStress) * area *
+	  tractionTpdt[1] / tractionShearMag;
+	
+	(*dLagrangeTpdt)[0] = dlp;
+	(*dLagrangeTpdt)[1] = dlq;
+	(*dLagrangeTpdt)[2] = 0.0;
+      } else {
+	(*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
+	(*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
+	(*dLagrangeTpdt)[2] = 0.0;
+      } // if/else	
       
     } else {
       // else friction exceeds value necessary, so stick

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -1826,8 +1826,8 @@
   topology::Field<topology::SubMesh>& area = _fields->get("area");
   const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
   area.newSection(slip, 1);
+  area.allocate();
   area.vectorFieldType(topology::FieldBase::SCALAR);
-  area.allocate();
   area.zero();
   const ALE::Obj<RealSection>& areaSection = area.section();
   assert(!areaSection.isNull());
@@ -1976,6 +1976,7 @@
   const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
   buffer.cloneSection(slip);
   buffer.zero();
+  assert(buffer.vectorFieldType() == topology::FieldBase::VECTOR);
 
   logger.stagePop();
 } // _allocateBufferVectorField
@@ -1999,6 +2000,7 @@
   const topology::Field<topology::SubMesh>& area = _fields->get("area");
   buffer.cloneSection(area);
   buffer.zero();
+  assert(buffer.vectorFieldType() == topology::FieldBase::SCALAR);
 
   logger.stagePop();
 } // _allocateBufferScalarField

Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveTract.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -302,10 +302,10 @@
     // Create section to hold initial tractions.
     _fields->add("initial traction", "initial_traction");
     topology::Field<topology::SubMesh>& traction = _fields->get("initial traction");
+    traction.newSection(topology::FieldBase::CELLS_FIELD, numQuadPts*spaceDim);
+    traction.allocate();
     traction.scale(pressureScale);
     traction.vectorFieldType(topology::FieldBase::MULTI_VECTOR);
-    traction.newSection(topology::FieldBase::CELLS_FIELD, numQuadPts*spaceDim);
-    traction.allocate();
     const ALE::Obj<RealSection>& tractionSection = traction.section();
     assert(!tractionSection.isNull());
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -400,6 +400,7 @@
   
   // Return tensor section to satisfy member function definition. Code
   // should never get here.
+  throw std::logic_error("Internal error.");
   topology::Field<topology::Mesh>& buffer = 
     _outputFields->get("buffer (tensor)");    
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -315,7 +315,8 @@
 // ----------------------------------------------------------------------
 // Get physical property or state variable field.
 void
-pylith::materials::Material::getField(topology::Field<topology::Mesh> *field, const char* name) const
+pylith::materials::Material::getField(topology::Field<topology::Mesh> *field,
+				      const char* name) const
 { // getField
   // Logging of allocation is handled by getField() caller since it
   // manages the memory for the field argument.

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -128,11 +128,6 @@
     _fieldAvg->newSection(fieldIn, fiberDim);
     _fieldAvg->allocate();
   } else if (_fieldAvg->sectionSize() != cells->size()*fiberDim) {
-#if 1 // :BUG: Avoid memory leak in section->clear()
-    _fieldAvg->clear();
-#else
-    delete _fieldAvg; _fieldAvg = new field_type(fieldIn.mesh());
-#endif
     _fieldAvg->newSection(fieldIn, fiberDim);
     _fieldAvg->allocate();
   } // else

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriter.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriter.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriter.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -56,8 +56,13 @@
 
   const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh = mesh.sieveMesh();
   assert(!sieveMesh.isNull());
-  _context = std::string("output_") + sieveMesh->getName() + 
-    ((label) ? std::string("_") + std::string(label) : std::string(""));
+
+  ostringstream s;
+  s << "output_"
+    << sieveMesh->getName();
+  if (label)
+    s << "_" << label << labelId;
+  _context = s.str();
 } // open
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -201,14 +201,12 @@
     err = VecDestroy(&elemVec); CHECK_PETSC_ERROR(err);
     delete[] tmpVertices; tmpVertices = 0;
 
-    if (!rank) {
-      hid_t h5 = -1;
-      err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
-      assert(h5 >= 0);
-      const int cellDim = mesh.dimension();
-      HDF5::writeAttribute(h5, "/topology/cells", "cell_dim", (void*)&cellDim,
-			   H5T_NATIVE_INT);
-    } // if
+    hid_t h5 = -1;
+    err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
+    assert(h5 >= 0);
+    const int cellDim = mesh.dimension();
+    HDF5::writeAttribute(h5, "/topology/cells", "cell_dim", (void*)&cellDim,
+			 H5T_NATIVE_INT);
   } catch (const std::exception& err) {
     std::ostringstream msg;
     msg << "Error while opening HDF5 file "
@@ -235,12 +233,16 @@
   _timesteps.clear();
   _tstampIndex = 0;
 
-  Xdmf metafile;
-  const std::string& hdf5filename = _hdf5Filename();
-  const int indexExt = hdf5filename.find(".h5");
-  std::string xdmfFilename = 
-    std::string(hdf5filename, 0, indexExt) + ".xmf";
-  metafile.write(xdmfFilename.c_str(), _hdf5Filename().c_str());
+  int rank = 0;
+  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
+  if (!rank) {
+    Xdmf metafile;
+    const std::string& hdf5filename = _hdf5Filename();
+    const int indexExt = hdf5filename.find(".h5");
+    std::string xdmfFilename = 
+      std::string(hdf5filename, 0, indexExt) + ".xmf";
+    metafile.write(xdmfFilename.c_str(), _hdf5Filename().c_str());
+  } // if
 } // close
 
 // ----------------------------------------------------------------------
@@ -312,7 +314,7 @@
     err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
     err = VecSetBlockSize(vector, blockSize); CHECK_PETSC_ERROR(err);
 
-    if (!rank && 0 == istep) {
+    if (0 == istep) {
       hid_t h5 = -1;
       err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
       assert(h5 >= 0);
@@ -366,7 +368,6 @@
     field.scatterSectionToVector(context);
     PetscVec vector = field.vector(context);
     assert(vector);
-
     const ALE::Obj<typename mesh_type::RealSection>& section = field.section();
     assert(!section.isNull());      
     assert(!sieveMesh->getLabelStratum(labelName, depth).isNull());
@@ -401,7 +402,7 @@
     err = PetscViewerHDF5PopGroup(_viewer); CHECK_PETSC_ERROR(err);
     err = VecSetBlockSize(vector, blockSize); CHECK_PETSC_ERROR(err);
 
-    if (!rank && 0 == istep) {
+    if (0 == istep) {
       hid_t h5 = -1;
       err = PetscViewerHDF5GetFileId(_viewer, &h5); CHECK_PETSC_ERROR(err);
       assert(h5 >= 0);
@@ -451,7 +452,7 @@
 { // _writeTimeStamp
   PetscErrorCode err = 0;
 
-  if (!rank) {
+  if (0 == rank) {
     err = VecSetValue(_tstamp, 0, t, INSERT_VALUES); CHECK_PETSC_ERROR(err);
   } // if
   err = VecAssemblyBegin(_tstamp); CHECK_PETSC_ERROR(err);

Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -261,11 +261,15 @@
   _tstampIndex = 0;
   deallocate();
 
-  Xdmf metafile;
-  const std::string& hdf5filename = _hdf5Filename();
-  const int indexExt = hdf5filename.find(".h5");
-  std::string xdmfFilename = std::string(hdf5filename, 0, indexExt) + ".xmf";
-  metafile.write(xdmfFilename.c_str(), _hdf5Filename().c_str());
+  int rank = 0;
+  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
+  if (!rank) {
+    Xdmf metafile;
+    const std::string& hdf5filename = _hdf5Filename();
+    const int indexExt = hdf5filename.find(".h5");
+    std::string xdmfFilename = std::string(hdf5filename, 0, indexExt) + ".xmf";
+    metafile.write(xdmfFilename.c_str(), _hdf5Filename().c_str());
+  } // if
 } // close
 
 // ----------------------------------------------------------------------

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Distributor.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -102,11 +102,11 @@
   // Setup and allocate field
   const int fiberDim = 1;
   topology::Field<topology::Mesh> partition(mesh);
+  partition.newSection(topology::FieldBase::CELLS_FIELD, fiberDim);
+  partition.allocate();
   partition.scale(1.0);
   partition.label("partition");
   partition.vectorFieldType(topology::FieldBase::SCALAR);
-  partition.newSection(topology::FieldBase::CELLS_FIELD, fiberDim);
-  partition.allocate();
   const ALE::Obj<RealSection>& partitionSection = partition.section();
   assert(!partitionSection.isNull());
 

Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/Field.cc	2011-07-22 20:35:08 UTC (rev 18797)
@@ -87,6 +87,7 @@
       err = VecDestroy(&s_iter->second.scatterVec);CHECK_PETSC_ERROR(err);
     } // if
   } // for
+  _scatters.clear();
 } // deallocate
 
 // ----------------------------------------------------------------------
@@ -146,6 +147,9 @@
 void
 pylith::topology::Field<mesh_type, section_type>::newSection(void)
 { // newSection
+  // Clear memory
+  clear();
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
 
@@ -167,8 +171,12 @@
 { // newSection
   typedef typename mesh_type::SieveMesh::point_type point_type;
 
+  // Clear memory
+  clear();
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
+
   if (fiberDim < 0) {
     std::ostringstream msg;
     msg << "Fiber dimension (" << fiberDim << ") for field '" << _metadata.label
@@ -203,6 +211,9 @@
 { // newSection
   typedef typename mesh_type::SieveMesh::point_type point_type;
 
+  // Clear memory
+  clear();
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
   if (fiberDim < 0) {
@@ -261,6 +272,9 @@
 pylith::topology::Field<mesh_type, section_type>::newSection(const Field& src,
 					       const int fiberDim)
 { // newSection
+  // Clear memory
+  clear();
+
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Field");
 
@@ -299,21 +313,22 @@
 void
 pylith::topology::Field<mesh_type, section_type>::cloneSection(const Field& src)
 { // cloneSection
+  std::string origLabel = _metadata.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");
 
-  deallocate();
-  std::string origLabel = _metadata.label;
   _metadata = src._metadata;
   label(origLabel.c_str());
 
-  const ALE::Obj<section_type>& srcSection = src.section();
-  if (!srcSection.isNull() && _section.isNull()) {
-    logger.stagePop();
-    newSection();
-    logger.stagePush("Field");
-  }
-
   if (!_section.isNull()) {
     if (!srcSection->sharedStorage()) {
       _section->setAtlas(srcSection->getAtlas());
@@ -1006,14 +1021,17 @@
   err = VecSetFromOptions(sinfo.vector); CHECK_PETSC_ERROR(err);  
 
 #if 0
-  std::cout << "CONTEXT: " << context 
+  std::cout << "["<<sieveMesh->commRank()<<"] CONTEXT: " << context 
 	    << ", orderLabel: " << orderLabel
 	    << ", section size w/BC: " << _section->sizeWithBC()
 	    << ", section size: " << _section->size()
 	    << ", section storage size: " << _section->getStorageSize()
 	    << ", global numbering size: " << numbering->getGlobalSize()
-	    << ", global size: " << order->getGlobalSize()
+	    << ", global order size: " << order->getGlobalSize()
+	    << ", local numbering size: " << numbering->getLocalSize()
+	    << ", local order size: " << order->getLocalSize()
 	    << ", scatter from size: " << sinfo.scatter->from_n
+	    << ", scatter: " << sinfo.scatter
 	    << std::endl;
 #endif
   
@@ -1139,8 +1157,33 @@
 { // _getScatter
   assert(context);
 
-  const bool isNewScatter = _scatters.find(context) == _scatters.end();
+  bool isNewScatter = _scatters.find(context) == _scatters.end();
 
+  // Synchronize creation of scatter (empty sections may have
+  // leftover, reusable scatters that need to be cleared out).
+  int numNewScatterLocal = (isNewScatter) ? 1 : 0;
+  int numNewScatter = 0;
+  MPI_Allreduce(&numNewScatterLocal, &numNewScatter, 1, MPI_INT, MPI_MAX,
+		_mesh.comm());
+  if (numNewScatter && !isNewScatter) {
+    // remove old scatter
+    ScatterInfo& sinfo = _scatters[context];
+    PetscErrorCode err = 0;
+    if (sinfo.vector) {
+      err = VecDestroy(&sinfo.vector);CHECK_PETSC_ERROR(err);
+    } // if
+    if (sinfo.scatter) {
+      err = VecScatterDestroy(&sinfo.scatter);CHECK_PETSC_ERROR(err);
+    } // if
+
+    if (sinfo.scatterVec) {
+      err = VecDestroy(&sinfo.scatterVec);CHECK_PETSC_ERROR(err);
+    } // if
+
+    _scatters.erase(context);
+    isNewScatter = true;
+  } // if
+
   if (isNewScatter && !createOk) {
     std::ostringstream msg;
     msg << "Scatter for context '" << context << "' does not exist.";

Modified: short/3D/PyLith/trunk/pylith/problems/ExplicitLumped.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/ExplicitLumped.py	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/pylith/problems/ExplicitLumped.py	2011-07-22 20:35:08 UTC (rev 18797)
@@ -132,10 +132,10 @@
     self._info.log("Creating lumped Jacobian matrix.")
     from pylith.topology.topology import MeshField
     jacobian = MeshField(self.mesh)
+    jacobian.newSection(jacobian.VERTICES_FIELD, dimension)
+    jacobian.allocate()
     jacobian.label("jacobian")
     jacobian.vectorFieldType(jacobian.VECTOR)
-    jacobian.newSection(jacobian.VERTICES_FIELD, dimension)
-    jacobian.allocate()
     self.jacobian = jacobian
     self._debug.log(resourceUsageString())
 

Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py	2011-07-22 20:20:23 UTC (rev 18796)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py	2011-07-22 20:35:08 UTC (rev 18797)
@@ -473,9 +473,9 @@
 
     lengthScale = normalizer.lengthScale()
     solution = self.fields.get("dispIncr(t->t+dt)")
+    solution.newSection(solution.VERTICES_FIELD, dimension)
     solution.vectorFieldType(solution.VECTOR)
     solution.scale(lengthScale.value)
-    solution.newSection(solution.VERTICES_FIELD, dimension)
     if self.splitFields:
       solution.splitDefault()
       for integrator in self.integratorsMesh + self.integratorsSubMesh:



More information about the CIG-COMMITS mailing list