[cig-commits] r20615 - short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble

[brad at geodynamics.org]

Author: brad
Date: 2012-08-21 10:00:25 -0700 (Tue, 21 Aug 2012)
New Revision: 20615

Modified:
   short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble/ElasticityExplicit.cc
   short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
Log:
Adjust assert so nondimensional artificial viscosity for numerical damping can be zero.

Modified: short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble/ElasticityExplicit.cc	2012-08-21 16:30:28 UTC (rev 20614)
+++ short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble/ElasticityExplicit.cc	2012-08-21 17:00:25 UTC (rev 20615)
@@ -223,7 +223,7 @@
 				    _normalizer->densityScale());
 
   const PylithScalar dt = _dt;
-  assert(_normViscosity > 0.0);
+  assert(_normViscosity >= 0.0);
   assert(dt > 0);
   const PylithScalar viscosity = dt*_normViscosity;
 
@@ -495,7 +495,7 @@
             _normalizer->densityScale());
 
   const PylithScalar dt = _dt;
-  assert(_normViscosity > 0.0);
+  assert(_normViscosity >= 0.0);
   assert(dt > 0);
   const PylithScalar viscosity = dt*_normViscosity;
 

Modified: short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc	2012-08-21 16:30:28 UTC (rev 20614)
+++ short/3D/PyLith/branches/v1.7-stable/libsrc/pylith/feassemble/ElasticityExplicitLgDeform.cc	2012-08-21 17:00:25 UTC (rev 20615)
@@ -185,17 +185,14 @@
   assert(!accSection.isNull());
   RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
 
-#if 0 // Numerical damping not yet implemented
   scalar_array velCell(numBasis*spaceDim);
   const ALE::Obj<RealSection>& velSection = 
     fields->get("velocity(t)").section();
   assert(!velSection.isNull());
   RestrictVisitor velVisitor(*velSection, velCell.size(), &velCell[0]);
 
-  scalar_array dispAdjCell(numBasis*spaceDim);
-#endif
-
   scalar_array dispCell(numBasis*spaceDim);
+  scalar_array dispAdjCell(numBasis*spaceDim);
   const ALE::Obj<RealSection>& dispSection = fields->get("disp(t)").section();
   assert(!dispSection.isNull());
   RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
@@ -217,7 +214,7 @@
 				    _normalizer->densityScale());
 
   const PylithScalar dt = _dt;
-  assert(_normViscosity > 0.0);
+  assert(_normViscosity >= 0.0);
   assert(dt > 0);
   const PylithScalar viscosity = dt*_normViscosity;
 
@@ -311,21 +308,17 @@
     } // for
     PetscLogFlops(numQuadPts*(2+numBasis*(1+numBasis*(2*spaceDim))));
 
-#if 0 // Numerical damping not yet implemented. Is small strain
-      // formulation compatible with numerical damping?
-
     // Numerical damping. Compute displacements adjusted by velocity
     // times normalized viscosity.
     dispAdjCell = dispCell + viscosity * velCell;
-#endif
 
     // Compute B(transpose) * sigma, first computing strains
-    _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCell,
+    _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispAdjCell,
 		     numBasis, numQuadPts, spaceDim);
     calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
     const scalar_array& stressCell = _material->calcStress(strainCell, true);
 
-    CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispCell);
+    CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispAdjCell);
 
     // Assemble cell contribution into field
     residualVisitor.clear();
@@ -426,17 +419,14 @@
   assert(!accSection.isNull());
   RestrictVisitor accVisitor(*accSection, accCell.size(), &accCell[0]);
 
-#if 0 // Numerical damping not yet implemented
   scalar_array velCell(numBasis*spaceDim);
   const ALE::Obj<RealSection>& velSection = 
     fields->get("velocity(t)").section();
   assert(!velSection.isNull());
   RestrictVisitor velVisitor(*velSection, velCell.size(), &velCell[0]);
 
-  scalar_array dispAdjCell(numBasis*spaceDim);
-#endif
-
   scalar_array dispCell(numBasis*spaceDim);
+  scalar_array dispAdjCell(numBasis*spaceDim);
   const ALE::Obj<RealSection>& dispSection = fields->get("disp(t)").section();
   assert(!dispSection.isNull());
   RestrictVisitor dispVisitor(*dispSection, dispCell.size(), &dispCell[0]);
@@ -457,6 +447,11 @@
     _normalizer->pressureScale() / (_normalizer->lengthScale() *
 				    _normalizer->densityScale());
 
+  const PylithScalar dt = _dt;
+  assert(_normViscosity >= 0.0);
+  assert(dt > 0);
+  const PylithScalar viscosity = dt*_normViscosity;
+
   // Get parameters used in integration.
   scalar_array valuesIJ(numBasis);
 
@@ -552,20 +547,17 @@
 	  accCell[iBasis*spaceDim+iDim];
     PetscLogFlops(numQuadPts*(4+numBasis*3));
 
-#if 0 // Numerical damping not yet implemented. Is small strain
-      // formulation compatible with numerical damping?
-
     // Numerical damping. Compute displacements adjusted by velocity
     // times normalized viscosity.
     dispAdjCell = dispCell + viscosity * velCell;
-#endif
+
     // Compute B(transpose) * sigma, first computing strains
-    _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispCell,
+    _calcDeformation(&deformCell, basisDeriv, coordinatesCell, dispAdjCell,
 		     numBasis, numQuadPts, spaceDim);
     calcTotalStrainFn(&strainCell, deformCell, numQuadPts);
     const scalar_array& stressCell = _material->calcStress(strainCell, true);
 
-    CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispCell);
+    CALL_MEMBER_FN(*this, elasticityResidualFn)(stressCell, dispAdjCell);
     
     // Assemble cell contribution into field
     residualVisitor.clear();