Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15437) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15438) @@ -6894,16 +6894,20 @@ 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,4 +6924,8 @@ 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,4 +6951,20 @@ for(i=0;ivalues[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;ivalues[i*numdof*4+3+j*4+3]+=Ke_temp_stab[i][j]; + } @@ -7781,57 +7805,4 @@ } - 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;ibegin();igend();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,9 +7825,10 @@ 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,4 +7837,5 @@ 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,4 +7843,5 @@ 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,4 +7852,7 @@ 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,4 +7873,15 @@ 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