dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Mar 12, 2009 at 09:51:38AM +0100, Johan Hake wrote:
> On Thursday 12 March 2009 09:37:58 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?
> 
> Fine with me. 
> 
> Could you add some unit tests for this?
> 
> Johan

Yes, they are currently under sandbox/misc/python. I'll move them to
the unit tests if/when we've agreed on this.

--
Anders

Attachment: signature.asc
Description: Digital signature