source: issm/oecreview/Archive/20545-21336/ISSM-20718-20719.diff@ 21337

../trunk-jpl/src/c/modules/Calvingx/Calvingx.cpp

@@ -26 +26 @@
                         femmodel->CalvingRateDevx();
                         femmodel->ElementOperationx(&Element::CalvingRateDev);
                         break;
+                case CalvingMinthicknessEnum:
+                        femmodel->CalvingRateMinthicknessx();
+                        break;
                 default:
                         _error_("Caving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
         }

../trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -2076 +2076 @@
                                 case CalvingDevEnum:
                                         this->CalvingRateDev();
                                         break;
+                                case CalvingMinthicknessEnum:
+                                        this->CalvingRateMinthickness();
+                                        break;
                                 default:
                                         _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
                         }

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

@@ -175 +175 @@
                 virtual void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part)=0;
                 virtual void        CalvingRateLevermann(void)=0;
                 virtual void       CalvingRateDev(void){_error_("not implemented yet");};
+                virtual void       CalvingRateMinthickness(void){_error_("not implemented yet");};
                 virtual void       WriteLevelsetSegment(DataSet* segments){_error_("not implemented yet");};
                 virtual void       ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){_error_("not implemented yet");};
                 virtual IssmDouble CharacteristicLength(void)=0;

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

@@ -337 +337 @@
         delete gauss;
 }
 /*}}}*/
+void       Tria::CalvingRateMinthickness(){/*{{{*/
+
+        IssmDouble  calvingratex[NUMVERTICES];
+        IssmDouble  calvingratey[NUMVERTICES];
+        IssmDouble  calvingrate[NUMVERTICES];
+        IssmDouble  vx,vy,H,minH;
+
+        /*Retrieve all inputs and parameters we will need*/
+        Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
+        Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
+        Input* H_input  = inputs->GetInput(ThicknessEnum); _assert_(H_input);
+        parameters->FindParam(&minH,CalvingMinthicknessEnum);
+
+        /* Start looping on the number of vertices: */
+        GaussTria* gauss=new GaussTria();
+        for(int iv=0;iv<NUMVERTICES;iv++){
+                gauss->GaussVertex(iv);
+
+                /*Get velocity components and thickness*/
+                vx_input->GetInputValue(&vx,gauss);
+                vy_input->GetInputValue(&vy,gauss);
+                H_input->GetInputValue(&H,gauss);
+
+                /*Assign values*/
+                if(H>minH){
+                        calvingratex[iv]=0.;
+                        calvingratey[iv]=0.;
+                        calvingrate[iv]=0.;
+                }
+                else{
+                        calvingratex[iv]=vx+1e-2; //counter balance advance + add some retreat
+                        calvingratey[iv]=vy+1e-2;
+                        calvingrate[iv]=sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
+                }
+        }
+
+        /*Add input*/
+        this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+        this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+        /*Clean up and return*/
+        delete gauss;
+}
+/*}}}*/
 void       Tria::WriteLevelsetSegment(DataSet* segments){/*{{{*/
 
         if(!this->IsZeroLevelset(MaskIceLevelsetEnum)) return;

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

@@ -53 +53 @@
                 void        AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
                 void                    CalvingRateLevermann();
                 void                    CalvingRateDev();
+                void                    CalvingRateMinthickness();
                 void                    WriteLevelsetSegment(DataSet* segments);
                 void        ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments);
                 IssmDouble  CharacteristicLength(void);

../trunk-jpl/src/c/classes/FemModel.h

@@ -94 +94 @@
                 void DeviatoricStressx();
                 void CalvingRateLevermannx();
                 void CalvingRateDevx();
+                void CalvingRateMinthicknessx();
                 void ResetLevelset();
                 #ifdef  _HAVE_DAKOTA_
                 void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);

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

@@ -2165 +2165 @@
         }
 }
 /*}}}*/
+void FemModel::CalvingRateMinthicknessx(){/*{{{*/
+
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                element->CalvingRateMinthickness();
+        }
+}
+/*}}}*/
 void FemModel::StrainRateparallelx(){/*{{{*/
 
         for(int i=0;i<elements->Size();i++){

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

@@ -197 +197 @@
         switch(calvinglaw){
                 case DefaultCalvingEnum:
                 case CalvingDevEnum:
+                case CalvingMinthicknessEnum:
                         lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
                         if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
                         calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
@@ -250 +251 @@
                 /*Get calving speed*/
                 switch(calvinglaw){
                         case DefaultCalvingEnum:
+                        case CalvingDevEnum:
+                        case CalvingMinthicknessEnum:
                                 lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
                                 if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
                                 calvingrate_input->GetInputValue(&calvingrate,gauss);
@@ -270 +273 @@
                                 break;
 
                         case CalvingLevermannEnum:
-                        case CalvingDevEnum:
                                 calvingratex_input->GetInputValue(&c[0],gauss);
                                 if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
                                 meltingrate_input->GetInputValue(&meltingrate,gauss);