Changeset 10628
- Timestamp:
- 11/14/11 10:59:49 (13 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r10597 r10628 822 822 /*return PentaRef field*/ 823 823 return this->element_type; 824 } 825 /*}}}*/ 826 /*FUNCTION Penta::GetElementSizes{{{1*/ 827 void Penta::GetElementSizes(double* hx,double* hy,double* hz){ 828 829 double xyz_list[NUMVERTICES][3]; 830 double xmin,ymin,zmin; 831 double xmax,ymax,zmax; 832 833 /*Get xyz list: */ 834 GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES); 835 xmin=xyz_list[0][0]; xmax=xyz_list[0][0]; 836 ymin=xyz_list[0][1]; ymax=xyz_list[0][1]; 837 zmin=xyz_list[0][2]; zmax=xyz_list[0][2]; 838 839 for(int i=1;i<NUMVERTICES;i++){ 840 if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0]; 841 if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0]; 842 if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1]; 843 if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1]; 844 if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2]; 845 if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2]; 846 } 847 848 *hx=xmax-xmin; 849 *hy=ymax-ymin; 850 *hz=zmax-zmin; 824 851 } 825 852 /*}}}*/ … … 3180 3207 K[0][0]=D_scalar_stab*pow(u,2); K[0][1]=D_scalar_stab*fabs(u)*fabs(v); 3181 3208 K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2); 3182 3183 GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); 3209 _error_("TO BE RECODED"); 3210 3211 //GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); 3184 3212 3185 3213 TripleMultiply(&B_stab[0][0],2,numdof,1, … … 3304 3332 int stabilization; 3305 3333 int i,j,ig,found=0; 3306 double Jdet,u,v,w,um,vm,wm ;3307 double epsvel=2.220446049250313e-16;3334 double Jdet,u,v,w,um,vm,wm,vel; 3335 double h,hx,hy,hz,vx,vy,vz; 3308 3336 double gravity,rho_ice,rho_water; 3309 3337 double heatcapacity,thermalconductivity,dt; … … 3321 3349 double D_scalar_trans,D_scalar_stab; 3322 3350 double D[3][3]; 3323 double K[ 2][2]={0.0};3351 double K[3][3]={0.0}; 3324 3352 Tria* tria=NULL; 3325 3353 GaussPenta *gauss=NULL; … … 3374 3402 GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 3375 3403 3376 vx_input->GetInputValue(&u, gauss); 3377 vy_input->GetInputValue(&v, gauss); 3378 vz_input->GetInputValue(&w, gauss); 3379 vxm_input->GetInputValue(&um,gauss); 3380 vym_input->GetInputValue(&vm,gauss); 3381 vzm_input->GetInputValue(&wm,gauss); 3404 vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um; 3405 vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm; 3406 vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm; 3382 3407 3383 3408 D_scalar_advec=gauss->weight*Jdet; 3384 3409 if(dt) D_scalar_advec=D_scalar_advec*dt; 3385 3410 3386 D[0][0]=D_scalar_advec* (u-um);D[0][1]=0; D[0][2]=0;3387 D[1][0]=0; D[1][1]=D_scalar_advec* (v-vm);D[1][2]=0;3388 D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_advec* (w-wm);3411 D[0][0]=D_scalar_advec*vx; D[0][1]=0; D[0][2]=0; 3412 D[1][0]=0; D[1][1]=D_scalar_advec*vy; D[1][2]=0; 3413 D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_advec*vz; 3389 3414 3390 3415 TripleMultiply(&B_advec[0][0],3,numdof,1, … … 3410 3435 if(stabilization==1){ 3411 3436 /*Build K: */ 3412 D_scalar_stab=gauss->weight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel); 3437 GetElementSizes(&hx,&hy,&hz); 3438 vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14; 3439 h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.)); 3440 3441 K[0][0]=h/(2*vel)*vx*vx; K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz; 3442 K[1][0]=h/(2*vel)*vy*vx; K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz; 3443 K[2][0]=h/(2*vel)*vz*vx; K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz; 3444 3445 D_scalar_stab=gauss->weight*Jdet; 3413 3446 if(dt) D_scalar_stab=D_scalar_stab*dt; 3414 K[0][0]=D_scalar_stab*pow(u,2); K[0][1]=D_scalar_stab*fabs(u)*fabs(v); 3415 K[1][0]=D_scalar_stab*fabs(u)*fabs(v);K[1][1]=D_scalar_stab*pow(v,2); 3416 3417 GetBArtdiff(&B_stab[0][0],&xyz_list[0][0],gauss); 3418 3419 TripleMultiply(&B_stab[0][0],2,numdof,1, 3420 &K[0][0],2,2,0, 3421 &B_stab[0][0],2,numdof,0, 3447 for(i=0;i<3;i++) for(j=0;j<3;j++) K[i][j] = D_scalar_stab*K[i][j]; 3448 3449 GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 3450 3451 TripleMultiply(&Bprime_advec[0][0],3,numdof,1, 3452 &K[0][0],3,3,0, 3453 &Bprime_advec[0][0],3,numdof,0, 3422 3454 &Ke->values[0],1); 3423 3455 } 3424 3456 else if(stabilization==2){ 3425 3426 3457 GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss); 3427 3428 3458 tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity); 3429 3459 -
issm/trunk/src/c/objects/Elements/Penta.h
r10571 r10628 170 170 void GetSidList(int* sidlist); 171 171 int GetElementType(void); 172 void GetElementSizes(double* hx,double* hy,double* hz); 172 173 void GetInputListOnVertices(double* pvalue,int enumtype); 173 174 void GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue); -
issm/trunk/src/c/objects/Elements/PentaRef.cpp
r10514 r10628 384 384 } 385 385 386 }387 /*}}}*/388 /*FUNCTION PentaRef::GetBArtdiff {{{1*/389 void PentaRef::GetBArtdiff(double* B_stab, double* xyz_list, GaussPenta* gauss){390 /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.391 * For node i, Bi' can be expressed in the actual coordinate system392 * by:393 * Bi_stab=[ dh/dx ]394 * [ dh/dy ]395 * where h is the interpolation function for node i.396 *397 * We assume B has been allocated already, of size: 2x(NDOF1*NUMNODESP1)398 */399 400 /*Same thing in the actual coordinate system: */401 double dbasis[3][NUMNODESP1];402 403 /*Get dbasis in actual coordinates system : */404 GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);405 406 /*Build B': */407 for (int i=0;i<NUMNODESP1;i++){408 *(B_stab+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i];409 *(B_stab+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i];410 }411 386 } 412 387 /*}}}*/ -
issm/trunk/src/c/objects/Elements/PentaRef.h
r10514 r10628 45 45 void GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss); 46 46 void GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss); 47 void GetBArtdiff(double* B_stab, double* xyz_list, GaussPenta* gauss);48 47 void GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss); 49 48 void GetBVert(double* B, double* xyz_list, GaussPenta* gauss);
Note:
See TracChangeset
for help on using the changeset viewer.