Changeset 25597

Timestamp:

09/27/20 13:12:43 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: trying to simplify GetGrounded ratio

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (2 diffs)

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

@@ -1670 +1670 @@
         int               domaintype,index1,index2;
         const IssmPDouble epsilon        = 1.e-15;
-        IssmDouble        phi,s1,s2,area_init,area_grounded;
+        IssmDouble        phi,s1,s2;
         IssmDouble        gl[NUMVERTICES];
         IssmDouble        xyz_bis[3][3];
@@ -1708 +1708 @@
 
                         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_("Tria::GetGroundedPortion: levelset case not allowed for element #" << this->Sid()+1 << ":" << gl[0] << "|" << gl[1] << "|" << gl[2]);
-
-                        /*Compute fraction of grounded element*/
-                        if(domaintype==Domain3DsurfaceEnum){ //hack, need to be implemented in a Tria 3D
-                                GetJacobianDeterminant3D(&area_init, xyz_list,NULL);
-                                GetJacobianDeterminant3D(&area_grounded, &xyz_bis[0][0],NULL);
+                        if(mainlyfloating){
+                                phi = (1-s1*s2);
                         }
                         else{
-                                GetJacobianDeterminant(&area_init, xyz_list,NULL);
-                                GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
-                        }
-                        if(mainlyfloating==true) area_grounded=area_init-area_grounded;
-                        phi=area_grounded/area_init;
+                                phi = s1*s2;
+                        }
                 }
         }
         else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");
 
-        if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");
-
+        _assert_(phi<=1. && phi>=0.);
         return phi;
 }