Changeset 3013

Timestamp:

May 30, 2006, 11:27:54 AM (18 years ago)

Author:

duncan

Message:

fixing the black screen of death bug at the mesh generation end.

Location:

inundation

Files:

• fit_interpolate/fit.py (modified) (1 diff)
• pmesh/mesh.py (modified) (3 diffs)
• pmesh/test_mesh.py (modified) (7 diffs)
• pmesh/test_mesh_interface.py (modified) (1 diff)
• triangle/test_triangmodule.py (modified) (1 diff)
• triangle/triang.c (modified) (8 diffs)

inundation/fit_interpolate/fit.py

@@ -318 +318 @@
 
         self._build_coefficient_matrix_B(verbose)
-        indexes = range(0,m)
         loners = self.mesh.get_lone_vertices()
-        #print "loners",loners
-        loners.sort()
-        loners.reverse()
-        #print "loners",loners
-        for lone_index in loners:
-            indexes.pop(lone_index)
-        #print "fit self.Atz",take(self.Atz,indexes)
+        # FIXME  - finish this.
+        #if len(loners)>0:
+        #    msg = 'ERROR (least_squares): Too few data points\n'
+        #    msg += 'There are only %d data points and alpha == 0. ' %n
+        #    msg += 'Need at least %d\n' %m
+        #    msg += 'Alternatively, set smoothing parameter alpha to a small '
+        #    msg += 'positive value,\ne.g. 1.0e-3.'
+        #    raise ToFewPointsError(msg)
+        
         
         return conjugate_gradient(self.B, self.Atz, self.Atz,

inundation/pmesh/mesh.py

@@ -1308 +1308 @@
         triangle attribute list: [(T1att), (T2att), ...]
             (list of a list of strings)
+        lone point list:[point1, ...] (list of integers)
         """
         #Clear the current generated mesh values
@@ -1375 +1376 @@
                                                    ObjectNeighbor[2])
             index += 1
-
+        genDict['lonepointlist'].sort() 
+        genDict['lonepointlist'].reverse()       
+        for loner in genDict['lonepointlist']:
+            # Remove the loner vertex
+            #print "Removing the loner", loner
+            self.meshVertices.pop(loner)
 
     def setMesh(self, genDict):
@@ -1925 +1931 @@
         load_mesh.loadASCII.export_mesh_file(ofile,dict)
 
+    # FIXME(DSG-DSG):Break this into two functions.
+    #One for the outline points.
+    #One for the mesh points.
+    # Note: this function is not in the gui
     def exportPointsFile(self,ofile):
         """

inundation/pmesh/test_mesh.py

@@ -231 +231 @@
 
 
-    def testisUserSegmentNew (self):
+    def testisUserSegmentNewII (self):
         mesh = Mesh()
         a = mesh.addUserVertex(0.0, 0.0)
@@ -943 +943 @@
         d = Vertex (1,2)
         e = Vertex (3,1)
-        f = Vertex (3,1)
+        #f = Vertex (3,1)
       
         s1 = Segment(a,b, tag = 50)
@@ -991 +991 @@
                         lFile[2] == "0.0,3.0" and
                         lFile[3] == "3.0,3.0" and
-                        lFile[4] == "1.0,2.0" and
-                        lFile[5] == "3.0,1.0"
+                        lFile[4] == "1.0,2.0"
                         ,
                         'exported Ascii xya file is wrong')
@@ -1906 +1905 @@
         
  
-    def not_ready_test_duplicat_verts_do_what(self):
+    def test_duplicat_verts_are_removed(self):
     
       
@@ -1913 +1912 @@
         c = Vertex (4.0,4.0)
         d = Vertex (4.0,0.0)
+        e = Vertex (4.0,0.0) # duplicate point
     
         s1 = Segment(a,b, tag = "50")
         s2 = Segment(b,c, tag = "40")
         s3 = Segment(c,d, tag = "30")
