dolfin team mailing list archive

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

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.

--
Martin