1#!/usr/bin/python 2"""2dreantrantcorner-generator.py -- generate a mesh for the reentrant corner 3 4SYNOPSIS: 5./2dreantrantcorner-generator.py [-s] [-t TRIANGLE_OPTS] NODES OPENING 6 7PARAMETERS: 8-s Call showme on the generated mesh 9 10NODES Number of nodes on the circumference 11 12OPENING the opening angle in degrees 13""" 14 15from math import pi, cos, sin 16from os import system, chdir 17from shutil import rmtree 18from sys import argv, stderr, stdout 19from tempfile import mkdtemp 20 21 22showme = False 23triangle_opts = "" 24params = [] 25i = 1 26while i < len(argv): 27 if argv[i] == "-s": 28 showme = True 29 elif argv[i] == "-t": 30 i+=1 31 triangle_opts = argv[i] 32 elif argv[i][0:2] == "-t": 33 triangle_opts = argv[i][2:] 34 elif argv[i][0] == "-": 35 print __doc__ 36 exit(1) 37 else: 38 params += [argv[i]] 39 i+=1 40 41if(len(params) != 2): 42 print __doc__ 43 exit(1) 44 45nodes = int(params[0]) 46angle = float(params[1])/360*2*pi 47 48dir = mkdtemp() 49chdir(dir) 50stderr.write("Using Temporary directory %s\n" % dir) 51 52poly = file("mesh.poly", "w") 53#First line: <# of vertices> <dimension (must be 2)> <# of attributes> <# of boundary markers (0 or 1)> 54poly.write("%d 2 0 0\n" % (nodes+1)) 55#Following lines: <vertex #> <x> <y> [attributes] [boundary marker] 56for i in range(nodes): 57 used = 2*pi - angle 58 piece_angle = used / (nodes-1) 59 my_theta = piece_angle * i 60 x = cos(my_theta) 61 y = sin(my_theta) 62 poly.write("%d %.17g %.17g\n" % (i, x, y)) 63poly.write("%d 0 0\n" % nodes) 64#One line: <# of segments> <# of boundary markers (0 or 1)> 65poly.write("%d 0\n" % (nodes+1)) 66#Following lines: <segment #> <endpoint> <endpoint> [boundary marker] 67for i in range(nodes+1): 68 poly.write("%d %g %g\n" % (i, i, (i+1)%(nodes+1))) 69#One line: <# of holes> 70poly.write("0\n") 71#Following lines: <hole #> <x> <y> 72#Optional line: <# of regional attributes and/or area constraints> 73#Optional following lines: <region #> <x> <y> <attribute> <max area> 74poly.close() 75 76system("triangle -p %s mesh.poly 1>&2" % triangle_opts) 77 78if showme: 79 system("showme mesh.1.ele 1>&2") 80 81stdout.write("DGF\n") 82stdout.write("%% generated by %s\n" % " ".join(argv)) 83 84node = file("mesh.1.node") 85nodes = node.readline().split()[0]; 86stdout.write("Vertex\n") 87for i in range(int(nodes)): 88 (x, y) = node.readline().split()[1:3] 89 stdout.write("%s\t%s\n" % (x, y)) 90node.close() 91stdout.write("#\n") 92 93ele = file("mesh.1.ele") 94elems = ele.readline().split()[0] 95stdout.write("Simplex\n") 96for i in range(int(elems)): 97 nodes = ele.readline().split()[1:4] 98 stdout.write("\t".join(nodes)+"\n") 99ele.close() 100stdout.write("#\n") 101 102stdout.write("Boundarydomain\n") 103stdout.write("default 1\n") 104stdout.write("#\n") 105 106rmtree(dir) 107