Changeset 5715


Ignore:
Timestamp:
09/09/10 09:23:48 (15 years ago)
Author:
Mathieu Morlighem
Message:

Separated MacAyeal2d and MacAyeal3d in icefront (easier to handle Gauss points)

Location:
issm/trunk/src/c
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r5707 r5715  
    115115        RiftfrontEnum,
    116116        SegmentRiftfrontEnum,
    117         MacAyealIceFrontEnum,
     117        MacAyeal2dIceFrontEnum,
     118        MacAyeal3dIceFrontEnum,
    118119        PattynIceFrontEnum,
    119120        StokesIceFrontEnum,
  • issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

    r5707 r5715  
    100100                case RiftfrontEnum : return "Riftfront";
    101101                case SegmentRiftfrontEnum : return "SegmentRiftfront";
    102                 case MacAyealIceFrontEnum : return "MacAyealIceFront";
     102                case MacAyeal2dIceFrontEnum : return "MacAyeal2dIceFront";
     103                case MacAyeal3dIceFrontEnum : return "MacAyeal3dIceFront";
    103104                case PattynIceFrontEnum : return "PattynIceFront";
    104105                case StokesIceFrontEnum : return "StokesIceFront";
  • issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

    r5707 r5715  
    9898        else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    9999        else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
    100         else if (strcmp(name,"MacAyealIceFront")==0) return MacAyealIceFrontEnum;
     100        else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum;
     101        else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum;
    101102        else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum;
    102103        else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

    r5607 r5715  
    6060
    6161                /*Create and  add load: */
    62                 if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum)){
    63                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyealIceFrontEnum,DiagnosticHorizAnalysisEnum));
     62                if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum) && iomodel->dim==2){
     63                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal2dIceFrontEnum,DiagnosticHorizAnalysisEnum));
     64                        count++;
     65                }
     66                else if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum) && iomodel->dim==3){
     67                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
    6468                        count++;
    6569                }
     
    7377                }
    7478                else if ((int)*(iomodel->elements_type+element)==(MacAyealPattynApproximationEnum)){
    75                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyealIceFrontEnum,DiagnosticHorizAnalysisEnum));
     79                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
    7680                        count++;
    7781                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));
  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5714 r5715  
    6868        icefront_mparid=iomodel->numberofelements+1; //matlab indexing
    6969
    70         if (in_icefront_type==MacAyealIceFrontEnum){
     70        if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
    7171                icefront_node_ids[0]=iomodel->nodecounter+(int)*(iomodel->pressureload+segment_width*i+0);
    7272                icefront_node_ids[1]=iomodel->nodecounter+(int)*(iomodel->pressureload+segment_width*i+1);
     
    330330                int type;
    331331                inputs->GetParameterValue(&type,TypeEnum);
    332                 if (type==MacAyealIceFrontEnum){
    333                         CreatePVectorDiagnosticMacAyeal( pg);
     332                if (type==MacAyeal2dIceFrontEnum){
     333                        CreatePVectorDiagnosticMacAyeal2d( pg);
     334                }
     335                else if (type==MacAyeal3dIceFrontEnum){
     336                        CreatePVectorDiagnosticMacAyeal3d( pg);
    334337                }
    335338                else if (type==PattynIceFrontEnum){
     
    419422
    420423/*Icefront numerics: */
    421 /*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal{{{1*/
    422 void Icefront::CreatePVectorDiagnosticMacAyeal( Vec pg){
     424/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{1*/
     425void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg){
    423426
    424427        int       i;
     
    450453        Input*  bed_input=NULL;
    451454        Tria*   tria=NULL;
    452         Penta*  penta=NULL;
    453455
    454456        BeamRef* beam=NULL;
    455457
    456         //check that the element is onbed (collapsed formulation) otherwise:pe=0
    457         if(element->Enum()==PentaEnum){
    458                 if(!element->GetOnBed()) return;
    459                 thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
    460                 bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
    461         }
    462         else{
    463                 thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);
    464                 bed_input      =((Tria*)element)->inputs->GetInput(BedEnum);
    465         }
    466 
    467458        /*Recover inputs: */
     459        thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);
     460        bed_input      =((Tria*)element)->inputs->GetInput(BedEnum);
    468461        inputs->GetParameterValue(&fill,FillEnum);
    469462
     
    514507                }
    515508                else{
    516                         ISSMERROR("fill type %i not supported yet",fill);
     509                        ISSMERROR("fill type %s not supported yet",EnumToString(fill));
    517510                }
    518511                pressure = ice_pressure + water_pressure + air_pressure;
     
    528521
    529522                       
     523        /*Plug pe_g into vector: */
     524        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
     525
     526        /*Free ressources:*/
     527        xfree((void**)&segment_gauss_coord);
     528        xfree((void**)&gauss_weights);
     529        xfree((void**)&doflist);
     530
     531}
     532/*}}}*/
     533/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{1*/
     534void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){
     535
     536        int       i;
     537        const int numnodes= NUMVERTICESSEG;
     538        const int numdofs = numnodes *NDOF2;
     539        int*      doflist=NULL;
     540        double    xyz_list[numnodes][3];
     541
     542        /*Objects: */
     543        double    pe_g[numdofs] = {0.0};
     544
     545        /*Segment*/
     546        double   normal[2];
     547        int      num_gauss;
     548        double*  segment_gauss_coord=NULL;
     549        double*  gauss_weights=NULL;
     550        double   gauss_weight;
     551        double   gauss_coord;
     552        int      ig;
     553        double   Jdet;
     554        double   L[2];
     555        double   thickness,bed,pressure;
     556        double   ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity;
     557        double   surface_under_water,base_under_water;
     558        int      fill;
     559
     560        /*Inputs*/
     561        Input*  thickness_input=NULL;
     562        Input*  bed_input=NULL;
     563        Penta*  penta=NULL;
     564
     565        BeamRef* beam=NULL;
     566
     567        //check that the element is onbed (collapsed formulation) otherwise:pe=0
     568        if(!element->GetOnBed()) return;
     569
     570        /*Recover inputs: */
     571        thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
     572        bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
     573        inputs->GetParameterValue(&fill,FillEnum);
     574
     575        /* Get dof list and node coordinates: */
     576        GetDofList(&doflist,MacAyealApproximationEnum);
     577        GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
     578
     579        /*Get Matpar parameters*/
     580        rho_water=matpar->GetRhoWater();
     581        rho_ice  =matpar->GetRhoIce();
     582        gravity  =matpar->GetG();
     583
     584        /*Get normal*/
     585        GetNormal(&normal[0],xyz_list);
     586
     587        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
     588        num_gauss=3;
     589        gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
     590
     591        for(ig=0;ig<num_gauss;ig++){
     592
     593                /*Pick up the gaussian point: */
     594                gauss_weight=gauss_weights[ig];
     595                gauss_coord=segment_gauss_coord[ig];
     596
     597                //compute Jacobian of segment
     598                beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord);
     599
     600                /*Get thickness and bed*/
     601                this->GetParameterValue(&thickness,gauss_coord,thickness_input);
     602                this->GetParameterValue(&bed,      gauss_coord,bed_input);
     603
     604                /*Compute total pressure applied to ice front*/
     605                if (fill==WaterEnum){
     606                        //icefront ends in water:
     607                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
     608                        air_pressure=0;
     609
     610                        //Now deal with water pressure
     611                        surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
     612                        base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
     613                        water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
     614                }
     615                else if (fill==AirEnum){
     616                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
     617                        air_pressure=0;
     618                        water_pressure=0;
     619                }
     620                else{
     621                        ISSMERROR("fill type %i not supported yet",fill);
     622                }
     623                pressure = ice_pressure + water_pressure + air_pressure;
     624
     625                /*Get L*/
     626                beam->GetNodalFunctions(&L[0],gauss_coord);
     627
     628                for (i=0;i<numnodes;i++){
     629                        pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i];
     630                        pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i];
     631                }
     632        } //for(ig=0;ig<num_gauss;ig++)
     633
     634
    530635        /*Plug pe_g into vector: */
    531636        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
     
    818923               
    819924        /*How many nodes? :*/
    820         if(type==MacAyealIceFrontEnum)numberofnodes=2;
    821         else numberofnodes=4;
     925        if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
     926         numberofnodes=2;
     927        else
     928         numberofnodes=4;
    822929       
    823930        /*Figure out size of doflist: */
  • issm/trunk/src/c/objects/Loads/Icefront.h

    r5714 r5715  
    7676                /*}}}*/
    7777                /*Load management: {{{1*/
    78                 void  CreatePVectorDiagnosticMacAyeal( Vec pg);
    79                 void  CreatePVectorDiagnosticPattyn( Vec pg);
    80                 void  CreatePVectorDiagnosticStokes( Vec pg);
     78                void  CreatePVectorDiagnosticMacAyeal2d(Vec pg);
     79                void  CreatePVectorDiagnosticMacAyeal3d(Vec pg);
     80                void  CreatePVectorDiagnosticPattyn(Vec pg);
     81                void  CreatePVectorDiagnosticStokes(Vec pg);
    8182                void  GetDofList(int** pdoflist,int approximation_enum=0);
    8283                void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
Note: See TracChangeset for help on using the changeset viewer.