Changeset 4946


Ignore:
Timestamp:
08/03/10 10:31:37 (15 years ago)
Author:
Mathieu Morlighem
Message:

improved Ice Front 2d

Location:
issm/trunk/src/c/objects/Loads
Files:
2 edited

Legend:

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

    r4839 r4946  
    397397void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg){
    398398
    399         int i,j;
    400        
    401         const int numgrids=2;
    402         const int NDOF2=2;
    403         int   numberofdofspernode;
    404         const int numdofs=numgrids*NDOF2;
    405         int   doflist[numdofs];
    406         double xyz_list[numgrids][3];
    407 
    408         double x1,y1,x2,y2;
    409 
    410         /* input parameters: */
    411         double thickness_list_element[6];
    412         double thickness_list[2];
    413         double bed_list_element[6];
    414         double bed_list[2];
    415        
    416         int    grid1,grid2;
    417         double normal[2];
    418         double length;
    419         int    element_type;
     399        int       i;
     400        const int numgrids = 2;
     401        const int NDOF2    = 2;
     402        const int numdofs = numgrids *NDOF2;
     403        int       numberofdofspernode;
     404        int       doflist[numdofs];
     405        double    xyz_list[numgrids][3];
    420406
    421407        /*Objects: */
    422         double  pe_g[numdofs]={0.0};
    423         Matpar*  matpar=NULL;
    424         Node**   element_nodes=NULL;
    425         Node**   nodes=NULL;
    426         Element* element=NULL;
     408        double    pe_g[numdofs] = {0.0};
     409        Matpar   *matpar        = NULL;
     410        Node    **nodes         = NULL;
     411        Element  *element       = NULL;
     412
     413        /*Segment*/
     414        double   normal[2];
     415        int      num_gauss;
     416        double*  segment_gauss_coord=NULL;
     417        double*  gauss_weights=NULL;
     418        double   gauss_weight;
     419        double   gauss_coord;
     420        int      ig;
     421        double   Jdet;
     422        double   L[2];
     423        double   thickness,bed,pressure;
     424        double   ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity;
     425        double   surface_under_water,base_under_water;
     426        int      fill;
     427
     428        /*Inputs*/
     429        Input*  thickness_input=NULL;
     430        Input*  bed_input=NULL;
     431        Tria*   tria=NULL;
     432        Penta*  penta=NULL;
    427433
    428434        /*Recover hook objects: */
    429         matpar=(Matpar*)hmatpar->delivers();
     435        matpar =(Matpar*)hmatpar->delivers();
    430436        element=(Element*)helement->delivers();
    431         nodes=(Node**)hnodes->deliverp();
    432 
    433         element_type=element->Enum();
     437        nodes  =(Node**)hnodes->deliverp();
     438        BeamRef* beam=NULL;
    434439
    435440        //check that the element is onbed (collapsed formulation) otherwise:pe=0
    436         if(element_type==PentaEnum){
    437                 if  (!element->GetOnBed()){
    438                         return;
    439                 }
    440         }
    441                
    442         /*Identify which grids are comprised in the segment. */
    443         if(element_type==TriaEnum)element_nodes=(Node**)xmalloc(3*sizeof(Node*));
    444         if(element_type==PentaEnum)element_nodes=(Node**)xmalloc(6*sizeof(Node*));
    445         element->GetNodes((void**)element_nodes);
    446 
    447         /*go through first 3 grids (all grids for tria, bed grids for penta) and identify 1st and 2nd grid: */
    448         for(i=0;i<3;i++){
    449                 if (nodes[0]==element_nodes[i])grid1=i;
    450                 if (nodes[1]==element_nodes[i])grid2=i;
    451         }
    452        
     441        if(element->Enum()==PentaEnum){
     442                if(!element->GetOnBed()) return;
     443                thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
     444                bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
     445        }
     446        else{
     447                thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);
     448                bed_input      =((Tria*)element)->inputs->GetInput(BedEnum);
     449        }
     450
     451        /*Recover inputs: */
     452        inputs->GetParameterValue(&fill,FillEnum);
     453
    453454        /* Get dof list and node coordinates: */
    454455        GetDofList(&doflist[0],&numberofdofspernode);
    455         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    456        
    457         /*Now build thickness_list_element and bed_list_element: */
    458         element->GetThicknessList(&thickness_list_element[0]);
    459         element->GetBedList(&bed_list_element[0]);
    460                
    461         /*Build thickness_list and bed_list: */
    462         thickness_list[0]=thickness_list_element[grid1];
    463         thickness_list[1]=thickness_list_element[grid2];
    464         bed_list[0]=bed_list_element[grid1];
    465         bed_list[1]=bed_list_element[grid2];
    466 
    467         //Recover grid coordinates
    468         x1=xyz_list[0][0];
    469         y1=xyz_list[0][1];
    470         x2=xyz_list[1][0];
    471         y2=xyz_list[1][1];
    472 
    473         //Compute length and normal of segment
    474         normal[0]=cos(atan2(x1-x2,y2-y1));
    475         normal[1]=sin(atan2(x1-x2,y2-y1));
    476         length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
    477 
    478         //Compute load contribution for this segment:
    479         SegmentPressureLoad(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list,bed_list,normal,length);
     456        GetVerticesCoordinates(&xyz_list[0][0],nodes,numgrids);
     457       
     458        /*Get Matpar parameters*/
     459        rho_water=matpar->GetRhoWater();
     460        rho_ice=matpar->GetRhoIce();
     461        gravity=matpar->GetG();
     462
     463        /*Get normal*/
     464        GetNormal(&normal[0],xyz_list);
     465
     466        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
     467        num_gauss=3;
     468        GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
     469
     470        for(ig=0;ig<num_gauss;ig++){
     471
     472                /*Pick up the gaussian point: */
     473                gauss_weight=gauss_weights[ig];
     474                gauss_coord=segment_gauss_coord[ig];
     475
     476                //compute Jacobian of segment
     477                beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord);
     478
     479                /*Get thickness and bed*/
     480                this->GetParameterValue(&thickness,gauss_coord,thickness_input);
     481                this->GetParameterValue(&bed,      gauss_coord,bed_input);
     482
     483                /*Compute total pressure applied to ice front*/
     484                if (fill==WaterEnum){
     485                        //icefront ends in water:
     486                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
     487                        air_pressure=0;
     488
     489                        //Now deal with water pressure
     490                        surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
     491                        base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
     492                        water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
     493                }
     494                else if (fill==AirEnum){
     495                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
     496                        air_pressure=0;
     497                        water_pressure=0;
     498                }
     499                else{
     500                        ISSMERROR("fill type %i not supported yet",fill);
     501                }
     502                pressure = ice_pressure + water_pressure + air_pressure;
     503
     504                /*Get L*/
     505                beam->GetNodalFunctions(&L[0],gauss_coord);
     506
     507                for (i=0;i<numgrids;i++){
     508                        pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i];
     509                        pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i];
     510                }
     511        } //for(ig=0;ig<num_gauss;ig++)
     512
     513        xfree((void**)&segment_gauss_coord);
     514        xfree((void**)&gauss_weights);
    480515               
    481516        /*Plug pe_g into vector: */
    482517        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
    483 
    484         /*Free ressources:*/
    485         xfree((void**)&element_nodes);
    486518
    487519}
     
    13221354}
    13231355/*}}}*/
    1324 /*FUNCTION Icefront::SegmentPressureLoad {{{1*/
    1325 
    1326 void Icefront::SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length){
    1327 
    1328         double   nx,ny;
    1329         double   h1,h2,b1,b2;
    1330         int      num_gauss;
    1331         double*  segment_gauss_coord=NULL;
    1332         double*  gauss_weights=NULL;
    1333         int      ig;
    1334         double   Jdet;
    1335         double   thickness,bed;
    1336         double   ice_pressure,water_pressure,air_pressure;
    1337         double   surface_under_water,base_under_water;
    1338         double   pressure;
    1339         int    fill;
    1340 
    1341         /*Recover inputs: */
    1342         inputs->GetParameterValue(&fill,FillEnum);
    1343 
    1344         nx=normal[0];
    1345         ny=normal[1];
    1346 
    1347 
    1348         //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
    1349         num_gauss=3;
    1350         GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
    1351 
    1352         //recover thickness and bed at two grids
    1353         h1=thickness_list[0];
    1354         h2=thickness_list[1];
    1355         b1=bed_list[0];
    1356         b2=bed_list[1];
    1357 
    1358         //compute Jacobian of segment
    1359         Jdet=1./2.*length;
    1360 
    1361         for(ig=0;ig<num_gauss;ig++){
    1362 
    1363                 thickness=h1*(1+segment_gauss_coord[ig])/2+h2*(1-segment_gauss_coord[ig])/2;
    1364                 bed=b1*(1+segment_gauss_coord[ig])/2+b2*(1-segment_gauss_coord[ig])/2;
    1365 
    1366                 if (fill==WaterEnum){
    1367                         //icefront ends in water:
    1368                         ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    1369                         air_pressure=0;
    1370 
    1371                         //Now deal with water pressure
    1372                         surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
    1373                         base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
    1374                         water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
    1375                 }
    1376                 else if (fill==AirEnum){
    1377                         ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
    1378                         air_pressure=0;
    1379                         water_pressure=0;
    1380                 }
    1381                 else{
    1382                         ISSMERROR("fill type %i not supported yet",fill);
    1383                 }
    1384 
    1385                 pressure = ice_pressure + water_pressure + air_pressure;
    1386 
    1387                 pe_g[2*0+0]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*nx;
    1388                 pe_g[2*0+1]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*ny;
    1389                 pe_g[2*1+0]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*nx;
    1390                 pe_g[2*1+1]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*ny;
    1391 
    1392         } //for(ig=0;ig<num_gauss;ig++)
    1393 
    1394         xfree((void**)&segment_gauss_coord);
    1395         xfree((void**)&gauss_weights);
    1396 }
    1397 /*}}}*/
     1356/*FUNCTION Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/
     1357void Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){
     1358
     1359        /*output: */
     1360        double value;
     1361
     1362        /*dynamic objects pointed to by hooks: */
     1363        Element* element=NULL;
     1364        Node**   nodes=NULL;
     1365
     1366        /*recover objects from hooks: */
     1367        element=(Element*)helement->delivers();
     1368        nodes=(Node**)hnodes->deliverp();
     1369
     1370        /*Get value on Element 1*/
     1371        element->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in);
     1372
     1373        /*Assign output pointers:*/
     1374        *pvalue=value;
     1375}
     1376/*}}}*/
     1377/*FUNCTION Icefront::GetNormal {{{1*/
     1378void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){
     1379
     1380        /*Build unit outward pointing vector*/
     1381        const int numgrids=2;
     1382        double vector[2];
     1383        double norm;
     1384
     1385        vector[0]=xyz_list[1][0] - xyz_list[0][0];
     1386        vector[1]=xyz_list[1][1] - xyz_list[0][1];
     1387
     1388        norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
     1389
     1390        normal[0]= + vector[1]/norm;
     1391        normal[1]= - vector[0]/norm;
     1392}
     1393/*}}}*/
  • issm/trunk/src/c/objects/Loads/Icefront.h

    r4575 r4946  
    7373                void  CreatePVectorDiagnosticStokes( Vec pg);
    7474                void  GetDofList(int* doflist,int* pnumberofdofs);
    75                 void  SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length);
    7675                void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    7776                                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
    7877                void  QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,
    7978                                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
     79                void GetParameterValue(double* pvalue, double gauss_coord,Input* input_in);
     80                void GetNormal(double* normal,double xyz_list[2][3]);
    8081                /*}}}*/
    8182};
Note: See TracChangeset for help on using the changeset viewer.