Changeset 18277

Timestamp:

07/21/14 17:43:07 (11 years ago)

Author:

cborstad

Message:

CHG: working on new damage calculation as a requested output of a Stressbalance solution

Location:

issm/trunk-jpl/src/c

Files:

• analyses/DamageEvolutionAnalysis.cpp (modified) (1 diff)
• analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)
• classes/Elements/Element.cpp (modified) (3 diffs)

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

@@ -548 +548 @@
 
         /*Compute stress tensor: */
-        element->ComputeDeviatoricStressTensor();
+        element->ComputeStrainRate();
 
         /*retrieve what we need: */

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

@@ -89 +89 @@
         int     numoutputs;
         char**  requestedoutputs = NULL;
+        int     materials_type;
 
         parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
@@ -114 +115 @@
                 parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
                 parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
+        }
+
+        iomodel->Constant(&materials_type,MaterialsEnum);
+        if(materials_type==MatdamageiceEnum){
+                parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum));
+                parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum));
         }
 

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

@@ -46 +46 @@
         IssmDouble  eps_xx,eps_xy,eps_yy,eps_xz,eps_yz,eps_zz,eps_eff;
         IssmDouble  epsmin=1.e-27;
-        int         dim;
-
-        /*Retrieve all inputs we will be needing: */
-        /* TODO: retrieve parameters eps_0 and eps_f and input DamageD(bar?) */
+        IssmDouble  eps_0,eps_f,sigma_0,B,D,n;
+        int         dim,counter=0;
+        IssmDouble  k1,k2,threshold=1.e-12;
+
+        /* Retrieve parameters */
         this->GetVerticesCoordinates(&xyz_list);
         this->ComputeStrainRate();
         parameters->FindParam(&dim,DomainDimensionEnum);
+        parameters->FindParam(&eps_f,DamageC1Enum);
+        parameters->FindParam(&sigma_0,DamageStressThresholdEnum);
+
+        /* Retrieve inputs */
         Input* eps_xx_input=this->GetInput(StrainRatexxEnum); _assert_(eps_xx_input);
         Input* eps_yy_input=this->GetInput(StrainRateyyEnum); _assert_(eps_yy_input);
@@ -65 +70 @@
         }
 
-        /* Start looping on the number of vertices: */
+        /* Fetch number of nodes and allocate output*/
+   int numnodes = this->GetNumberOfNodes();
+   IssmDouble* newD = xNew<IssmDouble>(numnodes);
+
+        /* Retrieve domain-dependent inputs */
+        Input* n_input=this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
+   Input* damage_input = NULL;
+   Input* B_input = NULL;
+        int domaintype;
+   parameters->FindParam(&domaintype,DomainTypeEnum);
+        if(domaintype==Domain2DhorizontalEnum){
+      damage_input = this->GetInput(DamageDbarEnum);  _assert_(damage_input);
+      B_input=this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
+   }
+   else{
+      damage_input = this->GetInput(DamageDEnum);   _assert_(damage_input);
+      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: */
         Gauss* gauss=this->NewGauss();
-        int numvertices = this->GetNumberOfVertices();
-        for (int iv=0;iv<numvertices;iv++){
-                gauss->GaussVertex(iv);
+        for (int i=0;i<numnodes;i++){
+                gauss->GaussNode(this->GetElementType(),i);
 
                 eps_xx_input->GetInputValue(&eps_xx,gauss);
@@ -82 +115 @@
 
                 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
-                eps_eff=sqrt(eps_xx*eps_xx+eps_yy*eps_yy+eps_xy*eps_xy+eps_xz*eps_xz+eps_yz*eps_yz+eps_xx*eps_xx*epsmin*epsmin);
-
-                /*TODO: compute kappa from initial D, then compute new D */
-
-        }
-
-        /* TODO: add newdamage input to DamageEnum and NewDamageEnum */
+                eps_eff=sqrt(eps_xx*eps_xx+eps_yy*eps_yy+eps_xy*eps_xy+eps_xz*eps_xz+eps_yz*eps_yz+eps_xx*eps_yy+epsmin*epsmin);
+
+                B_input->GetInputValue(&B,gauss);
+      n_input->GetInputValue(&n,gauss);
+      damage_input->GetInputValue(&D,gauss);
+
+                /* 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 */
+                
+                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));
+
+        if( fabs(k2-k1)/(fabs(k2))<threshold ){
+                break;
+        }
+        else{
+           k1=k2;
+           counter++;
+        }
+
+        if(counter>50) break;
+        }
+
+                /* 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));
+                }
+                else newD[i]=D;
+        }
+
+        /* Add new damage input to DamageEnum and NewDamageEnum */
+        this->AddInput(NewDamageEnum,newD,this->GetElementType());
+        if(domaintype==Domain2DhorizontalEnum){
+                this->AddInput(DamageDbarEnum,newD,this->GetElementType());
+        }
+        else{
+                this->AddInput(DamageDEnum,newD,this->GetElementType());
+        }
 
         /*Clean up and return*/
         xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(newD);
         delete gauss;