[cig-commits] r14201 - long/3D/SNAC/trunk/Snac/snac2vtk

[echoi at geodynamics.org]

Author: echoi
Date: 2009-03-03 07:46:40 -0800 (Tue, 03 Mar 2009)
New Revision: 14201

Modified:
   long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c
Log:
* Fixed a bug that can potentially lead to nan or inf failurePotential.



Modified: long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c
===================================================================
--- long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c	2009-03-03 03:06:53 UTC (rev 14200)
+++ long/3D/SNAC/trunk/Snac/snac2vtk/snac2vtk.c	2009-03-03 15:46:40 UTC (rev 14201)
@@ -1046,10 +1046,20 @@
      *  Calculate the failure potential for hillslope angle
      */
 /*     elementStressMeasures->failurePotential= ( (fabs(elementStressMeasures->maxShearStress)-1e6) */
-/* 					       /(-2*elementStressMeasures->pressure/3.0) ); */
-    elementStressMeasures->failurePotential= fabs(-elementStressMeasures->slopeShearStress/elementStressMeasures->slopeNormalStress);
-					       
+/* 					       /(-2*elementStressMeasures->pessure/3.0) ); */
 
+	/* If slopeNormalStress is 0, failurePotential is not defined. Assign some indicative value: -1 for now. */
+	/* If the computed slopeNormalStress is infinite, assign -1 again. */
+	/* In either case, a warning message would be desirable. Or a switch to turn off failure potential calculations might be better. */
+	/* -EChoi 2009/03/03 */
+	if( elementStressMeasures->slopeNormalStress == 0.0 )
+		elementStressMeasures->failurePotential = -1.0; 
+	else {
+		elementStressMeasures->failurePotential= fabs(-elementStressMeasures->slopeShearStress/elementStressMeasures->slopeNormalStress);
+		if( isinf( elementStressMeasures->failurePotential ) )
+			elementStressMeasures->failurePotential = -1.0;
+	}
+
 }