Changeset 19105 for issm/trunk/src/c/analyses/LsfReinitializationAnalysis.cpp

Timestamp:

02/12/15 16:48:40 (10 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 19103

Location:

issm/trunk

Files:

• . (modified) (2 props)
• src (modified) (1 prop)
• src/c (modified) (1 prop)
• src/c/analyses/LsfReinitializationAnalysis.cpp (modified) (8 diffs)

issm/trunk

@@ +1 @@
+build-fw
+build-ad
 nightlylog
 configure.sh

@@ +1 @@
+build-fw
+build-ad
 nightlylog
 configure.sh

issm/trunk/src/c

@@ -5 +5 @@
 *.obj
 *.exe
+issm
+kriging
 appscan.*
 ouncemake*
@@ -14 +16 @@
 probe.results
 stXXXX*
-.deps
-.dirstamp
+*.deps
+*.dirstamp
 .libs

@@ -5 +5 @@
 *.obj
 *.exe
+issm
+kriging
 appscan.*
 ouncemake*
@@ -14 +16 @@
 probe.results
 stXXXX*
-.deps
-.dirstamp
+*.deps
+*.dirstamp
 .libs

issm/trunk/src/c/analyses/LsfReinitializationAnalysis.cpp

@@ -10 +10 @@
 
 /*Model processing*/
+void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+        /* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+        /* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+        int finiteelement=P1Enum;
+        if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+        ::CreateNodes(nodes,iomodel,LsfReinitializationAnalysisEnum,finiteelement);
+        iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+}/*}}}*/
 int  LsfReinitializationAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
         return 1;
-}/*}}}*/
-void LsfReinitializationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-        /* Do nothing for now */
 }/*}}}*/
 void LsfReinitializationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
@@ -34 +43 @@
         iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
 }/*}}}*/
-void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-        int finiteelement=P1Enum;
-        if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-        ::CreateNodes(nodes,iomodel,LsfReinitializationAnalysisEnum,finiteelement);
-        iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-}/*}}}*/
-void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+void LsfReinitializationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
         /* Do nothing for now */
 }/*}}}*/
-void LsfReinitializationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-        /* Do nothing for now */
-}/*}}}*/
 
 /*Finite element Analysis*/
-void  LsfReinitializationAnalysis::Core(FemModel* femmodel){/*{{{*/
+void           LsfReinitializationAnalysis::Core(FemModel* femmodel){/*{{{*/
 
         /*parameters: */
@@ -249 +249 @@
         return pe;
         }/*}}}*/
-void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           LsfReinitializationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+         * For node i, Bi can be expressed in the actual coordinate system
+         * by: 
+         *       Bi=[ N ]
+         *          [ N ]
+         * where N is the finiteelement function for node i.
+         *
+         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions*/
+        IssmDouble* basis=xNew<IssmDouble>(numnodes);
+        element->NodalFunctions(basis,gauss);
+
+        /*Build B: */
+        for(int i=0;i<numnodes;i++){
+                B[numnodes*0+i] = basis[i];
+                B[numnodes*1+i] = basis[i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(basis);
+}/*}}}*/
+void           LsfReinitializationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+         * For node i, Bi' can be expressed in the actual coordinate system
+         * by: 
+         *       Bi_prime=[ dN/dx ]
+         *                [ dN/dy ]
+         * where N is the finiteelement function for node i.
+         *
+         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+        /*Build B': */
+        for(int i=0;i<numnodes;i++){
+                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
+                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
+
+}/*}}}*/
+IssmDouble     LsfReinitializationAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
+        // returns distance d of point q to straight going through points s0, s1
+        // d=|a x b|/|b|
+        // with a=q-s0, b=s1-s0
+        
+        /* Intermediaries */
+        const int dim=2;
+        int i;
+        IssmDouble a[dim], b[dim];
+        IssmDouble norm_b;
+
+        for(i=0;i<dim;i++){
+                a[i]=q[i]-s0[i];
+                b[i]=s1[i]-s0[i];
+        }
+        
+        norm_b=0.;
+        for(i=0;i<dim;i++)
+                norm_b+=b[i]*b[i];
+        norm_b=sqrt(norm_b);
+        _assert_(norm_b>0.);
+
+        return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
+}/*}}}*/
+void           LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
         IssmDouble   lsf;
@@ -284 +363 @@
 
 }/*}}}*/
-void LsfReinitializationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           LsfReinitializationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
         _error_("Not implemented yet");
 }/*}}}*/
-void LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
         int domaintype;
@@ -301 +380 @@
         }
 }/*}}}*/
