...

On Wed, Jun 13, 2012 at 02:50:56PM -0000, Yahya wrote:
> New question #200319 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/200319
>
> from dolfin import*
>
> mesh=UnitSquare(4,4)
>
> U=FunctionSpace(mesh, "Lagrange", 1)
> V=FunctionSpace(mesh, "Lagrange", 1)
> W=FunctionSpace(mesh, "Lagrange", 1)
>
> M=MixedFunctionSpace([U,V,W])
> boundary_parts = MeshFunction("uint", mesh, 1)
>
> tol = 1E-14
> class UpperBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and abs(x[1]-1) < tol
>
> Gamma_o = UpperBoundary()
> Gamma_o.mark(boundary_parts, 0)
>
> class LowerBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and abs(x[1]) < tol
>
> Gamma_c = LowerBoundary()
> Gamma_c.mark(boundary_parts, 1)
>
> class RestBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and abs(x[0]) < tol and abs(x[1]-1) < tol and abs(x[0]-1) < tol
>
> Gamma_t_c = RestBoundary()
> Gamma_t_c.mark(boundary_parts, 2)
>
> (u,q, lam)=TrialFunctions(M)
> (v_u,v_q, v_lam)=TestFunctions(M)
>
> u0=Constant("1.0")
> u_0=interpolate(u0, U)
> alpha=(1.0/1000)
>
> a=inner(grad(lam),grad(v_u))*dx + inner(grad(u),grad(v_lam))*dx+\
> inner(lam,v_u)*dx +inner(u,v_lam)*dx-q*v_lam*dx+\
> alpha*(q*v_q)*dx+(lam*v_q)*dx -u*v_u*ds(0)
> L= -u0*v_u*ds(0)
>
> #Question. HOW I CAN assemble THE VARIATIONAL EQUATION. HOW I CAN DEFINE solve(a==L,solution) IN THIS CASE.
> solution=Function(M)
> solve(a==L,solution)
> u,q,lam=solution.split()
>
> from math import cosh,sinh
>
> Ue=FunctionSpace(mesh, "Lagrange", 1)
> u_e=Expression("(u0/(v1-alph))*(((1+(cosh(1)/sinh(1)))*(-cosh(x[1])/(2*sinh(1))))+(x[1]*sinh(x[1])/(2*sinh(1))))",u0=1.0,v1=(((((1+((exp(1)+exp(-1))/(exp(1)-exp(-1))))*((-exp(1)-exp(-1))/(2*exp(1)-2*exp(-1)))))+(1*(exp(1)-exp(-1))/(2*exp(1)-2*exp(-1))))),alph=(1.0/1000))
> u_ext=interpolate(u_e,Ue)
> print'u:',u_ext.vector().array()
>
> Ve=FunctionSpace(mesh, "Lagrange", 1)
> q_e=Expression("(-u0/(v1-alph))*(cosh(x[1])/sinh(1))",u0=1.0,v1=(((1+(cosh(1)/sinh(1)))*(-cosh(1)/(2*sinh(1))))+(1*sinh(1)/(2*sinh(1)))),alph=(1.0/1000))
> q_ext=interpolate(q_e,Ve)
> print'q:',q_ext.vector().array()
> print'q.max:',max( q_ext.vector().array())
> print'q.min:',min( q_ext.vector().array())
>
> We=FunctionSpace(mesh, "Lagrange", 1)
> lam_e=Expression("(alph*u0/(v1-alph))*cosh(x[1])/sinh(1)",u0=1.0,v1=(((1+(cosh(1)/sinh(1)))*(-cosh(1)/(2*sinh(1))))+(1*sinh(1)/(2*sinh(1)))),alph=(1.0/1000))
> lam_ext=interpolate(lam_e,We)
> print'lam:',lam_ext.vector().array()
>
> #Question.HOW CAN SEE (u,q,lam) VALUE WITHOUT USING range() . I USED u.vector().array(), q.vector().array(), lam.vecto().arra() #DOES NOT WORK
> print 'u.e:', solution.vector().array()[range(0,25)]
> print 'q.e:', solution.vector().array()[range(25,50)]
> print 'lam.e:', solution.vector().array()[range(50,75)]
>
> #Question.ERROR IS NOT DECREASING WITH FACTOR 1/4 AS MESH(4,4) CHANGED MESH(8,8), MESH(16,16)
> error_u = u*u*dx + u_ext*u_ext*dx-2*u*u_ext*dx
> E_u = sqrt(assemble(error_u))
> print '|| u - u.e ||_2 : ', E_u
>
> error_q = q*q*dx + q_ext*q_ext*dx-2*q*q_ext*dx
> E_q = sqrt(assemble(error_q))
> print '|| q - q.e ||_2 : ', E_q
>
> error_lam = lam*lam*dx + lam_ext*lam_ext*dx-2*lam*lam_ext*dx
> E_lam = sqrt(assemble(error_lam))
> print '|| lam - lam.e ||_2 : ', E_lam
>
> #Question. HOW I CAN assemble THE VARIATIONAL EQUATION. HOW I CAN DEFINE solve(a==L,solution) IN THIS CASE.
> #Question.HOW CAN SEE (u,q,lam) VALUE WITHOUT USING range() . I USED u.vector().array(), q.vector().array(), lam.vecto().arra() DOES NOT WORK
> #Question.ERROR IS NOT DECREASING WITH FACTOR 1/4 AS MESH(4,4) CHANGED MESH(8,8), MESH(16,16)

