On Thu, Mar 12, 2009 at 10:25:36AM +0100, Martin Sandve Alnæs wrote: > I have serious problems with the idea of letting user-defined > functions change nature to discrete functions as a side effect > of other operations, effectively hiding the user-defined implementation > of eval. > > (I know this concept wasn't introduced in this discussion, but it's related).
You mean by calling u.vector()? Yes, I agree it's problematic. But I don't know how to handle it otherwise. Consider the following example: u = Function(V) solve(A, u.vector(), b) First u is a user-defined function since it doesn't have a vector. Then it becomes discrete when we ask for the vector. The problem I think is that there is no way for us to check whether a user has overloaded eval() without trying to call it. -- Anders > Martin > > > > On Thu, Mar 12, 2009 at 10:00 AM, Garth N. Wells <[email protected]> wrote: > > > > > > Garth N. Wells wrote: > >> > >> Anders Logg wrote: > >>> On Wed, Mar 11, 2009 at 11:15:44PM +0100, Anders Logg wrote: > >>>> On Wed, Mar 11, 2009 at 11:01:56PM +0100, DOLFIN wrote: > >>>>> One or more new changesets pushed to the primary dolfin repository. > >>>>> A short summary of the last three changesets is included below. > >>>>> > >>>>> changeset: 5853:f1ef6132a568d5a56e5c70b17ce118c19bfa961c > >>>>> tag: tip > >>>>> user: Anders Logg <[email protected]> > >>>>> date: Wed Mar 11 23:01:49 2009 +0100 > >>>>> files: ChangeLog demo/pde/poisson/cpp/Poisson.h > >>>>> dolfin/function/Function.cpp sandbox/misc/Poisson.form > >>>>> sandbox/misc/Poisson.h sandbox/misc/README sandbox/misc/SConstruct > >>>>> sandbox/misc/cpp/Poisson.form sandbox/misc/cpp/Poisson.h > >>>>> sandbox/misc/cpp/SConstruct sandbox/misc/cpp/main.cpp > >>>>> sandbox/misc/main.cpp > >>>>> description: > >>>>> Automatically interpolate user-defined functions on assignment > >>>> This is something we discussed a while back but we didn't agree on > >>>> whether or not it was a good idea. I think it is and it's easy enough > >>>> to remove if there is strong enough pressure against it. > >>>> > >>>> Here are two examples of assignment: > >>>> > >>>> Case 1: Time-stepping with user-defined initial data > >>>> > >>>> # Initializations > >>>> mesh = UnitSquare(32, 32) > >>>> V = FunctionSpace(mesh, "Lagrange", 1) > >>>> u0 = Function(V, "sin(x[0])") > >>>> u1 = Function(V) > >>>> > >>>> # Time stepping > >>>> for i in range(10): > >>>> > >>>> print i > >>>> > >>>> # Solve for u1 > >>>> u1.vector() > >>>> > >>>> # Assign u0 = u1 > >>>> u0.assign(u1) > >>>> > >>>> > >>>> This works fine since u1 is defined by a vector of dofs so the > >>>> assignment is allowed. > >>>> > >>>> Case 2: Time-stepping with user-defined coefficient > >>>> > >>>> # Initializations > >>>> mesh = UnitSquare(32, 32) > >>>> V = FunctionSpace(mesh, "Lagrange", 1) > >>>> w0 = Function(V) > >>>> w1 = Function(V, "sin(t*x[0])") > >>>> > >>>> # Time stepping > >>>> for i in range(10): > >>>> > >>>> print i > >>>> > >>>> # Update w1 > >>>> w1.t = float(i) > >>>> > >>>> # Solve for u > >>>> > >>>> # Assign w0 = w1 (does not work) > >>>> w0.assign(w1) > >>>> #w0 = interpolate(w1, V) > >>>> #w0 = project(w1, V) > >>>> > >>>> This breaks since assignment is not allowed from the user-defined > >>>> Function w1. Interpolation or projection helps, but each of these > >>>> return a new function, which will confuse the JIT compiler (at least > >>>> the current FFC JIT compiler) and lead to excessive generation of code > >>>> (no cache reuse). > >>>> > >>>> The new version of the assignment operator allows this kind of > >>>> assignment and automatically interpolates when necessary. It also > >>>> prints out the following message: > >>>> > >>>> Assignment from user-defined function, interpolating. > >>>> > >>>> So it should be clear what happens. Any objections? > >>> Any more thoughts on this? > >>> > >> > >> Below is how I would prefer to have it: > >> > >> Case 1: Time-stepping with user-defined initial data > >> > >> # Initializations > >> mesh = UnitSquare(32, 32) > >> V = FunctionSpace(mesh, "Lagrange", 1) > >> u0 = Function(V, "sin(x[0])") > >> u1 = Function(V) > >> > >> # Interpolate the initial data > >> u0.interpolate() > >> > >> # Time stepping > >> for i in range(10): > >> > >> print i > >> > >> # Solve for u1 > >> u1.vector() > >> > >> # Assign u0 = u1 > >> u0.assign(u1) > >> > > > > I've changed my mind for the above example. If u0 is defined as > > > > u0 = Function(V, "sin(x[0])") > > > > (it has a FunctionSpace) then interpolation upon assignment is OK. What > > I would object to is a user-defined Function which does not yet have a > > FunctionSpace being interpolated (which might be prevented by the > > current design anyway). > > > > Garth > > > >> We now have interpolation upon assignment, but why not for example > >> projection upon assignment? Either way the Function being assigned to is > >> not the same as the Function on the RHS of the assignment operator. > >> > >> Garth > >> > >>> > >>> ------------------------------------------------------------------------ > >>> > >>> _______________________________________________ > >>> DOLFIN-dev mailing list > >>> [email protected] > >>> http://www.fenics.org/mailman/listinfo/dolfin-dev > >> > >> > >> _______________________________________________ > >> DOLFIN-dev mailing list > >> [email protected] > >> http://www.fenics.org/mailman/listinfo/dolfin-dev > > > > > > _______________________________________________ > > DOLFIN-dev mailing list > > [email protected] > > http://www.fenics.org/mailman/listinfo/dolfin-dev > > > _______________________________________________ > DOLFIN-dev mailing list > [email protected] > http://www.fenics.org/mailman/listinfo/dolfin-dev
signature.asc
Description: Digital signature
_______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
