[cig-commits] r14823 - in short/3D/PyLith/branches/pylith-swig: libsrc/faults libsrc/topology pylith/faults pylith/meshio

brad at geodynamics.org brad at geodynamics.org
Wed Apr 29 21:52:58 PDT 2009


Author: brad
Date: 2009-04-29 21:52:57 -0700 (Wed, 29 Apr 2009)
New Revision: 14823

Modified:
   short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh
   short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh
   short/3D/PyLith/branches/pylith-swig/pylith/faults/Fault.py
   short/3D/PyLith/branches/pylith-swig/pylith/meshio/OutputFaultKin.py
   short/3D/PyLith/branches/pylith-swig/pylith/meshio/OutputManagerSubMesh.py
Log:
Worked on fixing output for fault stuff.

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc	2009-04-30 00:46:29 UTC (rev 14822)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.cc	2009-04-30 04:52:57 UTC (rev 14823)
@@ -653,11 +653,12 @@
 				  const char* name,
 				  const topology::SolutionFields* fields)
 { // vertexField
-#if 0
-  assert(0 != fields);
+  assert(0 != _faultMesh);
+  assert(0 != _quadrature);
   assert(0 != _normalizer);
+  assert(0 != _fields);
 
-  const int cohesiveDim = _faultMesh.dimension();
+  const int cohesiveDim = _faultMesh->dimension();
   const int spaceDim = _quadrature->spaceDim();
 
   const int slipStrLen = strlen("final_slip");
@@ -667,68 +668,64 @@
   int fiberDim = 0;
   if (0 == strcasecmp("slip", name)) {
     const topology::Field<topology::SubMesh>& cumSlip = 
-      fields->get("cumulative slip");
+      _fields->get("cumulative slip");
     return cumSlip;
 
   } else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
     const ALE::Obj<RealSection>& orientationSection =
-      fields->get("orientation").section();
+      _fields->get("orientation").section();
     assert(!orientationSection.isNull());
-    strikeSection = orientationSection->getFibration(0);
+    const ALE::Obj<RealSection>& dirSection =
+      orientationSection->getFibration(0);
+    assert(!dirSection.isNull());
     _allocateBufferVectorField();
     assert(0 != _bufferVectorField);
-    const ALE::Obj<RealSection>& bufferSection =
-      _bufferVectorField->section();
-    bufferSection.values(strikeSection.values());
-    scale = 0.0;
-    fiberDim = spaceDim;
+    _bufferVectorField->copy(dirSection);
+    return *_bufferVectorField;
 
   } else if (2 == cohesiveDim && 0 == strcasecmp("dip_dir", name)) {
-    *fieldType = VECTOR_FIELD;
-    _bufferTmp = _orientation->getFibration(1);
-    scale = 0.0;
-    fiberDim = spaceDim;
+    const ALE::Obj<RealSection>& orientationSection =
+      _fields->get("orientation").section();
+    assert(!orientationSection.isNull());
+    const ALE::Obj<RealSection>& dirSection =
+      orientationSection->getFibration(1);
+    _allocateBufferVectorField();
+    assert(0 != _bufferVectorField);
+    _bufferVectorField->copy(dirSection);
+    return *_bufferVectorField;
 
   } else if (0 == strcasecmp("normal_dir", name)) {
-    *fieldType = VECTOR_FIELD;
+    const ALE::Obj<RealSection>& orientationSection =
+      _fields->get("orientation").section();
+    assert(!orientationSection.isNull());
     const int space = 
       (0 == cohesiveDim) ? 0 : (1 == cohesiveDim) ? 1 : 2;
-    _bufferTmp = _orientation->getFibration(space);
-    scale = 0.0;
-    fiberDim = spaceDim;
+    const ALE::Obj<RealSection>& dirSection =
+      orientationSection->getFibration(space);
+    assert(!dirSection.isNull());
+    _allocateBufferVectorField();
+    assert(0 != _bufferVectorField);
+    _bufferVectorField->copy(dirSection);
+    return *_bufferVectorField;
 
   } else if (0 == strncasecmp("final_slip_X", name, slipStrLen)) {
     const std::string value = std::string(name).substr(slipStrLen+1);
 
-    *fieldType = VECTOR_FIELD;
     const srcs_type::const_iterator s_iter = _eqSrcs.find(value);
     assert(s_iter != _eqSrcs.end());
-    _allocateBufferVertexVector();
-    topology::FieldOps::copyValues(_bufferVertexVector, 
-				   s_iter->second->finalSlip());
-    _bufferTmp = _bufferVertexVector;
-    scale = _normalizer->lengthScale();
-    fiberDim = spaceDim;
+    return s_iter->second->finalSlip();
 
   } else if (0 == strncasecmp("slip_time_X", name, timeStrLen)) {
-    *fieldType = SCALAR_FIELD;
     const std::string value = std::string(name).substr(timeStrLen+1);
     const srcs_type::const_iterator s_iter = _eqSrcs.find(value);
     assert(s_iter != _eqSrcs.end());
-    _allocateBufferVertexScalar();
-    topology::FieldOps::copyValues(_bufferVertexScalar, 
-				   s_iter->second->slipTime());
-    _bufferTmp = _bufferVertexScalar;
-    scale = _normalizer->timeScale();
-    fiberDim = 1;
+    return s_iter->second->slipTime();
 
   } else if (0 == strcasecmp("traction_change", name)) {
-    *fieldType = VECTOR_FIELD;
-    const ALE::Obj<RealSection>& solution = fields->getSolution();
-    _calcTractionsChange(&_bufferVertexVector, mesh, solution);
-    _bufferTmp = _bufferVertexVector;
-    scale = _normalizer->pressureScale();
-    fiberDim = spaceDim;
+    assert(0 != fields);
+    const topology::Field<topology::Mesh>& solution = fields->solution();
+    _allocateBufferVectorField();
+    _calcTractionsChange(_bufferVectorField, solution);
     
   } else {
     std::ostringstream msg;
@@ -737,24 +734,8 @@
     throw std::runtime_error(msg.str());
   } // else
 
