dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

2008/10/23 Garth N. Wells <gnw20@xxxxxxxxx>:
>
>
> Martin Sandve Alnæs wrote:
>>
>> 2008/10/22 Anders Logg <logg@xxxxxxxxx>:
>>>
>>> On Wed, Oct 22, 2008 at 10:00:07PM +0200, Martin Sandve Alnæs wrote:
>>>>
>>>> 2008/10/22 Anders Logg <logg@xxxxxxxxx>:
>>>>>
>>>>> On Wed, Oct 22, 2008 at 07:45:03PM +0100, Garth N. Wells wrote:
>>>>>>
>>>>>> Anders Logg wrote:
>>>>>>>
>>>>>>> On Wed, Oct 22, 2008 at 08:24:18PM +0200, Anders Logg wrote:
>>>>>>>>
>>>>>>>> On Wed, Oct 22, 2008 at 07:05:09PM +0100, Garth N. Wells wrote:
>>>>>>>>>
>>>>>>>>> Anders Logg wrote:
>>>>>>>>>>
>>>>>>>>>> On Wed, Oct 22, 2008 at 12:08:43PM +0100, Garth N. Wells wrote:
>>>>>>>>>>>
>>>>>>>>>>> The thread on the new Function design has digressed from the
>>>>>>>>>>> immediate
>>>>>>>>>>> issue, so I'm restarting it.
>>>>>>>>>>>
>>>>>>>>>>> The issue is how to deal with user defined functions in C++. What
>>>>>>>>>>> if we
>>>>>>>>>>> have a design such that:
>>>>>>>>>>>
>>>>>>>>>>> - All Functions must have a FunctionSpace
>>>>>>>>>>
>>>>>>>>>> Yes.
>>>>>>>>>>
>>>>>>>>>>> - A FunctionSpace does not have to be complete (a complete
>>>>>>>>>>> FunctionSpace
>>>>>>>>>>> having a Mesh, a FiniteElement and a DofMap). As a minimum
>>>>>>>>>>> requirement,
>>>>>>>>>>> and FunctionSpace must have a Mesh.
>>>>>>>>>>
>>>>>>>>>> I think it should always need to be complete. We can avoid a lot
>>>>>>>>>> of
>>>>>>>>>> trouble if we make firm requirements: a Function is always
>>>>>>>>>> associated
>>>>>>>>>> with a FunctionSpace and the FunctionSpace is always completely
>>>>>>>>>> defined.
>>>>>>>>>>
>>>>>>>>>>> - A FiniteElement and/or DofMap can be attached to a
>>>>>>>>>>> FunctionSpace after
>>>>>>>>>>> its creation
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> Related to the functions Function::interpolate(double*
>>>>>>>>>>> coefficients, ..)
>>>>>>>>>>> for interpolating Functions on cells
>>>>>>>>>>>
>>>>>>>>>>> - The new Function::interpolate functions do not take a
>>>>>>>>>>> FiniteElement as
>>>>>>>>>>> an argument, so it is not possible to interpolate a function in a
>>>>>>>>>>> different space. Is it desirable to allow Functions from one
>>>>>>>>>>> space to be
>>>>>>>>>>> interpolated in another?
>>>>>>>>>>
>>>>>>>>>> Yes, it's desirable. This can be done globally when needed:
>>>>>>>>>>
>>>>>>>>>>  Function v(V);
>>>>>>>>>>  Vector coefficients;
>>>>>>>>>>  v.interpolate(coefficients, W);
>>>>>>>>>>
>>>>>>>>>> This will compute the global coefficient vector for v on W.
>>>>>>>>>>
>>>>>>>>>> This is currently implemented with the assumption that the meshes
>>>>>>>>>> for
>>>>>>>>>> V and W are the same but could quite easily be extended to
>>>>>>>>>> non-matching meshes now that functions can be evaluated at
>>>>>>>>>> arbitrary
>>>>>>>>>> points (using GTS).
>>>>>>>>>>
>>>>>>>>>>> If we do this, would:
>>>>>>>>>>>
>>>>>>>>>>> - The above allow FiniteElement types to be checked at runtime
>>>>>>>>>>> for
>>>>>>>>>>> consistency (the FiniteElement passed to Function::interpolate
>>>>>>>>>>> should be
>>>>>>>>>>> the same as the Functions own FiniteElement for discrete
>>>>>>>>>>> Functions. This
>>>>>>>>>>> is what we did with the old design.)>
>>>>>>>>>>>
>>>>>>>>>>> - The above deal with the issue of user-defined functions which
>>>>>>>>>>> have a a
>>>>>>>>>>> FunctionSpace but no FiniteElement?
>>>>>>>>>>>
>>>>>>>>>>> - The above deal with special functions, like the mesh size h?
>>>>>>>>>>
>>>>>>>>>> Here's my suggestion for how to handle initialization of Functions
>>>>>>>>>> in C++.
>>>>>>>>>> There is no need for a circular dependency. First some simple
>>>>>>>>>> facts:
>>>>>>>>>>
>>>>>>>>>>  1. The constructor of a Form may require one or more Functions
>>>>>>>>>>  2. The constructor of a Function requires a FunctionSpace
>>>>>>>>>>
>>>>>>>>>> From this it follows that we must do something like
>>>>>>>>>>
>>>>>>>>>>  Function f(V);
>>>>>>>>>>  Poisson::BilinearForm a;
>>>>>>>>>>  Poisson::LinearForm L(f);
>>>>>>>>>>
>>>>>>>>>> The question is now how to initialize the FunctionSpace V. My
>>>>>>>>>> suggestion would be to extend the code generation to generate code
>>>>>>>>>> for
>>>>>>>>>> creating the FunctionSpace(s), just like we do already when we
>>>>>>>>>> generate UFC code + some extra code for defining the Form classes
>>>>>>>>>> (when using -l dolfin in FFC). This does not make FFC
>>>>>>>>>> DOLFIN-specific
>>>>>>>>>> or DOLFIN FFC-specific. One can still use other form compilers,
>>>>>>>>>> but the
>>>>>>>>>> interface will not be as nice.
>>>>>>>>>>
>>>>>>>>>> So, one would do something like this:
>>>>>>>>>>
>>>>>>>>>>  #include "Poisson.h"
>>>>>>>>>>
>>>>>>>>>>  int main()
>>>>>>>>>>  {
>>>>>>>>>>    Mesh mesh("mesh.xml");
>>>>>>>>>>    Poisson::TestSpace V(mesh);
>>>>>>>>>
>>>>>>>>> Why 'TestSpace'?
>>>>>>>>>
>>>>>>>>> Also, when many functions are present, how would we
>>>>>>>>> identify each one? Could FFC create a class
>>>>>>>>> Poisson::FunctionSpace::foo
>>>>>>>>> when 'foo' is the name of the function in the FFC input?
>>>>>>>>
>>>>>>>> We could let FFC create a number of function space classes:
>>>>>>>>
>>>>>>>>  Poisson::TestSpace
>>>>>>>>  Poisson::TrialSpace
>>>>>>>>  Poisson::FunctionSpace_0
>>>>>>>>  Poisson::FunctionSpace_1
>>>>>>>>  ...
>>>>>>>>
>>>>>>>> or even
>>>>>>>>
>>>>>>>>  Poisson::v::FunctionSpace
>>>>>>>>  Poisson::u::FunctionSpace
>>>>>>>>  Poisson::f::FunctionSpace
>>>>>>>>  Poisson::g::FunctionSpace
>>>>>>>>
>>>>>>>> (but I'm not sure it looks very nice).
>>>>>>>>
>>>>>>>> A complication I didn't think of before is how to handle the case
>>>>>>>> when
>>>>>>>> all spaces are the same (to get reuse of dofmaps). One option would
>>>>>>>> be
>>>>>>>> to either generate just one class if all spaces are the same, and
>>>>>>>> otherwise generate separate clases like above. So either one does
>>>>>>>>
>>>>>>>>  Poisson::FunctionSpace V(mesh);
>>>>>>>>  Function f(V);
>>>>>>>>  Function g(V);
>>>>>>>>  Poisson::LinearForm(f, g);
>>>>>>>>
>>>>>>>> or
>>>>>>>>
>>>>>>>>  Poisson::FunctionSpace_0 V(mesh);
>>>>>>>>  Poisson::FunctionSpace_1 W(mesh);
>>>>>>>>  Function f(V);
>>>>>>>>  Function g(W);
>>>>>>>>  Poisson::LinearForm(f, g);
>>>>>>>>
>>>>>>>>>>    Function f(V);
>>>>>>>>>>    Poisson::BilinearForm a;
>>>>>>>>>>    Poisson::LinearForm L(f);
>>>>>>>>>>    ...
>>>>>>>>>>  }
>>>>>>>>>>
>>>>>>>>>> First V, then f, then L.
>>>>>>>>>>
>>>>>>>>> Would this work for all the Functions in SpecialFunction.h?
>>>>>>>>
>>>>>>>> I haven't thought about it, but yes I guess it would.
>>>>>>>>
>>>>>>>> We can provide a set of predefined FunctionSpaces (like DG on
>>>>>>>> triangles and tetrahedra).
>>>>>>>
>>>>>>> I'm moving the discussion of the removal of DofMapSet over here since
>>>>>>> it's related.
>>>>>>>
>>>>>>> If a Form should own a set of FunctionSpaces (replacing the current
>>>>>>> DofMapSet), that would prevent reuse of FunctionSpaces (and DofMaps)
>>>>>>> across multiple Forms.
>>>>>>>
>>>>>>> Another option would be to always require that a Form is initialized
>>>>>>> with one or more FunctionSpaces. This makes sense and is in agreement
>>>>>>> with the requirement of a FunctionSpace when initializing a Function.
>>>>>>>
>>>>>>> It would then be
>>>>>>>
>>>>>>>  PoissonFunctionSpace V(mesh);
>>>>>>>  Function f(V);
>>>>>>>  Function g(V);
>>>>>>>  PoissonBilinearForm a(V, V);
>>>>>>>  PoissonLinearForm L(V, f, g);
>>>>>>>
>>>>>>> Then it's possible to reuse V across multiple forms (even forms not
>>>>>>> included in Poisson.h).
>>>>>>>
>>>>>> Could work. To keep things clear cut, I would use
>>>>>>
>>>>>>    PoissonFunctionSpace V(mesh);
>>>>>>    Array<FunctionSpace*> VxV(V, V);
>>>>>>    PoissonBilinearForm a(VxV);  // Might also nee a shared_ptr version
>>>>>>    PoissonLinearForm L(V, f, g);
>>>>>>
>>>>>> The first argument to a Form defines the test and/or trial spaces, the
>>>>>> other arguments are the Functions in the Form.
>>>>>>
>>>>>> Garth
>>>>>
>>>>> That looks like a complication. The constructor of a Form takes a
>>>>> variable number of arguments anyway, so it's no big deal to require
>>>>> one FunctionSpace for a linear form and two for a bilinear form.
>>>>> It also looks nicer (avoiding Array and pointers).
>>>>>
>>>> The good thing about it is that it allows building the argument list
>>>> dynamically depending on input parameters.
>>>
>>> When is that necessary?
>>>
>>> If it's needed it's simple enough to generate code for both: one
>>> constructor that takes (V, V) and one that takes VxV.
>>
>> An example would be changing material models in an otherwise
>> fixed application. Of course, one could probably handle that with
>> an application-specific "form factory" function or something.
>> But since it's trivial to do, having both would be nice.
>>
>
> It seems that we've converged on some details of the new Function design. I
> don't think we need to address the FFC generation of helper code for
> creating FunctionSpaces first up because we can patch a FunctionSpace
> together with what we have. Same goes for extracting finite elements and dof
> maps using nice names. What would be nice straight away is an extension of
> the generated dolfin::Form code that takes FunctionSpaces as arguments. Just
> need to figure out to work with  vectors or arrays of shared_ptr in Form.h
> for storing the FunctionSpaces.
>
> Garth

Since a FunctionSpace is likely to be shared between forms or other
parts of the application, it seems a good candidate for using shared_ptr.
However, then the syntax "Form a(V, V)" won't work.

--
Martin