Changeset 89

Timestamp:

04/28/09 15:15:45 (16 years ago)

Author:

Eric.Larour

Message:

New viscosity 3d for Pattyn elements

File:

• issm/trunk/src/c/objects/Matice.cpp (modified) (1 diff)

issm/trunk/src/c/objects/Matice.cpp

@@ -240 +240 @@
         *pviscosity2=viscosity2;
 }
+
+#undef __FUNCT__ 
+#define __FUNCT__ "MatIce::GetViscosity3d"
+void  Matice::GetViscosity3d(double* pviscosity3d, double* epsilon){
+
+        /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: 
+         *
+         *                                 2*B
+         * viscosity3d= -------------------------------------------------------------------
+         *     2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
+         *
+         *     where mu is the viscotiy, 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 Emg is being run, and we 
+         * return g, initial viscosity.
+         */
+        
+        /*output: */
+        double viscosity3d;
+
+        /*input strain rate: */
+        double exx,eyy,exy,exz,eyz;
+
+        /*Intermediary value A and exponent e: */
+        double A,e;
+
+        if (n==1){
+                /*Viscous behaviour! viscosity3d=B: */
+                viscosity3d=B;
+        }
+        else{
+                if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
+                                (epsilon[3]==0) && (epsilon[4]==0)){
+                        viscosity3d=pow(10,14);
+                }
+                else{
+
+                        /*Retrive strain rate components: */
+                        exx=*(epsilon+0);
+                        eyy=*(epsilon+1);
+                        exy=*(epsilon+2);
+                        exz=*(epsilon+3);
+                        eyz=*(epsilon+4);
+
+                        /*Build viscosity: viscosity3d=2*B/(2*A^e) */
+                        A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy;
+                        if(A==0){
+                                /*Maxiviscosity3dm viscosity for 0 shear areas: */
+                                viscosity3d=4.5*pow(10,17);
+                        }
+                        else{
+                                e=(n-1)/2/n;
+                        
+                                viscosity3d=2*B/(2*pow(A,e));
+                        }
+                }
+        }
+        #ifdef _DEBUG_
+        _printf_("Viscosity %lf\n",viscosity3d);
+        #endif
+
+        /*Assign output pointers:*/
+        *pviscosity3d=viscosity3d;
+}