On 29 January 2013 06:55, Anders Logg <l...@simula.no> wrote:
> 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).
>
And it should be documented along the lines of, e.g.
http://fenicsproject.org/documentation/dolfin/1.1.0/python/demo/pde/poisson/python/documentation.html
Garth
> --
> 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@lists.launchpad.net
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help : https://help.launchpad.net/ListHelp
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : dolfin@lists.launchpad.net
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
