Changeset 1217 for inundation/ga/storm_surge

Timestamp:

Apr 11, 2005, 6:23:47 PM (20 years ago)

Author:

prow

Message:

Nearly working polygon from triangles.

Location:

inundation/ga/storm_surge/pmesh

Files:

• mesh.py (modified) (11 diffs)
• pmesh.py (modified) (2 diffs)

inundation/ga/storm_surge/pmesh/mesh.py

@@ -2117 +2117 @@
                           colour = colour)
 
-
     def weed(self,Vertices,Segments):
         #weed out existing duplicates
@@ -2129 +2128 @@
             point = (vertex.x,vertex.y)
             point_keys[point]=vertex
+        #inlined would looks very ugly
 
         line_keys = {}
@@ -2151 +2151 @@
         return Vertices,Segments
 
+    def segs_to_dict(self,segments):
+        dict={}
+        for segment in segments:
+            vertex1 = segment.vertices[0]
+            vertex2 = segment.vertices[1]
+            point1 = (vertex1.x,vertex1.y)
+            point2 = (vertex2.x,vertex2.y)
+            line = (point1,point2)
+            dict[line]=segment
+        return dict
+
+    def seg2line(self,s):
+        return ((s.vertices[0].x,s.vertices[0].y,)\
+                (s.vertices[1].x,s.vertices[1].y))
+
+    def line2seg(self,line,tag=None):
+        point0 = self.ver2point(line[0])
+        point1 = self.ver2point(line[1])
+        return Segment(point0,point1,tag=tag)
+
+    def ver2point(self,vertex):
+        return (vertex.x,vertex.y)
+
+    def point2ver(self,point):
+        return Vertex(point[0],point[1])
+
     def triangles_to_polySet(self,setName):
+        seg2line = self.seg2line
+        ver2point = self.ver2point
+        line2seg = self.line2seg
+        point2ver = self.point2ver
+
+        from Numeric import array,allclose
         #turn the triangles into a set
         Triangles = self.sets[self.setID[setName]]
@@ -2158 +2190 @@
             Triangles_dict[triangle]=None
 
