The FEniCS Project (2003-current, www.fenicsproject.org) is a collaborative open-source scientific computing software project targeting automated solution of differential equation via finite element methods. Here, we will illustrate how to take advantage of it for rapid implementation and evaluation of an adaptive finite element method for the Poisson equation.
For simplicity, we will consider the Poisson equation with homogeneous Dirichlet boundary conditions. For a computational domain $\Omega$ with boundary $\partial \Omega$ and given data $f$ and $A$ find $u$ such that
\begin{align} \quad -\mathrm{div} \, A \, \mathrm{grad} \, u &= f \quad \quad in \, \Omega, \\ \quad u &= 0 \quad \quad on \, \partial \Omega. \end{align}In variational form, this problem reads: find $u \in H^1_0(\Omega)$ such that
\begin{equation} \quad \langle A \, \mathrm{grad} \, u, \mathrm{grad} \, v \rangle = \langle f, v \rangle \end{equation}for all $v \in H^1_0(\Omega)$, where $\langle \cdot, \cdot \rangle$ denotes the $L^2(\Omega)$ inner product.
Next, assume that you have a mesh $\mathcal{T}_h$ of the domain and a finite element space $V_h$ defined relative to this mesh. The discrete variational problem then reads: find $u_h \in V_h$ such that
\begin{equation} \quad \langle A \, \mathrm{grad} \, u_h, \mathrm{grad} \, v \rangle = \langle f, v \rangle \end{equation}for all $v \in V_h$.
FEniCS provides both a Python and a C++ interface: here we use the Python interface (obviously). To use FEniCS from Python, import it like this (assuming that you have a working FEniCS installation):
from fenics import *
Next usual step is to think about which domain we would like to consider. Let's start by considering a standard unit square: $\Omega = [0, 1]^2$ tessellated into $n \times n$ boxes and each box divided into two triangles:
# Define mesh and finite element space
n = 3
mesh = UnitCubeMesh(n, n, n)
Next, what finite elements would we like to use? Let's start with Lagrange elements (continuous piecewise polynomials) for instance of order $1$:
V = FunctionSpace(mesh, "Lagrange", 1)
DEBUG:FFC:Reusing form from cache.
Next step is to define the unknown trial function, test function and the data. We do this as follows:
# Define trial function (the unknown)
u = TrialFunction(V)
# Define test function
v = TestFunction(V)
# Define the source function
f = Expression("sin(4*pi*x[2])")
# Define A:
A = Constant(1.0)
# Define the boundary condition
bc = DirichletBC(V, 0.0, "on_boundary")
Ok, now we have all that we need to define the variational problem:
rhs = inner(A*grad(u), grad(v))*dx
lhs = inner(f, v)*dx
Next, we can solve it and plug result into a new Function u:
u = Function(V)
solve(rhs == lhs, u, bc)
print u.vector().array()
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
[ 0. 0. 0. 0. 0. 0. 0.0096225 0. 0. 0. 0. 0. 0. 0. 0. -0.0096225 0. 0.0096225 0. 0. 0. 0. 0.0096225 0. 0. 0. 0.0096225 0. 0. 0. 0. 0. 0. -0.0096225 0. 0. 0. -0.0096225 0. 0. 0. -0.0096225 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]
For plotting inline in this IPython Notebook, we can include this handy code from Chris Richardson:
%matplotlib inline
from dolfin import *
import matplotlib.pyplot as plt
import matplotlib.tri as tri
def mesh2triang(mesh):
xy = mesh.coordinates()
return tri.Triangulation(xy[:, 0], xy[:, 1], mesh.cells())
def plot_inline(obj):
plt.gca().set_aspect('equal')
if isinstance(obj, Function):
mesh = obj.function_space().mesh()
if (mesh.geometry().dim() != 2):
raise(AttributeError)
if obj.vector().size() == mesh.num_cells():
C = obj.vector().array()
plt.tripcolor(mesh2triang(mesh), C)
else:
C = obj.compute_vertex_values(mesh)
plt.tripcolor(mesh2triang(mesh), C, shading='gouraud')
elif isinstance(obj, Mesh):
if (obj.geometry().dim() != 2):
raise(AttributeError)
plt.triplot(mesh2triang(obj), color='k')
Now we can look at the solution and the mesh inline (in this notebook):
plot_inline(u)
Let's consider a version of the example of Nochetto and Veeser, Primer of Adaptive Finite Element Methods, page 14:
Let $\Omega$ be an L-shaped domain with coordinates $x$ and consider the Poisson problem (with $A = 1$) with solution
\begin{equation} u(r, \theta) = r^{2/3} \sin(2 \theta/3) \end{equation}where $r^2 = x[0]^2 + x[1]^2$ and $\theta = \arctan(x[1]/x[0])$.
Let's define a function that given a mesh solves this test case on that mesh first...
k = 1
parameters["form_compiler"]["quadrature_degree"] = 2*k + 1
def poisson(mesh):
# Define your favorite finite element space
V = FunctionSpace(mesh, "CG", k)
# Define the exact solution and compute the
# right hand side f symbolically
x = SpatialCoordinate(mesh)
r = sqrt(x[0]**2 + x[1]**2)
theta = atan_2(x[1], x[0])
u_exact = pow(r, 2./3)*sin(2./3*theta)
f = - div(grad(u_exact))
# Use the exact solution as the boundary condition
bc = DirichletBC(V, u_exact, "on_boundary")
# Set-up the variational problem and solve!
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
L = f*v*dx
u_h = Function(V)
solve(a == L, u_h, bc)
# Project the exact solution into V
P_u = project(u_exact, V)
return (u_h, P_u, f)
Ok, now we need a mesh. Either we have one lying about, or we could create one using CSG (constructive solid geometry) combined with a mesh generator. This is what the separate package "mshr" offers:
from mshr import *
big_box = Rectangle(Point(-1.0, -1.0), Point(1.0, 1.0))
small_box = Rectangle(Point(-1.0, -1.0), Point(0.0, 0.0))
mesh = generate_mesh(big_box - small_box, 1)
plot_inline(mesh)
Ok, let's try solving our Poisson problem on this mesh:
mesh = refine(mesh)
(u_h, P_u, f) = poisson(mesh)
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
plot_inline(u_h)
plot_inline(u_h.function_space().mesh())
That looks ok in the eyenorm, how about in some other norms? Let's try computing relevant Sobolev norms while refining the mesh a couple of times.
e_h = []
f_h = []
h = []
for i in range(4):
print "i = ", i
mesh = refine(mesh)
(u_h, P_h, f) = poisson(mesh)
e_h.append(errornorm(P_h, u_h, "L2"))
f_h.append(errornorm(P_h, u_h, "H10"))
h.append(mesh.hmax())
print "mesh sizes = ", h
print "L^2 errors = ", e_h
print "H^1_0 errors = ", f_h
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
i = 0 i = 1 i =
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
2 i =
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
3 mesh sizes = [0.35355339059327634, 0.17677669529663997, 0.0883883476483292, 0.04419417382417978] L^2 errors = [0.003505321993714835, 0.001383262553617792, 0.0005383586308775164, 0.00020958136026271565] H^1_0 errors = [0.03417998469550842, 0.021704653439925496, 0.013702196236018609, 0.008638029167554793]
We could plot the errors too, look here:
import pylab
pylab.loglog(h, e_h, '*-', label="$L^2$")
pylab.loglog(h, f_h, 'o-', label="$H^1_0$")
pylab.grid(True)
pylab.xlabel('h')
pylab.ylabel('E')
pylab.legend()
pylab.show()
Ok, those errors seem to decrease, but how quickly? Let's compute the convergence rates.
# Let's compute some convergence rates
import math
def convergence_rates(hs, errors):
rates = [(math.log(errors[i+1]/errors[i]))/(math.log(hs[i+1]/hs[i]))
for i in range(len(hs)-1)]
return rates
print "L2 rates = ", convergence_rates(h, e_h)
print "H^1_0 rates = ", convergence_rates(h, f_h)
L2 rates = [1.3414719573346296, 1.3614355586057703, 1.3610571408758334] H^1_0 rates = [0.6551473638919928, 0.6635972352031397, 0.6656330588015127]
Now, let's continue by trying the adaptive approach for the same problem. Let's start with our original mesh:
mesh = generate_mesh(big_box - small_box, 1)
mesh = refine(mesh)
plot_inline(mesh)
(u_h, P_h, f) = poisson(mesh)
plot_inline(u_h)
plot_inline(mesh)
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
Define the interior residual and the jump residual as follows: \begin{align} \quad r &= f + \mathrm{div} \, A \, \mathrm{grad} \, u_h \\ \quad j &= \mathrm{jump}(A \, \mathrm{grad} u_h, n) \end{align}
Now, we would like to compute our favorite error indicator, for instance this residual based one:
\begin{equation} E_T^2 = h_T^2 || r ||_{L^2(T)}^2 + h_T || j ||^2_{L^2(\partial T \backslash \partial \Omega)} \end{equation}# Get the cell sizes:
h = CellSize(mesh)
# Define the element residual:
r = f + div(grad(u_h))
# Define the jump residual
n = FacetNormal(mesh)
j = jump(grad(u_h), n)
Combining these with a set of piecewise constant basis functions that each are 1 on one cell and zero on the rest, we can define the error indicator as one variational form as follows. (This can be viewed as a neat trick.)
# Create a DG0 function space (piecewise constants, no continuity)
DG0 = FunctionSpace(mesh, "DG", 0)
# Create a test function on this space
w = TestFunction(DG0)
# Define the error indicator form
E_T = h**2*r**2*w*dx + 2*avg(w)*avg(h)*j*j*dS
# Define a function to hold the result of the evaluation
indicators = Function(DG0)
# Assemble the error indicator form into the coefficient vector
assemble(E_T, tensor=indicators.vector())
# Look at them!
plot_inline(indicators)
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
Let's combine this into an 'estimate' function:
def estimate(u_h, f):
# Get the cell sizes:
mesh = u_h.function_space().mesh()
h = CellSize(mesh)
# Define the element and jump residual:
r = f + div(grad(u_h))
n = FacetNormal(mesh)
j = jump(grad(u_h), n)
# Define the error indicator form
DG0 = FunctionSpace(mesh, "DG", 0)
w = TestFunction(DG0)
E_T = h**2*r*r*w*dx #+ 2*avg(w)*avg(h)*j*j*dS
# Evaluate the form
indicators = Function(DG0)
assemble(E_T, tensor=indicators.vector())
return indicators
So, now we can do this and still get exactly the same result:
mesh = generate_mesh(big_box - small_box, 1)
mesh = refine(mesh)
(u_h, P_h, f) = poisson(mesh)
indicators = estimate(u_h, f)
plot_inline(indicators)
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Level 25:FFC:Calling FFC just-in-time (JIT) compiler, this may take some time.
INFO:FFC:Compiling form ffc_form_cc24019d78b701a0baf517e241b1d5e1542dc95e
INFO:FFC:Compiler stage 1: Analyzing form(s)
INFO:FFC:-----------------------------------
INFO:FFC:
INFO:FFC: Geometric dimension: 2
Number of cell subdomains: 0
Rank: 1
Arguments: '(v_0)'
Number of coefficients: 1
Coefficients: '[f_4810]'
Unique elements: 'DG0(?), CG1(?)'
Unique sub elements: 'DG0(?), CG1(?)'
INFO:FFC: Extracting monomial form representation from UFL form
DEBUG:FFC: Monomial extraction failed: No handler defined for expression Division.
DEBUG:FFC: Estimated cost of tensor representation: -1
INFO:FFC: representation: auto --> quadrature
INFO:FFC: quadrature_degree: 3
INFO:FFC: quadrature_rule: auto --> default
INFO:FFC:
INFO:FFC:Compiler stage 1 finished in 0.037158 seconds.
INFO:FFC:Compiler stage 2: Computing intermediate representation
INFO:FFC:-------------------------------------------------------
INFO:FFC: Computing representation of 2 elements
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
INFO:FFC: Computing representation of 2 dofmaps
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
INFO:FFC: Computing representation of integrals
INFO:FFC: Computing quadrature representation
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC:
QG-utils, psi_tables:
{6: {FiniteElement('Lagrange', Domain(Cell('triangle', 2), label='dolfin_mesh_with_id_4794', data='<data with id 4794>'), 1, quad_scheme=None): {None: {None: {(0, 1): array([[ -1.00000000e+00, -1.00000000e+00, -1.00000000e+00,
-1.00000000e+00, -1.00000000e+00, -1.00000000e+00],
[ 8.88178420e-16, 8.88178420e-16, 8.88178420e-16,
8.88178420e-16, 8.88178420e-16, 8.88178420e-16],
[ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
1.00000000e+00, 1.00000000e+00, 1.00000000e+00]]), (0, 0): array([[ 0.10903901, 0.23193337, 0.10903901, 0.65902762, 0.23193337,
0.65902762],
[ 0.65902762, 0.65902762, 0.23193337, 0.23193337, 0.10903901,
0.10903901],
[ 0.23193337, 0.10903901, 0.65902762, 0.10903901, 0.65902762,
0.23193337]]), (0, 2): array([[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.]]), (2, 0): array([[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.]]), (1, 0): array([[-1., -1., -1., -1., -1., -1.],
[ 1., 1., 1., 1., 1., 1.],
[ 0., 0., 0., 0., 0., 0.]]), (1, 1): array([[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.]])}}}, VectorElement('Lagrange', Domain(Cell('triangle', 2), label='dolfin_mesh_with_id_4794', data='<data with id 4794>'), 1, dim=2, quad_scheme=None): {None: {None: {(0, 1): array([[[ -1.00000000e+00, -1.00000000e+00, -1.00000000e+00,
-1.00000000e+00, -1.00000000e+00, -1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],
[[ 8.88178420e-16, 8.88178420e-16, 8.88178420e-16,
8.88178420e-16, 8.88178420e-16, 8.88178420e-16],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],
[[ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
1.00000000e+00, 1.00000000e+00, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]],
[[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, -1.00000000e+00, -1.00000000e+00,
-1.00000000e+00, -1.00000000e+00, -1.00000000e+00]],
[[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 8.88178420e-16, 8.88178420e-16, 8.88178420e-16,
8.88178420e-16, 8.88178420e-16, 8.88178420e-16]],
[[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
1.00000000e+00, 1.00000000e+00, 1.00000000e+00]]]), (1, 0): array([[[-1., -1., -1., -1., -1., -1.],
[ 0., 0., 0., 0., 0., 0.]],
[[ 1., 1., 1., 1., 1., 1.],
[ 0., 0., 0., 0., 0., 0.]],
[[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.]],
[[ 0., 0., 0., 0., 0., 0.],
[-1., -1., -1., -1., -1., -1.]],
[[ 0., 0., 0., 0., 0., 0.],
[ 1., 1., 1., 1., 1., 1.]],
[[ 0., 0., 0., 0., 0., 0.],
[ 0., 0., 0., 0., 0., 0.]]]), (0, 0): array([[[ 0.10903901, 0.23193337, 0.10903901, 0.65902762, 0.23193337,
0.65902762],
[ 0. , 0. , 0. , 0. , 0. ,
0. ]],
[[ 0.65902762, 0.65902762, 0.23193337, 0.23193337, 0.10903901,
0.10903901],
[ 0. , 0. , 0. , 0. , 0. ,
0. ]],
[[ 0.23193337, 0.10903901, 0.65902762, 0.10903901, 0.65902762,
0.23193337],
[ 0. , 0. , 0. , 0. , 0. ,
0. ]],
[[ 0. , 0. , 0. , 0. , 0. ,
0. ],
[ 0.10903901, 0.23193337, 0.10903901, 0.65902762, 0.23193337,
0.65902762]],
[[ 0. , 0. , 0. , 0. , 0. ,
0. ],
[ 0.65902762, 0.65902762, 0.23193337, 0.23193337, 0.10903901,
0.10903901]],
[[ 0. , 0. , 0. , 0. , 0. ,
0. ],
[ 0.23193337, 0.10903901, 0.65902762, 0.10903901, 0.65902762,
0.23193337]]])}}}, FiniteElement('Discontinuous Lagrange', Domain(Cell('triangle', 2), label='dolfin_mesh_with_id_4794', data='<data with id 4794>'), 0, quad_scheme=None): {None: {None: {(0, 0): array([[ 1., 1., 1., 1., 1., 1.]])}}}}}
DEBUG:FFC:
QG-utils, psi_tables, flat_tables:
{'FE0_D10': array([[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.]]), 'FE0_D02': array([[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]]), 'FE0_D01': array([[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00]]), 'FE0_D11': array([[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]]), 'FE2_C0_D01': array([[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]), 'FE2_C0_D10': array([[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.]]), 'FE0_D20': array([[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]]), 'FE2_C1': array([[ 0. , 0. , 0. , 0.10903901, 0.65902762,
0.23193337],
[ 0. , 0. , 0. , 0.23193337, 0.65902762,
0.10903901],
[ 0. , 0. , 0. , 0.10903901, 0.23193337,
0.65902762],
[ 0. , 0. , 0. , 0.65902762, 0.23193337,
0.10903901],
[ 0. , 0. , 0. , 0.23193337, 0.10903901,
0.65902762],
[ 0. , 0. , 0. , 0.65902762, 0.10903901,
0.23193337]]), 'FE2_C0': array([[ 0.10903901, 0.65902762, 0.23193337, 0. , 0. ,
0. ],
[ 0.23193337, 0.65902762, 0.10903901, 0. , 0. ,
0. ],
[ 0.10903901, 0.23193337, 0.65902762, 0. , 0. ,
0. ],
[ 0.65902762, 0.23193337, 0.10903901, 0. , 0. ,
0. ],
[ 0.23193337, 0.10903901, 0.65902762, 0. , 0. ,
0. ],
[ 0.65902762, 0.10903901, 0.23193337, 0. , 0. ,
0. ]]), 'FE2_C1_D01': array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00]]), 'FE0': array([[ 0.10903901, 0.65902762, 0.23193337],
[ 0.23193337, 0.65902762, 0.10903901],
[ 0.10903901, 0.23193337, 0.65902762],
[ 0.65902762, 0.23193337, 0.10903901],
[ 0.23193337, 0.10903901, 0.65902762],
[ 0.65902762, 0.10903901, 0.23193337]]), 'FE1': array([[ 1.],
[ 1.],
[ 1.],
[ 1.],
[ 1.],
[ 1.]]), 'FE2_C1_D10': array([[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.]])}
DEBUG:FFC:
tables: {'FE0_D10': array([[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.]]), 'FE0_D02': array([[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]]), 'FE0_D01': array([[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00]]), 'FE2_C0_D01': array([[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ -1.00000000e+00, 8.88178420e-16, 1.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]), 'FE2_C0_D10': array([[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.]]), 'FE2_C1': array([[ 0. , 0. , 0. , 0.10903901, 0.65902762,
0.23193337],
[ 0. , 0. , 0. , 0.23193337, 0.65902762,
0.10903901],
[ 0. , 0. , 0. , 0.10903901, 0.23193337,
0.65902762],
[ 0. , 0. , 0. , 0.65902762, 0.23193337,
0.10903901],
[ 0. , 0. , 0. , 0.23193337, 0.10903901,
0.65902762],
[ 0. , 0. , 0. , 0.65902762, 0.10903901,
0.23193337]]), 'FE2_C0': array([[ 0.10903901, 0.65902762, 0.23193337, 0. , 0. ,
0. ],
[ 0.23193337, 0.65902762, 0.10903901, 0. , 0. ,
0. ],
[ 0.10903901, 0.23193337, 0.65902762, 0. , 0. ,
0. ],
[ 0.65902762, 0.23193337, 0.10903901, 0. , 0. ,
0. ],
[ 0.23193337, 0.10903901, 0.65902762, 0. , 0. ,
0. ],
[ 0.65902762, 0.10903901, 0.23193337, 0. , 0. ,
0. ]]), 'FE2_C1_D01': array([[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
-1.00000000e+00, 8.88178420e-16, 1.00000000e+00]]), 'FE0': array([[ 0.10903901, 0.65902762, 0.23193337],
[ 0.23193337, 0.65902762, 0.10903901],
[ 0.10903901, 0.23193337, 0.65902762],
[ 0.65902762, 0.23193337, 0.10903901],
[ 0.23193337, 0.10903901, 0.65902762],
[ 0.65902762, 0.10903901, 0.23193337]]), 'FE1': array([[ 1.],
[ 1.],
[ 1.],
[ 1.],
[ 1.],
[ 1.]]), 'FE2_C1_D10': array([[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.]])}
DEBUG:FFC:
name_map: {'FE0_D02': ['FE0_D11', 'FE0_D20']}
DEBUG:FFC:
inv_name_map: {'FE0_D10': 'FE0_D10', 'FE0_D11': 'FE0_D02', 'FE0_D01': 'FE0_D01', 'FE0_D02': 'FE0_D02', 'FE2_C0': 'FE2_C0', 'FE1': 'FE1', 'FE0_D20': 'FE0_D02', 'FE2_C1': 'FE2_C1', 'FE2_C0_D01': 'FE2_C0_D01', 'FE2_C1_D01': 'FE2_C1_D01', 'FE0': 'FE0', 'FE2_C0_D10': 'FE2_C0_D10', 'FE2_C1_D10': 'FE2_C1_D10'}
DEBUG:FFC:
QG-utils, psi_tables, unique_tables:
{'FE0_D10': array([[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.],
[-1., 1., 0.]]), 'FE0_D02': array([[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]]), 'FE0_D01': array([[-1., 0., 1.],
[-1., 0., 1.],
[-1., 0., 1.],
[-1., 0., 1.],
[-1., 0., 1.],
[-1., 0., 1.]]), 'FE2_C0_D01': array([[-1., 0., 1., 0., 0., 0.],
[-1., 0., 1., 0., 0., 0.],
[-1., 0., 1., 0., 0., 0.],
[-1., 0., 1., 0., 0., 0.],
[-1., 0., 1., 0., 0., 0.],
[-1., 0., 1., 0., 0., 0.]]), 'FE2_C0_D10': array([[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.],
[-1., 1., 0., 0., 0., 0.]]), 'FE2_C1': array([[ 0. , 0. , 0. , 0.10903901, 0.65902762,
0.23193337],
[ 0. , 0. , 0. , 0.23193337, 0.65902762,
0.10903901],
[ 0. , 0. , 0. , 0.10903901, 0.23193337,
0.65902762],
[ 0. , 0. , 0. , 0.65902762, 0.23193337,
0.10903901],
[ 0. , 0. , 0. , 0.23193337, 0.10903901,
0.65902762],
[ 0. , 0. , 0. , 0.65902762, 0.10903901,
0.23193337]]), 'FE2_C0': array([[ 0.10903901, 0.65902762, 0.23193337, 0. , 0. ,
0. ],
[ 0.23193337, 0.65902762, 0.10903901, 0. , 0. ,
0. ],
[ 0.10903901, 0.23193337, 0.65902762, 0. , 0. ,
0. ],
[ 0.65902762, 0.23193337, 0.10903901, 0. , 0. ,
0. ],
[ 0.23193337, 0.10903901, 0.65902762, 0. , 0. ,
0. ],
[ 0.65902762, 0.10903901, 0.23193337, 0. , 0. ,
0. ]]), 'FE2_C1_D01': array([[ 0., 0., 0., -1., 0., 1.],
[ 0., 0., 0., -1., 0., 1.],
[ 0., 0., 0., -1., 0., 1.],
[ 0., 0., 0., -1., 0., 1.],
[ 0., 0., 0., -1., 0., 1.],
[ 0., 0., 0., -1., 0., 1.]]), 'FE0': array([[ 0.10903901, 0.65902762, 0.23193337],
[ 0.23193337, 0.65902762, 0.10903901],
[ 0.10903901, 0.23193337, 0.65902762],
[ 0.65902762, 0.23193337, 0.10903901],
[ 0.23193337, 0.10903901, 0.65902762],
[ 0.65902762, 0.10903901, 0.23193337]]), 'FE1': array([[ 1.],
[ 1.],
[ 1.],
[ 1.],
[ 1.],
[ 1.]]), 'FE2_C1_D10': array([[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.],
[ 0., 0., 0., -1., 1., 0.]])}
DEBUG:FFC:
QG-utils, psi_tables, name_map:
{'FE0_D10': ('FE0_D10', (), False, False), 'FE0_D11': ('FE0_D02', (), True, False), 'FE0_D01': ('FE0_D01', (), False, False), 'FE0_D02': ('FE0_D02', (), True, False), 'FE2_C0': ('FE2_C0', (), False, False), 'FE1': ('FE1', (), False, True), 'FE0_D20': ('FE0_D02', (), True, False), 'FE2_C1': ('FE2_C1', (), False, False), 'FE2_C0_D01': ('FE2_C0_D01', (), False, False), 'FE2_C1_D01': ('FE2_C1_D01', (), False, False), 'FE0': ('FE0', (), False, False), 'FE2_C0_D10': ('FE2_C0_D10', (), False, False), 'FE2_C1_D10': ('FE2_C1_D10', (), False, False)}
INFO:FFC: Transforming cell integral
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
DEBUG:FFC: Reusing element from cache
INFO:FFC: Computing representation of forms
INFO:FFC:
INFO:FFC:Compiler stage 2 finished in 0.0937929 seconds.
INFO:FFC:Compiler stage 3: Optimizing intermediate representation
INFO:FFC:--------------------------------------------------------
INFO:FFC: Skipping optimizations, add -O to optimize
INFO:FFC:
INFO:FFC:Compiler stage 3 finished in 0.00133109 seconds.
INFO:FFC:Compiler stage 4: Generating code
INFO:FFC:---------------------------------
INFO:FFC: Generating code for 2 element(s)
DEBUG:FFC: Removing unused variable: tmp2
DEBUG:FFC: Removing unused variable: tmp1
DEBUG:FFC: Removing unused variable: tmp7
DEBUG:FFC: Removing unused variable: tmp6
DEBUG:FFC: Removing unused variable: tmp5
DEBUG:FFC: Removing unused variable: tt
DEBUG:FFC: Removing unused variable: ss
DEBUG:FFC: Removing unused variable: rr
DEBUG:FFC: Removing unused variable: tmp2
DEBUG:FFC: Removing unused variable: tmp1
DEBUG:FFC: Removing unused variable: tmp7
DEBUG:FFC: Removing unused variable: tmp6
DEBUG:FFC: Removing unused variable: tmp5
DEBUG:FFC: Removing unused variable: tt
DEBUG:FFC: Removing unused variable: ss
DEBUG:FFC: Removing unused variable: rr
DEBUG:FFC: Removing unused variable: tmp2
DEBUG:FFC: Removing unused variable: tmp1
DEBUG:FFC: Removing unused variable: tmp7
DEBUG:FFC: Removing unused variable: tmp6
DEBUG:FFC: Removing unused variable: tmp5
DEBUG:FFC: Removing unused variable: tt
DEBUG:FFC: Removing unused variable: ss
DEBUG:FFC: Removing unused variable: rr
DEBUG:FFC: Removing unused variable: tmp2
DEBUG:FFC: Removing unused variable: tmp1
DEBUG:FFC: Removing unused variable: tmp7
DEBUG:FFC: Removing unused variable: tmp6
DEBUG:FFC: Removing unused variable: tmp5
DEBUG:FFC: Removing unused variable: tt
DEBUG:FFC: Removing unused variable: ss
DEBUG:FFC: Removing unused variable: rr
DEBUG:FFC: Removing unused variable: tmp2
DEBUG:FFC: Removing unused variable: tmp1
DEBUG:FFC: Removing unused variable: tmp7
DEBUG:FFC: Removing unused variable: tmp6
DEBUG:FFC: Removing unused variable: tmp5
DEBUG:FFC: Removing unused variable: tt
DEBUG:FFC: Removing unused variable: ss
DEBUG:FFC: Removing unused variable: rr
DEBUG:FFC: Removing unused variable: tmp2
DEBUG:FFC: Removing unused variable: tmp1
DEBUG:FFC: Removing unused variable: tmp7
DEBUG:FFC: Removing unused variable: tmp6
DEBUG:FFC: Removing unused variable: tmp5
DEBUG:FFC: Removing unused variable: tt
DEBUG:FFC: Removing unused variable: ss
DEBUG:FFC: Removing unused variable: rr
DEBUG:FFC: Removing unused variable: tmp2
DEBUG:FFC: Removing unused variable: tmp1
DEBUG:FFC: Removing unused variable: tmp0
DEBUG:FFC: Removing unused variable: tmp7
DEBUG:FFC: Removing unused variable: tmp6
DEBUG:FFC: Removing unused variable: tmp5
DEBUG:FFC: Removing unused variable: tt
DEBUG:FFC: Removing unused variable: ss
DEBUG:FFC: Removing unused variable: rr
DEBUG:FFC: Removing unused variable: Y
DEBUG:FFC: Removing unused variable: X
DEBUG:FFC: Removing unused variable: C1
DEBUG:FFC: Removing unused variable: C0
INFO:FFC: Generating code for 2 dofmap(s)
INFO:FFC: Generating code for integrals
INFO:FFC: Generating code for forms
INFO:FFC:
INFO:FFC:Compiler stage 4 finished in 0.057447 seconds.
INFO:FFC:Compiler stage 4.1 finished in 1.90735e-06 seconds.
INFO:FFC:Compiler stage 5: Formatting code
INFO:FFC:---------------------------------
INFO:FFC: Output written to ./ffc_form_cc24019d78b701a0baf517e241b1d5e1542dc95e.h.
INFO:FFC:
INFO:FFC:Compiler stage 5 finished in 0.00183201 seconds.
INFO:FFC:FFC finished in 0.193662 seconds.
DEBUG:FFC:Compiling and linking Python extension module, this may take some time.
The next step is to use these error estimators for a marking strategy. Dorfler is a good choice:
import numpy
def mark(alpha, indicators):
# Sort eta_T in decreasing order and keep track of the cell numbers
etas = indicators.vector().array()
indices = etas.argsort()[::-1]
sorted = etas[indices]
# Compute sum and fraction of indicators
total = sum(sorted)
fraction = alpha*total
# Define cell function to hold markers
mesh = indicators.function_space().mesh()
markers = CellFunction("bool", mesh, False)
# Iterate over the cells
v = 0.0
for i in indices:
# Stop if we have marked enough
if v >= fraction:
break
# Otherwise
markers[i] = True
v += sorted[i]
return markers
# Let's look at these:
markers = mark(0.3, indicators)
print markers.array()
[False True True False False False True False False True False False True False True True False True True True False True False True]
Ok, last thing is to actually refine the mesh, that is easy:
mesh = refine(mesh, markers)
plot_inline(mesh)
import time
parameters["refinement_algorithm"] = "regular_cut"
#parameters["refinement_algorithm"] = "recursive_bisection"
mesh = generate_mesh(big_box - small_box, 1)
mesh = refine(mesh)
tolerance = 1.e-3
max_iterations = 10
alpha = 0.1
errors = []
etas = []
Ns = []
i = 1
while (True):
# Solve
(u_h, P_h, f) = poisson(mesh)
E_h = errornorm(u_h, P_h, "H10")
errors.append(E_h)
# Estimate
indicators = estimate(u_h, f)
# Check stopping criterion
eta = sum(indicators.vector().array())
etas.append(eta)
Ns.append(mesh.num_cells())
if eta < tolerance or i >= max_iterations:
break
# Mark
markers = mark(alpha, indicators)
# Refine
mesh = refine(mesh, markers)
#
i += 1
DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache. DEBUG:FFC:Reusing form from cache.
plot_inline(mesh)
plot_inline(u_h)
print "Ns = ", Ns
print "etas = ", etas
pylab.loglog(Ns, etas, '*-', label="Estimate")
pylab.loglog(Ns, errors, 'o-', label="Actual")
pylab.legend()
pylab.grid(True)
pylab.show()
Ns = [24, 64, 104, 152, 266, 376, 490, 769, 1187, 1725] etas = [1.8170118921801477, 0.9645955851256065, 0.59235414631234007, 0.38542080313965621, 0.22203444405926043, 0.16265809387480953, 0.12420782606645292, 0.081372939392402394, 0.051849152384037048, 0.036444731215497982]