Hi everyone,
i´m trying to solve a system, depending on the solution from an other system. It seems to work fine, there are no error-messages. But when i want to evaluate the solution, it is NaN everywere and i have no idea why.
from dolfin import *
# Set Mesh and Functionspaces:
mesh = CircleMesh(Point(0.0, 0.0), 1.0, 0.2)
V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q
class Inlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] < -0.9 + DOLFIN_EPS
class Outlet(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] > 0.9 - DOLFIN_EPS
class Shape(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0]) <= 0.9 - DOLFIN_EPS
Gamma_in = Inlet()
Gamma_out = Outlet()
Gamma_shape = Shape()
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
Gamma_in.mark(boundaries, 1)
Gamma_out.mark(boundaries, 2)
Gamma_shape.mark(boundaries, 3)
# Get Solution u:
inflow = Expression(('-x[0]-0.9', '0.0'))
bc0 = DirichletBC(W.sub(0), inflow, boundaries, 1)
shapeflow = Constant((0.0,0.0))
bc1 = DirichletBC(W.sub(0), shapeflow, boundaries, 3)
zero = Constant(0.0)
bc2 = DirichletBC(W.sub(1), zero, boundaries, 2)
bcs = [bc0, bc1, bc2]
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0.0, 0.0))
g = Constant(0.0)
a = inner(nabla_grad(u),nabla_grad(v))*dx + inner(grad(p),v)*dx + div(u)*q*dx
L = inner(f, v)*dx + g*q*dx
U = Function(W)
solve(a == L, U, bcs)
u, p = U.split()
# Get Solution lamb, depending on u:
ad0 = DirichletBC(W.sub(0), f, boundaries, 1)
ad1 = DirichletBC(W.sub(0), f, boundaries, 3)
ad2 = DirichletBC(W.sub(1), g, boundaries, 2)
ads = [ad0, ad1, ad2]
(lamb, lambp) = TrialFunctions(W)
(v, q) = TestFunctions(W)
n = FacetNormal(mesh)
ds = Measure("ds")[boundaries]
C = as_matrix([[u[0]*n[0], u[1]*n[0]], [u[0]*n[1], u[1]*n[1]]])
h = as_vector([2*curl(u).dx(1), -2*curl(u).dx(0)])
a = inner(nabla_grad(lamb),nabla_grad(v))*dx + inner(nabla_grad(lamb)*u,v)*dx + inner(nabla_grad(lamb).T*u,v)*dx + inner(lambp,div(v))*dx + inner(inner(u,n)*lamb,v)*ds(2) + inner(C*lamb,v)*ds(2)
L = inner(h, v)*dx + g*q*dx
La = Function(W)
solve(a == L, La, ads)
lamb, lambp = La.split()
print('Lambda at Center:')
print(lamb(0,0))
