Changeset 11341

Timestamp:

02/07/12 09:07:26 (13 years ago)

Author:

Mathieu Morlighem

Message:

Almost done with Jacobian Stokes

File:

• issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (modified) (7 diffs)

issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

@@ -7437 +7437 @@
                 case StokesApproximationEnum:
                         return CreateJacobianDiagnosticStokes();
+                case NoneApproximationEnum:
+                        return NULL;
                 default:
                         _error_("Approximation %s not supported yet",EnumToStringx(approximation));
@@ -7473 +7475 @@
         /*Intermediaries */
         int        i,j,ig;
-        int        approximation;
         double     xyz_list[NUMVERTICES][3];
         double     Jdet;
@@ -7489 +7490 @@
 
         /*Retrieve all inputs and parameters*/
-        inputs->GetInputValue(&approximation,ApproximationEnum);
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
@@ -7535 +7535 @@
 ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
 
-        _error_("Not implemented yet");
         /*Constants*/
-        const int    numdof=NDOF2*NUMVERTICES;
+        const int    numdof=NDOF4*NUMVERTICES;
 
         /*Intermediaries */
         int        i,j,ig;
-        int        approximation;
         double     xyz_list[NUMVERTICES][3];
         double     Jdet;
         double     eps1dotdphii,eps1dotdphij;
         double     eps2dotdphii,eps2dotdphij;
+        double     eps3dotdphii,eps3dotdphij;
         double     mu_prime;
         double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
-        double     eps1[3],eps2[3];
+        double     eps1[3],eps2[3],eps3[3];
         double     phi[NUMVERTICES];
         double     dphi[3][NUMVERTICES];
         GaussPenta *gauss=NULL;
 
-        /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
-        ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
+        /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
+        ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
 
         /*Retrieve all inputs and parameters*/
-        inputs->GetInputValue(&approximation,ApproximationEnum);
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
+        Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
 
         /* Start  looping on the number of gaussian points: */
@@ -7573 +7572 @@
                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
                 matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
-                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
-                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
-                eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
+                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
+                eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
+                eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
 
                 for(i=0;i<6;i++){
@@ -7583 +7582 @@
                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
-
-                                Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
-                                Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
-                                Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
-                                Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
+                                eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
+                                eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
+
+                                Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
+                                Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
+                                Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
+
+                                Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
+                                Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
+                                Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
+
+                                Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
+                                Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
+                                Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
                         }
                 }
@@ -7593 +7601 @@
 
         /*Transform Coordinate System*/
-        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
+        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
 
         /*Clean up and return*/