source: issm/oecreview/Archive/26740-27031/ISSM-27018-27019.diff@ 27032

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

@@ -651 +651 @@
         IssmDouble newlevelset[6];
         femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
 
-        if(calvinglaw==CalvingMinthicknessEnum){
+        /*Apply minimum thickness criterion*/
+        if(calvinglaw==CalvingMinthicknessEnum || calvinglaw==CalvingVonmisesEnum || calvinglaw==CalvingParameterizationEnum){
 
                 IssmDouble mig_max = femmodel->parameters->FindParam(MigrationMaxEnum);
                 IssmDouble dt      = femmodel->parameters->FindParam(TimesteppingTimeStepEnum);
@@ -690 +691 @@
                                 dis_input->GetInputValue(&distance,gauss);
 
                                 if(thickness<min_thickness && bed<sealevel && fabs(distance)<mig_max*dt && levelset<0){
-                                        newlevelset[in] = +400.; //Arbitrary > 0 number
+                                        newlevelset[in] = +400.; //Arbitrary > 0 number (i.e. deactivate this node)
                                 }
                                 else{
                                         newlevelset[in] = levelset;
@@ -713 +714 @@
    InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
    femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
 
-        if(calvinglaw==CalvingVonmisesEnum){
+   if(calvinglaw==CalvingHabEnum){
 
                 /*Intermediaries*/
-                IssmDouble thickness,bed,sealevel,distance,levelset;
-                IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum);
-
-                /*Loop over all elements of this partition*/
-                for(Object* & object : femmodel->elements->objects){
-                        Element* element  = xDynamicCast<Element*>(object);
-
-                        int      numnodes     = element->GetNumberOfNodes();
-                        Gauss*   gauss        = element->NewGauss();
-
-                        Input *H_input        = element->GetInput(ThicknessEnum);              _assert_(H_input);
-                        Input *b_input        = element->GetInput(BedEnum);                    _assert_(b_input);
-                        Input *sl_input       = element->GetInput(SealevelEnum);               _assert_(sl_input);
-                        Input *dis_input      = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input);
-                        Input *levelset_input = element->GetInput(MaskIceLevelsetEnum);        _assert_(levelset_input);
-
-                        /*Potentially constrain nodes of this element*/
-                        for(int in=0;in<numnodes;in++){
-                                gauss->GaussNode(element->GetElementType(),in);
-                                Node* node=element->GetNode(in);
-                                if(!node->IsActive()) continue;
-
-                                H_input->GetInputValue(&thickness,gauss);
-                                b_input->GetInputValue(&bed,gauss);
-                                sl_input->GetInputValue(&sealevel,gauss);
-                                dis_input->GetInputValue(&distance,gauss);
-                                levelset_input->GetInputValue(&levelset,gauss);
-                                if(thickness<min_thickness && bed<sealevel && fabs(distance)<mig_max*dt && levelset<0){
-                                        node->ApplyConstraint(0,+1.);
-                                }
-                                else{
-                                        /* no ice, set no spc */
-                                        node->DofInFSet(0);
-                                }
-                        }
-                        delete gauss;
-                }
-        }
-        else if(calvinglaw==CalvingParameterizationEnum){
-
-      /*Intermediaries*/
-      IssmDouble thickness,bed,sealevel,distance;
-      IssmDouble min_thickness = femmodel->parameters->FindParam(CalvingMinthicknessEnum);
-
-                /*Loop over all elements of this partition*/
-                for(Object* & object : femmodel->elements->objects){
-                        Element* element  = xDynamicCast<Element*>(object);
-
-                        int      numnodes = element->GetNumberOfNodes();
-                        Gauss*   gauss    = element->NewGauss();
-                        Input*   H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
-                        Input*   b_input = element->GetInput(BedEnum); _assert_(b_input);
-                        Input*   sl_input = element->GetInput(SealevelEnum); _assert_(sl_input);
-                        Input*   dis_input           = element->GetInput(DistanceToCalvingfrontEnum); _assert_(dis_input);
-
-                        /*Potentially constrain nodes of this element*/
-                        for(int in=0;in<numnodes;in++){
-                                gauss->GaussNode(element->GetElementType(),in);
-                                Node* node=element->GetNode(in);
-                                if(!node->IsActive()) continue;
-
-                                H_input->GetInputValue(&thickness,gauss);
-                                b_input->GetInputValue(&bed,gauss);
-                                sl_input->GetInputValue(&sealevel,gauss);
-                                dis_input->GetInputValue(&distance,gauss);
-                                if(thickness<min_thickness && bed<sealevel && fabs(distance)<mig_max*dt){
-                                        node->ApplyConstraint(0,+1.);
-                                }
-                                else {
-                                        /* no ice, set no spc */
-                                        node->DofInFSet(0);
-                                }
-                        }
-                        delete gauss;
-                }
-        }
-   else if(calvinglaw==CalvingHabEnum){
-
-                /*Intermediaries*/
                 IssmDouble  thickness,water_depth,distance,hab_fraction;
 
                 /*Loop over all elements of this partition*/