SyFi
0.3
|
00001 """This demo program solves Poisson's equation 00002 00003 - div grad u(x, y) = f(x, y) 00004 00005 on the unit square with source f given by 00006 00007 f(x, y) = 500*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02) 00008 00009 and boundary conditions given by 00010 00011 u(x, y) = 0 for x = 0 or x = 1 00012 """ 00013 00014 # Copyright (C) 2007-2008 Anders Logg 00015 # 00016 # This file is part of SyFi. 00017 # 00018 # SyFi is free software: you can redistribute it and/or modify 00019 # it under the terms of the GNU General Public License as published by 00020 # the Free Software Foundation, either version 2 of the License, or 00021 # (at your option) any later version. 00022 # 00023 # SyFi is distributed in the hope that it will be useful, 00024 # but WITHOUT ANY WARRANTY; without even the implied warranty of 00025 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00026 # GNU General Public License for more details. 00027 # 00028 # You should have received a copy of the GNU General Public License 00029 # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 00030 # 00031 # Modified to work with SFC by Kent-Andre Mardal 00032 # 00033 # First added: 2007-08-16 00034 # Last changed: 2008-12-13 00035 00036 import sys 00037 try: 00038 from dolfin import * 00039 parameters["form_compiler"]["name"] = "sfc" 00040 except: 00041 print "Dolfin is not found. Demo can not run" 00042 sys.exit(2) 00043 00044 # Create mesh and define function space 00045 mesh = UnitSquare(32, 32) 00046 V = FunctionSpace(mesh, "CG", 1) 00047 00048 # Define Dirichlet boundary (x = 0 or x = 1) 00049 class DirichletBoundary(SubDomain): 00050 def inside(self, x, on_boundary): 00051 return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS 00052 00053 # Define boundary condition 00054 u0 = Constant(0.0) 00055 bc = DirichletBC(V, u0, DirichletBoundary()) 00056 00057 # Define variational problem 00058 v = TestFunction(V) 00059 u = TrialFunction(V) 00060 f = Expression("500.0 * exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)") 00061 ff = Function(V) 00062 ff.interpolate(f) 00063 a = dot(grad(v), grad(u))*dx 00064 L = v*ff*dx 00065 00066 # Compute solution 00067 u = Function(V) 00068 solve(a == L, u, bc) 00069 00070 ok = (u.vector().norm("l2") - 142.420764968) < 10**-4 00071 00072 sys.exit(0 if ok else 1) 00073