[petsc-dev] Making FE fields in the beginning

[Geoffrey Irving]

On Wed, Nov 27, 2013 at 11:38 AM, Matthew Knepley <knepley at gmail.com> wrote:
> On Wed, Nov 27, 2013 at 12:43 PM, Geoffrey Irving <irving at naml.us> wrote:
>>
>> Currently, FE fields on a DMPlex are created with the function
>>
>>     DMPlexProjectFunction
>>
>> This function has three issues:
>>
>> 1. It doesn't take a void* argument, and therefore requires global
>> variables if you need to pass in additional data.
>
> I am not opposed to DMPlex functions taking contexts. However, my
> question is: how is intended to receive this context? So far, I am
> unconvinced that PetscFE needs an outside context. It just seems like
> bad design to me.

Currently the PetscFE functions take raw function pointers with no
context.  The new one would take a function pointer and a context.  If
it's bad design for these functions to take a context, why isn't it
bad design for them to take raw function pointers?

> The original intention is to pass evaluation data in as a mesh field. How
> does this break down?

Who makes the first mesh field?  God?

Less sarcastically, I am trying to shoehorn data into petsc from some
other source.  In this case, the other source is a Python script, but
it doesn't matter: the important thing is that the data does not
magically start out as a mesh field.  It also does not start out as a
parameter free C function pointer.

If you think it's best, I can hardwire understanding of the layout of
mesh fields into my code and just build them directly, but this is
rather annoying to have to do for anything higher than linear
elements.

>> 2. It operates one quadrature point and a time, so is inefficient in
>> situations where evaluation overhead can be amortized across multiple
>> points (especially from scripts).  For the workhorse FE routines this
>> is inconvenient but acceptable since the goal is probably OpenCL, but
>> it'd would be nice to have access to python and numpy in the setup
>> phase.
>
> I think this is oversimplified. DMPlex does operate in batches of quadrature
> points, in fact in batches of cells (although the batch size is 1 now). It is
> PetscDualSpace that calls things one quad point at a time. I would think that
> he Python interface would come in at the level of the PetscDualSpace evaluator. How
> does this differ from what you had planned?

I realize that part of the code is batched already, and don't really
understand the design motivation well enough to have an opinion.  Do
you mean that we'd make a new type of PetscDualSpace that would
internally store the context that matches the function pointers passed
in from DMPlexProjectFunctionLocal?  PetscDualSpace doesn't currently
have any control over the shape of the function pointers that it
receives from the DMPlex routines.

>> 3. It doesn't pass a cell argument, so data in an existing mesh fields
>> requires a hierarchy traversal to access.
>
> I think I do not understand something fundamental about how you are using this
> function. This is intended to be an orthogonal projection of the function f onto the
> finite element space V. Why would we need a cell index or additional information?
> Is it because the function f you have is not analytic, but another mesh field? I will
> think about optimizing this case.

Correct, I'm imagining functions that are not analytic, but are not
already petsc mesh fields.  I know what to do once I make them mesh
fields, but they are not there yet.  I don't need this support at the
moment, but I am going to in the not very distant future.  Simple
examples of nontrivial piecewise analytic functions are NURBS in 2D,
splines in 1D, or subdivision surfaces.

Geoffrey