[cig-commits] r16268 - short/3D/PyLith/trunk/libsrc/meshio

[brad at geodynamics.org]

Author: brad
Date: 2010-02-16 15:46:48 -0800 (Tue, 16 Feb 2010)
New Revision: 16268

Modified:
   short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc
Log:
Fixed bug in getting cells for cell filter (mimic depth/label in DataWriterVTK to account for fault and boundary meshes.

Modified: short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc	2010-02-16 20:17:44 UTC (rev 16267)
+++ short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc	2010-02-16 23:46:48 UTC (rev 16268)
@@ -93,10 +93,18 @@
   
   const ALE::Obj<SieveMesh>& sieveMesh = fieldIn.mesh().sieveMesh();
   assert(!sieveMesh.isNull());
+  const int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
+  const int depth = (0 == label) ? cellDepth : labelId;
+  const std::string labelName = (0 == label) ?
+    ((sieveMesh->hasLabel("censored depth")) ?
+     "censored depth" : "depth") : label;
 
-  const ALE::Obj<label_sequence>& cells = (0 == label) ?
-    sieveMesh->heightStratum(0) :
-    sieveMesh->getLabelStratum(label, labelId);
+  std::cout << "labelName: " << labelName
+	    << ", depth: " << depth
+	    << std::endl;
+
+  const ALE::Obj<label_sequence>& cells = 
+    sieveMesh->getLabelStratum(labelName, depth);
   assert(!cells.isNull());
   const typename label_sequence::iterator cellsBegin = cells->begin();
   const typename label_sequence::iterator cellsEnd = cells->end();
@@ -109,6 +117,11 @@
   const int fiberDim = totalFiberDim / numQuadPts;
   assert(fiberDim * numQuadPts == totalFiberDim);
 
+  fieldIn.view("INPUT FIELD");
+  std::cout << "CELLS SIZE: " << cells->size() << std::endl;
+  std::cout << "FIBERDIM: " << fiberDim << std::endl;
+  std::cout << "CELLS BEGIN: " << *cellsBegin << std::endl;
+
   // Allocate field if necessary
   ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
   logger.stagePush("Output");
@@ -168,6 +181,7 @@
        c_iter != cellsEnd;
        ++c_iter) {
     const double* values = sectionIn->restrictPoint(*c_iter);
+    assert(totalFiberDim == sectionIn->getFiberDimension(*c_iter));
     
     fieldAvgCell = 0.0;
     for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
@@ -178,6 +192,8 @@
   } // for
   PetscLogFlops( cells->size() * numQuadPts*fiberDim*3 );
 
+  _fieldAvg->view("AVERAGED FIELD");
+
   return *_fieldAvg;
 } // filter