source: issm/oecreview/Archive/16554-17801/ISSM-17193-17194.diff@ 17802

Index: ../trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.h       (revision 17193)
+++ ../trunk-jpl/src/c/classes/Elements/Element.h       (revision 17194)
@@ -220,6 +220,8 @@
                virtual int    PressureInterpolation()=0;
                virtual bool   IsZeroLevelset(int levelset_enum)=0;
                virtual void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
+               virtual void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;
+               virtual void   GetNormalFromLSF(IssmDouble *pnormal)=0;
 
                #ifdef _HAVE_RESPONSES_
                virtual void   AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp        (revision 17193)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp        (revision 17194)
@@ -890,6 +890,69 @@
        *pxyz_zero= xyz_zero;
 }
 /*}}}*/
+void  Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
+       
+       /* Intermediaries */
+       int i, dir,nrfrontnodes;
+       IssmDouble  levelset[NUMVERTICES];
+
+       /*Recover parameters and values*/
+       GetInputListOnVertices(&levelset[0],levelsetenum);
+
+       int* indicesfront = xNew<int>(NUMVERTICES);
+       /* Get nodes where there is no ice */
+       nrfrontnodes=0;
+       for(i=0;i<NUMVERTICES;i++){
+               if(levelset[i]>=0.){
+                       indicesfront[nrfrontnodes]=i;
+                       nrfrontnodes++;
+               }
+       }
+
+       IssmDouble* xyz_front = xNew<IssmDouble>(3*nrfrontnodes);
+       /* Return nodes */
+       for(i=0;i<nrfrontnodes;i++){
+               for(dir=0;dir<3;dir++){
+                       xyz_front[3*i+dir]=xyz_list[indicesfront[i]+dir];
+               }
+       }
+
+       *pxyz_front=xyz_front;
+
+       xDelete<int>(indicesfront);
+}/*}}}*/
+void  Tria::GetNormalFromLSF(IssmDouble *pnormal){/*{{{*/
+
+       /* Intermediaries */
+       const int dim=2;
+       int i,counter;
+       IssmDouble* xyz_list = NULL;
+       IssmDouble  dlevelset[dim], norm_dlevelset;
+       IssmDouble  normal[dim]={0.};
+
+       /*Retrieve all inputs and parameters*/
+       Input* levelset_input=this->GetInput(MaskIceLevelsetEnum); _assert_(levelset_input);
+       this->GetVerticesCoordinates(&xyz_list);
+       
+       counter=0;
+       Gauss* gauss = this->NewGauss(2);
+       for(int ig=gauss->begin();ig<gauss->end();ig++){
+               /* Get normal on node */
+               levelset_input->GetInputDerivativeValue(&dlevelset[0],xyz_list,gauss);
+               norm_dlevelset=0.;
+               for(i=0;i<dim;i++) norm_dlevelset+=dlevelset[i]*dlevelset[i]; 
+               norm_dlevelset=sqrt(norm_dlevelset)+1.e-14;
+               for(i=0;i<dim;i++) normal[i]+=dlevelset[i]/norm_dlevelset;
+               counter++;
+       }
+       _assert_(counter>0);
+       for(i=0;i<dim;i++) normal[i]/counter;
+       
+       pnormal=&normal[0];
+
+       delete gauss;
+       xDelete<IssmDouble>(xyz_list);
+}/*}}}*/
 /*FUNCTION Tria::GetNodeIndex {{{*/
 int Tria::GetNodeIndex(Node* node){
 
Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.h  (revision 17193)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.h  (revision 17194)
@@ -130,6 +130,8 @@
                int         VertexConnectivity(int vertexindex);
                void   VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
+               void        GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
+           void        GetNormalFromLSF(IssmDouble *pnormal);
                bool        IsZeroLevelset(int levelset_enum);
 
                #ifdef _HAVE_RESPONSES_
Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp       (revision 17193)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp       (revision 17194)
@@ -3165,7 +3165,7 @@
        /*Retrieve all inputs and parameters*/
        GetInputListOnVertices(&ls[0],levelset_enum);
 
-       /*If the level set is awlays <=0, there is no ice front here*/
+       /*If the level set has always same sign, there is no ice front here*/
        iszerols = false;
        if(IsIceInElement()){
                if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
Index: ../trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17193)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17194)
@@ -114,6 +114,8 @@
                int    PressureInterpolation();
                bool   IsZeroLevelset(int levelset_enum);
                void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
+               void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
+               void   GetNormalFromLSF(IssmDouble *pnormal){_error_("not implemented yet");};
                void   PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
                void   ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
                void   ResetCoordinateSystem(void);
Index: ../trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Seg.h   (revision 17193)
+++ ../trunk-jpl/src/c/classes/Elements/Seg.h   (revision 17194)
@@ -159,6 +159,8 @@
                void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
                bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
+               void            GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};
+               void        GetNormalFromLSF(IssmDouble *pnormal){_error_("not implemented yet");};
 
                #ifdef _HAVE_HYDROLOGY_
                void    GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode){_error_("not implemented yet");};