-        s4 = Segment(d,a, tag = "20")
-
-        e = Vertex (4.0,0.0) # duplicate point
+        s4 = Segment(d,e, tag = "no where seg")
+        s5 = Segment(e,a, tag = "20")
+
         
         m = Mesh(userVertices=[a,b,c,d,e],
-                 userSegments=[s1,s2,s3,s4])
+                 userSegments=[s1,s2,s3,s4,s5])
 
         seg = m.getUserSegments()
@@ -1928 +1928 @@
         holes = m.getHoles()
         regions = m.getRegions()
-        fileName = tempfile.mktemp(".tsh")
+        #fileName = tempfile.mktemp(".tsh")
         #fileName = "badmesh.tsh"
-        m.export_mesh_file(fileName)
+        #m.export_mesh_file(fileName)
         #print "***************************fileName", fileName
         #new_m = importMeshFromFile(fileName)
-        os.remove(fileName)
-        
-        m.generateMesh("Q", maxArea = 2.1 )
-
-        m.export_mesh_file("from_test_mesh.tsh")
+        #os.remove(fileName)
+        
+        m.generateMesh("Q", maxArea = 2000.1 )
+
+        #m.export_mesh_file("from_test_mesh.tsh")
         seg = m.getMeshSegments()
-        #print "seg",seg
-        #print "seg[0].tag"
-        #print seg[0].tag
-        #print "seg[0].tag"
-        
-        self.failUnless(seg[0].tag == 5 and
-                        seg[1].tag == 7 and
-                        seg[2].tag == 9 and
-                        seg[3].tag == 7 and
-                        seg[4].tag == 9,
-                        'seg tags are wrong')
-            
-        
-        #print '**@@@@@******'
-        #print "new_m",new_m 
-        #print '**@@@@@******'
-        #print "m",m 
-        #print '**@@@@@******'
-        
-
-        #self.failUnless( new_m == m,
-        #                 'loadASCIITestCase failed. test new 1')      
+        self.failUnless(4==len(seg),
+                        'FAILED!') 
+
+        vert = m.getMeshVertices() 
+        self.failUnless(4==len(vert),
+                        'FAILED!')
+ 
+    def test_duplicat_verts_are_removedII(self):
+    
+      
+        a = Vertex ( 0.0 ,0.0)
+        b = Vertex (0.0, 4.0)
+        c = Vertex (4.0,4.0)
+        d = Vertex (4.0,0.0)
+        e = Vertex (4.0,0.0) # duplicate point
+        f = Vertex (49.0,0.0) # unused point
+    
+        s1 = Segment(a,b, tag = "50")
+        s2 = Segment(b,c, tag = "40")
+        s3 = Segment(c,d, tag = "30")
+        s4 = Segment(d,e, tag = "no where seg")
+        s5 = Segment(e,a, tag = "20")
+
+        
+        m = Mesh(userVertices=[a,b,c,d,e,f],
+                 userSegments=[s1,s2,s3,s4,s5])
+
+        seg = m.getUserSegments()
+        points = m.getUserVertices()
+        holes = m.getHoles()
+        regions = m.getRegions()
+        #fileName = tempfile.mktemp(".tsh")
+        #fileName = "badmesh.tsh"
+        #m.export_mesh_file(fileName)
+        #print "***************************fileName", fileName
+        #new_m = importMeshFromFile(fileName)
+        #os.remove(fileName)
+        
+        m.generateMesh("Q", maxArea = 2000.1 )
+
+        #m.export_mesh_file("from_test_mesh.tsh")
+        seg = m.getMeshSegments()
+        self.failUnless(4==len(seg),
+                        'FAILED!') 
+
+        vert = m.getMeshVertices() 
+        self.failUnless(4==len(vert),
+                        'FAILED!')
+        
 def list_comp(A,B):
     yes = len(A)==len(B)
@@ -1973 +2000 @@
 if __name__ == "__main__":
     suite = unittest.makeSuite(meshTestCase,'test')
-    #suite = unittest.makeSuite(meshTestCase,'not_ready_test_duplicat_verts_do_what')
+    #suite = unittest.makeSuite(meshTestCase,'test_duplicat_verts_do_what')
     runner = unittest.TextTestRunner() #verbosity=2)
     runner.run(suite)

