%matplotlib inline
from ipywidgets import interact
import numpy as np
import ufl
from dolfin import *
parameters["plotting_backend"] = "matplotlib"
import logging
logging.getLogger("FFC").setLevel(logging.WARNING)
import matplotlib.pyplot as plt
def create_line_mesh(vertices):
"Given list of vertex coordinate tuples, build and return a mesh of intervals."
# Get dimensions
gdim = 1 if isinstance(vertices[0], float) else len(vertices[0])
tdim = 1
# Automatic choice of cellname for simplices
cellname = "interval"
# Indirect error checking and determination of tdim via ufl
ufl_cell = ufl.Cell(cellname, gdim)
assert tdim == ufl_cell.topological_dimension()
# Create mesh to return
mesh = Mesh()
# Open mesh in editor
me = MeshEditor()
me.open(mesh, cellname, tdim, gdim)
# Add vertices to mesh
nv = len(vertices)
me.init_vertices(nv)
if gdim == 1:
for i, v in enumerate(vertices):
me.add_vertex(i, v)
else:
for i, v in enumerate(vertices):
me.add_vertex(i, *v)
# Add cells to mesh
me.init_cells(nv-1)
for i in range(nv-1):
c = (i, i+1)
me.add_cell(i, *c)
me.close()
return mesh
def line1d():
n = 100
us = [i/float(n-1) for i in range(n)]
vertices = [u**3 for u in us]
return create_line_mesh(vertices)
def line2d():
n = 100
us = [i/float(n-1) for i in range(n)]
vertices = [(cos(1.5*DOLFIN_PI*u),
sin(1.5*DOLFIN_PI*u))
for u in us]
return create_line_mesh(vertices)
def line3d():
n = 100
us = [i/float(n-1) for i in range(n)]
vertices = [(cos(4.0*DOLFIN_PI*u),
sin(4.0*DOLFIN_PI*u),
2.0*u)
for u in us]
return create_line_mesh(vertices)
plot(UnitSquareMesh(5,5))
[<matplotlib.lines.Line2D at 0x7f496cf22dd0>, <matplotlib.lines.Line2D at 0x7f496cf22e90>]
plot(UnitCubeMesh(5,5,5))
mesh = UnitIntervalMesh(5)
mesh.coordinates()[:] = mesh.coordinates()[:]**2
plot(mesh)
[<matplotlib.lines.Line2D at 0x7f496b844e10>]
plot(line2d())
[<matplotlib.lines.Line2D at 0x7f496b7ced90>]
plot(line3d())
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x7f496b724b90>]
mesh = UnitSquareMesh(5,5)
cf = CellFunctionSizet(mesh)
for i in range(len(cf)):
cf[i] = i
plot(cf)
<matplotlib.collections.PolyCollection at 0x7f496a14f8d0>
mesh = UnitIntervalMesh(8)
V = FunctionSpace(mesh, "CG", 1)
f = interpolate(Expression("x[0]*x[0]"), V)
plot(f, title="f")
[<matplotlib.lines.Line2D at 0x7f4969ae2e50>]
mesh = UnitIntervalMesh(8)
V = FunctionSpace(mesh, "DG", 0)
f = interpolate(Expression("x[0]*x[0]"), V)
plot(f, title="f")
[<matplotlib.lines.Line2D at 0x7f496986ef10>]
mesh = UnitCubeMesh(5,5,5)
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.interpolate(Expression("x[1]"))
plot(f)
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7f496928c8d0>
mesh = UnitSquareMesh(32, 32)
e = Expression("sin(k*DOLFIN_PI*x[0]*x[1] + q*DOLFIN_PI*x[0] + r*DOLFIN_PI*x[0])", k=8, q=2, r=2)
V0 = FunctionSpace(mesh, "DG", 0)
f0 = Function(V0)
f0.interpolate(e)
V1 = FunctionSpace(mesh, "CG", 1)
f1 = Function(V1)
f1.interpolate(e)
V2 = FunctionSpace(mesh, "CG", 2)
f2 = Function(V2)
f2.interpolate(e)
plot(f0, title="piecewise constant function")
<matplotlib.collections.PolyCollection at 0x7f4968840650>
plot(f1)
<matplotlib.collections.TriMesh at 0x7f496ce25090>
plot(f2)
<matplotlib.collections.TriMesh at 0x7f49686227d0>
def display(n=(4, 80), k=(0,5), d=(0,2)):
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "CG" if d else "DG", d)
e = Expression("sin(k*DOLFIN_PI*x[0]*x[1])", k=4)
e.k = k
f = interpolate(e, V)
plot(f, title="n={n}, k={k}, d={d}".format(**locals()))
interact(display)
<function __main__.display>
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.interpolate(Expression("x[0]*x[0] + x[1]*x[1]*(x[1]-0.5)"))
plot(f, title="$f(x)$", mode="surface")
<matplotlib.collections.TriMesh at 0x7f496651bfd0>
plot(f, title="$f(x)$", mode="wireframe")
[<matplotlib.lines.Line2D at 0x7f49691b2650>, <matplotlib.lines.Line2D at 0x7f49691b4090>]
plot(f, title="$f(x)$", mode="contour")
<matplotlib.tri.tricontour.TriContourSet instance at 0x7f4966d33bd8>
# Not working:
#plot(f, title="$f(x)$", mode="warp")
mesh = UnitSquareMesh(16, 16)
ev = Expression(("x[1]", "-x[0]"))
U0 = VectorFunctionSpace(mesh, "DG", 0)
g0 = Function(U0)
g0.interpolate(ev)
plot(g0)
<matplotlib.quiver.Quiver at 0x7f496699c450>
mesh = UnitSquareMesh(10, 10)
ev = Expression(("x[1]-0.5", "x[0]-0.5"))
U2 = VectorFunctionSpace(mesh, "CG", 2)
g2 = Function(U2)
g2.interpolate(ev)
plot(g2)
<matplotlib.quiver.Quiver at 0x7f496570f090>
mesh = UnitSquareMesh(20, 20)
ev = Expression(("sin(k*3.14159*(x[1]-0.5))", "cos(k*3.14159*(x[0]-0.5))"), k=2)
ev.k = 1
U1 = VectorFunctionSpace(mesh, "CG", 1)
g1 = Function(U1)
def display(k=(1,5)):
ev.k = k
g1.interpolate(ev)
plot(g1)
interact(display)
<function __main__.display>
mesh = UnitCubeMesh(8, 8, 8)
V = VectorFunctionSpace(mesh, "CG", 1)
e = Expression(("x[1]", "-x[0]", "x[2]"))
f = interpolate(e, V)
plot(f)
<mpl_toolkits.mplot3d.art3d.Line3DCollection at 0x7f49650dfa10>
mesh = UnitSquareMesh(9, 9)
ev = Expression(("0.1*(x[1]-0.5)*(x[1]-0.5)", "0.2*x[0] - 0.1*x[0]*x[0]"))
U1 = VectorFunctionSpace(mesh, "CG", 1)
g = Function(U1)
g.interpolate(ev)
plot(g, mode="deformation")
<matplotlib.collections.TriMesh at 0x7f496489d250>