Changeset 24368

Timestamp:

11/20/19 09:58:23 (5 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added GL flux to 3d models

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Penta.cpp (modified) (1 diff)
• Elements/Penta.h (modified) (1 diff)
• Elements/Tria.cpp (modified) (1 diff)
• FemModel.cpp (modified) (1 diff)

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

@@ -1928 +1928 @@
         /*Clean up and return*/
         return groundedarea;
+}
+/*}}}*/
+IssmDouble Penta::GroundinglineMassFlux(bool scaled){/*{{{*/
+
+        /*Make sure there is a grounding line here*/
+        if(!IsOnBase()) return 0;
+        if(!IsIceInElement()) return 0;
+        if(!IsZeroLevelset(MaskGroundediceLevelsetEnum)) return 0;
+
+        /*Scaled not implemented yet...*/
+        _assert_(!scaled);
+
+        int               domaintype,index1,index2;
+        const IssmPDouble epsilon = 1.e-15;
+        IssmDouble        s1,s2;
+        IssmDouble        gl[NUMVERTICES];
+        IssmDouble        xyz_front[2][3];
+
+        IssmDouble *xyz_list = NULL;
+        this->GetVerticesCoordinates(&xyz_list);
+
+        /*Recover parameters and values*/
+        GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+
+        /*Be sure that values are not zero*/
+        if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+        if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+        if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+        int pt1 = 0;
+        int pt2 = 1;
+        if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+
+                /*Portion of the segments*/
+                s1=gl[2]/(gl[2]-gl[1]);
+                s2=gl[2]/(gl[2]-gl[0]);
+                if(gl[2]<0.){
+                        pt1 = 1; pt2 = 0;
+                }
+                xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+                xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+                xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
+                xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+                xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+                xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
+        }
+        else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+
+                /*Portion of the segments*/
+                s1=gl[0]/(gl[0]-gl[1]);
+                s2=gl[0]/(gl[0]-gl[2]);
+                if(gl[0]<0.){
+                        pt1 = 1; pt2 = 0;
+                }
+
+                xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
+                xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
+                xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
+                xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
+                xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
+                xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
+        }
+        else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+
+                /*Portion of the segments*/
+                s1=gl[1]/(gl[1]-gl[0]);
+                s2=gl[1]/(gl[1]-gl[2]);
+                if(gl[1]<0.){
+                        pt1 = 1; pt2 = 0;
+                }
+
+                xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
+                xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
+                xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
+                xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
+                xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
+                xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
+        }
+        else{
+                _error_("case not possible");
+        }
+
+
+        /*Some checks in debugging mode*/
+        _assert_(s1>=0 && s1<=1.);
+        _assert_(s2>=0 && s2<=1.);
+
+        /*Get normal vector*/
+        IssmDouble normal[3];
+        this->NormalSectionBase(&normal[0],&xyz_front[0][0]);
+        normal[0] = -normal[0];
+        normal[1] = -normal[1];
+
+        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+        /*Get inputs*/
+        IssmDouble flux = 0.;
+        IssmDouble vx,vy,thickness,Jdet;
+        IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
+        Input2 *thickness_input = this->GetInput2(ThicknessEnum); _assert_(thickness_input);
+        Input2 *vx_input        = this->GetInput2(VxAverageEnum); _assert_(vx_input);
+        Input2 *vy_input        = this->GetInput2(VyAverageEnum); _assert_(vy_input);
+
+        /*Start looping on Gaussian points*/
+        Gauss* gauss=this->NewGaussBase(xyz_list,&xyz_front[0][0],3);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+                thickness_input->GetInputValue(&thickness,gauss);
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                this->JacobianDeterminantLine(&Jdet,&xyz_front[0][0],gauss);
+
+                flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
+        }
+
+        /*Clean up and return*/
+        delete gauss;
+        return flux;
+
 }
 /*}}}*/

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

@@ -104 +104 @@
                 void           GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
                 IssmDouble     GroundedArea(bool scaled);
+                IssmDouble     GroundinglineMassFlux(bool scaled);
                 IssmDouble              IcefrontMassFluxLevelset(bool scaled);
                 IssmDouble     IceVolume(bool scaled);

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

@@ -2872 +2872 @@
 IssmDouble Tria::GroundinglineMassFlux(bool scaled){/*{{{*/
 
-        /*Make sure there is an ice front here*/
+        /*Make sure there is a grounding line here*/
         if(!IsIceInElement()) return 0;
         if(!IsZeroLevelset(MaskGroundediceLevelsetEnum)) return 0;

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

@@ -1620 +1620 @@
 }/*}}}*/
 void FemModel::GroundinglineMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/
+
+        /*First we need to depth average the velocities*/
+        int domaintype;
+        this->parameters->FindParam(&domaintype,DomainTypeEnum);
+        this->parameters->SetParam(VxEnum,InputToDepthaverageInEnum);
+        this->parameters->SetParam(VxAverageEnum,InputToDepthaverageOutEnum);
+        depthaverage_core(this);
+        if(domaintype!=Domain2DverticalEnum){
+                this->parameters->SetParam(VyEnum,InputToDepthaverageInEnum);
+                this->parameters->SetParam(VyAverageEnum,InputToDepthaverageOutEnum);
+                depthaverage_core(this);
+        }
 
         IssmDouble local_mass_flux = 0;