Changeset 17210 for issm/trunk-jpl/src

Timestamp:

02/04/14 16:04:28 (11 years ago)

Author:

jbondzio

Message:

revert to r17207

Location:

issm/trunk-jpl/src/c

Files:

• analyses/EnthalpyAnalysis.cpp (modified) (1 diff)
• analyses/LevelsetAnalysis.cpp (modified) (2 diffs)
• analyses/LevelsetAnalysis.h (modified) (1 diff)
• analyses/StressbalanceAnalysis.cpp (modified) (5 diffs)
• classes/Elements/Element.h (modified) (1 diff)
• classes/Elements/Penta.h (modified) (1 diff)
• classes/Elements/Seg.h (modified) (1 diff)
• classes/Elements/Tria.cpp (modified) (7 diffs)
• classes/Elements/Tria.h (modified) (1 diff)
• classes/FemModel.cpp (modified) (1 diff)
• cores/transient_core.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp ¶

@@ -826 +826 @@
         int solution_type, i;
         bool computebasalmeltingrates=true;
-        bool isdrainage=true;
+        bool isdrainage=false;
         bool updatebasalconstraints=true;
 

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp ¶

@@ -67 +67 @@
         bool save_results;
 
+        /* extrapolate */
+        Analysis* analysis = new ExtrapolationAnalysis();
+        femmodel->parameters->SetParam(VxEnum,ExtrapolationVariableEnum);
+        analysis->Core(femmodel);
+        delete analysis;
+
         /*activate formulation: */
         femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
 
+        /*recover parameters: */
+        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
+
         if(VerboseSolution()) _printf0_("call computational core:\n");
         solutionsequence_linear(femmodel);
 
-        /*recover parameters: */
-        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
         if(save_results){
                 if(VerboseSolution()) _printf0_("   saving results\n");
@@ -317 +324 @@
 }/*}}}*/
 
-/* Update of constraints */
-void LevelsetAnalysis::UpdateNoIceConstraints(FemModel* femmodel){/*{{{*/
-
-        IssmDouble* mask_ice = GetMaskOfIce(femmodel->elements, femmodel->nodes); 
-
-        for (int i=0;i<femmodel->nodes->Size();i++){
-                Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
-                if(node->InAnalysis(LevelsetAnalysisEnum)){
-                        if(mask_ice[node->Sid()]==1.){//FIXME: what should be done with actual spcs to ice model?
-//                              node->DofInFSet(0); /*remove spc*/ 
-                        }
-                        else{
-                                IssmDouble defval=0.;
-                                node->ApplyConstraint(1,defval); /*apply spc*/ 
-                        }
-                }
-        }
-}/*}}}*/
-IssmDouble* LevelsetAnalysis::GetMaskOfIce(Elements* elements, Nodes* nodes){/*{{{*/
-
-        int                 i;
-        IssmDouble*         mask_ice      = NULL;
-        Vector<IssmDouble>* vec_mask_ice  = NULL;
-        Element*            element       = NULL;
-
-        /*Initialize vector with number of vertices*/
-        IssmDouble numnodes=nodes->NumberOfNodes(LevelsetAnalysisEnum);
-        vec_mask_ice=new Vector<IssmDouble>(numnodes); //nodes at ice front that have ice at next time step
-        for(i=0;i<numnodes;i++)
-                vec_mask_ice[i]=0.;
-        /*Fill vector vertices that have no contact to ice: */
-        for(i=0;i<elements->Size();i++){
-                element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                SetMaskOfIceElement(vec_mask_ice, element);
-        }
-
-        /*Assemble vector and serialize */
-        vec_mask_ice->Assemble();
-        mask_ice=vec_mask_ice->ToMPISerial();
-        delete vec_mask_ice;
-        return mask_ice;
-
-}/*}}}*/
-void LevelsetAnalysis::SetMaskOfIceElement(Vector<IssmDouble>* vec_mask_ice, Element* element){/*{{{*/
-
-        /* Intermediaries */
-        int numnodes = element->GetNumberOfNodes();
-        
-        if(element->IsIceInElement()){
-                for(int i = 0;i<numnodes;i++){
-                        vec_mask_ice->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
-                }
-        }
-}/*}}}*/

issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h ¶

@@ -31 +31 @@
         void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 
-        /* Updating constraints */
-        void UpdateNoIceConstraints(FemModel* femmodel);
-        IssmDouble* GetMaskOfIce(Elements* elements, Nodes* nodes);
-        void SetMaskOfIceElement(Vector<IssmDouble>* vec_mask_ice, Element* element);
 };
 #endif

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp ¶

@@ -1368 +1368 @@
 
         /*If no front, return NULL*/
-        if(!element->IsIcefront()) return NULL;
+        if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
 
         /*Intermediaries*/
@@ -1391 +1391 @@
         IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
         element->GetVerticesCoordinates(&xyz_list);
-
-//      element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
-        element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
-
+        element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
         element->NormalSection(&normal[0],xyz_list_front);
-//      element->GetNormalFromLSF(&normal[0]);
 
         /*Start looping on Gaussian points*/
@@ -1785 +1781 @@
 
         /*If no front, return NULL*/
-        if(!element->IsIcefront()) return NULL;
+        if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
 
         /*Intermediaries*/
@@ -2208 +2204 @@
 
         /*If no front, return NULL*/
-        if(!element->IsIcefront()) return NULL;
+        if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
 
         /*Intermediaries*/
@@ -2974 +2970 @@
 
         /*If no front, return NULL*/
