dolfin team mailing list archive

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

2008/10/23 Garth N. Wells <gnw20@xxxxxxxxx>:
>
>
> Anders Logg wrote:
>> On Thu, Oct 23, 2008 at 10:22:54AM +0200, Martin Sandve Alnæs wrote:
>>> 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.
>>
>> Why not? I thought this would be exactly the thing that shared
>> pointers were good at.
>>
>
> I was planning to use to dual approach which we have elsewhere: if a
> reference is passed to the Form constructor, the Form  does not own the
> FunctionSpace; if a shared_ptr is passed to the Form constructor, the
> FunctionSpace may be shared. Simple usage would look like
>
>   FunctionSpace V;
>   Form a(V, V);
>
> and advanced usage
>
>   shared_ptr<FunctionSpace> V(new FunctionSpace);
>   Form a(V, V);
>
> Garth

Looks good.

Anders: In my comment about "Form a(V,V)" I was assuming
that V were FunctionSpace references, not shared_ptrs.
Consider this code:

FunctionSpace * V = new FunctionSpace(...);
Form a(*V, *V); // using reference arguments
delete V; // a doesn't know V is deleted

Using "a" later will cause segmentation faults (if you're lucky!).

--
Martin