[petsc-dev] Making FE fields in the beginning

[Matthew Knepley]

On Wed, Nov 27, 2013 at 1:58 PM, Geoffrey Irving <irving at naml.us> wrote:

> 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 bad part here is that C has few ways for me to tell it what kinds of
functions
are admissible. However, the convention is that the function depends on no
other information than the real space coordinate x, the field values u(x)
and
derivatives gradU(x), and auxiliary fields a(x) and gradA(x). With this, I
can
optimize the evaluation. If I allow arbitrary information into the
evaluation, I can't
do any optimization, and I don't think it adds anything.


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

I have missed the point of the whole discussion. This is about
bootstrapping data
into the system, Sorry, sometimes it takes a while for me to see the point.

Okay, forget what I said above. Yes, you really want f to be a closure, so
pulling
along the pointer makes perfect sense. We have to change the
PetscDualSpaceApply()
signature too.

I still don't like passing the cell number in since its so specific.

   Matt


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



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20131127/91317730/attachment.html>