diff --git a/Python FEA Simulations/__pycache__/pitFinalMaterials.cpython-38.pyc b/Python FEA Simulations/__pycache__/pitFinalMaterials.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..02d6d0709fdf550fcfb41ca80ed35936e6d9b749
Binary files /dev/null and b/Python FEA Simulations/__pycache__/pitFinalMaterials.cpython-38.pyc differ
diff --git a/Python FEA Simulations/__pycache__/pitFinalPhysics.cpython-38.pyc b/Python FEA Simulations/__pycache__/pitFinalPhysics.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..796df3cb7fa33ff5b882f3e4ff31be96fde30257
Binary files /dev/null and b/Python FEA Simulations/__pycache__/pitFinalPhysics.cpython-38.pyc differ
diff --git a/Python FEA Simulations/old/pitFinalFEA.py b/Python FEA Simulations/old/pitFinalFEA.py
new file mode 100644
index 0000000000000000000000000000000000000000..7ccf8064383dcbf360e967c7239fea3def4147fc
--- /dev/null
+++ b/Python FEA Simulations/old/pitFinalFEA.py
@@ -0,0 +1,116 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.patches as patches
+import matplotlib.colors as colors
+
+# PHYSICAL PARAMETERS
+u0 = 4e-7*np.pi
+m = np.array([0,1]) # every dipole will have the same magnetization vector, Bohn Magneton is 9.274e-24 J/T for context
+Idl = np.array((0,0,-10)) # current flow vector for biot-savart
+
+def calc_biot_H(r):
+ r_mag = np.sqrt(r.dot(r))
+ r_hat = r/r_mag
+ biot_H = 1/(4*np.pi*r_mag**2)*(np.cross(Idl,r_hat))
+ #print(biot_H)
+ return biot_H # numpy will automatically resolve the z componenet, so no need to delete it or append it to r_hat initially
+
+def calc_dipo_H(r, internalCheck):
+ r_mag = np.sqrt(r.dot(r)) # the magnitude of r (faster than np.linalg.norm(r))
+ r_hat = r/r_mag # unit vector is the vector divided by its magnitude
+ dipo_H = 1/(4*np.pi*r_mag**3)*(3*np.dot(r_hat,m)*r_hat - m)
+ if internalCheck == 4: # dirac delta is essentially all of the if statements that led us here...
+ dipo_H = dipo_H - 1/(3)*m # ***artificially swapped sign of this term
+ return dipo_H
+
+if __name__ == '__main__':
+# BUILD 2D SIMULATION
+ x_size = 30
+ y_size = 30
+
+ elements = np.zeros((x_size,y_size,4))
+ '''
+ elements[0] is a tristate where 0=air, 1=dipole, and 2=wire
+ elements[1] is the x component of the H field
+ elements[2] is the y component of the H field
+ elements[3] is the magnitude of that H vector
+ '''
+
+ # Now let's assign a region of permanent magnet
+ mag_w = 6
+ mag_h = 10
+ mag_o = [int(x_size/2-mag_w/2), int(y_size/2-mag_h/2)]
+ elements[mag_o[0]:mag_o[0]+mag_w, mag_o[1]:mag_o[1]+mag_h] = [1,0,0,0]
+
+ # Now let's assign a region of current carrying wire
+ wire_w = 2
+ wire_h = 2
+ wire_o = [6, int(x_size/2)]#[int(x_size/2-wire_w/2), int(y_size/2-wire_h/2)]
+ elements[wire_o[0]:wire_o[0]+wire_w, wire_o[1]:wire_o[1]+wire_h] = [2,0,0,0]
+
+# Now let's iterate through every element and for each one, find the impact of every PM in the array
+# Switch to make sure you are only checking magnet cells
+ for i in range(len(elements)):
+ for j in range(len(elements[i])):
+ pos = np.array([i,j])
+ for im in range(len(elements)):
+ for jm in range(len(elements[i])):
+ materialSpecial = elements[im][jm][0]
+ if materialSpecial >= 1 and (i,j) != (im,jm):
+ m_pos = np.array([im,jm])
+ r = pos - m_pos
+ # First determine the vector r pointing from the dipole to the current pos
+ internalCheck = 0
+ if i >= mag_o[0] and i <= mag_o[0]+mag_w-1: # if we are inside the magnet
+ internalCheck += 1
+ if j >= mag_o[1] and j <= mag_o[1]+mag_h-1:
+ internalCheck += 1
+ if im >= mag_o[0] and im <= mag_o[0]+mag_w-1: # if we are inside the magnet
+ internalCheck += 1
+ if jm >= mag_o[1] and jm <= mag_o[1]+mag_h-1:
+ internalCheck += 1
+ if materialSpecial == 1:
+ temp = calc_dipo_H(r,internalCheck)
+ elif materialSpecial == 2:
+ temp = calc_biot_H(r)
+ #temp = calc_biot_H(r)
+ elements[i,j][1] += temp[1] # artificially swapped order here to make plotting work
+ elements[i,j][2] += temp[0]
+# Now let's convert those mx and my vectors to magnitude for easier plotting
+ magnitudes = np.zeros((x_size,y_size))
+ for i in range(len(elements)):
+ for j in range(len(elements[i])):
+ mag_vector = np.array([elements[i,j][1],elements[i,j][2]])
+ elements[i,j][3] = np.sqrt(mag_vector.dot(mag_vector))
+ magnitudes[i][j] = elements[i,j][3]
+
+# PLOTTING
+ fig, ax = plt.subplots()
+ #plt.imshow(magnitudes)#, interpolation='nearest')
+ x,y = np.meshgrid(np.linspace(0,x_size,x_size),np.linspace(0,y_size,y_size))
+ u = np.zeros((x_size,y_size))
+ v = np.zeros((x_size,y_size))
+ m = np.zeros((x_size,y_size))
+ for i in range(len(elements)):
+ for j in range(len(elements[i])):
+ r = np.array([elements[i,j][1],elements[i,j][2]])
+ u[i,j] = elements[i,j][1]
+ v[i,j] = elements[i,j][2]
+ m[i,j] = elements[i,j][3]
+ # Normalize here, only applied to the quiver vectors
+ u[i,j] = u[i,j]/(m[i,j]**.95) # name the exponent <1 to scale favorably
+ v[i,j] = v[i,j]/(m[i,j]**.95)
+ m_log = np.log(m)
+ qq=plt.quiver(x,y,u,v,m,cmap=plt.cm.jet,norm=colors.LogNorm(vmin=m.min(), vmax=m.max()),width=.0025) # magnitude is still being used for color which is why this works
+ plt.colorbar(qq, cmap=plt.cm.jet)
+
+ # Add borders for Magnets and Wires
+ mag1 = patches.Rectangle((np.flip(mag_o)), mag_h, mag_w, linewidth=1, edgecolor='r', facecolor='none')
+ wire1 = patches.Rectangle((np.flip(wire_o-np.array([.25,.0]))), wire_h, wire_w, linewidth=1, edgecolor='g', facecolor='none')
+ ax.add_patch(mag1)
+ ax.add_patch(wire1)
+ #ax = plt.gca()
+ #ax.set_ylim(ax.get_ylim()[::-1])
+ #plt.gca().invert_yaxis()
+ plt.show()
\ No newline at end of file
diff --git a/Python FEA Simulations/old/pitFinalFEAold.py b/Python FEA Simulations/old/pitFinalFEAold.py
new file mode 100644
index 0000000000000000000000000000000000000000..ae937184d8eb62d81f34651f5180ccf9a8acd0c9
--- /dev/null
+++ b/Python FEA Simulations/old/pitFinalFEAold.py
@@ -0,0 +1,117 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.patches as patches
+import matplotlib.colors as colors
+
+# PHYSICAL PARAMETERS
+u0 = 4e-7*np.pi
+m = np.array([0,1]) # every dipole will have the same magnetization vector, Bohn Magneton is 9.274e-24 J/T for context
+Idl = np.array((0,0,-10)) # current flow vector for biot-savart
+
+def calc_biot_H(r):
+ r_mag = np.sqrt(r.dot(r))
+ r_hat = r/r_mag
+ biot_H = 1/(4*np.pi*r_mag**2)*(np.cross(Idl,r_hat))
+ #print(biot_H)
+ return biot_H # numpy will automatically resolve the z componenet, so no need to delete it or append it to r_hat initially
+
+def calc_dipo_H(r, internalCheck):
+ r_mag = np.sqrt(r.dot(r)) # the magnitude of r (faster than np.linalg.norm(r))
+ r_hat = r/r_mag # unit vector is the vector divided by its magnitude
+ dipo_H = 1/(4*np.pi*r_mag**3)*(3*np.dot(r_hat,m)*r_hat - m)
+ if internalCheck == 4: # dirac delta is essentially all of the if statements that led us here...
+ dipo_H = dipo_H - 1/(3)*m # ***artificially swapped sign of this term
+ return dipo_H
+
+if __name__ == '__main__':
+# BUILD 2D SIMULATION
+ x_size = 5
+ y_size = 5
+
+ elements = np.zeros((x_size,y_size,4))
+ '''
+ elements[0] is a tristate where 0=air, 1=dipole, and 2=wire
+ elements[1] is the x component of the H field
+ elements[2] is the y component of the H field
+ elements[3] is the magnitude of that H vector
+ '''
+
+ # Now let's assign a region of permanent magnet
+ mag_w = 6
+ mag_h = 10
+ mag_o = [int(x_size/2-mag_w/2), int(y_size/2-mag_h/2)]
+ #elements[mag_o[0]:mag_o[0]+mag_w, mag_o[1]:mag_o[1]+mag_h] = [1,0,0,0]
+
+ # Now let's assign a region of current carrying wire
+ wire_w = 2
+ wire_h = 2
+ wire_o = [2, int(x_size/2)]#[int(x_size/2-wire_w/2), int(y_size/2-wire_h/2)]
+ elements[wire_o[0]:wire_o[0]+wire_w, wire_o[1]:wire_o[1]+wire_h] = [2,0,0,0]
+
+# Now let's iterate through every element and for each one, find the impact of every PM in the array
+# Switch to make sure you are only checking magnet cells
+ for i in range(len(elements)):
+ for j in range(len(elements[i])):
+ pos = np.array([i,j])
+ for im in range(len(elements)):
+ for jm in range(len(elements[i])):
+ materialSpecial = elements[im][jm][0]
+ if materialSpecial >= 1 and (i,j) != (im,jm):
+ m_pos = np.array([im,jm])
+ r = pos - m_pos
+ print(r)
+ # First determine the vector r pointing from the dipole to the current pos
+ internalCheck = 0
+ if i >= mag_o[0] and i <= mag_o[0]+mag_w-1: # if we are inside the magnet
+ internalCheck += 1
+ if j >= mag_o[1] and j <= mag_o[1]+mag_h-1:
+ internalCheck += 1
+ if im >= mag_o[0] and im <= mag_o[0]+mag_w-1: # if we are inside the magnet
+ internalCheck += 1
+ if jm >= mag_o[1] and jm <= mag_o[1]+mag_h-1:
+ internalCheck += 1
+ if materialSpecial == 1:
+ temp = calc_dipo_H(r,internalCheck)
+ elif materialSpecial == 2:
+ temp = calc_biot_H(r)
+ #temp = calc_biot_H(r)
+ elements[i,j][1] += temp[1] # artificially swapped order here to make plotting work
+ elements[i,j][2] += temp[0]
+# Now let's convert those mx and my vectors to magnitude for easier plotting
+ magnitudes = np.zeros((x_size,y_size))
+ for i in range(len(elements)):
+ for j in range(len(elements[i])):
+ mag_vector = np.array([elements[i,j][1],elements[i,j][2]])
+ elements[i,j][3] = np.sqrt(mag_vector.dot(mag_vector))
+ magnitudes[i][j] = elements[i,j][3]
+
+# PLOTTING
+ fig, ax = plt.subplots()
+ #plt.imshow(magnitudes)#, interpolation='nearest')
+ x,y = np.meshgrid(np.linspace(0,x_size,x_size),np.linspace(0,y_size,y_size))
+ u = np.zeros((x_size,y_size))
+ v = np.zeros((x_size,y_size))
+ m = np.zeros((x_size,y_size))
+ for i in range(len(elements)):
+ for j in range(len(elements[i])):
+ r = np.array([elements[i,j][1],elements[i,j][2]])
+ u[i,j] = elements[i,j][1]
+ v[i,j] = elements[i,j][2]
+ m[i,j] = elements[i,j][3]
+ # Normalize here, only applied to the quiver vectors
+ u[i,j] = u[i,j]/(m[i,j]**.95) # name the exponent <1 to scale favorably
+ v[i,j] = v[i,j]/(m[i,j]**.95)
+ m_log = np.log(m)
+ qq=plt.quiver(x,y,u,v,m,cmap=plt.cm.jet,norm=colors.LogNorm(vmin=m.min(), vmax=m.max()),width=.0025) # magnitude is still being used for color which is why this works
+ plt.colorbar(qq, cmap=plt.cm.jet)
+
+ # Add borders for Magnets and Wires
+ mag1 = patches.Rectangle((np.flip(mag_o)), mag_h, mag_w, linewidth=1, edgecolor='r', facecolor='none')
+ wire1 = patches.Rectangle((np.flip(wire_o-np.array([.25,.0]))), wire_h, wire_w, linewidth=1, edgecolor='g', facecolor='none')
+ ax.add_patch(mag1)
+ ax.add_patch(wire1)
+ #ax = plt.gca()
+ #ax.set_ylim(ax.get_ylim()[::-1])
+ #plt.gca().invert_yaxis()
+ plt.show()
\ No newline at end of file
diff --git a/Python FEA Simulations/pitFinalFEA.py b/Python FEA Simulations/pitFinalFEA.py
new file mode 100644
index 0000000000000000000000000000000000000000..cd7d72d91181094c7b3ec12834464704f7a1a3f1
--- /dev/null
+++ b/Python FEA Simulations/pitFinalFEA.py
@@ -0,0 +1,118 @@
+from pickle import FALSE
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import matplotlib.patches as patches
+import matplotlib.colors as colors
+from pitFinalMaterials import *
+from pitFinalPhysics import *
+
+out_path = "output/magnet_b_field"
+iterations = 10
+
+# PHYSICAL PARAMETERS
+u0 = 4e-7*np.pi
+u_steel = u0*5000 # relative permeability of steel = 5000, permeability = 5000*u0
+ur = u_steel / u0
+Xm = ur - 1
+
+ferriteON = False
+
+if __name__ == '__main__':
+# BUILD 2D SIMULATION
+ x_size = 31
+ y_size = 31
+ elements = np.zeros((x_size,y_size,4))
+ ''' elements[0] is a tristate where 0=air, 1=dipole, and 2=wire (no longer used)
+ elements[1] is the x component of the H field
+ elements[2] is the y component of the H field
+ elements[3] is the magnitude of that H vector '''
+ magnets, wires, ferritics = materials(x_size, y_size)
+
+ # Now let's iterate through every element and for each one, find the impact of every PM in the array
+ for i in range(x_size):
+ for j in range(y_size):
+ pos = np.array([i,j])
+ # Here we want to calculate the effect of an entire magnet and/or an entire wire
+ for wire in wires:
+ wireElements = wire[4]
+ for wireElement in wireElements:
+ r = pos - wireElement # r is the distanc between the current point [i,j] and the wire element
+ if np.linalg.norm(r) != 0: # don't try to caclulate for r = 0 (there's no curl contributed)
+ elements[i,j][1] += calc_biot_B(wire[0], r)[0]
+ elements[i,j][2] += calc_biot_B(wire[0], r)[1]
+ for magnet in magnets:
+ ox, oy = magnet[3] # center of the magnet
+ x_pos, y_pos = pos - [ox,oy] # current x and y positions
+ mx, my = magnet[0] # magnetization vector
+ hx, hy = magnet[1], magnet[2] # height of the magnet
+ elements[i,j][1] += calc_magnet_B(x_pos, y_pos, mx, my, hx, hy, ox, oy)[0]
+ elements[i,j][2] += calc_magnet_B(x_pos, y_pos, mx, my, hx, hy, ox, oy)[1]
+
+ # Now lets store all of our initially induced B fields permanenetly
+ B_induced = elements
+
+ if len(ferritics) == 0:
+ iterations = 1
+ for iteration in range(iterations):
+ if ferriteON == True:
+ # Now let's assign our ferritic regions and update their B fields with a new u
+ for ferrite in ferritics:
+ i = ferrite[3][0]
+ j = ferrite[3][1]
+ ferrite[0] = np.array([elements[i,j][1], elements[i,j][2]]) * ur # set the m fields here which is B * ur
+
+ for i in range(x_size):
+ for j in range(y_size):
+ pos = np.array([i,j])
+ for ferrite in ferritics:
+ ox, oy = ferrite[3] # center of the magnet
+ x_pos, y_pos = pos - [ox,oy] # current x and y positions
+ mx, my = ferrite[0] # magnetization vector
+ hx, hy = 1, 1 # height of the magnet
+ elements[i,j][1] = B_induced[i,j][1] + calc_magnet_B(x_pos, y_pos, mx, my, hx, hy, ox, oy)[0]
+ elements[i,j][2] = B_induced[i,j][2] + calc_magnet_B(x_pos, y_pos, mx, my, hx, hy, ox, oy)[1]
+ ferriteON = True
+ # Now let's put every x and y value into its own array for easier plotting with quiver
+ x_samples = []
+ y_samples = []
+ x_vals = []
+ y_vals = []
+ storage = []
+ for i in range(x_size):
+ for j in range(y_size):
+ x_samples.append(i)
+ y_samples.append(j)
+ x_vals.append(elements[i,j][1])
+ y_vals.append(elements[i,j][2])
+ if i == round(x_size/2) and j >= 5:
+ print(elements[i,j][2])
+ storage.append(elements[i,j][2])
+ magnitudes = np.sqrt(np.square(x_vals) + np.square(y_vals))
+ # Plotting
+ fig, ax = plt.subplots()
+ ax.set_aspect(1.0)
+
+ norm = colors.LogNorm()
+ norm.autoscale(magnitudes)
+ cm = plt.cm.jet
+ sm = plt.cm.ScalarMappable(cmap=cm, norm=norm)
+ sm.set_array([])
+ abc = ax.quiver(x_samples, y_samples, x_vals/magnitudes, y_vals/magnitudes, color=cm(norm(magnitudes)))
+ fig.colorbar(sm)
+
+ # Add borders for Magnets and Wires, and Ferrites-?
+ for wire in wires:
+ wirePatch = patches.Rectangle(wire[3] - np.array([0.5,0.5]), wire[1], wire[2], linewidth=1, edgecolor='g', facecolor='none')
+ ax.add_patch(wirePatch)
+ for magnet in magnets:
+ magnetPatch = patches.Rectangle(magnet[3] - np.array([magnet[1]/2, magnet[2]/2]), magnet[1], magnet[2], linewidth=1, edgecolor='r', facecolor='none')
+ ax.add_patch(magnetPatch)
+ for ferrite in ferritics:
+ ferritePatch = patches.Rectangle(ferrite[3] - np.array([ferrite[1]/2, ferrite[2]/2]), ferrite[1], ferrite[2], linewidth=.5, edgecolor='b', facecolor='none')
+ ax.add_patch(ferritePatch)
+
+ # Save the figure.
+ fig.savefig(out_path + str(iteration) + ".png", dpi=150)
+ plt.close(fig)
+ #plt.show()
\ No newline at end of file
diff --git a/Python FEA Simulations/pitFinalMaterials.py b/Python FEA Simulations/pitFinalMaterials.py
new file mode 100644
index 0000000000000000000000000000000000000000..e7ee5510c4095edf2f6744a0d168ca10a550626e
--- /dev/null
+++ b/Python FEA Simulations/pitFinalMaterials.py
@@ -0,0 +1,66 @@
+import numpy as np
+
+def materials(x_size, y_size):
+ """Assign our materials of interests
+ """
+ # CREATE MAGNETS
+ magnets = [] # Now let's assign regions of permanent magnetic material
+ def magnetMaker(a, b, c, d):
+ m = np.array(a) # Bohr Magneton is 9.274e-24 J/T for context
+ mag_w = b # x dimension - must be odd
+ mag_h = c # y dimension - must be odd
+ mag_o = d # origin
+ magnets.append([m, mag_w, mag_h, mag_o,])
+ # syntax is: ferriteMaker([mx,my], mag_w, mag_h, [mag_ox, mag_oy])
+ offset = -10
+ magnetMaker([0,1],3,9,[round(x_size/2)+offset, round(y_size/2)-1])
+
+ # CREATE WIRES
+ wires = [] # Now let's assign a region of current carrying wire
+ def wireMaker(a, b, c, d):
+ Idl = np.array((0,0,a)) # current flow vector for biot-savart
+ wire_w = b
+ wire_h = c
+ wire_o = [d[0],d[1]]#[int(x_size/2-wire_w/2), int(y_size/2-wire_h/2)]
+ coords = []
+ for i in range(wire_w):
+ for j in range(wire_h):
+ wirePos = [wire_o[0] + i, wire_o[1] + j]
+ coords.append(wirePos)
+ wires.append([Idl, wire_w, wire_h, wire_o, coords])
+ # syntax is: wireMaker(Idl,wire_w,wire_h,[x_center, y_center])
+ z = 1
+ #wireMaker(z,2,1,[x_size/2 - 24, 5])
+ #wireMaker(z,2,1,[x_size/2 - 20, 5])
+ #wireMaker(z,2,1,[x_size/2 - 16, 5])
+ wireMaker(z,2,1,[x_size/2 - 12, 5])
+ wireMaker(z,2,1,[x_size/2 - 8, 5])
+ wireMaker(z,2,1,[x_size/2 - 4, 5])
+
+ #wireMaker(-z,2,1,[x_size/2 + 24, 5])
+ #wireMaker(-z,2,1,[x_size/2 + 20, 5])
+ #wireMaker(-z,2,1,[x_size/2 + 16, 5])
+ wireMaker(-z,2,1,[x_size/2 + 12, 5])
+ wireMaker(-z,2,1,[x_size/2 + 8, 5])
+ wireMaker(-z,2,1,[x_size/2 + 4, 5])
+ for wire in wires:
+ wire[3][0] -= 1
+
+ # CREATE FERITIC MATERIALS
+ ferritics = [] # Now let's assign a region of current carrying wire
+ def ferriteMaker(a, b, c, d):
+ M = np.array([0,0]) # current flow vector for biot-savart
+ ferr_w = b
+ ferr_h = c
+ ferr_o = [d[0],d[1]]#[int(x_size/2-wire_w/2), int(y_size/2-wire_h/2)]
+ for i in range(ferr_w):
+ for j in range(ferr_h):
+ ferrPos = [ferr_o[0] + i, ferr_o[1] + j]
+ ferritics.append([M, 1, 1, ferrPos])
+ # syntax is: ferriteMaker(m,1,1,[x_pos, y_pos])
+ #ferriteMaker(0,4,15,[round(x_size/2) + 0, round(y_size/2)-1 + 2 - 9])
+
+
+
+ # Now return everything
+ return magnets, wires, ferritics
\ No newline at end of file
diff --git a/Python FEA Simulations/pitFinalPhysics.py b/Python FEA Simulations/pitFinalPhysics.py
new file mode 100644
index 0000000000000000000000000000000000000000..026bda7ce11249a1fad54f22f85c5e89854a7ead
--- /dev/null
+++ b/Python FEA Simulations/pitFinalPhysics.py
@@ -0,0 +1,61 @@
+import numpy as np
+import math
+
+def current_sheet_field_2d(x_pos, y_pos, jz, y_lo, y_hi):
+ """Computes the B field (in teslas) of a sheet of current in the YZ plane.
+
+ The sheet extends from y = y_lo to y_hi in meters, and is positioned at x = 0.
+ jz gives the linear current density per meter along z (positive is out of the page)
+ """
+ scale = 1.0e-7 * jz
+ x_sqr = np.square(x_pos)
+ y_diff_hi = y_pos - y_hi
+ y_diff_lo = y_pos - y_lo
+ r_hi = np.sqrt(x_sqr + np.square(y_diff_hi))
+ r_lo = np.sqrt(x_sqr + np.square(y_diff_lo))
+
+ bx = -scale * (1 / r_hi - 1 / r_lo)
+ #print(bx)
+ by = scale * x_pos * (1 / (r_hi * (y_diff_hi + r_hi)) - 1 / (r_lo * (y_diff_lo + r_lo)))
+ if math.isnan(bx) == True or math.isnan(by) == True:
+ print(x_pos, y_pos)
+ return bx, by
+
+def calc_magnet_B(x_pos, y_pos, mx, my, hx, hy, ox, oy):
+ """Computes the B field teslas of a rectangle in the XY plane with uniform magnetization.
+
+ The rectangular region extends from x = -hx/2 to hx/2, and from y = -hy/2 to hy/2.
+ The magnetization is given by the components mx and my.
+ The answer is computed by finding the fields generated by the resulting bound current sheets
+ along the boundary of the region.
+ """
+ # right (bound current is -my)
+ bx, by = current_sheet_field_2d(x_pos - 0.5 * hx, y_pos, -my, -0.5 * hy, 0.5 * hy)
+
+ # left (bound current is my)
+ x_tmp, y_tmp = current_sheet_field_2d(x_pos + 0.5 * hx, y_pos, my, -0.5 * hy, 0.5 * hy)
+ bx += x_tmp
+ by += y_tmp
+
+ # top (bound current is mx)
+ # We rotate 90 degrees by sending y to x and x to -y.
+ x_tmp, y_tmp = current_sheet_field_2d(y_pos - 0.5 * hy, -x_pos, mx, -0.5 * hx, 0.5 * hx)
+ # Here we undo the rotation by sending x to y and y to -x.
+ bx -= y_tmp
+ by += x_tmp
+
+ # bottom (bound current -mx)
+ x_tmp, y_tmp = current_sheet_field_2d(y_pos + 0.5 * hy, -x_pos, -mx, -0.5 * hx, 0.5 * hx)
+ bx -= y_tmp
+ by += x_tmp
+ #return bx, by
+ #print(bx)
+ #input("input then enter")
+ return np.array([bx, by])
+
+def calc_biot_B(Idl, r):
+ u0 = 1e-7 # the 4*pi gets canceled later so elminiating it for speed
+ r_mag = np.sqrt(r.dot(r))
+ r_hat = r/r_mag
+ biot_B = u0/(r_mag**2)*(np.cross(Idl,r_hat))
+ return biot_B # numpy will automatically resolve the z componenet, so no need to delete it or append it to r_hat initially
diff --git a/README.md b/README.md
index 7d1548d44de7a2be9f68ca8b53d27d6fab65cb83..c70653a72df594e32338424f005268a4ab52709c 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
# Electromagnetic FEA Simulations
-We can create rudimentary 2D and 3D FEA simulations of current carrying wires and permanent magnets in free space using only two equations:
+In python we can create rudimentary 2D and 3D FEA simulations of current carrying wires and permanent magnets in free space using only two equations:
For current carrying wires we use [Biot-Savart](https://en.wikipedia.org/wiki/Biot%E2%80%93Savart_law):
@@ -30,7 +30,7 @@ Examining the results above, there are some inaccuracies and limitations worth n
## Attempt #2 - Magnets as Current Sheets
-Lots of changes to make after several very productive white-boarding sessions with Erik, who entirely came up with the solutions outlined below. First is that in order to address problem #2 above, and to properly calculate the field within a magnetized material, we need to switch from the dipole equation to something that will properly calculate a magnetic material's self-contribution. To do this, we can think of making the switch to effectively modeling a magnet with Biot-Savart, and rather than using the "free current" associated with charge flowing through a wire, we will use "bound current" associated with electron orbits. We can relate this bound current to the magnetization vector we are familiar with by taking its curl.
+Lots of changes to make after several very productive white-boarding sessions with Erik, who came up with the solution to the dipole problem outlined below. First is that in order to address problem #2 above, and to properly calculate the field within a magnetized material, we need to switch from the dipole equation to something that will properly calculate a magnetic material's self-contribution. To do this, we can think of making the switch to effectively modeling a magnet with Biot-Savart, and rather than using the "free current" associated with charge flowing through a wire, we will use "bound current" associated with electron orbits. We can relate this bound current to the magnetization vector we are familiar with by taking its curl.
@@ -84,22 +84,43 @@ For designing motors (assuming we stick with magnetostatics), the final piece re
So in the end more work to do here. We expect to see most of the permanent magnet's flux to recirculate through the steel, which would also reduce the flux outside of both materials (aside from the path into and out of the steel).
-## More things to try / consider
+## More Things
+
+### Running the Code
+
+All of the code is stored in this repository. From the command line you can call pitFinalFEA.py, which plots the B vector field contributed to an elements by a magnetic material at some distance away.
+
+The relevant physics functions are imported from pitFinalPhysics.py, where calc_magnet is used for both permanent magnets and magnetized ferritic materials.
+```python
+calc_biot_B(Idl, r) # returns Bx and By for a a current*length Idl at some distance and orientation r
+calc_magnet_B(x_pos, y_pos, mx, my, hx, hy, ox, oy) # returns Bx and By for some magnetized material with magnetization vector [mx,my]
+```
+
+To create materials we import pitFinalMaterials.py where we can append lines of code with the origin, height, width, and magnetic properties of an element or region of elements.
+
+```python
+wireMaker(Idl,wire_w,wire_h,[x_center, y_center]) # Idl sign will dictate right hand rule
+magnetMaker([mx,my], mag_w, mag_h, [mag_ox, mag_oy])
+ferriteMaker(m,ferr_w1,ferr_h,[x_pos, y_pos]) # note that this will then create w*h 1x1 elements rather than monolithic block
+```
### 2D Simulations to 3D
-In the validation section we talked about the limitations of looking at the simulation results in an absolute sense. If we wanted to, we should be able to extrapolate some symmetric 3D cases from these 2D results.
+In the validation section we talked about the limitations of looking at the simulation results in an absolute sense. If we wanted to, we should be able to extrapolate some symmetric 3D cases from these 2D results. Comsol has a very nice interface for this where you can actually revolve a 2D simulation to integrate a 3rd dimension without much computational overhead. I think
-For example, Biot-Savart tells us the contribution of the infinitesimal wire length dl, but what if the wire extends 5cm into the screen, and an adjacent magnet is only 2cm thick? In that case the sum of contributions from each element dl all contribute to the field strength at point A, but with increasing effective radius, and decreasing angle theta, which results in the cross product dl x r_hat also decreasing in magnitude. It turns out that if we integrate the contribution of an infinite wire on a single point A, we scale Biot-Savart by a factor of 2r. From my understanding, if we assume all magnetic materials scale infinitely into the page, this would result in linear scaling and could still provide useful relative approximations.
+Where I think these simulations can get you into trouble is if you are modeling materials with different length scales. In that case the sum of contributions from each element dl all contribute to the field strength at point A, but with increasing effective radius, and decreasing angle theta, which results in the cross product dl x r_hat also decreasing in magnitude. From one of the psets I was under the impression that if we integrate the contribution of an infinite wire on a single point A, we scale Biot-Savart by a factor of 2r. From my understanding, if we assume all magnetic materials scale infinitely into the page, this would result in linear scaling and could still provide useful relative approximations.
-In the case of a coil, center axis of the coil is at a fixed distance from the circumference of wire (assuming it isn't spiraling), and therefore we can scale Biot-Savart by 2*pi*r, but as we move towards the perimeter of the coil radially, the 3 dimensionality of the coil becomes more complex and nonlinear. I believe that unlike the infinite wire approximation, this would not result in linear scaling and perhaps would result in error unless accounted for?
+In the case of a revolved coil for example, the center axis of the coil is at a fixed distance from the circumference of wire (assuming it isn't spiraling), and therefore we can scale Biot-Savart by 2*pi*r, but as we move towards the perimeter of the coil radially, the 3 dimensionality of the coil becomes more complex and nonlinear. I believe that unlike the infinite wire approximation, this would not result in linear scaling and perhaps would result in error unless accounted for?
### Calculating Forces
We have a path to calculating forces via the [Lorentz Force Law](https://en.wikipedia.org/wiki/Lorentz_force) but I didn't get around to this.
-### Differences between this and a commercial FEA package
+### Differences between this and a commercial FEA software
+* Software this project is trying to emulate:
+ * [FEMM - Finite Element Method Magnetics](https://www.femm.info/wiki/HomePage)
+ * [Comsol](https://www.comsol.com/support/learning-center/article/Modeling-Electromagnetic-Coils-8251)
* Erik pointed out that rather than go stright to B fields, most commercial solvers will likely calculate H and then... This explains why commercial solvers are so concerned with enforcing an infinite air boundary around their models...
* Real FEA packages do a lot of front-end work to mesh their models (triangular meshes rather than just enforcing a uniform grid). This allows for users to get higher resolution results in regions of interest, and avoid enforcing an element size dictated by the simulations smallest regions of interest. With that said, the meshing problem can be looked at as completely seperate from the physics, and I'd imagine at Comsol there is a meshing department and then various physics departments that don't interact too much with one another...
* Run on the GPU
diff --git a/images/diySim1.png b/images/diySim1.png
index 13f07287c8bf0b36ab9965a47edae1620bc07205..f468602f930b0ec0dfd017d0bfaca383bdc8e951 100644
Binary files a/images/diySim1.png and b/images/diySim1.png differ