Changeset 5719


Ignore:
Timestamp:
09/09/10 12:00:59 (15 years ago)
Author:
Mathieu Morlighem
Message:

Fixed IceFront with new Gauss

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

Legend:

Unmodified
Added
Removed

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

    r5660 r5719  
    3232                virtual void   GetSolutionFromInputs(Vec solution)=0;
    3333                virtual void   GetNodes(void** nodes)=0;
     34                virtual int    GetNodeIndex(Node* node)=0;
    3435                virtual bool   GetShelf()=0;
    3536                virtual bool   GetOnBed()=0;

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

    r5707 r5719  
    916916                pnodes[i]=nodes[i];
    917917        }
     918}
     919/*}}}*/
     920/*FUNCTION Penta::GetNodeIndex {{{1*/
     921int Penta::GetNodeIndex(Node* node){
     922
     923        int i;
     924
     925        for(i=0;i<NUMVERTICES;i++){
     926                if(node==nodes[i])
     927                 return i;
     928        }
     929        ISSMERROR("Node provided not found among element nodes");
     930
    918931}
    919932/*}}}*/

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

    r5660 r5719  
    7979                void   DeleteResults(void);
    8080                void   GetNodes(void** nodes);
     81                int    GetNodeIndex(Node* node);
    8182                bool   GetOnBed();
    8283                bool   GetShelf();

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

    r5669 r5719  
    820820}
    821821/*}}}*/
     822/*FUNCTION Tria::GetNodeIndex {{{1*/
     823int Tria::GetNodeIndex(Node* node){
     824
     825        int i;
     826
     827        for(i=0;i<NUMVERTICES;i++){
     828                if(node==nodes[i])
     829                 return i;
     830        }
     831        ISSMERROR("Node provided not found among element nodes");
     832}
     833/*}}}*/
    822834/*FUNCTION Tria::GetOnBed {{{1*/
    823835bool Tria::GetOnBed(){

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

    r5669 r5719  
    7575                void   CreatePVector(Vec pg);
    7676                void   GetNodes(void** nodes);
     77                int    GetNodeIndex(Node* node);
    7778                bool   GetOnBed();
    7879                bool   GetShelf();

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

    r5638 r5719  
    432432        *(J+NDOF2*0+1)=0.5*(y2-y1);
    433433        *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
     434}
     435/*}}}*/
     436/*FUNCTION TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/
     437void TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss){
     438        /*The Jacobian is constant over the element, discard the gaussian points.
     439         * J is assumed to have been allocated.*/
     440
     441        double x1,y1,x2,y2;
     442
     443        x1=*(xyz_list+3*0+0);
     444        y1=*(xyz_list+3*0+1);
     445        x2=*(xyz_list+3*1+0);
     446        y2=*(xyz_list+3*1+1);
     447
     448        *J=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.));
     449}
     450/*}}}*/
     451/*FUNCTION TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss) {{{1*/
     452void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){
     453        /*The Jacobian determinant is constant over the element, discard the gaussian points.
     454         * J is assumed to have been allocated*/
     455
     456        GetSegmentJacobian(Jdet,xyz_list,gauss);
     457        if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
     458
    434459}
    435460/*}}}*/
    … …  
    557582}
    558583/*}}}*/
     584/*FUNCTION TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss) {{{1*/
     585void TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss,int index1,int index2){
     586        /*This routine returns the values of the nodal functions  at the gaussian point.*/
     587
     588        double BasisFunctions[3];
     589
     590        GetNodalFunctions(&BasisFunctions[0],gauss);
     591
     592        ISSMASSERT(index1>=0 && index1<3);
     593        ISSMASSERT(index2>=0 && index2<3);
     594        l1l2[0]=BasisFunctions[index1];
     595        l1l2[1]=BasisFunctions[index2];
     596}
     597/*}}}*/
    559598/*FUNCTION TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss) {{{1*/
    560599void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){

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

    r5638 r5719  
    3737                void GetJacobian(double* J, double* xyz_list,double* gauss);
    3838                void GetJacobian(double* J, double* xyz_list,GaussTria* gauss);
     39                void GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss);
     40                void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss);
    3941                void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss);
    4042                void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss);
    … …  
    4547                void GetNodalFunctions(double* l1l2l3, double* gauss);
    4648                void GetNodalFunctions(double* l1l2l3,GaussTria* gauss);
     49                void GetSegmentNodalFunctions(double* l1l2l3,GaussTria* gauss, int index1,int index2);
    4750                void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss);
    4851                void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss);

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

    r5641 r5719  
    3636        coord3=UNDEF;
    3737
     38}
     39/*}}}*/
     40/*FUNCTION GaussTria::GaussTria(int index1,int index2,int order) {{{1*/
     41GaussTria::GaussTria(int index1,int index2,int order){
     42
     43        /*Intermediaties*/
     44        double *seg_coords  = NULL;
     45        double *seg_weights = NULL;
     46        int     i,index3;
     47
     48        /*Get Segment gauss points*/
     49        numgauss=order;
     50        GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);
     51
     52        /*Allocate GaussTria fields*/
     53        coords1=(double*)xmalloc(numgauss*sizeof(double));
     54        coords2=(double*)xmalloc(numgauss*sizeof(double));
     55        coords3=(double*)xmalloc(numgauss*sizeof(double));
     56        weights=(double*)xmalloc(numgauss*sizeof(double));
     57
     58        /*Reverse grid1 and 2 if necessary*/
     59        if (index1>index2){
     60                index3=index1; index1=index2; index2=index3;
     61                for(i=0;i<numgauss;i++) seg_coords[i]=-seg_coords[i];
     62        }
     63
     64        /*Build Triangle Gauss point*/
     65        if (index1==0 && index2==1){
     66                for(i=0;i<numgauss;i++) coords1[i]=  0.5*(1-seg_coords[i]);
     67                for(i=0;i<numgauss;i++) coords2[i]=1-0.5*(1.-seg_coords[i]);
     68                for(i=0;i<numgauss;i++) coords3[i]=0;
     69                for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
     70        }
     71        else if (index1==0 && index2==2){
     72                for(i=0;i<numgauss;i++) coords1[i]=  0.5*(1-seg_coords[i]);
     73                for(i=0;i<numgauss;i++) coords2[i]=0;
     74                for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
     75                for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
     76        }
     77        else if (index1==1 && index2==2){
     78                for(i=0;i<numgauss;i++) coords1[i]=0;
     79                for(i=0;i<numgauss;i++) coords2[i]=  0.5*(1-seg_coords[i]);
     80                for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
     81                for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
     82        }
     83        else
     84         ISSMERROR("The 2 indices provided are not supported yet (user provided %i and %i)",index1,index2);
     85
     86        /*Initialize static fields as undefined*/
     87        weight=UNDEF;
     88        coord1=UNDEF;
     89        coord2=UNDEF;
     90        coord3=UNDEF;
     91
     92        /*clean up*/
     93        xfree((void**)&seg_coords);
     94        xfree((void**)&seg_weights);
    3895}
    3996/*}}}*/

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

    r5641 r5719  
    3131                GaussTria();
    3232                GaussTria(int order);
     33                GaussTria(int index1,int index2,int order);
    3334                ~GaussTria();
    3435

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

    r5715 r5719  
    425425void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg){
    426426
    427         int       i;
     427        /*Constants*/
    428428        const int numnodes= NUMVERTICESSEG;
    429429        const int numdofs = numnodes *NDOF2;
    430         int*      doflist=NULL;
    431         double    xyz_list[numnodes][3];
    432 
    433         /*Objects: */
    434         double    pe_g[numdofs] = {0.0};
    435 
    436         /*Segment*/
    437         double   normal[2];
    438         int      num_gauss;
    439         double*  segment_gauss_coord=NULL;
    440         double*  gauss_weights=NULL;
    441         double   gauss_weight;
    442         double   gauss_coord;
    443         int      ig;
    444         double   Jdet;
    445         double   L[2];
    446         double   thickness,bed,pressure;
    447         double   ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity;
    448         double   surface_under_water,base_under_water;
    449         int      fill;
    450 
    451         /*Inputs*/
    452         Input*  thickness_input=NULL;
    453         Input*  bed_input=NULL;
    454         Tria*   tria=NULL;
    455 
    456         BeamRef* beam=NULL;
    457 
    458         /*Recover inputs: */
    459         thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);
    460         bed_input      =((Tria*)element)->inputs->GetInput(BedEnum);
    461         inputs->GetParameterValue(&fill,FillEnum);
     430
     431        /*Intermediary*/
     432        int        ig,index1,index2,fill;
     433        double     Jdet;
     434        double     thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity;
     435        double     water_pressure,air_pressure,surface_under_water,base_under_water;
     436        int       *doflist = NULL;
     437        double     xyz_list[numnodes][3];
     438        double     pe_g[numdofs] = {0.0};
     439        double     normal[2];
     440        double     L[2];
     441        GaussTria *gauss;
     442
     443        Tria* tria=((Tria*)element);
    462444
    463445        /* Get dof list and node coordinates: */
    464446        GetDofList(&doflist,MacAyealApproximationEnum);
    465447        GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
    466        
    467         /*Get Matpar parameters*/
     448
     449        /*Recover inputs and parameters */
     450        Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
     451        Input* bed_input      =tria->inputs->GetInput(BedEnum);       ISSMASSERT(bed_input);
     452        inputs->GetParameterValue(&fill,FillEnum);
     453
    468454        rho_water=matpar->GetRhoWater();
    469455        rho_ice  =matpar->GetRhoIce();
    470456        gravity  =matpar->GetG();
    471457
    472         /*Get normal*/
    473458        GetNormal(&normal[0],xyz_list);
    474459
    475         //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
    476         num_gauss=3;
    477         gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
    478 
    479         for(ig=0;ig<num_gauss;ig++){
    480 
    481                 /*Pick up the gaussian point: */
    482                 gauss_weight=gauss_weights[ig];
    483                 gauss_coord=segment_gauss_coord[ig];
    484 
    485                 //compute Jacobian of segment
    486                 beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord);
    487 
    488                 /*Get thickness and bed*/
    489                 this->GetParameterValue(&thickness,gauss_coord,thickness_input);
    490                 this->GetParameterValue(&bed,      gauss_coord,bed_input);
    491 
    492                 /*Compute total pressure applied to ice front*/
    493                 if (fill==WaterEnum){
    494                         //icefront ends in water:
    495                         ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    496                         air_pressure=0;
    497 
    498                         //Now deal with water pressure
    499                         surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
    500                         base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
    501                         water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
     460        index1=tria->GetNodeIndex(nodes[0]);
     461        index2=tria->GetNodeIndex(nodes[1]);
     462        gauss=new GaussTria(index1,index2,3);
     463
     464        for(ig=gauss->begin();ig<gauss->end();ig++){
     465
     466                gauss->GaussPoint(ig);
     467
     468                thickness_input->GetParameterValue(&thickness,gauss);
     469                bed_input->GetParameterValue(&bed,gauss);
     470
     471                switch(fill){
     472                        case WaterEnum:
     473                                surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
     474                                base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
     475                                water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
     476                                break;
     477                        case AirEnum:
     478                                water_pressure=0;
     479                                break;
     480                        default:
     481                                ISSMERROR("fill type %s not supported yet",EnumToString(fill));
    502482                }
    503                 else if (fill==AirEnum){
    504                         ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    505                         air_pressure=0;
    506                         water_pressure=0;
     483                ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
     484                air_pressure=0;
     485                pressure = ice_pressure + water_pressure + air_pressure;
     486
     487                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
     488                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
     489
     490                for (int i=0;i<numnodes;i++){
     491                        pe_g[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*L[i];
     492                        pe_g[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*L[i];
    507493                }
    508                 else{
    509                         ISSMERROR("fill type %s not supported yet",EnumToString(fill));
    510                 }
    511                 pressure = ice_pressure + water_pressure + air_pressure;
    512 
    513                 /*Get L*/
    514                 beam->GetNodalFunctions(&L[0],gauss_coord);
    515 
    516                 for (i=0;i<numnodes;i++){
    517                         pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i];
    518                         pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i];
    519                 }
    520         } //for(ig=0;ig<num_gauss;ig++)
    521 
    522                        
    523         /*Plug pe_g into vector: */
     494        }
     495       
    524496        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
    525497
    526498        /*Free ressources:*/
    527         xfree((void**)&segment_gauss_coord);
    528         xfree((void**)&gauss_weights);
     499        delete gauss;
    529500        xfree((void**)&doflist);
    530 
    531501}
    532502/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.