Changeset 11373

Timestamp:

02/08/12 16:49:03 (13 years ago)

Author:

cborstad

Message:

matice changes for new damage variable

Location:

issm/trunk-jpl-damage/src/c/objects/Materials

Files:

• Matice.cpp (modified) (17 diffs)
• Matice.h (modified) (1 diff)

issm/trunk-jpl-damage/src/c/objects/Materials/Matice.cpp

@@ -218 +218 @@
 }
 /*}}}*/
+/*FUNCTION Matice::GetZ {{{1*/
+double Matice::GetZ(){
+
+        /*Output*/
+        double Z;
+
+        inputs->GetInputAverage(&Z,MaterialsRheologyZEnum);
+        return Z;
+}
+/*}}}*/
+/*FUNCTION Matice::GetZbar {{{1*/
+double Matice::GetZbar(){
+
+        /*Output*/
+        double Zbar;
+
+        inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum);
+        return Zbar;
+}
+/*}}}*/
 /*FUNCTION Matice::GetN {{{1*/
 double Matice::GetN(){
@@ -263 +283 @@
 void  Matice::GetViscosity2d(double* pviscosity, double* epsilon){
         /*From a string tensor and a material object, return viscosity, using Glen's flow law.
-                                                                                                    B
+                                                                                                  Z * B
           viscosity= -------------------------------------------------------------------
                                                   2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
 
-          where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity 
-          vector, and n the flow law exponent.
+          where viscosity is the viscosity, Z the damage effect variable (1-D), 
+          B the flow law parameter, (u,v) the velocity vector, and n the flow law exponent.
 
           If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we 
@@ -282 +302 @@
         /*Intermediary: */
         double A,e;
-        double B,n;
-
-        /*Get B and n*/
-        B=GetBbar();
+        double Btmp,B,n,Z;
+
+        /*Get B, Z and n*/
+        Btmp=GetBbar();
+        Z=GetZbar();
         n=GetN();
+        B=Z*Btmp;
 
         if (n==1){
@@ -318 +340 @@
         if(viscosity<=0) _error_("Negative viscosity");
         _assert_(B>0);
+        _assert_(Z>0);
         _assert_(n>0);
 
@@ -329 +352 @@
         /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
          *
-         *                                               B
+         *                                             Z * B
          * viscosity3d= -------------------------------------------------------------------
          *                      2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
@@ -348 +371 @@
         /*Intermediaries: */
         double A,e;
-        double B,n;
+        double Btmp,B,n,Z;
 
         /*Get B and n*/
-        B=GetB();
+        Btmp=GetB();
+        Z=GetZ();
         n=GetN();
+        B=Z*Btmp;
 
         if (n==1){
@@ -389 +414 @@
         if(viscosity3d<=0) _error_("Negative viscosity");
         _assert_(B>0);
+        _assert_(Z>0);
         _assert_(n>0);
 
@@ -399 +425 @@
         /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
          *
-         *                                          B
+         *                                        Z * B
          * viscosity3d= -------------------------------------------------------------------
          *                   2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
@@ -418 +444 @@
         /*Intermediaries: */
         double A,e;
-        double B,n;
+        double Btmp,B,Z,n;
         double eps0;
 
         /*Get B and n*/
         eps0=pow((double)10,(double)-27);
-        B=GetB();
-        n=GetN();
+        Btmp=GetB();
+        Z=GetZ();
+        n=GetB();
+        B=Z*Btmp;
         
         if (n==1){
@@ -461 +489 @@
         if(viscosity3d<=0) _error_("Negative viscosity");
         _assert_(B>0);
+        _assert_(Z>0);
         _assert_(n>0);
 
@@ -490 +519 @@
 
         /*Get B and n*/
-        B=GetBbar();
+        B=GetBbar(); /* why is this needed in this function? */
         n=GetN();
 
@@ -508 +537 @@
                 
                         viscosity_complement=1/(2*pow(A,e));
+                }
+        }
+        else{
+                viscosity_complement=4.5*pow((double)10,(double)17);
+        }
+
+        /*Checks in debugging mode*/
+        _assert_(B>0); /* shouldn't be necessary as B is never used! */
+        _assert_(n>0);
+        _assert_(viscosity_complement>0);
+                
+        /*Return: */
+        *pviscosity_complement=viscosity_complement;
+}
+/*}}}*/
+/*FUNCTION Matice::GetViscosityZComplement {{{1*/
+void  Matice::GetViscosityZComplement(double* pviscosity_complement, double* epsilon){
+        /*Return viscosity derivative for control method d(mu)/dZ: 
+         *
+         *                                                                                             B 
+         * dviscosity= -------------------------------------------------------------------
+         *                                2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
+         *
+         * If epsilon is NULL, it means this is the first time Gradjb is being run, and we 
+         * return mu20, initial viscosity.
+         */
+        
+        /*output: */
+        double viscosity_complement;
+
+        /*input strain rate: */
+        double exx,eyy,exy;
+
+        /*Intermediary value A and exponent e: */
+        double A,e;
+        double B,n;
+
+        /*Get B and n*/
+        B=GetBbar();
+        n=GetN();
+
+        if(epsilon){
+                exx=*(epsilon+0);
+                eyy=*(epsilon+1);
+                exy=*(epsilon+2);
+
+                /*Build viscosity: mu2=B/(2*A^e) */
+                A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
+                if(A==0){
+                        /*Maximum viscosity_complement for 0 shear areas: */
+                        viscosity_complement=2.25*pow((double)10,(double)17);
+                }
+                else{
+                        e=(n-1)/(2*n);
+                
+                        viscosity_complement=B/(2*pow(A,e));
                 }
         }
@@ -685 +770 @@
                         this->inputs->AddInput(new TriaP1Input(MaterialsRheologyBbarEnum,nodeinputs));
                 }
+                
+                /*Get Z*/
+                if (iomodel->Data(MaterialsRheologyZEnum)) {
+                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
+                        this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs));
+                }
 
                 /*Get n*/
