How can I reproduce the mesh used in the magnetostatic example?

mesh
1

I have tried to reproduce the nice mesh shown in Fig. 4.3 from
Solving PDEs in Python The FEniCS Tutorial I - p. 103 (109):

original
original822×850 670 KB

I have used the provided code:

from fenics import *
from mshr import *
from math import sin, cos, pi

a = 1.0 # inner radius of iron cylinder
b = 1.2 # outer radius of iron cylinder
c_1 = 0.8 # radius for inner circle of copper wires
c_2 = 1.4 # radius for outer circle of copper wires
r = 0.1 # radius of copper wires
R = 5.0 # radius of domain
n = 10 # number of windings

# Define geometry for background
domain = Circle(Point(0, 0), R)

# Define geometry for iron cylinder
cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)

# Define geometry for wires (N = North (up), S = South (down))
angles_N = [i*2*pi/n for i in range(n)]
angles_S = [(i + 0.5)*2*pi/n for i in range(n)]
wires_N = [Circle(Point(c_1*cos(v), c_1*sin(v)), r) for v in angles_N]
wires_S = [Circle(Point(c_2*cos(v), c_2*sin(v)), r) for v in angles_S]

# Set subdomain for iron cylinder
domain.set_subdomain(1, cylinder)

# Set subdomains for wires
for (i, wire) in enumerate(wires_N):
    domain.set_subdomain(2 + i, wire)
for (i, wire) in enumerate(wires_S):
    domain.set_subdomain(2 + n + i, wire)

# Create mesh
mesh = generate_mesh(domain, 32)

and visualised the besh:

plot(mesh)
import matplotlib.pyplot as plt
plt.xlim(-1.6, 1.6)
plt.ylim(-1.6, 1.6)

but results are disappointing:
Mshr.3

How can I reproduce the mesh? Also, is it possible to generalise that approach to 3D meshes?

2

You can try the following

from fenics import *
from mshr import *
from math import sin, cos, pi
import matplotlib.pyplot as plt

a = 1.0 # inner radius of iron cylinder
b = 1.2 # outer radius of iron cylinder
c_1 = 0.8 # radius for inner circle of copper wires
c_2 = 1.4 # radius for outer circle of copper wires
r = 0.1 # radius of copper wires
R = 5.0 # radius of domain
n = 10 # number of windings

# Define geometry for background
domain = Circle(Point(0, 0), R)

# Define geometry for iron cylinder
cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)

# Define geometry for wires (N = North (up), S = South (down))
angles_N = [i*2*pi/n for i in range(n)]
angles_S = [(i + 0.5)*2*pi/n for i in range(n)]
wires_N = [Circle(Point(c_1*cos(v), c_1*sin(v)), r, 100) for v in angles_N]
wires_S = [Circle(Point(c_2*cos(v), c_2*sin(v)), r, 100) for v in angles_S]

# Set subdomain for iron cylinder
domain.set_subdomain(1, cylinder)

# Set subdomains for wires
for (i, wire) in enumerate(wires_N):
    domain.set_subdomain(2 + i, wire)
for (i, wire) in enumerate(wires_S):
    domain.set_subdomain(2 + n + i, wire)

# Create mesh
mesh = generate_mesh(domain, 32)

plot(mesh)
plt.show()

The additional parameter in the Circle constructor defines the number of segments for the boundary approximation. The default seems to be set to 32 though, so i do not know why it does not work properly when you do not set the parameter.

3

To generalize this approach to 3D, i would strongly suggest using external sotware such as pygmsh/gmsh. This is mostly due to the fact that the support for mshr in the later years have become very limited.

1 Like
4

Thanks, it works.

In my version of dolfin (2018.1.0) segments defaults to 0. By setting it to 32 I got:

Mshr.3

5

On the visualization part check out vtkplotter shorcut function plot():

591×590 118 KB

1 Like