On 01/09/2014 11:19 AM, Jan Blechta wrote:
>>>> I would like to solve the following equation (which does not
>>>> directly come from a PDE):
>>>>
>>>> \int dV f(x) * \partial_r g(x) = \int dA u(x) * \partial_n g(x)
>>>>
>>>> f(x) is known, u(x) is unknown, and the equation should hold for
>>>> any g(x) that satisfy Laplace's equation.
>>>>
>>>> In other words, I'm looking for a weight function u(x), such that
>>>> the surface integral of the normal derivative of any g(x)
>>>> (weighted by u) gives the same result as the volume integral of
>>>> the radial derivative of
>>>> g(x) (weighted by the known function f(x)).
>>>>
>>>> Is it possible to do this with FEniCS?
>>>> It seems that the equation itself is easy to express in UFL, but I
>>>> am not sure how do deal with the fact that there are no boundary
>>>> conditions, and that any trial function g(x) needs to satisfy
>>>> Laplace's equation.
>>>
>>> Yes. Look at the demo
>>>
>>> demo/documented/neumann-poisson
>>>
[...]
>> However, even after going through the example, I'm not sure how I can
>> tell FEniCS to use only trial functions satisfying Laplace's equation.
>> In the example, the constraint on the test functions seems to fix just
>> the constant offset. It's not clear to me to extend this to something
>> more complicated, where the constraint itself takes the form of a PDE.
>> Do you think you could explain in a bit more detail?
>
> You change the constraint to Laplace problem, and the space R for
> Lagrange multiplier and corresponding test function to some appropriate
> FE subspace of Sobolev space W_0^{1,2}. Hint:
> http://en.wikipedia.org/wiki/Lagrange_multipliers_on_Banach_spaces
I have started trying to implement this, but I'm running into some
problems. I have written the following code (which does not yet impose
the extra constraint on the test functions):
from __future__ import division, print_function, absolute_import
from dolfin import *
set_log_level(PROGRESS)
mesh = UnitSquareMesh(30,30)
V = FunctionSpace(mesh, 'Lagrange', 1)
u = TrialFunction(V)
v = TestFunction(V)
jt = Expression('sqrt(pow(x[0], 2) + pow(x[1], 2))')
a = jt * grad(v)[0] * dx
L = u * dot(FacetNormal(mesh), grad(v)) * ds
bopt = Function(V)
solve(a == L, bopt,
solver_parameters = {'linear_solver': 'cg',
'preconditioner': 'ilu'})
but the call to solve() fails with:
> *** Error: Unable to define linear variational problem a(u, v) == L(v) for
> all v.
> *** Reason: Expecting the left-hand side to be a bilinear form (not rank 1).
> *** Where: This error was encountered inside LinearVariationalProblem.cpp.
> *** Process: 0
> ***
> *** DOLFIN version: 1.3.0
> *** Git changeset:
The message is pretty clear, but I don't know what to do about it. The
equation that I want to solve simply doesn't have a TrialFunction in L..
I guess I need some trick here?
Best,
Nikolaus
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
