source: issm/oecreview/Archive/22819-23185/ISSM-23172-23173.diff

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

@@ -55 +55 @@
                 xDelete<IssmDouble>(mask);
         }
 }/*}}}*/
-void GetMaskOfIceVerticesLSMx(FemModel* femmodel){/*{{{*/
 
+void GetMaskOfIceVerticesLSMx0(FemModel* femmodel){/*{{{*/
+
         /*Initialize vector with number of vertices*/
         int numvertices=femmodel->vertices->NumberOfVertices();
         if(numvertices==0)  return;
@@ -82 +83 @@
         delete vec_mask_ice;
         xDelete<IssmDouble>(mask_ice);
 }/*}}}*/
+void GetMaskOfIceVerticesLSMx(FemModel* femmodel){/*{{{*/
+
+        femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
+
+        /*Create vector on gset*/
+        int configuration_type;
+        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        int gsize=femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
+        if(gsize==0)  return;
+        Vector<IssmDouble>* vec_mask_ice=new Vector<IssmDouble>(gsize);
+
+        /*Fill vector with values: */
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+
+                if(element->IsIceInElement()){
+                        int numnodes = element->GetNumberOfNodes();
+                        int  gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum);
+                        int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum);
+                        IssmDouble* ones = xNew<IssmDouble>(gsize_local);
+                        for(int n=0;n<gsize_local;n++) ones[n] = 1.;
+                        vec_mask_ice->SetValues(gsize_local,glist_local,ones,INS_VAL);
+                        xDelete<IssmDouble>(ones);
+                        xDelete<int>(glist_local);
+                }
+        }
+
+        /*Assemble vector and serialize */
+        vec_mask_ice->Assemble();
+        IssmDouble* mask_ice=vec_mask_ice->ToMPISerial();
+        delete vec_mask_ice;
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                element->InputUpdateFromSolutionOneDof(mask_ice,IceMaskNodeActivationEnum);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(mask_ice);
+}/*}}}*/

../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h

@@ -8 +8 @@
 #include "../../classes/classes.h"
 
 void SetActiveNodesLSMx(FemModel* femmodel);
+void GetMaskOfIceVerticesLSMx0(FemModel* femmodel);
 void GetMaskOfIceVerticesLSMx(FemModel* femmodel);
-#endif  /* _UPDATESPCSX_H */
+#endif  /* _SETACTIVENODESLSMX_H*/

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

@@ -404 +404 @@
                         VerticesDofx(vertices,parameters);
 
                         if(VerboseMProcessor()) _printf0_("      detecting active vertices\n");
-                        GetMaskOfIceVerticesLSMx(this);
+                        GetMaskOfIceVerticesLSMx0(this);
                 }
 
                 if(VerboseMProcessor()) _printf0_("      resolving node constraints\n");
@@ -2632 +2632 @@
         delete this->constraints;       this->constraints       = new_constraints;
         delete this->materials;         this->materials = new_materials;
 
-        GetMaskOfIceVerticesLSMx(this);
+        GetMaskOfIceVerticesLSMx0(this);
 
         /*Insert MISMIP+ bed topography FIXME it could be stay in another place*/
         this->parameters->FindParam(&basalforcing_model,BasalforcingsEnum);

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

@@ -3958 +3958 @@
 
         _assert_(fieldvalue==0.); //field value != 0 not implemented yet
 
+        /*Get field on vertices (we do not allow for higher order elements!!)*/
+        IssmDouble lsf[NUMVERTICES];
+        this->GetInputListOnVertices(&lsf[0],fieldenum);
+
         /*1. check that we do cross fieldvalue in this element*/
-        Input* input = inputs->GetInput(fieldenum);
-        if(!input) _error_("Cannot calculate distance to "<<EnumToStringx(fieldenum)<<", input not found");
-        IssmDouble minvalue = input->Min();
+        IssmDouble minvalue = lsf[0];
+        IssmDouble maxvalue = lsf[0];
+        for(int i=1;i<NUMVERTICES;i++){
+                if(lsf[i]>maxvalue) maxvalue = lsf[i];
+                if(lsf[i]<minvalue) minvalue = lsf[i];
+        }
         if(minvalue>fieldvalue) return;
-        IssmDouble maxvalue = input->Max();
         if(maxvalue<fieldvalue) return;
 
         /*2. Find coordinates of where levelset crosses 0*/
@@ -3987 +3993 @@
                 }
         }
         else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
-                IssmDouble lsf[NUMVERTICES];
-                this->GetInputListOnVertices(&lsf[0],fieldenum);
+
                 for(int i=0;i<NUMVERTICES;i++){
                         if(lsf[indices[i]]==0.){
                                 x[counter]=xyz_list[indices[i]][0];

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

@@ -348 +348 @@
         /*Grab max of all cpus: */
         ISSM_MPI_Allreduce((void*)&flag,(void*)&allflag,1,ISSM_MPI_INT,ISSM_MPI_MAX,IssmComm::GetComm());
 
-        if(allflag){
-                return true;
-        }
-        else{
-                return false;
-        }
+        if(allflag) return true;
+        else        return false;
 }
 /*}}}*/