"""Classes implementing general 2D geometrical mesh of triangles. Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou Geoscience Australia, 2004 """ from general_mesh import General_mesh class Mesh(General_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(coordinates, triangles) where coordinates is either a list of 2-tuples or an Mx2 Numeric array of floats representing all x, y coordinates in the mesh. triangles 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] triangles = [ [1,0,2], [1,2,3] ] #bac, bce mesh = Mesh(points, triangles) #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. """ #FIXME: Maybe rename coordinates to points (as in a poly file) #But keep 'vertex_coordinates' #FIXME: Put in check for angles less than a set minimum def __init__(self, coordinates, triangles, boundary = None, tagged_elements = None): """ Build triangles from x,y coordinates (sequence of 2-tuples or Mx2 Numeric array of floats) and triangles (sequence of 3-tuples or Nx3 Numeric array of non-negative integers). """ from Numeric import array, zeros, Int, Float, maximum, sqrt, sum General_mesh.__init__(self, coordinates, triangles) N = self.number_of_elements #Allocate space for geometric quantities self.centroid_coordinates = zeros((N, 2), Float) self.radii = zeros(N, 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) #Get x,y coordinates for all triangles and store V = self.vertex_coordinates #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.centroid_coordinates[i] = centroid #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) #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 surrogate neighbour structure self.build_surrogate_neighbour_structure() #Build boundary dictionary mapping (id, edge) to symbolic tags self.build_boundary_dictionary(boundary) #Build tagged element dictionary mapping (tag) to array of elements self.build_tagged_elements_dictionary(tagged_elements) #Update boundary indices FIXME: OBSOLETE #self.build_boundary_structure() def __repr__(self): return 'Mesh: %d triangles, %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) N = self.number_of_elements neighbourdict = {} for i in range(N): #Register all segments as keys mapping to current triangle #and segment id a = self.triangles[i, 0] b = self.triangles[i, 1] c = self.triangles[i, 2] neighbourdict[a,b] = (i, 2) #(id, edge) neighbourdict[b,c] = (i, 0) #(id, edge) neighbourdict[c,a] = (i, 1) #(id, edge) #Step 2: #Go through triangles again, but this time #reverse direction of segments and lookup neighbours. for i in range(N): a = self.triangles[i, 0] b = self.triangles[i, 1] c = self.triangles[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 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 check 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 #FIXME: This assert violates internal boundaries (delete it) #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 tag 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 #FIXME: Make this situation an error in the future #and make another function which will #enable default boundary-tags where #tags a not specified boundary[ (vol_id, edge_id) ] =\ default_boundary_tag self.boundary = boundary def build_tagged_elements_dictionary(self, tagged_elements = None): """Build the dictionary of element tags. self.tagged_elements is a dictionary of element arrays, keyed by tag: { (tag): [e1, e2, e3..] } Postconditions: self.element_tag is defined """ from Numeric import array, Int if tagged_elements is None: tagged_elements = {} else: #Check that all keys in given boundary exist for tag in tagged_elements.keys(): tagged_elements[tag] = array(tagged_elements[tag]).astype(Int) msg = 'Not all elements exist. ' assert max(tagged_elements[tag]) < self.number_of_elements, msg #print "tagged_elements", tagged_elements self.tagged_elements = tagged_elements 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) Note: If a segment is listed in the boundary dictionary it *will* become a boundary - even if there is a neighbouring triangle. This would be the case for internal boundaries """ #FIXME: Now Obsolete - maybe use some comments from here in #domain.set_boundary 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: #FIXME: One would detect internal boundaries as follows #if self.neighbours[id, edge] > -1: # print 'Internal boundary' 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 that all vertices have been registered for v in self.vertexlist: msg = 'Some points do not belong to an element.\n' msg += 'Make sure all points appear as element vertices!' if v is None: print 'WARNING (mesh.py): %s' %msg break #Should this warning be an exception? #assert v is not None, msg #Check integrity of neighbour structure for i in range(N): for v in self.triangles[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. #FIXME (Ole): In pyvolution mark 3 this is OK again #NOTE: No longer works because neighbour structure is modified by # domain set_boundary. #for id, edge in self.boundary: # assert self.neighbours[id,edge] < 0 def get_centroid_coordinates(self): """Return all centroid coordinates. Return all centroid coordinates for all triangles as an Nx2 array (ordered as x0, y0 for each triangle) """ return self.centroid_coordinates