Commit 2124f847 authored by Siamak Pourkeivannour's avatar Siamak Pourkeivannour
Browse files

Still AC not working

parent fd94fb85
......@@ -5,3 +5,5 @@
Ccore GMSH/
From Andrea/
2D_MH_lin.py
Ask Joost/
from Mitrofan/
import nutils, numpy
import matplotlib.pyplot as plt
import matplotlib.collections
from nutils import cli, mesh, export
import time
__matrix__ = nutils.matrix.MKL() or nutils.matrix.Scipy() or nutils.matrix.Numpy()
__matrix__.__enter__()
def main():
pi = numpy.pi
btype = 'std'
degree = 2
m = 2 # number of winding layers
L = 50e-3
nelems = 40
nelems_x = tuple([20])
nelems_y = tuple([20])
lengths_x = tuple([50e-3])
lengths_y = tuple([50e-3])
subnames_x = 'CondX'
subnames_y = 'CondY'
subverts_x = numpy.cumsum([0, *lengths_x])
verts_x = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_x[:-1], subverts_x[1:], nelems_x)] + [
subverts_x[-1:]])
subverts_y = numpy.cumsum([0, *lengths_y])
verts_y = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_y[:-1], subverts_y[1:], nelems_y)] + [
subverts_y[-1:]])
domain, geom = nutils.mesh.rectilinear([verts_x, verts_y])
#
# yverts = numpy.cumsum([0, *nelems_y])
# domain = domain.withgroups(vgroups={k: domain[:, l:r] for k, l, r in zip(subnames_y, yverts[:-1], yverts[1:])})
# xverts = numpy.cumsum([0, *nelems_x])
# domain = domain.withgroups(
# vgroups={k: domain['CondY'][l:r] for k, l, r in zip(subnames_x, xverts[:-1], xverts[1:])})
ns = nutils.function.Namespace()
ns.x = geom
# with export.mplfigure('geometry.png') as fig:
# # Add geometry
# ax = fig.add_subplot(1, 1, 1)
# #ax.set_aspect('equal')
# smpl_gauss = domain.sample('gauss',degree*2)
# smpl_bezier = domain.sample('bezier',15)
# gauss = smpl_gauss.eval(geom)
# bezier = smpl_bezier.eval(geom)
# ax.add_collection(matplotlib.collections.LineCollection(bezier[smpl_bezier.hull], colors='k', linewidths=0.25, alpha=.25))
# # Add Gauss points
# #for i, q in enumerate(gauss):
# # ax.text(q[0] - 1e-3, q[1] + 0.25e-3, i, fontsize=4)
# # ax.plot(q[0], q[1], 'k.', markersize=1)
#
# ax.set_xlim([0, 50e-3])
# ax.set_ylim([0, 50e-3])
# Define helper function for creating values that are constant per subdomain.
# subdomconst = lambda *, default=0, **values: sum(v * domain.indicator(k) for k, v in values.items()) + default
start_time = time.time()
freq = 0
µ_0 = 4 * pi * 1e-7
σ_C = 5.998 * 1e7 # Copper conductivity
omega = 2 * pi * freq
I = 1 #coil current
Cond_Area = L * L
J_dc = I/Cond_Area
ns.µ0 = µ_0
ns.ϑ0 = 1 / µ_0
ns.σC = σ_C
ns.omega = omega
ns.J = J_dc
ns.mu = µ_0
ns.nu = 1 / ns.µ0
ns.sigma = ns.σC
ns.basis = domain.basis(btype, degree=degree).vector(2) #
ns.D = numpy.array([[0, 1], [-1, 0]])
ns.Az_i = 'basis_ni ?lhs_n'
ns.J_k = '<J, 0>_k'
ns.B_jk = 'D_ij Az_k,i'
# ns.curlubasis_ni = 'eps_ij ucbasis_ni,j'
# # Form the basis
# ns.ucbasis, ns.usbasis = nutils.function.chain([
# domain.basis('spline', degree=degree, knotvalues=[verts_x, verts_y]),
# domain.basis('spline', degree=degree, knotvalues=[verts_x, verts_y])])
#
# ns.eps = numpy.array([[0, 1], [-1, 0]]) # 2D
# ns.curlucbasis_ni = 'eps_ij ucbasis_n,j' # <<< ?
# ns.curlusbasis_ni = 'eps_ij usbasis_n,j'
#
# ns.uc = 'ucbasis_n ?lhs_n'
# ns.us = 'usbasis_n ?lhs_n'
#
# ns.curluc_i = 'eps_ij uc_,j'
# ns.curlus_i = 'eps_ij us_,j'
#
# res = domain.integral('(nu curlucbasis_ni curluc_i + nu curlusbasis_ni curlus_i) d:x' @ ns, degree=degree * 2)
#
# res -= domain.integral('(sigma omega ucbasis_n us - sigma omega usbasis_n uc) d:x' @ ns, degree=degree * 2)
#
# res -= domain.integral('(ucbasis_n J ) d:x' @ ns, degree=degree * 2) # careful with parenthesis
# sqr = domain.boundary['right'].integral('( curluc_1^2 + curlus_1^2 ) d:x' @ ns, degree=degree * 2)
# sqr += domain.boundary['left'].integral('( curluc_1^2 + curlus_1^2 ) d:x' @ ns, degree=degree * 2)
# sqr += domain.boundary['top'].integral('( curluc_0^2 + curlus_0^2 ) d:x' @ ns, degree=degree * 2)
# sqr += domain.boundary['bottom'].integral('( curluc_0^2 + curlus_0^2 ) d:x' @ ns, degree=degree * 2)
sqr = domain.boundary['right'].integral('( B_00^2 + B_01^2 ) d:x' @ ns, degree=degree * 2)
sqr += domain.boundary['left'].integral('( B_00^2 + B_01^2) d:x' @ ns, degree=degree * 2)
sqr += domain.boundary['top'].integral('( B_10^2 + B_11^2 ) d:x' @ ns, degree=degree * 2)
sqr += domain.boundary['bottom'].integral('( B_10^2 + B_11^2 ) d:x' @ ns, degree=degree * 2)
sqr += domain.integral(' B_00 d:x ' @ ns, degree=degree * 2)
sqr += domain.integral(' B_01 d:x ' @ ns, degree=degree * 2)
sqr += domain.integral(' B_10 d:x ' @ ns, degree=degree * 2)
sqr += domain.integral(' B_11 d:x ' @ ns, degree=degree * 2)
# sqr += domain.integral(' Az_1 d:x ' @ ns, degree=degree * 2)
cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)
res = domain.integral('( basis_nk,i Az_k,i + mu omega sigma basis_nk D_ik Az_i ) d:x' @ ns, degree=degree * 2)
res -= domain.integral('( mu basis_nk J_k ) d:x' @ ns, degree=degree * 2) # careful with parenthesis
# Solve the problem
lhs = nutils.solver.solve_linear('lhs', res, constrain=cons)
#ns.normAz = ' sqrt( Az_0^2 + Az_1^2 )'
#ns.Je = ' sigma omega Az_0 '
bezierC = domain.sample('bezier', 10)
print("--- %s seconds ---" % (time.time() - start_time))
xC, J0, nu, mu, Az_r, Az_i , Bx_r , Bx_i , By_r , By_i = bezierC.eval(['x_i', 'J_0', 'nu', 'mu', 'Az_0', 'Az_1', 'B_00', 'B_01', 'B_10', 'B_11'] @ ns, lhs=lhs)
#Plot
plt.figure(1)
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, Az_r, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Real Component of Magnetic Vector Potential [Wb/m]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
plt.figure(2)
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, Az_i, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Imaginary Component of Magnetic Vector Potential [Wb/m]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
plt.figure(3)
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, Bx_r, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Real Component of Bx [T]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
plt.figure(4)
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, Bx_i, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Imaginary Component of Bx [T]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
plt.figure(5)
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, By_r, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Real Component of By [T]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
plt.figure(6)
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, By_i, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Imaginary Component of By [T]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
plt.show()
#Define a 2d list or numpy matrix with coordinates
# x_end = 50e-3
# y_end = 50e-3
#
# x_coordinates = numpy.linspace(0,x_end,100)
# y_coordinates = numpy.linspace(0,y_end,100)
#
# coordinates = numpy.array(numpy.meshgrid(x_coordinates, y_coordinates)).T.reshape(-1, 2)
# # Define a sampling object
#
# located_points = domain.locate(geom, coordinates, tol=1e-10)
#
# # Evaluate the sampling object for a desired quantity
# Az = numpy.reshape(located_points.eval('uc' @ ns, lhs=lhs),(-1,100)).T
# Az_comsol = numpy.loadtxt('real_Az_10kHz.txt')
#
# #error = (Az-Az_comsol)/Az_
# error = 100 * numpy.divide((Az-Az_comsol), Az_comsol, out=numpy.zeros_like((Az-Az_comsol)), where=Az_comsol != 0)
# print(error)
#
# error_max = numpy.max(error)
# print(error_max)
# error_min = numpy.min(error)
# print(error_min)
# fig, ax = plt.subplots()
# c = ax.pcolormesh(x_coordinates, y_coordinates, error, shading='auto', cmap='RdBu',vmax=-95 , vmin=-105)
# ax.set_title('Deviation of Real Component of Az [Nutils - COMSOL] @10kHz')
# # set the limits of the plot to the limits of the data
# ax.axis([x_coordinates.min(), x_coordinates.max(), y_coordinates.min(), y_coordinates.max()])
# fig.colorbar(c, ax=ax)
# plt.show()
if __name__ == '__main__':
nutils.cli.run(main)
import nutils, numpy
import matplotlib.pyplot as plt
import matplotlib.collections
from nutils import cli, mesh, export
import time
__matrix__ = nutils.matrix.MKL() or nutils.matrix.Scipy() or nutils.matrix.Numpy()
__matrix__.__enter__()
def main():
pi = numpy.pi
btype = 'std'
degree = 2
L = 50e-3
nelems = 10
nelems_x = nelems,nelems,nelems
nelems_y = nelems,nelems,nelems
lengths_x = L, L, L
lengths_y = L, L, L
subnames_x = 'Air', 'CondX', 'Air'
subnames_y = 'Air', 'CondY', 'Air'
subverts_x = numpy.cumsum([0, *lengths_x])
verts_x = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_x[:-1], subverts_x[1:], nelems_x)] + [
subverts_x[-1:]])
subverts_y = numpy.cumsum([0, *lengths_y])
verts_y = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_y[:-1], subverts_y[1:], nelems_y)] + [
subverts_y[-1:]])
domain, geom = nutils.mesh.rectilinear([verts_x, verts_y])
yverts = numpy.cumsum([0, *nelems_y])
domain = domain.withgroups(vgroups={k: domain[:, l:r] for k, l, r in zip(subnames_y, yverts[:-1], yverts[1:])})
xverts = numpy.cumsum([0, *nelems_x])
domain = domain.withgroups(
vgroups={k: domain['CondY'][l:r] for k, l, r in zip(subnames_x, xverts[:-1], xverts[1:])})
ns = nutils.function.Namespace()
ns.x = geom
# with export.mplfigure('geometry.png') as fig:
# # Add geometry
# ax = fig.add_subplot(1, 1, 1)
# #ax.set_aspect('equal')
# smpl_gauss = domain.sample('gauss',degree*2)
# smpl_bezier = domain.sample('bezier',15)
# gauss = smpl_gauss.eval(geom)
# bezier = smpl_bezier.eval(geom)
# ax.add_collection(matplotlib.collections.LineCollection(bezier[smpl_bezier.hull], colors='k', linewidths=0.25, alpha=.25))
# # Add Gauss points
# #for i, q in enumerate(gauss):
# # ax.text(q[0] - 1e-3, q[1] + 0.25e-3, i, fontsize=4)
# # ax.plot(q[0], q[1], 'k.', markersize=1)
#
# ax.set_xlim([0, 175e-3])
# ax.set_ylim([0, 150e-3])
# Define helper function for creating values that are constant per subdomain.
subdomconst = lambda *, default=0, **values: sum(v * domain.indicator(k) for k, v in values.items()) + default
start_time = time.time()
freq = 100
µ_0 = 1.257*1e-6
σ_C = 5.998* 1e7 # Copper conductivity
ns.SigmaCond = σ_C
omega = 2 * pi * freq
I = 1 #coil current
Cond_Area = L * L
J_dc = I/Cond_Area
ns.µ0 = µ_0
ns.ϑ0 = 1 / µ_0
ns.σC = σ_C
ns.omega = omega
ns.J = J_dc
ns.mu = subdomconst(default= (1 * ns.µ0))
ns.nu = subdomconst(default= (1 / ns.µ0))
ns.sigma = subdomconst(CondX=ns.σC, default=0)
ns.Jdc = subdomconst( CondX=ns.J, default=0)
ns.delta = ns.mu * ns.sigma * ns.omega
spline_km_y = numpy.ones([sum(nelems_y) + 1], int)
spline_km_y[numpy.cumsum([0, *nelems_y])] = degree
ns.ucbasis, ns.usbasis = nutils.function.chain([
domain.basis(btype, degree=degree),
domain.basis(btype, degree=degree)])
# Setup magnetics problem
ns.eps = numpy.array([[0, 1], [-1, 0]]) # 2D
ns.curlucbasis_ni = 'eps_ij ucbasis_n,j' # <<< ?
ns.curlusbasis_ni = 'eps_ij usbasis_n,j'
ns.uc = 'ucbasis_n ?lhs_n'
ns.us = 'usbasis_n ?lhs_n'
ns.curluc_i = 'eps_ij uc_,j'
ns.curlus_i = 'eps_ij us_,j'
sqr = domain.boundary['bottom'].integral('(uc^2 + us^2) d:x' @ ns, degree=degree)
sqr += domain.boundary['top'].integral('(uc^2 + us^2) d:x' @ ns, degree=degree)
cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)
res = domain.integral('( curlucbasis_ni curluc_i + curlusbasis_ni curlus_i ) d:x' @ ns, degree=degree)
res -= domain.integral('( delta ucbasis_n us - delta usbasis_n uc ) d:x' @ ns, degree=degree)
res += domain.integral('( mu ucbasis_n Jdc ) d:x' @ ns, degree=degree) # careful with parenthesis
# Solve the problem
lhs = nutils.solver.solve_linear('lhs', res, constrain=cons)
ns.A = ' sqrt( uc uc + us us )'
bezierC = domain.sample('bezier', 10)
xC, Az_c, Az_s, A , J_0 = bezierC.eval(['x_i','uc', 'us', 'A', 'Jdc'] @ ns, lhs=lhs)
print("--- %s seconds ---" % (time.time() - start_time))
# Plot
plt.figure(2)
plt.tripcolor(xC[:, 0] * 1000, xC[:, 1] * 1000, bezierC.tri, A, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Magnetic Vector Potential [Wb/m]')
plt.xlabel('x [mm]')
plt.ylabel('y [mm]')
matplotlib.rcParams["figure.dpi"] = 0.1
plt.show()
# # Define a 2d list or numpy matrix with coordinates
# x_end = 5*L
# y_end = 3*L
# x_coordinates = numpy.linspace(0,x_end,250)
# y_coordinates = numpy.linspace(0,y_end,150)
# coordinates = numpy.array(numpy.meshgrid(x_coordinates, y_coordinates)).T.reshape(-1, 2)
# # Define a sampling object
# located_points = domain.locate(geom, coordinates, tol=1e-12)
# # Evaluate the sampling object for a desired quantity
# Az = numpy.reshape(located_points.eval('u' @ ns, lhs=lhs), (-1,150)).T
# print(len(Az))
# Az_comsol = numpy.loadtxt('Az_COMSOL.txt')
# print(len(Az_comsol))
#
# error = 100*numpy.true_divide((Az - Az_comsol),Az_comsol)
# error_max = numpy.max(error)
# print(error_max)
# error_min = numpy.min(error)
# print(error_min)
#
# fig, ax = plt.subplots()
# c = ax.pcolormesh(x_coordinates, y_coordinates, error, cmap='RdBu', vmin=error_min, vmax=error_max)
# ax.set_title('Error (%)')
# # set the limits of the plot to the limits of the data
# ax.axis([x_coordinates.min(), x_coordinates.max(), y_coordinates.min(), y_coordinates.max()])
# fig.colorbar(c, ax=ax)
# plt.show()
if __name__ == '__main__':
nutils.cli.run(main)
import nutils, numpy
import matplotlib.pyplot as plt
import matplotlib.collections
from nutils import cli, mesh, export
import time
__matrix__ = nutils.matrix.MKL() or nutils.matrix.Scipy() or nutils.matrix.Numpy()
__matrix__.__enter__()
def main():
pi = numpy.pi
btype = 'std'
degree = 2
L = 50e-3
nelems = 30
nelems_x = nelems, nelems, nelems
nelems_y = nelems, nelems, nelems
lengths_x = L, L, L
lengths_y = L, L, L
subnames_x = 'Air', 'CondX', 'Air'
subnames_y = 'Air', 'CondY', 'Air'
subverts_x = numpy.cumsum([0, *lengths_x])
verts_x = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_x[:-1], subverts_x[1:], nelems_x)] + [
subverts_x[-1:]])
subverts_y = numpy.cumsum([0, *lengths_y])
verts_y = numpy.concatenate(
[numpy.linspace(l, r, n + 1)[:-1] for l, r, n in zip(subverts_y[:-1], subverts_y[1:], nelems_y)] + [
subverts_y[-1:]])
domain, geom = nutils.mesh.rectilinear([verts_x, verts_y])
yverts = numpy.cumsum([0, *nelems_y])
domain = domain.withgroups(vgroups={k: domain[:, l:r] for k, l, r in zip(subnames_y, yverts[:-1], yverts[1:])})
xverts = numpy.cumsum([0, *nelems_x])
domain = domain.withgroups(
vgroups={k: domain['CondY'][l:r] for k, l, r in zip(subnames_x, xverts[:-1], xverts[1:])})
ns = nutils.function.Namespace()
ns.x = geom
# with export.mplfigure('geometry.png') as fig:
# # Add geometry
# ax = fig.add_subplot(1, 1, 1)
# #ax.set_aspect('equal')
# smpl_gauss = domain.sample('gauss',degree*2)
# smpl_bezier = domain.sample('bezier',15)
# gauss = smpl_gauss.eval(geom)
# bezier = smpl_bezier.eval(geom)
# ax.add_collection(matplotlib.collections.LineCollection(bezier[smpl_bezier.hull], colors='k', linewidths=0.25, alpha=.25))
# # Add Gauss points
# #for i, q in enumerate(gauss):
# # ax.text(q[0] - 1e-3, q[1] + 0.25e-3, i, fontsize=4)
# # ax.plot(q[0], q[1], 'k.', markersize=1)
#
# ax.set_xlim([0, 175e-3])
# ax.set_ylim([0, 150e-3])
# Define helper function for creating values that are constant per subdomain.
subdomconst = lambda *, default=0, **values: sum(v * domain.indicator(k) for k, v in values.items()) + default
freq = 1
ns.mu0 = pi * 4e-7
ns.sigmaC = 5.988e7 # Copper conductivity
ns.omega = 2 * pi * freq
I = 1 # coil current
Cond_Area = domain['CondX'].integral('1 d:x' @ ns, degree=degree * 2).eval()
#nutils.log.user('Conductor area: {:.4e}'.format(Cond_Area))
ns.Jdc = I / Cond_Area
ns.nu = 1 / ns.mu0
ns.sigma = subdomconst(CondX= ns.sigmaC, default=0)
ns.Js = subdomconst(CondX= ns.Jdc, default=0)
ns.delta = ns.sigma * ns.omega
ns.basis = domain.basis(btype, degree=degree).vector(2)
ns.D = numpy.array([[0, 1], [-1, 0]])
ns.Az_i = 'basis_ni ?lhs_n'
ns.J_k = '<Js, 0>_k'
ns.B_jk = 'D_ij Az_k,i'
ns.normA = ' sqrt( Az_0^2 + Az_1^2 )'
ns.Je_i = ' delta D_ij Az_j'
sqr = domain.boundary['right'].integral('( Az_0^2 + Az_1^2 ) d:x' @ ns, degree=degree * 2)
sqr += domain.boundary['left'].integral('( Az_0^2 + Az_1^2) d:x' @ ns, degree=degree * 2)
sqr += domain.boundary['top'].integral('( Az_0^2 + Az_1^2 ) d:x' @ ns, degree=degree * 2)
sqr += domain.boundary['bottom'].integral('( Az_0^2 + Az_1^2 ) d:x' @ ns, degree=degree * 2)
sqr += domain['CondX'].integral(' Je_0 d:x ' @ ns, degree=degree * 2)
sqr += domain['CondX'].integral(' Je_1 d:x ' @ ns, degree=degree * 2)
cons = nutils.solver.optimize('lhs', sqr, droptol=1e-15)
res = domain.integral('( nu basis_nk,i Az_k,i + delta basis_nk D_ik Az_i ) d:x' @ ns, degree=degree * 2)
res -= domain.integral('( basis_nk J_k ) d:x' @ ns, degree=degree * 2) # careful with parenthesis
# Solve the problem
lhs = nutils.solver.solve_linear('lhs', res, constrain=cons)
bezierC = domain.sample('bezier', 8)
# xC, Je20, Je21, divA, divJe, divJe20 = bezierC.eval(['x_i', 'Je_0', 'Je2_1', 'A_,0 + A_,1', 'Je_,0 + Je_,1','Je2_0,0 + Je2_0,1'] @ ns, lhs=lhs)
xC, Az_r, Az_i, Je_r, Je_i , Bx_r , Bx_i , By_r , By_i = bezierC.eval(['x_i', 'Az_0' , 'Az_1', 'Je_0', 'Je_1', 'B_00', 'B_01', 'B_10', 'B_11'] @ ns, lhs=lhs)
#IntdivJe = domain['CondX'].integral('(Je_,0 + Je_,1) d:x' @ ns, degree=degree * 2).eval(lhs=lhs)
#nutils.log.user('IntdivJe equal to: {:.4e} Watts'.format(IntdivJe))
plt.figure(1)
plt.tripcolor(xC[:, 0], xC[:, 1], bezierC.tri, Je_r, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Je real [A/m^2]')
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.figure(2)
plt.tripcolor(xC[:, 0], xC[:, 1], bezierC.tri, Je_i, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Je imag [A/m^2]')
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.figure(3)
plt.tripcolor(xC[:, 0], xC[:, 1], bezierC.tri, Az_r, shading='gouraud', cmap='cividis')
plt.colorbar(orientation="vertical")
plt.title('Magnetic Vector Potential real [Wb/m]')
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.figure(4)
plt.tripcolor(xC[:, 0], xC[:, 1], bezierC.tri, Az_i, shading='gouraud', cmap='cividis')