Changeset 18282

Timestamp:

07/22/14 13:16:55 (11 years ago)

Author:

cborstad

Message:

NEW: updated damage solution can now be called as a Requested Output of a Stressbalance solution. Requires that the matdamageice class be used with two extra parameters specified (stress_threshold and c1, a ductility parameter).

File:

• issm/trunk-jpl/src/c/classes/Elements/Element.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -88 +88 @@
       B_input=this->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
    }
-        /* Get B and D inputs appropriate for the type of domain */
-        /*n=this->material->GetN();*/
-   /*if(domaintype==Domain2DhorizontalEnum){
-                B=this->material->GetBbar();
-                D=this->material->GetDbar();
-   }
-        else{
-                B=this->material->GetB();
-                D=this->material->GetD();
-        }*/
         
         /* Start looping on the number of nodes: */
@@ -121 +111 @@
       damage_input->GetInputValue(&D,gauss);
 
+                /* Compute threshold strain rate from threshold stress */
+                eps_0=pow(sigma_0/B,n); 
+                _assert_(eps_f>eps_0);
+
                 /* Compute kappa (k) from pre-existing level of damage, using Newton-Raphson */
-                k1=exp(n*eps_0/(eps_f-eps_0)+n*log(1.-D)-log(eps_0)); /* initial guess */
-                
+                /* provide a reasonable initial guess */
+                if(D==0){
+                        k1=eps_0;
+                }
+                else{
+                        k1=exp(n*eps_0/(eps_f-eps_0)-n*log(1.-D)+log(eps_0)); /* initial guess */
+                }
+        
+                counter=0;
                 while(true){
         /*Newton step k2=k1-f(k1)/f'(k1) */
-        k2=k1-(1.-D-pow(eps_0/k1,1./n)*exp(-(k1-eps_0)/(eps_f-eps_0)))/(pow(eps_0/k1,1./n)*exp(-(k1-eps_0)/(eps_f-eps_0))*(-1./(eps_f-eps_0)-1./n/k1));
+                        k2=k1-(k1+(eps_f-eps_0)/n*log(k1)-eps_0+(eps_f-eps_0)*(log(1.-D)-1./n*log(eps_0)))/(1.+(eps_f-eps_0)/n/k1);
 
         if( fabs(k2-k1)/(fabs(k2))<threshold ){
@@ -132 +133 @@
         }
         else{
-           k1=k2;
+                                k1=k2;
            counter++;
         }
@@ -139 +140 @@
         }
 
-                /* Compute threshold strain rate from threshold stress */
-                eps_0=pow(sigma_0/B,n); 
-                _assert_(eps_f>eps_0);
-
                 if(eps_eff>k2){
-                        newD[i]=1.-pow(eps_0/k2,1./n)*exp(-(k2-eps_0)/(eps_f-eps_0));
+                        newD[i]=1.-pow(eps_0/eps_eff,1./n)*exp(-(eps_eff-eps_0)/(eps_f-eps_0));
                 }
                 else newD[i]=D;