dolfin team mailing list archive

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

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

The good thing about it is that it allows building the argument list
dynamically depending on input parameters.

A potential issue with generating FunctionSpace subclasses
is sharing a FunctionSpace between forms compiled separately.
That is easily handled though: FooForm need to take FunctionSpace
as input and check for compatibility with FooFunctionSpace.

--
Martin