dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Very nice. Can you add the __author__, __date__ strings etc as in the
other demos and then send it as an hg bundle (see appendix in manual)?
Put it under

  src/demo/pde/dielectric/python/demo.py

Also update the AUTHORS file.

-- 
Anders


On Thu, Feb 14, 2008 at 02:21:15PM +0100, Kristen Kaasbjerg wrote:
> 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)
> 
> 
> 
> 
> 
> 
> 

> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev