Changeset 25598

Timestamp:

09/27/20 20:53:38 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: simplifying Penta GL subelement grounded ratio calculation

Location:

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

Files:

• Penta.cpp (modified) (2 diffs)
• Tria.cpp (modified) (2 diffs)

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

@@ -1320 +1320 @@
         bool               mainlyfloating = true;
         const IssmPDouble  epsilon= 1.e-15;
-        IssmDouble         phi,s1,s2,area_init,area_grounded;
+        IssmDouble         phi,s1,s2;
         IssmDouble         gl[NUMVERTICES];
-        IssmDouble         xyz_bis[NUMVERTICES2D][3];
 
         /*Recover parameters and values*/
@@ -1344 +1343 @@
 
                 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];
-
-                        /*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]);
-
-                        /*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]);
                 }
                 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*/
-                        xyz_bis[0][0]=*(xyz_list+3*0+0);
-                        xyz_bis[0][1]=*(xyz_list+3*0+1);
-                        xyz_bis[0][2]=*(xyz_list+3*0+2);
-
-                        /*Portion of the segments*/
                         s1=gl[0]/(gl[0]-gl[1]);
                         s2=gl[0]/(gl[0]-gl[2]);
-
-                        /*New point 1*/
-                        xyz_bis[1][0]=*(xyz_list+3*0+0)+s1*(*(xyz_list+3*1+0)-*(xyz_list+3*0+0));
-                        xyz_bis[1][1]=*(xyz_list+3*0+1)+s1*(*(xyz_list+3*1+1)-*(xyz_list+3*0+1));
-                        xyz_bis[1][2]=*(xyz_list+3*0+2)+s1*(*(xyz_list+3*1+2)-*(xyz_list+3*0+2));
-
-                        /*New point 2*/
-                        xyz_bis[2][0]=*(xyz_list+3*0+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*0+0));
-                        xyz_bis[2][1]=*(xyz_list+3*0+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*0+1));
-                        xyz_bis[2][2]=*(xyz_list+3*0+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*0+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
-                        /*Coordinates of point 1: same as initial point 2*/
-                        xyz_bis[1][0]=*(xyz_list+3*1+0);
-                        xyz_bis[1][1]=*(xyz_list+3*1+1);
-                        xyz_bis[1][2]=*(xyz_list+3*1+2);
-
-                        /*Portion of the segments*/
                         s1=gl[1]/(gl[1]-gl[0]);
                         s2=gl[1]/(gl[1]-gl[2]);
-
-                        /*New point 0*/
-                        xyz_bis[0][0]=*(xyz_list+3*1+0)+s1*(*(xyz_list+3*0+0)-*(xyz_list+3*1+0));
-                        xyz_bis[0][1]=*(xyz_list+3*1+1)+s1*(*(xyz_list+3*0+1)-*(xyz_list+3*1+1));
-                        xyz_bis[0][2]=*(xyz_list+3*1+2)+s1*(*(xyz_list+3*0+2)-*(xyz_list+3*1+2));
-
-                        /*New point 2*/
-                        xyz_bis[2][0]=*(xyz_list+3*1+0)+s2*(*(xyz_list+3*2+0)-*(xyz_list+3*1+0));
-                        xyz_bis[2][1]=*(xyz_list+3*1+1)+s2*(*(xyz_list+3*2+1)-*(xyz_list+3*1+1));
-                        xyz_bis[2][2]=*(xyz_list+3*1+2)+s2*(*(xyz_list+3*2+2)-*(xyz_list+3*1+2));
                 }
                 else _error_("case not possible");
-
-                /*Compute fraction of grounded element*/
-                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(mainlyfloating){
+                        phi = (1-s1*s2);
+                }
+                else{
+                        phi = s1*s2;
+                }
+        }
+
+        _assert_(phi<=1. && phi>=0.);
         return phi;
 }

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

@@ -1672 +1672 @@
         IssmDouble        phi,s1,s2;
         IssmDouble        gl[NUMVERTICES];
-        IssmDouble        xyz_bis[3][3];
 
         /*Recover parameters and values*/
@@ -1719 +1718 @@
                                 s2=gl[1]/(gl[1]-gl[2]);
                         }
-                        else _error_("Tria::GetGroundedPortion: levelset case not allowed for element #" << this->Sid()+1 << ":" << gl[0] << "|" << gl[1] << "|" << gl[2]);
+                        else _error_("case not possible");
                         if(mainlyfloating){
                                 phi = (1-s1*s2);