+        point_keys = {}
+        userVertices={}
+        for vertex in self.getUserVertices():
+            point = ver2point(vertex)
+            if not point_keys.has_key(point):
+                point_keys[point]=vertex
+                userVertices[vertex]=True
+
+        userSegments = self.segs_to_dict(self.userSegments)
+        affine_lines = Affine_Linespace()
+        for line in userSegments.keys():
+            affine_lines.append(line)
+
+        alphaSegments = self.segs_to_dict(self.alphaUserSegments)
+        for line in alphaSegments.keys():
+            affine_lines.append(line)
+        
+        for triangle in Triangles:
+            for i in (0,1,2):
+                #for every triangles neighbour:
+                if not Triangles_dict.has_key(triangle.neighbors[i]):
+                #if the neighbour is not in the set:
+                    a = triangle.vertices[i-1]
+                    b = triangle.vertices[i-2]
+                    for vert in (a,b):
+                        point = ver2point(vert)
+                        if not point_keys.has_key(point):
+                            userVertices[vert]=True
+                            point_keys[point]=vert
+
+                    #Get possible matches:
+                    point_a = ver2point(a)
+                    point_b = ver2point(b)
+                    midpoint = ((a.x+b.x)/2,(a.y+b.y)/2)
+                    line = (point_a,point_b)
+                    tag = None
+                    possible_lines = affine_lines[line] 
+                    found = 0                            
+                    for user_line in possible_lines:
+                        if self.point_on_line(midpoint,user_line):
+                            if user_line in userSegments:
+                                parent_segment = userSegments.pop(user_line)
+                            else:
+                                parent_segment = alphaSegments.pop(user_line)
+                            tag = parent_segment.tag
+                            offspring = [line]
+                            offspring.extend(self.subtract_line(user_line,line))
+                            affine_lines.remove(user_line)
+                            for newline in offspring:
+                                line_vertices = []
+                                for point in newline:
+                                    if point_keys.has_key(point):
+                                        vert = point_keys[point]
+                                    else:
+                                        vert = Vertex(point[0],point[1])
+                                        userVertices[vert]=True
+                                        point_keys[point]=vert
+                                    line_vertices.append(vert)
+                                segment = Segment(line_vertices[0],line_vertices[1],tag)
+                                userSegments[newline]=segment
+                                affine_lines.append(newline)
+                            found=1
+                            break
+                    if not found:
+                        line_vertices = []
+                        for point in line:
+                            if point_keys.has_key(point):
+                                vert = point_keys[point]
+                            else:
+                                vert = Vertex(point[0],point[1])
+                                userVertices[vert]=True
+                                point_keys[point]=vert
+                            line_vertices.append(vert)
+                        segment = Segment(line_vertices[0],line_vertices[1],tag)
+                        userSegments[line]=segment
+                        affine_lines.append(line)
+        
+        self.userVerticies = []
+        self.userSegments = []
+        self.alphaSegments = []
+        #for key in userSegments.keys():
+        #    self.userSegments.append(userSegments[key])
+
+        #for key in alphaSegments.keys():
+        #    self.userSegments.append(alphaSegments[key])
+        #FIXME change alpha to user and sync with pmesh
+        #ie self.alphaSegments.append(alphaSegments[key])
+
+        return userVertices,userSegments,alphaSegments
+
+    def subtract_line(self,parent,child):
+        #Subtracts child from parent
+        #Requires that the child is a 
+        #subline of parent to work.
+
+        from Numeric import allclose,dot,array
+        A= parent[0]
+        B= parent[1]
+        a = child[0]
+        b = child[1]
+
+        A_array = array(parent[0])
+        B_array = array(parent[1])
+        a_array   = array(child[0])
+        b_array   = array(child[1])
+
+        assert not A == B
+        assert not a == b
+
+        answer = []
+
+        if not (allclose(A_array,a_array)\
+             or allclose(B_array,b_array)\
+             or allclose(A_array,b_array)\
+             or allclose(a_array,B_array)):
+            if dot(A_array-a_array,A_array-a_array) \
+            < dot(A_array-b_array,A_array-b_array):
+                sibling1 = (A,a)
+                sibling2 = (B,b)
+                return [sibling1,sibling2]
+            else:
+                sibling1 = (A,b)
+                sibling2 = (B,a)
+                return [sibling1,sibling2]
+
+        elif allclose(A_array,a_array):
+            if allclose(B_array,b_array):
+                return []
+            else:
+                sibling = (b,B)
+                return [sibling]
+        elif allclose(B_array,b_array):
+            sibling = (a,A)
+            return [sibling]
+
+        elif allclose(A_array,b_array):
+            if allclose(B,a):
+                return []
+            else:
+                sibling = (a,B)
+                return [sibling]
+        elif allclose(a_array,B_array):
+            sibling = (b,A)
+            return [sibling]
+
+
+    def point_on_line(self,point,line):
+        x=point[0]
+        y=point[1]
+        x0=line[0][0]
+        x1=line[1][0]
+        y0=line[0][1]
+        y1=line[1][1]
+        return point_on_line(x,y,x0,y0,x1,y1)
+        #a wrapper for the c func.
+
+    def line_length(self,line):
+        x0=line[0][0]
+        x1=line[1][0]
+        y0=line[0][1]
+        y1=line[1][1]
+        return ((x1-x0)**2-(y1-y0)**2)**0.5     
+
+    def triangles_to_polySet2(self,setName):
+        from Numeric import array,allclose
+
+        #turn the triangles into a set
+        Triangles = self.sets[self.setID[setName]]
+        Triangles_dict = {}
+        for triangle in Triangles:
+            Triangles_dict[triangle]=None
+
         userVertices = {}
-        userSegments = []
+        userSegments = {}
         point_keys = {}
+        affine_lines = {}
         for vertex in self.getUserVertices():
             point = (vertex.x,vertex.y)
@@ -2167 +2372 @@
         line_keys = {}
         for segment in self.getUserSegments():
+        #inlined would looks very ugly
             vertex1 = segment.vertices[0]
             vertex2 = segment.vertices[1]
@@ -2179 +2385 @@
             line1 = (point1,point2)
             line2 = (point2,point1)
