"""Library of standard meshes and facilities for reading various mesh file formats """ def rectangular(m, n, len1=1.0, len2=1.0, origin = (0.0, 0.0)): """Setup a rectangular grid of triangles with m+1 by n+1 grid points and side lengths len1, len2. If side lengths are omitted the mesh defaults to the unit square. len1: x direction (left to right) len2: y direction (bottom to top) Return to lists: points and elements suitable for creating a Mesh or FVMesh object, e.g. Mesh(points, elements) """ from config import epsilon #E = m*n*2 #Number of triangular elements #P = (m+1)*(n+1) #Number of initial vertices delta1 = float(len1)/m delta2 = float(len2)/n #Dictionary of vertex objects vertices = {} points = [] for i in range(m+1): for j in range(n+1): vertices[i,j] = len(points) points.append([i*delta1 + origin[0], j*delta2 + origin[1]]) #Construct 2 triangles per rectangular element and assign tags to boundary elements = [] boundary = {} for i in range(m): for j in range(n): v1 = vertices[i,j+1] v2 = vertices[i,j] v3 = vertices[i+1,j+1] v4 = vertices[i+1,j] #Update boundary dictionary and create elements if i == m-1: boundary[(len(elements), 2)] = 'right' if j == 0: boundary[(len(elements), 1)] = 'bottom' elements.append([v4,v3,v2]) #Lower element if i == 0: boundary[(len(elements), 2)] = 'left' if j == n-1: boundary[(len(elements), 1)] = 'top' elements.append([v1,v2,v3]) #Upper element return points, elements, boundary def oblique(m, n, lenx = 1.0, leny = 1.0, theta = 8.95, origin = (0.0, 0.0)): """Setup a oblique grid of triangles with m segments in the x-direction and n segments in the y-direction """ from Numeric import array import math from config import epsilon deltax = lenx/float(m) deltay = leny/float(n) a = 0.75*lenx*math.tan(theta/180.*math.pi) x1 = lenx y1 = 0 x2 = lenx y2 = leny x3 = 0.25*lenx y3 = leny x4 = x3 y4 = 0 a2 = a/(x1-x4) a1 = -a2*x4 a4 = ((a1 + a2*x3)/y3-(a1 + a2*x2)/y2)/(x2-x3) a3 = 1. - (a1 + a2*x3)/y3 - a4*x3 # Dictionary of vertex objects vertices = {} points = [] for i in range(m+1): x = deltax*i for j in range(n+1): y = deltay*j if x > 0.25*lenx: y = a1 + a2*x + a3*y + a4*x*y vertices[i,j] = len(points) points.append([x + origin[0], y + origin[1]]) # Construct 2 triangles per element elements = [] boundary = {} for i in range(m): for j in range(n): v1 = vertices[i,j+1] v2 = vertices[i,j] v3 = vertices[i+1,j+1] v4 = vertices[i+1,j] #Update boundary dictionary and create elements if i == m-1: boundary[(len(elements), 2)] = 'right' if j == 0: boundary[(len(elements), 1)] = 'bottom' elements.append([v4,v3,v2]) #Lower if i == 0: boundary[(len(elements), 2)] = 'left' if j == n-1: boundary[(len(elements), 1)] = 'top' elements.append([v1,v2,v3]) #Upper return points, elements, boundary #OLD FORMULA FOR SETTING BOUNDARY TAGS - OBSOLETE NOW # for id, face in M.boundary: # e = element_class.instances[id] # x0, y0, x1, y1, x2, y2 = e.get_instance_vertex_coordinates() # if face==2: # #Left or right# # if abs(x0-origin[0]) < epsilon: # M.boundary[(id,face)] = 'left' # elif abs(origin[0]+lenx-x0) < epsilon: # M.boundary[(id,face)] = 'right' # else: # print face, id, id%m, m, n # raise 'Left or Right Unknown boundary' # elif face==1: # #Top or bottom # if x0 > 0.25*lenx and abs(y0-a1-a2*x0-origin[1]) < epsilon or\ # x0 <= 0.25*lenx and abs(y0-origin[1]) < epsilon: # M.boundary[(id,face)] = 'bottom' # elif abs(origin[1]+leny-y0) < epsilon: # M.boundary[(id,face)] = 'top' # else: # print face, id, id%m, m, n # raise 'Top or Bottom Unknown boundary' # else: # print face, id, id%m, m, n # raise 'Unknown boundary' # # return M def circular(m, n, radius=1.0, center = (0.0, 0.0)): """Setup a circular grid of triangles with m concentric circles and with n radial segments. If radius is are omitted the mesh defaults to the unit circle radius. radius: radius of circle #FIXME: The triangles become degenerate for large values of m or n. """ from math import pi, cos, sin radius = float(radius) #Ensure floating point format #Dictionary of vertex objects and list of points vertices = {} points = [[0.0, 0.0]] #Center point vertices[0, 0] = 0 for i in range(n): theta = 2*i*pi/n x = cos(theta) y = sin(theta) for j in range(1,m+1): delta = j*radius/m vertices[i,j] = len(points) points.append([delta*x, delta*y]) #Construct 2 triangles per element elements = [] for i in range(n): for j in range(1,m): i1 = (i + 1) % n #Wrap around v1 = vertices[i,j+1] v2 = vertices[i,j] v3 = vertices[i1,j+1] v4 = vertices[i1,j] elements.append([v4,v2,v3]) #Lower elements.append([v1,v3,v2]) #Upper #Do the center v1 = vertices[0,0] for i in range(n): i1 = (i + 1) % n #Wrap around v2 = vertices[i,1] v3 = vertices[i1,1] elements.append([v1,v2,v3]) #center return points, elements def from_polyfile(name): """Read mesh from .poly file, an obj like file format listing first vertex coordinates and then connectivity """ from util import anglediff from math import pi import os.path root, ext = os.path.splitext(name) if ext == 'poly': filename = name else: filename = name + '.poly' fid = open(filename) points = [] #x, y values = [] #z ##vertex_values = [] #Repeated z triangles = [] #v0, v1, v2 lines = fid.readlines() keyword = lines[0].strip() msg = 'First line in .poly file must contain the keyword: POINTS' assert keyword == 'POINTS', msg offending = 0 i = 1 while keyword == 'POINTS': line = lines[i].strip() i += 1 if line == 'POLYS': keyword = line break fields = line.split(':') assert int(fields[0]) == i-1, 'Point indices not consecutive' #Split the three floats xyz = fields[1].split() x = float(xyz[0]) y = float(xyz[1]) z = float(xyz[2]) points.append([x, y]) values.append(z) k = i while keyword == 'POLYS': line = lines[i].strip() i += 1 if line == 'END': keyword = line break fields = line.split(':') assert int(fields[0]) == i-k, 'Poly indices not consecutive' #Split the three indices vvv = fields[1].split() i0 = int(vvv[0])-1 i1 = int(vvv[1])-1 i2 = int(vvv[2])-1 #Check for and exclude degenerate areas x0 = points[i0][0] y0 = points[i0][1] x1 = points[i1][0] y1 = points[i1][1] x2 = points[i2][0] y2 = points[i2][1] area = abs((x1*y0-x0*y1)+(x2*y1-x1*y2)+(x0*y2-x2*y0))/2 if area > 0: #Ensure 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) if a0 < pi and a1 < pi and a2 < pi: #all is well j0 = i0 j1 = i1 j2 = i2 else: #Swap two vertices j0 = i1 j1 = i0 j2 = i2 triangles.append([j0, j1, j2]) ##vertex_values.append([values[j0], values[j1], values[j2]]) else: offending +=1 print 'Removed %d offending triangles out of %d' %(offending, len(lines)) return points, triangles, values def strang_mesh(filename): """Read Strang generated mesh. """ from math import pi from util import anglediff fid = open(filename) points = [] # List of x, y coordinates triangles = [] # List of vertex ids as listed in the file for line in fid.readlines(): fields = line.split() if len(fields) == 2: # we are reading vertex coordinates points.append([float(fields[0]), float(fields[1])]) elif len(fields) == 3: # we are reading triangle point id's (format ae+b) triangles.append([int(float(fields[0]))-1, int(float(fields[1]))-1, int(float(fields[2]))-1]) else: raise 'wrong format in ' + filename elements = [] #Final list of elements for t in triangles: #Get vertex coordinates v0 = t[0] v1 = t[1] v2 = t[2] x0 = points[v0][0] y0 = points[v0][1] x1 = points[v1][0] y1 = points[v1][1] x2 = points[v2][0] y2 = points[v2][1] #Check that points are arranged in counter clock-wise order vec0 = [x1-x0, y1-y0] vec1 = [x2-x1, y2-y1] vec2 = [x0-x2, y0-y2] a0 = anglediff(vec1, vec0) a1 = anglediff(vec2, vec1) a2 = anglediff(vec0, vec2) if a0 < pi and a1 < pi and a2 < pi: elements.append([v0, v1, v2]) else: elements.append([v0, v2, v1]) return points, elements # def contracting_channel_mesh(m, n, x1 = 0.0, x2 = 1./3., x3 = 2./3., x4 = 1.0, # y1 =0.0, y4 = -1./4., y8 = 1.0, # origin = (0.0, 0.0), point_class=Point, # element_class=Triangle, mesh_class=Mesh): # """Setup a oblique grid of triangles # with m segments in the x-direction and n segments in the y-direction # Triangle refers to the actual class or subclass to be instantiated: # e.g. if Volume is a subclass of Triangle, # this function can be invoked with the keywords # oblique_mesh(...,Triangle=Volume, Mesh=Domain) # """ # from Numeric import array # from visual import rate# # import math # from config import epsilon # E = m*n*2 #Number of triangular elements # P = (m+1)*(n+1) #Number of initial vertices # initialise_consecutive_datastructure(P+E, E, # point_class, # element_class, # mesh_class) # deltax = (x4 - x1)/float(m) # deltay = (y8 - y1)/float(n) # a = y4 - y1 # if y8 - a <= y1 + a: # print a,y1,y4 # raise 'Constriction is too large reduce a' # y2 = y1 # y3 = y4 # x5 = x4 # y5 = y8 - a # x6 = x3 # y6 = y5 # x7 = x2 # y7 = y8 # x8 = x1 # a2 = a/(x3 - x2) # a1 = a - a*x3/(x3 - x2) # a4 = (-a + a2*(x7 - x6))/(x6 - x7)/y7 # a3 = (y7 - a1 - x7*a2 - a4*x7*y7)/y7 # # Dictionary of vertex objects # vertices = {} # for i in range(m+1): # x = deltax*i # for j in range(n+1): # y = deltay*j # if x > x2: # if x < x3: # y = a1 + a2*x + a3*y + a4*x*y # else: # y = a + y*(y5 - y4)/(y8 - y1) # vertices[i,j] = Point(x + origin[0],y + origin[1]) # # Construct 2 elements per element # elements = [] # for i in range(m): # for j in range(n): # v1 = vertices[i,j+1] # v2 = vertices[i,j] # v3 = vertices[i+1,j+1] # v4 = vertices[i+1,j] # elements.append(element_class(v4,v3,v2)) #Lower # elements.append(element_class(v1,v2,v3)) #Upper # M = mesh_class(elements) # #Set a default tagging # for id, face in M.boundary: # e = element_class.instances[id] # x_0, y_0, x_1, y_1, x_2, y_2 = e.get_instance_vertex_coordinates() # lenx = x4 - x1 # if face==2: # #Left or right# # if abs(x_0-origin[0]) < epsilon: # M.boundary[(id,face)] = 'left' # elif abs(origin[0]+lenx-x_0) < epsilon: # M.boundary[(id,face)] = 'right' # else: # print face, id, id%m, m, n # raise 'Left or Right Unknown boundary' # elif face==1: # #Top or bottom # if x_0 <= x2 and abs(y_0-y1-origin[1]) < epsilon or\ # x_0 > x3 and abs(y_0-y3-origin[1]) < epsilon or\ # x_0 > x2 and x_0 <= x3 and abs(y_0-(y2+(y3-y2)*(x_0-x2)/(x3-x2)+origin[1])) < epsilon: # M.boundary[(id,face)] = 'bottom' # elif x_0 <= x7 and abs(y_0-y8-origin[1]) < epsilon or\ # x_0 > x6 and abs(y_0-y6-origin[1]) < epsilon or\ # x_0 > x7 and x_0 <= x6 and abs(y_0-(y7+(y6-y7)*(x_0-x7)/(x6-x7)+origin[1])) < epsilon: # M.boundary[(id,face)] = 'top' # else: # print face, id, id%m, m, n # raise 'Top or Bottom Unknown boundary' # else: # print face, id, id%m, m, n # raise 'Unknown boundary' # # print id, face, M.boundary[(id,face)],x_0,y_0,x_1,y_1,x_2,y_2 # return M # #Map from edge number to indices of associated vertices # edge_map = ((1,2), (0,2), (0,1))