Hypercircle法による誤差評価

準備

In [1]:
from dolfin import *

#------------ mesh triangulation and data function f  --------------#
N=16
mesh = UnitSquareMesh(N, N)
f = Expression("2*x[0]*(1- x[0]) + 2*x[1]*(1- x[1])", degree=4, domain=mesh)  # C 言語で書かれた数式

適合有限要素解の計算

In [2]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
%matplotlib inline


#---------------   boundary condition  ----------------------------#
V = FunctionSpace(mesh, "Lagrange", 1)

def boundary(x):
    return near(x[0], 0) or near(x[0],1.0) or near(x[1], 0) or near(x[1],1.0)

u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

#---------------   Compute solution   ------------------------------#
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx
b = f*v*dx

conforming_fem_u = Function(V)
solve(a == b, conforming_fem_u, bc)


#-------------  Error computation ------------------------------------#
exact_solution = Expression("x[0]*(1- x[0])*x[1]*(1-x[1])", degree=4)
mesh2 = refine(mesh)
V4 = FunctionSpace(mesh2, "DG", 4);
fem_error = project(conforming_fem_u - exact_solution, V4);
print "Error in L2 norm: %f"% norm(fem_error,"L2")
print "Error in H1 norm: %f"% norm(fem_error,"H1")

/usr/lib/python2.7/dist-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

Calling FFC just-in-time (JIT) compiler, this may take some time.
Error in L2 norm: 0.000366
Error in H1 norm: 0.015185

混合法の有限要素解の計算

In [3]:
degree=1

RT_Element = FiniteElement("RT", mesh.ufl_cell(), degree)
DG_Element = FiniteElement("DG", mesh.ufl_cell(), degree-1)


MixedV = FunctionSpace( mesh , MixedElement ([RT_Element, DG_Element]) )

#---------------   Variational equation  ----------------------------#
(p,u) = TrialFunctions(MixedV)
(q,v) = TestFunctions(MixedV)

a = (inner(p, q)+ div(p)*v + div(q)*u) *dx
b = -f*v*dx

#---------------   Compute solution   ------------------------------#
mixed_fem_sol = Function(MixedV)
solve(a == b, mixed_fem_sol)

(p,u) = mixed_fem_sol.split()

Hypercircle法による誤差評価

誤差評価式：

$$ \| \nabla u - \nabla u_h \| \le C_0(h) \|f-f_h\| + \|\nabla u_h - p_h\| $$

以下、$Err1:= C_0(h) \|f-f_h\|$ と$Err2:= \|\nabla u_h - p_h\|$を計算します。

In [4]:
DG = FunctionSpace(mesh, "DG", degree-1)

C0_h = 1/pi/N
fh = project(f, DG);

Err1 = C0_h * sqrt( assemble( (f-fh)*(f-fh)*dx ) )

v1 = grad(conforming_fem_u)-p
Err2 = sqrt( assemble(  inner(v1,v1)*dx(mesh) ) )

print "Error estimation by using Hypercricle method:",  Err1 + Err2

Calling FFC just-in-time (JIT) compiler, this may take some time.
Error estimation by using Hypercricle method: 0.0182513641172