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

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

@@ -220 +220 @@
                 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;

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

@@ -890 +890 @@
         *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){
 

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

@@ -130 +130 @@
                 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_

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

@@ -3165 +3165 @@
         /*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.)){

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

@@ -114 +114 @@
                 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);

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

@@ -159 +159 @@
                 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");};