fenics team mailing list archive

[Harish Narayanan <harish.mlists@xxxxxxxxx>]

Dear fellow PDE solvers,

The sort of question I am about the phrase might seem strange to some,
but perhaps it will seem natural to others.

Last evening, I was curious about meshing quirky geometries like the
Möbius strip in DOLFIN. This morning, I was curious about solving a PDE
on such a domain, and so I came up with the attached python script.

The script purportedly sets up a Poisson equation (with homogeneous
boundary conditions and a constant source) on a Möbius strip (with four
twists) and "solves it." I see a pretty looking plot with some
psychedelic colours, but can someone who knows more about this stuff be
kind enough to clarify what I've ended up doing?

Harish

from numpy import *
from dolfin import *

# Mobius strip parameters
# Number of times you'd like your strip to twist (only even numbers allowed!)
num_twists = 4
# The width of the strip
width = 0.5

# Mesh parameters
# Number of nodes along the length of the strip
nl = 300
# Number of nodes along the width of the strip (>= 2)
nw = 14

# Generate suitable ranges for parameters
u_range = arange(nl, dtype='d')/(nl)*2*pi
v_range = arange(nw, dtype='d')/(nw - 1.0)*width

# Create the mesh
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh, "triangle", 2, 3)

# Create a list to store the vertices
editor.init_vertices(nl*nw)

# Populate the list of vertices
j = 0
for u in u_range:
    for v in v_range:
        editor.add_vertex(j,  cos(u) + v*cos(num_twists*u/2.0)*cos(u), \
                              sin(u) + v*cos(num_twists*u/2.0)*sin(u), \
                              v*sin(num_twists*u/2.0))
        j = j + 1

# Create a list to store the cells
editor.init_cells(nl*(nw - 1)*2)

# Populate the list of cells
k = 0
for i in range(nl - 1):
    for j in range(nw - 1):
        editor.add_cell(k    , i*nw + j, (i + 1)*nw + j + 1, i*nw + j + 1) 
        editor.add_cell(k + 1, i*nw + j, (i + 1)*nw + j    , (i + 1)*nw + j + 1)
        k = k + 2
# Close the geometry
for j in range(nw - 1):
    editor.add_cell(k    , (nl - 1)*nw + j, j + 1, (nl - 1)*nw + j + 1) 
    editor.add_cell(k + 1, (nl - 1)*nw + j, j    , j + 1)
    k = k + 2

editor.close()

# Setup and solve the Poisson's equation
V = FunctionSpace(mesh, "CG", 1)
f = Expression("1.0")
zero = Constant(0.0)

boundary = compile_subdomains(["on_boundary"])
bc = DirichletBC(V, zero, boundary) 

v = TestFunction(V)
u = TrialFunction(V)

a = inner(grad(v), grad(u))*dx
L = v*f*dx

problem = VariationalProblem(a, L, bc)
u = problem.solve()

plot(u, mesh = mesh, interactive = True)

file = File("mobius.pvd")
file << u