Constraining the solution and its derivative on the boundary

Hi, mathematically this is not a well-posed problem (in the sense of
Hademard). The Poisson problem should
have one BC for each point on the boundary. Furthermore, when you use
Lagrange elements you only
have one dof for each point on the boundary - making it impossible to set
two different BC there.

Kent

On 9 December 2012 12:50, Aravind <email address hidden>wrote:

> Question #216323 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/216323
>
> Description changed to:
> Hi,
>
> I'm trying to solve the Poisson equation on the unit interval with the
> following boundary conditions:
>
> u''(x) = 1 for 0 <= x <= 1 BC1: u(1) = 0 BC2: u'(1) = 0
>
> Mathematically, this is a well-posed problem, but it is not possible to
> apply BC2 as a natural boundary condition in the variational form,
> because the test function is defined to be zero at x=1 due to the
> constraint posed by BC1. I was able to solve this problem using a mixed
> formulation similar to mixed-poisson demo problem. However, this seems
> to be an inefficient method because it uses 2 degrees of freedom for
> every node in the mesh.
>
> I saw a similar question
> (https://answers.launchpad.net/fenics/+question/212434) where a boundary
> condition was applied using a scalar Lagrange multiplier. I tried to
> apply this method here, but without success. Basically, I introduced the
> value of the derivative at x=0 as an additional unknown, lambda, and I
> tried to constrain the value of lambda so that u'(1) = 0. Could someone
> tell me whether this possible at all and if so, what changes I could
> make to the following code:
>
> from dolfin import *
>
> # Create mesh and define function space
> mesh = UnitInterval(32)
> Wp = FunctionSpace(mesh, "CG", 1)
> R = FunctionSpace(mesh, "R", 0)
> W = Wp * R
>
> # Define trial and test functions
> (u, lmbda) = TrialFunctions(W)
> (v, w) = TestFunctions(W)
>
> # Define source function
> f = Constant(-1.0)
>
> # # Define variational form ***** CORRECT SOLUTION OBTAINED BY SPECIFYING
> u'(0) INSTEAD OF u'(1) *****
> # a = inner(grad(u), grad(v))*dx
> # L = f*v*dx + v*ds(0)
>
> # Define variational form ***** ATTEMPT USING LAGRANGE MULTIPLIER - DOES
> NOT WORK *****
> a = inner(grad(u), grad(v))*dx - lmbda*v*ds(0) + w*(lmbda - grad(u))*ds(0)
> L = f*v*dx
>
> # Define boundary condition
> def right(x):
> return x[0] > 1.0 - DOLFIN_EPS
> bc = DirichletBC(W.sub(0), Constant(0.0), right)
>
> # Compute solution
> U = Function(W)
> solve(a == L, U, bc)
> (u, lmbda) = U.split(True)
>
> grad_u = project(grad(u), FunctionSpace(mesh, 'Lagrange', 1))
>
> # Plot solution
> print u.vector().array()[[0,-1]]
> print lmbda.vector().array()
> plot(u, title='Solution')
> plot(grad_u, title='Gradient')
> interactive()
>
> Thanks for the help.
>
> Aravind
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.
>

Thanks for your reply, Kent.

Perhaps I'm mistaken in assuming that it is mathematically well-posed, but I can solve the problem analytically and I believe that the solution is unique. It is trivial to show that u(x) = (1-x)^2/2 is a solution to the PDE. I did not try to look further for a rigorous proof of well-posedness. Would it help to view this problem as an "initial"-value problem rather than as a boundary value problem?

I agree that Lagrange elements may not be the most appropriate for this problem. I tried Hermite elements, which are supposed to have an additional DOF for the slope at the boundary. However, I hit a dead end because I get the error 'Domain "interval" invalid for "Hermite" finite element'. If you can point me to some helpful documentation or a working example, that would be great. Perhaps I can address this issue in another thread after a more thorough investigation.

Meanwhile, sticking to Lagrange elements, it may not be possible to mathematically enforce two boundary conditions at a point, but wouldn't it be possible to impose the derivative constraint through a penalization approach (e.g. the Lagrange multiplier approach)? I'm able to solve the PDE with Lagrange elements using a mixed formulation by splitting the PDE into a system of first-order PDEs. This gives me some confidence that I may be on the right track. Please correct me if I'm wrong.

Aravind

I managed to solve the problem using the Lagrange multiplier method. I'm just posting the solution here so that someone else may benefit from it. The explanation may be excessive for most of you, but I hope it will benefit newbies like me.

The original PDE: u''(x) = -f(x) for 0 <= x <= 1 where f(x) = -1 BC1: u(1) = 0 BC2: u'(1) = 0 .... (1)

The problem cannot be solved using the variational method because it lacks a boundary condition at x=0 and has 2 B.C.'s at x = 1. We can transform the problem as follows:

u''(x) = -f(x) for 0 <= x <= 1 where f(x) = -1 BC1: u(1) = 0 BC2: u'(0) = lambda .... (2)

Note that we have replaced the original BC2 with a new condition that is applied at x = 0. In doing so, we have to introduce an additional constant, lambda, which serves as the unknown value of the gradient at x = 0. In order to determine lambda, we integrate the PDE once as follows:

u'(1) - u'(0) = - int_{x=0}^{x=1} f(x) dx .... (3) where int stands for the integral symbol

Here we apply the original BC2 which states that u'(1) = 0. We also replace u'(0) with lambda, using the definition of lambda, to get:

lambda = int_{x=0}^{x=1} f(x) dx .... (4)

Solving Equations (2) and (4) simultaneously yields the correct solution. This can be expressed in variational form and solved in FEniCS. The Python program in my initial comment stays the same except for the lines where bilinear and linear forms, 'a' and 'L', are defined. The correct version should be:

a = inner(grad(u), grad(v))*dx + lmbda*v*ds(0) + lmbda*w*dx # (I'm not sure why lmbda*w*ds is incorrect and lmbda*w*dx is)
L = f*v*dx + f*w*dx

If you are wondering why I am interested in solving this problem, Eqn. (1) describes the variation of bending moment due to a distributed load on a cantilever beam described using Euler-Bernoulli theory. Eqn. (4), in effect, states that the shear force at the anchor (left end) is equal to the integral of the distributed load, which is required for equilibrium.

I was trying to solve for the displacement in the beam, which is a 4th order PDE: u''''(x) = f(x). It has 4 boundary conditions u(0) = u'(0) = 0 and u'''(1) = u''''(1) = 0. To solve this 4th order PDE, I used a mixed formulation to convert it into a system of 2 second-order PDEs: u''(x) = -M(x) and M''(x) = -f(x), where M(x) is the bending moment in the beam. Using Lagrange multipliers as above, I can solve the entire system using 2 DOFs at each node in the mesh (to represent u(x) and M(x)) and 2 constant Lagrange multipliers (to represent the unknown values of the gradient at the boundary). Without Lagrange multipliers, I have to solve a system of four 1st order PDEs representing u(x) and its first 3 derivatives. I can post the complete code for the Euler-Bernoulli beam if anyone is interested.

I would be particular interested about that EB beam solution if you have one?