# PDE Homogenization

Registered by Pedro Guarderas

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()

## Blueprint information

Status:
Not started
Approver:
None
Priority:
Undefined
Drafter:
None
Direction:
Needs approval
Assignee:
Definition:
Review
Series goal:
None
Implementation:
Unknown
Milestone target:
None

(?)

### Work Items

This blueprint contains Public information
Everyone can see this information.

No subscribers.