dolfin team mailing list archive

[Xavier Raynaud <question151961@xxxxxxxxxxxxxxxxxxxxx>]

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

Xavier Raynaud posted a new comment:
I actually found a way to deal with this question but it does not seem
natural: I apply zero dirichlet boundary condition on the interior (!)
of the domain. Is there any more natural way in Fenics of defining
spaces of functions defined on the boundary?

You can find here a description of the problem
http://dl.dropbox.com/u/4003624/testcase.pdf

Here is the code:

from dolfin import *


mesh = UnitSquare(32,32)


movie=False

# Time step, Time interval
dt=0.1
T=80

# Choose finite element vector space
U = FunctionSpace(mesh, "Lagrange", 1)
P = FunctionSpace(mesh, "Lagrange", 1)

W=U*P

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


def boundary(x):
   if x[0]<0.0001 or x[0]>1-0.0001 or x[1]<0.0001 or x[1]>1-0.0001:
      return False
   else:
      return True

bc = DirichletBC(W.sub(1), Constant(1.0), boundary)


#dissipation coefficient
alpha=0.001

uprev=interpolate(Constant(0.0),U)
pprev=interpolate(Constant(0.0),P)

# visualization parameters
h=plot(uprev)
# h.rescale=True
h.rescale=False
h.set_min_max(0,100)
h.init_writer(".tmp")


# Construct bilinear and linear form
(u,p) = TrialFunction(W)
(v,q) = TestFunction(W)

# k1=0
k2=1
gamma=1

# a = u*v*dx+dt*alpha*inner(grad(u),
grad(v))*dx-(-k1*gamma*u*p*v+k2*gamma*(1-p))*ds+(p-dt*(-k1*gamma*u*p*v+k2*gamma*(1-p)))*w*ds

a = u*v*dx+dt*alpha*inner(grad(u),grad(v))*dx+k2*gamma*p*v*ds+(1+dt*k2)*p*q*ds
L = uprev*v*dx+k2*gamma*v*ds+(pprev+dt*k2)*q*ds


A=assemble(a)

w=Function(W)


# debug_here()

# Time integration
t=0
i=1
while t<=T:

   b=assemble(L)
   bc.apply(A,b)
   solve(A,w.vector(),b)

   (u,p)=w.split()
   h.update(data=u)
   if movie:
      h.write_png()

   i+=1
   t+=dt
   uprev.assign(u)
   pprev.assign(p)

if movie:
   h.movie("movie2.avi", fps=10, cleanup=True)

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