Quarter annulus¶

In this notebook, we present a set of solutions for 2D using Elementally.

The Geoemetry¶

We'll consider a 2D quarter annulus starting at $\theta=0$ going to $\theta=\pi$, with an inner diameter of 2.0 and an outer diameter of 3.0

Setting up the Python environment¶

Throughout this example we'll use the %pylab magic (populates the environment with e.g. numpy as np and pyplot as plt) and we'll get the modules from elementally.

In [1]:

%pylab
%matplotlib inline
import elementally.meshers as meshers
import elementally.assemblers as assemblers
import elementally.boundaries as boundaries
# sys.path.insert(0, '../elementally')

Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

/usr/lib/python2.7/dist-packages/pkg_resources.py:1031: UserWarning: /home/tarjei/.python-eggs is writable by group/others and vulnerable to attack when used with get_resource_filename. Consider a more secure location (set with .set_extraction_path or the PYTHON_EGG_CACHE environment variable).
  warnings.warn(msg, UserWarning)

Creating the mesh¶

The mesh itself is created using Andreas Kloeckner terrific (Meshpy)[http://mathema.tician.de/software/meshpy/] package, wrapped in an Elementally front-end.

In [2]:

# Creating a mesh on the form of a quarter annulus, centred in origo, with an outer diameter of 4. and an inner of 3.
mesh = meshers.quarter_annulus_2d(0.3, 0, 3.14, (0, 0), 3, 4)
# Remapping the C++ arrayed points from the Meshpy mesh over to a Numpy array
points = np.array(mesh.points)
# Creating plotable trianfles on the form of a (n,3) matrix, where n is the number of elements
triangles = np.array(mesh.elements)

In [3]:

plt.triplot(points[:, 0], points[:,1], triangles)
plt.show()

The Poisson Equation¶

Firstly, let's solve the Poisson equation on our established mesh.

Boundary Conditions¶

With the equation written as $\Delta u = f,$ we'll here solve for $u$ using a simple "Poisson function" $f$ of $f=1$. Further, we'll impose two kinds of boundary conditions, Dirichlet and Neuman boundaries. We can sum this up as follows:

In [4]:

# Poisson function:
def f(x):
    return 1.

# Dirichlet function:
def g(x):
    return 0.

# Neumann function:
def h(x):
    return 0.

The boundary conditions themselves are set up in a dictionary called boundary_dict, which is a nested dictionary (here) containing two main keys (vis. dir and neu) where the "facet markers" of each facet-zone (set in quarter_annulus_2d) is set to a corresponding function handle.

In [5]:

boundary_dict = {'dir': {2: g, 3: g}, 'neu': {1: h}}

Assembly of the matrices¶

There's a separate module for the assembling of matrices in elementally. We can get the desired matrices through the following lines.

In [6]:

assembly = assemblers.Poisson_2d(mesh, f)
A = assembly.A
b = assembly.b
points = assembly.points

Imposing the boundary conditions¶

We impose the boundary conditions on the assembled matrices through the boundaries module.

In [7]:

A,b = boundaries.impose_dirichlet(boundary_dict, mesh, b, A)

Solving the equations¶

We can now find $u$ using numpy's eminent linalg module.

In [8]:

import numpy.linalg as la
u = la.solve(A, b)

Plotting the solution¶

Finally, we can visualize the resulting solution through the following matplotlib objects.

In [9]:

import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

fig1 = plt.figure(1)
ax = fig1.gca(projection='3d')

ax.plot_trisurf(points[:,0], points[:,1], u,
                triangles = mesh.elements,
                cmap=cm.jet, linewidth=0.2)

plt.show(1)