Changeset 18956

Timestamp:

12/11/14 01:20:21 (10 years ago)

Author:

jbondzio

Message:

ADD: Computation of Tria-area covered by ice. Computation of determination of ice distribution in element.

Location:

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

Files:

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

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

@@ -879 +879 @@
 }
 /*}}}*/
+IssmDouble Tria::GetAreaIce(void){/*{{{*/
+
+        /*return area of element covered by ice*/
+        /*Intermediaries*/
+        int numiceverts;
+        IssmDouble area_fraction;
+        IssmDouble s[2]; // s:fraction of intersected triangle edges that lie inside ice
+        int* indices=NULL;
+
+        this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
+
+        switch (numiceverts){
+                case 0: // no vertex has ice: element is ice free
+                        area_fraction=0.;
+                        break;
+                case 1: // one vertex has ice: get area of triangle
+                        area_fraction=s[0]*s[1];
+                        break;
+                case 2: // two vertices have ice: get area of quadrangle
+                        area_fraction=s[0]+s[1]-s[0]*s[1];
+                        break;
+                case NUMVERTICES: // all vertices have ice: return triangle area
+                        area_fraction=1.;
+                        break;
+                default:
+                        _error_("Wrong number of ice vertices in Tria::GetAreaIce!");
+                        break;
+        }
+        _assert_((area_fraction>=0.) && (area_fraction<=1.));
+
+        xDelete<int>(indices);
+        return area_fraction*this->GetArea();
+}/*}}}*/
 void       Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/
         /*Computeportion of the element that is grounded*/ 
@@ -1184 +1217 @@
         xDelete<int>(indicesfront);
 }/*}}}*/
+void                    Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
+        
+        /* GetLevelsetIntersection computes: 
+         * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
+         * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/
+
+        /*Intermediaries*/
+        int i, numiceverts, numnoiceverts;
+        int ind0, ind1, lastindex;
+        int indices_ice[NUMVERTICES],indices_noice[NUMVERTICES];
+        IssmDouble lsf[NUMVERTICES];
+        int* indices = xNew<int>(NUMVERTICES);
+
+        /*Retrieve all inputs and parameters*/
+        GetInputListOnVertices(&lsf[0],levelset_enum);
+
+        /* Determine distribution of ice over element.
+         * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
+        lastindex=0;
+        for(i=0;i<NUMVERTICES;i++){ // go backwards along vertices, and check for sign change
+                ind0=(NUMVERTICES-i)%NUMVERTICES;
+                ind1=(ind0-1+NUMVERTICES)%NUMVERTICES;
+                if((lsf[ind0]-level)*(lsf[ind1]-level)<=0.){ // levelset has been crossed, find last index belonging to segment
+                        if(lsf[ind1]==level) //if levelset intersects 2nd vertex, choose this vertex as last
+                                lastindex=ind1;
+                        else
+                                lastindex=ind0;
+                        break;
+                }
+        }
+
+        numiceverts=0;
+        numnoiceverts=0;
+        for(i=0;i<NUMVERTICES;i++){
+                ind0=(lastindex+i)%NUMVERTICES;
+                if(lsf[i]<=level){
+                        indices_ice[numiceverts]=i;
+                        numiceverts++;
+                }
+                else{
+                        indices_noice[numnoiceverts]=i;
+                        numnoiceverts++;
+                }
+        }
+        //merge indices 
+        for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
+        for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
+
+        switch (numiceverts){
+                case 0: // no vertex has ice: element is ice free, no intersection
+                        for(i=0;i<2;i++)
+                                fraction[i]=0.;
+                        break;
+                case 1: // one vertex has ice:
+                        for(i=0;i<2;i++){
+                                fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
+                        }
+                        break;
+                case 2: // two vertices have ice: get area of quadrangle
+                        for(i=0;i<2;i++){
+                                fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
+                        }
+                        break;
+                case NUMVERTICES: // all vertices have ice: return triangle area
+                        for(i=0;i<2;i++)
+                                fraction[i]=1.;
+                        break;
+                default:
+                        _error_("Wrong number of ice vertices in Tria::GetLevelsetIntersection!");
+                        break;
+        }
+
+        *pindices=indices;
+        *pnumiceverts=numiceverts;
+}
+/*}}}*/
 void       Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
         

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

@@ -143 +143 @@
                 void           AddInput(int input_enum, IssmDouble* values, int interpolation_enum);
                 IssmDouble     GetArea(void);
+                IssmDouble      GetAreaIce(void);
                 void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
+                void            GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level);
                 int            GetElementType(void);
                 void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);