On 01/22/2014 06:56 AM, Nico Schlömer wrote:
> Here's a complete example:
>
> =============== *snip* ===============
> from dolfin import *
>
> mesh = UnitSquareMesh(20, 20)
> V = FunctionSpace(mesh, 'CG', 1)
>
> f = Expression('sin(x[0]) * sin(x[1])')
>
> u = TrialFunction(V)
> v = TestFunction(V)
> a = dot(grad(u), grad(v)) * dx
> L = f * v * dx
>
> A = assemble(a)
> b = assemble(L)
>
> # Make sure the system is consistent, i.e., the right-hand side is in the
> image
> # of A. This is equivalent to the rhs not having a component in the kernel of
> # A (since A is symmetric).
> e = Function(V)
> e.interpolate(Constant(1.0))
> evec = e.vector()
> evec /= norm(evec)
> alpha = b.inner(evec)
> b -= alpha * evec
Could you explain how this enforces that b is in the image of A? It
seems to me that you're instead making sure that b has no component in
the nullspace of A (which is spanned by evec). But that doesn't make
sense to me, the nullspace of A and the image of A could be subspaces of
entirely different spaces.
What I would have done is this:
# Make sure that Neumann boundary conditions satisfy
# compatibility condition
net_flux = assemble(f * dx)
area = assemble(interpolate(Constant(1), V) * ds)
g = Constant(net_flux) / Constant(area)
L = f * v * dx + g * v * ds
b = assemble(L)
Best,
Nikolaus
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
