dolfin team mailing list archive

[Chris Richardson <question105900@xxxxxxxxxxxxxxxxxxxxx>]

Question #105900 on DOLFIN changed:
https://answers.launchpad.net/dolfin/+question/105900

    Status: Answered => Open

Chris Richardson is still having a problem:
OK, thanks that clears up one thing...

My question now becomes how to set the flux BC on the other boundaries when they are not so simple and 
straight.

e.g. I thought this would work, but is giving real problems at the
corners:

mesh = UnitSquare(20,20)
for x in mesh.coordinates():
     x[1] += 0.5*x[0]
n= FacetNormal(mesh)

R = FunctionSpace(mesh,"BDM",1)
W = FunctionSpace(mesh,"DG",0)
V = R+W

(u,p) = TrialFunctions(V)
(v,q) = TestFunctions(V)

# boundary conditions                                                                                                               
def boundary0(x,on_bound):
    return x[0]>DOLFIN_EPS and x[0]<1.0-DOLFIN_EPS and on_bound

u0 = Constant((0.0,0.0))
bcs = DirichletBC(V.sub(0), u0, boundary0)


f=Expression("-1.0*(x[0]==0.0)+1.0*(x[0]==1.0)")
# f=Constant(0.0)                                                                                                                   

# 'Permeability' k                                                                                                                  
k=Expression("1.0+exp(-40*pow(x[1]-0.7,2))")

a = (k*dot(v,u) - div(v)*p + q*div(u))*dx
L = dot(v,f*n)*ds

# Compute solution                                                                                                                  
problem = VariationalProblem(a, L,bcs)
(u, p) = problem.solve().split()

# Plot solution                                                                                                                     
plot(p, interactive=True)
plot(u, interactive=True)

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.