PDE Homogenization

I was cheking the dolfin examples and I saw that there is no example of PDE
Homogenization, therefore I decided to construct an 2D example, try it please.

Let me know if everything is fine or if we can improve the example, expecially
for the construction of the "fast" functions w(x/e).

Here the source:

#----------------------------------------------------------------
# PDE Homogenization
#
# author: Pedro Guarderas
# date: 29/01/2013
#
#----------------------------------------------------------------

from dolfin import *
from math import floor

#----------------------------------------------------------------
N = 5.0
M = 5.0
alpha = 1.0
beta = 0.01
g = -9.8

class BExp( Expression ):
 def eval( self, values, x ):
  values[0] = 0.0

class FExp( Expression ):
 def eval( self, values, x ):
  values[0] = g * x[1]

class AExp( Expression ):
 def eval( self, values, x ):
  values[0] = alpha * floor( 2 * x[0] * N + 0.5 ) % 2
  values[1] = beta
  values[2] = beta
  values[3] = alpha * floor( 2 * x[1] * N + 0.5 ) % 2

 def value_shape(self):
         return (2,2)

#----------------------------------------------------------------
# Problem
mesh = UnitSquareMesh( 50, 50 )
V = FunctionSpace( mesh, "Lagrange", 2 )

u = TrialFunction( V )
v = TestFunction( V )
f = FExp()

A = AExp()

u0 = BExp()
def u0_boundary( x, on_boundary ):
 return on_boundary

bc = DirichletBC( V, u0, u0_boundary )

a = inner( A * grad( u ), grad( v ) ) * dx
L = f * v * dx

u = Function( V )
solve( a == L, u, bc )

plot( u, wireframe = True,
 title = 'Solution',
 axes = True )

fileH = File( "solution.pvd" )
fileH << u

#----------------------------------------------------------------
# Homogenization method
e = 0.05

class AExps( Expression ):
 def eval( self, values, x ):
  values[0] = alpha * floor( 2 * x[0] + 0.5 ) % 2
  values[1] = beta
  values[2] = beta
  values[3] = alpha * floor( 2 * x[1] + 0.5 ) % 2

 def value_shape(self):
         return (2,2)

class AExpc( Expression ):
 def eval( self, values, x ):
  values[0] = alpha * floor( 2 * x[0] / e + 0.5 ) % 2
  values[1] = beta
  values[2] = beta
  values[3] = alpha * floor( 2 * x[1] / e + 0.5 ) % 2

 def value_shape(self):
         return (2,2)

# Cell problems
# Tricky part: fast and slow solutions
e1 = Constant( (1.0, 0.0) )
e2 = Constant( (0.0, 1.0) )

Acs = AExps()
Ac = AExpc()

meshc = UnitSquareMesh( 50, 50 )
W = FunctionSpace( meshc, "Lagrange", 2 )

# first cell
w1 = TrialFunction( W ) # fast w1(x) = w1s(x/e)
w1s = TrialFunction( W ) # slow
v1 = TestFunction( W )
b1 = DirichletBC( W, u0, u0_boundary )

a1 = e * inner( Ac * grad( w1 ), grad( v1 ) ) * dx
L1 = -inner( Ac * e1, grad(v1) ) * dx

a1s = inner( Acs * grad( w1s ), grad( v1 ) ) * dx
L1s = -inner( Acs * e1, grad(v1) ) * dx

w1 = Function( W )
w1s = Function( W )
solve( a1 == L1, w1, b1 )
solve( a1s == L1s, w1s, b1 )

plot( w1s, wireframe = True,
 title = "cell 1",
 axes = True )

file1 = File( "cell_solution_1.pvd" )
file1 << w1

# second cell
w2 = TrialFunction( W ) # fast w2(x) = w2s(x/e)
w2s = TrialFunction( W ) # slow
v2 = TestFunction( W )
b2 = DirichletBC( W, u0, u0_boundary )

a2 = e * inner( Ac * grad( w2 ), grad( v2 ) ) * dx
L2 = -inner( Ac * e2, grad(v2) ) * dx

a2s = inner( Acs * grad( w2s ), grad( v2 ) ) * dx
L2s = -inner( Acs * e2, grad(v2) ) * dx

w2 = Function( W )
w2s = Function( W )
solve( a2 == L2, w2, b2 )
solve( a2s == L2s, w2s, b2 )

plot( w2s, wireframe = True,
 title = "cell 2",
 axes = True )

file2 = File( "cell_solution_2.pvd" )
file2 << w2

#----------------------------------------------------------------
# Homogenized equation
# u0

# Determine A* homogenization of A
# here, we use the slow functions
A11 = inner( A * ( e1 + grad(w1s) ), e1 ) * dx
A12 = inner( A * ( e2 + grad(w2s) ), e1 ) * dx
A21 = inner( A * ( e1 + grad(w1s) ), e2 ) * dx
A22 = inner( A * ( e2 + grad(w2s) ), e1 ) * dx

a11 = assemble(A11)
a12 = assemble(A12)
a21 = assemble(A21)
a22 = assemble(A22)

print "a list = ", a11, a12, a21, a22

class AExph( Expression ):
 def eval( self, values, x ):
  values[0] = a11
  values[1] = a12
  values[2] = a21
  values[3] = a22

 def value_shape(self):
         return (2,2)

Ah = AExph()

uh = TrialFunction( W )
vh = TestFunction( W )
bh = DirichletBC( W, u0, u0_boundary )

ah = inner( Ah * grad( uh ), grad( vh ) ) * dx
Lh = f * vh * dx

uh = Function( W )
solve( ah == Lh, uh, bh )

plot( uh, wireframe = True,
 title = 'Homogenization u_0',
 axes = True )

fileh = File( "homgenization.pvd" )
fileh << uh

#----------------------------------------------------------------
# Homogenization, approximation of first order
# ue = u0 + e * inner( gra(u_0), [w1,w2]' )
vhe = TestFunction(W)
uhe = TrialFunction(W)
bhe = DirichletBC( W, u0, u0_boundary )

ae = uhe * vhe * dx
Le = ( uh + e * ( inner( grad( uh ), e1 ) * w1
 + inner( grad( uh ), e2 ) * w2 ) ) * vhe * dx

uhe = Function(W)
solve( ae == Le, uhe, bhe )

plot( uhe, wireframe = True,
 title = 'Homogenization u_e',
 axes = True )

filehe = File( "homgenization_e.pvd" )
filehe << uhe

interactive()