FEniCS Project

Defining a Dirichlet boundary condition over the vector coefficients

Asked by Camille Gouttebroze

Hi,

I would like to solve a 1D problem where the solution is 2D vector field. One of the boundary condition is of the following kind:
a * u[0] + b * u[1] = c
where a, b and c are constant. To achieve this, I use the following penalization method.

# Create mesh and define function space
mesh = UnitInterval(32)
V = VectorFunctionSpace(mesh, "Lagrange", 1, dim=2)

# Define left boundary domain
def left(x):
    return x[0] < DOLFIN_EPS

# Define left boundary condition
bc = DirichletBC(V, Constant((0.,0.)), left)

# Penalization for the right boundary condition
a = 1.
b = 2.
c = 4.
p = Expression("x[0] == 1.0 ? 1.e6 : 0.")

# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement

# Total potential energy
Pi = inner( grad(u), grad(u) )*dx + p*(a*u[0]+b*u[1]-c)**2*dx

# Compute first variation of Pi (directional derivative about u in the direction of v)
f = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(f, u, du)

# Solve variational problem
solve(f == 0, u, bc, J=J)

# Plot solution
plot(u[0])
plot(u[1])
interactive()

This works fine on this small example, but i would prefer to use Lagrange multipliers for more complex cases. Is there any way to do this with DirichletBC or any other class ?

Best regards,

Camille

Question information

Language:
English Edit question
Status:
Solved
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Solved by:
Anders Logg
Solved:
Last query:
Last reply:
Revision history for this message
Anders Logg (logg) said : 		#1

You can try defining a scalar Lagrange multiplier and use it in combination with *ds to impose the boundary condition. For an example of Lagrange multipliers, see the demo

http://bazaar.launchpad.net/~dolfin-core/dolfin/trunk/view/head:/demo/undocumented/neumann-poisson/python/demo_neumann-poisson.py

Revision history for this message
Camille Gouttebroze (camille-gouttebroze) said : 		#2

I tried to transpose the demo to my case. I did the following :

# Create mesh and define function spaces
mesh = UnitInterval(32)
V = VectorFunctionSpace(mesh, "Lagrange", 1, dim=2)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Define functions
(u,c) = TrialFunction(W)

wich ended with this error :

(u,c) = TrialFunction(W)
ValueError: too many values to unpack

Revision history for this message
Martin Sandve Alnæs (martinal) said : 		#3

On 30 October 2012 10:21, Camille Gouttebroze
<email address hidden> wrote:
> Question #212434 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/212434
>
> Status: Answered => Open
>
> Camille Gouttebroze is still having a problem:
> I tried to transpose the demo to my case. I did the following :
>
> # Create mesh and define function spaces
> mesh = UnitInterval(32)
> V = VectorFunctionSpace(mesh, "Lagrange", 1, dim=2)
> R = FunctionSpace(mesh, "R", 0)
> W = V * R
>
> # Define functions
> (u,c) = TrialFunction(W)
>
> wich ended with this error :
>
> (u,c) = TrialFunction(W)
> ValueError: too many values to unpack

You need to do
(u,c) = TrialFunctions(W) # Note the s!

or just
u,c = TrialFunctions(W)

or equivalently
uc = TrialFunction(W)
u,c = split(uc)

Martin

Revision history for this message
Kristian B. Ølgaard (k.b.oelgaard) said : 		#4

On 30 October 2012 10:21, Camille Gouttebroze
<email address hidden> wrote:
> Question #212434 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/212434
>
> Status: Answered => Open
>
> Camille Gouttebroze is still having a problem:
> I tried to transpose the demo to my case. I did the following :
>
> # Create mesh and define function spaces
> mesh = UnitInterval(32)
> V = VectorFunctionSpace(mesh, "Lagrange", 1, dim=2)
> R = FunctionSpace(mesh, "R", 0)
> W = V * R
>
> # Define functions
> (u,c) = TrialFunction(W)

