Changeset 18957

Timestamp:

12/11/14 01:32:48 (10 years ago)

Author:

jbondzio

Message:

CHG: compute Tria ice volume now also over elements intersected by 0-levelset

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (1 diff)

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

@@ -1504 +1504 @@
 IssmDouble Tria::IceVolume(void){/*{{{*/
 
-        /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
-        IssmDouble base,surface,bed;
-        IssmDouble xyz_list[NUMVERTICES][3];
-
-        if(!IsIceInElement())return 0;
-
-        /*First get back the area of the base*/
-        base=this->GetArea();
-
-        /*Now get the average height*/
-        Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
-        Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
-        surface_input->GetInputAverage(&surface);
-        base_input->GetInputAverage(&bed);
-
-        /*Return: */
+        /*The volume of a truncated prism is area_base * 1/numedges sum(length of edges)*/
+
+        /*Intermediaries*/
+        int i, numiceverts;
+        IssmDouble area_base,surface,base,Haverage;
+        IssmDouble Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
+        IssmDouble s[2]; // s:fraction of intersected triangle edges, that lies inside ice
+        int* indices=NULL;
+        IssmDouble* H=NULL;
+
+        if(!IsIceInElement())return 0.;
+
         int domaintype;
         parameters->FindParam(&domaintype,DomainTypeEnum);
+
+        if(IsIcefront()){
+                area_base=this->GetAreaIce();
+                //Assumption: linear ice thickness profile on element. 
+                //Hence ice thickness at intersection of levelset function with triangle edge is linear interpolation of ice thickness at vertices.
+                this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
+                GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+                GetInputListOnVertices(&bases[0],BaseEnum);
+                for(i=0;i<NUMVERTICES;i++) Haux[i]= surfaces[indices[i]]-bases[indices[i]]; //sort thicknesses in ice/noice
+                int numthk=numiceverts+2;
+                H=xNew<IssmDouble>(numthk);
+                switch(numiceverts){
+                        case 1: // average over triangle 
+                                H[0]=Haux[0];
+                                H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
+                                H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
+                                break;
+                        case 2: // average over quadrangle
+                                H[0]=Haux[0];
+                                H[1]=Haux[1];
+                                H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
+                                H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
+                                break;
+                        default:
+                                _error_("Number of ice covered vertices wrong in Tria::IceVolume()");
+                                break;
+                }
+                Haverage=0.;
+                for(i=0;i<numthk;i++)   Haverage+=H[i];
+                Haverage/=IssmDouble(numthk);
+        }
+        else{
+                /*First get back the area of the base*/
+                area_base=this->GetArea();
+
+                /*Now get the average height*/
+                Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
+                Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
+                surface_input->GetInputAverage(&surface);
+                base_input->GetInputAverage(&base);
+                Haverage=surface-base;
+        }
+
+        /*Cleanup & return: */
+        xDelete<int>(indices);
+        xDelete<IssmDouble>(H);
+
         if(domaintype==Domain2DverticalEnum){
-          return base;
+          return area_base;
         }
         else{
-          return base*(surface-bed);
+          return area_base*Haverage;
         }
 }