Gluing together two boundaries

Asked by Artur Palha on 2013-03-01

I am new to FEniCS so probably this has a very simple solution.

The problem is that I wish to solve Poisson's equation in 2D with Dirichlet boundary conditions on a section of a hollow cylinder of outer radius R and inner radius r.

I thought of starting with the Creating More Complex Domains tutorial and generate my mesh. So I take a square and map it to the hollow cylinder. With this mapping, the top and bottom boundaries of the square will be mapped (up to machine precision) to the same nodes.

My question is, how can I glue them together in order to say that the degrees of freedom are the same? I tried using periodic boundary conditions but it did not work. Any help?

This is how I tried to do it:

from dolfin import *
import numpy

# The mesh generation part follow the Creating More Complex Domains tutorial

# this part just generates the mesh in a square and them maps it to a hollow
# cylinder

# ---------------------------------------------------------------------------

# define the inflow and the radius of the cylinder
r = 1 # the radius of the cylinder
R = 5 # the outer radius of the mesh
v_inf = [1,0]

# define the transformation of the mesh into the cylinder
Theta = 2*pi
nr = 4 # divisions in r direction
nt = 8 # divisions in theta direction
mesh = RectangleMesh(0, 0, 1, 1, nr, nt, 'crossed')

# First make a denser mesh towards r=a
x = mesh.coordinates()[:,0]
y = mesh.coordinates()[:,1]
s = 1.3

def denser(x, y):
    return [x**s, y]

x_bar, y_bar = denser(x, y)
xy_bar_coor = numpy.array([x_bar, y_bar]).transpose()
mesh.coordinates()[:] = xy_bar_coor

def cylinder(r, s, rMin, rMax):
    return [(r*(rMax-rMin)+rMin)*numpy.cos(Theta*s), (r*(rMax-rMin)+rMin)*numpy.sin(Theta*s)]

x_hat, y_hat = cylinder(x_bar, y_bar, r, R)
xy_hat_coor = numpy.array([x_hat, y_hat]).transpose()
mesh.coordinates()[:] = xy_hat_coor

plot(mesh, title='hollow cylinder')

# ---------------------------------------------------------------------------

# Now that the mesh is generate, generate the function spaces and the
# analytical solution

# define the function space associated to the mesh
V = FunctionSpace(mesh, 'Lagrange', 1)

# Define analytical solution and boundary conditions
phi_exact = Expression('v_inf_x*x[0]*(1+r*r/(x[0]*x[0] + x[1]*x[1])) + v_inf_y*x[1]*(1+r*r/(x[0]*x[0] + x[1]*x[1]))', v_inf_x=v_inf[0], v_inf_y=v_inf[1], r=r)

# ---------------------------------------------------------------------------

# Generate the subdomains for the Dirichlet boundary conditions (the outer
# and inner surfaces of the cylinder)

# Su domain for Dirichlet boundary condition
class DirichletBoundaryOuter(SubDomain):
    def inside(self, x, on_boundary):
        #return bool((abs(x[0] - 1) < DOLFIN_EPS) and on_boundary)
        return bool((abs(x[0]*x[0]+x[1]*x[1] - R*R) < DOLFIN_EPS) and on_boundary)

class DirichletBoundaryInner(SubDomain):
    def inside(self, x, on_boundary):
        return bool((abs(x[0]*x[0]+x[1]*x[1] - r*r) < DOLFIN_EPS) and on_boundary)
        #return bool((abs(x[0] + 1) < DOLFIN_EPS) and on_boundary)

# ---------------------------------------------------------------------------

# Generate the subdomain for periodic boundary conditions

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool((abs(x[1]) < 100*DOLFIN_EPS) and (x[0] > -0.1) and on_boundary)

    def map(self, x, y):
        y[0] = x[0]
        y[1] = x[1]

# ---------------------------------------------------------------------------

# Create the boundary conditions

# Create Dirichlet boundary condition
dbc = DirichletBoundaryOuter()
bc0 = DirichletBC(V, phi_exact, dbc)

dbc = DirichletBoundaryInner()
bc1 = DirichletBC(V, phi_exact, dbc)

# Create periodic boundary condition
pbc = PeriodicBoundary()
bc2 = PeriodicBC(V, pbc)

# Collect boundary conditions
bcs = [bc0, bc1, bc2]

# ---------------------------------------------------------------------------

# Simple Poisson problem

# Define variational problem
w = TrialFunction(V) # the trial functions
v = TestFunction(V) # the test functions
f = Constant(0)
a = inner(nabla_grad(w), nabla_grad(v))*dx
L = inner(f,v)*dx

# Compute solution
phi = Function(V) # the solution is a function defined as a linear combination of V functions

solve(a == L, phi, bcs)

interactive()

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Artur Palha
Solved:
2013-03-01
Last query:
2013-03-01
Last reply:
2013-03-01
Anders Logg (logg) said : #1

If you are using version 1.1, you should be able to generate your mesh
using CSG (Constructive Solid Geometry). Look at the CSG demos in
DOLFIN. Essentially, you can do

c = Cylinder(..., r)
C = Cylinder(..., R)
g = C - c
m = Mesh(g)

--
Anders

