Changeset 5726


Ignore:
Timestamp:
09/09/10 15:02:52 (15 years ago)
Author:
Mathieu Morlighem
Message:

Use GaussPenta for Hutter 3d

Location:
issm/trunk/src/c/objects/Elements
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r5724 r5726  
    27402740        int     ig;
    27412741        GaussPenta *gauss=NULL;
    2742 
    2743         /* specific gaussian point: */
    2744         double  gauss_weight,area_gauss_weight,vert_gauss_weight;
    2745         double  gauss_coord[4];
    2746         double  gauss_coord_tria[3];
    2747         int     area_order=5;
    2748         int       num_vert_gauss=5;
    2749 
    27502742        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    27512743        double  viscosity;
     
    34163408        double    pe_g[numdofs]={0.0};
    34173409        double    xyz_list[NUMVERTICES][3];
     3410        double    xyz_list_segment[2][3];
    34183411        double    z_list[NUMVERTICES];
    34193412        double    z_segment[2];
    34203413        double    slope2,constant_part;
    3421         double    gauss[numdofs][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
    34223414        int       node0,node1;
    3423         BeamRef*  beam=NULL;
    3424 
    3425         /*gaussian points: */
    3426         int      num_gauss;
    3427         double*  segment_gauss_coord=NULL;
    3428         double   gauss_coord;
    3429         double*  gauss_weights=NULL;
    3430         double   gauss_weight;
     3415        GaussPenta* gauss=NULL;
    34313416        int      ig;
    34323417        double   slope[2];
     
    34443429        int  connectivity[2];
    34453430
    3446         /*Inputs*/
    3447         Input* thickness_input;
    3448         Input* surface_input;
    3449         Input* slopex_input;
    3450         Input* slopey_input;
    3451 
    34523431        /*recover some inputs: */
    34533432        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
     
    34693448        B=matice->GetB();
    34703449
    3471         //Initilaize beam
    3472         beam=new BeamRef();
    3473 
    34743450        //recover extra inputs
    3475         thickness_input=inputs->GetInput(ThicknessEnum);
    3476         surface_input=inputs->GetInput(SurfaceEnum);
    3477         slopex_input=inputs->GetInput(SurfaceSlopeXEnum);
    3478         slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
     3451        Input* thickness_input=inputs->GetInput(ThicknessEnum);
     3452        Input* surface_input=inputs->GetInput(SurfaceEnum);
     3453        Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum);
     3454        Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum);
    34793455
    34803456        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    34813457        for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2];
    34823458
    3483         //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
    3484         num_gauss=3;
    3485         gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
    3486 
    34873459        /*Loop on the three segments*/
    34883460        for(i=0;i<3;i++){
    3489                 slopex_input->GetParameterValue(&slope[0], &gauss[i][0]);
    3490                 slopey_input->GetParameterValue(&slope[1], &gauss[i][0]);
    3491                 surface_input->GetParameterValue(&surface, &gauss[i][0]);
    3492                 thickness_input->GetParameterValue(&thickness, &gauss[i][0]);
    3493 
    3494                 //compute slope2 slopex and slopey
    3495                 slope2=pow(slope[0],2)+pow(slope[1],2);
    3496 
    3497                 //%compute constant_part
    3498                 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
    3499 
    3500                 z_segment[0]=z_list[i];
    3501                 z_segment[1]=z_list[i+3];
    35023461
    35033462                node0=i;
    35043463                node1=i+3;
    35053464
     3465                for(j=0;j<3;j++){
     3466                        xyz_list_segment[0][j]=xyz_list[node0][j];
     3467                        xyz_list_segment[1][j]=xyz_list[node1][j];
     3468                }
     3469
    35063470                connectivity[0]=nodes[node0]->GetConnectivity();
    35073471                connectivity[1]=nodes[node1]->GetConnectivity();
    35083472
    35093473                /*Loop on the Gauss points: */
    3510                 for(ig=0;ig<num_gauss;ig++){
    3511 
    3512                         //Pick up the gaussian point and its weight:
    3513                         gauss_weight=gauss_weights[ig];
    3514                         gauss_coord=segment_gauss_coord[ig];
    3515 
    3516                         beam->GetParameterValue(&z_g, &z_segment[0],gauss_coord);
    3517 
    3518                         //Get Jacobian determinant:
    3519                         beam->GetJacobianDeterminant(&Jdet, &z_segment[0],gauss_coord);
     3474                gauss=new GaussPenta(node0,node1,3);
     3475                for(ig=gauss->begin();ig<gauss->end();ig++){
     3476                        gauss->GaussPoint(ig);
     3477
     3478                        slopex_input->GetParameterValue(&slope[0],gauss);
     3479                        slopey_input->GetParameterValue(&slope[1],gauss);
     3480                        surface_input->GetParameterValue(&surface,gauss);
     3481                        thickness_input->GetParameterValue(&thickness,gauss);
     3482
     3483                        slope2=pow(slope[0],2)+pow(slope[1],2);
     3484                        constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
     3485
     3486                        PentaRef::GetParameterValue(&z_g,&z_list[0],gauss);
     3487                        GetSegmentJacobianDeterminant(&Jdet,&xyz_list_segment[0][0],gauss);
    35203488
    35213489                        if (onsurface){
    3522                                 for(j=0;j<NDOF2;j++){
    3523                                         pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight/(double)connectivity[1];
    3524                                 }
     3490                                for(j=0;j<NDOF2;j++) pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight/(double)connectivity[1];
    35253491                        }
    35263492                        else{//connectivity is too large, should take only half on it
    3527                                 for(j=0;j<NDOF2;j++){
    3528                                         pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight*2/(double)connectivity[1];
    3529                                 }
     3493                                for(j=0;j<NDOF2;j++) pe_g[2*node1+j]+=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss->weight*2/(double)connectivity[1];
    35303494                        }
    35313495                }
     3496                delete gauss;
    35323497
    35333498                //Deal with lower surface
    35343499                if (onbed){
    3535 
    3536                         //compute ub
    35373500                        constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness;
    35383501                        ub=constant_part*slope[0];
    35393502                        vb=constant_part*slope[1];
    35403503
    3541                         //Add to pe:
    35423504                        pe_g[2*node0]+=ub/(double)connectivity[0];
    35433505                        pe_g[2*node0+1]+=vb/(double)connectivity[0];
     
    35493511
    35503512        /*Clean up */
    3551         delete beam;
    3552         xfree((void**)&gauss_weights);
    3553         xfree((void**)&segment_gauss_coord);
    35543513        xfree((void**)&doflist);
    35553514}
  • issm/trunk/src/c/objects/Elements/PentaRef.cpp

    r5722 r5726  
    16991699}
    17001700/*}}}*/
     1701/*FUNCTION PentaRef::GetSegmentJacobianDeterminant{{{1*/
     1702void PentaRef::GetSegmentJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss){
     1703        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     1704         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
     1705
     1706        double x1,x2,y1,y2,z1,z2;
     1707
     1708        x1=*(xyz_list+3*0+0);
     1709        y1=*(xyz_list+3*0+1);
     1710        z1=*(xyz_list+3*0+2);
     1711        x2=*(xyz_list+3*1+0);
     1712        y2=*(xyz_list+3*1+1);
     1713        z2=*(xyz_list+3*1+2);
     1714
     1715        *Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.) + pow(z2-z1,2.));
     1716        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
     1717
     1718}
     1719/*}}}*/
    17011720/*FUNCTION PentaRef::GetJacobianInvert {{{1*/
    17021721void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,GaussPenta* gauss){
  • issm/trunk/src/c/objects/Elements/PentaRef.h

    r5722 r5726  
    5959                void GetJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
    6060                void GetTriaJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
     61                void GetSegmentJacobianDeterminant(double*  Jdet, double* xyz_list,GaussPenta* gauss);
    6162                void GetJacobianInvert(double*  Jinv, double* xyz_list,GaussPenta* gauss);
    6263                void GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss);
Note: See TracChangeset for help on using the changeset viewer.