Changeset 18136

Timestamp:

06/10/14 16:29:46 (11 years ago)

Author:

cborstad

Message:

CHG: working on a new damage evolution source term

Location:

issm/trunk-jpl/src/c/analyses

Files:

• DamageEvolutionAnalysis.cpp (modified) (5 diffs)
• DamageEvolutionAnalysis.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

@@ -266 +266 @@
 
         /*Intermediaries*/
-        int      domaintype;
+        int      domaintype,damagelaw;
         Element* basalelement;
         IssmDouble  Jdet,dt;
@@ -285 +285 @@
         }
 
+
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes = basalelement->GetNumberOfNodes();
@@ -295 +296 @@
         basalelement->GetVerticesCoordinates(&xyz_list);
         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-        this->CreateDamageFInput(basalelement);
+        basalelement->FindParam(&damagelaw,DamageLawEnum);
+        if(damagelaw==1 | damagelaw==2){
+                this->CreateDamageFInputPralong(basalelement);
+        }
+        else if(damagelaw==3){
+                this->CreateDamageFInputExp(basalelement);
+        }
         Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input);
         Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input);
@@ -440 +447 @@
 
 /*Intermediaries*/
-void DamageEvolutionAnalysis::CreateDamageFInput(Element* element){/*{{{*/
+void DamageEvolutionAnalysis::CreateDamageFInputPralong(Element* element){/*{{{*/
 
         /*Intermediaries */
@@ -527 +534 @@
         delete gauss;
 }/*}}}*/
+void DamageEvolutionAnalysis::CreateDamageFInputExp(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        IssmDouble epsf,stress_threshold,eps0;
+        IssmDouble damage,B,n,epseff;
+        IssmDouble eps_xx,eps_yy,eps_xy,eps1,eps2,epstmp;
+        int domaintype,damagelaw;
+
+        /*Fetch number of vertices and allocate output*/
+        int numnodes = element->GetNumberOfNodes();
+        IssmDouble* f   = xNew<IssmDouble>(numnodes);
+
+        /*retrieve parameters:*/
+        element->FindParam(&epsf,DamageC1Enum);
+        element->FindParam(&stress_threshold,DamageStressThresholdEnum);
+        element->FindParam(&domaintype,DomainTypeEnum);
+        element->FindParam(&damagelaw,DamageLawEnum);
+
+        /*Compute stress tensor: */
+        element->ComputeDeviatoricStressTensor();
+
+        /*retrieve what we need: */
+        Input* eps_xx_input  = element->GetInput(StrainRatexxEnum);     _assert_(eps_xx_input);
+        Input* eps_xy_input  = element->GetInput(StrainRatexyEnum);     _assert_(eps_xy_input);
+        Input* eps_yy_input  = element->GetInput(StrainRateyyEnum);     _assert_(eps_yy_input);
+        Input*  n_input=element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
+        Input* damage_input = NULL;
+        Input* B_input = NULL;
+        if(domaintype==Domain2DhorizontalEnum){
+                damage_input = element->GetInput(DamageDbarEnum);       _assert_(damage_input);
+                B_input=element->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
+        }
+        else{
+                damage_input = element->GetInput(DamageDEnum);   _assert_(damage_input);
+                B_input=element->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
+        }
+
+        /*Calculate damage evolution source term: */
+        Gauss* gauss=element->NewGauss();
+        for (int i=0;i<numnodes;i++){
+                gauss->GaussNode(element->GetElementType(),i);
+                
+                eps_xx_input->GetInputValue(&eps_xx,gauss);
+                eps_xy_input->GetInputValue(&eps_xy,gauss);
+                eps_yy_input->GetInputValue(&eps_yy,gauss);
+                B_input->GetInputValue(&B,gauss);
+                n_input->GetInputValue(&n,gauss);
+                damage_input->GetInputValue(&damage,gauss);
+        
+                /*Calculate principal effective strain rates*/
+                eps1=(eps_xx+eps_yy)/2.+sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
+                eps2=(eps_xx+eps_yy)/2.-sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
+                if(fabs(eps2)>fabs(eps1)){epstmp=eps2; eps2=eps1; eps1=epstmp;}
+
+                /*Calculate effective strain rate and threshold strain rate*/
+                epseff=1./sqrt(2.)*sqrt(eps1*eps1-eps1*eps2+eps2*eps2);
+                eps0=pow(stress_threshold/B,n);
+
+                if(epseff>eps0){
+                        f[i]=1.-pow(eps0/epseff,1./n)*exp(-(epseff-eps0)/(epsf-eps0))-damage;
+                }
+                else f[i]=0;
+        }
+
+        /*Add input*/
+        element->AddInput(DamageFEnum,f,element->GetElementType());
+        
+        /*Clean up and return*/
+        xDelete<IssmDouble>(f);
+        delete gauss;
+}/*}}}*/

issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h

@@ -34 +34 @@
 
                 /*Intermediaries*/
-                void CreateDamageFInput(Element* element);
+                void CreateDamageFInputPralong(Element* element);
+                void CreateDamageFInputExp(Element* element);
 };
 #endif