"""
This demo program illustrates how to solve Poisson's equation
- div grad u(x, y) = f(x, y)
on the unit square with pure Neumann boundary conditions:
du/dn(x, y) = -sin(5*x)
and source f given by
f(x, y) = 10*exp(-((x - 0.5)^2 + (y - 0.5)^2) / 0.02)
Since only Neumann conditions are applied, u is only determined up to
a constant c by the above equations. An addition constraint is thus
required, for instance
\int u = 0
This can be accomplished by introducing the constant c as an
additional unknown (to be sought in the space of real numbers)
and the above constraint.
"""
# Copyright (C) 2010 Marie E. Rognes
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see .
#
# Modified by Anders Logg 2011
#
# First added: 2010-05-10
# Last changed: 2012-11-12
# Begin demo
from dolfin import *
# Create mesh and define function space
mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R
# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("-sin(5*x[0])")
a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = f*v*dx + g*v*ds
# Compute solution
w = Function(W)
solve(a == L, w)
(u, c) = w.split()
# Plot solution
plot(u, interactive=True)