On Wed, Nov 27, 2013 at 1:58 PM, Geoffrey Irving <[email protected]> wrote:
> On Wed, Nov 27, 2013 at 11:38 AM, Matthew Knepley <[email protected]> > wrote: > > On Wed, Nov 27, 2013 at 12:43 PM, Geoffrey Irving <[email protected]> > 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
