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)) _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
