A package to solve second order elliptic problems in 2D
Project description
ellipt2d
A package for solving second order elliptic problems in 2D
Free software: MIT license
Documentation: https://ellipt2d.readthedocs.io.
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.
History
2.0.0 (2020-12-15)
First release on PyPI.
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.