DOLFIN

Beginners question: applying boundary conditions

Asked by Thomas Ward

I am struggling with a very simple problem

unit cube
div(grad(u))=0
du/dn=-1.4 on x=0
du/dn=1.4 on x=1
du/dn+1.4*1.4/9.81*u=0 on z=1
du/dn=0 on y=0,y=1,z=0

below is my attempt with pydolfin, but it gives u=zero everywhere
any help will be much appreciated, thanks, Tom

from dolfin import *
mesh = UnitCube(5, 5, 5)
V = FunctionSpace(mesh, "CG", 1)

subdomains = MeshFunction("uint",mesh, mesh.topology().dim())
left, right = compile_subdomains(["(fabs(x[0]) < DOLFIN_EPS) && on_boundary", "(fabs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary"])
btm, top = compile_subdomains(["(fabs(x[2]) < DOLFIN_EPS) && on_boundary", "(fabs(x[2] - 1.0) < DOLFIN_EPS) && on_boundary"])

top.mark(subdomains, 1)
right.mark(subdomains, 2)
left.mark(subdomains, 3)

v = TestFunction(V);
u = TrialFunction(V);

g = Expression("1.4")
alpha=Expression("1.4*1.4/9.81")

a = inner(grad(v), grad(u))*dx +(alpha*v*u)*ds(1)
L = g*v*ds(2)-g*v*ds(3)

problem = VariationalProblem(a, L)
u = problem.solve()
plot(u, interactive=True)

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Hake
Solved:
Last query:
Last reply:
Revision history for this message
Best Johan Hake (johan-hake) said : 		#1

You are almost there:

Try this slightly corrected code:

Johan

*******************************************
from dolfin import *
mesh = UnitCube(5, 5, 5)
V = FunctionSpace(mesh, "CG", 1)

# You want the facets not the cells
subdomains = MeshFunction("uint", mesh, mesh.topology().dim()-1)

left, right = compile_subdomains(["(fabs(x[0]) < DOLFIN_EPS) && on_boundary",
"(fabs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary"])
btm, top = compile_subdomains(["(fabs(x[2]) < DOLFIN_EPS) && on_boundary",
"(fabs(x[2] - 1.0) < DOLFIN_EPS) && on_boundary"])

# Mark all subdomains with a number you are not using in the form
subdomains.set_all(4)
top.mark(subdomains, 1)
right.mark(subdomains, 2)
left.mark(subdomains, 3)

v = TestFunction(V);
u = TrialFunction(V);

# Use Constants for, well, constants
g = Constant(1.4)
alpha=Constant(1.4*1.4/9.81)

a = inner(grad(v), grad(u))*dx +(alpha*v*u)*ds(1)
L = g*v*ds(2)-g*v*ds(3)

# Pass your subdomains to VariationalProblem
problem = VariationalProblem(a, L, exterior_facet_domains=subdomains)
u = problem.solve()
plot(u, interactive=True)

