dolfin team mailing list archive

[Pedro Guarderas <namstdpfgo@xxxxxxxxx>]

Hello

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

Pedro


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