Changeset 18958

Timestamp:

12/12/14 01:01:37 (10 years ago)

Author:

jbondzio

Message:

minor: code cleanup (aka killing a code-monster)

File:

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

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

@@ -1275 +1275 @@
                         }
                         break;
-                case 2: // two vertices have ice: get area of quadrangle
+                case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
                         for(i=0;i<2;i++){
                                 fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
@@ -3197 +3197 @@
 /*}}}*/
 void       Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
-
-        int         normal_orientation=0;
-        IssmDouble  s1,s2;
-        IssmDouble  levelset[NUMVERTICES];
-
-        /*Recover parameters and values*/
-        IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
-        GetInputListOnVertices(&levelset[0],levelsetenum);
-
-        if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
-                /*Portion of the segments*/
-                s1=levelset[2]/(levelset[2]-levelset[1]);
-                s2=levelset[2]/(levelset[2]-levelset[0]);
-
-                if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle depending on distribution of levelsetfunction
-                /*New point 1*/
-                xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
-                xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
-                xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
-
-                /*New point 0*/
-                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
-                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
-                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
-        }
-        else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
-                /*Portion of the segments*/
-                s1=levelset[0]/(levelset[0]-levelset[2]);
-                s2=levelset[0]/(levelset[0]-levelset[1]);
-
-                if(levelset[0]<0.) normal_orientation=1;
-                /*New point 1*/
-                xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
-                xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
-                xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
-
-                /*New point 2*/
-                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
-                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
-                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
-        }
-        else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
-                /*Portion of the segments*/
-                s1=levelset[1]/(levelset[1]-levelset[0]);
-                s2=levelset[1]/(levelset[1]-levelset[2]);
-
-                if(levelset[1]<0.) normal_orientation=1;
-                /*New point 0*/
-                xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
-                xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
-                xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
-
-                /*New point 2*/
-                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
-                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
-                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
-        }
-        else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
-                xyz_zero[3*0+0]=xyz_list[0*3+0];
-                xyz_zero[3*0+1]=xyz_list[0*3+1];
-                xyz_zero[3*0+2]=xyz_list[0*3+2];
-
-                /*New point 2*/
-                xyz_zero[3*1+0]=xyz_list[1*3+0];
-                xyz_zero[3*1+1]=xyz_list[1*3+1];
-                xyz_zero[3*1+2]=xyz_list[1*3+2];
-        }
-        else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
-                xyz_zero[3*0+0]=xyz_list[2*3+0];
-                xyz_zero[3*0+1]=xyz_list[2*3+1];
-                xyz_zero[3*0+2]=xyz_list[2*3+2];
-
-                /*New point 2*/
-                xyz_zero[3*1+0]=xyz_list[0*3+0];
-                xyz_zero[3*1+1]=xyz_list[0*3+1];
-                xyz_zero[3*1+2]=xyz_list[0*3+2];
-        }
-        else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
-                xyz_zero[3*0+0]=xyz_list[1*3+0];
-                xyz_zero[3*0+1]=xyz_list[1*3+1];
-                xyz_zero[3*0+2]=xyz_list[1*3+2];
-
-                /*New point 2*/
-                xyz_zero[3*1+0]=xyz_list[2*3+0];
-                xyz_zero[3*1+1]=xyz_list[2*3+1];
-                xyz_zero[3*1+2]=xyz_list[2*3+2];
-        }
-        else _error_("Case not covered");
-
-        /*Assign output pointer*/
-        *pxyz_zero= xyz_zero;
+        /* Return coordinates where levelset intersects element edges.
+         * Attention: In case that no intersection exists, NULL pointer is returned.*/
+
+        /*Intermediaries*/
+        const int dim=3;
+        int numiceverts;
+        int i, n, e, counter;
+        IssmDouble s[2];
+        int* indices=NULL;
+        IssmDouble* xyz_zero=NULL;
+
+        this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
+        
+        //TODO: check if for 2 iceverts front segment is oriented in CCW way
+        
+        if(numiceverts>0) xyz_zero=xNew<IssmDouble>(2*dim);
+        if((numiceverts>0)&&(numiceverts<NUMVERTICES)){
+                counter=0;
+                for(i=0;i<numiceverts;i++){     // iterate over ice vertices
+                        for(n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
+                                for(e=0;e<dim;e++){ // spatial direction
+                                        int ind_ice             =dim*indices[i]+e;
+                                        int ind_noice   =dim*indices[n]+e;
+                                        int ind                 =dim*counter+e;
+                                        xyz_zero[ind]=xyz_list[ind_ice]+s[counter]*(xyz_list[ind_noice]-xyz_list[ind_ice]);
+                                }
+                                counter++;
+                        }
+                }
+        }
+        else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
+                IssmDouble lsf[NUMVERTICES];
+                this->GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+                counter=0;
+                for(i=0;i<NUMVERTICES;i++){
+                        if(lsf[indices[i]]==0.){
+                                for(e=0;e<dim;e++)      xyz_zero[dim*counter+e]=xyz_list[dim*indices[i]+e];
+                                counter++;
+                        }
+                        if(counter==2) break;
+                }
+        }
+        _assert_(counter==2);
+
+        /*Cleanup & return*/
+        xDelete<int>(indices);
+        *pxyz_zero=xyz_zero;
 }
 /*}}}*/