-  if (0 != scale) {
-    // dimensionalize values
-    double_array values(fiberDim);
-    const ALE::Obj<Mesh::label_sequence>& vertices = _faultMesh->depthStratum(0);
-    assert(!vertices.isNull());
-    const Mesh::label_sequence::iterator verticesEnd = vertices->end();
-    for (Mesh::label_sequence::iterator v_iter=vertices->begin(); 
-	 v_iter != verticesEnd;
-	 ++v_iter) {
-      assert(fiberDim == _bufferTmp->getFiberDimension(*v_iter));
-      _bufferTmp->restrictPoint(*v_iter, &values[0], values.size());
-      _normalizer->dimensionalize(&values[0], values.size(), scale);
-      _bufferTmp->updatePointAll(*v_iter, &values[0]);
-    } // for
-  } // if
-
-  return _bufferTmp;
-#endif
+  assert(0 != _bufferScalarField);
+  return *_bufferScalarField;
 } // vertexField
 
 // ----------------------------------------------------------------------
@@ -771,7 +752,8 @@
   throw std::runtime_error(msg.str());
 
   // Return generic section to satisfy member function definition.
-  return *_bufferVectorField;
+  assert(0 != _bufferScalarField);
+  return *_bufferScalarField;
 } // cellField
 
 // ----------------------------------------------------------------------
@@ -1218,5 +1200,42 @@
 #endif
 } // _calcTractionsChange
 
