dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Thanks for the example.

This would be a nice example to include in DOLFIN. The following steps
would be necessary in order to do this:

1. It needs to follow the style of other demos wrt copyright notice
etc.

2. We need a corresponding C++ version which should look as much as
the Python demo as possible.

3. It needs more commenting (look at other examples).

4. You need to sign and send the licensing statement (see web page).

--
Anders


On Tue, Jan 29, 2013 at 01:37:10AM +0100, Pedro Guarderas wrote:
> 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()
> #----------------------------------------------------------------
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp