Changeset 23189 for issm/trunk/src/c/analyses/LevelsetAnalysis.cpp

Timestamp:

08/28/18 09:45:51 (7 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 23187

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c/analyses/LevelsetAnalysis.cpp (modified) (20 diffs)

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

@@ -12 +12 @@
 
 void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-        IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,P1Enum);
+        int finiteelement;
+        iomodel->FindConstant(&finiteelement,"md.levelset.fe");
+        IoModelToConstraintsx(constraints,iomodel,"md.levelset.spclevelset",LevelsetAnalysisEnum,finiteelement);
 }
 /*}}}*/
@@ -19 +21 @@
 }/*}}}*/
 void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-        int finiteelement=P1Enum;
+        int finiteelement;
+        iomodel->FindConstant(&finiteelement,"md.levelset.fe");
         if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
         ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
@@ -32 +35 @@
 
         /*Finite element type*/
-        int finiteelement = P1Enum;
+        int finiteelement;
+        iomodel->FindConstant(&finiteelement,"md.levelset.fe");
 
         /*Update elements: */
@@ -43 +47 @@
                 }
         }
-        
+
         iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
         iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
@@ -177 +181 @@
 
         /*Calving threshold*/
-        
+
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes    = basalelement->GetNumberOfNodes();
@@ -326 +330 @@
                                  }
                                 break;
-                                
+
                         case CalvingLevermannEnum:
                                 calvingratex_input->GetInputValue(&c[0],gauss);
@@ -357 +361 @@
                                  }
                                 break;
-                        
+
                         case CalvingHabEnum:
                                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
@@ -378 +382 @@
                                  }
                                 break;
-                        
+
                         case CalvingCrevasseDepthEnum:
                                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
@@ -385 +389 @@
 
                                 if(groundedice<0) meltingrate = 0.;
-                                
+
                                 norm_dlsf=0.;
                                 for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
@@ -516 +520 @@
 }/*}}}*/
 ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
-        
+
         if(!element->IsOnBase()) return NULL;
         Element* basalelement = element->SpawnBasalElement();
@@ -525 +529 @@
         IssmDouble  lsf;
         IssmDouble* xyz_list = NULL;
-        
+
         /*Fetch number of nodes and dof for this finite element*/
         int numnodes = basalelement->GetNumberOfNodes();
@@ -532 +536 @@
         ElementVector* pe = basalelement->NewElementVector();
         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
-        
+
         if(dt!=0.){
                 /*Initialize basis vector*/
@@ -623 +627 @@
         // d=|a x b|/|b|
         // with a=q-s0, b=s1-s0
-        
+
         /* Intermediaries */
         const int dim=2;
@@ -634 +638 @@
                 b[i]=s1[i]-s0[i];
         }
-        
+
         norm_b=0.;
         for(i=0;i<dim;i++)
@@ -662 +666 @@
                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
         }
-}/*}}}*/
-void           LevelsetAnalysis::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++)
-                sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
-
-        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(xIsNan<IssmDouble>(lsf_old) || 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>(xyz_list);
-        xDelete<IssmDouble>(xyz_list_zero);
 }/*}}}*/
 void           LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
@@ -724 +675 @@
         IssmDouble  bed,water_depth;
         IssmDouble  levelset;
-        
+
         femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
 
@@ -756 +707 @@
                 }
         }
-        
+
         if(calvinglaw==CalvingHabEnum){
 
@@ -762 +713 @@
                 femmodel->elements->InputDuplicate(MaskIceLevelsetEnum, DistanceToCalvingfrontEnum);
                 femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
-                
+
                 /*Loop over all elements of this partition*/
                 for(int i=0;i<femmodel->elements->Size();i++){
                         Element* element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                        
+
                         rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
                         rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
@@ -797 +748 @@
                 }
         }
-        
+
         if(calvinglaw==CalvingCrevasseDepthEnum){
-        
+
                 int                 nflipped,local_nflipped;
                 Vector<IssmDouble>* vec_constraint_nodes = NULL;
                 IssmDouble* constraint_nodes = NULL;
-                
+
                 /*Get the DistanceToCalvingfront*/
                 femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
@@ -895 +846 @@
                                 gauss->GaussNode(element->GetElementType(),in);
                                 Node* node=element->GetNode(in);
-                                
+
                                 if(constraint_nodes[node->Sid()]>0.){
                                         node->ApplyConstraint(0,+1.);