Changeset 4946

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:

• Icefront.cpp (modified) (2 diffs)
• Icefront.h (modified) (1 diff)

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

@@ -397 +397 @@
 void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg){
 
-        int i,j;
-        
-        const int numgrids=2;
-        const int NDOF2=2;
-        int   numberofdofspernode;
-        const int numdofs=numgrids*NDOF2;
-        int   doflist[numdofs];
-        double xyz_list[numgrids][3];
-
-        double x1,y1,x2,y2;
-
-        /* input parameters: */
-        double thickness_list_element[6];
-        double thickness_list[2];
-        double bed_list_element[6];
-        double bed_list[2];
-        
-        int    grid1,grid2;
-        double normal[2];
-        double length;
-        int    element_type;
+        int       i;
+        const int numgrids = 2;
+        const int NDOF2    = 2;
+        const int numdofs = numgrids *NDOF2;
+        int       numberofdofspernode;
+        int       doflist[numdofs];
+        double    xyz_list[numgrids][3];
 
         /*Objects: */
-        double  pe_g[numdofs]={0.0};
-        Matpar*  matpar=NULL;
-        Node**   element_nodes=NULL;
-        Node**   nodes=NULL;
-        Element* element=NULL;
+        double    pe_g[numdofs] = {0.0};
+        Matpar   *matpar        = NULL;
+        Node    **nodes         = NULL;
+        Element  *element       = NULL;
+
+        /*Segment*/
+        double   normal[2];
+        int      num_gauss;
+        double*  segment_gauss_coord=NULL;
+        double*  gauss_weights=NULL;
+        double   gauss_weight;
+        double   gauss_coord;
+        int      ig;
+        double   Jdet;
+        double   L[2];
+        double   thickness,bed,pressure;
+        double   ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity;
+        double   surface_under_water,base_under_water;
+        int      fill;
+
+        /*Inputs*/
+        Input*  thickness_input=NULL;
+        Input*  bed_input=NULL;
+        Tria*   tria=NULL;
+        Penta*  penta=NULL;
 
         /*Recover hook objects: */
-        matpar=(Matpar*)hmatpar->delivers();
+        matpar =(Matpar*)hmatpar->delivers();
         element=(Element*)helement->delivers();
-        nodes=(Node**)hnodes->deliverp();
-
-        element_type=element->Enum();
+        nodes  =(Node**)hnodes->deliverp();
+        BeamRef* beam=NULL;
 
         //check that the element is onbed (collapsed formulation) otherwise:pe=0
-        if(element_type==PentaEnum){
-                if  (!element->GetOnBed()){
-                        return;
-                }
-        }
-                
-        /*Identify which grids are comprised in the segment. */
-        if(element_type==TriaEnum)element_nodes=(Node**)xmalloc(3*sizeof(Node*));
-        if(element_type==PentaEnum)element_nodes=(Node**)xmalloc(6*sizeof(Node*));
-        element->GetNodes((void**)element_nodes);
-
-        /*go through first 3 grids (all grids for tria, bed grids for penta) and identify 1st and 2nd grid: */
-        for(i=0;i<3;i++){
-                if (nodes[0]==element_nodes[i])grid1=i;
-                if (nodes[1]==element_nodes[i])grid2=i;
-        }
-        
+        if(element->Enum()==PentaEnum){
+                if(!element->GetOnBed()) return;
+                thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
+                bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
+        }
+        else{
+                thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);
+                bed_input      =((Tria*)element)->inputs->GetInput(BedEnum);
+        }
+
+        /*Recover inputs: */
+        inputs->GetParameterValue(&fill,FillEnum);
+
         /* Get dof list and node coordinates: */
         GetDofList(&doflist[0],&numberofdofspernode);
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-        
-        /*Now build thickness_list_element and bed_list_element: */
-        element->GetThicknessList(&thickness_list_element[0]);
-        element->GetBedList(&bed_list_element[0]);
-                
-        /*Build thickness_list and bed_list: */
-        thickness_list[0]=thickness_list_element[grid1];
-        thickness_list[1]=thickness_list_element[grid2];
-        bed_list[0]=bed_list_element[grid1];
-        bed_list[1]=bed_list_element[grid2];
-
-        //Recover grid coordinates
-        x1=xyz_list[0][0];
-        y1=xyz_list[0][1];
-        x2=xyz_list[1][0];
-        y2=xyz_list[1][1];
-
-        //Compute length and normal of segment
-        normal[0]=cos(atan2(x1-x2,y2-y1));
-        normal[1]=sin(atan2(x1-x2,y2-y1));
-        length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
-
-        //Compute load contribution for this segment:
-        SegmentPressureLoad(pe_g,matpar->GetRhoWater(),matpar->GetRhoIce(),matpar->GetG(),thickness_list,bed_list,normal,length);
+        GetVerticesCoordinates(&xyz_list[0][0],nodes,numgrids);
+        
+        /*Get Matpar parameters*/
+        rho_water=matpar->GetRhoWater();
+        rho_ice=matpar->GetRhoIce();
+        gravity=matpar->GetG();
+
+        /*Get normal*/
+        GetNormal(&normal[0],xyz_list);
+
+        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
+        num_gauss=3;
+        GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
+
+        for(ig=0;ig<num_gauss;ig++){
+
+                /*Pick up the gaussian point: */
+                gauss_weight=gauss_weights[ig];
+                gauss_coord=segment_gauss_coord[ig];
+
+                //compute Jacobian of segment
+                beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord);
+
+                /*Get thickness and bed*/
+                this->GetParameterValue(&thickness,gauss_coord,thickness_input);
+                this->GetParameterValue(&bed,      gauss_coord,bed_input);
+
+                /*Compute total pressure applied to ice front*/
+                if (fill==WaterEnum){
+                        //icefront ends in water: 
+                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
+                        air_pressure=0;
+
+                        //Now deal with water pressure
+                        surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
+                        base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
+                        water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
+                }
+                else if (fill==AirEnum){
+                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
+                        air_pressure=0;
+                        water_pressure=0;
+                }
+                else{
+                        ISSMERROR("fill type %i not supported yet",fill);
+                }
+                pressure = ice_pressure + water_pressure + air_pressure;
+
+                /*Get L*/
+                beam->GetNodalFunctions(&L[0],gauss_coord);
+
+                for (i=0;i<numgrids;i++){
+                        pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i];
+                        pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i];
+                }
+        } //for(ig=0;ig<num_gauss;ig++)
+
+        xfree((void**)&segment_gauss_coord);
+        xfree((void**)&gauss_weights);
                 
         /*Plug pe_g into vector: */
         VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