-void LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-        /* Do nothing for now */
-}/*}}}*/
-void LsfReinitializationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-         * For node i, Bi can be expressed in the actual coordinate system
-         * by: 
-         *       Bi=[ N ]
-         *          [ N ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions*/
-        IssmDouble* basis=xNew<IssmDouble>(numnodes);
-        element->NodalFunctions(basis,gauss);
-
-        /*Build B: */
-        for(int i=0;i<numnodes;i++){
-                B[numnodes*0+i] = basis[i];
-                B[numnodes*1+i] = basis[i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(basis);
-}/*}}}*/
-void LsfReinitializationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-         * For node i, Bi' can be expressed in the actual coordinate system
-         * by: 
-         *       Bi_prime=[ dN/dx ]
-         *                [ dN/dy ]
-         * where N is the finiteelement function for node i.
-         *
-         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
-         */
-
-        /*Fetch number of nodes for this finite element*/
-        int numnodes = element->GetNumberOfNodes();
-
-        /*Get nodal functions derivatives*/
-        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-        /*Build B': */
-        for(int i=0;i<numnodes;i++){
-                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
-        }
-
-        /*Clean-up*/
-        xDelete<IssmDouble>(dbasis);
-
-}/*}}}*/
-
-/* Other */
-void LsfReinitializationAnalysis::SetReinitSPCs(FemModel* femmodel){/*{{{*/
+bool           LsfReinitializationAnalysis::ReinitConvergence(Vector<IssmDouble>* lsfg,Vector<IssmDouble>* lsfg_old,IssmDouble reltol){/*{{{*/
+
+        /*Output*/
+        bool converged = true;
+
+        /*Intermediary*/
+        Vector<IssmDouble>* dlsfg    = NULL;
+        IssmDouble          norm_dlsf,norm_lsf;
+
+        /*compute norm(du)/norm(u)*/
+        dlsfg=lsfg_old->Duplicate(); lsfg_old->Copy(dlsfg); dlsfg->AYPX(lsfg,-1.0);
+        norm_dlsf=dlsfg->Norm(NORM_TWO); norm_lsf=lsfg_old->Norm(NORM_TWO);
+        if (xIsNan<IssmDouble>(norm_dlsf) || xIsNan<IssmDouble>(norm_lsf)) _error_("convergence criterion is NaN!");
+        if((norm_dlsf/norm_lsf)<reltol){
+                if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" < "<<reltol*100<<" %\n");
+        }
+        else{ 
+                if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" > "<<reltol*100<<" %\n");
+                converged=false;
+        }
+
+        /*Cleanup*/
+        delete dlsfg;
+
+        return converged;
+}/*}}}*/
+void           LsfReinitializationAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
+
+        /* Intermediaries */
+        int i;
+        IssmDouble dmaxp,dmaxm,val;
+        Element* element;
+
+        /*Initialize vector with number of vertices*/
+        int numvertices=femmodel->vertices->NumberOfVertices();
+
+        Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
+        GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
+        
+        /* set distance on elements intersected by zero levelset */
+        for(i=0;i<femmodel->elements->Size();i++){
+                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
+                        SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
+                }
+        }
+        vec_dist_zerolevelset->Assemble();
+
+        /* Get maximum distance to interface along vertices */
+        dmaxp=0.; dmaxm=0.;
+        for(i=0;i<numvertices;i++){
+                vec_dist_zerolevelset->GetValue(&val,i); 
+                if((val>0.) && (val>dmaxp))
+                         dmaxp=val;
+                else if((val<0.) && (val<dmaxm))
+                         dmaxm=val;
+        }
+        //wait until all values are computed
+
+        /* set all none intersected vertices to max/min distance */
+        for(i=0;i<numvertices;i++){
+                vec_dist_zerolevelset->GetValue(&val,i);
+                if(val==1.) //FIXME: improve check
+                        vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL);
+                else if(val==-1.)
+                        vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL);
+        }
+
+        /*Assemble vector and serialize */
+        vec_dist_zerolevelset->Assemble();
+        IssmDouble* dist_zerolevelset=vec_dist_zerolevelset->ToMPISerial();
+        InputUpdateFromVectorx(femmodel,dist_zerolevelset,MaskIceLevelsetEnum,VertexSIdEnum);
+
+        /*Clean up and return*/
+        delete vec_dist_zerolevelset;
+        delete dist_zerolevelset;
+
+}/*}}}*/
+void           LsfReinitializationAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
+
+        if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
+                return;
+
+        /* Intermediaries */
+        int dim=3;
+        int i,d;
+        IssmDouble dist,lsf_old;
+
+        int numvertices=element->GetNumberOfVertices();
+
+        IssmDouble* lsf = xNew<IssmDouble>(numvertices);
+        IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
+        IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
+        IssmDouble* s0 = xNew<IssmDouble>(dim);
+        IssmDouble* s1 = xNew<IssmDouble>(dim);
+        IssmDouble* v = xNew<IssmDouble>(dim);
+        IssmDouble* xyz_list = NULL;
+        IssmDouble* xyz_list_zero = NULL;
+
+        /* retrieve inputs and parameters */
+        element->GetVerticesCoordinates(&xyz_list);
+        element->GetInputListOnVertices(lsf,MaskIceLevelsetEnum);
+
+        /* get sign of levelset function */
+        for(i=0;i<numvertices;i++)
+                sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
+
+        element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
+        for(i=0;i<dim;i++){
+                s0[i]=xyz_list_zero[0+i];
+                s1[i]=xyz_list_zero[3+i];
+        }
+
+        /* get signed_distance of vertices to zero levelset straight */
+        for(i=0;i<numvertices;i++){
+                for(d=0;d<dim;d++)
+                        v[d]=xyz_list[dim*i+d];
+                dist=GetDistanceToStraight(&v[0],&s0[0],&s1[0]);
+                signed_dist[i]=sign_lsf[i]*dist;
+        }
+        
+        /* insert signed_distance into vec_signed_dist, if computed distance is smaller */
+        for(i=0;i<numvertices;i++){
+                vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
+                /* initial lsf values are +-1. Update those fields or if distance to interface smaller.*/
+                if(fabs(lsf_old)==1. || fabs(signed_dist[i])<fabs(lsf_old))
+                        vec_signed_dist->SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL);
+        }
+
+        xDelete<IssmDouble>(lsf);
+        xDelete<IssmDouble>(sign_lsf);
+        xDelete<IssmDouble>(signed_dist);
+        xDelete<IssmDouble>(s0);
+        xDelete<IssmDouble>(s1);
+        xDelete<IssmDouble>(v);
+
+}/*}}}*/
+void           LsfReinitializationAnalysis::SetReinitSPCs(FemModel* femmodel){/*{{{*/
 
         int i,k, numnodes;
@@ -369 +525 @@
         /* deactivate all spcs */
         for(i=0;i<femmodel->elements->Size();i++){
-                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
                 for(k=0;k<element->GetNumberOfNodes();k++){
                         node=element->GetNode(k);
@@ -382 +538 @@
         /* reactivate spcs on elements intersected by zero levelset */
         for(i=0;i<femmodel->elements->Size();i++){
-                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
                 if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
                         /*iterate over nodes and set spc */
@@ -399 +555 @@
 
 }/*}}}*/
-void LsfReinitializationAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
-
-        /* Intermediaries */
-        int i;
-        IssmDouble dmaxp,dmaxm,val;
-        Element* element;
-
-        /*Initialize vector with number of vertices*/
-        int numvertices=femmodel->vertices->NumberOfVertices();
-
-        Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
-        GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
-        
-        /* set distance on elements intersected by zero levelset */
-        for(i=0;i<femmodel->elements->Size();i++){
-                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
-                        SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
-                }
-        }
-        vec_dist_zerolevelset->Assemble();
-
-        /* Get maximum distance to interface along vertices */
-        dmaxp=0.; dmaxm=0.;
-        for(i=0;i<numvertices;i++){
-                vec_dist_zerolevelset->GetValue(&val,i); 
-                if((val>0.) && (val>dmaxp))
-                         dmaxp=val;
-                else if((val<0.) && (val<dmaxm))
-                         dmaxm=val;
-        }
-        //wait until all values are computed
-
-        /* set all none intersected vertices to max/min distance */
-        for(i=0;i<numvertices;i++){
-                vec_dist_zerolevelset->GetValue(&val,i);
-                if(val==1.) //FIXME: improve check
-                        vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL);
-                else if(val==-1.)
-                        vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL);
-        }
-
-        /*Assemble vector and serialize */
-        vec_dist_zerolevelset->Assemble();
-        IssmDouble* dist_zerolevelset=vec_dist_zerolevelset->ToMPISerial();
-        InputUpdateFromVectorx(femmodel,dist_zerolevelset,MaskIceLevelsetEnum,VertexSIdEnum);
-
-        /*Clean up and return*/
-        delete vec_dist_zerolevelset;
-        delete dist_zerolevelset;
-
-}/*}}}*/
-void LsfReinitializationAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
-
-        if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
-                return;
-
-        /* Intermediaries */
-        int dim=3;
-        int i,d;
-        IssmDouble dist,lsf_old;
-
-        int numvertices=element->GetNumberOfVertices();
-
-        IssmDouble* lsf = xNew<IssmDouble>(numvertices);
-        IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
-        IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
-        IssmDouble* s0 = xNew<IssmDouble>(dim);
-        IssmDouble* s1 = xNew<IssmDouble>(dim);
-        IssmDouble* v = xNew<IssmDouble>(dim);
-        IssmDouble* xyz_list = NULL;
-        IssmDouble* xyz_list_zero = NULL;
-
-        /* retrieve inputs and parameters */
-        element->GetVerticesCoordinates(&xyz_list);
-        element->GetInputListOnVertices(lsf,MaskIceLevelsetEnum);
-
-        /* get sign of levelset function */
-        for(i=0;i<numvertices;i++)
-                sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
-
-        element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
-        for(i=0;i<dim;i++){
-                s0[i]=xyz_list_zero[0+i];
-                s1[i]=xyz_list_zero[3+i];
-        }
-
-        /* get signed_distance of vertices to zero levelset straight */
-        for(i=0;i<numvertices;i++){
-                for(d=0;d<dim;d++)
-                        v[d]=xyz_list[dim*i+d];
-                dist=GetDistanceToStraight(&v[0],&s0[0],&s1[0]);
-                signed_dist[i]=sign_lsf[i]*dist;
-        }
-        
-        /* insert signed_distance into vec_signed_dist, if computed distance is smaller */
-        for(i=0;i<numvertices;i++){
-                vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
-                /* initial lsf values are +-1. Update those fields or if distance to interface smaller.*/
-                if(fabs(lsf_old)==1. || fabs(signed_dist[i])<fabs(lsf_old))
-                        vec_signed_dist->SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL);
-        }
-
-        xDelete<IssmDouble>(lsf);
-        xDelete<IssmDouble>(sign_lsf);
-        xDelete<IssmDouble>(signed_dist);
-        xDelete<IssmDouble>(s0);
-        xDelete<IssmDouble>(s1);
-        xDelete<IssmDouble>(v);
-
-}/*}}}*/
-IssmDouble LsfReinitializationAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
-        // returns distance d of point q to straight going through points s0, s1
-        // d=|a x b|/|b|
-        // with a=q-s0, b=s1-s0
-        
-        /* Intermediaries */
-        const int dim=2;
-        int i;
-        IssmDouble a[dim], b[dim];
-        IssmDouble norm_b;
-
-        for(i=0;i<dim;i++){
-                a[i]=q[i]-s0[i];
-                b[i]=s1[i]-s0[i];
-        }
-        
-        norm_b=0.;
-        for(i=0;i<dim;i++)
-                norm_b+=b[i]*b[i];
-        norm_b=sqrt(norm_b);
-        _assert_(norm_b>0.);
-
-        return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
-}/*}}}*/
-bool LsfReinitializationAnalysis::ReinitConvergence(Vector<IssmDouble>* lsfg,Vector<IssmDouble>* lsfg_old,IssmDouble reltol){/*{{{*/
-
-        /*Output*/
-        bool converged = true;
-
-        /*Intermediary*/
-        Vector<IssmDouble>* dlsfg    = NULL;
-        IssmDouble          norm_dlsf,norm_lsf;
-
-        /*compute norm(du)/norm(u)*/
-        dlsfg=lsfg_old->Duplicate(); lsfg_old->Copy(dlsfg); dlsfg->AYPX(lsfg,-1.0);
-        norm_dlsf=dlsfg->Norm(NORM_TWO); norm_lsf=lsfg_old->Norm(NORM_TWO);
-        if (xIsNan<IssmDouble>(norm_dlsf) || xIsNan<IssmDouble>(norm_lsf)) _error_("convergence criterion is NaN!");
-        if((norm_dlsf/norm_lsf)<reltol){
-                if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" < "<<reltol*100<<" %\n");
-        }
-        else{ 
-                if(VerboseConvergence()) _printf0_("\n"<<setw(50)<<left<<"   Velocity convergence: norm(du)/norm(u)"<<norm_dlsf/norm_lsf*100<<" > "<<reltol*100<<" %\n");
-                converged=false;
-        }
-
-        /*Cleanup*/
-        delete dlsfg;
-
-        return converged;
-}/*}}}*/
+void           LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+        /* Do nothing for now */
+}/*}}}*/

@@ +1 @@
+build-fw
+build-ad
 nightlylog
 configure.sh

@@ -5 +5 @@
 *.obj
 *.exe
+issm
+kriging
 appscan.*
 ouncemake*
@@ -14 +16 @@
 probe.results
 stXXXX*
-.deps
-.dirstamp
+*.deps
+*.dirstamp
 .libs