[cig-commits] r13447 - in short/3D/PyLith/trunk: libsrc/feassemble modulesrc/bc

[brad at geodynamics.org]

Author: brad
Date: 2008-12-03 11:22:41 -0800 (Wed, 03 Dec 2008)
New Revision: 13447

Modified:
   short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
   short/3D/PyLith/trunk/modulesrc/bc/bc.pyxe.src
Log:
Fixed nondimensionalization bugs.

Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc	2008-12-03 19:22:16 UTC (rev 13446)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc	2008-12-03 19:22:41 UTC (rev 13447)
@@ -19,6 +19,9 @@
 
 #include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
 #include "pylith/topology/FieldsManager.hh" // USES FieldsManager
+
+#include "spatialdata/units/Nondimensional.hh" // USES Nondimensional
+
 #include "pylith/utils/array.hh" // USES double_array
 
 #include <cstring> // USES memcpy()
@@ -309,6 +312,7 @@
   const ALE::Obj<real_section_type>& coordinates = 
     mesh->getRealSection("coordinates");
   assert(!coordinates.isNull());
+  const ALE::Obj<real_section_type>& disp = fields->getSolution();
   
   // Get cell geometry information that doesn't depend on cell
   const int numQuadPts = _quadrature->numQuadPts();
@@ -322,6 +326,10 @@
   const int totalFiberDim = numQuadPts * tensorSize;
   double_array totalStrain(totalFiberDim);
   totalStrain = 0.0;
+  double_array stress(totalFiberDim);
+
+  assert(0 != _normalizer);
+  const double pressureScale = _normalizer->pressureScale();
   
   // Allocate buffer for property field.
   if (field->isNull() || 
@@ -333,8 +341,6 @@
     mesh->allocate(*field);
   } // if
   
-  const ALE::Obj<real_section_type>& disp = fields->getSolution();
-  /// const int dispAtlasTag = fields->getSolutionAtlasTag(materialId);
   
   // Loop over cells
   int c_index = 0;
@@ -358,7 +364,9 @@
       (*field)->updatePoint(*c_iter, &totalStrain[0]);
     } else {
       _material->getPropertiesCell(*c_iter, numQuadPts);
-      const double_array& stress = _material->calcStress(totalStrain);
+      stress = _material->calcStress(totalStrain);
+      _normalizer->dimensionalize(&stress[0], stress.size(),
+				  pressureScale);
       (*field)->updatePoint(*c_iter, &stress[0]);	
     } // else
   } // for
@@ -395,9 +403,14 @@
   const int numQuadPts = _quadrature->numQuadPts();
   
   // Allocate vector for total strain
-  double_array totalStrain(numQuadPts*tensorSize);
+  const int totalFiberDim = numQuadPts*tensorSize;
+  double_array totalStrain(totalFiberDim);
   totalStrain = 0.0;
+  double_array stress(totalFiberDim);
   
+  assert(0 != _normalizer);
+  const double pressureScale = _normalizer->pressureScale();
+  
   // Allocate buffer for tensor field.
   if (field->isNull()) {
     const int fiberDim = numQuadPts * tensorSize;
@@ -414,7 +427,9 @@
        ++c_iter) {
     (*field)->restrictPoint(*c_iter, &totalStrain[0], totalStrain.size());
     _material->getPropertiesCell(*c_iter, numQuadPts);
-    const double_array& stress = _material->calcStress(totalStrain);
+    stress = _material->calcStress(totalStrain);
+    _normalizer->dimensionalize(&stress[0], stress.size(),
+				pressureScale);
     (*field)->updatePoint(*c_iter, &stress[0]);	
   } // for
 } // _calcStressFromStrain

Modified: short/3D/PyLith/trunk/modulesrc/bc/bc.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/bc.pyxe.src	2008-12-03 19:22:16 UTC (rev 13446)
+++ short/3D/PyLith/trunk/modulesrc/bc/bc.pyxe.src	2008-12-03 19:22:41 UTC (rev 13447)
@@ -12,8 +12,6 @@
 
 #header{
 #include "pylith/bc/BoundaryCondition.hh"
-#include "pylith/feassemble/Constraint.hh"
-#include "pylith/feassemble/Integrator.hh"
 #include "pylith/bc/AbsorbingDampers.hh"
 #include "pylith/bc/DirichletPoints.hh"
 #include "pylith/bc/DirichletBoundary.hh"
@@ -332,6 +330,37 @@
       free(fixedDOF)
 
 
+  property normalizer:
+    def __set__(self, value):
+      """
+      Set nondimensionalizer.
+      """
+      # create shim for method 'normalizer'
+      #embed{ void DirichletPoints_normalizer_set(void* objVptr, void* dimVptr)
+      try {
+        assert(0 != objVptr);
+        assert(0 != dimVptr);
+        spatialdata::units::Nondimensional* dim =
+          (spatialdata::units::Nondimensional*) dimVptr;
+        ((pylith::bc::DirichletPoints*) objVptr)->normalizer(*dim);
+      } catch (const std::exception& err) {
+        PyErr_SetString(PyExc_RuntimeError,
+                        const_cast<char*>(err.what()));
+      } catch (const ALE::Exception& err) {
+        PyErr_SetString(PyExc_RuntimeError,
+                        const_cast<char*>(err.msg().c_str()));
+      } catch (...) {
+        PyErr_SetString(PyExc_RuntimeError,
+                        "Caught unknown C++ exception.");
+      } // try/catch
+      #}embed
+      if not value.name == "spatialdata_units_Nondimensional":
+        raise TypeError, \
+              "Argument must be extension module type " \
+              "'spatialdata::units::Nondimensional'."
+      DirichletPoints_normalizer_set(self.thisptr, ptrFromHandle(value))
+
+
   def setConstraintSizes(self, field, mesh):
     """
     Set number of degrees of freedom that are constrained at points in field.
@@ -577,6 +606,37 @@
       free(fixedDOF)
 
 
+  property normalizer:
+    def __set__(self, value):
+      """
+      Set nondimensionalizer.
+      """
+      # create shim for method 'normalizer'
+      #embed{ void DirichletBoundary_normalizer_set(void* objVptr, void* dimVptr)
+      try {
+        assert(0 != objVptr);
+        assert(0 != dimVptr);
+        spatialdata::units::Nondimensional* dim =
+          (spatialdata::units::Nondimensional*) dimVptr;
+        ((pylith::bc::DirichletBoundary*) objVptr)->normalizer(*dim);
+      } catch (const std::exception& err) {
+        PyErr_SetString(PyExc_RuntimeError,
+                        const_cast<char*>(err.what()));
+      } catch (const ALE::Exception& err) {
+        PyErr_SetString(PyExc_RuntimeError,
+                        const_cast<char*>(err.msg().c_str()));
+      } catch (...) {
+        PyErr_SetString(PyExc_RuntimeError,
+                        "Caught unknown C++ exception.");
+      } // try/catch
+      #}embed
+      if not value.name == "spatialdata_units_Nondimensional":
+        raise TypeError, \
+              "Argument must be extension module type " \
+              "'spatialdata::units::Nondimensional'."
+      DirichletBoundary_normalizer_set(self.thisptr, ptrFromHandle(value))
+
+
   def setConstraintSizes(self, field, mesh):
     """
     Set number of degrees of freedom that are constrained at points in field.