dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

I cannot spot anything right now.

You might want to write down the math of your form and maybe someone else 
might give you a hint.

Johan

On Monday January 10 2011 10:36:53 B. Emek Abali wrote:
> Question #140910 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/140910
> 
> B. Emek Abali posted a new comment:
> it is an iterative nonlinear problem, by cutting all the leaves away,
> this code leads still to NaN solutions
> 
> from dolfin import *
> import numpy
> 
> def FormF(v,w,p,q):
> 	F = w[i]*v[i].dx(j)*v[j]*dx
> 	return F
> 
> def FormG(v,w,dv,p,q,dp):
> 	G = w[i]*dv[i].dx(j)*v[j]*dx
> 	return G
> 
> mesh = Box(0.0, -0.5, -0.5, 2.0, 0.5, 0.5, 10, 5, 5)
> V = VectorFunctionSpace(mesh, 'CG', 1)
> Q = FunctionSpace(mesh, 'CG', 1)
> ME = MixedFunctionSpace([V,Q])
> 
> subdomains = MeshFunction("uint", mesh, mesh.topology().dim()-1)
> right, left, back, front, top, bottom = compile_subdomains(["(fabs(x[0] -
> 2.0) < DOLFIN_EPS) && on_boundary", "(fabs(x[0] + 0.0) < DOLFIN_EPS) &&
> on_boundary", "(fabs(x[1] - 0.5) < DOLFIN_EPS) && on_boundary",
> "(fabs(x[1] + 0.5) < DOLFIN_EPS) && on_boundary", "(fabs(x[2] - 0.5) <
> DOLFIN_EPS) && on_boundary", "(fabs(x[2] + 0.5) < DOLFIN_EPS) &&
> on_boundary"])
> 
> null = Constant(0.0)
> class TopDrag(Expression):
> 	def eval(self, vec, x):
> 		vec[0] = 1.0
> 		vec[1] = 0.0
> 		vec[2] = 0.0
> 	def dim(self):
> 		return 3
> 
> class BottomHold(Expression):
> 	def eval(self, vec, x):
> 		vec[0] = 0.0
> 		vec[1] = 0.0
> 		vec[2] = 0.0
> 	def dim(self):
> 		return 3
> 
> vOnTop = TopDrag()
> vOnBottom = BottomHold()
> pOut = 100.0 #right side pressure
> pIn = 100.0 #left side pressure
> 
> bc = [DirichletBC(ME.sub(0), vOnBottom, bottom), DirichletBC(ME.sub(0),
> vOnTop, top), DirichletBC(ME.sub(0).sub(1), null, front),
> DirichletBC(ME.sub(0).sub(1), null, back), DirichletBC(ME.sub(1), pIn,
> left), DirichletBC(ME.sub(1), pOut, right)]
> 
> u = Function(ME)
> v = u.split()[0]
> p = u.split()[1]
> omega = TestFunction(ME)
> w, q = split(omega)
> du = TrialFunction(ME)
> dv, dp = split(du)
> 
> F = FormF(v,w,p,q)
> G = FormG(v,w,dv,p,q,dp)
> a = lhs(F-G)
> L = rhs(F-G)
> A = assemble(a)
> b = assemble(L)
> for condition in bc:
> 	condition.apply(A, b)
> 
> du = Function(ME)
> solve(A, du.vector(), b)
> 
> print 'some value:', u((0.5,0.5,0.5))
> print 'some value:', du((0.5,0.5,0.5))