You need to post a more specific question. This is incomprehensible.

--
Anders

This problem is too complex for us to debug. It would take me a long
time to check that your analytic solution is correct and that you are
actually solving the correct problem.

One hint I can give you is to use the errornorm() function. Try

  help(errornorm)

That will at least simplify the computation of the error.

Also try starting with a simpler solution. It seems unnecessarily
complex.

--
Anders

On Sun, Jun 17, 2012 at 08:06:13PM -0000, Yahya wrote:
> Question #200319 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/200319
>
> Description changed to:
> #Question.ERROR IS NOT DECREASING WITH FACTOR 1/4 AS MESH(4,4) CHANGED
> MESH(8,8), MESH(16,16) Can you tell me reason
>
> from dolfin import*
>
> mesh=UnitSquare(4,4)
>
> U=FunctionSpace(mesh, "Lagrange", 1)
> V=FunctionSpace(mesh, "Lagrange", 1)
> W=FunctionSpace(mesh, "Lagrange", 1)
>
> M=MixedFunctionSpace([U,V,W])
> boundary_parts = MeshFunction("uint", mesh, 1)
>
> tol = 1E-14
> class UpperBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return on_boundary and abs(x[1]-1) < tol
>
> Gamma_o = UpperBoundary()
> Gamma_o.mark(boundary_parts, 0)
> (u,q, lam)=TrialFunctions(M)
> (v_u,v_q, v_lam)=TestFunctions(M)
>
> u0=Constant("1.0")
> u_0=interpolate(u0, U)
> alpha=(1.0/1000)
>
> a=inner(grad(lam),grad(v_u))*dx + inner(grad(u),grad(v_lam))*dx+\
> inner(lam,v_u)*dx +inner(u,v_lam)*dx-q*v_lam*dx+\
> alpha*(q*v_q)*dx+(lam*v_q)*dx -u*v_u*ds(0)
> L= -u0*v_u*ds(0)
>
> solution=Function(M)
> solve(a==L,solution)
> u,q,lam=solution.split()
>
> from math import cosh,sinh
>
> Ue=FunctionSpace(mesh, "Lagrange", 1)
> u_e=Expression("(u0/(v1-alph))*(((1+(cosh(1)/sinh(1)))*(-cosh(x[1])/(2*sinh(1))))+(x[1]*sinh(x[1])/(2*sinh(1))))",u0=1.0,v1=(((((1+((exp(1)+exp(-1))/(exp(1)-exp(-1))))*((-exp(1)-exp(-1))/(2*exp(1)-2*exp(-1)))))+(1*(exp(1)-exp(-1))/(2*exp(1)-2*exp(-1))))),alph=(1.0/1000))
> u_ext=interpolate(u_e,Ue)
> print'u:',u_ext.vector().array()
>
> Ve=FunctionSpace(mesh, "Lagrange", 1)
> q_e=Expression("(-u0/(v1-alph))*(cosh(x[1])/sinh(1))",u0=1.0,v1=(((1+(cosh(1)/sinh(1)))*(-cosh(1)/(2*sinh(1))))+(1*sinh(1)/(2*sinh(1)))),alph=(1.0/1000))
> q_ext=interpolate(q_e,Ve)
> print'q:',q_ext.vector().array()
>
> We=FunctionSpace(mesh, "Lagrange", 1)
> lam_e=Expression("(alph*u0/(v1-alph))*cosh(x[1])/sinh(1)",u0=1.0,v1=(((1+(cosh(1)/sinh(1)))*(-cosh(1)/(2*sinh(1))))+(1*sinh(1)/(2*sinh(1)))),alph=(1.0/1000))
>
>
> error_u = u*u*dx + u_ext*u_ext*dx-2*u*u_ext*dx
> E_u = sqrt(assemble(error_u))
> print '|| u - u.e ||_2 : ', E_u
>
> error_q = q*q*dx + q_ext*q_ext*dx-2*q*q_ext*dx
> E_q = sqrt(assemble(error_q))
> print '|| q - q.e ||_2 : ', E_q
>
> error_lam = lam*lam*dx + lam_ext*lam_ext*dx-2*lam*lam_ext*dx
> E_lam = sqrt(assemble(error_lam))
> print '|| lam - lam.e ||_2 : ', E_lam
>