Difficulty in reproducing the behavior of a non-linear magnet (from phenics to dolfinx)

dolfinx
1

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();
1 Like
2

Could you be more specific about what you are unsure about?
Is there a specific plot that doesn’t look right to you after simulation?

If the initial plots are fine, could you remove what is not needed from the posts to make them easier to digest?

Could you be more specific, what doubts are you having?

3

Your initial attempts look really good, but as @dokken says we need a little more context and precision in your question to be able to help.

4

Thank you for answering me, dokken and @nate

Below, I’ve posted the code snippet I’m trying to reproduce and the corresponding graphics. I think they might help you understand my problem.

Old code

# 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)

first
first539×738 112 KB

# 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)

B
B404×728 113 KB

norm_B
norm_B528×729 91.3 KB

New code

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()

A_
A_671×641 35.2 KB

# 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();

B1
B1712×509 48.2 KB

5

Dear @eegdo,

There are a few things that puzzles me with your DOLFINx code above.
In the provided code, you have not defined u or v, which means that the code can be executed.

Secondly, i would use another mu_r, i.e.

def mu_r(x):
    return np.interp(x[0], Bi_steel, mi_steel / mu0)


MU_space = functionspace(mesh, ("DG", 2))
mu_r_function = Function(MU_space)
mu_r_function.x.array[:] = mu0
mu_r_function.interpolate(mu_r, cells0=cell_tags.find(inner_tag))
mu_r_function.x.scatter_forward()

# Convert mu_r to a spatial function in Dolfinx

from dolfinx.io import VTXWriter

with VTXWriter(MPI.COMM_WORLD, "mu_r.bp", [mu_r_function], engine="BP5") as bp:
    bp.write(0.0)

which yields

image
image674×737 12.7 KB

I would also set solver options for your linear system:

problem = LinearProblem(
    a,
    L,
    u=A_,
    bcs=[bc],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
    },
)
problem.solve()

which yields:

image
image725×737 21.8 KB

image
image725×737 28.8 KB

as one now can observe, the boundary conditions are different.

If you then change the space you set bcs on

facets = locate_entities_boundary(mesh, tdim - 1, lambda x: np.full(x.shape[1], True))
dofs = locate_dofs_topological(V, tdim - 1, facets)
bc = dirichletbc(default_scalar_type(0), dofs, V)

you get

image
image725×737 15.1 KB

image
image725×737 27.3 KB
os
which is closer to what you have.
The full code I used

# 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

# Visualization
import pyvista
from mpi4py import MPI

# FEniCS
from ufl import (
    SpatialCoordinate,
    TestFunction,
    TrialFunction,
    dx,
    dot,
    grad,
    as_vector,
)
from dolfinx import default_scalar_type
from dolfinx.io import gmshio
from dolfinx.fem import (
    Constant,
    Expression,
    Function,
    functionspace,
    dirichletbc,
    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


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)

V = functionspace(mesh, ("Lagrange", 2))
tdim = mesh.topology.dim
u = TrialFunction(V)
v = TestFunction(V)
# 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(V, tdim - 1, facets)
bc = dirichletbc(default_scalar_type(0), dofs, V)

# 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):
    return np.interp(x[0], Bi_steel, mi_steel / mu0)


MU_space = functionspace(mesh, ("DG", 2))
mu_r_function = Function(MU_space)
mu_r_function.x.array[:] = mu0
mu_r_function.interpolate(mu_r, cells0=cell_tags.find(inner_tag))
mu_r_function.x.scatter_forward()

# Convert mu_r to a spatial function in Dolfinx

from dolfinx.io import VTXWriter

with VTXWriter(MPI.COMM_WORLD, "mu_r.bp", [mu_r_function], engine="BP5") as xdmf:
    xdmf.write(0.0)

# 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],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
    },
)
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)

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=0.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()
3 Likes
6

Thank you a lot, dokken. You helped me to understand better what was going on. The last question that I have is about the possibility to calculate on the boundary (more specifally on the right side). I tried to include ds in the code you provided, like it is done in the tutorial, but I don’t know how to integrate it in the rest of the code. I think it is the reason your plots were different from the first 3 I posted.

V = functionspace(mesh, ("Lagrange", 2))
tdim = mesh.topology.dim
u = TrialFunction(V)
v = TestFunction(V)
# 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))

