dolfin team mailing list archive

["Garth N. Wells" <gnw20@xxxxxxxxx>]

On 29 January 2013 08:12, Johan Hake <hake.dev@xxxxxxxxx> wrote:
> 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.
>

It is not an onerous burden to document a demo when making a
contribution.  We want to document the demos, but we go backwards if
undocumented demos are added faster than documented demos. Compared to
other projects, our documentation of demo programs is weak and should
be improved.

Garth

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