ellipt2d 2.1.0

Features

• Supports triangulated domain generated by the _pytriangle: https://github.com/pletzer/pytriangle package

• Stiffness is a tensor

• Can output solution to VTK file

• Simple, easy to use

How to solve an elliptic problem with Ellipt2d

Define the domain. We recommend to use the pytriangle package to triangule a domain. This involves specifying boundary points and segments. If there are holes in the domain, then you’ll need to specify those as well. As an example, we’ll triangulate an annulus:

import triangle
import numpy

# number of outer and inner poloidal points
nto, nti = 16, 8

dto, dti = 2*numpy.pi / nto, 2*numpy.pi / nti
to = numpy.linspace(0., 2*numpy.pi - dto, nto)
ti = numpy.linspace(0., 2*numpy.pi - dti, nti)

# outer boundary, points and segments go counterclockwise
bound_pts = [(numpy.cos(t), numpy.sin(t)) for t in to]
bound_seg = [(i, i + 1) for i in range(nto - 1)] + [(nto - 1, 0)] # close the contour

# add the inner boundary, points and segments go clockwise
bound_pts += [(0.3*numpy.cos(t), 0.3*numpy.sin(t)) for t in ti]
bound_seg += [(i, i+1) for i in range(nto, nto + nti - 1)] + [(nto + nti - 1, nto)]

# now create the triangulation
grid = triangle.Triangle()
grid.set_points(bound_pts)
grid.set_segments(bound_seg)
grid.set_holes([(0., 0.)])
grid.triangulate()

Create an Ellipt2d instance, here a Laplace operator:

import ellipt2d

# - div F . grad u + g = s
# F = [[fxx, fxy], [fxy, fyy]] and g are defined on cells and s on nodes
# here fxx = fyy = 1 and fxy = g = 0
nnodes = grid.get_num_points()
ncells = grid.get_num_triangles()
fxx = fyy = numpy.ones(ncells)
fxy = g = numpy.zeros(ncells)
s = numpy.zeros(nnodes)

equ = ellipt2d.Ellipt2d(grid=grid, fxx=fxx, fxy=fxy, fyy=fyy, g=g, s=s)

Set the boundary conditions. You only need to specify boundary conditions if they are different from the zero flux boundary conditions. Here we specify the Dirichlet boundary conditions on the outer contour (first nto - 1 points):

# gather all the points that are on the external boundary where the radius is larger than 0.31
nodes = grid.get_points() # [(x, y, 0 or 1), ...]
ext_pts = [(i, nodes[i][0][0], nodes[i][0][1]) for i in range(nnodes) if nodes[i][1] == 1 and nodes[i][0][0]**2 + nodes[i][0][1]**2 > 0.31**2]
# Dirichlet BC is cos of angle
db = {bp[0]: numpy.cos(numpy.arctan2(bp[2], bp[1])) for bp in ext_pts}
equ.setDirichletBoundaryConditions(db)

Solve the linear system:

u = equ.solve()

Save the solution in a VTK file for plotting with Paraview or VisIt:

equ.saveVTK(filename='sol.vtk', solution=u, sol_name='u')

Credits

This package was created with Cookiecutter and the audreyr/cookiecutter-pypackage project template.