Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15472) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15473) @@ -2308,5 +2308,5 @@ /*Compute the length of this edge and compare it to the minimal length*/ - length=pow(pow(xyz_list[node0][0]-xyz_list[node1][0],2.0)+pow(xyz_list[node0][1]-xyz_list[node1][1],2.0)+pow(xyz_list[node0][2]-xyz_list[node1][2],2.0),0.5); + length=sqrt(pow(xyz_list[node0][0]-xyz_list[node1][0],2)+pow(xyz_list[node0][1]-xyz_list[node1][1],2)+pow(xyz_list[node0][2]-xyz_list[node1][2],2)); if(lengthvalues, &Ke_temp[0][0]); + //Ke->Echo(); + //_error_("stop"); /*Transform Coordinate System*/ @@ -6915,4 +6917,15 @@ GaussPenta *gauss=NULL; + /*Stabilization*/ + bool stabilization = false; + IssmDouble dbasis[3][6]; + IssmDouble dmu[3]; + IssmDouble mu; + IssmDouble dnodalbasis[6][6][3]; + IssmDouble SU[6][4][4]; + IssmDouble SW[6][4][4]; + int p,q,ii; + int c=3; //index of pressure + /*If on water or not Stokes, skip stiffness: */ inputs->GetInputValue(&approximation,ApproximationEnum); @@ -6932,5 +6945,5 @@ /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); + gauss=new GaussPenta(3,2); for(int ig=gauss->begin();igend();ig++){ @@ -6951,27 +6964,91 @@ &D[0][0],8,8,0, &B_prime[0][0],8,24,0, - &Ke_temp[0][0],1); + &Ke_temp[0][0],0); 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,6,1, - &D_stab[0][0],3,3,0, - &B_stab[0][0],3,6,0, - &Ke_temp_stab[0][0],1); - - for(i=0;ivalues[numdof*(i*4+3)+j*4+3]+=Ke_temp_stab[i][j]; - - } - + if(stabilization){ + GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + dmu[0]=0.; dmu[1]=0.; dmu[2]=0.; + mu = viscosity; + for(i=0;i<6;i++){ + for(p=0;p<6;p++){ + for(j=0;j<3;j++){ + dnodalbasis[i][p][j] = dbasis[j][p]; + } + } + } + for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){ + SU[p][i][j]=0.; + SW[p][i][j]=0.; + } + for(p=0;p<6;p++){ + for(i=0;i<3;i++){ + SU[p][i][c] += dbasis[i][p]; + SW[p][c][i] += dbasis[i][p]; + for(j=0;j<3;j++){ + SU[p][i][i] += -dmu[j]*dbasis[j][p]; + SU[p][i][j] += -dmu[i]*dbasis[i][p]; + for(ii=0;ii<6;ii++){ + SU[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; + SU[p][i][j] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; + } + SW[p][i][i] += -dmu[j]*dbasis[j][p]; + SW[p][j][i] += -dmu[j]*dbasis[i][p]; + for(ii=0;ii<6;ii++){ + SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; + SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; + } + } + } + } + IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity); + for(p=0;p<6;p++){ + for(q=0;q<6;q++){ + for(i=0;i<4;i++){ + for(j=0;j<4;j++){ + Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][0]*SU[q][0][j]; + Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][1]*SU[q][1][j]; + Ke->values[(p*4+i)*numdof+q*4+j] += gauss->weight*Jdet*tau*SW[p][i][2]*SU[q][2][j]; + } + } + } + } + + } + + + //viscosity = 1000*rigidity; + //D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning; + //if(id==24){ + //printf("tau %g\n",1./3.*pow(diameter,2)/(8.*2.*viscosity)); + //} + //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,6,1, + // &D_stab[0][0],3,3,0, + // &B_stab[0][0],3,6,0, + // &Ke_temp_stab[0][0],0); + + //for(i=0;ivalues[numdof*(i*4+3)+j*4+3]+=Ke_temp_stab[i][j]; + + } + +// if(id==24){ +// _error_(""); +// } + + //Ke->Echo(); + //_error_("stop"); /*Transform Coordinate System*/ TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum); + if(id==24){ + //printarray(Ke->values,24,24); + } /*Clean up and return*/ @@ -7790,7 +7867,7 @@ for(i=0;iweight*l1l7[i]; - Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*l1l7[i]; - Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*l1l7[i]; - Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*l1l7[i]; + Pe_gaussian[i*NDOF4+0]+=rho_ice*forcex*Jdet*gauss->weight*l1l7[i]; + Pe_gaussian[i*NDOF4+1]+=rho_ice*forcey*Jdet*gauss->weight*l1l7[i]; + Pe_gaussian[i*NDOF4+2]+=rho_ice*forcez*Jdet*gauss->weight*l1l7[i]; } @@ -7828,5 +7905,5 @@ int i,j; int approximation; - IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab; + IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab,viscosity; IssmDouble forcex,forcey,forcez,diameter,stokesreconditioning; IssmDouble xyz_list[NUMVERTICES][3]; @@ -7834,6 +7911,15 @@ 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; + + /*Stabilization*/ + bool stabilization = false; + IssmDouble dbasis[3][6]; + IssmDouble dmu[3]; + IssmDouble mu; + IssmDouble dnodalbasis[6][6][3]; + IssmDouble SW[6][4][4]; + int p,q,ii; + int c=3; //index of pressure /*Initialize Element vector and return if necessary*/ @@ -7851,4 +7937,7 @@ Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); + 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); /*Find minimal length*/ @@ -7856,5 +7945,5 @@ /* Start looping on the number of gaussian points: */ - gauss=new GaussPenta(5,5); + gauss=new GaussPenta(2,3); for(int ig=gauss->begin();igend();ig++){ @@ -7863,4 +7952,6 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); GetNodalFunctionsP1(&l1l6[0], gauss); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + material->GetViscosity3dStokes(&viscosity,&epsilon[0]); loadingforcex_input->GetInputValue(&forcex, gauss); @@ -7869,26 +7960,67 @@ for(i=0;iweight*l1l6[i]; - Pe_gaussian[i*NDOF4+0]+=forcex*Jdet*gauss->weight*l1l6[i]; - Pe_gaussian[i*NDOF4+1]+=forcey*Jdet*gauss->weight*l1l6[i]; - Pe_gaussian[i*NDOF4+2]+=forcez*Jdet*gauss->weight*l1l6[i]; + pe->values[i*NDOF4+2]+=-rho_ice*gravity*Jdet*gauss->weight*l1l6[i]; + pe->values[i*NDOF4+0]+=rho_ice*forcex*Jdet*gauss->weight*l1l6[i]; + pe->values[i*NDOF4+1]+=rho_ice*forcey*Jdet*gauss->weight*l1l6[i]; + pe->values[i*NDOF4+2]+=rho_ice*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;ivalues[i]+=Pe_gaussian[i]; + if(stabilization){ + GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss); + dmu[0]=0.; dmu[1]=0.; dmu[2]=0.; + mu = viscosity; + for(i=0;i<6;i++){ + for(p=0;p<6;p++){ + for(j=0;j<3;j++){ + dnodalbasis[i][p][j] = dbasis[j][p]; + } + } + } + for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){ + SW[p][i][j]=0.; + } + for(p=0;p<6;p++){ + for(i=0;i<3;i++){ + SW[p][c][i] += dbasis[i][p]; + for(j=0;j<3;j++){ + SW[p][i][i] += -dmu[j]*dbasis[j][p]; + SW[p][j][i] += -dmu[j]*dbasis[i][p]; + for(ii=0;ii<6;ii++){ + SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii]; + SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii]; + } + } + } + } + IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity); + for(p=0;p<6;p++){ + for(j=0;j<4;j++){ + pe->values[p*4+j] += gauss->weight*Jdet*tau*forcex*SW[p][j][0]; + pe->values[p*4+j] += gauss->weight*Jdet*tau*forcey*SW[p][j][1]; + pe->values[p*4+j] += gauss->weight*Jdet*tau*(forcez-rho_ice*gravity)*SW[p][j][2]; + } + } + + } + ///*Add stabilization*/ + //D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning; + //GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0],gauss); + + //for(i=0;ivalues,24,1); + _error_("STOP"); + } /*Clean up and return*/