dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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 <gnw20@xxxxxxxxx> 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 <logg@xxxxxxxxx>
> >>>>> 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
> >>> DOLFIN-dev@xxxxxxxxxx
> >>> http://www.fenics.org/mailman/listinfo/dolfin-dev
> >>
> >>
> >> _______________________________________________
> >> DOLFIN-dev mailing list
> >> DOLFIN-dev@xxxxxxxxxx
> >> http://www.fenics.org/mailman/listinfo/dolfin-dev
> >
> >
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/mailman/listinfo/dolfin-dev
> >
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev

Attachment: signature.asc
Description: Digital signature