SyFi  0.3
demo.py
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines