Changeset 3160


Ignore:
Timestamp:
03/02/10 14:34:14 (15 years ago)
Author:
Mathieu Morlighem
Message:

minor

Location:
issm/trunk/src/c/objects
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Penta.cpp

    r3159 r3160  
    30493049        double z1,z2,z3,z4,z5,z6;
    30503050
    3051         double sqrt3=SQRT3;
    3052 
    30533051        /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
    30543052        A1=gauss_coord[0];
    … …  
    30573055
    30583056        xi=A2-A1;
    3059         eta=sqrt3*A3;
     3057        eta=SQRT3*A3;
    30603058        zi=gauss_coord[3];
    30613059
    … …  
    30823080
    30833081
    3084         *(J+NDOF3*0+0)=1.0/4.0*(x1-x2-x4+x5)*zi+1.0/4.0*(-x1+x2-x4+x5);
    3085         *(J+NDOF3*1+0)=sqrt3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+sqrt3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
    3086         *(J+NDOF3*2+0)=sqrt3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +1.0/4.0*(-x1+x5-x2+x4);
    3087 
    3088         *(J+NDOF3*0+1)=1.0/4.0*(y1-y2-y4+y5)*zi+1.0/4.0*(-y1+y2-y4+y5);
    3089         *(J+NDOF3*1+1)=sqrt3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+sqrt3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
    3090         *(J+NDOF3*2+1)=sqrt3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+1.0/4.0*(y1-y2-y4+y5)*xi+1.0/4.0*(y4-y1+y5-y2);
    3091 
    3092         *(J+NDOF3*0+2)=1.0/4.0*(z1-z2-z4+z5)*zi+1.0/4.0*(-z1+z2-z4+z5);
    3093         *(J+NDOF3*1+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+sqrt3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
    3094         *(J+NDOF3*2+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+1.0/4.0*(z1-z2-z4+z5)*xi+1.0/4.0*(-z1+z5-z2+z4);
     3082        *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
     3083        *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
     3084        *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
     3085
     3086        *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
     3087        *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
     3088        *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
     3089
     3090        *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
     3091        *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
     3092        *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
    30953093
    30963094#ifdef _ISSM_DEBUG_
    … …  
    34613459        const int numgrids=6;
    34623460        double A1,A2,A3,z;
    3463         double sqrt3=SQRT3;
    34643461
    34653462        A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
    … …  
    34703467
    34713468        /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
    3472         *(dl1dl6+numgrids*0+0)=-1.0/2.0*(1.0-z)/2.0;
    3473         *(dl1dl6+numgrids*1+0)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
    3474         *(dl1dl6+numgrids*2+0)=-1.0/2.0*A1;
     3469        *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;
     3470        *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
     3471        *(dl1dl6+numgrids*2+0)=-0.5*A1;
    34753472
    34763473        /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
    3477         *(dl1dl6+numgrids*0+1)=1.0/2.0*(1.0-z)/2.0;
    3478         *(dl1dl6+numgrids*1+1)=-1.0/2.0/sqrt3*(1.0-z)/2.0;
    3479         *(dl1dl6+numgrids*2+1)=-1.0/2.0*A2;
     3474        *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;
     3475        *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
     3476        *(dl1dl6+numgrids*2+1)=-0.5*A2;
    34803477
    34813478        /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
    34823479        *(dl1dl6+numgrids*0+2)=0.0;
    3483         *(dl1dl6+numgrids*1+2)=1.0/sqrt3*(1.0-z)/2.0;
    3484         *(dl1dl6+numgrids*2+2)=-1.0/2.0*A3;
     3480        *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
     3481        *(dl1dl6+numgrids*2+2)=-0.5*A3;
    34853482
    34863483        /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
    3487         *(dl1dl6+numgrids*0+3)=-1.0/2.0*(1.0+z)/2.0;
    3488         *(dl1dl6+numgrids*1+3)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
    3489         *(dl1dl6+numgrids*2+3)=1.0/2.0*A1;
     3484        *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;
     3485        *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
     3486        *(dl1dl6+numgrids*2+3)=0.5*A1;
    34903487
    34913488        /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
    3492         *(dl1dl6+numgrids*0+4)=1.0/2.0*(1.0+z)/2.0;
    3493         *(dl1dl6+numgrids*1+4)=-1.0/2.0/sqrt3*(1.0+z)/2.0;
    3494         *(dl1dl6+numgrids*2+4)=1.0/2.0*A2;
     3489        *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;
     3490        *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
     3491        *(dl1dl6+numgrids*2+4)=0.5*A2;
    34953492
    34963493        /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
    34973494        *(dl1dl6+numgrids*0+5)=0.0;
    3498         *(dl1dl6+numgrids*1+5)=1.0/sqrt3*(1.0+z)/2.0;
    3499         *(dl1dl6+numgrids*2+5)=1.0/2.0*A3;
     3495        *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
     3496        *(dl1dl6+numgrids*2+5)=0.5*A3;
    35003497}
    35013498/*}}}*/
    … …  
    35083505         * natural coordinate system) at the gaussian point. */
    35093506
    3510         double sqrt3=SQRT3;
    35113507        int    numgrids=7; //six plus bubble grids
    35123508
    35133509        double r=gauss_coord[1]-gauss_coord[0];
    3514         double s=-3.0/sqrt3*(gauss_coord[0]+gauss_coord[1]-2.0/3.0);
     3510        double s=-3.0/SQRT3*(gauss_coord[0]+gauss_coord[1]-2.0/3.0);
    35153511        double zeta=gauss_coord[3];
    35163512
    35173513        /*First nodal function: */
    3518         *(dl1dl7+numgrids*0+0)=-1.0/2.0*(1.0-zeta)/2.0;
    3519         *(dl1dl7+numgrids*1+0)=-sqrt3/6.0*(1.0-zeta)/2.0;
    3520         *(dl1dl7+numgrids*2+0)=-1.0/2.0*(-1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
     3514        *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;
     3515        *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
     3516        *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
    35213517
    35223518        /*Second nodal function: */
    3523         *(dl1dl7+numgrids*0+1)=1.0/2.0*(1.0-zeta)/2.0;
    3524         *(dl1dl7+numgrids*1+1)=-sqrt3/6.0*(1.0-zeta)/2.0;
    3525         *(dl1dl7+numgrids*2+1)=-1.0/2.0*(1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
     3519        *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;
     3520        *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
     3521        *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
    35263522
    35273523        /*Third nodal function: */
    35283524        *(dl1dl7+numgrids*0+2)=0;
    3529         *(dl1dl7+numgrids*1+2)=sqrt3/3.0*(1.0-zeta)/2.0;
    3530         *(dl1dl7+numgrids*2+2)=-1.0/2.0*(sqrt3/3.0*s+1.0/3.0);
     3525        *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
     3526        *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
    35313527
    35323528        /*Fourth nodal function: */
    3533         *(dl1dl7+numgrids*0+3)=-1.0/2.0*(1.0+zeta)/2.0;
    3534         *(dl1dl7+numgrids*1+3)=-sqrt3/6.0*(1.0+zeta)/2.0;
    3535         *(dl1dl7+numgrids*2+3)=1.0/2.0*(-1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
     3529        *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;
     3530        *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
     3531        *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
    35363532
    35373533        /*Fith nodal function: */
    3538         *(dl1dl7+numgrids*0+4)=1.0/2.0*(1.0+zeta)/2.0;
    3539         *(dl1dl7+numgrids*1+4)=-sqrt3/6.0*(1.0+zeta)/2.0;
    3540         *(dl1dl7+numgrids*2+4)=1.0/2.0*(1.0/2.0*r-sqrt3/6.0*s+1.0/3.0);
     3534        *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;
     3535        *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
     3536        *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
    35413537
    35423538        /*Sixth nodal function: */
    35433539        *(dl1dl7+numgrids*0+5)=0;
    3544         *(dl1dl7+numgrids*1+5)=sqrt3/3.0*(1.0+zeta)/2.0;
    3545         *(dl1dl7+numgrids*2+5)=1.0/2.0*(sqrt3/3.0*s+1.0/3.0);
     3540        *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
     3541        *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
    35463542
    35473543        /*Seventh nodal function: */
    3548         *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(sqrt3*s+1.0);
    3549         *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(sqrt3*pow(s,2.0)-2.0*s-sqrt3*pow(r,2.0));
     3544        *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
     3545        *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
    35503546        *(dl1dl7+numgrids*2+6)=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(-2.0*zeta);
    35513547
    … …  
    37113707        epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2);
    37123708
    3713         epsilon_eff=1/pow(2,1.0/2.0)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),1.0/2.0);
     3709        epsilon_eff=1/pow(2,0.5)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),0.5);
    37143710        *phi=2*pow(epsilon_eff,2.0)*viscosity;
    37153711}

  • issm/trunk/src/c/objects/Tria.cpp

    r3159 r3160  
    388388        if(numpar->artdiff){
    389389                //Get the Jacobian determinant
    390                 gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
     390                gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
    391391                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
    392392
    393393                //Build K matrix (artificial diffusivity matrix)
    394                 v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
    395                 v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
     394                v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
     395                v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
    396396
    397397                K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
    … …  
    611611        if(numpar->artdiff){
    612612                //Get the Jacobian determinant
    613                 gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
     613                gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
    614614                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
    615615
    616616                //Build K matrix (artificial diffusivity matrix)
    617                 v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
    618                 v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
     617                v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
     618                v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
    619619
    620620                K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!!
    … …  
    14081408        if(numpar->artdiff==1){
    14091409                //Get the Jacobian determinant
    1410                 gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
     1410                gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD;
    14111411                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
    14121412
    14131413                //Build K matrix (artificial diffusivity matrix)
    1414                 v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
    1415                 v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
     1414                v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
     1415                v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
    14161416
    14171417                K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]);
    … …  
    32623262
    32633263
    3264         *(J+NDOF2*0+0)=1.0/2.0*(x2-x1);
     3264        *(J+NDOF2*0+0)=0.5*(x2-x1);
    32653265        *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
    3266         *(J+NDOF2*0+1)=1.0/2.0*(y2-y1);
     3266        *(J+NDOF2*0+1)=0.5*(y2-y1);
    32673267        *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
    32683268}
    … …  
    34683468
    34693469        /*First nodal function: */
    3470         *(dl1dl3+numgrids*0+0)=-1.0/2.0;
     3470        *(dl1dl3+numgrids*0+0)=-0.5;
    34713471        *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
    34723472
    34733473        /*Second nodal function: */
    3474         *(dl1dl3+numgrids*0+1)=1.0/2.0;
     3474        *(dl1dl3+numgrids*0+1)=0.5;
    34753475        *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
    34763476
    … …  
    42914291
    42924292        mass_flux= rho_ice*length*( 
    4293                                 (1.0/3.0*(h1-h2)*(vx1-vx2)+1.0/2.0*h2*(vx1-vx2)+1.0/2.0*(h1-h2)*vx2+h2*vx2)*normal[0]+
    4294                                 (1.0/3.0*(h1-h2)*(vy1-vy2)+1.0/2.0*h2*(vy1-vy2)+1.0/2.0*(h1-h2)*vy2+h2*vy2)*normal[1]
     4293                                (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
     4294                                (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
    42954295                                );
    42964296        return mass_flux;
Note: See TracChangeset for help on using the changeset viewer.