dolfin team mailing list archive

[Kristen Kaasbjerg <cosby@xxxxxxxxx>]

Hi,

Here is an example of a solution to Laplace equation
with a discontinuous dielectric function.
The obtained solution is compared to the exact solution
by calculating the L2-norm of the difference.

Thought you might want to include this in the demo dir
sometime ?

Kristen

# This demo program finds the electrostatic potential
# in the unit square by solving Laplace's equation
#
#      div \epsilon_r grad u(x, y) = 0
#
# The lower half is filled with a dielectric material
# with dielectric constant \epsilon_r.
#
# The boundary conditions are:
#
#     u(x, y)     = 0  for y > 0
#     u(x,y)      = V  for y = 0

from dolfin import *

l = 1. ; h = 1. ; h_ = .5 ; e_r = 10 ; V = 1.

# Create mesh and finite element
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh,"triangle", 2,2) # topological_dimension (cells) and geometric_dimension (vertices)
editor.initVertices(4)
editor.initCells(2)
editor.addVertex(0, 0.0, 0.0)
editor.addVertex(1, 0.0, h)
editor.addVertex(2, l, h)
editor.addVertex(3, l,0)
editor.addCell(0, 0, 1, 2)
editor.addCell(1, 0, 2, 3)
editor.close()

for refine in range(4):
    mesh.refine()

element = FiniteElement("Lagrange", "triangle", 2) # last argument: order of the polynomial
DG0 = FiniteElement("DG", "triangle", 0) # 0'th order are possible for DG elements

# Source term
class Source(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        values[0] = 0.0

# Exact solution
class Exact(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        phi = 0
        N = 20
        pi = DOLFIN_PI
        if x[1]<=h_:
            for n in range(N):
                n_ = 2*n+1 #n>0!
                X = (1.-exp(2*n_*pi*(1-h_)))/(1+exp(2*n_*pi*(1-h_)))
                phi -= sin(n_*pi*x[0])*4*V/(pi*n_) * ( (1+e_r*X)*exp(n_*pi*(x[1]-h_))+
                                                       (-1.+e_r*X)*exp(-n_*pi*(x[1]-h_)) ) /\
                                                       ( (1.-e_r*X)*exp(n_*pi*h_) -
                                                         (1.+e_r*X)*exp(-n_*pi*h_))
        else:
            for n in range(N):
                n_ = 2*n+1 #n>0!
                X = (1.-exp(2*n_*pi*(1-h_)))/(1+exp(2*n_*pi*(1-h_)))
                X1 = 1./(1+exp(2*n_*pi*(1-h_)))
                X2 = exp(2*n_*pi*(1-h_))/(1+exp(2*n_*pi*(1-h_)))
                phi += sin(n_*pi*x[0])*4*V/(pi*n_) * ( 2*e_r*X1*exp(n_*pi*(x[1]-h_)) -
                                                    2*e_r*X2*exp(-n_*pi*(x[1]-h_))) /\
                                                    ( (1.+e_r*X)*exp(-n_*pi*h_) - (1.-e_r*X)*exp(n_*pi*h_) )
        values[0] = phi 
        
# Coefficient
class Coefficient(Function):
    def __init__(self, element, mesh):
        Function.__init__(self, element, mesh)
    def eval(self, values, x):
        if x[1] <= h_:
            values[0] = e_r
        else:
            values[0] = 1.0

# Dirichlet boundary condition
class DirichletFunction(Function):
    def __init__(self, element,mesh):
        Function.__init__(self, element,mesh)
    def eval(self, values, x):
        if x[1] < DOLFIN_EPS:
            values[0] = V
        else:
            values[0] = 0.0
            
# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return bool(on_boundary and True)

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

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

# Define boundary condition
u0 = DirichletFunction(element,mesh)
boundary = DirichletBoundary()
bc = DirichletBC(u0, mesh, boundary)

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

plot(u)

#Calculate difference between exact and FEM solution
# Use higher order element for exact solution,
# because it is an interpolation of the exact solution
# in the finite element space!
Pk = FiniteElement("Lagrange", "triangle", 5)
exact = Exact(Pk,mesh) 

e = u - exact
L2_norm = e*e*dx
norm = assemble(L2_norm,mesh)

print "L2-norm is: ", norm
#plot(exact)