dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Feb 12, 2008 at 03:42:48PM +0100, cosby@xxxxxxxxx wrote:

> Hi again
> 
> Sorry for having been absent from the discussion, but for some
> reason that I'm still to figure out, I only receive the daily
> Dolfin digest mails.
> So what conclusion did you reach ? Can one just go ahead and
> multiply a discontinuous coefficient to the variational form:
>      \epsilon * grad(v)*grad(u) = 0
> 
> Or will there be an additional term from the boundary integral
> along the interface where \epsilon is discontinuous ?
> Anders, could I ask you to send me your modified poisson script
> where you had added the discontinuous coefficients ?
> Cheers
> Kristen

Sure, here it is.

-- 
Anders

# This demo program solves Poisson's equation
#
#     - div grad u(x, y) = f(x, y)
#
# on the unit square with source f given by
#
#     f(x, y) = 500*exp(-((x-0.5)^2 + (y-0.5)^2)/0.02)
#
# and boundary conditions given by
#
#     u(x, y)     = 0  for x = 0
#     du/dn(x, y) = 1  for x = 1
#     du/dn(x, y) = 0  otherwise

__author__ = "Anders Logg (logg@xxxxxxxxx)"
__date__ = "2007-08-16 -- 2007-08-21"
__copyright__ = "Copyright (C) 2007 Anders Logg"
__license__  = "GNU LGPL Version 2.1"

from dolfin import *

# Create mesh and finite element
mesh = UnitSquare(32, 32)
element = FiniteElement("Lagrange", "triangle", 1)
DG0 = FiniteElement("DG", "triangle", 0)

# Source term
class Source(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        dx = x[0] - 0.5
        dy = x[1] - 0.5
        values[0] = 500.0*exp(-(dx*dx + dy*dy)/0.02)

# Neumann boundary condition
class Flux(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        if x[0] > DOLFIN_EPS:
            values[0] = 25.0*sin(5.0*DOLFIN_PI*x[1])
        else:
            values[0] = 0.0

# Coefficient
class Coefficient(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        if x[1] > 0.5:
            values[0] = 1.0
        else:
            values[0] = 0.1

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(on_boundary and x[0] < DOLFIN_EPS)

# Define variational problem
v = TestFunction(element)
u = TrialFunction(element)
f = Source(element, mesh)
g = Flux(element, mesh)
b = Coefficient(DG0, mesh)

a = b*dot(grad(v), grad(u))*dx
L = v*f*dx + v*g*ds

# Define boundary condition
u0 = Function(mesh, 0.0)
boundary = DirichletBoundary()
bc = DirichletBC(u0, mesh, boundary)

# Solve PDE and plot solution
pde = LinearPDE(a, L, mesh, bc)
u = pde.solve()
plot(u)

# Save solution to file
file = File("poisson.pvd")
file << u