Changeset 22988

Timestamp:

07/20/18 14:27:12 (7 years ago)

Author:

Mathieu Morlighem

Message:

BUG: reverted change for FemModel and Tria

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/Tria.cpp (modified) (1 diff)
• FemModel.cpp (modified) (1 diff)

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

@@ -3301 +3301 @@
 }
 /*}}}*/
-void       Tria::ResetLevelsetFromSegmentlist(IssmDouble* distances){/*{{{*/
-
+void       Tria::ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments){/*{{{*/
+
+        /*Intermediaries*/
+        IssmDouble d,xn,yn;
 
         /*Get current levelset and vertex coordinates*/
         IssmDouble ls[NUMVERTICES];
+        IssmDouble  xyz_list[NUMVERTICES][3];
+        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
         GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
 
         /*Get distance from list of segments and reset ls*/
         for(int j=0;j<NUMVERTICES;j++){
+                IssmDouble dmin = 1.e+50;
+                for(int i=0;i<numsegments;i++){
+                        IssmDouble x = xyz_list[j][0];
+                        IssmDouble y = xyz_list[j][1];
+                        IssmDouble l2 = (segments[4*i+2]-segments[4*i+0])*(segments[4*i+2]-segments[4*i+0]) + (segments[4*i+3]-segments[4*i+1])*(segments[4*i+3]-segments[4*i+1]);
+
+                        /*Segment has a length of 0*/
+                        if(l2==0.){
+                                d = (x-segments[4*i+0])*(x-segments[4*i+0])+(y-segments[4*i+1])*(y-segments[4*i+1]);
+                                if(d<dmin) dmin = d;
+                                continue;
+                        }
+
+                        /*Consider the line extending the segment, parameterized as v + t (w - v).
+                         *We find projection of point p onto the line.
+                         *It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
+                        IssmDouble t = ((x-segments[4*i+0])*(segments[4*i+2]-segments[4*i+0]) + (y-segments[4*i+1])*(segments[4*i+3]-segments[4*i+1]))/l2;
+                        if(t < 0.0){
+                                // Beyond the 'v' end of the segment
+                                d = (x-segments[4*i+0])*(x-segments[4*i+0])+(y-segments[4*i+1])*(y-segments[4*i+1]);
+                        }
+                        else if (t > 1.0){
+                                // Beyond the 'w' end of the segment
+                                d = (x-segments[4*i+2])*(x-segments[4*i+2])+(y-segments[4*i+3])*(y-segments[4*i+3]);
+                        }
+                        else{
+                                // Projection falls on the segment
+                                xn = segments[4*i+0] + t * (segments[4*i+2] - segments[4*i+0]);
+                                yn = segments[4*i+1] + t * (segments[4*i+3] - segments[4*i+1]);
+                                d = (x-xn)*(x-xn)+(y-yn)*(y-yn);
+                        }
+
+                        if(d<dmin) dmin = d;
+                }
+
+                /*Update signed distance*/
+                dmin = sqrt(dmin);
                 if(dmin>10000) dmin=10000;
                 if(ls[j]>0){
-                        ls[j] = distances[this->vertices->Pid()];
+                        ls[j] = dmin;
                 }
                 else{
-                        ls[j] = - distances[this->vertices->Pid()];
+                        ls[j] = - dmin;
                 }
         }

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1050 +1050 @@
 
         /*3: Add distance input to all elements*/
-        IssmDouble* distances = xNew<IssmDouble>(vertices->Size());
-        IssmDouble  d,xn,yn;
-        for(int i=0;i<vertices->Size();i++){
-                Vertex* vertex=dynamic_cast<Vertex*>(this->vertices->GetObjectByOffset(i));
-                IssmDouble dmin = 1.e+50;
-
-                for(int i=0;i<numseg;i++){
-                        IssmDouble x = vertex->x;
-                        IssmDouble y = vertex->y;
-                        IssmDouble l2 = (allsegmentlist[4*i+2]-allsegmentlist[4*i+0])*(allsegmentlist[4*i+2]-allsegmentlist[4*i+0]) + (allsegmentlist[4*i+3]-allsegmentlist[4*i+1])*(allsegmentlist[4*i+3]-allsegmentlist[4*i+1]);
-
-                        /*Segment has a length of 0*/
-                        if(l2==0.){
-                                d = (x-allsegmentlist[4*i+0])*(x-allsegmentlist[4*i+0])+(y-allsegmentlist[4*i+1])*(y-allsegmentlist[4*i+1]);
-                                if(d<dmin) dmin = d;
-                                continue;
-                        }
-
-                        /*Consider the line extending the segment, parameterized as v + t (w - v).
-                         *We find projection of point p onto the line.
-                         *It falls where t = [(p-v) . (w-v)] / |w-v|^2*/
-                        IssmDouble t = ((x-allsegmentlist[4*i+0])*(allsegmentlist[4*i+2]-allsegmentlist[4*i+0]) + (y-allsegmentlist[4*i+1])*(allsegmentlist[4*i+3]-allsegmentlist[4*i+1]))/l2;
-                        if(t < 0.0){
-                                // Beyond the 'v' end of the segment
-                                d = (x-allsegmentlist[4*i+0])*(x-allsegmentlist[4*i+0])+(y-allsegmentlist[4*i+1])*(y-allsegmentlist[4*i+1]);
-                        }
-                        else if (t > 1.0){
-                                // Beyond the 'w' end of the segment
-                                d = (x-allsegmentlist[4*i+2])*(x-allsegmentlist[4*i+2])+(y-allsegmentlist[4*i+3])*(y-allsegmentlist[4*i+3]);
-                        }
-                        else{
-                                // Projection falls on the segment
-                                xn = allsegmentlist[4*i+0] + t * (allsegmentlist[4*i+2] - allsegmentlist[4*i+0]);
-                                yn = allsegmentlist[4*i+1] + t * (allsegmentlist[4*i+3] - allsegmentlist[4*i+1]);
-                                d = (x-xn)*(x-xn)+(y-yn)*(y-yn);
-                        }
-
-                        if(d<dmin) dmin = d;
-                }
-
-                /*Update signed distance*/
-                distances[vertex->pid] = sqrt(dmin);
-        }
-
-        //element->CreateDistanceInputFromSegmentlist(allsegmentlist,numseg,distanceenum);
-        InputUpdateFromVectorx(this,distances,distanceenum,VertexPIdEnum);
-        printf("-------------- file: FemModel.cpp line: %i\n",__LINE__); 
+        for(int i=0;i<elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+                if(!element->IsOnBase()) continue;
+                element->CreateDistanceInputFromSegmentlist(allsegmentlist,numseg,distanceenum);
+        }
+
+        /*Extrude if necessary*/
+        int elementtype;
+        this->parameters->FindParam(&elementtype,MeshElementtypeEnum);
+        if(elementtype==PentaEnum){
+                InputExtrudex(this,distanceenum,-1);
+        }
+        else if(elementtype==TriaEnum){
+                /*no need to extrude*/
+        }
+        else{
+                _error_("not implemented yet");
+        }
 
         /*Clean up and return*/
-        xDelete<IssmDouble>(distances);
         xDelete<IssmDouble>(allsegmentlist);
 }/*}}}*/