+// ----------------------------------------------------------------------
+// Allocate buffer for vector field.
+void
+pylith::faults::FaultCohesiveKin::_allocateBufferVectorField(void)
+{ // _allocateBufferVectorField
+  if (0 != _bufferVectorField)
+    return;
 
+  // Create vector field; use same shape/chart as cumulative slip field.
+  assert(0 != _faultMesh);
+  assert(0 != _fields);
+  _bufferVectorField = new topology::Field<topology::SubMesh>(*_faultMesh);
+  const topology::Field<topology::SubMesh>& slip = 
+    _fields->get("cumulative slip");
+  _bufferVectorField->newSection(slip);
+  _bufferVectorField->allocate();
+  _bufferVectorField->zero();
+} // _allocateBufferVectorField
+
+// ----------------------------------------------------------------------
+// Allocate buffer for scalar field.
+void
+pylith::faults::FaultCohesiveKin::_allocateBufferScalarField(void)
+{ // _allocateBufferScalarField
+  if (0 != _bufferScalarField)
+    return;
+
+  // Create vector field; use same shape/chart as area field.
+  assert(0 != _faultMesh);
+  assert(0 != _fields);
+  _bufferScalarField = new topology::Field<topology::SubMesh>(*_faultMesh);
+  const topology::Field<topology::SubMesh>& area = _fields->get("area");
+  _bufferScalarField->newSection(area);
+  _bufferScalarField->allocate();
+  _bufferScalarField->zero();
+} // _allocateBufferScalarField
+
+
 // End of file 

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh	2009-04-30 00:46:29 UTC (rev 14822)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/faults/FaultCohesiveKin.hh	2009-04-30 04:52:57 UTC (rev 14823)
@@ -209,6 +209,12 @@
   void _calcTractionsChange(topology::Field<topology::SubMesh>* tractions,
 			    const topology::Field<topology::Mesh>& solution);
 
+  /// Allocate buffer for vector field.
+  void _allocateBufferVectorField(void);
+
+  /// Allocate buffer for scalar field.
+  void _allocateBufferScalarField(void);
+
   // NOT IMPLEMENTED ////////////////////////////////////////////////////
 private :
 

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc	2009-04-30 00:46:29 UTC (rev 14822)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.cc	2009-04-30 04:52:57 UTC (rev 14823)
@@ -286,6 +286,47 @@
 } // copy
 
 // ----------------------------------------------------------------------
+// Copy field values.
+template<typename mesh_type>
+void
+pylith::topology::Field<mesh_type>::copy(const ALE::Obj<typename mesh_type::RealSection>& osection)
+{ // copy
+  // Check compatibility of sections
+  const int srcSize = (!osection.isNull()) ? osection->size() : 0;
+  const int dstSize = (!_section.isNull()) ? _section->size() : 0;
+  if (srcSize != dstSize) {
+    std::ostringstream msg;
+
+    msg << "Cannot copy values from Sieve section "
+	<< _label << "'. Sections are incompatible.\n"
+	<< "  Source section:\n"
+	<< "    size: " << srcSize
+	<< "  Destination section:\n"
+	<< "    space dim: " << spaceDim() << "\n"
+	<< "    vector field type: " << _vecFieldType << "\n"
+	<< "    scale: " << _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));
+      _section->updatePoint(*c_iter, osection->restrictPoint(*c_iter));
+    } // for
+  } // if
+} // copy
+
+// ----------------------------------------------------------------------
 // Add two fields, storing the result in one of the fields.
 template<typename mesh_type>
 void

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh	2009-04-30 00:46:29 UTC (rev 14822)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Field.hh	2009-04-30 04:52:57 UTC (rev 14823)
@@ -189,6 +189,12 @@
    */
   void copy(const Field& field);
 
