← Back to team overview

dolfin team mailing list archive

Re: [HG DOLFIN] Automatically interpolate user-defined functions on assignment

 

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).

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
>


Follow ups

References