In Python, I am trying to solve finite element methods for optimal control problems on Neumann conditions. But assemble does not work properly getting some error in line writing likes this given below
error_y = (y - y_e)(y - y_e)dx
error = sqrt(assemble(error_y))
print"y_error:",error
Similiar, Same way getting error in line given below
error = sqrt(assemble(error_p))
error = sqrt(assemble(error_u))
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
Here Python Code for Finite element methods for optimal control problems for further clarification.
from dolfin import *
from math import *
import sys
Create mesh and define function space
nx = 32
ny = 32
y_0=0
mesh = UnitSquare (nx,ny)
X = FunctionSpace (mesh,'Lagrange',1)
Y = FunctionSpace (mesh,'Lagrange',1)
Z = X*Y
(y,p) = TrialFunctions(Z)
(si,phi) = TestFunctions(Z)
Define Sourse values
f = Expression('x[1]x[2]/2-x[1]')
y0 = Constant ('1')
d = Constant('1') # regularize parameter
a = inner(grad(phi),grad(p))dx+phipdx-phiyds(1)+inner(grad(y),grad(si))dx+ysidx-(1/d)psids(0)
L = fsidx-y0phids(1)
boundary
boundary_parts = MeshFunction ('uint', mesh, 1)
Mark lower boundary facets as subdomain 0
class LowerNeumannBoundary(SubDomain):
def inside(self,x,on_boundary):
tol = 1E-14
return on_boundary and abs(x[1]) < tol
tolerance for coordinate comparisons
L = LowerNeumannBoundary()
L.mark(boundary_parts, 0)
Mark upper boundary facets as subdomain 1
class UpperNeumannBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and abs(x[1] - 1) < tol
tolerance for coordinate comparisons
U = UpperNeumannBoundary()
U.mark(boundary_parts, 1)
all of the Rest boundaries
class RestNeumannBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and (abs(x[0]) < tol or abs (x[0] - 1) < tol)
tolerance for coordinate comparisons
Verification
u=0
y_exact = Expression('(cosh(x[1] - 1)/(1 - dsinh(1)sinh(1)))(y_0+1/sinh(1) - .5) - (cosh(x[1] - 1)/sinh(1)) + x[1] * x[1]/2 - x[1] + 1' ,y_0 = 1,d = 1)
p_exact = Expression('((dsinh(1)cosh(x[1]))/(1-dsinh(1)sinh(1)))(y_0 + (1.0/sinh(1)) - .5)',y_0 = 1,d = 1)
u_exact = Expression('-sinh(1)/(1-psinh(1)sinh(1))*(u0+(1/sinh(1))-0.5)',u0=1,p=1)
Xe = FunctionSpace(mesh, 'Lagrange', 5)
y_e = interpolate(y_exact,Xe)
Ye = FunctionSpace(mesh, 'Lagrange', 5)
p_e = interpolate (p_exact, Ye,)
ze = FunctionSpace(mesh, 'Lagrange', 5)
u_e = interpolate(u_exact ,ze)
error_y = (y - y_e)(y - y_e)dx
error = sqrt(assemble(error_y))
print"y_error:",error
error_p = (p - p_e)(p - p_e)dx
error = sqrt(assemble(error_p))
print"p_error:",error
error_u = (u - u_e)(u - u_e)dx
error = sqrt(assemble(error_u))
print"u_error:",error
R = RestNeumannBoundary()
R.mark(boundary_parts,23)
M= assemble(y_0phids(1),exterior_facet_domains=boundary_parts)
Compute solution
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
s = Function(Z)
solve(A,s.vector(),b)
(y,p) = s.split()
print'(y,p) :',s.vector().array()
print'n'
print'(y) :',y.vector().array()
print'n'
print'(p) :',p.vector().array()
print'(M) :',M.array()
B= (.5,.5)
print'y_e at the centor:',y_e(B)
print'y at the centor:',y(B)
print'p_e at the centor:',p_e(B)
print'p at the centor:',p(B)
cost=inner(y-y_0,y-y_0)ds(1)+(1.0/d)inner(p,p)ds(0)
J_h=(1.0/2)assemble(cost, exterior_facet_domains=boundary_parts)
coste=inner(y_e_Ve-y_0,y_e_Ve-y_0)ds(1)+(1.0/d)inner(p_e_Ve,p_e_Ve)ds(0)
J_ex=(1.0/2)assemble(coste, exterior_facet_domains=boundary_parts)
E5=abs(J_h-J_ex)
t_j=f_pv+s_pv
print'(y_m):',max (y.vector().array())
print'(y_em):',max (y_e_Ve.vector().array())
print'(y_n):',min (y.vector().array())
print'(y_en):',min (y_e_Ve.vector().array())
print"n"
print'f_pv:',f_pv
print's_pv:',s_pv
print'tot:',t_j
print'J_h(y,u):',J_h
print'J_e(y,u):',J_ex
print'(p_m):', max(p.vector().array())
print'(p_em):', max(p_e_Ve.vector().array())
print'(p_n):', min(p.vector().array())
print'(p_en):', min(p_e_Ve.vector().array())
Plot solution
plot(y,title="yplot")
plot(p,title="pplot")
plot(mesh,title="mesh")
interactive()
