Changeset 17307

Timestamp:

02/19/14 13:30:46 (11 years ago)

Author:

jbondzio

Message:

BUG: removing error messages

Location:

issm/trunk-jpl/src/c/analyses

Files:

• LsfReinitializationAnalysis.cpp (modified) (1 diff)
• LsfReinitializationAnalysis.h (modified) (1 diff)

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

@@ -6 +6 @@
 
 /*Model processing*/
-int  LsfReinitializationAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){
-        _error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){
-        _error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){
-        _error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){
-        _error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){
-        _error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){
-        _error_("not implemented yet");
-}
+int  LsfReinitializationAnalysis::DofsPerNode(int** doflist,int meshtype,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){/*{{{*/
+        /* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+        /* Do nothing for now */
+}/*}}}*/
+void LsfReinitializationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+        /* 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){/*{{{*/
         _error_("not implemented yet");
-}
-ElementVector* LsfReinitializationAnalysis::CreateDVector(Element* element){
+}/*}}}*/
+ElementVector* LsfReinitializationAnalysis::CreateDVector(Element* element){/*{{{*/
         _error_("not implemented yet");
-}
-ElementMatrix* LsfReinitializationAnalysis::CreateJacobianMatrix(Element* element){
+}/*}}}*/
+ElementMatrix* LsfReinitializationAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
         _error_("not implemented yet");
-}
-ElementMatrix* LsfReinitializationAnalysis::CreateKMatrix(Element* element){
+}/*}}}*/
+ElementMatrix* LsfReinitializationAnalysis::CreateKMatrix(Element* element){/*{{{*/
         _error_("not implemented yet");
-}
-ElementVector* LsfReinitializationAnalysis::CreatePVector(Element* element){
+}/*}}}*/
+ElementVector* LsfReinitializationAnalysis::CreatePVector(Element* element){/*{{{*/
         _error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){
+}/*}}}*/
+void LsfReinitializationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
         _error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){
+}/*}}}*/
+void LsfReinitializationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
         _error_("not implemented yet");
-}
-void LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){
+}/*}}}*/
+void LsfReinitializationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
         _error_("not implemented yet");
-}
+}/*}}}*/
 
+/* Other */
+void LsfReinitializationAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
+
+        /* Intermediaries */
+        int i;
+
+        /*Initialize vector with number of vertices*/
+        int numvertices=femmodel->vertices->NumberOfVertices();
+        Vector<IssmDouble>* vec_dist_zerolevelset=new Vector<IssmDouble>(numvertices); //vertices that have ice at next time step
+
+        /*TODO: Fill vector with values of old level set function: */
+        for(i=0;i<numvertices;i++){
+                vec_dist_zerolevelset->SetValue(i,0.,INS_VAL);
+        }
+        for(i=0;i<femmodel->elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                if(element->IsZeroLevelset(MaskIceLevelsetEnum))
+                        SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
+        }
+
+        /*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 */
+        const int dim=3;
+        int i,d;
+        int numvertices=element->GetNumberOfVertices();
+        IssmDouble s0[dim], s1[dim], v[dim];
+        IssmDouble dist,lsf_old;
+
+        IssmDouble* lsf = xNew<IssmDouble>(numvertices);
+        IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
+        IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
+        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++){
+                if(lsf[i]==0.)
+                        sign_lsf[i]=1.;
+                else
+                        sign_lsf[i]=lsf[i]/fabs(lsf[i]);
+        }
+
+        element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
+        for(d=0;d<dim;d++){
+                s0[d]=xyz_list_zero[0+d];
+                s1[d]=xyz_list_zero[3+d];
+        }
+
+        /* 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[3*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());
+                if(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);
+}/*}}}*/
+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;
+}/*}}}*/
+

issm/trunk-jpl/src/c/analyses/LsfReinitializationAnalysis.h

@@ -29 +29 @@
         void InputUpdateFromSolution(IssmDouble* solution,Element* element);
         void UpdateConstraints(FemModel* femmodel);
+        
+        /* Other */
+        void SetDistanceOnIntersectedElements(FemModel* femmodel);
+        void SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_dist, Element* element);
+        IssmDouble GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1);
 };
 #endif