dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

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)