+  /** Copy field values.
+   *
+   * @param field Field to copy.
+   */
+  void copy(const ALE::Obj<typename mesh_type::RealSection>& field);
+
   /** Add two fields, storing the result in one of the fields.
    *
    * @param field Field to add.

Modified: short/3D/PyLith/branches/pylith-swig/pylith/faults/Fault.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/faults/Fault.py	2009-04-30 00:46:29 UTC (rev 14822)
+++ short/3D/PyLith/branches/pylith-swig/pylith/faults/Fault.py	2009-04-30 04:52:57 UTC (rev 14823)
@@ -161,10 +161,10 @@
     ModuleFault.initialize(self, 
                            self.mesh, self.upDir, self.normalDir, self.matDB)
 
-    #if None != self.output:
-    #  self.output.initialize(normalizer, self.faultQuadrature)
-    #  self.output.writeInfo()
-    #  self.output.open(totalTime, numTimeSteps)
+    if None != self.output:
+      self.output.initialize(normalizer, self.faultQuadrature)
+      self.output.writeInfo()
+      self.output.open(totalTime, numTimeSteps)
 
     self._logger.eventEnd(logEvent)
     return

Modified: short/3D/PyLith/branches/pylith-swig/pylith/meshio/OutputFaultKin.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/meshio/OutputFaultKin.py	2009-04-30 00:46:29 UTC (rev 14822)
+++ short/3D/PyLith/branches/pylith-swig/pylith/meshio/OutputFaultKin.py	2009-04-30 04:52:57 UTC (rev 14823)
@@ -25,39 +25,37 @@
   Python object for managing output of finite-element information for
   faults with kinematic ruptures.
 
+  Inventory
+
+  @class Inventory
+  Python object for managing OutputFaultKin facilities and properties.
+  
+  \b Properties
+  @li \b vertex_info_fields Names of vertex info fields to output.
+  @li \b vertex_data_fields Names of vertex data fields to output.
+  
+  \b Facilities
+  @li None
+
   Factory: output_manager
   """
 
   # INVENTORY //////////////////////////////////////////////////////////
 
-  class Inventory(OutputManagerSubMesh.Inventory):
-    """
-    Python object for managing OutputFaultKin facilities and properties.
-    """
+  import pyre.inventory
 
-    ## @class Inventory
-    ## Python object for managing OutputFaultKin facilities and properties.
-    ##
-    ## \b Properties
-    ## @li \b vertex_info_fields Names of vertex info fields to output.
-    ## @li \b vertex_data_fields Names of vertex data fields to output.
-    ##
-    ## \b Facilities
-    ## @li None
+  vertexInfoFields = pyre.inventory.list("vertex_info_fields",
+                                         default=["normal_dir",
+                                                  "final_slip_rupture",
+                                                  "slip_time_rupture"])
+  vertexInfoFields.meta['tip'] = "Names of vertex info fields to output."
 
-    import pyre.inventory
+  vertexDataFields = pyre.inventory.list("vertex_data_fields", 
+                                         default=["slip",
+                                                  "traction_change"])
+  vertexDataFields.meta['tip'] = "Names of vertex data fields to output."
 
-    vertexInfoFields = pyre.inventory.list("vertex_info_fields",
-                                           default=["normal_dir",
-                                                    "final_slip_rupture",
-                                                    "slip_time_rupture"])
-    vertexInfoFields.meta['tip'] = "Names of vertex info fields to output."
 
-    vertexDataFields = pyre.inventory.list("vertex_data_fields", 
-                                           default=["slip",
-                                                    "traction_change"])
-    vertexDataFields.meta['tip'] = "Names of vertex data fields to output."
-
   # PUBLIC METHODS /////////////////////////////////////////////////////
 
   def __init__(self, name="outputfaultkin"):

Modified: short/3D/PyLith/branches/pylith-swig/pylith/meshio/OutputManagerSubMesh.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/meshio/OutputManagerSubMesh.py	2009-04-30 00:46:29 UTC (rev 14822)
+++ short/3D/PyLith/branches/pylith-swig/pylith/meshio/OutputManagerSubMesh.py	2009-04-30 04:52:57 UTC (rev 14823)
@@ -107,6 +107,7 @@
     ModuleOutputManager.appendVertexField(self, t, field)
     return
 
+
   def _appendCellField(self, t, field):
     """
     Call C++ appendCellField();



More information about the CIG-COMMITS mailing list