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