Révision | 5f3b4941f4485411daafa5c3a42017647da24277 (tree) |
---|---|
l'heure | 2007-11-06 18:03:33 |
Auteur | iselllo |
Commiter | iselllo |
I added the code dolfin_demo.py which solves Poisson equation on a 2D
box using the finite-element code dolfin.
@@ -0,0 +1,54 @@ | ||
1 | +#! /usr/bin/env python | |
2 | + | |
3 | +from dolfin import * | |
4 | + | |
5 | +# Create mesh and finite element | |
6 | +mesh = UnitSquare(32, 32) | |
7 | +element = FiniteElement("Lagrange", "triangle", 1) | |
8 | + | |
9 | +# Source term | |
10 | +class Source(Function): | |
11 | + def __init__(self, element, mesh): | |
12 | + Function.__init__(self, element, mesh) | |
13 | + def eval(self, values, x): | |
14 | + dx = x[0] - 0.5 | |
15 | + dy = x[1] - 0.5 | |
16 | + values[0] = 500.0*exp(-(dx*dx + dy*dy)/0.02) | |
17 | + | |
18 | +# Neumann boundary condition | |
19 | +class Flux(Function): | |
20 | + def __init__(self, element, mesh): | |
21 | + Function.__init__(self, element, mesh) | |
22 | + def eval(self, values, x): | |
23 | + if x[0] > DOLFIN_EPS: | |
24 | + values[0] = 25.0*sin(5.0*DOLFIN_PI*x[1]) | |
25 | + else: | |
26 | + values[0] = 0.0 | |
27 | + | |
28 | +# Sub domain for Dirichlet boundary condition | |
29 | +class DirichletBoundary(SubDomain): | |
30 | + def inside(self, x, on_boundary): | |
31 | + return bool(on_boundary and x[0] < DOLFIN_EPS) | |
32 | + | |
33 | +# Define variational problem | |
34 | +v = TestFunction(element) | |
35 | +u = TrialFunction(element) | |
36 | +f = Source(element, mesh) | |
37 | +g = Flux(element, mesh) | |
38 | + | |
39 | +a = dot(grad(v), grad(u))*dx | |
40 | +L = v*f*dx + v*g*ds | |
41 | + | |
42 | +# Define boundary condition | |
43 | +u0 = Function(mesh, 0.0) | |
44 | +boundary = DirichletBoundary() | |
45 | +bc = DirichletBC(u0, mesh, boundary) | |
46 | + | |
47 | +# Solve PDE and plot solution | |
48 | +pde = LinearPDE(a, L, mesh, bc) | |
49 | +u = pde.solve() | |
50 | +plot(u) | |
51 | + | |
52 | +# Save solution to file | |
53 | +file = File("poisson.pvd") | |
54 | +file << u |