Skip to content

Instantly share code, notes, and snippets.

@kingoflolz

kingoflolz/femm3.py

Created January 4, 2019 03:13
Show Gist options
  • Save kingoflolz/1c5b6960ab84f0dcc1651a842af39cdf to your computer and use it in GitHub Desktop.
Save kingoflolz/1c5b6960ab84f0dcc1651a842af39cdf to your computer and use it in GitHub Desktop.
Download ZIP
Raw
femm3.py
import random
import subprocess
import numpy as np
import scipy.optimize
import os
import time
import copy
import multiprocessing
class Simulation:
def __init__(self):
self.commands = ["newdocument(0)", 'mi_probdef(0, "millimeters", "planar", 1E-8)']
self.commands += ['mi_addboundprop("dirichlet", 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)']
self.points = {}
self.materials = {}
self.circuits = 0
self.boundaries = 0
def run_sim(self, id):
self.commands += ['mi_saveas("Untitled{}.fem")'.format(id), "mi_analyse(1)", "mi_loadsolution()"]
def save(self, varname, id):
self.commands += ['f = openfile("femmoutfile{}", "a+")'.format(id), 'write(f, {})'.format(varname),
'write(f, "\\n")', 'flush()', 'closefile(f)']
def calculate_force(self, blocks, id):
self.commands += ['mo_hidedensityplot()']
self.commands += ['mo_hidecontourplot()']
self.commands += ['mo_seteditmode("area")']
for i in blocks:
self.commands += ['mo_selectblock({}, {})'.format(i[0], i[1])]
# self.commands += ["pause()"]
self.commands += ['force = mo_blockintegral(18)',
'mo_close()']
self.save("force", id)
def zoom_useful(self):
x = []
y = []
for k, v in self.points.items():
x.append(k[0])
y.append(k[1])
self.commands += ["mo_zoom({}, {}, {}, {})".format(min(x), min(y), max(x), max(y))]
self.commands += ['mo_showdensityplot(1,0,2,0,"bmag")']
self.commands += ['mo_showcontourplot(19,-0.01,0.01,real)']
def take_cap(self, n, id):
self.commands += ['mo_savebitmap("id{}n{}.bmp")'.format(id, n)]
def move_group(self, group, x, y):
self.commands += ["mi_clearselected()"]
# self.commands += ["pause()"]
self.commands += ['mi_selectgroup({})'.format(group)]
# self.commands += ["pause()"]
self.commands += ['mi_movetranslate({}, {}, 4)'.format(x, y)]
# self.commands += ["pause()"]
self.commands += ["mi_clearselected()"]
def add_point(self, x, y, group=0):
if (x, y) not in self.points:
self.commands += ["mi_addnode({}, {})".format(x, y)]
self.points[(x, y)] = 0
def draw_line(self, x1, y1, x2, y2, group=0, propname=""):
self.add_point(x1, y1)
self.add_point(x2, y2)
self.commands += ["mi_addsegment({}, {}, {}, {})".format(x1, y1, x2, y2)]
self.commands += ["mi_selectsegment({}, {})".format((x1 + x2) / 2, (y1 + y2) / 2)]
self.commands += ['mi_setsegmentprop("{}", 0, 0, 0, {})'.format(propname, group)]
self.commands += ["mi_clearselected()"]
def rect_with_mat(self, x1, y1, x2, y2, name, mag_angle=0, current=0, group=0):
self.draw_contour(0, 0, [
(x1, y1),
(x1, y2),
(x2, y2),
(x2, y1),
], group=group)
self.set_materials((x1 + x2) / 2, (y1 + y2) / 2, name, mag_angle, current)
def draw_contour(self, x, y, l, group=0):
for iind, i in enumerate(l):
j = l[iind - 1]
self.draw_line(x + i[0], y + i[1], x + j[0], y + j[1], group=group)
def create_boundaries(self):
x = []
y = []
for k, v in self.points.items():
x.append(k[0])
y.append(k[1])
x_center = (max(x) - min(x)) / 2 + min(x)
y_center = (max(y) - min(y)) / 2 + min(y)
r = max(max(x) - min(x), max(y) - min(y))
self.commands += ["mi_makeABC(7, {}, {}, {}, {})".format(r, x_center, y_center, 0)]
def create_apb(self, x1, y1, x2, y2, x3, y3, x4, y4, group=0):
self.boundaries += 1
self.commands += ['mi_addboundprop("apb{}", 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0)'.format(self.boundaries)]
self.draw_line(x1, y1, x2, y2, propname="apb{}".format(self.boundaries), group=group)
self.draw_line(x3, y3, x4, y4, propname="apb{}".format(self.boundaries), group=group)
def create_dirichlet_line(self, x1, y1, x2, y2, group=0):
self.draw_line(x1, y1, x2, y2, propname="dirichlet".format(self.boundaries), group=group)
def add_materials(self, name):
if name not in self.materials:
self.materials[name] = 0
self.commands += ['mi_getmaterial("{}")'.format(name)]
def set_materials(self, x, y, name, mag_angle=0, current=0, group=0):
self.add_materials(name)
current_str = ""
if current != 0:
self.circuits += 1
self.commands += ['mi_addcircprop("circuit{}", {}, 1)'.format(self.circuits, current, 1)]
current_str = "circuit" + str(self.circuits)
self.commands += ["mi_clearselected()"]
self.commands += ["mi_addblocklabel({}, {})".format(x, y)]
self.commands += ["mi_selectlabel({}, {})".format(x, y)]
self.commands += ['mi_setblockprop("{}", 1, 0, "{}", {}, {}, 1)'.format(name, current_str, mag_angle, group)]
self.commands += ["mi_clearselected()"]
def output(self):
return "\n\r".join(self.commands)
def quit(self):
self.commands += ["quit()"]
class Motor:
def __init__(self, id, stator_dist=9.0):
self.teeth_height = 2
self.notch_height = 6
self.teeth_width = 8
self.notch_width = 2
self.air_gap = 0.5
self.back_iron_thickness = 3
self.stator_dist = stator_dist
self.magnet_width = self.stator_dist * 6 / 7 - 1
self.magnet_thickness = 2
self.offset = 0
self.sim = Simulation()
self.id = id
def draw_outer_magnets(self, x, a):
if a:
a1 = 90
else:
a1 = 270
self.sim.rect_with_mat(x, (self.notch_height + self.teeth_height * 2) / 2 + self.air_gap,
x + self.magnet_width,
(self.notch_height + self.teeth_height * 2) / 2 + self.air_gap + self.magnet_thickness,
"NdFeB 40 MGOe", mag_angle=a1)
self.sim.rect_with_mat(x, -(self.notch_height + self.teeth_height * 2) / 2 - self.air_gap,
x + self.magnet_width,
-(self.notch_height + self.teeth_height * 2) / 2 - self.air_gap - self.magnet_thickness,
"NdFeB 40 MGOe", mag_angle=a1)
def draw_stator(self, x, y, cur):
self.sim.draw_contour(x, y, [(self.notch_width / 2, self.notch_height / 2),
(self.notch_width / 2, -self.notch_height / 2),
(self.teeth_width / 2, -self.notch_height / 2),
(self.teeth_width / 2, -self.notch_height / 2 - self.teeth_height),
(-self.teeth_width / 2, -self.notch_height / 2 - self.teeth_height),
(-self.teeth_width / 2, -self.notch_height / 2),
(-self.notch_width / 2, -self.notch_height / 2),
(-self.notch_width / 2, self.notch_height / 2),
(-self.teeth_width / 2, self.notch_height / 2),
(-self.teeth_width / 2, self.notch_height / 2 + self.teeth_height),
(self.teeth_width / 2, self.notch_height / 2 + self.teeth_height),
(self.teeth_width / 2, self.notch_height / 2),
], group=1)
self.sim.draw_line(self.teeth_width / 2 + x, self.notch_height / 2 + y, self.teeth_width / 2 + x,
-self.notch_height / 2 + y,
group=1)
self.sim.draw_line(-self.teeth_width / 2 + x, self.notch_height / 2 + y, -self.teeth_width / 2 + x,
-self.notch_height / 2 + y,
group=1)
copper_width = (self.teeth_width - self.notch_width) / 2
current = 50 * (copper_width * self.notch_height) * cur
self.sim.set_materials(x, y, "M-22 Steel", group=1)
self.sim.set_materials(self.notch_width / 2 + 0.01 + x, y, "Copper", current=int(current), group=1)
self.sim.set_materials(-self.notch_width / 2 - 0.01 + x, y, "Copper", current=-int(current), group=1)
def build(self, steps, pic):
for iind, i in enumerate([1, -1, 0.5, -0.5, 0.5, -0.5]):
self.draw_stator(((iind - 2.5) * self.stator_dist), 0, i)
magnet_space = self.stator_dist * 6 / 7 - self.magnet_width
for i in range(7):
self.draw_outer_magnets((i - 3.5) * (self.magnet_width + magnet_space) + magnet_space / 2, i % 2)
min_x = -3.5 * (self.magnet_width + magnet_space)
max_x = 3.5 * (self.magnet_width + magnet_space)
stator_end_y = (self.notch_height + self.teeth_height * 2) / 2
back_iron_start_y = stator_end_y + self.air_gap + self.magnet_thickness
back_iron_end_y = stator_end_y + self.air_gap + self.magnet_thickness + self.back_iron_thickness
air_end_y = back_iron_end_y + 2
magnet_start_y = stator_end_y + self.air_gap
self.sim.rect_with_mat(min_x,
back_iron_start_y,
max_x,
back_iron_end_y,
"1018 Steel")
self.sim.rect_with_mat(min_x,
-back_iron_start_y,
max_x,
-back_iron_end_y,
"1018 Steel")
self.sim.set_materials(0, (self.notch_height + self.air_gap / 2) / 2 + self.teeth_height, "Air")
self.sim.set_materials(0, back_iron_end_y + 0.1, "Air")
self.sim.set_materials(0, -back_iron_end_y - 0.1, "Air")
self.sim.create_apb(min_x, back_iron_start_y, min_x, back_iron_end_y, max_x, back_iron_start_y, max_x,
back_iron_end_y)
self.sim.create_apb(min_x, -back_iron_start_y, min_x, -back_iron_end_y, max_x, -back_iron_start_y, max_x,
-back_iron_end_y)
self.sim.create_apb(min_x, magnet_start_y, min_x, stator_end_y, max_x, magnet_start_y, max_x,
stator_end_y)
self.sim.create_apb(min_x, -magnet_start_y, min_x, -stator_end_y, max_x, -magnet_start_y, max_x,
-stator_end_y)
self.sim.create_apb(min_x, magnet_start_y, min_x, back_iron_start_y, max_x, magnet_start_y, max_x,
back_iron_start_y)
self.sim.create_apb(min_x, -magnet_start_y, min_x, -back_iron_start_y, max_x, -magnet_start_y, max_x,
-back_iron_start_y)
self.sim.create_apb(min_x, air_end_y, min_x, back_iron_end_y, max_x, air_end_y, max_x, back_iron_end_y)
self.sim.create_apb(min_x, -air_end_y, min_x, -back_iron_end_y, max_x, -air_end_y, max_x, -back_iron_end_y)
self.sim.create_dirichlet_line(min_x, air_end_y, max_x, air_end_y)
self.sim.create_dirichlet_line(min_x, -air_end_y, max_x, -air_end_y)
self.sim.create_apb(min_x, stator_end_y, min_x, -stator_end_y, max_x, stator_end_y, max_x, -stator_end_y,
group=1)
b = 0
m = (magnet_space + self.magnet_width) / steps
self.sim.move_group(1, b, 0)
for i in range(steps):
k = b + m * i
self.sim.run_sim(self.id)
if pic:
print("taking image")
self.sim.zoom_useful()
self.sim.take_cap(i, self.id)
force_blocks = []
for i in range(6):
o = (i - 2.5) * self.stator_dist
force_blocks += [(o + k, 0), (o + k + self.teeth_width / 2 - 0.01, 0),
(o + k - self.teeth_width / 2 + 0.01, 0)]
self.sim.calculate_force(force_blocks, self.id)
self.sim.move_group(1, m, 0)
self.sim.quit()
open("sim{}.lua".format(self.id), mode="w+").write(self.sim.output())
def run(self):
a = subprocess.Popen(["/home/ben/.wine/drive_c/femm42/bin/femm.exe",
"-lua-script=/home/ben/PycharmProjects/axial_motor2/sim{}.lua".format(self.id)])
return a
def analyze(self):
torques = []
for line in open("femmoutfile{}".format(self.id), mode="r+").readlines():
torques.append(abs(float(line.strip())))
os.remove("femmoutfile{}".format(self.id))
os.remove("Untitled{}.fem".format(self.id))
os.remove("sim{}.lua".format(self.id))
os.remove("Untitled{}.ans".format(self.id))
return torques
def ef(f, epsilon, xk, f0, k, args):
ei = np.zeros((len(xk),), float)
ei[k] = 1.0
d = epsilon * ei
return (f(*((xk + d,) + args)) - f0) / d[k]
if __name__ == '__main__':
def fitness(vars, steps=1, pic=False):
start = time.time()
r = int(time.time() * 1000) * 1000 + random.randint(0, 1000)
time.sleep(0.01)
m1 = Motor(r)
m1.teeth_height = vars[0]
m1.notch_height = vars[1]
m1.notch_width = vars[2]
m1.air_gap = vars[3]
m1.back_iron_thickness = vars[4]
m1.magnet_thickness = vars[5]
m2 = copy.deepcopy(m1)
m2.id = int(time.time() * 1000) * 1000 + random.randint(0, 1000)
time.sleep(0.01)
m2.stator_dist = 12.251
m3 = copy.deepcopy(m2)
m3.id = int(time.time() * 1000) * 1000 + random.randint(0, 1000)
time.sleep(0.01)
m3.teeth_width = 11.215
m3.notch_width += 11.215 - 8
m3.magnet_width = 6 * 11.215 / 7 - 1
m4 = copy.deepcopy(m3)
m4.stator_dist = 15.24
m4.id = int(time.time() * 1000) * 1000 + random.randint(0, 1000)
time.sleep(0.01)
m4.stator_dist = 15.248
m1.build(steps, pic)
m2.build(steps, pic)
m3.build(steps, pic)
m4.build(steps, pic)
a = m1.run()
b = m2.run()
c = m3.run()
d = m4.run()
a.wait()
b.wait()
c.wait()
d.wait()
t1 = np.array(m1.analyze())
t2 = np.array(m2.analyze())
t3 = np.array(m3.analyze())
t4 = np.array(m4.analyze())
print(t1)
print(t2)
print(t3)
print(t4)
torque = 2 * (t1 * 3 * 0.0185 + t2 * 3 * 0.0215 + t3 * 2.8 * 0.0245 + t4 * 2.8 * 0.0245)
actual = -torque.max() / (
2 * m1.teeth_height + m1.notch_height + 2 * m1.air_gap + 2 * m1.back_iron_thickness + 2 * m1.magnet_thickness + 20)
print(vars)
print(actual)
print(torque)
print("thickness: {}".format(
2 * m1.teeth_height + m1.notch_height + 2 * m1.air_gap + 2 * m1.back_iron_thickness + 2 * m1.magnet_thickness + 20))
print("simulation took {} s".format(time.time() - start))
return actual
bounds = [(0, 8), (1, 15), (0.5, 3.5), (0.5, 2), (0.5, 3), (1, 5)]
constraints = [
scipy.optimize.LinearConstraint(np.array([0, 0, 1, 0, 0, 0]), 0.1, 4)
]
pool = multiprocessing.Pool(2)
def jac(xk, epsilon=0.01, args=(), f0=None):
"""
See ``approx_fprime``. An optional initial function value arg is added.
"""
if f0 is None:
f0 = fitness(*((xk,) + args))
grad = pool.starmap(ef, [(fitness, epsilon, xk, f0, i, args) for i in range(len(xk))])
print(grad)
return grad
start = np.array([0.8092011 , 8.52388944, 2.39536697, 0.5, 2.43616628,2.37874529])
fitness(start, steps=1, pic=True)
print(scipy.optimize.minimize(fitness, start, bounds=bounds,
constraints=constraints,
options={'eps': 0.01, "maxiter": 200}, jac=jac))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment