dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

Hello!

I have added a python version of the discontinuous galerkin advection 
diffusion demo. I used the newly added functionality to compile cpp function 
on the fly. As I only more or less copied someone elses code I am not sure I 
got the copyright message right.

Please have a look!

Johan

# HG changeset patch
# User "Johan Hake <hake@xxxxxxxxx>"
# Date 1213866469 -7200
# Node ID a9bc99897bb6deb200316c5c1985ac9414f0b4e4
# Parent  4bddacadb8b7e14138a52d223c37a6c17427c8be
Added python version of advection diffusion demo

diff -r 4bddacadb8b7 -r a9bc99897bb6 demo/pde/dg/advection_diffusion/python/README
--- a/demo/pde/dg/advection_diffusion/python/README	Wed Jun 18 19:30:31 2008 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,2 +0,0 @@
-A Python version of this demo is currently missing.
-Please consider contributing the missing code.
diff -r 4bddacadb8b7 -r a9bc99897bb6 demo/pde/dg/advection_diffusion/python/demo.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demo/pde/dg/advection_diffusion/python/demo.py	Thu Jun 19 11:07:49 2008 +0200
@@ -0,0 +1,116 @@
+""" Steady state advection-diffusion equation,
+discontinuous formulation using full upwinding.
+
+Constant velocity field with homogeneous Dirichlet boundary conditions
+on all boundaries.
+
+Implemented in python from cpp demo by Johan Hake
+
+"""
+
+__author__ = "Johan Hake (hake@xxxxxxxxx)"
+__date__ = "2008-06-19"
+__copyright__ = "Copyright (C) 2008 Johan Hake"
+__license__  = "GNU LGPL Version 2.1"
+
+from dolfin import *
+
+def forms(scalar, vector, constant):
+    # Import local ffc to be able to define forms
+    # Note: This is needed because the compiled functions that is used below are
+    # only cpp_Functions not ffc_Functions, and cannot be used to define forms
+    import ffc
+    
+    # Test and trial functions
+    v  = TestFunction(scalar)
+    u  = TrialFunction(scalar)
+    
+    # Coeffisient functions
+    f  = ffc.Function(scalar)
+    b  = ffc.Function(vector)    # Note: b('+') == b('-')
+    of = ffc.Function(constant)  # This function is a switch that determines if a
+                                 # facet is an outflow facet or not 1.0 or 0.0
+    n  = ffc.FacetNormal("triangle")
+    h  = ffc.MeshSize("triangle")
+                            
+    kappa = Constant("triangle")
+    alpha = Constant("triangle")
+    
+    def upwind(u, b):
+        return [b[i]('+')*(of('+')*u('+') + of('-')*u('-')) for i in range(len(b))]
+    
+    # Bilinear form
+    a_int = dot( grad(v), mult(kappa, grad(u)) - mult(b,u) )*dx
+    
+    a_fac = kappa('+')*alpha('+')/h('+')*dot(jump(v, n), jump(u, n))*dS \
+            - kappa('+')*dot(avg(grad(v)), jump(u, n))*dS \
+            - kappa('+')*dot(jump(v, n), avg(grad(u)))*dS
+    
+    a_gd = kappa*alpha/h*v*u*ds \
+           - kappa*dot(grad(v), mult(u,n))*ds \
+           - kappa*dot(mult(v,n), grad(u))*ds
+    
+    a_vel = dot(jump(v, n), upwind(u, b))*ffc.dS + dot(mult(v, n), mult(b, of*u))*ds
+    
+    a = a_int + a_fac + a_gd + a_vel
+    
+    # Linear form
+    L = v*f*dx
+    
+    return a, L
+
+mesh =  UnitSquare (64, 64)
+
+# Defining the finite element room
+scalar_DG = FiniteElement("Discontinuous Lagrange", "triangle", 2)
+scalar_CG = FiniteElement("Lagrange", "triangle", 2)
+vector_CG = VectorElement("Lagrange", "triangle", 2)
+constant  = FiniteElement("Discontinuous Lagrange", "triangle", 0)
+
+# Defining dolfin coeffisient functions
+file_string = open('functions2D.h').read()
+coeffisients = compile_functions(file_string,mesh)
+
+# Extracting the compiled functions
+source   = coeffisients.pop(0)
+velocity = coeffisients[0]
+outflow_facet = coeffisients[1]
+
+# Setting members of the compiled functions
+outflow_facet.velocity = velocity
+source.C               = 0.0
+
+coeffisients.append(FacetNormal("triangle", mesh))
+coeffisients.append(AvgMeshSize("triangle", mesh)) 
+coeffisients.append(Function(constant,mesh,0.0))   # The diffusivity
+coeffisients.append(Function(constant,mesh,20.0))  # The penalty term
+
+a, L = forms(scalar_DG, vector_CG, constant)
+
+# Assembly and solve system
+A = assemble(a, mesh, coeffisients)
+b = assemble(L, mesh, [source])
+
+uh = Function(scalar_DG, mesh, Vector())
+
+solve(A, uh.vector(), b)
+
+# Defining projection forms 
+vp = TestFunction(scalar_CG)
+up = TrialFunction(scalar_CG)
+
+u0 = Function(scalar_DG)
+
+ap = dot(vp,up)*dx
+Lp = dot(vp,uh)*dx
+
+# Define and solve PDE
+pde = LinearPDE(ap, Lp, mesh)
+
+up = pde.solve()
+
+file = File("temperature.pvd")
+file << up
+
+# Plot solution
+plot(up, interactive=True)
diff -r 4bddacadb8b7 -r a9bc99897bb6 demo/pde/dg/advection_diffusion/python/functions2D.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demo/pde/dg/advection_diffusion/python/functions2D.h	Thu Jun 19 11:07:49 2008 +0200
@@ -0,0 +1,87 @@
+// Source term
+class Source2D : public dolfin::Function
+{
+public:
+    
+  Source2D(dolfin::Mesh& mesh) : dolfin::Function(mesh) 
+  {
+    C = 1.0;
+  }
+
+  real eval(const real* x) const
+  {
+
+    real vx = -exp(x[0])*(x[1]*cos(x[1]) + sin(x[1]));
+    real vy = exp(x[0])*(x[1]*sin(x[1]));
+
+    real ux = 5.0*DOLFIN_PI*cos(5.0*DOLFIN_PI*x[0])*sin(5.0*DOLFIN_PI*x[1]);
+    real uy = 5.0*DOLFIN_PI*sin(5.0*DOLFIN_PI*x[0])*cos(5.0*DOLFIN_PI*x[1]);
+    real uxx = -25.0*DOLFIN_PI*DOLFIN_PI*sin(5.0*DOLFIN_PI*x[0])*sin(5.0*DOLFIN_PI*x[1]);
+    real uyy = -25.0*DOLFIN_PI*DOLFIN_PI*sin(5.0*DOLFIN_PI*x[0])*sin(5.0*DOLFIN_PI*x[1]);
+
+    return vx*ux + vy*uy - C*(uxx + uyy);
+  }
+  dolfin::real C;
+};
+
+// Velocity
+class Velocity2D : public dolfin::Function
+{
+public:
+    
+  Velocity2D(dolfin::Mesh& mesh) : dolfin::Function(mesh) {}
+
+  void eval(real* values, const real* x) const
+  {
+    values[0] = -exp(x[0])*(x[1]*cos(x[1]) + sin(x[1]));
+    values[1] = exp(x[0])*(x[1]*sin(x[1]));
+  }
+};
+
+class OutflowFacet2D : public dolfin::Function
+{
+public:
+
+  OutflowFacet2D(dolfin::Mesh& mesh) : dolfin::Function(mesh)
+  {
+    velocity = 0;
+  }
+
+  real eval(const real* x) const
+  {
+    // If there is no facet (assembling on interior), return 0.0
+    if (facet() < 0)
+      return 0.0;
+    else
+    {
+      if (!velocity)
+	error("Attach a velocity function.");
+      real normal_vector[2];
+      real velocities[2] = {0.0, 0.0};
+      
+      // Compute facet normal
+      for (dolfin::uint i = 0; i < cell().dim(); i++)
+        normal_vector[i] = cell().normal(facet(), i);
+      
+      // Get velocities
+      velocity->eval(velocities, x);
+      
+      // Compute dot product of the facet outward normal and the velocity vector
+      real dot = 0.0;
+      for (dolfin::uint i = 0; i < cell().dim(); i++)
+        dot += normal_vector[i]*velocities[i];
+      
+      // If dot product is positive the facet is an outflow facet, 
+      // meaning the contribution from this cell is on the upwind side.
+      if (dot > DOLFIN_EPS)
+      {
+        return 1.0;
+      }
+      else
+      {
+        return 0.0;
+      }
+    }
+  }
+  dolfin::Function* velocity;
+};