dolfin team mailing list archive

[Marie Rognes <meg@xxxxxxxxx>]

Chris Richardson wrote:
> 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:
>
>   

I've attached a revised version. Take a look and see if it makes sense.

Note that with (the previous):

   def boundary0(x,on_bound):
   (*)    return x[0]>DOLFIN_EPS and x[0]<1.0-DOLFIN_EPS and on_bound
   bcs = DirichletBC(V.sub(0), u0, boundary0)
     

only the facets where both endpoints (and the midpoint) satisfy (*) will 
be marked. (Somebody correct me if I'm wrong here.) In particular, the 
left-most "bottom" edge will not be marked. 

Further, then v*n will not vanish on this edge, and hence you will get 
contributions from this edge in:

	L = dot(v,f*n)*ds

Then, the definition of f on the edges starts mattering etc.

The moral is that one should be really careful with the boundary 
conditions ;)

--
Marie



from dolfin import *

n = 16
mesh = UnitSquare(n,n)
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
#V =  R * W #in later versions of DOLFIN

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

# boundary conditions
class PressureBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] < DOLFIN_EPS or x[0] > (1.0 - DOLFIN_EPS))

class VelocityBoundary(SubDomain):
    def inside(self, x, on_boundary):
        left = x[0] < DOLFIN_EPS and (x[1] > DOLFIN_EPS
                                      and x[1] < 1.0 - DOLFIN_EPS)
        right = abs(x[0] - 1.0) < DOLFIN_EPS and (abs(x[1] - .5) > DOLFIN_EPS
                                                  and x[1] < 1.5 - DOLFIN_EPS)
        return on_boundary and not (left or right)

pressure_boundary = PressureBoundary()
velocity_boundary = VelocityBoundary()

markers = MeshFunction("uint", mesh, 1)
markers.set_all(2)
pressure_boundary.mark(markers, 1)
velocity_boundary.mark(markers, 0)

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

f = Expression("-1.0*(x[0]==0.0)+1.0*(x[0]==1.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, n)*f*ds(1)

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

# Plot solution
plot(u)
plot(p)
interactive()