dolfin team mailing list archive

["B. Emek Abali" <question140910@xxxxxxxxxxxxxxxxxxxxx>]

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))

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.