dolfin team mailing list archive

[Johan Hake <hake.dev@xxxxxxxxx>]

On 01/29/2013 08:57 AM, Garth N. Wells wrote:
> On 29 January 2013 06:55, Anders Logg <logg@xxxxxxxxx> 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

Is the documentation policy required for all new demos added? That is a
pretty strong criteria.

I agree with all of the bullet points Anders wrote. I also think that
the documentation is _very_ nice and should be the gold standard, but I
am not sure it is wise to reject demos without it.

Johan

> 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@xxxxxxxxxxxxxxxxxxx
>>> Unsubscribe : https://launchpad.net/~dolfin
>>> More help   : https://help.launchpad.net/ListHelp
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~dolfin
>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help   : https://help.launchpad.net/ListHelp
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>