Changeset 15438

Timestamp:

07/05/13 12:50:39 (12 years ago)

Author:

seroussi

Message:

FIX: fixed stokes

File:

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

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

@@ -6894 +6894 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSViscous(void){
 
-        int        numdof = NDOF4*NUMVERTICES;
+        int        numdof  = NUMVERTICES*NDOF4;
+        int        numvert = NUMVERTICES;
 
         /*Intermediaries */
         int        i,j,approximation;
-        IssmDouble Jdet,viscosity,stokesreconditioning;
+        IssmDouble Jdet,viscosity,stokesreconditioning,diameter,rigidity;
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         IssmDouble B[8][24];
         IssmDouble B_prime[8][24];
-        IssmDouble D_scalar;
+        IssmDouble B_stab[3][numvert];
+        IssmDouble D_scalar,D_scalar_stab;
         IssmDouble D[8][8]={0.0};
+        IssmDouble D_stab[3][3]={0.0};
         IssmDouble Ke_temp[24][24]={0.0}; //for the six nodes
+        IssmDouble Ke_temp_stab[6][6]={0.0}; //for the six nodes
         GaussPenta *gauss=NULL;
 
@@ -6920 +6924 @@
         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
 
+        /*Find minimal length and B*/
+        rigidity=material->GetB();
+        diameter=MinEdgeLength(xyz_list);
+
         /* Start  looping on the number of gaussian points: */
         gauss=new GaussPenta(5,5);
@@ -6943 +6951 @@
 
                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_temp[i][j];
+
+                /*Add stabilization*/
+                D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2.)/(4.*1000*rigidity)*stokesreconditioning;
+                GetBConduct(&B_stab[0][0],&xyz_list[0][0],gauss); 
+
+                D_stab[0][0]=D_scalar_stab; D_stab[0][1]=0;             D_stab[0][2]=0;
+                D_stab[1][0]=0;             D_stab[1][1]=D_scalar_stab; D_stab[1][2]=0;
+                D_stab[2][0]=0;             D_stab[2][1]=0;             D_stab[2][2]=D_scalar_stab;
+
+                TripleMultiply(&B_stab[0][0],3,numvert,1,
+                                        &D_stab[0][0],3,3,0,
+                                        &B_stab[0][0],3,numvert,0,
+                                        &Ke_temp_stab[0][0],1);
+
+                for(i=0;i<numvert;i++) for(j=0;j<numvert;j++) Ke->values[i*numdof*4+3+j*4+3]+=Ke_temp_stab[i][j];
+
         }
 
@@ -7781 +7805 @@
         }
 
-        if(IsOnBed() || IsOnSurface()){
-                IssmDouble      xyz_list_tria[NUMVERTICES2D][3];
-                IssmDouble  pi=3.141592653589793;
-                IssmDouble  x,y,z;
-                IssmDouble  basis[6]; //for the six nodes of the penta
-                GaussPenta *gauss=NULL;
-                Input* x_input=inputs->GetInput(MeshXEnum);   _assert_(x_input);
-                Input* y_input=inputs->GetInput(MeshYEnum);   _assert_(y_input);
-                Input* z_input=inputs->GetInput(MeshZEnum);   _assert_(z_input);
-
-
-                if(IsOnBed()){ 
-                        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
-                        gauss=new GaussPenta(0,1,2,2);
-                }
-                if(IsOnSurface()){
-                        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i+3][j];
-                        gauss=new GaussPenta(3,4,5,2);
-                }
-
-
-                for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-                        gauss->GaussPoint(ig);
-
-                        GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0], gauss);
-                        GetNodalFunctionsP1(basis, gauss);
-                        x_input->GetInputValue(&x, gauss);
-                        y_input->GetInputValue(&y, gauss);
-                        z_input->GetInputValue(&z, gauss);
-                        
-
-                        forcex=-((cos(2*pi*z)-1)*2*pi*sin(2*pi*y)*cos(2*pi*x)-2*(cos(2*pi*x)-1)*2*pi*sin(2*pi*y)*cos(2*pi*z))*1;
-                        forcey=-((cos(2*pi*z)-1)*2*pi*sin(2*pi*x)*cos(2*pi*y)+(cos(2*pi*y)-1)*2*pi*sin(2*pi*x)*cos(2*pi*z))*1;
-                        forcez=-(-2*pi*2*sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z)+sin(2*pi*x)*sin(2*pi*y)*sin(2*pi*z));
-
-                        if(IsOnBed()){
-                                for(i=0;i<3;i++){
-                                        Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*basis[i]*-1.;
-                                        Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*basis[i]*-1.;
-                                        Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*basis[i]*-1.;
-                                }
-                        }
-                        if(IsOnSurface()){
-                                for(i=3;i<6;i++){
-                                        Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*basis[i]*1.;
-                                        Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*basis[i]*1.;
-                                        Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*basis[i]*1.;
-                                }
-                        }
-                }
-        }
-
         /*Condensation*/
         ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]);
@@ -7854 +7825 @@
         int        i,j;
         int        approximation;
-        IssmDouble Jdet,gravity,rho_ice;
-        IssmDouble forcex,forcey,forcez;
+        IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab;
+        IssmDouble forcex,forcey,forcez,diameter,stokesreconditioning;
         IssmDouble xyz_list[NUMVERTICES][3];
         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
         IssmDouble l1l6[6]; //for the six nodes and the bubble 
+        IssmDouble dh1dh6[3][NUMVERTICES];
         IssmDouble Pe_gaussian[numdof]={0.0}; //for the six nodes and the bubble 
         GaussPenta *gauss=NULL;
@@ -7865 +7837 @@
         inputs->GetInputValue(&approximation,ApproximationEnum);
         if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
+        parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
 
@@ -7870 +7843 @@
         rho_ice=matpar->GetRhoIce();
         gravity=matpar->GetG();
+        B=material->GetB();
         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         Input* vx_input=inputs->GetInput(VxEnum);   _assert_(vx_input);
@@ -7878 +7852 @@
         Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
 
+        /*Find minimal length*/
+        diameter=MinEdgeLength(xyz_list);
+
         /* Start  looping on the number of gaussian points: */
         gauss=new GaussPenta(5,5);
@@ -7896 +7873 @@
                         Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*l1l6[i];
                         Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*l1l6[i];
+                }
+
+                /*Add stabilization*/
+                D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2.)/(4.*1000*B)*stokesreconditioning;
+                GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0],gauss);
+
+                for(i=0;i<NUMVERTICES;i++){
+                        Pe_gaussian[i*NDOF4+3]+=-rho_ice*gravity*D_scalar_stab*dh1dh6[2][i];
+                        Pe_gaussian[i*NDOF4+3]+=forcex*D_scalar_stab*dh1dh6[0][i];
+                        Pe_gaussian[i*NDOF4+3]+=forcey*D_scalar_stab*dh1dh6[1][i];
+                        Pe_gaussian[i*NDOF4+3]+=forcez*D_scalar_stab*dh1dh6[2][i];
                 }
         }