dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Mar 12, 2009 at 09:00:35AM +0000, Garth N. Wells 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

Yes, the first line of the assignment operator is

  if (!v.has_function_space())
    error("Cannot copy Functions which do not have a FunctionSpace.");

-- 
Anders

Attachment: signature.asc
Description: Digital signature