← Back to branch summary

~ubuntu-branches/ubuntu/maverick/dolfin/maverick

« back to all changes in this revision

Viewing changes to demo/pde/dielectric/python/demo.py

  • Committer: Bazaar Package Importer
  • Author(s): Johannes Ring
  • Date: 2008-09-16 08:41:20 UTC
  • Revision ID: james.westby@ubuntu.com-20080916084120-i8k3u6lhx3mw3py3
Tags: upstream-0.9.2
Import upstream version 0.9.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/pde/dielectric/python/demo.py
 
1
# This demo program calculates the electrostatic potential
 
2
# in the unit square by solving Laplace's equation
 
3
#
 
4
#      div \epsilon_r grad u(x, y) = 0
 
5
#
 
6
# The lower half is filled with a dielectric material
 
7
# with dielectric constant \epsilon_r.
 
8
#
 
9
# The boundary conditions are:
 
10
#
 
11
#     u(x, y)     = 0  for y > 0
 
12
#     u(x,y)      = V  for y = 0
 
13
 
 
14
__author__ = "Kristen Kaasbjerg (cosby@fys.ku.dk)"
 
15
__date__ = "2008-02-14 -- 2008-12-13"
 
16
__copyright__ = ""
 
17
__license__  = "GNU LGPL Version 2.1"
 
18
 
 
19
# Modified by Kristian Oelgaard 2008
 
20
 
 
21
from dolfin import *
 
22
 
 
23
l = 1.0; h = 1.0 # unit square 
 
24
h_ = 0.5         # position of the dielectric interface
 
25
e_r = 10         # dielectric constant for y<h_ (e_r=1 for y>h_
 
26
V = 1.0          # applied voltage at the y=0 boundary
 
27
 
 
28
# Create mesh and finite element
 
29
mesh = Mesh()
 
30
editor = MeshEditor()
 
31
editor.open(mesh, "triangle", 2, 2) 
 
32
editor.initVertices(4)
 
33
editor.initCells(2)
 
34
editor.addVertex(0, 0.0, 0.0)
 
35
editor.addVertex(1, 0.0, h)
 
36
editor.addVertex(2, l, h)
 
37
editor.addVertex(3, l,0)
 
38
editor.addCell(0, 0, 1, 2)
 
39
editor.addCell(1, 0, 2, 3)
 
40
editor.close()
 
41
 
 
42
# Refine mesh
 
43
for refine in range(4):
 
44
    mesh.refine()
 
45
 
 
46
# Function spaces
 
47
V2 = FunctionSpace(mesh, "CG", 2) # last argument: order of the polynomial
 
48
 
 
49
# Discontinuous element needed for the dielectric function
 
50
V0 = FunctionSpace(mesh, "DG", 0) # 0'th order are possible for DG elements
 
51
 
 
52
# Exact solution
 
53
class Exact(Function):
 
54
    def eval(self, values, x):
 
55
        u = 0
 
56
        N = 20 # N needs to be rather large (~200) in other to have u=V for y=0 !
 
57
        pi = DOLFIN_PI
 
58
        if x[1]<=h_:
 
59
            for n in range(N):
 
60
                n_ = 2*n+1 #n>0!
 
61
                X = (1.-exp(2*n_*pi*(1-h_)))/(1+exp(2*n_*pi*(1-h_)))
 
62
                u -= sin(n_*pi*x[0])*4*V/(pi*n_) * ( (1+e_r*X)*exp(n_*pi*(x[1]-h_))+
 
63
                                                     (-1.+e_r*X)*exp(-n_*pi*(x[1]-h_)) ) /\
 
64
                                                     ( (1.-e_r*X)*exp(n_*pi*h_) -
 
65
                                                       (1.+e_r*X)*exp(-n_*pi*h_))
 
66
        else:
 
67
            for n in range(N):
 
68
                n_ = 2*n+1 #n>0!
 
69
                X = (1.-exp(2*n_*pi*(1-h_)))/(1+exp(2*n_*pi*(1-h_)))
 
70
                X1 = 1./(1+exp(2*n_*pi*(1-h_)))
 
71
                X2 = exp(2*n_*pi*(1-h_))/(1+exp(2*n_*pi*(1-h_)))
 
72
                u += sin(n_*pi*x[0])*4*V/(pi*n_) * ( 2*e_r*X1*exp(n_*pi*(x[1]-h_)) -
 
73
                                                     2*e_r*X2*exp(-n_*pi*(x[1]-h_))) /\
 
74
                                                     ( (1.+e_r*X)*exp(-n_*pi*h_) -
 
75
                                                       (1.-e_r*X)*exp(n_*pi*h_) )
 
76
        values[0] = u 
 
77
        
 
78
# Dielectric constant
 
79
class Coefficient(Function):
 
80
    def eval(self, values, x):
 
81
        if x[1] <= h_:
 
82
            values[0] = e_r
 
83
        else:
 
84
            values[0] = 1.0
 
85
 
 
86
# Dirichlet boundary condition
 
87
class DirichletFunction(Function):
 
88
    def eval(self, values, x):
 
89
        if x[1] < DOLFIN_EPS:
 
90
            values[0] = V
 
91
        else:
 
92
            values[0] = 0.0
 
93
            
 
94
# Sub domain for Dirichlet boundary condition
 
95
class DirichletBoundary(SubDomain):
 
96
    def inside(self, x, on_boundary):
 
97
        return on_boundary and True
 
98
 
 
99
# Define variational problem
 
100
v = TestFunction(V2)
 
101
u = TrialFunction(V2)
 
102
f = Constant(mesh, 0)
 
103
b = Coefficient(V0)
 
104
a = b*dot(grad(v), grad(u))*dx
 
105
L = v*f*dx
 
106
 
 
107
# Define boundary condition
 
108
u0 = DirichletFunction(V2)
 
109
bc = DirichletBC(V2, u0, DirichletBoundary())
 
110
 
 
111
# Solve PDE and plot solution
 
112
problem = VariationalProblem(a, L, bc)
 
113
u = problem.solve()
 
114
 
 
115
plot(u)
 
116
 
 
117
# Calculate difference between exact and FEM solution
 
118
# Use higher order element for exact solution,
 
119
# because it is an interpolation of the exact solution
 
120
# in the finite element space!
 
121
Pk = FunctionSpace(mesh, "CG", 5)
 
122
exact = Exact(Pk)
 
123
 
 
124
e = u - exact
 
125
L2_norm = e*e*dx
 
126
norm = sqrt(assemble(L2_norm))
 
127
 
 
128
print "L2-norm of error is: ", norm