Changeset 26626

Timestamp:

11/15/21 14:39:11 (3 years ago)

Author:

vverjans

Message:

CHG: IceVolume() in Tria.cpp and Penta.cpp can account for partial ice coverage of frontal elements (currently disabled using false condition)

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Penta.cpp (modified) (3 diffs)
• Penta.h (modified) (1 diff)
• Tria.cpp (modified) (3 diffs)

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

@@ -1236 +1236 @@
         /*return PentaRef field*/
         return this->element_type;
+}
+/*}}}*/
+void       Penta::GetFractionGeometry2D(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){/*{{{*/
+
+  /*Compute portion of element that is grounded based on levelset at the 3 lower vertices of the Penta element*/
+   bool               trapezeisnegative=true;
+   int                point;
+   const IssmPDouble  epsilon= 1.e-15;
+   IssmDouble         f1,f2,phi;
+
+   /*Weights*/
+   Gauss* gauss = NULL;
+   IssmDouble loadweights_g[NUMVERTICES2D];
+   IssmDouble total_weight = 0;
+
+   _assert_(!xIsNan<IssmDouble>(gl[0]));
+   _assert_(!xIsNan<IssmDouble>(gl[1]));
+   _assert_(!xIsNan<IssmDouble>(gl[2]));
+
+   /*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 positive or negative: */
+   if(gl[0]>0 && gl[1]>0 && gl[2]>0){
+      point = 0;
+      f1    = 1.;
+      f2    = 1.;
+   }
+   else if(gl[0]<0 && gl[1]<0 && gl[2]<0){
+      point = 0;
+      f1    = 0.;
+      f2    = 0.;
+   }
+        else{
+                if(gl[0]*gl[1]*gl[2]<0) trapezeisnegative = false;
+
+                /*Find the similar nodes*/
+                if(gl[0]*gl[1]>0){ 
+                        point = 2;
+                        f1    = gl[2]/(gl[2]-gl[0]);
+                        f2    = gl[2]/(gl[2]-gl[1]);
+                }
+                else if(gl[1]*gl[2]>0){ 
+                        point = 0;
+                        f1    = gl[0]/(gl[0]-gl[1]);
+                        f2    = gl[0]/(gl[0]-gl[2]);
+                }
+                else if(gl[0]*gl[2]>0){ 
+                        point = 1;
+                        f1    = gl[1]/(gl[1]-gl[2]);
+                        f2    = gl[1]/(gl[1]-gl[0]);
+                }
+                else _error_("case not possible");
+        }
+        if(trapezeisnegative) phi = 1-f1*f2;
+        else                  phi = f1*f2;
+        
+        /*Compute weights*/
+        gauss = this->NewGauss(point,f1,f2,trapezeisnegative,2);
+
+        total_weight = 0.0;
+        for(int i=0;i<NUMVERTICES2D;i++)weights[i] = 0;
+        while(gauss->next()){
+                GetNodalFunctions(&loadweights_g[0],gauss,P1Enum);
+                for(int i=0;i<NUMVERTICES2D;i++)weights[i] += loadweights_g[i]*gauss->weight;
+                total_weight += gauss->weight;
+        }
+
+        /*Normalize to phi*/
+        if(total_weight>0.){
+                for(int i=0;i<NUMVERTICES2D;i++) weights[i] = weights[i]*phi/total_weight;
+        }
+        else for(int i=0;i<NUMVERTICES2D;i++) weights[i] = 0.0;
+
+        /*Cleanup*/
+        delete gauss;
+        
+        /*Assign output pointers*/
+        *pphi               = phi;
+        *ppoint1            = point;
+        *pfraction1         = f1;
+        *pfraction2         = f2;
+        *ptrapezeisnegative = trapezeisnegative;
 }
 /*}}}*/
@@ -2099 +2184 @@
         IssmDouble base,height,scalefactor;
         IssmDouble xyz_list[NUMVERTICES][3];
+        IssmDouble lsf[NUMVERTICES];
 
         if(!IsIceInElement())return 0;