inundation/pmesh/test_mesh_interface.py

@@ -221 +221 @@
                         'FAILED!')
 
+
+    def test_create_mesh_from_regions_with_duplicate_verts(self):
+
+        # These are the absolute values
+        
+        polygon_absolute = [[0.0,0.0],
+                            [0,4.0],
+                            [4.0,4.0],
+                            [4.0,0.0],
+                            [4.0,0.0]]
+        
+        x_p = -10
+        y_p = -40
+        zone = 808
+        geo_ref_poly = Geo_reference(zone, x_p, y_p)
+        polygon = geo_ref_poly.change_points_geo_ref(polygon_absolute)
+
+        boundary_tags = {'50':[0],
+                         '40':[1],
+                         '30':[2],
+                         'no where seg':[3],
+                         '20':[4]
+                         }
+       
+        m = create_mesh_from_regions(polygon,
+                                     boundary_tags,
+                                     10000000,
+                                     poly_geo_reference=geo_ref_poly)
+        
+
+        fileName = "badmesh.tsh"
+        #m.export_mesh_file(fileName)
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":

inundation/triangle/test_triangmodule.py

@@ -298 +298 @@
                         'triangleattributelist is wrong!')
         
-        
+     
+    def test_lone_verts(self):
+
+        points = []
+        seglist = []
+        holelist = []
+        regionlist = []
+
+        points = [(0.0,0.0),(0.0,10.0),(3.0,0.0),(3.0,10.0),(0.0,10.0)]
+        pointattlist = [[],[],[],[],[]]
+        regionlist.append( (1.2,1.2,5.0) )
+        seglist = [(0,4),(4,1),(1,3),(3,2),(2,0)]
+        segattlist = [0,0,0,0,0]
+        trilist = []
+        mode = "Qzp"
+        data = triang.genMesh(points,seglist,holelist,regionlist,
+                              pointattlist,segattlist,trilist, mode)
+        #print "data['generatedtrianglelist']", data['generatedtrianglelist']
+        self.failUnless(data['generatedtrianglelist'] ==[(4, 0, 2), (2, 3, 4)],
+                        'trianglelist is wrong!')
+        self.failUnless(data['generatedsegmentlist'] ==[(0, 4), (4, 3),
+                                                        (3, 2), (2, 0)],
+                        'segmentlist is wrong!')
+        self.failUnless(data['generatedpointlist'] ==[(0.0, 0.0), (0.0, 10.0),
+                                                      (3.0, 0.0), (3.0, 10.0),
+                                                      (0.0,10.0)],
+                        ' is wrong!')
+        self.failUnless(data['lonepointlist'] ==[1],
+                        'lonepointlist is wrong!')
+        #print "", data['lonepointlist']
+   
+    def test_lone_vertsII(self):
+
+        points = []
+        seglist = []
+        holelist = []
+        regionlist = []
+
+        points = [(0.0,0.0),(0.0,0.0),(0.0,10.0),(0.0,10.0),(10.0,10.0),
+                  (10.0,10.0),(0.0,10.0),(10.0,0.0)]
+        
+        pointattlist = []
+        for point in points:
+            pointattlist.append([])
+        seglist = [(0,1),(1,2),(2,3),(3,4),(4,5),(5,7),(7,0)]       
+        segattlist = []
+        for seg in seglist:
+            segattlist.append(0)
+        trilist = []
+        mode = "Qzp"
+        data = triang.genMesh(points,seglist,holelist,regionlist,
+                              pointattlist,segattlist,trilist, mode)
+        #print "data['generatedtrianglelist']", data['generatedtrianglelist']
+        self.failUnless(data['generatedtrianglelist'] ==[(6, 1, 7), (7, 5, 6)],
+                        'trianglelist is wrong!')
+        self.failUnless(data['generatedsegmentlist'] ==[(1, 6), (6, 5),
+                                                        (5, 7), (7, 1)],
+                        'segmentlist is wrong!')
+        self.failUnless(data['generatedpointlist'] ==[(0.0,0.0),(0.0,0.0),
+                                                      (0.0,10.0),(0.0,10.0),
+                                                      (10.0,10.0),
+                                                      (10.0,10.0),(0.0,10.0),
+                                                      (10.0,0.0)],
+                        ' is wrong!')
+        self.failUnless(data['lonepointlist'] ==[0,2,3,4],
+                        'lonepointlist is wrong!')
+
 if __name__ == "__main__":
 

