dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Oct 22, 2008 at 11:17:07PM +0100, Garth N. Wells wrote:
> 
> 
> 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.

Yes! It's starting to look good.

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

I'll probably try adding some simple support for this in FFC right
away. It should just be a few lines of code.

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

ok. I see you've started.

Martin: Note that we can make this work with SFC just as well as for
FFC. You just need to add a flag -l dolfin that adds some sugar on top
of the UFC. One may still use "pure UFC" code (and in particular the
-l dolfin flag won't be used in the Python interface) but then one
needs to write a few extra lines of code. Feel free to copy any
relevant code from FFC (dolfinformat.py).

-- 
Anders

Attachment: signature.asc
Description: Digital signature