Changeset 5722


Ignore:
Timestamp:
09/09/10 14:17:56 (15 years ago)
Author:
Mathieu Morlighem
Message:

Deleted some functions and added more GaussPenta

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

Legend:

Unmodified
Added
Removed

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

    r5719 r5722  
    508508        int  dofp[1]={3};
    509509        double Jdet2d;
    510         Tria* tria=NULL;
    511510
    512511        /*Gauss*/
    513         int     num_gauss,ig;
    514         double* first_gauss_area_coord  =  NULL;
    515         double* second_gauss_area_coord =  NULL;
    516         double* third_gauss_area_coord  =  NULL;
    517         double* gauss_weights           =  NULL;
    518         double  gauss_weight;
    519         double  gauss_coord[4];
     512        int     ig;
     513        GaussPenta* gauss;
    520514
    521515        /*Output*/
    … …  
    529523        bool onbed;
    530524        int  approximation;
    531         Input* pressure_input=NULL;
    532         Input* vx_input=NULL;
    533         Input* vy_input=NULL;
    534         Input* vz_input=NULL;
    535525
    536526        /*parameters: */
    … …  
    569559        }
    570560
    571         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    572         GaussLegendreTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2);
    573 
    574561        /*Retrieve all inputs we will be needing: */
    575         pressure_input=inputs->GetInput(PressureEnum);
    576         vx_input=inputs->GetInput(VxEnum);
    577         vy_input=inputs->GetInput(VyEnum);
    578         vz_input=inputs->GetInput(VzEnum);
     562        Input* pressure_input=inputs->GetInput(PressureEnum);
     563        Input* vx_input=inputs->GetInput(VxEnum);
     564        Input* vy_input=inputs->GetInput(VyEnum);
     565        Input* vz_input=inputs->GetInput(VzEnum);
    579566
    580567        /* Start  looping on the number of gaussian points: */
     568        gauss=new GaussPenta(0,1,2,2);
    581569        for (ig=0; ig<num_gauss; ig++){
    582570
    583                         /*Pick up the gaussian point: */
    584                         gauss_weight=*(gauss_weights+ig);
    585                         gauss_coord[0]=*(first_gauss_area_coord+ig);
    586                         gauss_coord[1]=*(second_gauss_area_coord+ig);
    587                         gauss_coord[2]=*(third_gauss_area_coord+ig);
    588                         gauss_coord[3]=-1.0; //take the base
    589 
    590                         /*Compute strain rate viscosity and pressure: */
    591                         this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
    592                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    593                         pressure_input->GetParameterValue(&pressure, &gauss_coord[0]);
    594 
    595                         /*Compute Stress*/
    596                         sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure
    597                         sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning;
    598                         sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning;
    599                         sigma_xy=2*viscosity*epsilon[3];
    600                         sigma_xz=2*viscosity*epsilon[4];
    601                         sigma_yz=2*viscosity*epsilon[5];
    602 
    603                         /*Get normal vector to the bed */
    604                         BedNormal(&bed_normal[0],xyz_list_tria);
    605 
    606                         /*basalforce*/
    607                         basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2];
    608                         basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2];
    609                         basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2];
    610 
    611                         /*Get the Jacobian determinant */
    612                         tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0],gauss_coord);
    613                         value+=sigma_zz*Jdet2d*gauss_weight;
    614                         surface+=Jdet2d*gauss_weight;
     571                gauss->GaussPoint(ig);
     572
     573                /*Compute strain rate viscosity and pressure: */
     574                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     575                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     576                pressure_input->GetParameterValue(&pressure, &gauss_coord[0]);
     577
     578                /*Compute Stress*/
     579                sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure
     580                sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning;
     581                sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning;
     582                sigma_xy=2*viscosity*epsilon[3];
     583                sigma_xz=2*viscosity*epsilon[4];
     584                sigma_yz=2*viscosity*epsilon[5];
     585
     586                /*Get normal vector to the bed */
     587                BedNormal(&bed_normal[0],xyz_list_tria);
     588
     589                /*basalforce*/
     590                basalforce[0] += sigma_xx*bed_normal[0] + sigma_xy*bed_normal[1] + sigma_xz*bed_normal[2];
     591                basalforce[1] += sigma_xy*bed_normal[0] + sigma_yy*bed_normal[1] + sigma_yz*bed_normal[2];
     592                basalforce[2] += sigma_xz*bed_normal[0] + sigma_yz*bed_normal[1] + sigma_zz*bed_normal[2];
     593
     594                /*Get the Jacobian determinant */
     595                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
     596                value+=sigma_zz*Jdet2d*gauss_weight;
     597                surface+=Jdet2d*gauss_weight;
    615598        }
    616599        value=value/surface;
    … …  
    27352718        const int numdof2d=NUMVERTICES2D*NDOF4;
    27362719
    2737         /*Collapsed formulation: */
    2738         Tria*  tria=NULL;
    2739 
    27402720        /*Grid data: */
    27412721        double     xyz_list[NUMVERTICES][3];
    … …  
    27612741        int     ig;
    27622742        GaussPenta *gauss=NULL;
    2763         GaussTria  *gauss_tria=NULL;
    2764 
    27652743
    27662744        /* specific gaussian point: */
    … …  
    27862764        bool shelf;
    27872765
    2788         /*inputs: */
    2789         Input* vx_input=NULL;
    2790         Input* vy_input=NULL;
    2791         Input* vz_input=NULL;
    2792 
    27932766        /*retrive parameters: */
    27942767        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    28082781        GetDofList(&doflist,StokesApproximationEnum);
    28092782
    2810         /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    2811                 get tria gaussian points as well as segment gaussian points. For tria gaussian
    2812                 points, order of integration is 2, because we need to integrate the product tB*D*B'
    2813                 which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian
    2814                 points, same deal, which yields 3 gaussian points.*/
    2815 
    28162783        /*Retrieve all inputs we will be needing: */
    2817         vx_input=inputs->GetInput(VxEnum);
    2818         vy_input=inputs->GetInput(VyEnum);
    2819         vz_input=inputs->GetInput(VzEnum);
     2784        Input* vx_input=inputs->GetInput(VxEnum);
     2785        Input* vy_input=inputs->GetInput(VyEnum);
     2786        Input* vz_input=inputs->GetInput(VzEnum);
    28202787
    28212788        /* Start  looping on the number of gaussian points: */
    … …  
    28532820                for(i=0;i<27;i++) for(j=0;j<27;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
    28542821        }
     2822        delete gauss; //gauss of previous loop
    28552823
    28562824        if((onbed==1) && (shelf==0)){
    … …  
    28662834
    28672835                /* Start  looping on the number of gaussian points: */
    2868                 delete gauss; //gauss of previous loop
    2869                 gauss=new GaussPenta();
    2870                 gauss->GaussFaceTria(0,1,2,2);
    2871                 gauss_tria=new GaussTria();
     2836                gauss=new GaussPenta(0,1,2,2);
    28722837                for (ig=gauss->begin();ig<gauss->end();ig++){
    28732838
    28742839                        gauss->GaussPoint(ig);
    2875                         gauss->SynchronizeGaussTria(gauss_tria);
    28762840
    28772841                        /*Get the Jacobian determinant */
    2878                         tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria);
     2842                        GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    28792843
    28802844                        /*Get L matrix if viscous basal drag present: */
    … …  
    29202884                /*Free ressources:*/
    29212885                delete friction;
     2886                delete gauss;
    29222887
    29232888        } //if ( (onbed==1) && (shelf==0))
    … …  
    29302895
    29312896        /*Free ressources:*/
    2932         delete gauss;
    2933         delete gauss_tria;
    29342897        xfree((void**)&doflist);
    2935 
    2936         return;
    29372898}
    29382899/*}}}*/
    … …  
    37423703        int     ig;
    37433704        GaussPenta *gauss=NULL;
    3744         GaussTria  *gauss_tria=NULL;
    37453705
    37463706        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    … …  
    38543814                for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j];
    38553815        }
     3816        delete gauss;
    38563817
    38573818        /*Deal with 2d friction at the bedrock interface: */
    … …  
    38653826
    38663827                /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    3867                 delete gauss; //gauss of previous loop
    3868                 gauss=new GaussPenta();
    3869                 gauss->GaussFaceTria(0,1,2,2);
    3870                 gauss_tria=new GaussTria();
     3828                gauss=new GaussPenta(0,1,2,2)
    38713829                for (ig=gauss->begin();ig<gauss->end();ig++){
    38723830
    38733831                        gauss->GaussPoint(ig);
    3874                         gauss->SynchronizeGaussTria(gauss_tria);
    38753832
    38763833                        /*Get the Jacobian determinant */
    3877                         tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria);
     3834                        GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    38783835
    38793836                        /* Get bed at gaussian point */
    … …  
    38953852                        }
    38963853                }
     3854                /*Free ressources:*/
     3855                delete gauss;
    38973856        } //if ( (onbed==1) && (shelf==1))
    38983857
    … …  
    39053864        /*Free ressources:*/
    39063865        delete gauss;
    3907         delete gauss_tria;
    39083866        xfree((void**)&doflist);
    39093867

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

    r5677 r5722  
    16771677}
    16781678/*}}}*/
     1679/*FUNCTION PentaRef::GetTriaJacobianDeterminant{{{1*/
     1680void PentaRef::GetTriaJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss){
     1681        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     1682         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     1683
     1684        double x1,x2,x3,y1,y2,y3,z1,z2,z3;
     1685
     1686        x1=*(xyz_list+3*0+0);
     1687        y1=*(xyz_list+3*0+1);
     1688        z1=*(xyz_list+3*0+2);
     1689        x2=*(xyz_list+3*1+0);
     1690        y2=*(xyz_list+3*1+1);
     1691        z2=*(xyz_list+3*1+2);
     1692        x3=*(xyz_list+3*2+0);
     1693        y3=*(xyz_list+3*2+1);
     1694        z3=*(xyz_list+3*2+2);
     1695
     1696        *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
     1697        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
     1698
     1699}
     1700/*}}}*/
    16791701/*FUNCTION PentaRef::GetJacobianInvert {{{1*/
    16801702void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){

  • issm/trunk/src/c/objects/Elements/PentaRef.h

    r5677 r5722  
    5858                void GetJacobian(double* J, double* xyz_list,GaussPenta* gauss);
    5959                void GetJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
     60                void GetTriaJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
    6061                void GetJacobianInvert(double*  Jinv, double* xyz_list,GaussPenta* gauss);
    6162                void GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss);

  • issm/trunk/src/c/objects/Gauss/GaussPenta.cpp

    r5675 r5722  
    8080}
    8181/*}}}*/
     82/*FUNCTION GaussPenta::GaussPenta{{{1*/
     83GaussPenta::GaussPenta(int index1, int index2, int index3, int order){
     84
     85        /*Basal Tria*/
     86        if(index1==0 && index2==1 && index3==2){
     87
     88                /*Get GaussTria*/
     89                GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);
     90
     91                /*compute z coordinate*/
     92                coords4=(double*)xmalloc(numgauss*sizeof(double));
     93                for(int i=0;i<numgauss;i++) coords4[i]=-1.0;
     94        }
     95        else{
     96                ISSMERROR("Tria not supported yet");
     97        }
     98
     99}
     100/*}}}*/
    82101/*FUNCTION GaussPenta::~GaussPenta(){{{1*/
    83102GaussPenta::~GaussPenta(){

  • issm/trunk/src/c/objects/Gauss/GaussPenta.h

    r5675 r5722  
    3434                GaussPenta();
    3535                GaussPenta(int order_horiz,int order_vert);
     36                GaussPenta(int index1, int index2, int index3, int order);
    3637                ~GaussPenta();
    3738

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5720 r5722  
    13541354}
    13551355/*}}}*/
    1356 /*FUNCTION Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/
    1357 void Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){
    1358 
    1359         /*output: */
    1360         double value;
    1361 
    1362         /*Get value on Element 1*/
    1363         element->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in);
    1364 
    1365         /*Assign output pointers:*/
    1366         *pvalue=value;
    1367 }
    1368 /*}}}*/
    13691356/*FUNCTION Icefront::GetNormal {{{1*/
    13701357void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){

  • issm/trunk/src/c/objects/Loads/Icefront.h

    r5715 r5722  
    8585                void  QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    8686                                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
    87                 void GetParameterValue(double* pvalue, double gauss_coord,Input* input_in);
    8887                void GetNormal(double* normal,double xyz_list[2][3]);
    8988                /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.