Changeset 15517

Timestamp:

07/18/13 15:55:23 (12 years ago)

Author:

seroussi

Message:

CHG: startin g new icefront

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Penta.cpp (modified) (2 diffs)
• Elements/Tria.cpp (modified) (5 diffs)
• Elements/Tria.h (modified) (2 diffs)
• gauss/GaussTria.cpp (modified) (1 diff)
• gauss/GaussTria.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

@@ -6949 +6949 @@
 
         /* Start  looping on the number of gaussian points: */
-        gauss=new GaussPenta(3,2);
+        gauss=new GaussPenta(5,5);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
 
@@ -7927 +7927 @@
 
         /* Start  looping on the number of gaussian points: */
-        gauss=new GaussPenta(3,2);
+        gauss=new GaussPenta(5,5);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
 

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -692 +692 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GetAreaCoordinates{{{*/
+void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints){
+        /*Computeportion of the element that is grounded*/ 
+
+        int         i,j,k;
+        IssmDouble  area_init,area_portion;
+        IssmDouble  xyz_bis[3][3];
+
+        GetJacobianDeterminant(&area_init, &xyz_list[0][0],NULL);
+
+        /*Initialize xyz_list with original xyz_list of triangle coordinates*/
+        for(j=0;j<3;j++){ 
+                for(k=0;k<3;j++){
+                        xyz_bis[j][k]=xyz_list[j][k];
+                }
+        }
+        for(i=0;i<numpoints;i++){
+                for(j=0;j<3;j++){ 
+                        for(k=0;k<3;j++){
+                                /*Change appropriate line*/
+                                xyz_bis[j][k]=xyz_zero[i][k];
+                        }
+
+                        /*Compute area fraction*/
+                        GetJacobianDeterminant(&area_portion, &xyz_bis[0][0],NULL);
+                        *(area_coordinates+3*i+j)=area_portion/area_init;
+
+                        /*Reinitialize xyz_list*/
+                        for(k=0;k<3;j++){
+                                /*Reinitialize xyz_list with original coordinates*/
+                                xyz_bis[j][k]=xyz_list[j][k];
+                        }
+                }
+        }
+}
+/*}}}*/
 /*FUNCTION Tria::GetDofList {{{*/
 void  Tria::GetDofList(int** pdoflist, int approximation_enum,int setenum){
@@ -724 +760 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GetGroundedPart{{{*/
+void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){
+        /*Computeportion of the element that is grounded*/ 
+
+        bool               floating=true;
+        int                point;
+        const IssmPDouble  epsilon= 1.e-15;
+        IssmDouble         gl[3];
+        IssmDouble         f1,f2;
+
+        /*Recover parameters and values*/
+        GetInputListOnVertices(&gl[0],GLlevelsetEnum);
+
+        /*Be sure that values are not zero*/
+        if(gl[0]==0) gl[0]=gl[0]+epsilon;
+        if(gl[1]==0) gl[1]=gl[1]+epsilon;
+        if(gl[2]==0) gl[2]=gl[2]+epsilon;
+
+        /*Check that not all nodes are grounded or floating*/
+        if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
+                point=0;
+                f1=1.;
+                f2=1.;
+        }
+        else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
+                point=0;
+                f1=0.;
+                f2=0.;
+        }
+        else{
+                if(gl[0]*gl[1]*gl[2]<0) floating=false;
+
+                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+                        point=2;
+                        f1=gl[2]/(gl[2]-gl[0]);
+                        f2=gl[2]/(gl[2]-gl[1]);
+                }
+                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+                        point=0;
+                        f1=gl[0]/(gl[0]-gl[1]);
+                        f2=gl[0]/(gl[0]-gl[2]);
+                }
+                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+                        point=1;
+                        f1=gl[1]/(gl[1]-gl[2]);
+                        f2=gl[1]/(gl[1]-gl[0]);
+                }
+        }
+        *point1=point;
+        *fraction1=f1;
+        *fraction2=f2;
+        *mainlyfloating=floating;
+}
+/*}}}*/
 /*FUNCTION Tria::GetGroundedPortion{{{*/
 IssmDouble Tria::GetGroundedPortion(IssmDouble* xyz_list){
@@ -826 +916 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetGroundedPart{{{*/
-void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){
+/*FUNCTION Tria::GetSegmentNormal {{{*/
+void Tria:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]){
+
+        /*Build unit outward pointing vector*/
+        IssmDouble vector[2];
+        IssmDouble 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;
+}
+/*}}}*/
+/*FUNCTION Tria::GetZeroLevelsetCoordinates{{{*/
+void Tria::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum){
         /*Computeportion of the element that is grounded*/ 
 
-        bool               floating=true;
-        int                point;
-        const IssmPDouble  epsilon= 1.e-15;
-        IssmDouble         gl[3];
-        IssmDouble         f1,f2;
+        int         normal_orientation;
+        IssmDouble  s1,s2;
+        IssmDouble  levelset[3];
 
         /*Recover parameters and values*/
-        GetInputListOnVertices(&gl[0],GLlevelsetEnum);
-
-        /*Be sure that values are not zero*/
-        if(gl[0]==0) gl[0]=gl[0]+epsilon;
-        if(gl[1]==0) gl[1]=gl[1]+epsilon;
-        if(gl[2]==0) gl[2]=gl[2]+epsilon;
-
-        /*Check that not all nodes are grounded or floating*/
-        if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
-                point=0;
-                f1=1.;
-                f2=1.;
-        }
-        else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
-                point=0;
-                f1=0.;
-                f2=0.;
-        }
-        else{
-                if(gl[0]*gl[1]*gl[2]<0) floating=false;
-
-                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
-                        point=2;
-                        f1=gl[2]/(gl[2]-gl[0]);
-                        f2=gl[2]/(gl[2]-gl[1]);
-                }
-                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
-                        point=0;
-                        f1=gl[0]/(gl[0]-gl[1]);
-                        f2=gl[0]/(gl[0]-gl[2]);
-                }
-                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
-                        point=1;
-                        f1=gl[1]/(gl[1]-gl[2]);
-                        f2=gl[1]/(gl[1]-gl[0]);
-                }
-        }
-        *point1=point;
-        *fraction1=f1;
-        *fraction2=f2;
-        *mainlyfloating=floating;
+        GetInputListOnVertices(&levelset[0],levelsetenum);
+
+        if(levelset[0]*levelset[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+                /*Portion of the segments*/
+                s1=levelset[2]/(levelset[2]-levelset[1]);
+                s2=levelset[2]/(levelset[2]-levelset[0]);
+
+                if(levelset[2]>0) normal_orientation=0;
+                /*New point 1*/
+                *(xyz_zero+3*normal_orientation+0)=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
+                *(xyz_zero+3*normal_orientation+1)=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
+                *(xyz_zero+3*normal_orientation+2)=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
+
+                /*New point 0*/
+                *(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
+                *(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
+                *(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
+        }
+        else if(levelset[1]*levelset[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+                /*Portion of the segments*/
+                s1=levelset[0]/(levelset[0]-levelset[2]);
+                s2=levelset[0]/(levelset[0]-levelset[1]);
+
+                if(levelset[0]>0) normal_orientation=0;
+                /*New point 1*/
+                *(xyz_zero+3*normal_orientation+0)=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
+                *(xyz_zero+3*normal_orientation+1)=xyz_list[0][1]+s1*(xyz_list[2][1]-xyz_list[0][1]);
+                *(xyz_zero+3*normal_orientation+2)=xyz_list[0][2]+s1*(xyz_list[2][2]-xyz_list[0][2]);
+
+                /*New point 2*/
+                *(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[0][0]+s2*(xyz_list[1][0]-xyz_list[0][0]);
+                *(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[0][1]+s2*(xyz_list[1][1]-xyz_list[0][1]);
+                *(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[0][2]+s2*(xyz_list[1][2]-xyz_list[0][2]);
+        }
+        else if(levelset[0]*levelset[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+                /*Portion of the segments*/
+                s1=levelset[1]/(levelset[1]-levelset[0]);
+                s2=levelset[1]/(levelset[1]-levelset[2]);
+
+                if(levelset[1]>0) normal_orientation=0;
+                /*New point 0*/
+                *(xyz_zero+3*normal_orientation+0)=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
+                *(xyz_zero+3*normal_orientation+1)=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
+                *(xyz_zero+3*normal_orientation+2)=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
+
+                /*New point 2*/
+                *(xyz_zero+3*(1-normal_orientation)+0)=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
+                *(xyz_zero+3*(1-normal_orientation)+1)=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
+                *(xyz_zero+3*(1-normal_orientation)+2)=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
+                }
+
 }
 /*}}}*/
@@ -3006 +3120 @@
 
         /*Intermediaries */
-        int            i,j;
-        IssmDouble     ls[3];
-        IssmDouble     xyz_list[NUMVERTICES][3];
-
-        /*Fetch number of nodes and dof for this finite element*/
-        int numnodes = this->NumberofNodes();
-        int numdof   = numnodes*NDOF2;
-        Icefront *icefront=NULL;
+        IssmDouble  ls[3];
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        bool        isfront;
 
         return NULL;
+
         /*Retrieve all inputs and parameters*/
-        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         GetInputListOnVertices(&ls[0],IcelevelsetEnum);
-        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-        GaussTria*     gauss  = new GaussTria(2);
-        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-        /*Create Ice Front if necessary*/
+
+        /*If the level set is awlays <0, there is no ice front here*/
+        isfront = false;
         if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
                 if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
-                        //icefront=new Icefront("2d",inputs,matpar,MacAyealApproximationEnum,analysis_type);
+                        isfront = true;
+                }
+        }
+
+        /*If no front, return NULL*/
+        if(!isfront) return NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int         numnodes = this->NumberofNodes();
+        int         numdof   = numnodes*NDOF2;
+        IssmDouble  rho_ice,rho_water,gravity;
+        IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure,air_pressure;
+        IssmDouble  surface_under_water,base_under_water,pressure;
+        GaussTria*  gauss;
+        IssmDouble* basis = xNew<IssmDouble>(numnodes);
+        IssmDouble  xyz_list_front[2][3];
+        IssmDouble  area_coordinates[2][3];
+        IssmDouble  normal[2];
+
+        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
+        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+        Input* bed_input      =inputs->GetInput(BedEnum);       _assert_(bed_input);
+        rho_water=matpar->GetRhoWater();
+        rho_ice  =matpar->GetRhoIce();
+        gravity  =matpar->GetG();
+        GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,IcelevelsetEnum);
+        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
+        GetSegmentNormal(&normal[0],xyz_list_front);
+
+        /*Start looping on Gaussian points*/
+        gauss=new GaussTria(area_coordinates,3);
+
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+                thickness_input->GetInputValue(&thickness,gauss);
+                bed_input->GetInputValue(&bed,gauss);
+
+                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));
+                ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
+                air_pressure=0;
+
+                pressure = ice_pressure + water_pressure + air_pressure;
+
+                GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
+                GetNodalFunctions(basis,gauss);
+
+                for (int i=0;i<numnodes;i++){
+                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
+                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
                 }
         }
@@ -3037 +3196 @@
         delete gauss;
         return pe;
+
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/Elements/Tria.h

@@ -195 +195 @@
                 ElementVector* CreatePVectorSlope(void);
                 IssmDouble     GetArea(void);
+                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints);
                 int            GetElementType(void);
                 void             GetDofList(int** pdoflist,int approximation_enum,int setenum);
@@ -202 +203 @@
                 IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
                 void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
+                void           GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[2][3]);
+                void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum);
                 void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
                 void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);

issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp

@@ -87 +87 @@
                 coords2[i]=0.5*(b1+b2) + 0.5*seg_coords[i]*(b2-b1);
                 coords3[i]=0.5*(c1+c2) + 0.5*seg_coords[i]*(c2-c1);
+                weights[i]=seg_weights[i];
+        }
+
+        /*Initialize static fields as undefined*/
+        weight=UNDEF;
+        coord1=UNDEF;
+        coord2=UNDEF;
+        coord3=UNDEF;
+
+        /*clean up*/
+        xDelete<double>(seg_coords);
+        xDelete<double>(seg_weights);
+}
+/*}}}*/
+/*FUNCTION GaussTria::GaussTria(IssmDouble area_coordinates,int order) {{{*/
+GaussTria::GaussTria(IssmDouble area_coordinates[2][3],int order){
+
+        /*Intermediaties*/
+        IssmPDouble *seg_coords  = NULL;
+        IssmPDouble *seg_weights = NULL;
+        int     i,index3;
+
+        /*Get Segment gauss points*/
+        numgauss=order;
+        GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);
+
+        /*Allocate GaussTria fields*/
+        coords1=xNew<IssmDouble>(numgauss);
+        coords2=xNew<IssmDouble>(numgauss);
+        coords3=xNew<IssmDouble>(numgauss);
+        weights=xNew<IssmDouble>(numgauss);
+
+        /*Build Triangle Gauss point*/
+        for(i=0;i<numgauss;i++){
+                coords1[i]=0.5*(area_coordinates[0][0]+area_coordinates[1][0]) + 0.5*seg_coords[i]*(area_coordinates[1][0]-area_coordinates[0][0]);
+                coords2[i]=0.5*(area_coordinates[0][1]+area_coordinates[1][1]) + 0.5*seg_coords[i]*(area_coordinates[1][1]-area_coordinates[0][1]);
+                coords3[i]=0.5*(area_coordinates[0][2]+area_coordinates[1][2]) + 0.5*seg_coords[i]*(area_coordinates[1][2]-area_coordinates[0][2]);
                 weights[i]=seg_weights[i];
         }

issm/trunk-jpl/src/c/classes/gauss/GaussTria.h

@@ -31 +31 @@
                 GaussTria(int index1,int index2,int order);
                 GaussTria(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
+                GaussTria(IssmDouble area_coordinates[3][3],int order);
                 ~GaussTria();