@@ -2104 +2190 @@
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
-        /*First calculate the area of the base (cross section triangle)
-         * http://en.wikipedia.org/wiki/Pentangle
-         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
-        base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
-
-        if(scaled==true){ //scale for area projection correction
-                Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
-                scalefactor_input->GetInputAverage(&scalefactor);
-                base=base*scalefactor;
-        }
-
-        /*Now get the average height*/
-        height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
+        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+        /*Deal with partially ice-covered elements*/
+        if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
+                bool        istrapneg;
+                int         point;
+                IssmDouble  f1,f2,phi;
+                IssmDouble* heights = xNew<IssmDouble>(NUMVERTICES2D);
+                IssmDouble  weights[NUMVERTICES2D];
+                IssmDouble  lsf2d[NUMVERTICES2D];
+                for(int i=0;i<NUMVERTICES2D;i++){
+                        heights[i] = xyz_list[i+NUMVERTICES2D][2]-xyz_list[i][2];
+                        lsf2d[i]   = lsf[i];
+                }
+                GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&istrapneg,lsf2d);
+                
+                IssmDouble basetot;
+                height = 0.0;
+                for(int i=0;i<NUMVERTICES2D;i++) height += weights[i]*heights[i];
+                basetot = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
+                base    = basetot*phi;  
+        
+                /*Account for scaling factor averaged over subelement 2D area*/
+                if(scaled==true){
+                        IssmDouble* scalefactor_vertices   = xNew<IssmDouble>(NUMVERTICES);
+                        Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
+                        /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
+                        scalefactor = 0.0;
+                        for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]*scalefactor_vertices[i];
+                        base = base*scalefactor;
+                        xDelete<IssmDouble>(scalefactor_vertices);
+                }
+                xDelete<IssmDouble>(heights); 
+        }
+
+        else{ 
+                /*First calculate the area of the base (cross section triangle)
+                 * http://en.wikipedia.org/wiki/Pentangle
+                 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
+                base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
+        
+                if(scaled==true){ //scale for area projection correction
+                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
+                        scalefactor_input->GetInputAverage(&scalefactor);
+                        base=base*scalefactor;
+                }
+        
+                /*Now get the average height*/
+                height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
+        }
 
         /*Return: */

issm/trunk-jpl/src/c/classes/Elements/Penta.h ¶

@@ -86 +86 @@
                 IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
                 void           GetFractionGeometry(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl){_error_("not implemented yet");};
+                void           GetFractionGeometry2D(IssmDouble* weights, IssmDouble* pphi, int* ppoint1,IssmDouble* pfraction1,IssmDouble* pfraction2, bool* ptrapezeisnegative, IssmDouble* gl);
                 void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area,  int levelsetenum){_error_("not implemented yet");};
                 void       GetNodalWeightsAndAreaAndCentroidsFromLeveset(IssmDouble* loadweights, IssmDouble* ploadarea, IssmDouble* platbar, IssmDouble* plongbar, IssmDouble late, IssmDouble longe, IssmDouble area, int levelset1enum, int levelset2enum){_error_("not implemented yet");};

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

@@ -3489 +3489 @@
         parameters->FindParam(&domaintype,DomainTypeEnum);
 
+        /*Relict code
         if(false && IsIcefront()){
                 //Assumption: linear ice thickness profile on element.
@@ -3543 +3544 @@
                 for(i=0;i<numthk;i++)   Haverage+=H[i];
                 Haverage/=IssmDouble(numthk);
-        }
+        }*/
+
+   IssmDouble* lsf = xNew<IssmDouble>(NUMVERTICES);
+   Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+   /*Deal with partially ice-covered elements*/
+        if(false && (lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0)){
+                bool istrapneg;
+      int point;
+      IssmDouble* weights  = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble* surfaces = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble* bases    = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble* Hice     = xNew<IssmDouble>(NUMVERTICES);
+      IssmDouble area_basetot,f1,f2,phi;
+      /*Average thickness over subelement*/
+                Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+      Element::GetInputListOnVertices(&bases[0],BaseEnum);
+      GetFractionGeometry(weights,&phi,&point,&f1,&f2,&istrapneg,lsf);
+      for(int i=0;i<NUMVERTICES;i++) Hice[i] = surfaces[i]-bases[i];
+      Haverage = 0.0;
+      for(int i=0;i<NUMVERTICES;i++) Haverage += weights[i]*Hice[i];
+                /*Get back area of ice-covered base*/
+                area_basetot = this->GetArea();
+                area_base    = phi*area_basetot;
+
+                /*Account for scaling factor averaged over subelement*/
+                if(scaled==true){
+                        IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
+                        Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
+                        scalefactor = 0.0;
+                        for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]*scalefactor_vertices[i];
+                        area_base = area_base*scalefactor;
+                        xDelete<IssmDouble>(scalefactor_vertices);
+                }
+
+                /*Cleanup*/
+                xDelete<IssmDouble>(weights);
+                xDelete<IssmDouble>(surfaces);
+                xDelete<IssmDouble>(bases);
+                xDelete<IssmDouble>(Hice);
+   }
         else{
                 /*First get back the area of the base*/
@@ -3561 +3601 @@
         }
 
+        
         /*Cleanup & return: */
         xDelete<int>(indices);
         xDelete<IssmDouble>(H);
         xDelete<IssmDouble>(SF);
+        xDelete<IssmDouble>(lsf);
 
         if(domaintype==Domain2DverticalEnum){