should be
(u,c) = TrialFunctions(W)

> wich ended with this error :
>
> (u,c) = TrialFunction(W)
> ValueError: too many values to unpack
>
> --
> You received this question notification because you are a member of
> FEniCS Team, which is an answer contact for FEniCS Project.

Revision history for this message
Camille Gouttebroze (camille-gouttebroze) said : 		#5

Thank you Kristian and Martin for your answers. This allowed me to transpose to Anders' demo to my case under the following form :

# Create mesh and define function spaces
mesh = UnitInterval(32)
V = VectorFunctionSpace(mesh, "Lagrange", 1, dim=2)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Define left boundary domain
def left(x):
    return x[0] < DOLFIN_EPS

# Define left boundary condition
bc = DirichletBC(V, Constant((0.,0.)), left)

# Constants for the right boundary condition
a,b,c = 1.,2.,4.

# Define functions
(u,lmbda) = TrialFunctions(W) # Incremental displacement
(v,mu) = TestFunctions(W) # Test function
w = Function(W)

# Forms
a = inner( grad(u), grad(v) ) * dx + mu * ( a*u[0]+b*u[1]-c ) * ds + lmbda * ( a*v[0]+b*v[1]-c ) * ds
L = inner( Constant((0.,0.)), v ) * dx

# Solve variational problem
solve(a == L, w, bc)
(u,lmbda) = w.split()

# Plot solution
plot(u[0])
plot(u[1])
interactive()

It works for c = 0, as its the case in Anders' demo. But in the others cases, i get the following error :
Exception: All terms in form must have same rank.

How should I write the bilinear form for the cases where c is not 0 ?

Revision history for this message
Best Anders Logg (logg) said : 		#6

On Tue, Oct 30, 2012 at 11:11:14AM -0000, Camille Gouttebroze wrote:
> Question #212434 on FEniCS Project changed:
> https://answers.launchpad.net/fenics/+question/212434
>
> Status: Answered => Open
>
> Camille Gouttebroze is still having a problem:
> Thank you Kristian and Martin for your answers. This allowed me to
> transpose to Anders' demo to my case under the following form :
>
> # Create mesh and define function spaces
> mesh = UnitInterval(32)
> V = VectorFunctionSpace(mesh, "Lagrange", 1, dim=2)
> R = FunctionSpace(mesh, "R", 0)
> W = V * R
>
> # Define left boundary domain
> def left(x):
> return x[0] < DOLFIN_EPS
>
> # Define left boundary condition
> bc = DirichletBC(V, Constant((0.,0.)), left)
>
> # Constants for the right boundary condition
> a,b,c = 1.,2.,4.
>
> # Define functions
> (u,lmbda) = TrialFunctions(W) # Incremental displacement
> (v,mu) = TestFunctions(W) # Test function
> w = Function(W)
>
> # Forms
> a = inner( grad(u), grad(v) ) * dx + mu * ( a*u[0]+b*u[1]-c ) * ds + lmbda * ( a*v[0]+b*v[1]-c ) * ds
> L = inner( Constant((0.,0.)), v ) * dx
>
> # Solve variational problem
> solve(a == L, w, bc)
> (u,lmbda) = w.split()
>
> # Plot solution
> plot(u[0])
> plot(u[1])
> interactive()
>
> It works for c = 0, as its the case in Anders' demo. But in the others cases, i get the following error :
> Exception: All terms in form must have same rank.
>
> How should I write the bilinear form for the cases where c is not 0 ?

Move the c-term to the right-hand side.

Or better, write the whole thing as

F = inner(grad(u), grad(v))*dx ... - inner(Constant((0.,0.)), v)*dx

Then extract the left- and right-hand sides:

a = lhs(F)
L = rhs(F)

--
Anders

Revision history for this message
Camille Gouttebroze (camille-gouttebroze) said : 		#7

Thanks Anders Logg, that solved my question.