Skip to content
Snippets Groups Projects
Commit a41f13cb authored by asantosnet's avatar asantosnet
Browse files

Updated syntax

parent 0555da58
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
poisson/examples/graphene_MOS.py
+ 96
− 179
View file @ a41f13cb
'''
Example using a MOS structure
'''
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import cKDTree
import mayavi.mlab as mlp
from poisson.poisson import mesher
from poisson.poisson.discretized_poisson import (DiscretizePoisson,
get_value_from_points)
from poisson.poisson import solver
from poisson.poisson.poisson_problem import PoissonProblem
from poisson.poisson.repeated_values import remove_from_grid
def pot_points_interpolate(points_out, points_input, discretized_obj,
coordinates, deg_interpolation=[1, 1], dim=[0,1]):
regions = discretized_obj.regions
voronoi_prop = discretized_obj.voronoi_prop
pos_dirichlet = np.concatenate((regions.convex_voltage.index_of_closed_surface_points,
regions.other_voltage.index_of_closed_surface_points))
pos_dirichlet = pos_dirichlet.astype(int)
pos_neuman = regions.charge.index_of_closed_points
pos_all = np.concatenate((pos_dirichlet, pos_neuman))
# only a temporary fix. *len(pos_all) so it trows an error
# if the point we are looking at does not exists in pos_all
# which can be easily tested.
mapping = np.ones(len(discretized_obj.mesh), dtype=int)*(len(pos_all))
mapping[pos_all] = np.arange(len(pos_all))
mesh = discretized_obj.mesh[pos_all]
tree = cKDTree(mesh)
# Charge variation normalized by vaccum dielectric constant
points_charge = np.concatenate((points_out[pos_dirichlet],
points_input[pos_neuman])).astype(int)
eta_vide = 1+ 0*8.85e-12
points_charge= (points_charge / eta_vide)
charge_val = get_value_from_points(coordinates, points_charge, mapping,
pos_all, voronoi_prop, tree, vector=True,
interpolate=deg_interpolation[0])
print(charge_val[len(pos_dirichlet):-1])
# Voltage variation
points_voltage = np.concatenate((points_input[pos_dirichlet],
points_out[pos_neuman]))
voltage_val = get_value_from_points(coordinates, points_voltage,
mapping, pos_all,
voronoi_prop, tree, vector=True,
interpolate=deg_interpolation[1])
fig1 = plt.figure()
ax1 = fig1.add_subplot(111)
voltageplot = ax1.imshow(np.reshape(voltage_val, (int(np.sqrt(len(voltage_val))),
int(np.sqrt(len(voltage_val))))),
extent=(-10,10,-20,20),
interpolation='nearest', cmap='viridis')
plt.axes().set_aspect('equal', 'datalim')
plt.colorbar(voltageplot)
fig1 = plt.figure()
ax1 = fig1.add_subplot(111)
charge_plot = ax1.imshow(np.reshape(charge_val, (int(np.sqrt(len(charge_val))),
int(np.sqrt(len(charge_val))))),
extent=(-10,10,-20,20),
interpolation='nearest', cmap='viridis')
#seismic
plt.axes().set_aspect('equal', 'datalim')
plt.colorbar(charge_plot)
plt.show()
from matplotlib import pyplot as plt
from poisson.continuous import shapes
from poisson import (DiscretePoisson, GridBuilder, ContinuousGeometry,
LinearProblem)
### Geomtry
# make bridge gate
bridge = mesher.rectangle([20, 4, 2], corner=[-10, 0, 0])
bridge = shapes.Rectangle([20, 4, 2], corner=[-10, 0, 0])
bridge_bbox = [-20, 20, 0, 4, 0, 3]
# insulator1
insulator = mesher.rectangle([20,40,2], corner=[-10,-20,-2.1])
insulator = shapes.Rectangle([20,40,2], corner=[-10,-20,-2.1])
insulator_box = [-20,20,-20,20,-2,0]
# make Graphene
plane_1 = mesher.rectangle([20, 20, 2], corner=[-10, -20, -4.4])
plane_2 = mesher.rectangle([20, 20, 2], corner=[-10, 0, -4.4])
hole1 = mesher.rectangle([5, 14, 2], corner=[-10.1, -5, -4.4])
hole2 = mesher.rectangle([5, 14, 2], corner=[5.1, -5, -4.4])
hole3 = mesher.rectangle([10, 14, 2], corner=[-5, -5, -4.4])
graphene_elc1 = mesher.function_difference(plane_1, mesher.function_union(hole1,
hole2,
hole3))
graphene_elc2 = mesher.function_difference(plane_2, mesher.function_union(hole1,
hole2,
hole3))
plane_1 = shapes.Rectangle([20, 20, 2], corner=[-10, -20, -4.4])
plane_2 = shapes.Rectangle([20, 20, 2], corner=[-10, 0, -4.4])
hole1 = shapes.Rectangle([5, 14, 2], corner=[-10.1, -5, -4.4])
hole2 = shapes.Rectangle([5, 14, 2], corner=[5.1, -5, -4.4])
hole3 = shapes.Rectangle([10, 14, 2], corner=[-5, -5, -4.4])
graphene_elc1 = plane_1 - (hole1 + hole2 + hole3)
graphene_elc2 = plane_2 - (hole1 + hole2 + hole3)
graphene_bbox = [-10, 10, -20, 20, -4.5, -2]
# mesh
mesh = mesher.Mesher()
mesh.add_mesh_square(bridge_bbox, 1, bridge, 1.5)
mesh.add_mesh_square(graphene_bbox, 1, graphene_elc1, 1)
mesh.add_mesh_square(graphene_bbox, 1, graphene_elc2, 1)
mesh.add_mesh_square(insulator_box, 1, insulator, 0.5)
mesh.add_mesh_square(graphene_bbox, 1, hole1, 0.5)
mesh.add_mesh_square(graphene_bbox, 1, hole2, 0.5)
mesh.add_mesh_square(graphene_bbox, 1, hole3, 0.7)
mesh_points, l = remove_from_grid(mesh.mesh_points)
print(l)
## 3D Mesh visualization
#x, y, z = mesh.mesh_points.T
#mlp.points3d(x, y, z, mesh.point_label,
# scale_factor=0.25)
#mlp.show()
### Define Dirichelt, Neuman, Mixte and Dielectric
space = [bridge, insulator, plane_1, plane_2]
convex_dirichlet = [bridge,
graphene_elc1,
graphene_elc2]
other_dirichlet = []
neuman = [insulator,
hole1,
hole2,
hole3]
mixed_BC = []
die_co = 8.85e-12
dielectric = [mesher.wrapper_func(insulator, value=die_co),
mesher.wrapper_func(hole1, value=die_co),
mesher.wrapper_func(hole2, value=die_co),
mesher.wrapper_func(hole3, value=die_co),
mesher.wrapper_func(graphene_elc1, value=die_co),
mesher.wrapper_func(graphene_elc2, value=die_co),
mesher.wrapper_func(bridge, value=die_co)]
poissonpro = PoissonProblem(space, dielectric,
convex_dirichlet, other_dirichlet, neuman,
mixed_BC)
### Build the capacitance matrix and discretize the system
discretized_obj = DiscretizePoisson(poissonpro, mesh.mesh_points)
discretized_obj.make_voronoi()
discretized_obj.discretize(refinement_convex='Voronoi',
refinement_other='Neuman-Dirichlet')
discretized_obj.construct_capacitance_mat(distance_normalization=1)
### Define boundary conditions
dirichlet_val = [mesher.wrapper_func(bridge, value=10.0),
mesher.wrapper_func(graphene_elc1, value=5.0),
mesher.wrapper_func(graphene_elc2, value=-5.0)]
neuman_val = [mesher.wrapper_func(insulator, value=0.0),
mesher.wrapper_func(hole1, value=0),
mesher.wrapper_func(hole2, value=0),
mesher.wrapper_func(hole3, value=10e17)]
### Solve the linear system of equations
points_out, points_input = solver.solv_syst(discretized_obj,
dirichlet_val,
neuman_val)
### Plot result
# Plot result from a random set of coordinates - cut along z
z = -3.5
ymax = 20
ymin = -20
xmax = 10
xmin = -10
ypoints = np.linspace(ymin, ymax, 200)
xpoints = np.linspace(xmin, xmax, 200)
coord = np.dstack([x.flatten() for x in np.meshgrid(xpoints, ypoints, z)])[0]
pot_points_interpolate(points_out, points_input, discretized_obj,
coordinates=coord, deg_interpolation=[0, 0], dim=[0, 1])
# GRID
grid = GridBuilder()
grid.add_mesh_square(bridge_bbox, 1, bridge, 1.5)
grid.add_mesh_square(graphene_bbox, 1, graphene_elc1, 1)
grid.add_mesh_square(graphene_bbox, 1, graphene_elc2, 1)
grid.add_mesh_square(insulator_box, 0.25, insulator, 0.5)
grid.add_mesh_square(graphene_bbox, 1, hole1, 0.5)
grid.add_mesh_square(graphene_bbox, 1, hole2, 0.5)
grid.add_mesh_square(graphene_bbox, 0.1, hole3, 0.7)
# Geometry
die_insulator = 12
die_bridge = 3
die_graphene = 1.5
geometry = ContinuousGeometry(space = bridge + insulator + plane_1 + plane_2,
voltage={'bridge':bridge, 'Elec_1':graphene_elc1,
'Elec_2':graphene_elc2},
dielectric=[(insulator, 12),
(bridge, die_bridge),
(hole3, die_graphene)])
bbox = [-20,20, -20,20]
plot_geometry = geometry.plot(direction=2, bbox=bbox, plot_type='2D')
plot_geometry(variable=0)
plot_geometry(variable=-2)
plot_geometry(variable=4)
plot_geometry(variable=-4.5)
geometry.plot(points=grid.points, plot_type='3D', scale_factor = 0.2)
plt.show()
# Discrete system
discrete_poi = DiscretePoisson(
geometry, grid=grid,
selection = {'Neuman-Dirichlet':[['voltage', '*']]})
# Build the linear problem and solve it
linear_prob_inst = LinearProblem(discrete_poi,
voltage_val = [(bridge, 10.0),
(graphene_elc1, 5.0),
(graphene_elc2, -5.0)],
charge_val = [(hole3, 10e7)],
is_charge_density=True)
# Save data to vtk file
current = '/'.join((os.path.dirname(os.path.abspath(__file__)),
'example_sphere3D'))
linear_prob_inst.save_to_vtk(filename=current)
#### Plotting
# Plot 2D cut
plot_volt, plot_charge = linear_prob_inst.plot_cut_2d(direction=2,
npoints=(1000,1000))
plot_volt(variable=0, colorbar_label='Voltage (V)')
plot_charge(variable=0, colorbar_label='Charge')
plt.show()
plot_voltage, plot_charge = linear_prob_inst.plot_cut_1d(
directions=(1, 2), bbox='default', npoints=2000)
fig2 = plt.figure()
ax_voltage = fig2.add_subplot(111)
t, data_volt = plot_voltage(
(0, 0), ax = ax_voltage, marker='.', color='k',
label='Voltage simulation', linestyle='None')
ax_voltage.set_xlabel('{0}(nm)'.format('y'))
ax_voltage.set_ylabel('Voltage(V)')
ax_voltage.legend(loc='upper center')
ax_charge = ax_voltage.twinx()
tt, data_charge = plot_charge(
(0, 0), ax = ax_charge, marker='.', color='b',
label='Charge simulation', linestyle='None')
ax_charge.set_xlabel('{0}(nm)'.format('y'))
ax_charge.set_ylabel(r'Charge density $(\#.nm^{{{-2}}})$')
ax_charge.legend(loc='lower center')
plt.show()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment