Guest User

Untitled

a guest
Apr 10th, 2018
1,145
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #! /usr/bin/env python
  2. """
  3. Demo for Nitsche-type free-slip boundary conditions
  4. """
  5.  
  6. __author__ = "Christian Waluga (waluga@ma.tum.de)"
  7. __copyright__ = "Copyright (c) 2013 %s" % __author__
  8.  
  9. from dolfin import *
  10.  
  11. # avoid outputs like 'solving variational problem' by FEniCS
  12. set_log_level(ERROR)
  13.  
  14. mesh = refine(Mesh("cylinderbump.xml"))
  15.  
  16. # naive equal order element
  17. V = VectorFunctionSpace(mesh, 'CG', 1)
  18. Q = FunctionSpace(mesh, 'CG', 1)
  19.  
  20. mu = Constant(1.0)
  21. F = Constant((0.0,0.0))
  22.  
  23. W = V*Q
  24.  
  25. # variational problem
  26. u, p = TrialFunctions(W)
  27. v, q = TestFunctions(W)
  28.  
  29. # boundary conditions
  30. leftexp = Expression(("x[1]*(3-x[1])/10", "0.0"))
  31. left = DirichletBC(W.sub(0), leftexp, "near(x[0], -5)")
  32. top = DirichletBC(W.sub(0), (0,0), "near(x[1], 3)")
  33. bcs = [left, top]
  34.  
  35. class Bottom(SubDomain):
  36.   def inside(self, x, on_boundary):
  37.     return on_boundary and not (near(x[0], -5) or near(x[0], 5) or near(x[1], 3))
  38. bottom = Bottom()
  39.  
  40. boundaries = FacetFunction("size_t", mesh)
  41. boundaries.set_all(0)
  42. bottom.mark(boundaries, 1)
  43. ds = Measure("ds")[boundaries]
  44.  
  45. alpha = Constant(1./10)
  46. beta = Constant(10)
  47. h = CellSize(mesh)
  48. n = FacetNormal(mesh)
  49.  
  50. # (bi)linear forms
  51. def a(u,v): return inner(mu*grad(u),grad(v))*dx
  52. def b(v,q): return - div(v)*q*dx
  53. def f(v):   return dot(F, v)*dx
  54. def t(u,p): return dot(2*mu*sym(grad(u)),n) - p*n
  55.  
  56. stokes = a(u,v) + b(v,p) + b(u,q) - f(v) \
  57.        - dot(n,t(u,p))*dot(v,n)*ds(1) - dot(u,n)*dot(n,t(v,q))*ds(1) \
  58.        + beta/h*dot(u,n)*dot(v,n)*ds(1) \
  59.        + alpha*h**2*dot(F - grad(p), grad(q))*dx
  60.  
  61. # solve variational problem
  62. wh = Function(W)
  63. print 'size:', wh.vector().size()
  64.  
  65. solve(lhs(stokes) == rhs(stokes), wh, bcs, \
  66.   solver_parameters = { \
  67.     "linear_solver": "tfqmr", \
  68.     "preconditioner": "ml_amg"} )
  69.  
  70. uh, ph = wh.split()
  71. plot(uh, interactive = True)
  72. File("cylinderbump-velocity.pvd") << uh
  73. File("cylinderbump-pressure.pvd") << ph
  74.  
  75. exit()
RAW Paste Data