Getting error on assemble

Asked by K.Manickam on 2013-09-12

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

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+phi*p*dx-phi*y*ds(1)+inner(grad(y),grad(si))*dx+y*si*dx-(1/d)*p*si*ds(0)
L = f*si*dx-y0*phi*ds(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 - d*sinh(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('((d*sinh(1)*cosh(x[1]))/(1-d*sinh(1)*sinh(1)))*(y_0 + (1.0/sinh(1)) - .5)',y_0 = 1,d = 1)
u_exact = Expression('-sinh(1)/(1-p*sinh(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_0*phi*ds(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()

Some one help to rectify error.
Thanks
Manickam.K

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
2013-09-12
Last reply:
2013-09-23
Johannes Ring (johannr) said : #1

FEniCS no longer uses Launchpad for Questions & Answers. Please consult the documentation on the FEniCS web page for where and how to (re)post your question: http://fenicsproject.org/support/

Can you help with this problem?

Provide an answer of your own, or ask K.Manickam for more information if necessary.

To post a message you must log in.