Hi, guys, what’s up?
I’m an undergraduate student and FEniCS enthusiast. I’m still learning how to use everything it offers. At the moment, I’m trying to refactor the code below for a version that uses dolfinx.
The first graphics are not a problem. I start to have doubts when creating function spaces.
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d, splev, splrep
from mshr import *
mi0 = 4*np.pi*1e-7
Hc = 1000
# mesh parameters
nmesh = 30 # mesh density
tol = 0 # mesh tolerance
# geometry parameters
# air box
r_air = 80e-3
z_air = 80e-3
air = Rectangle(Point(tol,-z_air), Point(r_air,z_air))
# Magnet
r_mag = 30e-3
z_mag = 30e-3
core = Rectangle(Point(tol,-z_mag), Point(r_mag,z_mag))
# domain
domain = air + core
# Set subdomain for wire
domain.set_subdomain(1,core)
# Create mesh
mesh = generate_mesh(domain, nmesh)
# Plot mesh
fig = plt.figure(figsize=(10,10))
plot(mesh,title='Fig. 1 - Problem geometry');
## Interpolation - Define non-linear permeability
# Steel
Bi_steel = np.array([0, 0.221950, 0.245515, 0.344303, 0.375573, 0.454417, 0.627981,
0.670134, 0.861453, 1.075180, 1.241074, 1.423388, 1.656238, 1.686626,
1.813505, 1.964422, 1.979083, 2.012433, 2.021337, 2.033503, 2.050973,
2.052071, 2.191983, 2.197328, 2.240825, 2.309729, 2.327795, 2.435784])
Hi_steel = np.array([tol, 225.366667, 237.316667, 291.793333, 310.450000, 358.730000,
483.890000, 520.136667, 723.673333, 1071.333333, 1570.566667,
2775.500000, 6290.533333, 7049.866667, 12338.666667, 26304.666667,
28581.000000, 36287.000000, 39022.333333, 43292.333333, 50590.000000,
51118.333333, 134313.333333, 138566.666667, 168803.333333, 223476.666667,
237853.333333, 321480.000000])
mi_steel = np.divide(Bi_steel, Hi_steel)
mi_steel[0] = mi_steel[1]/2
#mi_steel[0] = mi0
plt.plot(Bi_steel, mi_steel/mi0,'-.');
plt.xlabel('Bi');
plt.ylabel('$\mu$');
f = interp1d(Bi_steel, mi_steel)
f2 = interp1d(Bi_steel, mi_steel, kind='quadratic')
spl = splrep(Bi_steel, mi_steel)
xnew = np.linspace(0, Bi_steel[-1], num = 101)
sp2 = splev(xnew, spl)
fig = plt.figure(figsize=(10,10))
plt.plot(Bi_steel, mi_steel, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
plt.plot(xnew, sp2, '-.')
plt.legend(['data', 'linear', 'quadratic', 'spline'], loc='best')
plt.show()
# Define function space
V = FunctionSpace(mesh, 'P', 2)
# Define boundary condition
bc = DirichletBC(V, Constant(0), 'on_boundary')
r = Expression('x[0]', degree=1)
W = VectorFunctionSpace(mesh, 'P', 2)
phi = Function(V)
def mi(phi):
B = project(as_vector((-(1/r)*phi.dx(1), (1/r)*phi.dx(0))), W)
norm_B = norm(B)
return f2(norm_B)
# Define subdomain markers and integration measure
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)
# Define current sheet
boundary_markers = MeshFunction('size_t', mesh, 1)
class IMA(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return near(x[0], r_mag, tol) and x[1] <= z_mag and x[1] >= -z_mag
sheet = IMA()
sheet.mark(boundary_markers, 1)
dS = Measure('dS', domain=mesh, subdomain_data=boundary_markers)
# Define magnetic permeability
class Permeability(UserExpression):
def __init__(self, markers, **kwargs):
super().__init__(**kwargs)
self.markers = markers
self.mi_0 = mi0 #coil and air
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = mi(phi)
else:
values[0] = self.mi_0
mu = Permeability(markers, degree=1)
fig = plt.figure(figsize=(10,10));
plot(markers,title='Fig. 2 - Subdomínios do problema');
# Define variational problem
# Define current densities
f = 0 #
# Define variational problem
phi = Function(V) # rho * A_theta
v = TestFunction(V)
# Poisson equation in axisymmetric cylindrical coordinates
r = Expression('x[0]', degree=1)
F = (1/mu)*(1/r)*(Dx(v,0)*Dx(phi,0) + Dx(v,1)*Dx(phi,1))*dx - Hc*v('+')*dS(1)
# Compute solution
solve(F==0, phi, bc)
fig = plt.figure(figsize=(10,10))
auxfig = plot(phi);
plt.colorbar(auxfig)
# Compute magnetic field (B = curl A)
W = VectorFunctionSpace(mesh, 'P', 1)
B = project(as_vector((-(1/r)*phi.dx(1),(1/r)*phi.dx(0))), W);
norm_B = sqrt(inner(B,B));
fig = plt.figure(figsize=(10,10))
auxfig = plot(norm_B);
plt.colorbar(auxfig)
fig = plt.figure(figsize=(10,10))
plot(B)
# Plot norm_B
z = np.linspace(-z_air , z_air, 101)
v_line = [(r_mag/2, z_) for z_ in z]
norm_B_v_line = np.array([norm_B(point) for point in v_line])
fig = plt.figure(figsize=(8,5))
plt.plot(z, norm_B_v_line)
plt.grid(True)
plt.xlabel('$z$')
plt.ylabel('$|B| (T)$')
plt.plot(z, f2(norm_B_v_line))
I’ve already done the code below, but the graphs don’t make sense. I’d like to know what you think. Any help would be very welcome.
# Basic Libraries
import numpy as np # For vector and matrix manipulation
import gmsh # Interface with the software GMSH
from scipy.interpolate import interp1d, splev, splrep
from petsc4py import PETSc
# Visualization
import pyvista
from mpi4py import MPI
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import plotly.express as px
# FEniCS
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dx, dot, grad, as_vector, sqrt, inner
from dolfinx import default_scalar_type
from dolfinx.io import gmshio, XDMFFile, VTKFile
from dolfinx.fem import Constant, Expression, Function, functionspace, dirichletbc, locate_dofs_geometrical, locate_dofs_topological
from dolfinx.mesh import locate_entities_boundary
from dolfinx.fem.petsc import LinearProblem
from dolfinx.plot import vtk_mesh
rank = MPI.COMM_WORLD.rank
gmsh.initialize()
mi0 = 4*np.pi*1e-7
Hc = 1000
inch = 25.4e-3
# Domínio com ar
r_air = 80e-3
z_air = 80e-3
# Ímã
r_mag = 30e-3
z_mag = 30e-3
# Dimensões do ímã e model rank
gdim = 2 # Geometric dimension of the mesh
model_rank = 0
mesh_comm = MPI.COMM_WORLD
if mesh_comm.rank == model_rank:
# Tags
inner_tag = 1 # interior do ímã
outer_tag = 2 # domínio com ar
# Cria o retângulo de ar
air_rectangle = gmsh.model.occ.addRectangle(0, -z_air, 0, r_air, 2*z_air, tag=outer_tag)
gmsh.model.occ.synchronize()
# Cria o retângulo do ímã
magnet_rectangle = gmsh.model.occ.addRectangle(0, -z_mag, 0, r_mag, 2*z_mag,tag = inner_tag)
gmsh.model.occ.synchronize()
whole_domain = gmsh.model.occ.fragment([(2, air_rectangle)], [(2, magnet_rectangle)])
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(gdim, [inner_tag], tag = 1, name = "inner") # Elementos internos
gmsh.model.addPhysicalGroup(gdim, [ outer_tag], tag = 2, name = "outer") # Elementos externos
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", inch/8)
gmsh.model.mesh.generate(gdim)
gmsh.write("ima.msh")
# Converte a malha criada pelo gmsh para um formato que o FEniCS compreende
mesh,cell_tags,facet_tags = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
# Finaliza o GMSH
gmsh.finalize()
# Steel
Bi_steel = np.array([0, 0.221950, 0.245515, 0.344303, 0.375573, 0.454417, 0.627981,
0.670134, 0.861453, 1.075180, 1.241074, 1.423388, 1.656238, 1.686626,
1.813505, 1.964422, 1.979083, 2.012433, 2.021337, 2.033503, 2.050973,
2.052071, 2.191983, 2.197328, 2.240825, 2.309729, 2.327795, 2.435784])
Hi_steel = np.array([0, 225.366667, 237.316667, 291.793333, 310.450000, 358.730000,
483.890000, 520.136667, 723.673333, 1071.333333, 1570.566667,
2775.500000, 6290.533333, 7049.866667, 12338.666667, 26304.666667,
28581.000000, 36287.000000, 39022.333333, 43292.333333, 50590.000000,
51118.333333, 134313.333333, 138566.666667, 168803.333333, 223476.666667,
237853.333333, 321480.000000])
mi_steel = np.divide(Bi_steel, Hi_steel)
mi_steel[0] = mi_steel[1]/2
#mi_steel[0] = mi0
# Plot 1
plt.plot(Bi_steel, mi_steel/mi0,'-.');
plt.xlabel('Bi');
plt.ylabel('$\mu$');
# Plot 2
fig = go.Figure()
x = np.arange(56)
fig.add_trace(go.Scatter(x = Bi_steel, y = mi_steel/mi0, name = 'permeabilidade', line = dict(dash = 'dash', color = 'royalblue', width = 4)))
fig.update_layout(title = 'mi x Bi', xaxis_title='Bi', yaxis_title='Permeabilidade Relativa (mi)', width=800, height=600)
fig.show()
f = interp1d(Bi_steel, mi_steel)
f2 = interp1d(Bi_steel, mi_steel, kind='quadratic')
spl = splrep(Bi_steel, mi_steel)
xnew = np.linspace(0, Bi_steel[-1], num = 101)
sp2 = splev(xnew, spl)
# Plot 1
fig = plt.figure(figsize=(10,10))
plt.plot(Bi_steel, mi_steel, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
plt.plot(xnew, sp2, '-.')
plt.legend(['data', 'linear', 'quadratic', 'spline'], loc='best')
plt.show()
# Plot 2
# Criação do gráfico
fig = go.Figure()
# Adicionando os dados originais
fig.add_trace(go.Scatter(x=Bi_steel, y=mi_steel, mode='markers', name='Dados originais', marker=dict(color='black', size=8)))
# Adicionando interpolação linear
fig.add_trace(go.Scatter(x=xnew, y=f(xnew), mode='lines', name='Linear', line=dict(color='blue', width=2)))
# Adicionando interpolação quadrática
fig.add_trace(go.Scatter(x=xnew, y=f2(xnew), mode='lines', name='Quadrática', line=dict(dash='dash', color='green', width=2)))
# Adicionando spline
fig.add_trace(go.Scatter(x=xnew, y=sp2, mode='lines', name='Spline', line=dict(dash='dot', color='red', width=2)))
# Configurando layout
fig.update_layout(
title='Interpolação de Dados',
xaxis_title='Campo Magnético (Bi)',
yaxis_title='Permeabilidade (mi)',
width=800,
height=600,
legend=dict(title='Métodos de Interpolação', orientation='h', x=0.5, xanchor='center')
)
fig.show()
Q = functionspace(mesh, ("DG", 0))
V = functionspace(mesh, ("Lagrange", 2))
tdim = mesh.topology.dim
# Define the radial coordinate r
# Define the radial coordinate r from x[0]
x = SpatialCoordinate(mesh)
r = x[0]
# Update the weak form to include the magnet
facets = locate_entities_boundary(mesh, tdim - 1, lambda x: np.full(x.shape[1], True))
dofs = locate_dofs_topological(Q, tdim - 1, facets)
bc = dirichletbc(default_scalar_type(0), dofs, Q)
# Update the weak shape to include the magnet
mur_air = 1 # Relative permeability in air
mu0 = 4 * np.pi * 1e-7 # Permeabilidade do vácuo
Hc = float(1000) # Vacuum permeability
# Function to define the relative permeability in the domainio
def mu_r(x):
# Define spatially variable relative permeability
in_magnet = (x[0] <= r_mag) & (np.abs(x[1]) <= z_mag)
return np.where(in_magnet, np.interp(x[0], Bi_steel, mi_steel / mu0), mur_air)
# Convert mu_r to a spatial function in Dolfinx
mu_r_function = Function(V)
mu_r_function.interpolate(mu_r)
# Variational equation for the magnet
a = (1 / mu_r_function) * (1 / r) * dot(grad(u), grad(v)) * dx
# Modify the source term to include Hc
M = Constant(mesh, Hc)
L = mu0 * M * v * dx
# Solve the linear problem
A_ = Function(V)
problem = LinearProblem(a, L, u=A_, bcs=[bc])
problem.solve()
# Calculate the magnetic field
W = functionspace(mesh, ("DG", 0, (mesh.geometry.dim, )))
B = Function(W)
B_expr = Expression(as_vector((-(1/r)*A_.dx(1), (1/r)*A_.dx(0))), W.element.interpolation_points())
B.interpolate(B_expr)
plotter = pyvista.Plotter()
# Converts the grid to an UnstructuredGrid
A_grid = pyvista.UnstructuredGrid(*vtk_mesh(V))
# Add the data from the A_z field
A_grid.point_data["A_"] = A_.x.array
A_grid.set_active_scalars("A_")
# Applies the deformation based on the A_z field
warp = A_grid.warp_by_scalar("A_")
# Add the deformed mesh to the plotter
actor = plotter.add_mesh(warp, show_edges=False)
plotter.view_xy()
plotter.show()
# Iterpolate B again to mach vtk_mesh DoF.
Wl = functionspace(mesh, ("Lagrange", 2, (mesh.geometry.dim, )))
Bl = Function(Wl)
Bl.interpolate(B)
pyvista.start_xvfb()
topology, cell_types, geo = vtk_mesh(V)
values = np.zeros((geo.shape[0], 3), dtype=np.float64)
values[:, :len(Bl)] = Bl.x.array.real.reshape((geo.shape[0], len(Bl)))
# Create a point cloud of glyphs
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geo)
function_grid["Bl"] = values
glyphs = function_grid.glyph(orient="Bl", factor=.5)
# Create a pyvista-grid for the mesh
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
grid = pyvista.UnstructuredGrid(*vtk_mesh(mesh, mesh.topology.dim))
# Create plotter
plotter = pyvista.Plotter()
#plotter.add_mesh(grid, style="wireframe", color="k")
plotter.add_mesh(glyphs)
plotter.view_xy()
#plotter.window_size = [1000, 250];
#plotter.camera.zoom(3)
plotter.show();