-        if(!element->IsIcefront()) return NULL;
+        if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
 
         /*Intermediaries*/

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

@@ -220 +220 @@
                 virtual int    PressureInterpolation()=0;
                 virtual bool   IsZeroLevelset(int levelset_enum)=0;
-                virtual bool   IsIcefront(void)=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;

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

@@ -114 +114 @@
                 int    PressureInterpolation();
                 bool   IsZeroLevelset(int levelset_enum);
-                bool   IsIcefront(void){_error_("not implemented yet");};
                 void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
                 void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};

issm/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");};
-                bool            IsIcefront(void){_error_("not implemented yet");};
                 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");};

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

@@ -909 +909 @@
                 }
         }
-        _assert_(nrfrontnodes==2);
-        
-        /* arrange order of frontnodes such that they are oriented counterclockwise */
-        if((NUMVERTICES+indicesfront[0]-indicesfront[1])%NUMVERTICES!=2){
-                int index=indicesfront[0];
-                indicesfront[0]=indicesfront[1];
-                indicesfront[1]=index;
-        }
+
         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[3*indicesfront[i]+dir];
-                }
-        }
-
-//      for(i=0;i<nrfrontnodes;i++){
-//              _printf0_("coords frontnode " << i << " :[");
-//              for(dir=0;dir<3;dir++){
-//                      xyz_front[3*i+dir]=xyz_list[indicesfront[i]+dir];
-//                      _printf0_(xyz_front[3*i+dir] << "; ");
-//              }
-//              _printf0_("]\n");
-//      }
+                        xyz_front[3*i+dir]=xyz_list[indicesfront[i]+dir];
+                }
+        }
 
         *pxyz_front=xyz_front;
@@ -938 +922 @@
         xDelete<int>(indicesfront);
 }/*}}}*/
-void  Tria::GetNormalFromLSF(IssmDouble *normal){/*{{{*/
+void  Tria::GetNormalFromLSF(IssmDouble *pnormal){/*{{{*/
 
         /* Intermediaries */
@@ -945 +929 @@
         IssmDouble* xyz_list = NULL;
         IssmDouble  dlevelset[dim], norm_dlevelset;
+        IssmDouble  normal[dim]={0.};
 
         /*Retrieve all inputs and parameters*/
@@ -951 +936 @@
         
         counter=0;
-        for(i=0;i<dim;i++) normal[i]=0.;
-
         Gauss* gauss = this->NewGauss(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -964 +947 @@
         }
         _assert_(counter>0);
-        for(i=0;i<dim;i++) normal[i]/=IssmDouble(counter);
+        for(i=0;i<dim;i++) normal[i]/counter;
+        
+        pnormal=&normal[0];
 
         delete gauss;
@@ -2601 +2586 @@
         GetInputListOnVertices(&ls[0],levelset_enum);
 
-        /* If levelset function changes sign, or is zero on at least two vertices, there is a zero level set here */
+        /*If the level set is awlays <0, there is no ice front here*/
         iszerols= false;
         if(IsIceInElement()){
@@ -2612 +2597 @@
 }
 /*}}}*/
-/*FUNCTION Tria::IsIcefront{{{*/
-bool Tria::IsIcefront(void){
-
-        bool isicefront;
-        int i,nrice;
-        IssmDouble ls[NUMVERTICES];
-
-        /*Retrieve all inputs and parameters*/
-        GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
-
-        /* If only one vertex has ice, there is an ice front here */
-        isicefront= false;
-        if(IsIceInElement()){
-                nrice=0;        
-                for(i=0;i<NUMVERTICES;i++)
-                        if(ls[i]<0.) nrice++;
-                if(nrice==1) isicefront= true;
-        }
-        return isicefront;
-}
-/*}}}*/
-
 
 #ifdef _HAVE_RESPONSES_

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

@@ -132 +132 @@
                 void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
                 void        GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
-                void        GetNormalFromLSF(IssmDouble *normal);
+            void        GetNormalFromLSF(IssmDouble *pnormal);
                 bool        IsZeroLevelset(int levelset_enum);
-                bool        IsIcefront(void);
 
                 #ifdef _HAVE_RESPONSES_

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

@@ -407 +407 @@
         if(VerboseModule()) _printf0_("   Updating constraints for time: " << time << "\n");
 
+        // analysis->UpdateConstraints();
+        
         /*Second, constraints might be time dependent: */
         SpcNodesx(nodes,constraints,parameters,analysis_type); 

issm/trunk-jpl/src/c/cores/transient_core.cpp ¶

@@ -162 +162 @@
                 if(islevelset){
                         if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
-
-                        /* extrapolate */
-                        Analysis* extanalysis = new ExtrapolationAnalysis();
-                        int vars[2] = {VxEnum, VyEnum};
-                        for(int iv=0; i<2;i++){
-                                femmodel->parameters->SetParam(vars[i],ExtrapolationVariableEnum);
-                                extanalysis->Core(femmodel);
-                        }
-                        delete extanalysis;
-
-                        LevelsetAnalysis* lsanalysis = new LevelsetAnalysis();
-                        /* solve level-set equation */
-                        lsanalysis->Core(femmodel);
-                        /* update spc boundary conditions for ice model */
-//                      lsanalysis->UpdateNoIceConstraints(femmodel);
-                        delete lsanalysis;
+                        analysis = new LevelsetAnalysis();
+                        analysis->Core(femmodel);
+                        delete analysis;
                 }