inkscape_python_figures/plot_mesh.py at main · otvam/inkscape_python_figures · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
"""
Example for PyVista plots with the following features:
    - Selection of the camera/zoom.
    - 3D plot (field pattern and geometry).
    - 2D plot (slice of the 3D geometry).

The EM problem solution has been obtained with PyPEEC (https://pypeec.otvam.ch).
PyPEEC is a 3D quasi-magnetostatic PEEC solver (voxel-based and FFT-accelerated).
"""

__author__ = "Thomas Guillod"
__copyright__ = "Thomas Guillod - Dartmouth College"
__license__ = "Mozilla Public License Version 2.0"

import os
import pyvista as pv
import utils_pv

# ###########################################################
# Plot global parameters
# ###########################################################

# set global parameters
utils_pv.set_global()

# create the output folder
os.makedirs("render", exist_ok=True)

# ###########################################################
# Load the 3D objects
# ###########################################################

# load the EM solution (VTK file)
solution = pv.read("plot_data/mesh_solution.vtk")

# load the STL geometry
coil = pv.read("plot_data/mesh_coil.stl")
core = pv.read("plot_data/mesh_core.stl")

# scale the units
coil = coil.scale(0.001)
core = core.scale(0.001)

# extract the winding potential
potential = solution.threshold(scalars="potential")

# extract the core field
field = solution.threshold(scalars="field")

# smooth the field patterns
potential = potential.cell_data_to_point_data()
field = field.cell_data_to_point_data()

# define the colormap
cmap_potential = "viridis"
cmap_field = "magma"

# get the colormap ticks
clim_potential = (0.35, 0.65)
clim_field = (0.0, 2.0)

# ###########################################################
# Find the camera angle
# ###########################################################

# create a plotter
pl = pv.Plotter(window_size=(1024, 768))

# add the geometry
pl.add_mesh(coil, color="orange", show_edges=True)
pl.add_mesh(core, color="gray", show_edges=True)
pl.add_title("Select the camera angle")

# plot and get the final camera angle
cpos = pl.show(return_cpos=True)

# ###########################################################
# Plot the winding potential / 3D
# ###########################################################

# create a plotter
pl = pv.Plotter(window_size=(1024, 768))

# add the geometry
pl.add_mesh(potential, scalars="potential", cmap=cmap_potential, clim=clim_potential)
pl.add_mesh(core, color="gray")

# plot, save, and crop margin
pl.show(screenshot="render/mesh_winding_3D.png", cpos=cpos)
utils_pv.get_crop("render/mesh_winding_3D.png", margin=25)

# ###########################################################
# Plot the core field / 3D
# ###########################################################

# create a plotter
pl = pv.Plotter(window_size=(1024, 768))

# add the geometry
pl.add_mesh(field, scalars="field", cmap=cmap_field, clim=clim_field)
pl.add_mesh(coil, color="orange")

# plot, save, and crop margin
pl.show(screenshot="render/mesh_core_3D.png", cpos=cpos)
utils_pv.get_crop("render/mesh_core_3D.png", margin=25)

# ###########################################################
# Plot the winding potential / 2D
# ###########################################################

# create a plotter
pl = pv.Plotter(window_size=(1024, 768), lighting="none")

# create a 2D slice of the 3D solution
potential = potential.slice(normal=(0.0, 1.0, 0.0), origin=(0.0, 0.0, 0.0))
core = core.clip(normal=(0.0, 1.0, 0.0), origin=(0.0, 0.0, 0.0))

# add the geometry
pl.add_mesh(potential, scalars="potential", cmap=cmap_potential, clim=clim_potential)
pl.add_mesh(core, color="gray")

# plot, save, and crop margin
pl.enable_2d_style()
pl.enable_parallel_projection()
pl.show(screenshot="render/mesh_winding_2D.png", cpos="xz")
utils_pv.get_crop("render/mesh_winding_2D.png", margin=25)

# ###########################################################
# Plot the core field / 2D
# ###########################################################

# create a plotter
pl = pv.Plotter(window_size=(1024, 768), lighting="none")

# create a 2D slice of the 3D solution
field = field.slice(normal=(0.0, 1.0, 0.0), origin=(0.0, 0.0, 0.0))
coil = coil.clip(normal=(0.0, 1.0, 0.0), origin=(0.0, 0.0, 0.0))

# add the geometry
pl.add_mesh(field, scalars="field", cmap=cmap_field, clim=clim_field)
pl.add_mesh(coil, color="orange")

# plot, save, and crop margin
pl.enable_2d_style()
pl.enable_parallel_projection()
pl.show(screenshot="render/mesh_core_2D.png", cpos="xz")
utils_pv.get_crop("render/mesh_core_2D.png", margin=25)