Changeset 20871

Timestamp:

07/08/16 12:12:08 (9 years ago)

Author:

adhikari

Message:

CHG: GetAreaSpherical added to compute spherical area

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Element.h (modified) (1 diff)
• Elements/Penta.h (modified) (1 diff)
• Elements/Seg.h (modified) (1 diff)
• Elements/Tetra.h (modified) (1 diff)
• Elements/Tria.cpp (modified) (5 diffs)
• Elements/Tria.h (modified) (1 diff)
• FemModel.cpp (modified) (2 diffs)

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

@@ -298 +298 @@
                 #ifdef _HAVE_SEALEVELRISE_
                 virtual IssmDouble    GetArea3D(void)=0;
+                virtual IssmDouble    GetAreaSpherical(void)=0;
                 virtual IssmDouble    OceanAverage(IssmDouble* Sg)=0;
                 virtual IssmDouble    OceanArea(void)=0;

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

@@ -64 +64 @@
                 void           FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
                 IssmDouble     GetArea3D(void){_error_("not implemented yet!");};
+                IssmDouble     GetAreaSpherical(void){_error_("not implemented yet!");};
                 void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
                 Element*       GetBasalElement(void);

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

@@ -163 +163 @@
                 void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
                 IssmDouble     GetArea3D(void){_error_("not implemented yet!");};
+                IssmDouble     GetAreaSpherical(void){_error_("not implemented yet!");};
 
 #ifdef _HAVE_GIA_

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

@@ -63 +63 @@
                 void        FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){_error_("not implemented yet");};
                 IssmDouble     GetArea3D(void){_error_("not implemented yet!");};
+                IssmDouble     GetAreaSpherical(void){_error_("not implemented yet!");};
                 Element*    GetBasalElement(void){_error_("not implemented yet");};
                 int         GetElementType(void);

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

@@ -912 +912 @@
         return area_fraction*this->GetArea();
 }/*}}}*/
+IssmDouble Tria::GetAreaSpherical(void){/*{{{*/
+
+        bool spherical=true;
+        IssmDouble llr_list[NUMVERTICES][3];
+        IssmDouble x1,y1,z1,x2,y2,z2,x3,y3,z3;
+        IssmDouble arc12,arc23,arc31,semi_peri,excess; 
+
+        /*retrieve coordinates: lat,long,radius */
+        ::GetVerticesCoordinates(&llr_list[0][0],vertices,NUMVERTICES,spherical);
+        x1=llr_list[0][0]/180*PI; y1=llr_list[0][1]/180*PI; z1=llr_list[0][2];
+        x2=llr_list[1][0]/180*PI; y2=llr_list[1][1]/180*PI; z2=llr_list[1][2];
+        x3=llr_list[2][0]/180*PI; y3=llr_list[2][1]/180*PI; z3=llr_list[2][2];
+
+        /*compute great circle distance between vertices */
+        arc12=2.*asin(sqrt(pow(sin((x2-x1)/2),2.0)+cos(x1)*cos(x2)*pow(sin((y2-y1)/2),2)));
+        arc23=2.*asin(sqrt(pow(sin((x3-x2)/2),2.0)+cos(x2)*cos(x3)*pow(sin((y3-y2)/2),2)));
+        arc31=2.*asin(sqrt(pow(sin((x1-x3)/2),2.0)+cos(x3)*cos(x1)*pow(sin((y1-y3)/2),2)));
+
+        /*semi parameter */ 
+        semi_peri=(arc12+arc23+arc31)/2; 
+
+        /*spherical excess */
+        excess=4.*atan(sqrt(tan(semi_peri/2)*tan((semi_peri-arc12)/2)*tan((semi_peri-arc23)/2)*tan((semi_peri-arc31)/2))); 
+
+        /*area = excess*radius^2 */
+        return excess*pow((z1+z2+z3)/3,2); 
+}
+/*}}}*/
 void       Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/
         /*Computeportion of the element that is grounded*/ 
@@ -3543 +3571 @@
 IssmDouble    Tria::OceanArea(void){ /*{{{*/
 
-        if(IsWaterInElement()) return GetArea3D();
+        if(IsWaterInElement()) return GetAreaSpherical();
         else return 0;
 
@@ -3555 +3583 @@
 
                 /*Compute area of element:*/
-                area=GetArea3D();
+                area=GetAreaSpherical();
 
                 /*Average Sg over vertices:*/
@@ -3658 +3686 @@
 
         /*Compute area of element:*/
-        area=GetArea3D();
+        area=GetAreaSpherical();
 
         /*Compute ice thickness change: */
@@ -3761 +3789 @@
 
         /*Compute area of element:*/
-        area=GetArea3D();
+        area=GetAreaSpherical();
 
         /* Where is the centroid of this element?:{{{*/

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

@@ -157 +157 @@
                 IssmDouble         GetArea3D(void);
                 IssmDouble         GetAreaIce(void);
+                IssmDouble         GetAreaSpherical(void);
                 void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
                 int            GetElementType(void);

issm/trunk-jpl/src/c/classes/FemModel.cpp ¶

@@ -2322 +2322 @@
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
                 oceanarea_cpu += element->OceanArea();
-                eartharea_cpu+= element->GetArea3D();
+                eartharea_cpu+= element->GetAreaSpherical();
         }
         ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
@@ -2382 +2382 @@
                 Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
                 oceanarea_cpu += element->OceanArea();
-                eartharea_cpu+= element->GetArea3D();
+                eartharea_cpu+= element->GetAreaSpherical();
         }
         ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );