source: issm/oecreview/Archive/15392-16133/ISSM-15813-15814.diff@ 16134

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

@@ -817 +817 @@
         bool               floating=true;
         int                point;
         const IssmPDouble  epsilon= 1.e-15;
-        IssmDouble         gl[3];
+        IssmDouble         gl[NUMVERTICES];
         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;
+        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
@@ -868 +868 @@
 IssmDouble Tria::GetGroundedPortion(IssmDouble* xyz_list){
         /*Computeportion of the element that is grounded*/ 
 
-        bool               mainlyfloating = true;
-        const IssmPDouble  epsilon= 1.e-15;
+        bool              mainlyfloating = true;
+        const IssmPDouble epsilon        = 1.e-15;
         IssmDouble         phi,s1,s2,area_init,area_grounded;
-        IssmDouble         gl[3];
+        IssmDouble         gl[NUMVERTICES];
         IssmDouble         xyz_bis[3][3];
 
         /*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;
+        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

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

@@ -1023 +1023 @@
         bool               mainlyfloating = true;
         const IssmPDouble  epsilon= 1.e-15;
         IssmDouble         phi,s1,s2,area_init,area_grounded;
-        IssmDouble         gl[3];
-        IssmDouble         xyz_bis[3][3];
+        IssmDouble         gl[NUMVERTICES];
+        IssmDouble         xyz_bis[NUMVERTICES2D][3];
 
         /*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;
+        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
@@ -1047 +1047 @@
 
                 if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
                         /*Coordinates of point 2: same as initial point 2*/
-                        xyz_bis[2][0]=*(xyz_list+3*2+0);
-                        xyz_bis[2][1]=*(xyz_list+3*2+1);
-                        xyz_bis[2][2]=*(xyz_list+3*2+2);
+                        xyz_bis[2][0]=xyz_list[3*2+0];
+                        xyz_bis[2][1]=xyz_list[3*2+1];
+                        xyz_bis[2][2]=xyz_list[3*2+2];
 
                         /*Portion of the segments*/
                         s1=gl[2]/(gl[2]-gl[1]);
                         s2=gl[2]/(gl[2]-gl[0]);
 
                         /*New point 1*/
-                        xyz_bis[1][0]=*(xyz_list+3*2+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*2+0));
-                        xyz_bis[1][1]=*(xyz_list+3*2+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*2+1));
-                        xyz_bis[1][2]=*(xyz_list+3*2+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*2+2));
+                        xyz_bis[1][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
+                        xyz_bis[1][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
+                        xyz_bis[1][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
 
                         /*New point 0*/
-                        xyz_bis[0][0]=*(xyz_list+3*2+0)+s2*(*(xyz_list+3*0+0)-*(xyz_list+3*2+0));
-                        xyz_bis[0][1]=*(xyz_list+3*2+1)+s2*(*(xyz_list+3*0+1)-*(xyz_list+3*2+1));
-                        xyz_bis[0][2]=*(xyz_list+3*2+2)+s2*(*(xyz_list+3*0+2)-*(xyz_list+3*2+2));
+                        xyz_bis[0][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
+                        xyz_bis[0][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
+                        xyz_bis[0][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
                 }
                 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
                         /*Coordinates of point 0: same as initial point 2*/
@@ -1107 +1107 @@
                 }
 
                 /*Compute fraction of grounded element*/
-                GetJacobianDeterminant(&area_init, xyz_list,NULL);
-                GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
+                GetTriaJacobianDeterminant(&area_init, xyz_list,NULL);
+                GetTriaJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
                 if(mainlyfloating==true) area_grounded=area_init-area_grounded;
                 phi=area_grounded/area_init;
         }
 
-        if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
+        if(phi>1. || phi<0.) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
 
         return phi;
 }
@@ -7659 +7659 @@
         int         i,j;
         int         analysis_type,migration_style;
         IssmDouble  xyz_list[NUMVERTICES][3];
-        IssmDouble  xyz_list_tria[NUMVERTICES2D][3]={0.};
+        IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
         IssmDouble  alpha2,Jdet;
         IssmDouble  phi=1.0;
         IssmDouble  DL_scalar;
@@ -7692 +7692 @@
         friction=new Friction("2d",inputs,matpar,analysis_type);
 
         /*Recover portion of element that is grounded*/
-        if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list[0][0]);
+        if(migration_style==SubelementMigrationEnum) phi=this->GetGroundedPortion(&xyz_list_tria[0][0]);
 
         /* Start  looping on the number of gaussian points: */
         gauss=new GaussPenta(0,1,2,2);