dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Oct 23, 2008 at 09:37:29AM +0100, Garth N. Wells wrote:
> 
> 
> 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

This is what I also had in mind.

And then I guess we need a member function

  const FunctionSpace& function_space(uint i) const;

-- 
Anders

Attachment: signature.asc
Description: Digital signature