Changeset 5720


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

Spawn tria for MacAyeal 3d ice front

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5719 r5720  
    504504void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){
    505505
    506         int       i;
    507         const int numnodes= NUMVERTICESSEG;
    508         const int numdofs = numnodes *NDOF2;
    509         int*      doflist=NULL;
    510         double    xyz_list[numnodes][3];
    511 
    512         /*Objects: */
    513         double    pe_g[numdofs] = {0.0};
    514 
    515         /*Segment*/
    516         double   normal[2];
    517         int      num_gauss;
    518         double*  segment_gauss_coord=NULL;
    519         double*  gauss_weights=NULL;
    520         double   gauss_weight;
    521         double   gauss_coord;
    522         int      ig;
    523         double   Jdet;
    524         double   L[2];
    525         double   thickness,bed,pressure;
    526         double   ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity;
    527         double   surface_under_water,base_under_water;
    528         int      fill;
    529 
    530         /*Inputs*/
    531         Input*  thickness_input=NULL;
    532         Input*  bed_input=NULL;
    533         Penta*  penta=NULL;
    534 
    535         BeamRef* beam=NULL;
    536 
    537         //check that the element is onbed (collapsed formulation) otherwise:pe=0
    538         if(!element->GetOnBed()) return;
    539 
    540         /*Recover inputs: */
    541         thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
    542         bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
    543         inputs->GetParameterValue(&fill,FillEnum);
    544 
    545         /* Get dof list and node coordinates: */
    546         GetDofList(&doflist,MacAyealApproximationEnum);
    547         GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
    548 
    549         /*Get Matpar parameters*/
    550         rho_water=matpar->GetRhoWater();
    551         rho_ice  =matpar->GetRhoIce();
    552         gravity  =matpar->GetG();
    553 
    554         /*Get normal*/
    555         GetNormal(&normal[0],xyz_list);
    556 
    557         //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
    558         num_gauss=3;
    559         gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
    560 
    561         for(ig=0;ig<num_gauss;ig++){
    562 
    563                 /*Pick up the gaussian point: */
    564                 gauss_weight=gauss_weights[ig];
    565                 gauss_coord=segment_gauss_coord[ig];
    566 
    567                 //compute Jacobian of segment
    568                 beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord);
    569 
    570                 /*Get thickness and bed*/
    571                 this->GetParameterValue(&thickness,gauss_coord,thickness_input);
    572                 this->GetParameterValue(&bed,      gauss_coord,bed_input);
    573 
    574                 /*Compute total pressure applied to ice front*/
    575                 if (fill==WaterEnum){
    576                         //icefront ends in water:
    577                         ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    578                         air_pressure=0;
    579 
    580                         //Now deal with water pressure
    581                         surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
    582                         base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
    583                         water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
    584                 }
    585                 else if (fill==AirEnum){
    586                         ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    587                         air_pressure=0;
    588                         water_pressure=0;
    589                 }
    590                 else{
    591                         ISSMERROR("fill type %i not supported yet",fill);
    592                 }
    593                 pressure = ice_pressure + water_pressure + air_pressure;
    594 
    595                 /*Get L*/
    596                 beam->GetNodalFunctions(&L[0],gauss_coord);
    597 
    598                 for (i=0;i<numnodes;i++){
    599                         pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i];
    600                         pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i];
    601                 }
    602         } //for(ig=0;ig<num_gauss;ig++)
    603 
    604 
    605         /*Plug pe_g into vector: */
    606         VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
    607 
    608         /*Free ressources:*/
    609         xfree((void**)&segment_gauss_coord);
    610         xfree((void**)&gauss_weights);
    611         xfree((void**)&doflist);
    612 
     506        Icefront* icefront=NULL;
     507        Penta*    penta=NULL;
     508        Tria*     tria=NULL;
     509        bool      onbed;
     510
     511        /*Cast element onto Penta*/
     512        penta   =(Penta*)this->element;
     513
     514        /*Return if not on bed*/
     515        penta->inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     516        if(!onbed) return;
     517
     518        /*Spawn Tria and call MacAyeal2d*/
     519        tria    =(Tria*)penta->SpawnTria(0,1,2);
     520        icefront=(Icefront*)this->copy();
     521        icefront->element=tria;
     522        icefront->inputs->AddInput(new IntInput(TypeEnum,MacAyeal2dIceFrontEnum));
     523        icefront->CreatePVectorDiagnosticMacAyeal2d(pg);
     524
     525        /*clean-up*/
     526        delete tria->matice;
     527        delete tria;
     528        delete icefront;
    613529}
    614530/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.