Changeset 7000

Timestamp:

01/07/11 16:34:33 (14 years ago)

Author:

seroussi

Message:

finished coupling viscous matrix MacAyeal/Stokes

Location:

issm/trunk/src/c/objects

Files:

• Elements/Penta.cpp (modified) (8 diffs)
• Elements/PentaRef.cpp (modified) (2 diffs)
• Node.cpp (modified) (1 diff)

issm/trunk/src/c/objects/Elements/Penta.cpp ¶

@@ -719 +719 @@
         const int    numdofm=NDOF2*NUMVERTICES2D;
         const int    numdofs=NDOF4*NUMVERTICES;
-        const int    numdoftotal=numdofm+numdofs;
+        const int    numdoftotal=2*numdofm+numdofs;
 
         /*Intermediaries */
         int         i,j,ig;
         double      Jdet;
-        double      viscosity; //viscosity
-        double      epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+        double      viscosity,stokesreconditioning; //viscosity
+        double      epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
         double      xyz_list[NUMVERTICES][3];
-        double      B[4][numdofs];
+        double      B[4][numdofs+3];
         double      Bprime[4][numdofm];
         double      B2[3][numdofm];
-        double      Bprime2[3][numdofs];
+        double      Bprime2[3][numdofs+3];
         double      D[4][4]={0.0};            // material matrix, simple scalar matrix.
         double      D2[3][3]={0.0};            // material matrix, simple scalar matrix.
@@ -736 +736 @@
         double      Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix 
         double      Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix 
-        double      Ke_gg_gaussian[numdofs][numdofm]; //stiffness matrix evaluated at the gaussian point.
-        double      Ke_gg_gaussian2[numdofm][numdofs]; //stiffness matrix evaluated at the gaussian point.
+        double      Ke_gg_gaussian[numdofs+3][numdofm]; //stiffness matrix evaluated at the gaussian point.
+        double      Ke_gg_gaussian2[numdofm][numdofs+3]; //stiffness matrix evaluated at the gaussian point.
         GaussPenta *gauss=NULL;
         GaussTria  *gauss_tria=NULL;
@@ -754 +754 @@
         /* Get node coordinates and dof list: */
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
@@ -777 +778 @@
                 D_scalar=2*viscosity*gauss->weight*Jdet;
                 for (i=0;i<3;i++) D[i][i]=D_scalar;
-                D[4][4]=(-1)*gauss->weight*Jdet;
+                D[3][3]=-gauss->weight*Jdet*stokesreconditioning;
                 for (i=0;i<3;i++) D2[i][i]=D_scalar;
 
-                TripleMultiply( &B[0][0],4,numdofs,1,
+                TripleMultiply( &B[0][0],4,numdofs+3,1,
                                         &D[0][0],4,4,0,
                                         &Bprime[0][0],4,numdofm,0,
@@ -787 +788 @@
                 TripleMultiply( &B2[0][0],3,numdofm,1,
                                         &D2[0][0],3,3,0,
-                                        &Bprime2[0][0],3,numdofs,0,
+                                        &Bprime2[0][0],3,numdofs+3,0,
                                         &Ke_gg_gaussian2[0][0],0);
 
@@ -794 +795 @@
         } 
         for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
-        for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[j][i];
-        _error_("coupling not finished yet");
+        for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j];
 
         /*Clean-up and return*/
@@ -990 +990 @@
         double      viscosity, oldviscosity, newviscosity, viscosity_overshoot;
         double      epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+        double      epsilons[6]; //6 for stokes
         double      xyz_list[NUMVERTICES][3];
         double      B[3][numdof2d];
@@ -1040 +1041 @@
                 }
                 else if (approximation==MacAyealStokesApproximationEnum){
-                        this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-                        matice->GetViscosity3dStokes(&newviscosity,&epsilon[0]);
+                        this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+                        matice->GetViscosity3dStokes(&newviscosity,&epsilons[0]);
                 }
                 else _error_("approximation %i (%s) not supported yet",approximation,EnumToString(approximation));

issm/trunk/src/c/objects/Elements/PentaRef.cpp ¶

@@ -134 +134 @@
                 *(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=l1l6[i];
         }
-
 }
 /*}}}*/
@@ -245 +244 @@
                 *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=dh1dh7[0][i];
                 *(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
-                *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=0;
-                *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=0;
-                *(Bprime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
         }
 

issm/trunk/src/c/objects/Node.cpp ¶

@@ -104 +104 @@
                         if (iomodel->vertices_type[io_index]==MacAyealStokesApproximationEnum && iomodel->gridonmacayeal[io_index]){
                                 if(!iomodel->gridonbed[io_index]){
-                                        for(k=1;k<=gsize;k++) this->FreezeDof(k);
+                                        for(k=1;k<=2;k++) this->FreezeDof(k);
                                 }
                         }