> 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.
For a symmetric linear operator A, the image is orthogonal on the kernel. Consequently, if I want something to be in the image of an operator, I have to make sure it has no component in the kernel, i.e., it is orthogonal to it. --Nico On Wed, Jan 22, 2014 at 5:34 PM, Nikolaus Rath <[email protected]> wrote: > 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
