Hans Petter Langtangen and I are working on a new, revised, and extended version of the FEniCS Tutorial. The plan is two volumes to cover both basic and advanced FEniCS programming. Here’s a sample from Chapter 3: A Gallery of PDE Solvers:

And here’s a bit of code to go along with the movie.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 |
""" FEniCS tutorial demo program: Incompressible Navier-Stokes equations for flow around a cylinder using the Incremental Pressure Correction Scheme (IPCS). rho (u' + u . grad(u)) - div(sigma) = f div(u) = 0 """ from __future__ import print_function from fenics import * from mshr import * import numpy as np T = 5.0 # final time num_steps = 5000 # number of time steps dt = T / num_steps # time step size rho = 1.0 # density nu = 0.001 # kinematic viscosity # Create mesh channel = Rectangle(Point(0, 0), Point(2.2, 0.41)) cylinder = Circle(Point(0.2, 0.2), 0.05) geometry = channel - cylinder mesh = generate_mesh(geometry, 64) # Define function spaces V = VectorFunctionSpace(mesh, 'P', 2) Q = FunctionSpace(mesh, 'P', 1) # Define boundaries inflow = 'near(x[0], 0)' outflow = 'near(x[0], 2.2)' walls = 'near(x[1], 0) || near(x[1], 0.41)' cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3' # Define inflow profile inflow_profile = ('4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)', '0') # Define boundary conditions bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow) bcu_walls = DirichletBC(V, Constant((0, 0)), walls) bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder) bcp_outflow = DirichletBC(Q, Constant(0), outflow) bcu = [bcu_inflow, bcu_walls, bcu_cylinder] bcp = [bcp_outflow] # Define trial and test functions u = TrialFunction(V) v = TestFunction(V) p = TrialFunction(Q) q = TestFunction(Q) # Define functions for solutions at previous and current time steps u0 = Function(V) u1 = Function(V) p0 = Function(Q) p1 = Function(Q) # Define expressions used in variational forms U = 0.5*(u0 + u) n = FacetNormal(mesh) f = Constant((0, 0)) k = Constant(dt) rho = Constant(rho) nu = Constant(nu) # Define symmetric gradient def epsilon(u): return sym(grad(u)) # Define stress tensor def sigma(u, p): return 2*nu*sym(grad(u)) - p*Identity(len(u)) # Define variational problem for step 1 F1 = rho*dot((u - u0) / k, v)*dx + rho*dot(grad(u0)*u0, v)*dx \ + inner(sigma(U, p0), epsilon(v))*dx \ + dot(p0*n, v)*ds - dot(nu*grad(U).T*n, v)*ds \ - dot(f, v)*dx a1 = lhs(F1) L1 = rhs(F1) # Define variational problem for step 2 a2 = dot(grad(p), grad(q))*dx L2 = dot(grad(p0), grad(q))*dx - (1/k)*div(u1)*q*dx # Define variational problem for step 3 a3 = dot(u, v)*dx L3 = dot(u1, v)*dx - k*dot(grad(p1 - p0), v)*dx # Assemble matrices A1 = assemble(a1) A2 = assemble(a2) A3 = assemble(a3) # Apply boundary conditions to matrices [bc.apply(A1) for bc in bcu] [bc.apply(A2) for bc in bcp] # Create VTK files for saving solution vtkfile_u = File('velocity.pvd') vtkfile_p = File('pressure.pvd') # Create progress bar progress = Progress('Time-stepping') set_log_level(PROGRESS) # Time-stepping t = 0 for n in xrange(num_steps): # Update current time t += dt # Step 1: Tentative velocity step b1 = assemble(L1) [bc.apply(b1) for bc in bcu] solve(A1, u1.vector(), b1, 'bicgstab') # Step 2: Pressure correction step b2 = assemble(L2) [bc.apply(b2) for bc in bcp] solve(A2, p1.vector(), b2, 'bicgstab') # Step 3: Velocity correction step b3 = assemble(L3) solve(A3, u1.vector(), b3, 'bicgstab') # Plot solution plot(u1, title='Velocity') plot(p1, title='Pressure') # Save solution to file vtkfile_u << u1 vtkfile_p << p1 # Update previous solution u0.assign(u1) p0.assign(p1) # Update progress progress.update(t / T) # Hold plot interactive() |

Hello Dr. Logg,

We are very grateful to you and your associates for making fenics open source with so many tutorials and demos. I am a MSc student trying to implement wave equation in elasoplastic media using fenics. I am finding it very difficult owing to the lack of detailed documentation and undocumented demos. I know the team stays very busy with so much work but it would be great if you could find some time for development and refinement of FSM app in future.

Thanks and Regards,

The new tutorial will contain many new examples. Maybe not exactly the one you are looking for… but hopefully it will help you figure out how to implement your particular problem in FEniCS.

[…] the velocity field, we use a previously computed solution of the incompressible Navier-Stokes equations which we can read/write using a FEniCS TimeSeries. The variational problem for the equations […]

Hello Dr. Logg,

What is the likely release date for the new tutorial?

Thank you,

I hope the first volume will be online within 1-2 weeks.

Hi Dr. Logg, I was wondering if there was an update on the release of the new tutorial. Thanks, Jonas

Yes, we’re currently proofreading it. It should be ready within 1-2 weeks.