| 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() |
|---|
