Changeset 23880

Timestamp:

04/22/19 10:18:19 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: merge mass flux calculations

Location:

issm/trunk-jpl/src/c

Files:

• classes/Elements/Tria.cpp (modified) (4 diffs)
• shared/Numerics/constants.h (modified) (1 diff)

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

@@ -2600 +2600 @@
         IssmDouble mass_flux=0.;
         IssmDouble xyz_list[NUMVERTICES][3];
-        IssmDouble normal[2];
-        IssmDouble length,rho_ice;
-        IssmDouble h1,h2;
-        IssmDouble vx1,vx2,vy1,vy2;
-        GaussTria* gauss_1=NULL;
-        GaussTria* gauss_2=NULL;
+        IssmDouble vx1,vx2,vy1,vy2,h1,h2;
 
         /*Get material parameters :*/
-        rho_ice=FindParam(MaterialsRhoIceEnum);
+        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
 
         /*First off, check that this segment belongs to this element: */
@@ -2617 +2612 @@
 
         /*get area coordinates of 0 and 1 locations: */
-        gauss_1=new GaussTria();
+        GaussTria* gauss_1=new GaussTria();
         gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
-        gauss_2=new GaussTria();
+        GaussTria* gauss_2=new GaussTria();
         gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
 
-        normal[0]=cos(atan2(x1-x2,y2-y1));
-        normal[1]=sin(atan2(x1-x2,y2-y1));
-
-        length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
-
+        /*Get segment length and normal (needs to be properly oriented)*/
+        IssmDouble nx=cos(atan2(x1-x2,y2-y1));
+        IssmDouble ny=sin(atan2(x1-x2,y2-y1));
+        IssmDouble length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
+
+        /*Get velocity and thickness*/
+        this->parameters->FindParam(&domaintype,DomainTypeEnum);
         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
-        this->parameters->FindParam(&domaintype,DomainTypeEnum);
         Input* vx_input=NULL;
         Input* vy_input=NULL;
@@ -2648 +2644 @@
 
         mass_flux= rho_ice*length*(  
-                                (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
-                                (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
+                                (1./3.*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*nx+
+                                (1./3.*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*ny
                                 );
 
@@ -2659 +2655 @@
 /*}}}*/
 IssmDouble Tria::MassFlux(IssmDouble* segment){/*{{{*/
-
-        int        domaintype;
-        IssmDouble mass_flux=0.;
-        IssmDouble xyz_list[NUMVERTICES][3];
-        IssmDouble normal[2];
-        IssmDouble length,rho_ice;
-        IssmDouble x1,y1,x2,y2,h1,h2;
-        IssmDouble vx1,vx2,vy1,vy2;
-        GaussTria* gauss_1=NULL;
-        GaussTria* gauss_2=NULL;
-
-        /*Get material parameters :*/
-        rho_ice=FindParam(MaterialsRhoIceEnum);
-
-        /*First off, check that this segment belongs to this element: */
-        if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
-
-        /*Recover segment node locations: */
-        x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
-
-        /*Get xyz list: */
-        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-        /*get area coordinates of 0 and 1 locations: */
-        gauss_1=new GaussTria();
-        gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
-        gauss_2=new GaussTria();
-        gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
-
-        normal[0]=cos(atan2(x1-x2,y2-y1));
-        normal[1]=sin(atan2(x1-x2,y2-y1));
-
-        length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
-
-        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
-        this->parameters->FindParam(&domaintype,DomainTypeEnum);
-        Input* vx_input=NULL;
-        Input* vy_input=NULL;
-        if(domaintype==Domain2DhorizontalEnum){
-                vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-                vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-        }
-        else{
-                vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
-                vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
-        }
-
-        thickness_input->GetInputValue(&h1, gauss_1);
-        thickness_input->GetInputValue(&h2, gauss_2);
-        vx_input->GetInputValue(&vx1,gauss_1);
-        vx_input->GetInputValue(&vx2,gauss_2);
-        vy_input->GetInputValue(&vy1,gauss_1);
-        vy_input->GetInputValue(&vy2,gauss_2);
-
-        mass_flux= rho_ice*length*(  
-                                (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
-                                (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
-                                );
-
-        /*clean up and return:*/
-        delete gauss_1;
-        delete gauss_2;
-        return mass_flux;
+        return this->MassFlux(segment[0],segment[1],segment[2],segment[3],reCast<int>(segment[4]));
 }
 /*}}}*/

issm/trunk-jpl/src/c/shared/Numerics/constants.h

@@ -7 +7 @@
 
 #define UNDEF -9999
-#define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
 #define SQRT2 1.414213562373095048801688724209698078569671875376948073176679738
 #define SQRT3 1.732050807568877293527446341505872366942805253810380628055806979