Changeset 15054

Timestamp:

05/21/13 09:28:30 (12 years ago)

Author:

seroussi

Message:

BUG: fixed grounding line problem

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• Tria.cpp (modified) (3 diffs)
• Tria.h (modified) (1 diff)

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

@@ -1240 +1240 @@
 
         return phi;
+}
+/*}}}*/
+/*FUNCTION Tria::GetGroundedPart{{{*/
+void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){
+        /*Computeportion of the element that is grounded*/ 
+
+        bool               floating=true;
+        int                point;
+        const IssmPDouble  epsilon= 1.e-15;
+        IssmDouble         gl[3];
+        IssmDouble         f1,f2;
+
+
+        /*Recover parameters and values*/
+        GetInputListOnVertices(&gl[0],GLlevelsetEnum);
+
+        /*Be sure that values are not zero*/
+        if(gl[0]==0) gl[0]=gl[0]+epsilon;
+        if(gl[1]==0) gl[1]=gl[1]+epsilon;
+        if(gl[2]==0) gl[2]=gl[2]+epsilon;
+
+        /*Check that not all nodes are grounded or floating*/
+        if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
+                point=0;
+                f1=1.;
+                f2=1.;
+        }
+        else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
+                point=0;
+                f1=0.;
+                f2=0.;
+        }
+        else{
+                if(gl[0]*gl[1]*gl[2]<0) floating=false;
+
+                if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+                        point=2;
+                        f1=gl[2]/(gl[2]-gl[0]);
+                        f2=gl[2]/(gl[2]-gl[1]);
+                }
+                else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+                        point=0;
+                        f1=gl[0]/(gl[0]-gl[1]);
+                        f2=gl[0]/(gl[0]-gl[2]);
+                }
+                else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+                        point=1;
+                        f1=gl[1]/(gl[1]-gl[2]);
+                        f2=gl[1]/(gl[1]-gl[0]);
+                }
+        }
+        *point1=point;
+        *fraction1=f1;
+        *fraction2=f2;
+        *mainlyfloating=floating;
 }
 /*}}}*/
@@ -3316 +3371 @@
 
         /*Intermediaries*/
+        bool       mainlyfloating;
         int        i,j;
         int        analysis_type,migration_style;
-        int        gausspoints=2;
+        int        point1;
         IssmDouble alpha2;
         IssmDouble Jdet;
+        IssmDouble fraction1,fraction2;
         IssmDouble phi=1.0;
         IssmDouble L[2][numdof];
@@ -3352 +3409 @@
         if(migration_style==SubelementMigration2Enum){
                 gllevelset_input=inputs->GetInput(GLlevelsetEnum); _assert_(gllevelset_input);
-                gausspoints=20;
+                this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+                gauss=new GaussTria(point1,fraction1,fraction2,mainlyfloating,2);
+        }
+        else{
+                gauss=new GaussTria(2);
         }
 
         /* Start  looping on the number of gaussian points: */
-        gauss=new GaussTria(gausspoints);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
 

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

@@ -204 +204 @@
                 void           GetConnectivityList(int* connectivity);
                 IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
+                void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
                 void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
                 void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);