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

Timestamp:

02/18/16 16:16:01 (9 years ago)

Author:

seroussi

Message:

NEW: compute exact total melting

File:

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

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

@@ -3073 +3073 @@
 
         /*The fbmb[kg yr-1] of one element is area[m2] * melting_rate [kg m^-2 yr^-1]*/
-        IssmDouble base,fbmb,rho_ice;
+        int        point1;
+        bool       mainlyfloating;
+        IssmDouble fbmb=0;
+        IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet;
         IssmDouble Total_Fbmb=0;
         IssmDouble xyz_list[NUMVERTICES][3];
+        Gauss*     gauss     = NULL;
+
+   if(!IsIceInElement())return 0;
 
         /*Get material parameters :*/
         rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
-
-   if(!IsIceInElement())return 0;
-
+        Input* floatingmelt_input = this->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input); 
+        Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
-        /*First calculate the area of the base (cross section triangle)
-         * http://en.wikipedia.org/wiki/Triangle
-         * 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])); // area of element in m2
-
-        /*Now get the average SMB over the element*/
-        Input* fbmb_input = inputs->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fbmb_input);
-        fbmb_input->GetInputAverage(&fbmb);                                                                                                                                                                                             // average melting rate on element in m ice s-1
-   Total_Fbmb=rho_ice*base*fbmb;                                                                                                                                                                                                                        // melting rate on element in kg s-1
+        this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+        /* Start  looping on the number of gaussian points: */
+        gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+                this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
+                floatingmelt_input->GetInputValue(&floatingmelt,gauss);
+                fbmb+=floatingmelt*Jdet*gauss->weight;
+        }
+
+   Total_Fbmb=rho_ice*fbmb;             // from volume to mass
 
         /*Return: */
+        delete gauss;
         return Total_Fbmb;
 }
@@ -3101 +3110 @@
 
         /*The gbmb[kg yr-1] of one element is area[m2] * gounded melting rate [kg m^-2 yr^-1]*/
-        IssmDouble base,gbmb,rho_ice;
+        int        point1;
+        bool       mainlyfloating;
+        IssmDouble gbmb=0;
+        IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet;
         IssmDouble Total_Gbmb=0;
         IssmDouble xyz_list[NUMVERTICES][3];
+        Gauss*     gauss     = NULL;
+
+   if(!IsIceInElement())return 0;
 
         /*Get material parameters :*/
         rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
-
-   if(!IsIceInElement())return 0;
-
+        Input* groundedmelt_input = this->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input); 
+        Input* gllevelset_input = this->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
-        /*First calculate the area of the base (cross section triangle)
-         * http://en.wikipedia.org/wiki/Triangle
-         * 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])); // area of element in m2
-
-        /*Now get the average grounded melting rate over the element*/
-        Input* gbmb_input = inputs->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gbmb_input);
-        gbmb_input->GetInputAverage(&gbmb);                                                                                                                                                                                             // average smb on element in m ice s-1
-   Total_Gbmb=rho_ice*base*gbmb;                                                                                                                                                                                                                        // smb on element in kg s-1
+        this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+        /* Start  looping on the number of gaussian points: */
+        gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+                this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
+                groundedmelt_input->GetInputValue(&groundedmelt,gauss);
+                gbmb+=groundedmelt*Jdet*gauss->weight;
+        }
+
+   Total_Gbmb=rho_ice*gbmb;             // from volume to mass
 
         /*Return: */
+        delete gauss;
         return Total_Gbmb;
 }