Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 10627) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 10628) @@ -822,4 +822,31 @@ /*return PentaRef field*/ return this->element_type; +} +/*}}}*/ +/*FUNCTION Penta::GetElementSizes{{{1*/ +void Penta::GetElementSizes(double* hx,double* hy,double* hz){ + + double xyz_list[NUMVERTICES][3]; + double xmin,ymin,zmin; + double xmax,ymax,zmax; + + /*Get xyz list: */ + GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); + xmin=xyz_list[0][0]; xmax=xyz_list[0][0]; + ymin=xyz_list[0][1]; ymax=xyz_list[0][1]; + zmin=xyz_list[0][2]; zmax=xyz_list[0][2]; + + for(int i=1;ixmax) xmax=xyz_list[i][0]; + if(xyz_list[i][1]ymax) ymax=xyz_list[i][1]; + if(xyz_list[i][2]zmax) zmax=xyz_list[i][2]; + } + + *hx=xmax-xmin; + *hy=ymax-ymin; + *hz=zmax-zmin; } /*}}}*/ @@ -3180,6 +3207,7 @@ K[0][0]=D_scalar_stab*pow(u,2); K[0][1]=D_scalar_stab*fabs(u)*fabs(v); K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2); - - GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); + _error_("TO BE RECODED"); + + //GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); TripleMultiply(&B_stab[0][0],2,numdof,1, @@ -3304,6 +3332,6 @@ int stabilization; int i,j,ig,found=0; - double Jdet,u,v,w,um,vm,wm; - double epsvel=2.220446049250313e-16; + double Jdet,u,v,w,um,vm,wm,vel; + double h,hx,hy,hz,vx,vy,vz; double gravity,rho_ice,rho_water; double heatcapacity,thermalconductivity,dt; @@ -3321,5 +3349,5 @@ double D_scalar_trans,D_scalar_stab; double D[3][3]; - double K[2][2]={0.0}; + double K[3][3]={0.0}; Tria* tria=NULL; GaussPenta *gauss=NULL; @@ -3374,17 +3402,14 @@ GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); - vx_input->GetInputValue(&u, gauss); - vy_input->GetInputValue(&v, gauss); - vz_input->GetInputValue(&w, gauss); - vxm_input->GetInputValue(&um,gauss); - vym_input->GetInputValue(&vm,gauss); - vzm_input->GetInputValue(&wm,gauss); + vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um; + vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm; + vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm; D_scalar_advec=gauss->weight*Jdet; if(dt) D_scalar_advec=D_scalar_advec*dt; - D[0][0]=D_scalar_advec*(u-um);D[0][1]=0; D[0][2]=0; - D[1][0]=0; D[1][1]=D_scalar_advec*(v-vm);D[1][2]=0; - D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_advec*(w-wm); + D[0][0]=D_scalar_advec*vx; D[0][1]=0; D[0][2]=0; + D[1][0]=0; D[1][1]=D_scalar_advec*vy; D[1][2]=0; + D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_advec*vz; TripleMultiply(&B_advec[0][0],3,numdof,1, @@ -3410,20 +3435,25 @@ if(stabilization==1){ /*Build K: */ - D_scalar_stab=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel); + GetElementSizes(&hx,&hy,&hz); + vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14; + h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.)); + + K[0][0]=h/(2*vel)*vx*vx; K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz; + K[1][0]=h/(2*vel)*vy*vx; K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz; + K[2][0]=h/(2*vel)*vz*vx; K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz; + + D_scalar_stab=gauss->weight*Jdet; if(dt) D_scalar_stab=D_scalar_stab*dt; - K[0][0]=D_scalar_stab*pow(u,2); K[0][1]=D_scalar_stab*fabs(u)*fabs(v); - K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2); - - GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); - - TripleMultiply(&B_stab[0][0],2,numdof,1, - &K[0][0],2,2,0, - &B_stab[0][0],2,numdof,0, + for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j]; + + GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); + + TripleMultiply(&Bprime_advec[0][0],3,numdof,1, + &K[0][0],3,3,0, + &Bprime_advec[0][0],3,numdof,0, &Ke->values[0],1); } else if(stabilization==2){ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); - tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity); Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 10627) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 10628) @@ -170,4 +170,5 @@ void GetSidList(int* sidlist); int GetElementType(void); + void GetElementSizes(double* hx,double* hy,double* hz); void GetInputListOnVertices(double* pvalue,int enumtype); void GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue); Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 10627) +++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp (revision 10628) @@ -384,29 +384,4 @@ } -} -/*}}}*/ -/*FUNCTION PentaRef::GetBArtdiff {{{1*/ -void PentaRef::GetBArtdiff(double* B_stab, double* xyz_list, GaussPenta* gauss){ - /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. - * For node i, Bi' can be expressed in the actual coordinate system - * by: - * Bi_stab=[ dh/dx ] - * [ dh/dy ] - * where h is the interpolation function for node i. - * - * We assume B has been allocated already, of size: 2x(NDOF1*NUMNODESP1) - */ - - /*Same thing in the actual coordinate system: */ - double dbasis[3][NUMNODESP1]; - - /*Get dbasis in actual coordinates system : */ - GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss); - - /*Build B': */ - for (int i=0;i