On Fri, Mar 01, 2013 at 10:56:06AM -0000, Artur Palha wrote:
> New question #223147 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/223147
>
> I am new to FEniCS so probably this has a very simple solution.
>
> The problem is that I wish to solve Poisson's equation in 2D with Dirichlet boundary conditions on a section of a hollow cylinder of outer radius R and inner radius r.
>
> I thought of starting with the Creating More Complex Domains tutorial and generate my mesh. So I take a square and map it to the hollow cylinder. With this mapping, the top and bottom boundaries of the square will be mapped (up to machine precision) to the same nodes.
>
> My question is, how can I glue them together in order to say that the degrees of freedom are the same? I tried using periodic boundary conditions but it did not work. Any help?
>
> This is how I tried to do it:
>
>
> from dolfin import *
> import numpy
>
> # The mesh generation part follow the Creating More Complex Domains tutorial
>
> # this part just generates the mesh in a square and them maps it to a hollow
> # cylinder
>
> # ---------------------------------------------------------------------------
>
> # define the inflow and the radius of the cylinder
> r = 1 # the radius of the cylinder
> R = 5 # the outer radius of the mesh
> v_inf = [1,0]
>
> # define the transformation of the mesh into the cylinder
> Theta = 2*pi
> nr = 4 # divisions in r direction
> nt = 8 # divisions in theta direction
> mesh = RectangleMesh(0, 0, 1, 1, nr, nt, 'crossed')
>
> # First make a denser mesh towards r=a
> x = mesh.coordinates()[:,0]
> y = mesh.coordinates()[:,1]
> s = 1.3
>
> def denser(x, y):
> return [x**s, y]
>
> x_bar, y_bar = denser(x, y)
> xy_bar_coor = numpy.array([x_bar, y_bar]).transpose()
> mesh.coordinates()[:] = xy_bar_coor
>
> def cylinder(r, s, rMin, rMax):
> return [(r*(rMax-rMin)+rMin)*numpy.cos(Theta*s), (r*(rMax-rMin)+rMin)*numpy.sin(Theta*s)]
>
> x_hat, y_hat = cylinder(x_bar, y_bar, r, R)
> xy_hat_coor = numpy.array([x_hat, y_hat]).transpose()
> mesh.coordinates()[:] = xy_hat_coor
>
> plot(mesh, title='hollow cylinder')
>
> # ---------------------------------------------------------------------------
>
> # Now that the mesh is generate, generate the function spaces and the
> # analytical solution
>
> # define the function space associated to the mesh
> V = FunctionSpace(mesh, 'Lagrange', 1)
>
> # Define analytical solution and boundary conditions
> phi_exact = Expression('v_inf_x*x[0]*(1+r*r/(x[0]*x[0] + x[1]*x[1])) + v_inf_y*x[1]*(1+r*r/(x[0]*x[0] + x[1]*x[1]))', v_inf_x=v_inf[0], v_inf_y=v_inf[1], r=r)
>
>
> # ---------------------------------------------------------------------------
>
> # Generate the subdomains for the Dirichlet boundary conditions (the outer
> # and inner surfaces of the cylinder)
>
> # Su domain for Dirichlet boundary condition
> class DirichletBoundaryOuter(SubDomain):
> def inside(self, x, on_boundary):
> #return bool((abs(x[0] - 1) < DOLFIN_EPS) and on_boundary)
> return bool((abs(x[0]*x[0]+x[1]*x[1] - R*R) < DOLFIN_EPS) and on_boundary)
>
> class DirichletBoundaryInner(SubDomain):
> def inside(self, x, on_boundary):
> return bool((abs(x[0]*x[0]+x[1]*x[1] - r*r) < DOLFIN_EPS) and on_boundary)
> #return bool((abs(x[0] + 1) < DOLFIN_EPS) and on_boundary)
>
> # ---------------------------------------------------------------------------
>
> # Generate the subdomain for periodic boundary conditions
>
> # Sub domain for Periodic boundary condition
> class PeriodicBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return bool((abs(x[1]) < 100*DOLFIN_EPS) and (x[0] > -0.1) and on_boundary)
>
> def map(self, x, y):
> y[0] = x[0]
> y[1] = x[1]
>
> # ---------------------------------------------------------------------------
>
> # Create the boundary conditions
>
> # Create Dirichlet boundary condition
> dbc = DirichletBoundaryOuter()
> bc0 = DirichletBC(V, phi_exact, dbc)
>
> dbc = DirichletBoundaryInner()
> bc1 = DirichletBC(V, phi_exact, dbc)
>
> # Create periodic boundary condition
> pbc = PeriodicBoundary()
> bc2 = PeriodicBC(V, pbc)
>
> # Collect boundary conditions
> bcs = [bc0, bc1, bc2]
>
> # ---------------------------------------------------------------------------
>
> # Simple Poisson problem
>
> # Define variational problem
> w = TrialFunction(V) # the trial functions
> v = TestFunction(V) # the test functions
> f = Constant(0)
> a = inner(nabla_grad(w), nabla_grad(v))*dx
> L = inner(f,v)*dx
>
> # Compute solution
> phi = Function(V) # the solution is a function defined as a linear combination of V functions
>
> solve(a == L, phi, bcs)
>
> interactive()
>
>
>

Artur Palha (artur-palha) said : #2

Thanks a lot Anders. This solves the problem of generating the mesh with the hole. But I wanted to generate a structured mesh. Do you know if that is possible within Dolfin, or should I use a mesh generator and import the mesh to Dolfin?

Once again thanks a lot.

-artur palha