Changeset 27035 for issm/trunk/src/c/modules/KillIcebergsx/KillIcebergsx.cpp

Timestamp:

06/01/22 05:01:48 (3 years ago)

Author:

Mathieu Morlighem

Message:

merged trunk-jpl and trunk for revision 27033

Location:

issm/trunk

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/c/modules/KillIcebergsx/KillIcebergsx.cpp (modified) (9 diffs)

issm/trunk/src/c/modules/KillIcebergsx/KillIcebergsx.cpp

@@ -9 +9 @@
 #include "../InputUpdateFromVectorx/InputUpdateFromVectorx.h"
 
-void KillIcebergsx(FemModel* femmodel){
+int KillIcebergsx(FemModel* femmodel){
 
         /*Intermediaries*/
@@ -18 +18 @@
         IssmDouble *local_mask   = xNewZeroInit<IssmDouble>(nbv_local);
         bool       *element_flag = xNewZeroInit<bool>(femmodel->elements->Size());
+        IssmDouble ice;
 
         /*Step 1, go through all elements and put 1 in local_mask if the element is grounded*/
@@ -27 +28 @@
                         element_flag[i] = true;
                 }
-                else{
-                        if(element->IsGrounded()){
+                else {
+                        if(element->IsAllGrounded()){
+                                /* only look at element with ice but not fully grounded */
                                 int numvertices = element->GetNumberOfVertices();
-                                for(int v=0;v<numvertices;v++) local_mask[element->vertices[v]->Lid()] = 1.;
+                                Gauss* gauss=element->NewGauss();
+                                Input* icelevel_input = element->GetInput(MaskIceLevelsetEnum);                         _assert_(icelevel_input);
+
+                                for(int v=0;v<numvertices;v++) {
+                                        gauss->GaussVertex(v);
+                                        icelevel_input->GetInputValue(&ice,gauss);
+                                        /* The initial mask is very strict, we look at all grounded elements and set the mask for ice nodes only. */
+                                        if (ice < 0) local_mask[element->vertices[v]->Lid()] = 1.;
+                                }
+                                delete gauss;
                         }
                 }
@@ -48 +59 @@
                 /*Local iterations on partition*/
                 bool keepgoing    = true;
-                int  didsomething = 0;
                 int  iter         = 1;
                 while(keepgoing){
@@ -54 +64 @@
 
                         keepgoing    = false;
-                        didsomething = 0;
                         int i=0;
                         for(Object* & object : femmodel->elements->objects){
@@ -62 +71 @@
                                         int numvertices = element->GetNumberOfVertices();
                                         bool found1 = false;
+                                        IssmDouble sumlocalmask = 0.;
+
                                         for(int j=0;j<numvertices;j++){
                                                 lid = element->vertices[j]->Lid();
-                                                if(local_mask[lid]>0.){
+                                                /*we need to have at least a sharing edge, to extend the mask*/
+                                                sumlocalmask += local_mask[lid];
+                                                if(sumlocalmask > 1.5){
                                                         found1 = true;
                                                         break;
@@ -73 +86 @@
                                                 for(int j=0;j<numvertices;j++){
                                                         lid = element->vertices[j]->Lid();
-                                                        if(local_mask[lid]==0.){
-                                                                local_mask[lid]=1.;
-                                                                keepgoing = true;
-                                                                didsomething = 1;
-                                                        }
+                                                        local_mask[lid]=1.;
                                                 }
+                                                keepgoing = true;
                                         }
                                 }
@@ -99 +109 @@
         xDelete<bool>(element_flag);
 
+        int killbergReinit = 0;
         /*OK, now deactivate iceberg and count the number of deactivated vertices*/
         for(Object* & object : femmodel->elements->objects){
@@ -120 +131 @@
                                 element->AddInput(MaskIceLevelsetEnum,values,P1Enum);
                                 xDelete<IssmDouble>(values);
+                                killbergReinit += 1;
                         }
                 }
         }
-
         /*cleanup*/
         xDelete<IssmDouble>(local_mask);
+
+        /*Recompute the sign distance for the levelset function*/
+        return killbergReinit;
 }