[cig-commits] r15039 - in short/3D/PyLith/trunk: libsrc/faults libsrc/topology pylith/problems

[brad at geodynamics.org]

Author: brad
Date: 2009-05-24 16:44:29 -0700 (Sun, 24 May 2009)
New Revision: 15039

Modified:
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
   short/3D/PyLith/trunk/libsrc/topology/SubMesh.cc
   short/3D/PyLith/trunk/pylith/problems/Formulation.py
   short/3D/PyLith/trunk/pylith/problems/Implicit.py
Log:
Fixed some nondimensionalization bugs.

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2009-05-24 20:23:46 UTC (rev 15038)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2009-05-24 23:44:29 UTC (rev 15039)
@@ -152,12 +152,6 @@
 
   // Compute tributary area for each vertex in fault mesh.
   _calcArea();
-
-  // Create empty tractions field for change in fault tractions.
-  _fields->add("tractions", "traction_change");
-  topology::Field<topology::SubMesh>& tractions = _fields->get("tractions");
-  tractions.vectorFieldType(topology::FieldBase::VECTOR);
-  tractions.scale(_normalizer->pressureScale());
 } // initialize
 
 // ----------------------------------------------------------------------
@@ -1165,8 +1159,8 @@
   assert(0 != _fields);
   assert(0 != _normalizer);
 
-  tractions->scale(_normalizer->pressureScale());
   tractions->label("traction_change");
+  tractions->scale(_normalizer->pressureScale());
 
   // Get vertices from mesh of domain.
   const ALE::Obj<SieveMesh>& sieveMesh = dispT.mesh().sieveMesh();
@@ -1237,6 +1231,7 @@
   assert(!tractionsSection.isNull());
   tractions->zero();
   
+  const double pressureScale = tractions->scale();
   for (SieveMesh::label_sequence::iterator v_iter = verticesBegin; 
        v_iter != verticesEnd;
        ++v_iter)
@@ -1255,7 +1250,11 @@
       const double* areaVertex = areaSection->restrictPoint(vertexFault);
       assert(0 != areaVertex);
 
-      const double scale = stiffnessVertex[0] / areaVertex[0];
+      // Account for conditioning and then "nondimensionalize" so that
+      // values appear to be nondimensionalized for future
+      // compatibility with nondimensionalization.
+      const double scale = 
+	stiffnessVertex[0] / (areaVertex[0] * pressureScale);
       for (int i=0; i < fiberDim; ++i)
 	tractionsVertex[i] = dispTVertex[i] * scale;
 

Modified: short/3D/PyLith/trunk/libsrc/topology/SubMesh.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/SubMesh.cc	2009-05-24 20:23:46 UTC (rev 15038)
+++ short/3D/PyLith/trunk/libsrc/topology/SubMesh.cc	2009-05-24 23:44:29 UTC (rev 15039)
@@ -71,6 +71,9 @@
   } // if
   _mesh->setRealSection("coordinates", 
 			meshSieveMesh->getRealSection("coordinates"));
+  if (meshSieveMesh->hasRealSection("coordinates_dimensioned"))
+  _mesh->setRealSection("coordinates_dimensioned", 
+			meshSieveMesh->getRealSection("coordinates_dimensioned"));
 
   // Create the parallel overlap
   const ALE::Obj<SieveMesh::sieve_type>& sieve = _mesh->getSieve();

Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py	2009-05-24 20:23:46 UTC (rev 15038)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py	2009-05-24 23:44:29 UTC (rev 15039)
@@ -213,6 +213,7 @@
     self.fields.solutionName(solnName)
     solution = self.fields.solution()
     solution.vectorFieldType(solution.VECTOR)
+    solution.scale(normalizer.lengthScale().value)
     solution.newSection(solution.VERTICES_FIELD, dimension)
     for constraint in self.constraints:
       constraint.setConstraintSizes(solution)

Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py	2009-05-24 20:23:46 UTC (rev 15038)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py	2009-05-24 23:44:29 UTC (rev 15039)
@@ -97,9 +97,12 @@
     self.fields.solutionName("dispIncr(t->t+dt)")
 
     # Set fields to zero
+    lengthScale = normalizer.lengthScale()
     disp = self.fields.get("disp(t)")
+    disp.scale(lengthScale.value)
     disp.zero()
     dispIncr = self.fields.get("dispIncr(t->t+dt)")
+    dispIncr.scale(lengthScale.value)
     dispIncr.zero()
     residual = self.fields.get("residual")
     residual.zero()