+            AB = affine_line(point1,point2)
+            flat_AB_up = ceilAB[0]+AB[1]+AB[2]
             if not (line_keys.has_key(line1) \
                  or line_keys.has_key(line2)):
@@ -2209 +2417 @@
                         b=point_keys[(b.x,b.y)]
                         userVertices[b]=point_keys[(b.x,b.y)]
-                    assert userVertices.has_key(b)  
+                    assert userVertices.has_key(b)
 
                     if not (line_keys.has_key(((a.x,a.y),(b.x,b.y)))\
@@ -2217 +2425 @@
                         assert ((a.x,a.y)!=(b.x,b.y))
                         assert a!=b
-                        userSegments.append(Segment(a,b))
+                        AB = affine_line(a,b)
+                        flat_AB = AB[0]+AB[1]+AB[2]
+                        if affine_lines.has_key(flat_AB):
+                            userSegments[Segment(a,b)]=None
                         line_keys[((a.x,a.y),(b.x,b.y))]=None
 
-        userVertices,userSegments = self.weed(userVertices.keys(),userSegments)
+        userVertices,userSegments = self.weed(userVertices.keys(),userSegments.keys())
         self.userVertices.extend(userVertices)
         self.userSegments.extend(userSegments)
@@ -2296 +2507 @@
         self.meshTriangles[i]=replacement
         assert replacement in self.meshTriangles
-
-
-
 
 def importUngenerateFile(ofile):
@@ -2834 +3042 @@
 
     def repairCaseOne(self):
-        r = self.r
+        r = self.rkey
 
 
@@ -2897 +3105 @@
         triangle1.neighbors[j]=triangle2
 
-         
-
+class Discretised_Tuple_Set:
+    """
+    if a={(0.0):[(0.01),(0.02)],(0.2):[(0.17)]}
+    a[(0.01)]=a[(0.0)]=[(0.01),(0.02)]
+    a[(10000)]=[] #NOT KEYERROR
+
+    a.append[(0.01)]
+    => {0.0:[(0.01),(0.02),(0.01)],0.2:[(0.17)]}
+
+    #NOT IMPLIMENTED
+    a.remove[(0.01)]
+    => {(0.0):[(0.02),(0.01)],0.2:[(0.17)]}
+
+    a.remove[(0.17)]
+    => {(0.0):[(0.02),(0.01)],0.2:[]}
+    #NOT IMPLIMENTED
+    at a.precision = 2:
+        a.round_up[0.0]=
+        a.round_flat[0.0]=
+        a.round_downp0.0]=
+
+        a.up((0.1,2.04))=
+
+    If tolerance = 0, nothing gets rounded into
+    two bins. If tolerance = 0.5, everything does.
+
+    Ideally, precision can be set high, so that multiple
+    entries are rarely in the same bin. And tolerance should
+    be low (<0.1 for 1 dimension!,<(0.1/n) for small n!!)
+    so that it is rare to put items in mutiple bins.
+
+    Ex bins per entry = product(a,b,c...,n)
+    a = 1 or 2 s.t. Ex(a) = 1+2*tolerance
+    b = 1 or 2 ... 
+
+    BUT!!! to avoid missing the right bin:
+    (-10)**(precision+1)*tolerance must be greater than the 
+    greatest possible variation that an identical element
+    can display.
+    """
+    def __init__(self,precision = 6,tolerance = 0.01):
+        self.__precision__ = precision
+        self.__tolerance__ = tolerance
+        assert tolerance <= 0.5
+        self.__items__ = {}
+        from math import frexp
+        self.frexp = frexp
+
+    def __repr__(self):
+        return '%s'%self.__items__
+
+    def rounded_keys(self,key):
+        key = tuple(key)
+        keys = [key]
+        keys = self.__rounded_keys__(key)
+        return (keys)
+
+    def __rounded_keys__(self,key):
+        keys = []
+        up_key=list(key)
+        for i in range(len(up_key)):
+            up_key[i]=self.round_up(key[i])
+        keys.append(tuple(up_key))
+
+        #round down the value
+        for i in range(len(up_key)):
+            if not keys[0][i]==self.round_down(key[i]):
+            #unless round_up = round_down 
+            #so the keys stay unique
+                for j in range(len(keys)):
+                    down_key = list(keys[j])
+                    down_key[i]=self.round_down(key[i])
+                    down_key=tuple(down_key)
+                    keys.append(down_key)
+        return keys
+
+    def append(self,item):
+        keys = self.rounded_keys(item)
+        for key in keys:
+            if self.__items__.has_key(key):
+                self.__items__[key].append(item)
+            else:
+                self.__items__[key]=[item]
+
+    def __getitem__(self,key):
+        answer = []
+        keys = self.rounded_keys(key)
+        for key in keys:
+            if self.__items__.has_key(key):
+                answer.extend(self.__items__[key])
+        #if len(answer)==0:
+        #    raise KeyError#FIXME or return KeyError
+        #                  #FIXME or just return []?
+        else:
+            return answer #FIXME or unique(answer)?
+
+    def __delete__(self,item):
+        keys = self.rounded_keys(item)
+        answer = False
+        #if any of the possible keys contains
+        #a list, return true
+        for key in keys:        
+            if self.__items__.has_key(key):
+                if item in self.__items__[key]:
+                    self.__items__[key].remove(item)
+
+    def remove(self,item):
+        self.__delete__(item)
+
+    def __contains__(self,item):
+
+        keys = self.rounded_keys(item)
+        answer = False
+        #if any of the possible keys contains
+        #a list, return true
+        for key in keys:        
+            if self.__items__.has_key(key):
+                if item in self.__items__[key]:
+                    answer = True
+        return answer
+
+
+    def has_item(self,item):
+        return self.__contains__(item)
+
+    def round_up(self,value):
+         tolerance=self.__tolerance__
+         #Rounding up the value
+         m,e = self.frexp(value)
+         #m is the mantissa, e the exponent
+         # 0.5 < |m| < 1.0
+         m = m*(1+tolerance*(10**-(self.__precision__+1))/((m**2)**0.5))
+         #bump m up, by 5% or whatever
+         #self.precision dictates
+         m = round(m,self.__precision__)
+     #round m off so that the bump up 
+     #will be discernable if and only if
+         return m*(2.**e)
+
+    def round_down(self,value):
+         tolerance=self.__tolerance__
+         #Rounding down the value
+         m,e = self.frexp(value)
+         #m is the mantissa, e the exponent
+         # 0.5 < m < 1.0
+         m = m*(1-tolerance*(10**-(self.__precision__+1))/((m**2)**0.5))
+         #bump the |m| down, by 5% or whatever
+         #self.precision dictates
+         m = round(m,self.__precision__)
+     #round m off so that the bump down 
+     #will be discernable if and only if
+         return m*(2.**e)
+
+    def round_flat(self,value):
+    #redundant
+         m,e = self.frexp(value)
+         m = round(m,self.__precision__)
+         return m*(2.**e)
+
+    def keys(self):
+        return self.__items__.keys()
+
+    
+
+
+class Mapped_Discretised_Tuple_Set(Discretised_Tuple_Set):
+    """
+    This is a discretised tuple set, but 
+    based on a mapping. The mapping MUST
+    return a sequence.
+
+    example: 
+    def weight(animal):
+        return [animal.weight]
+    
+    a = Mapped_Discretised_Tuple_Set(weight)
+    a.append[cow]
+    a.append[fox]
+    a.append[horse]
+
+    a[horse]    -> [cow,horse]
+    a[dog]      -> [fox]
+    a[elephant] -> []
+    """
+    def __init__(self,mapping,precision = 6, tolerance=0.01):
+        Discretised_Tuple_Set.__init__\
+         (self,precision,tolerance = tolerance)
+        self.mapping = mapping
+
+    def rounded_keys(self,key):
+        mapped_key = tuple(self.mapping(key))
+        keys = self.__rounded_keys__(mapped_key)
+        return keys
+
+class Affine_Linespace(Mapped_Discretised_Tuple_Set):
+    """
+    The affine linespace creates a way to record and compare lines.
+    Precision is a bit of a hack, but it creates a way to avoid 
+    misses caused by round offs (between two lines of different
+    lenghts, the short one gets rounded off more). 
+    I am starting to think that a quadratic search would be faster.
+    Nearly.
+    """
+    def __init__(self,precision = 4,tolerance=0.04):
+        Mapped_Discretised_Tuple_Set.__init__\
+            (self,self.affine_line,\
+            precision=precision,tolerance=tolerance)
+
+    def affine_line(self,line):
+        point_1 = line[0]
+        point_2 = line[1]
+        #returns the equation of a line
+        #between two points, in the from
+        #(a,b,-c), as in ax+by-c=0
+        #or line *dot* (x,y,1) = (0,0,0)
+
+        #Note that it normalises the line
+        #(a,b,-c) so Line*Line = 1.
+        #This does not change the mathematical
+        #properties, but it makes comparism
+        #easier.
+
+        #There are probably better algorithms.
+        x1 = point_1[0]
+        x2 = point_2[0]
+        y1 = point_1[1]
+        y2 = point_2[1]
+        dif_x = x1-x2
+        dif_y = y1-y2
+
+        if dif_x == dif_y == 0:
+            msg = 'points are the same'
+            raise msg
+        elif abs(dif_x)>=abs(dif_y):
+            alpha = (-dif_y)/dif_x
+            #a = alpha * b
+            b = 1.
+            c = (x1*alpha+x2*alpha+y1+y2)/2.
+            a = alpha*b
+        else:
+            beta = dif_x/(-dif_y)
+            #b = beta * a
+            a = 1.
+            c = (x1+x2+y1*beta+y2*beta)/2.
+            b = beta*a
+
+        mag = (a**2+b**2+c**2)**(0.5)
+
+        if a == 0:
+            sign_a = 1.
+        else:
+            sign_a = a/((a**2)**0.5)
+        if b == 0:
+            sign_b = 1.
+        else:
+            sign_b = b/((b**2)**0.5)
+        if c == 0:
+            sign_c = 1.
+        else:
+            sign_c = c/((c**2)**0.5)
+        a = a/mag*sign_a
+        b = b/mag*sign_b
+        c = c/mag*sign_c
+        return a,b,c
+
+    
 if __name__ == "__main__":
     #from mesh import *

inundation/ga/storm_surge/pmesh/pmesh.py

@@ -339 +339 @@
 
     def triangles_to_polySet(self,parent):
-        self.mesh.triangles_to_polySet(self.selSet)
-        #event = None
-        #for v in vertices.keys():
-        #    if vertices[v] is True:
-        #        x = v.x*self.SCALE
-        #        y = v.y*self.SCALE
-        #        vertices[v]=self.drawVertex(x,y,event)
-        #print 'segments'
-        #print segments
-        #print 'vertices.keys()'
-        #print vertices.keys()
-        #for s in segments:
-        #    v0 = vertices[s[0]]
-        #    v1 = vertices[s[1]]
-        #    self.drawSegment(v0,v1)
+        userVertices,userSegments,alphaSegments = \
+            self.mesh.triangles_to_polySet(self.selSet)
+
         self.canvas.delete(ALL)
+
+        self.Vertices = visualmesh.vPoints(mesh.Vertex)
+        self.Segments = visualmesh.vSegments(mesh.Segment)
+
+        event = None
+        print 'len(userVertices.keys())'
+        print len(userVertices.keys())
+        print 'len(userSegments.keys())'
+        print len(userSegments.keys())
+        print 'len(alphaSegments.keys())'
+        print len(alphaSegments.keys())
+        for v in userVertices.keys():
+            if userVertices[v] is True:
+                x = v.x*self.SCALE
+                y = v.y*self.SCALE
+                userVertices[(v.x,v.y)]=self.drawVertex(x,y,event)
+
+        for line in userSegments.keys():
+            v0 = userVertices[line[0]]
+            v1 = userVertices[line[1]]
+            segment = self.drawSegment(v0,v1)
+            segment.set_tag(userSegments[line].tag)
+
+        for line in alphaSegments.keys():
+            v0 = userVertices[line[0]]
+            v1 = userVertices[line[1]]
+            segment = self.drawSegment(v0,v1)
+            segment.set_tag(alphaSegments[line].tag)
         self.visualiseMesh(self.mesh)
 
@@ -458 +474 @@
         region = self.drawRegion(x_origin, y_origin, 0)
         region.setTag("setheight5")
+
+        
         
         x_origin = 30