source: inundation/ga/storm_surge/pyvolution/mesh.py @ 262

"""Classes implementing general 2D geometrical mesh of triangles.

   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
   Geoscience Australia, 2004   
"""


class Mesh:
    """Collection of triangular elements (purely geometric)

    A triangular element is defined in terms of three vertex ids,
    ordered counter clock-wise, 
    each corresponding to a given coordinate set.
    Vertices from different elements can point to the same
    coordinate set.

    Coordinate sets are implemented as an N x 2 Numeric array containing
    x and y coordinates. 
    

    To instantiate:
       Mesh(points, vertices)

    where

      points is either a list of 2-tuples or an Mx2 Numeric array of
      floats representing all x, y coordinates in the mesh.

      vertices is either a list of 3-tuples or an Nx3 Numeric array of
      integers representing indices of all vertices in the mesh.
      Each vertex is identified by its index i in [0, M-1].


    Example:
        a = [0.0, 0.0]
        b = [0.0, 2.0]
        c = [2.0,0.0]
        e = [2.0, 2.0]
        
        points = [a, b, c, e]
        vertices = [ [1,0,2], [1,2,3] ]   #bac, bce 
        mesh = Mesh(points, vertices)
    
        #creates two triangles: bac and bce


    Mesh takes the optional third argument boundary which is a
    dictionary mapping from (element_id, edge_id) to boundary tag.
    The default value is None which will assign the defualt_boundary_tag
    as specified in config.py to all boundary edges.
    """

    def __init__(self, coordinates, vertices, boundary = None):
        """
        Build triangles from x,y coordinates (sequence of 2-tuples or
        Mx2 Numeric array of floats) and vertices (sequence of 3-tuples
        or Nx3 Numeric array of non-negative integers).
        """

        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum

        #Only one of them
        self.vertices = array(vertices).astype(Int)
        self.coordinates = array(coordinates)


        #Input checks
        msg = 'Vertices must an Nx2 Numeric array or a sequence of 2-tuples'
        assert len(self.vertices.shape) == 2, msg

        msg = 'Coordinates must an Mx2 Numeric array or a sequence of 2-tuples'
        assert len(self.coordinates.shape) == 2, msg        

        msg = 'Vertex indices reference non-existing coordinate sets'
        assert max(max(self.vertices)) <= self.coordinates.shape[0], msg


        #Register number of elements (N)
        self.number_of_elements = N = self.vertices.shape[0]

        #Allocate space for geometric quantities
        self.centroids = zeros((N, 2), Float)
        #FIXME: Should be renamed to centroid_coordinates
        
        self.areas = zeros(N, Float)
        self.radii = zeros(N, Float)
        self.edgelengths = zeros((N, 3), Float)

        self.neighbours = zeros((N, 3), Int)
        self.neighbour_edges = zeros((N, 3), Int)
        self.number_of_boundaries = zeros(N, Int) 
        self.surrogate_neighbours = zeros((N, 3), Int)        
        
        self.normals = zeros((N, 6), Float)

        #Get x,y coordinates for all vertices for all triangles
        #and store
        self.vertex_coordinates = V = self.compute_vertex_coordinates()
        

        ##print 'Initialise mesh'
        #Initialise each triangle
        for i in range(N):
            #if i % (N/10) == 0: print '(%d/%d)' %(i, N)
            
            x0 = V[i, 0]; y0 = V[i, 1]
            x1 = V[i, 2]; y1 = V[i, 3]
            x2 = V[i, 4]; y2 = V[i, 5]            

            #Compute centroid
            centroid = array([(x0 + x1 + x2)/3, (y0 + y1 + y2)/3])
            self.centroids[i] = centroid

            #Area
            self.areas[i] = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2

            msg = 'Triangle (%f,%f), (%f,%f), (%f, %f)' %(x0,y0,x1,y1,x2,y2)
            msg += ' is degenerate:  area == %f' %self.areas[i]
            assert self.areas[i] > 0.0, msg
        
            #Midpoints        
            m0 = array([(x1 + x2)/2, (y1 + y2)/2])
            m1 = array([(x0 + x2)/2, (y0 + y2)/2])
            m2 = array([(x1 + x0)/2, (y1 + y0)/2])

            #The radius is the distance from the centroid of
            #a triangle to the midpoint of the side of the triangle
            #closest to the centroid
            d0 = sqrt(sum( (centroid-m0)**2 )) 
            d1 = sqrt(sum( (centroid-m1)**2 ))        
            d2 = sqrt(sum( (centroid-m2)**2 ))
        
            self.radii[i] = min(d0, d1, d2)

            #Normals
            #The normal vectors
            # - point outward from each edge
            # - are orthogonal to the edge
            # - have unit length
            # - Are enumerated according to the opposite corner:
            #   (First normal is associated with the edge opposite
            #    the first vertex, etc)
            # - Stored as six floats n0x,n0y,n1x,n1y,n2x,n2y per triangle

            n0 = array([x2 - x1, y2 - y1])
            l0 = sqrt(sum(n0**2))

            n1 = array([x0 - x2, y0 - y2])
            l1 = sqrt(sum(n1**2))

            n2 = array([x1 - x0, y1 - y0])
            l2 = sqrt(sum(n2**2))

            #Normalise
            n0 /= l0
            n1 /= l1
            n2 /= l2

            #Compute and store
            self.normals[i, :] = [n0[1], -n0[0],
                                  n1[1], -n1[0],
                                  n2[1], -n2[0]]             

            #Edgelengths
            self.edgelengths[i, :] = [l0, l1, l2]        

            #Initialise Neighbours (-1 means that it is a boundary neighbour)
            self.neighbours[i, :] = [-1, -1, -1]

            #Initialise edge ids of neighbours
            #In case of boundaries this slot is not used
            self.neighbour_edges[i, :] = [-1, -1, -1]


        #Build neighbour structure    
        self.build_neighbour_structure()

        #Build vertex list
        self.build_vertexlist()        

        #Build surrogate neighbour structure    
        self.build_surrogate_neighbour_structure()        

        #Build boundary dictionary mapping (id, edge) to symbolic tags        
        self.build_boundary_dictionary(boundary)

        #Update boundary indices 
        self.build_boundary_structure()        

        
    def __len__(self):
        return self.number_of_elements

    def __repr__(self):
        return 'Mesh: %d vertices, %d elements, %d boundary segments'\
               %(self.coordinates.shape[0], len(self), len(self.boundary))

            
    def build_neighbour_structure(self):
        """Update all registered triangles to point to their neighbours.
        
        Also, keep a tally of the number of boundaries for each triangle

        Postconditions:
          neighbours and neighbour_edges is populated
          number_of_boundaries integer array is defined.
        """
        
        #Step 1:
        #Build dictionary mapping from segments (2-tuple of points)
        #to left hand side edge (facing neighbouring triangle)
        #Also build list of vertices containing indices of triangles  

        N = self.number_of_elements        
        neighbourdict = {}
        vertexlist = [None]*len(self.coordinates)
        for i in range(N):

            #Register all segments as keys mapping to current triangle
            #and segment id
            a = self.vertices[i, 0]
            b = self.vertices[i, 1]
            c = self.vertices[i, 2]            
            neighbourdict[a,b] = (i, 2) #(id, edge)
            neighbourdict[b,c] = (i, 0) #(id, edge)
            neighbourdict[c,a] = (i, 1) #(id, edge)

            #Register the vertices as keys mapping to 
            #(triangle, edge) tuples associated with them
            #This is used for smoothing
            for vertex_id, v in enumerate([a,b,c]):
                if vertexlist[v] is None:
                    vertexlist[v] = []

                vertexlist[v].append( (i, vertex_id) )    

        #Step 2:
        #Go through triangles again, but this time
        #reverse direction of segments and lookup neighbours.
        for i in range(N):
            a = self.vertices[i, 0]
            b = self.vertices[i, 1]
            c = self.vertices[i, 2]

            self.number_of_boundaries[i] = 3
            if neighbourdict.has_key((b,a)):
                self.neighbours[i, 2] = neighbourdict[b,a][0]
                self.neighbour_edges[i, 2] = neighbourdict[b,a][1]
                self.number_of_boundaries[i] -= 1
                
            if neighbourdict.has_key((c,b)):
                self.neighbours[i, 0] = neighbourdict[c,b][0]
                self.neighbour_edges[i, 0] = neighbourdict[c,b][1]
                self.number_of_boundaries[i] -= 1               
                
            if neighbourdict.has_key((a,c)):
                self.neighbours[i, 1] = neighbourdict[a,c][0]
                self.neighbour_edges[i, 1] = neighbourdict[a,c][1]
                self.number_of_boundaries[i] -= 1              


        self.vertexlist = vertexlist
    
    def build_vertexlist(self):
        """Build vertexlist index by vertex ids and for each entry (point id)
        build a list of (triangles, vertex_id) pairs that use the point
        as vertex.

        Preconditions:
          self.coordinates and self.vertices are defined  
        
        Postcondition:
          self.vertexlist is build
        """

        vertexlist = [None]*len(self.coordinates)
        for i in range(self.number_of_elements):

            a = self.vertices[i, 0]
            b = self.vertices[i, 1]
            c = self.vertices[i, 2]            

            #Register the vertices as keys mapping to 
            #(triangle, edge) tuples associated with them
            #This is used for smoothing
            for vertex_id, v in enumerate([a,b,c]):
                if vertexlist[v] is None:
                    vertexlist[v] = []

                vertexlist[v].append( (i, vertex_id) )    

        self.vertexlist = vertexlist
    

    def build_surrogate_neighbour_structure(self):
        """Build structure where each triangle edge points to its neighbours
        if they exist. Otherwise point to the triangle itself.

        The surrogate neighbour structure is useful for computing gradients
        based on centroid values of neighbours.

        Precondition: Neighbour structure is defined
        Postcondition:
          Surrogate neighbour structure is defined:
          surrogate_neighbours: i0, i1, i2 where all i_k >= 0 point to
          triangles.

        """

        N = self.number_of_elements
        for i in range(N):
            #Find all neighbouring volumes that are not boundaries
            for k in range(3):
                if self.neighbours[i, k] < 0:
                    self.surrogate_neighbours[i, k] = i #Point this triangle
                else:
                    self.surrogate_neighbours[i, k] = self.neighbours[i, k]



    def build_boundary_dictionary(self, boundary = None):
        """Build or input the dictionary of boundary tags.
         self.boundary is a dictionary of tags,
         keyed by volume id and edge:
         { (id, edge): tag, ... }
         
         Postconditions:
            self.boundary is defined.
        """

        from config import default_boundary_tag
        
        if boundary is None:
            boundary = {}
            for vol_id in range(self.number_of_elements):
                for edge_id in range(0, 3):
                    if self.neighbours[vol_id, edge_id] < 0:
                        boundary[(vol_id, edge_id)] = default_boundary_tag
        else:
            #Check that all keys in given boundary exist
            for vol_id, edge_id in boundary.keys():
                msg = 'Segment (%d, %d) does not exist' %(vol_id, edge_id)
                a, b = self.neighbours.shape
                assert vol_id < a and edge_id < b, msg
                
                msg = 'Segment (%d, %d) is not a boundary' %(vol_id, edge_id)
                assert self.neighbours[vol_id, edge_id] < 0, msg

            #Check that all boundary segments are assigned a value
            for vol_id in range(self.number_of_elements):
                for edge_id in range(0, 3):
                    if self.neighbours[vol_id, edge_id] < 0:
                        if not boundary.has_key( (vol_id, edge_id) ):
                            msg = 'WARNING: Given boundary does not contain '
                            msg += 'tags for edge (%d, %d). '\
                                   %(vol_id, edge_id)
                            msg += 'Assigning default tag (%s).'\
                                   %default_boundary_tag

                            #FIXME: Print only as per verbosity
                            #print msg
                            boundary[ (vol_id, edge_id) ] =\
                                      default_boundary_tag
            

            
        self.boundary = boundary
            

    def build_boundary_structure(self):
        """Traverse boundary and 
        enumerate neighbour indices from -1 and
        counting down.

        Precondition:
            self.boundary is defined.
        Post condition:
            neighbour array has unique negative indices for boundary
            boundary_segments array imposes an ordering on segments
            (not otherwise available from the dictionary)             
        """

        if self.boundary is None:
            msg = 'Boundary dictionary must be defined before '
            msg += 'building boundary structure'
            raise msg

        
        self.boundary_segments = self.boundary.keys()
        self.boundary_segments.sort()

        index = -1        
        for id, edge in self.boundary_segments:
            self.neighbours[id, edge] = index
            index -= 1

    
    def get_boundary_tags(self):
        """Return list of available boundary tags
        """
        
        tags = {}
        for v in self.boundary.values():
            tags[v] = 1

        return tags.keys()    
            
    

    def check_integrity(self):
        """Check that triangles are internally consistent e.g.
        that area corresponds to edgelengths, that vertices
        are arranged in a counter-clockwise order, etc etc
        Neighbour structure will be checked by class Mesh
        """
        
        from config import epsilon
        from math import pi
        from util import anglediff

        N = self.number_of_elements
        
        #Get x,y coordinates for all vertices for all triangles
        V = self.get_vertex_coordinates()

        #Check each triangle
        for i in range(N):
            x0 = V[i, 0]; y0 = V[i, 1]
            x1 = V[i, 2]; y1 = V[i, 3]
            x2 = V[i, 4]; y2 = V[i, 5]            

            #Check that area hasn't been compromised
            area = self.areas[i]
            ref = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2
            msg = 'Wrong area for vertex coordinates: %f %f %f %f %f %f'\
                  %(x0,y0,x1,y1,x2,y2)
            assert abs((area - ref)/area) < epsilon, msg
        
            #Check that points are arranged in counter clock-wise order
            v0 = [x1-x0, y1-y0]
            v1 = [x2-x1, y2-y1]
            v2 = [x0-x2, y0-y2]

            a0 = anglediff(v1, v0)
            a1 = anglediff(v2, v1)
            a2 = anglediff(v0, v2)        

            msg = '''Vertices (%s,%s), (%s,%s), (%s,%s) are not arranged
            in counter clockwise order''' %(x0, y0, x1, y1, x2, y2)
            assert a0 < pi and a1 < pi and a2 < pi, msg

            #Check that normals are orthogonal to edge vectors
            #Note that normal[k] lies opposite vertex k

            normal0 = self.normals[i, 0:2]
            normal1 = self.normals[i, 2:4]
            normal2 = self.normals[i, 4:6]            

            for u, v in [ (v0, normal2), (v1, normal0), (v2, normal1) ]:
                assert u[0]*v[0] + u[1]*v[1] < epsilon


        #Check integrity of neighbour structure
        for i in range(N):
            for v in self.vertices[i, :]:
                #Check that all vertices have been registered
                assert self.vertexlist[v] is not None

                #Check that this triangle is listed with at least one vertex
                assert (i, 0) in self.vertexlist[v] or\
                       (i, 1) in self.vertexlist[v] or\
                       (i, 2) in self.vertexlist[v]
                
                

            #Check neighbour structure
            for k, neighbour_id in enumerate(self.neighbours[i,:]):

                #Assert that my neighbour's neighbour is me
                #Boundaries need not fulfill this
                if neighbour_id >= 0:
                    edge = self.neighbour_edges[i, k]
                    assert self.neighbours[neighbour_id, edge] == i



        #Check that all boundaries have
        # unique, consecutive, negative indices                    
        L = len(self.boundary)
        for i in range(L):
            id, edge = self.boundary_segments[i]
            assert self.neighbours[id, edge] == -i-1


        #NOTE: This assert doesn't hold true if there are internal boundaries
        #FIXME: Look into this further.
        #for id, edge in self.boundary:
        #    assert self.neighbours[id,edge] < 0
        

    

    def get_normals(self):
        """Return all normal vectors.
        Return normal vectors for all triangles as an Nx6 array
        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
        """
        return self.normals


    def get_normal(self, i, j):
        """Return normal vector j of the i'th triangle.
        Return value is the numeric array slice [x, y]
        """
        return self.normals[i, 2*j:2*j+2]     
    

            
    def get_vertex_coordinates(self):
        """Return all vertex coordinates.
        Return all vertex coordinates for all triangles as an Nx6 array
        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)        
        """
        return self.vertex_coordinates


    def get_vertex_coordinate(self, i, j):
        """Return coordinates for vertex j of the i'th triangle.
        Return value is the numeric array slice [x, y]
        """
        return self.vertex_coordinates[i, 2*j:2*j+2]     


    def compute_vertex_coordinates(self):
        """Return vertex coordinates for all triangles as an Nx6 array
        (ordered as x0, y0, x1, y1, x2, y2 for each triangle)
        """

        from Numeric import zeros, Float

        N = self.number_of_elements
        vertex_coordinates = zeros((N, 6), Float) 

        for i in range(N):
            for j in range(3):
                k = self.vertices[i,j]  #Index of vertex 0
                v_k = self.coordinates[k]
                vertex_coordinates[i, 2*j+0] = v_k[0]
                vertex_coordinates[i, 2*j+1] = v_k[1]

        return vertex_coordinates
  



def rectangular_mesh(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)):
    from mesh_factory import rectangular
    points, vertices, boundary = rectangular(m, n, len1, len2, origin)

    return Mesh(points, vertices, boundary)


#FIXME: Get oblique and contracting and circular meshes from Chris