On 01/29/2013 08:57 AM, Garth N. Wells wrote: > On 29 January 2013 06:55, Anders Logg <[email protected]> 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 : [email protected] >>> Unsubscribe : https://launchpad.net/~dolfin >>> More help : https://help.launchpad.net/ListHelp >> >> _______________________________________________ >> Mailing list: https://launchpad.net/~dolfin >> Post to : [email protected] >> Unsubscribe : https://launchpad.net/~dolfin >> More help : https://help.launchpad.net/ListHelp > > _______________________________________________ > Mailing list: https://launchpad.net/~dolfin > Post to : [email protected] > Unsubscribe : https://launchpad.net/~dolfin > More help : https://help.launchpad.net/ListHelp > _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