-
-        /*Free ressources:*/
-        xfree((void**)&element_nodes);
 
 }
@@ -1322 +1354 @@
 }
 /*}}}*/
-/*FUNCTION Icefront::SegmentPressureLoad {{{1*/
-
-void Icefront::SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length){
-
-        double   nx,ny;
-        double   h1,h2,b1,b2;
-        int      num_gauss;
-        double*  segment_gauss_coord=NULL;
-        double*  gauss_weights=NULL;
-        int      ig;
-        double   Jdet;
-        double   thickness,bed;
-        double   ice_pressure,water_pressure,air_pressure;
-        double   surface_under_water,base_under_water;
-        double   pressure;
-        int    fill;
-
-        /*Recover inputs: */
-        inputs->GetParameterValue(&fill,FillEnum);
-
-        nx=normal[0];
-        ny=normal[1];
-
-
-        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
-        num_gauss=3;
-        GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
-
-        //recover thickness and bed at two grids
-        h1=thickness_list[0];
-        h2=thickness_list[1];
-        b1=bed_list[0];
-        b2=bed_list[1];
-
-        //compute Jacobian of segment
-        Jdet=1./2.*length;
-
-        for(ig=0;ig<num_gauss;ig++){
-
-                thickness=h1*(1+segment_gauss_coord[ig])/2+h2*(1-segment_gauss_coord[ig])/2;
-                bed=b1*(1+segment_gauss_coord[ig])/2+b2*(1-segment_gauss_coord[ig])/2;
-
-                if (fill==WaterEnum){
-                        //icefront ends in water: 
-                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
-                        air_pressure=0;
-
-                        //Now deal with water pressure
-                        surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
-                        base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
-                        water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
-                }
-                else if (fill==AirEnum){
-                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
-                        air_pressure=0;
-                        water_pressure=0;
-                }
-                else{
-                        ISSMERROR("fill type %i not supported yet",fill);
-                }
-
-                pressure = ice_pressure + water_pressure + air_pressure;
-
-                pe_g[2*0+0]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*nx;
-                pe_g[2*0+1]+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*ny;
-                pe_g[2*1+0]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*nx;
-                pe_g[2*1+1]+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*ny;
-
-        } //for(ig=0;ig<num_gauss;ig++)
-
-        xfree((void**)&segment_gauss_coord);
-        xfree((void**)&gauss_weights);
-}
-/*}}}*/
+/*FUNCTION Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/
+void Icefront::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){
+
+        /*output: */
+        double value;
+
+        /*dynamic objects pointed to by hooks: */
+        Element* element=NULL;
+        Node**   nodes=NULL;
+
+        /*recover objects from hooks: */
+        element=(Element*)helement->delivers();
+        nodes=(Node**)hnodes->deliverp();
+
+        /*Get value on Element 1*/
+        element->GetParameterValue(&value,nodes[0],nodes[1],gauss_coord,input_in);
+
+        /*Assign output pointers:*/
+        *pvalue=value;
+}
+/*}}}*/
+/*FUNCTION Icefront::GetNormal {{{1*/
+void Icefront:: GetNormal(double* normal,double xyz_list[2][3]){
+
+        /*Build unit outward pointing vector*/
+        const int numgrids=2;
+        double vector[2];
+        double norm;
+
+        vector[0]=xyz_list[1][0] - xyz_list[0][0];
+        vector[1]=xyz_list[1][1] - xyz_list[0][1];
+
+        norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
+
+        normal[0]= + vector[1]/norm;
+        normal[1]= - vector[0]/norm;
+}
+/*}}}*/

issm/trunk/src/c/objects/Loads/Icefront.h

@@ -73 +73 @@
                 void  CreatePVectorDiagnosticStokes( Vec pg);
                 void  GetDofList(int* doflist,int* pnumberofdofs);
-                void  SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length);
                 void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 
                                               double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
                 void  QuadPressureLoadStokes(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 
                                               double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
+                void GetParameterValue(double* pvalue, double gauss_coord,Input* input_in);
+                void GetNormal(double* normal,double xyz_list[2][3]);
                 /*}}}*/
 };