[cig-commits] r16495 - in long/3D/SNAC/trunk/Snac/plugins: plasticSPR viscoplasticSPR

[echoi at geodynamics.org]

Author: echoi
Date: 2010-04-06 12:35:30 -0700 (Tue, 06 Apr 2010)
New Revision: 16495

Modified:
   long/3D/SNAC/trunk/Snac/plugins/plasticSPR/Constitutive.c
   long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.c
Log:

Print some more useful info when frictionAngle is zero or negative. 
	- Do a quick fix if it's because of remeshed plastic strain is negative (simply set to zero).
	- Otherwise print and die.



Modified: long/3D/SNAC/trunk/Snac/plugins/plasticSPR/Constitutive.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plasticSPR/Constitutive.c	2010-04-03 00:05:23 UTC (rev 16494)
+++ long/3D/SNAC/trunk/Snac/plugins/plasticSPR/Constitutive.c	2010-04-06 19:35:30 UTC (rev 16495)
@@ -171,8 +171,19 @@
 				}
 			}
 			else {
-				/* frictionAngle < 0.0 violates second law of thermodynamics */
-				abort();
+				if( plasticStrain < 0.0 ) {
+					plasticElement->plasticStrain[tetra_I] = 0.0;
+					plasticStrain = plasticElement->plasticStrain[tetra_I];
+					fprintf(stderr,"Warning: negative plastic strain. Setting to zero, but check if remesher is on and this happended for an external tet. rank:%d elem:%d tet:%d plasticStrain=%e frictionAngle=%e\n",context->rank,element_lI,tetra_I,plasticStrain,frictionAngle);
+					frictionAngle = material->frictionAngle[0];
+					dilationAngle = material->dilationAngle[0];
+					cohesion = material->cohesion[0];
+				}
+				else {
+					/* frictionAngle < 0.0 violates second law of thermodynamics */
+					fprintf(stderr,"Error due to an unknown reason: rank:%d elem:%d tet:%d plasticStrain=%e frictionAngle=%e\n",context->rank,element_lI,tetra_I,plasticStrain,frictionAngle);
+					assert(0);
+				}
 			}
 			
 			
@@ -225,6 +236,9 @@
 				depm = ( dep1 + dep2 + dep3 ) / 3.0f;
 				plasticElement->plasticStrain[tetra_I] += sqrt( 0.5f * ((dep1-depm) * (dep1-depm) + (dep2-depm) * (dep2-depm) + (dep3-depm) * (dep3-depm) + depm*depm) );
 
+				/* linear healing */
+				plasticElement->plasticStrain[tetra_I] *= (1.0/(1.0+context->dt/3.0e+13));
+
 				memset( stress, 0, sizeof((*stress)) );
 				/* Resolve back to global axes  */
 				for( m = 0; m < 3; m++ ) {

Modified: long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.c	2010-04-03 00:05:23 UTC (rev 16494)
+++ long/3D/SNAC/trunk/Snac/plugins/viscoplasticSPR/Constitutive.c	2010-04-06 19:35:30 UTC (rev 16495)
@@ -286,9 +286,19 @@
 				}
 			}
 			else {
-				/* frictionAngle < 0.0 violates second law of thermodynamics */
-				fprintf(stderr,"rank:%d elem:%d tet:%d plasticStrain=%e\n",element_lI,tetra_I,plasticStrain);
-				assert(0);
+				if( plasticStrain < 0.0 ) {
+					viscoplasticElement->plasticStrain[tetra_I] = 0.0;
+					plasticStrain = viscoplasticElement->plasticStrain[tetra_I];
+					fprintf(stderr,"Warning: negative plastic strain. Setting to zero, but check if remesher is on and this happended for an external tet. rank:%d elem:%d tet:%d plasticStrain=%e frictionAngle=%e\n",context->rank,element_lI,tetra_I,plasticStrain,frictionAngle);
+					frictionAngle = material->frictionAngle[0];
+					dilationAngle = material->dilationAngle[0];
+					cohesion = material->cohesion[0];
+				}
+				else {
+					/* frictionAngle < 0.0 violates second law of thermodynamics */
+					fprintf(stderr,"Error due to an unknown reason: rank:%d elem:%d tet:%d plasticStrain=%e frictionAngle=%e\n",context->rank,element_lI,tetra_I,plasticStrain,frictionAngle);
+					assert(0);
+				}
 			}
 
 			if( yieldcriterion == mohrcoulomb )
@@ -344,6 +354,9 @@
 						depm = ( dep1 + dep2 + dep3 ) / 3.0f;
 						viscoplasticElement->plasticStrain[tetra_I] += sqrt( 0.5f * ((dep1-depm) * (dep1-depm) + (dep2-depm) * (dep2-depm) + (dep3-depm) * (dep3-depm) + depm*depm) );
 
+						/* linear healing */
+						viscoplasticElement->plasticStrain[tetra_I] *= (1.0/(1.0+context->dt/3.0e+13));
+
 						/* Stress projection back to euclidean coordinates */
 						memset( stress, 0, sizeof((*stress)) );
 						/* Resolve back to global axes  */