1 | # ---------------------------------------------------------------------- |
---|
2 | # EXAMPLE: Rappture <field> elements |
---|
3 | # ====================================================================== |
---|
4 | # AUTHOR: Martin Hunt, Purdue University |
---|
5 | # Copyright (c) 2015 HUBzero Foundation, LLC |
---|
6 | # |
---|
7 | # See the file "license.terms" for information on usage and |
---|
8 | # redistribution of this file, and for a DISCLAIMER OF ALL WARRANTIES. |
---|
9 | # ====================================================================== |
---|
10 | |
---|
11 | import Rappture |
---|
12 | import sys |
---|
13 | import numpy as np |
---|
14 | import sympy |
---|
15 | from sympy.abc import _clash |
---|
16 | from sympy.utilities.lambdify import lambdify |
---|
17 | |
---|
18 | # uncomment these for debugging |
---|
19 | # sys.stderr = open('tool.err', 'w') |
---|
20 | # sys.stdout = open('tool.out', 'w') |
---|
21 | |
---|
22 | rx = Rappture.PyXml(sys.argv[1]) |
---|
23 | |
---|
24 | # Rather than using eval() which is unsafe and can collide with |
---|
25 | # local variables, we use sympy to parse the symbolic expression. |
---|
26 | # lambdify turns it into a callable function. |
---|
27 | formula = rx['input.string(formula).current'].value |
---|
28 | try: |
---|
29 | formula = sympy.sympify(formula, _clash) |
---|
30 | except: |
---|
31 | raise ValueError("Formula is not a valid expression.") |
---|
32 | formula = lambdify(formula.free_symbols, formula, modules=['numpy', 'mpmath', 'sympy']) |
---|
33 | |
---|
34 | xmin, xmax = 0, 4 |
---|
35 | ymin, ymax = 0, 4 |
---|
36 | num_steps = 5 |
---|
37 | |
---|
38 | # Generate a uniform grid 2D mesh and field values... |
---|
39 | m2d = rx['output.mesh(m2d)'] |
---|
40 | m2d['about.label'] = "2D Mesh" |
---|
41 | m2d['dim'] = 2 |
---|
42 | m2d['units'] = "um" |
---|
43 | m2d['hide'] = "yes" |
---|
44 | m2d['grid.xaxis.min'] = xmin |
---|
45 | m2d['grid.xaxis.max'] = xmax |
---|
46 | m2d['grid.xaxis.numpoints'] = num_steps |
---|
47 | m2d['grid.yaxis.min'] = ymin |
---|
48 | m2d['grid.yaxis.max'] = ymax |
---|
49 | m2d['grid.yaxis.numpoints'] = num_steps |
---|
50 | |
---|
51 | f2d = rx['output.field(f2d)'] |
---|
52 | f2d['about.label'] = "2D Field" |
---|
53 | f2d['component.mesh'] = "output.mesh(m2d)" |
---|
54 | |
---|
55 | x = np.linspace(xmin, xmax, num_steps) |
---|
56 | y = np.linspace(ymin, ymax, num_steps) |
---|
57 | |
---|
58 | # Important. meshgrid defaults to matlab compatibility |
---|
59 | # mode which transposes things. Tell it to use |
---|
60 | # normal 'ij' indexing instead. Or use mgrid() |
---|
61 | xx, yy = np.meshgrid(x, y, indexing='ij') |
---|
62 | pts = formula(xx, yy, 1) |
---|
63 | f2d['component.values'] = pts |
---|
64 | |
---|
65 | vizmethod = rx['input.choice(3D).current'].value |
---|
66 | |
---|
67 | if vizmethod == "grid": |
---|
68 | # |
---|
69 | # Generate a uniform grid 3D mesh and field values... |
---|
70 | # |
---|
71 | m3d = rx['output.mesh(m3d)'] |
---|
72 | m3d['about.label'] = "3D Uniform Mesh" |
---|
73 | m3d['dim'] = 3 |
---|
74 | m3d['units'] = "um" |
---|
75 | m3d['hide'] = "yes" |
---|
76 | m3d['grid.xaxis.min'] = xmin |
---|
77 | m3d['grid.xaxis.max'] = xmax |
---|
78 | m3d['grid.xaxis.numpoints'] = 5 |
---|
79 | m3d['grid.yaxis.min'] = ymin |
---|
80 | m3d['grid.yaxis.max'] = ymax |
---|
81 | m3d['grid.yaxis.numpoints'] = 5 |
---|
82 | m3d['grid.zaxis.min'] = 0.0 |
---|
83 | m3d['grid.zaxis.max'] = 1.0 |
---|
84 | m3d['grid.zaxis.numpoints'] = 2 |
---|
85 | |
---|
86 | f3d = rx['output.field(f3d)'] |
---|
87 | f3d['about.label'] = "3D Field" |
---|
88 | f3d['component.mesh'] = "output.mesh(m3d)" |
---|
89 | |
---|
90 | xx, yy, zz = np.mgrid[xmin:xmax:5j, ymin:ymax:5j, 0:1:2j] |
---|
91 | pts = formula(xx, yy, zz) |
---|
92 | f3d['component.values'] = pts |
---|
93 | |
---|
94 | if vizmethod == "unstructured": |
---|
95 | # |
---|
96 | # Generate an unstructured grid 3D mesh and field values... |
---|
97 | # |
---|
98 | m3d = rx['output.mesh(m3d)'] |
---|
99 | |
---|
100 | m3d['about.label'] = "3D Unstructured Mesh" |
---|
101 | m3d['dim'] = 3 |
---|
102 | m3d['units'] = "um" |
---|
103 | m3d['hide'] = "yes" |
---|
104 | |
---|
105 | f3d = rx['output.field(f3d)'] |
---|
106 | f3d['about.label'] = "3D Field" |
---|
107 | f3d['component.mesh'] = "output.mesh(m3d)" |
---|
108 | |
---|
109 | x = np.linspace(xmin, xmax, 5) |
---|
110 | y = np.linspace(ymin, ymax, 5) |
---|
111 | z = np.linspace(0, 1, 2) |
---|
112 | xx, yy, zz = np.meshgrid(x, y, z, indexing='ij') |
---|
113 | |
---|
114 | # create list on index points with x changing fastest |
---|
115 | upts = '' |
---|
116 | for c in z: |
---|
117 | for b in y: |
---|
118 | for a in x: |
---|
119 | upts += "%s %s %s\n" % (a, b, c) |
---|
120 | |
---|
121 | # or you could use a list comprehension |
---|
122 | # upts = ' '.join(["%s %s %s\n" % (a, b, c) for c in z for b in y for a in x]) |
---|
123 | m3d['unstructured.points'] = upts |
---|
124 | |
---|
125 | pts = formula(xx, yy, zz) |
---|
126 | f3d['component.values'] = pts |
---|
127 | |
---|
128 | m3d['unstructured.hexahedrons'] = """ |
---|
129 | 0 1 6 5 25 26 31 30 |
---|
130 | 1 2 7 6 26 27 32 31 |
---|
131 | 2 3 8 7 27 28 33 32 |
---|
132 | 3 4 9 8 28 29 34 33 |
---|
133 | 5 6 11 10 30 31 36 35 |
---|
134 | 8 9 14 13 33 34 39 38 |
---|
135 | 10 11 16 15 35 36 41 40 |
---|
136 | 13 14 19 18 38 39 44 43 |
---|
137 | 15 16 21 20 40 41 46 45 |
---|
138 | 16 17 22 21 41 42 47 46 |
---|
139 | 17 18 23 22 42 43 48 47 |
---|
140 | 18 19 24 23 43 44 49 48 |
---|
141 | """ |
---|
142 | |
---|
143 | rx.close() |
---|