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 言語で書かれた数式
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")
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()
誤差評価式:
$$ \| \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\|$を計算します。
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