Changeset 5844


Ignore:
Timestamp:
09/16/10 11:58:59 (15 years ago)
Author:
seroussi
Message:

clean tria

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5843 r5844  
    24222422        double     DLprime[2][2]                    = {0.0};
    24232423        double     DL_scalar;
    2424         double     Ke_gg[numdof][numdof]            = {0.0};
    24252424        double     Ke_gg_gaussian[numdof][numdof]   = {0.0};
    24262425        double     Ke_gg_thickness1[numdof][numdof] = {0.0};
    … …  
    24282427        GaussTria *gauss                            = NULL;
    24292428
    2430            /*Initialize Element matrix and return if necessary*/
     2429        /*Initialize Element matrix and return if necessary*/
    24312430        if(IsOnWater()) return NULL;
    24322431        ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
    … …  
    24472446        }
    24482447
    2449         //Create Artificial diffusivity once for all if requested
     2448        /*Create Artificial diffusivity once for all if requested*/
    24502449        if(artdiff){
    24512450                gauss=new GaussTria();
    … …  
    24602459        }
    24612460
    2462         /* Start  looping on the number of gaussian points: */
     2461        /*Start looping on the number of gaussian points:*/
    24632462        gauss=new GaussTria(2);
    24642463        for (ig=gauss->begin();ig<gauss->end();ig++){
    … …  
    25192518ElementMatrix* Tria::CreateKMatrixBalancedthickness_DG(void){
    25202519
    2521         /* local declarations */
    2522         const int    numdof=NDOF1*NUMVERTICES;
    2523         int             i,j;
    2524         int     ig;
    2525         double       xyz_list[NUMVERTICES][3];
    2526         GaussTria *gauss=NULL;
    2527         double B[2][NUMVERTICES];
    2528         double Bprime[2][NUMVERTICES];
    2529         double DL[2][2]={0.0};
    2530         double DLprime[2][2]={0.0};
    2531         double DL_scalar;
    2532         double Ke_gg[numdof][numdof]={0.0};
    2533         double Ke_gg2[numdof][numdof]={0.0};
    2534         double Jdettria;
    2535         double  vx,vy;
    2536         int     dim;
     2520        /*Constants*/
     2521        const int  numdof=NDOF1*NUMVERTICES;
     2522
     2523        /*Intermediaries*/
     2524        int        i,j,ig,dim;
     2525        double     vx,vy,Jdettria;
     2526        double     xyz_list[NUMVERTICES][3];
     2527        double     B[2][NUMVERTICES];
     2528        double     Bprime[2][NUMVERTICES];
     2529        double     DL[2][2]={0.0};
     2530        double     DL_scalar;
     2531        double     Ke_gg[numdof][numdof]={0.0};
     2532        GaussTria  *gauss=NULL;
    25372533
    25382534        /*Initialize Element matrix and return if necessary*/
    … …  
    25462542        Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
    25472543
    2548         /* Start  looping on the number of gaussian points: */
     2544        /*Start looping on the number of gaussian points:*/
    25492545        gauss=new GaussTria(2);
    25502546        for (ig=gauss->begin();ig<gauss->end();ig++){
    … …  
    25522548                gauss->GaussPoint(ig);
    25532549
     2550                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    25542551                /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
    25552552                GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    … …  
    25592556                vy_input->GetParameterValue(&vy,gauss);
    25602557
    2561                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    25622558                DL_scalar=-gauss->weight*Jdettria;
    2563 
    2564                 DLprime[0][0]=DL_scalar*vx;
    2565                 DLprime[1][1]=DL_scalar*vy;
     2559                DL[0][0]=DL_scalar*vx;
     2560                DL[1][1]=DL_scalar*vy;
    25662561
    25672562                TripleMultiply( &B[0][0],2,numdof,1,
    2568                                         &DLprime[0][0],2,2,0,
     2563                                        &DL[0][0],2,2,0,
    25692564                                        &Bprime[0][0],2,numdof,0,
    2570                                         &Ke_gg2[0][0],0);
    2571 
    2572                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg2[i][j];
     2565                                        &Ke_gg[0][0],0);
     2566
     2567                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg[i][j];
    25732568        }
    25742569
    … …  
    25812576ElementMatrix* Tria::CreateKMatrixBalancedvelocities(void){
    25822577
    2583         /* local declarations */
     2578        /*Constants*/
    25842579        const int    numdof=NDOF1*NUMVERTICES;
    2585         int             i,j;
    2586         double       xyz_list[NUMVERTICES][3];
    2587         int     ig;
    2588         GaussTria *gauss=NULL;
    2589         double L[NUMVERTICES];
    2590         double B[2][NUMVERTICES];
    2591         double Bprime[2][NUMVERTICES];
    2592         double DL[2][2]={0.0};
    2593         double DLprime[2][2]={0.0};
    2594         double DL_scalar;
    2595         double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix
    2596         double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    2597         double Ke_gg_velocities1[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    2598         double Ke_gg_velocities2[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
    2599         double Jdettria;
     2580
     2581        /*Intermediaries */
     2582        bool    artdiff;
     2583        int     i,j,ig,dim;
     2584        double  nx,ny,norm,Jdettria;
     2585        double  dvx[2],dvy[2];
     2586        double  vx,vy,dvxdx,dvydy;
     2587        double  v_gauss[2]={0.0};
    26002588        double  surface_normal[3];
    26012589        double  surface_list[3];
    2602         double  nx,ny,norm;
    2603         double  dvx[2];
    2604         double  dvy[2];
    2605         double  vx,vy;
    2606         double  dvxdx,dvydy;
    2607         double  v_gauss[2]={0.0};
     2590        double  xyz_list[NUMVERTICES][3];
     2591        double  B[2][NUMVERTICES];
     2592        double  Bprime[2][NUMVERTICES];
    26082593        double  K[2][2]={0.0};
    26092594        double  KDL[2][2]={0.0};
    2610         int     dim;
    2611         bool artdiff;
     2595        double  DLprime[2][2]={0.0};
     2596        double  DL_scalar;
     2597        double  Ke_gg_gaussian[numdof][numdof]   = {0.0};
     2598        double  Ke_gg_velocities[numdof][numdof] = {0.0};
     2599        GaussTria *gauss=NULL;
    26122600
    26132601        /*Initialize Element matrix and return if necessary*/
    … …  
    26442632        }
    26452633        if(nx==0 && ny==0){
    2646                 nx=0;
    2647                 ny=1;
     2634                nx=0; ny=1;
    26482635        }
    26492636        norm=pow( pow(nx,2)+pow(ny,2) , (double).5);
    2650         nx=nx/norm;
    2651         ny=ny/norm;
     2637        nx=nx/norm; ny=ny/norm;
    26522638
    26532639        //Create Artificial diffusivity once for all if requested
    … …  
    26602646                vxaverage_input->GetParameterAverage(&v_gauss[0]);
    26612647                vyaverage_input->GetParameterAverage(&v_gauss[1]);
    2662 
    26632648                K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!!
    26642649                K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]);
    … …  
    26712656                gauss->GaussPoint(ig);
    26722657
     2658                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    26732659                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
    26742660                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
    2675                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    26762661
    26772662                vxaverage_input->GetParameterValue(&vx,gauss);
    … …  
    26792664                vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
    26802665                vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     2666
    26812667                dvxdx=dvx[0];
    26822668                dvydy=dvy[1];
    2683 
    26842669                DL_scalar=gauss->weight*Jdettria;
    26852670
    … …  
    26902675                                        &DLprime[0][0],2,2,0,
    26912676                                        &Bprime[0][0],2,numdof,0,
    2692                                         &Ke_gg_velocities2[0][0],0);
    2693 
    2694                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_velocities2[i][j];
     2677                                        &Ke_gg_velocities[0][0],0);
     2678
     2679                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_velocities[i][j];
    26952680
    26962681                if(artdiff){
    … …  
    27042689
    27052690                        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
    2706 
    27072691                }
    2708 
    27092692        }
    27102693
    … …  
    27382721        int        i,j,ig;
    27392722        double     xyz_list[NUMVERTICES][3];
    2740         GaussTria *gauss = NULL;
    27412723        double     viscosity,newviscosity,oldviscosity;
    2742         double     viscosity_overshoot,thickness;
    2743         double     epsilon[3];    /* epsilon=[exx,eyy,exy];    */
    2744         double     oldepsilon[3]; /* oldepsilon=[exx,eyy,exy]; */
     2724        double     viscosity_overshoot,thickness,Jdet;
     2725        double     epsilon[3],oldepsilon[3];    /* epsilon=[exx,eyy,exy];    */
    27452726        double     B[3][numdof];
    27462727        double     Bprime[3][numdof];
    … …  
    27482729        double     D_scalar;
    27492730        double     Ke_g[numdof][numdof];
    2750         double     Jdet;
     2731        GaussTria *gauss = NULL;
    27512732
    27522733        /*Initialize Element matrix and return if necessary*/
    … …  
    27692750                gauss->GaussPoint(ig);
    27702751
    2771                 thickness_input->GetParameterValue(&thickness, gauss);
     2752                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2753                GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
     2754                GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
     2755
    27722756                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    27732757                this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    27742758                matice->GetViscosity2d(&viscosity, &epsilon[0]);
    27752759                matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
    2776                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     2760                thickness_input->GetParameterValue(&thickness, gauss);
    27772761
    27782762                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    27792763                D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
    27802764                for (i=0;i<3;i++) D[i][i]=D_scalar;
    2781 
    2782                 GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
    2783                 GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
    27842765
    27852766                TripleMultiply(&B[0][0],3,numdof,1,
Note: See TracChangeset for help on using the changeset viewer.