inundation/triangle/triang.c

@@ -76 +76 @@
   PyObject *Attrlist;
   REAL Attr;
-  int i, j, iatt, n;
+  int i, j, iatt, n, write_here;
   int a,b,c;
   int marker;
@@ -90 +90 @@
   int listsize,attsize ;
   PyObject  *ii;
-     
+  
+  int *points_connected;
+       
   if(!PyArg_ParseTuple(args,(char *)"OOOOOOOO",&hull,&seglist,&holelist,&regionlist,&pointattributelist,&segmarkerlist,&trianglelist,&mode)){
     return NULL;
@@ -369 +371 @@
   out.regionlist = (REAL *)NULL;
     
+  
   /*printf("\n\nTriangulate input args: %s \n\n", mod); */
   triangulate(mod, &in, &out, (struct triangulateio *)NULL );
@@ -379 +382 @@
   holder = PyDict_New();
      
+  /* list of int's, used to keep track of which verts are connected to
+     triangles. */
+  points_connected = (int *)malloc(out.numberofpoints*sizeof(int));
+  
   /* Add pointlist */
   listsize = out.numberofpoints;
@@ -387 +394 @@
                                     out.pointlist[i*2  ], out.pointlist[i*2+1]);  
     PyList_SetItem(holderlist,i, mlist);
+    points_connected[i] = 0;
   }  
   ii=PyString_FromString("generatedpointlist");
@@ -419 +427 @@
   ii=PyString_FromString("generatedpointattributelist");
   PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);  
+ 
   
   /* Add triangle list */
@@ -427 +436 @@
                                     out.trianglelist[i*3  ], out.trianglelist [i*3+1], out.trianglelist [i*3+2]);    
     PyList_SetItem(holderlist,i, mlist);
+    /* printf(" A vert index %i\n",out.trianglelist[i*3] );
+    printf(" A vert index %i\n",out.trianglelist[i*3+1] );
+    printf(" A vert index %i\n",out.trianglelist[i*3+2] ); */
+    points_connected[out.trianglelist[i*3]] = 1;
+    points_connected[out.trianglelist[i*3+1]] = 1;
+    points_connected[out.trianglelist[i*3+2]] = 1;
   }    
   ii=PyString_FromString("generatedtrianglelist");
   PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);
-     
+  
+  /* convert the points_connected vector from a true(1) false(0) vector, where
+     index is the vert, to a vector of the lone verts, at the beggining
+     of the vector. */
+  write_here = 0;   
+  for(i=0; i<out.numberofpoints;i++){
+    if (0 == points_connected[i]) {
+      points_connected[write_here] = i;
+      write_here ++;
+    }
+  }   
+  /* printf(" ******************** \n" );
+  for(i=0; i<write_here;i++){
+    printf(" A vert index %i\n",points_connected[i] );
+    } */
+  
+  listsize = write_here;
+  holderlist = PyList_New(listsize);
+  for(i=0; i<listsize;i++){
+   PyObject *mlist = Py_BuildValue((char *)"i", points_connected[i]);    
+    PyList_SetItem(holderlist,i, mlist);
+  }
+  ii=PyString_FromString("lonepointlist");
+  PyDict_SetItem(holder, ii, holderlist); Py_DECREF(ii); Py_DECREF(holderlist);  
      
   /* Add triangle attribute list */
@@ -531 +569 @@
     free(out.regionlist); out.regionlist=NULL;
   }
-    
+  if(!points_connected ){
+    free(points_connected ); points_connected =NULL;
+  }
+    
+  
+  
   /*  IN */