@@ -706 +797 @@
                                                 }
                                                 break;
+                                        case MaterialsRheologyZbarEnum:
+                                                if (iomodel->Data(MaterialsRheologyZEnum)){
+                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
+                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
+                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
+                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
+                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                                }
+                                                break;
                                 }
                         }
@@ -726 +826 @@
                         for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyBEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
                         this->inputs->AddInput(new PentaP1Input(MaterialsRheologyBEnum,nodeinputs));
+                }
+                
+                /*Get Z*/
+                if (iomodel->Data(MaterialsRheologyZEnum)) {
+                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
+                        this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs));
                 }
 
@@ -748 +854 @@
                                                 }
                                                 break;
+                                        case MaterialsRheologyZbarEnum:
+                                                if (iomodel->Data(MaterialsRheologyZEnum)){
+                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
+                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
+                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
+                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
+                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                                }
+                                                break;
+
                                 }
                         }
@@ -766 +882 @@
                                 name==MaterialsRheologyBEnum ||
                                 name==MaterialsRheologyBbarEnum ||
+                                name==MaterialsRheologyZEnum ||
+                                name==MaterialsRheologyZbarEnum ||
                                 name==MaterialsRheologyNEnum
                 ){

issm/trunk-jpl-damage/src/c/objects/Materials/Matice.h

@@ -68 +68 @@
                 void   GetViscosity3dStokes(double* pviscosity3d, double* epsilon);
                 void   GetViscosityComplement(double* pviscosity_complement, double* pepsilon);
+                void   GetViscosityZComplement(double* pviscosity_complement, double* pepsilon);
                 double GetB();
                 double GetBbar();
+                double GetZ();
+                double GetZbar();
                 double GetN();
                 bool   IsInput(int name);