Changeset 5715

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:

• EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• EnumDefinitions/EnumToString.cpp (modified) (1 diff)
• EnumDefinitions/StringToEnum.cpp (modified) (1 diff)
• modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (modified) (2 diffs)
• objects/Loads/Icefront.cpp (modified) (7 diffs)
• objects/Loads/Icefront.h (modified) (1 diff)

issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

@@ -115 +115 @@
         RiftfrontEnum,
         SegmentRiftfrontEnum,
-        MacAyealIceFrontEnum,
+        MacAyeal2dIceFrontEnum,
+        MacAyeal3dIceFrontEnum,
         PattynIceFrontEnum,
         StokesIceFrontEnum,

issm/trunk/src/c/EnumDefinitions/EnumToString.cpp

@@ -100 +100 @@
                 case RiftfrontEnum : return "Riftfront";
                 case SegmentRiftfrontEnum : return "SegmentRiftfront";
-                case MacAyealIceFrontEnum : return "MacAyealIceFront";
+                case MacAyeal2dIceFrontEnum : return "MacAyeal2dIceFront";
+                case MacAyeal3dIceFrontEnum : return "MacAyeal3dIceFront";
                 case PattynIceFrontEnum : return "PattynIceFront";
                 case StokesIceFrontEnum : return "StokesIceFront";

issm/trunk/src/c/EnumDefinitions/StringToEnum.cpp

@@ -98 +98 @@
         else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
         else if (strcmp(name,"SegmentRiftfront")==0) return SegmentRiftfrontEnum;
-        else if (strcmp(name,"MacAyealIceFront")==0) return MacAyealIceFrontEnum;
+        else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum;
+        else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum;
         else if (strcmp(name,"PattynIceFront")==0) return PattynIceFrontEnum;
         else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;

issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

@@ -60 +60 @@
 
                 /*Create and  add load: */
-                if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum)){
-                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyealIceFrontEnum,DiagnosticHorizAnalysisEnum));
+                if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum) && iomodel->dim==2){
+                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal2dIceFrontEnum,DiagnosticHorizAnalysisEnum));
+                        count++;
+                }
+                else if ((int)*(iomodel->elements_type+element)==(MacAyealApproximationEnum) && iomodel->dim==3){
+                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
                         count++;
                 }
@@ -73 +77 @@
                 }
                 else if ((int)*(iomodel->elements_type+element)==(MacAyealPattynApproximationEnum)){
-                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyealIceFrontEnum,DiagnosticHorizAnalysisEnum));
+                        loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,MacAyeal3dIceFrontEnum,DiagnosticHorizAnalysisEnum));
                         count++;
                         loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,PattynIceFrontEnum,DiagnosticHorizAnalysisEnum));

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

@@ -68 +68 @@
         icefront_mparid=iomodel->numberofelements+1; //matlab indexing
 
-        if (in_icefront_type==MacAyealIceFrontEnum){
+        if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
                 icefront_node_ids[0]=iomodel->nodecounter+(int)*(iomodel->pressureload+segment_width*i+0);
                 icefront_node_ids[1]=iomodel->nodecounter+(int)*(iomodel->pressureload+segment_width*i+1);
@@ -330 +330 @@
                 int type;
                 inputs->GetParameterValue(&type,TypeEnum);
-                if (type==MacAyealIceFrontEnum){
-                        CreatePVectorDiagnosticMacAyeal( pg);
+                if (type==MacAyeal2dIceFrontEnum){
+                        CreatePVectorDiagnosticMacAyeal2d( pg);
+                }
+                else if (type==MacAyeal3dIceFrontEnum){
+                        CreatePVectorDiagnosticMacAyeal3d( pg);
                 }
                 else if (type==PattynIceFrontEnum){
@@ -419 +422 @@
 
 /*Icefront numerics: */
-/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal{{{1*/
-void Icefront::CreatePVectorDiagnosticMacAyeal( Vec pg){
+/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{1*/
+void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg){
 
         int       i;
@@ -450 +453 @@
         Input*  bed_input=NULL;
         Tria*   tria=NULL;
-        Penta*  penta=NULL;
 
         BeamRef* beam=NULL;
 
-        //check that the element is onbed (collapsed formulation) otherwise:pe=0
-        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: */
+        thickness_input=((Tria*)element)->inputs->GetInput(ThicknessEnum);
+        bed_input      =((Tria*)element)->inputs->GetInput(BedEnum);
         inputs->GetParameterValue(&fill,FillEnum);
 
@@ -514 +507 @@
                 }
                 else{
-                        ISSMERROR("fill type %i not supported yet",fill);
+                        ISSMERROR("fill type %s not supported yet",EnumToString(fill));
                 }
                 pressure = ice_pressure + water_pressure + air_pressure;
@@ -528 +521 @@
 
                         
+        /*Plug pe_g into vector: */
+        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
+
+        /*Free ressources:*/
+        xfree((void**)&segment_gauss_coord);
+        xfree((void**)&gauss_weights);
+        xfree((void**)&doflist);
+
+}
+/*}}}*/
+/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{1*/
+void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){
+
+        int       i;
+        const int numnodes= NUMVERTICESSEG;
+        const int numdofs = numnodes *NDOF2;
+        int*      doflist=NULL;
+        double    xyz_list[numnodes][3];
+
+        /*Objects: */
+        double    pe_g[numdofs] = {0.0};
+
+        /*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;
+        Penta*  penta=NULL;
+
+        BeamRef* beam=NULL;
+
+        //check that the element is onbed (collapsed formulation) otherwise:pe=0
+        if(!element->GetOnBed()) return;
+
+        /*Recover inputs: */
+        thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
+        bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
+        inputs->GetParameterValue(&fill,FillEnum);
+
+        /* Get dof list and node coordinates: */
+        GetDofList(&doflist,MacAyealApproximationEnum);
+        GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
+
+        /*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<numnodes;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++)
+
+
         /*Plug pe_g into vector: */
         VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
@@ -818 +923 @@
                 
         /*How many nodes? :*/
-        if(type==MacAyealIceFrontEnum)numberofnodes=2;
-        else numberofnodes=4;
+        if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
+         numberofnodes=2;
+        else 
+         numberofnodes=4;
         
         /*Figure out size of doflist: */

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

@@ -76 +76 @@
                 /*}}}*/
                 /*Load management: {{{1*/
-                void  CreatePVectorDiagnosticMacAyeal( Vec pg);
-                void  CreatePVectorDiagnosticPattyn( Vec pg);
-                void  CreatePVectorDiagnosticStokes( Vec pg);
+                void  CreatePVectorDiagnosticMacAyeal2d(Vec pg);
+                void  CreatePVectorDiagnosticMacAyeal3d(Vec pg);
+                void  CreatePVectorDiagnosticPattyn(Vec pg);
+                void  CreatePVectorDiagnosticStokes(Vec pg);
                 void  GetDofList(int** pdoflist,int approximation_enum=0);
                 void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list,