facet_markers = np.ones(np.shape(facets),dtype=np.int32)
facet_tag = meshtags(mesh, gdim-1, facets, facet_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

dofs = locate_dofs_topological(V, tdim - 1, facets)
bc = dirichletbc(default_scalar_type(0), dofs, V)

# 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
7

Dear @eegdo,

It seems like you are almost there. However, there are a few deviations compared to what you did in legacy FEniCS

You’d next like to apply Hc to one side of this interface (Which side specifically, the core, tagged with 1?)

You can do this by following: Add one-sided integration example as example of custom integration · Issue #158 · jorgensd/dolfinx-tutorial · GitHub

Here is a fully working example:

# Basic Libraries
import ufl
import dolfinx
from dolfinx.io import VTXWriter
import numpy as np  # For vector and matrix manipulation
import gmsh  # Interface with the software GMSH
from scipy.interpolate import interp1d, splev, splrep

# Visualization
import pyvista
from mpi4py import MPI

# FEniCS
from ufl import (
    SpatialCoordinate,
    TestFunction,
    TrialFunction,
    dx,
    dot,
    grad,
    as_vector,
)
from dolfinx import default_scalar_type
from dolfinx.io import gmshio
from dolfinx.fem import (
    Constant,
    Expression,
    Function,
    functionspace,
    dirichletbc,
    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


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)

V = functionspace(mesh, ("Lagrange", 2))
tdim = mesh.topology.dim
u = TrialFunction(V)
v = TestFunction(V)
# 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(V, tdim - 1, facets)
bc = dirichletbc(default_scalar_type(0), dofs, V)

# 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


def ima(x):
    return np.isclose(x[0], r_mag) & (x[1] <= z_mag) & (x[1] >= -z_mag)


interal_ima_facets = dolfinx.mesh.locate_entities(
    mesh, mesh.topology.dim - 1, ima)


integration_entities = []
facet_map = mesh.topology.index_map(mesh.topology.dim-1)
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
f_to_c = mesh.topology.connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-1)

for i, facet in enumerate(interal_ima_facets):
    # Only loop over facets owned by the process to avoid duplicate integration
    if facet >= facet_map.size_local:
        continue
    # Find cells connected to facet
    cells = f_to_c.links(facet)
    # Get value of cells
    marked_cells = cell_tags.values[cells]
    # Get the cell marked with 1
    correct_cell = np.flatnonzero(marked_cells == 1)

    assert len(correct_cell) == 1
    # Get local index of facet
    local_facets = c_to_f.links(cells[correct_cell[0]])
    local_index = np.flatnonzero(local_facets == facet)
    assert len(local_index) == 1

    # Append integration entities
    integration_entities.append(cells[correct_cell[0]])
    integration_entities.append(local_index[0])

ds = ufl.Measure("ds", domain=mesh, subdomain_data=[
                 (8, np.asarray(integration_entities, dtype=np.int32))])

# Function to define the relative permeability in the domainio


def mu_r(x):
    return np.interp(x[0], Bi_steel, mi_steel / mu0)


MU_space = functionspace(mesh, ("DG", 2))
mu_r_function = Function(MU_space)
mu_r_function.x.array[:] = mu0
mu_r_function.interpolate(mu_r, cells0=cell_tags.find(inner_tag))
mu_r_function.x.scatter_forward()

# Convert mu_r to a spatial function in Dolfinx


with VTXWriter(MPI.COMM_WORLD, "mu_r.bp", [mu_r_function], engine="BP5") as xdmf:
    xdmf.write(0.0)

# 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
L += Hc * v * ds(8)

# Solve the linear problem
A_ = Function(V)
problem = LinearProblem(
    a,
    L,
    u=A_,
    bcs=[bc],
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
        "ksp_error_if_not_converged": True,
    },
)
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)

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=0.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()

yielding:

image
image1132×802 60.3 KB

image
image1818×1943 488 KB

You would have to define which side of the interior facet is correct, i.e.

    # Get the cell marked with 1
    correct_cell = np.flatnonzero(marked_cells == 1)

is it those marked with 1 or 2?

It would also be good to get some justification as to where the term:

comes from.
Is it that \frac{1}{\mu r}(\nabla u\vert_1 - \nabla u \vert_2) \cdot n_1 = Hc ?