# 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)"
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
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) \
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