On Tuesday May 4 2010 11:11:16 Thomas Ward wrote:
> New question #109629 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/109629
>
> I am struggling with a very simple problem
>
> unit cube
> div(grad(u))=0
> du/dn=-1.4 on x=0
> du/dn=1.4 on x=1
> du/dn+1.4*1.4/9.81*u=0 on z=1
> du/dn=0 on y=0,y=1,z=0
>
> below is my attempt with pydolfin, but it gives u=zero everywhere
> any help will be much appreciated, thanks, Tom
>
> from dolfin import *
> mesh = UnitCube(5, 5, 5)
> V = FunctionSpace(mesh, "CG", 1)
>
> subdomains = MeshFunction("uint",mesh, mesh.topology().dim())
> left, right = compile_subdomains(["(fabs(x[0]) < DOLFIN_EPS) &&
> on_boundary", "(fabs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary"]) btm, top
> = compile_subdomains(["(fabs(x[2]) < DOLFIN_EPS) && on_boundary",
> "(fabs(x[2] - 1.0) < DOLFIN_EPS) && on_boundary"])
>
> top.mark(subdomains, 1)
> right.mark(subdomains, 2)
> left.mark(subdomains, 3)
>
> v = TestFunction(V);
> u = TrialFunction(V);
>
> g = Expression("1.4")
> alpha=Expression("1.4*1.4/9.81")
>
> a = inner(grad(v), grad(u))*dx +(alpha*v*u)*ds(1)
> L = g*v*ds(2)-g*v*ds(3)
>
> problem = VariationalProblem(a, L)
> u = problem.solve()
> plot(u, interactive=True)

Revision history for this message
Garth Wells (garth-wells) said : 		#2

On 04/05/10 19:45, Johan Hake wrote:
> Question #109629 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/109629
>
> Status: Open => Answered
>
> Johan Hake proposed the following answer:
> You are almost there:
>
> Try this slightly corrected code:
>

Also, try applying some Dirichlet bcs to the problem.

Garth

> Johan
>
> *******************************************
> from dolfin import *
> mesh = UnitCube(5, 5, 5)
> V = FunctionSpace(mesh, "CG", 1)
>
> # You want the facets not the cells
> subdomains = MeshFunction("uint", mesh, mesh.topology().dim()-1)
>
> left, right = compile_subdomains(["(fabs(x[0])< DOLFIN_EPS)&& on_boundary",
> "(fabs(x[0] - 1.0)< DOLFIN_EPS)&& on_boundary"])
> btm, top = compile_subdomains(["(fabs(x[2])< DOLFIN_EPS)&& on_boundary",
> "(fabs(x[2] - 1.0)< DOLFIN_EPS)&& on_boundary"])
>
> # Mark all subdomains with a number you are not using in the form
> subdomains.set_all(4)
> top.mark(subdomains, 1)
> right.mark(subdomains, 2)
> left.mark(subdomains, 3)
>
> v = TestFunction(V);
> u = TrialFunction(V);
>
>
> # Use Constants for, well, constants
> g = Constant(1.4)
> alpha=Constant(1.4*1.4/9.81)
>
> a = inner(grad(v), grad(u))*dx +(alpha*v*u)*ds(1)
> L = g*v*ds(2)-g*v*ds(3)
>
> # Pass your subdomains to VariationalProblem
> problem = VariationalProblem(a, L, exterior_facet_domains=subdomains)
> u = problem.solve()
> plot(u, interactive=True)
>
> On Tuesday May 4 2010 11:11:16 Thomas Ward wrote:
>> New question #109629 on DOLFIN:
>> https://answers.launchpad.net/dolfin/+question/109629
>>
>> I am struggling with a very simple problem
>>
>> unit cube
>> div(grad(u))=0
>> du/dn=-1.4 on x=0
>> du/dn=1.4 on x=1
>> du/dn+1.4*1.4/9.81*u=0 on z=1
>> du/dn=0 on y=0,y=1,z=0
>>
>> below is my attempt with pydolfin, but it gives u=zero everywhere
>> any help will be much appreciated, thanks, Tom
>>
>> from dolfin import *
>> mesh = UnitCube(5, 5, 5)
>> V = FunctionSpace(mesh, "CG", 1)
>>
>> subdomains = MeshFunction("uint",mesh, mesh.topology().dim())
>> left, right = compile_subdomains(["(fabs(x[0])< DOLFIN_EPS)&&
>> on_boundary", "(fabs(x[0] - 1.0)< DOLFIN_EPS)&& on_boundary"]) btm, top
>> = compile_subdomains(["(fabs(x[2])< DOLFIN_EPS)&& on_boundary",
>> "(fabs(x[2] - 1.0)< DOLFIN_EPS)&& on_boundary"])
>>
>> top.mark(subdomains, 1)
>> right.mark(subdomains, 2)
>> left.mark(subdomains, 3)
>>
>> v = TestFunction(V);
>> u = TrialFunction(V);
>>
>> g = Expression("1.4")
>> alpha=Expression("1.4*1.4/9.81")
>>
>> a = inner(grad(v), grad(u))*dx +(alpha*v*u)*ds(1)
>> L = g*v*ds(2)-g*v*ds(3)
>>
>> problem = VariationalProblem(a, L)
>> u = problem.solve()
>> plot(u, interactive=True)
>

Revision history for this message
Thomas Ward (tomwardathome-gmail) said : 		#3

Hey Johan, ... many thanks that did the trick....
I am off to figure out what the difference between a facet and a cell is :)
Tom

Revision history for this message
Thomas Ward (tomwardathome-gmail) said : 		#4

Thanks Johan Hake, that solved my question.