The thing about the "special element" finite element spaces is that they can always be embedded in a DG space of some kind. So maybe it's always enough to project to the DG space, and then do interpolation from there?
--cjc On 8 September 2015 at 09:39, Martin Sandve Alnæs <marti...@simula.no> wrote: > As Jan says. > > For some "special elements", the dof evaluation is not point evaluation > but integration over a cell entity (e.g. facet), which is done via > quadrature. For point evaluation dofs, the "quadrature rule" is just a > single point (the dof coordinate) and weight 1.0 for a scalar element, or > weights 1.0 for one component and 0.0 for the other components for vector > elements. > > So a more generic model for evaluation of dofs than what we have today > with evaluate_dofs would be something like: > - for each mesh entity there is: > - a set of evaluation points > - a set of dofs > - a small (sparse or dense) matrix of weights such that dofs = matrix * > function components evaluated in points > This way the interpolation can directly tabulate the points (for each mesh > entity, for each point on entity) and do the optimal number of function > evaluations, and there is no ufc::function callback involved which > simplifies the design. > > However this still ignores Piola mappings, so that needs to be handled. > > And I'm not volunteering for the implementation :) > > Martin > > > On 8 September 2015 at 10:23, Jan Blechta <blec...@karlin.mff.cuni.cz> > wrote: > >> On Tue, 8 Sep 2015 08:58:54 +0200 >> Mikael Mortensen <mikael.morten...@gmail.com> wrote: >> >> > >> > > On 07 Sep 2015, at 15:35, Martin Sandve Alnæs <marti...@simula.no> >> > > wrote: >> > > >> > > What if the ufc finite element can return a quadrature rule for >> > > evaluation of dofs instead of evaluate_dofs taking a callback? Then >> > > dolfin can handle evaluation and taking the weighted sum without >> > > involving ufc::function at all. >> > > >> > >> > Not quite sure I follow. evaluate_dofs is performing point >> > evaluations and then some additional work for special elements. I >> > don’t think quadrature rules could help out with this? >> >> I think Martin suggests that instead of execution of DOF code, a >> recipe/formula for evaluating DOF should be passed to DOLFIN. From some >> reason Martin calls the DOF formula by "quadrature rule". >> >> Jan >> >> > >> > M >> > >> > > 7. sep. 2015 09.53 skrev "Mikael Mortensen" >> > > <mikael.morten...@gmail.com <mailto:mikael.morten...@gmail.com>>: Hi >> > > >> > > As you all rememeber the LagrangeInterpolator class can already do >> > > interpolation in parallel on non-matching meshes for Lagrange >> > > elements. I had another look at what it would take to get >> > > interpolation working in parallel for any element and it turns out >> > > I do not need to change all that much. I have implemented a first >> > > solution in the branch >> > > >> https://bitbucket.org/fenics-project/dolfin/branch/mikael/interpolation-parallel >> > > < >> https://bitbucket.org/fenics-project/dolfin/branch/mikael/interpolation-parallel >> > >> > > just to test in a function called >> > > “LagrangeInterpolator::interpolateall". This is what the algorithm >> > > looks like to get interpolation in parallel working for any element >> > > >> > > Parallel interpolation from Function u0 to Function u living on >> > > different meshes with different partitioning: >> > > >> > > 1) Create a set of all different points of u that will require an >> > > eval (like tabulate_all_coordinates but with a set of coordinates) >> > > 2) Create global bounding boxes for the partitioned mesh of u0 and >> > > distribute to all processors. >> > > >> > > 3) Using bounding boxes, compute the processes that *may* own the >> > > points in need of eval. >> > > >> > > 4) Distribute interpolation points to potential owners who >> > > subsequently tries to evaluate u0. If successful, return values >> > > (result of eval) of u0 to owner. 5) Now all the results of u0.eval >> > > will be on the processes that needs it. It remains to run over all >> > > cells locally (for u) and call evaluate_dofs using these values. >> > > This is a bit tricky since evaluate_dofs needs to take a >> > > ufc::function as argument. To this end I have a solution where I >> > > wrap all results of u0.eval in an Expression, calling a map from >> > > point x to result >> > > >> > > >> > > // Create map from coordinate to global result of eval >> > > static std::map<std::vector<double>, std::vector<double>, >> > > lt_coordinate> coords_to_values(lt_coordinate(1.0e-12)); >> > > >> > > …. Fill map coords_to_values using MPI ++ >> > > >> > > // Wrap coords_to_values in an Expression >> > > class WrapperFunction : public Expression >> > > { >> > > public: >> > > >> > > WrapperFunction(int value_shape) : Expression(value_shape) {}; >> > > >> > > mutable std::vector<double> _xx; >> > > >> > > void eval(Array<double>& values, const Array<double>& x) const >> > > { >> > > for (uint j = 0; j < x.size(); j++) >> > > _xx[j] = x[j]; >> > > >> > > const std::vector<double>& v = coords_to_values[_xx]; // <— >> > > Map from x to u0.eval result for (std::size_t j = 0; j < v.size(); >> > > j++) values[j] = v[j]; // Put >> > > u0.eval in values of Expression } >> > > }; >> > > >> > > 6) Run over local mesh calling evaluate_dofs with wrapped function >> > > as the ufc function. >> > > >> > > 7) Update coefficients of local vector with results from >> > > evaluate_dofs >> > > >> > > Finished:-) >> > > >> > > I have tested it for Nedelec elements of order 1 and bubble >> > > elements. Higher order Nedelec elements do not have a >> > > tabulate_coordinates function implemented, which makes it a bit >> > > more difficult to create the list of all interpolation points. I >> > > think this can be solved quite easily, though, with running over >> > > the local mesh once and collecting all x’s used in the evals. >> > > >> > > I do not have any strong opinion on ufc::cell vs dolfin::cell or >> > > collective vs non-collective eval, but I certainly think both these >> > > should be available to the user. The collective eval is not used in >> > > the above algorithm, because it is more efficient to collect all >> > > interpolation points in one large structure and distribute, than it >> > > is to do each point for itself. >> > > >> > > The parallel interpolation routine should probably be put in the >> > > Function class and not the specific LagrangeInterpolator class, but >> > > for now I’ve just been testing. >> > > >> > > M >> > > >> > > >> > >> On 04 Sep 2015, at 11:02, Chris Richardson <ch...@bpi.cam.ac.uk >> > >> <mailto:ch...@bpi.cam.ac.uk>> wrote: >> > >> >> > >> On 03/09/2015 15:54, Martin Sandve Alnæs wrote: >> > >>> On 3 September 2015 at 16:50, Garth N. Wells <gn...@cam.ac.uk >> > >>> <mailto:gn...@cam.ac.uk>> wrote: >> > >>>> On 3 September 2015 at 14:42, Chris Richardson >> > >>>> <ch...@bpi.cam.ac.uk <mailto:ch...@bpi.cam.ac.uk>> wrote: >> > >>>> On 03/09/2015 14:10, Martin Sandve Alnæs wrote: >> > >>>> Part of the complexity in this chain steps from the numerous >> > >>>> dolfin signatures. As discussed before, these also need cleanup >> > >>>> to clarify collective/non-collective operations. While at it, we >> > >>>> could also vectorize on x by allowing a number of coordinates >> > >>>> wherever there's an >> > >>>> x argument in the code below. >> > >>>> For a moment I'd like to ignore the ufc/ffc parts and figure out >> > >>>> how >> > >>>> many eval signatures we really need in dolfin? In particular to >> > >>>> resolve the collective/non-collective part. Then we can design >> > >>>> the callback mechanism to match the ideal dolfin design. >> > >>>> Here's one take (dodging the type of 'cell'): >> > >>>> Is it as simple as just having a "collective" and >> > >>>> "non_collective" eval? >> > >>>> Typically, processes will iterate over some set of "local" >> > >>>> points, asking for evaluations ( - OK, they may be "local" on >> > >>>> one mesh, but not another). >> > >>>> When the point is not "on process", the process needs to get it, >> > >>>> somehow. >> > >>>> But all the other processes are unaware - they are thinking about >> > >>>> different points. >> > >>>> It is quite unusual for all processes to want to evaluate the >> > >>>> same point at the same time, for a simple collective operation. >> > >>>> My thought was that we would have to store up a list of points, >> > >>>> and then do a collective operation to resolve all the >> > >>>> cross-process evaluations. >> > >>>> Or maybe I have missed something... >> > >>> My impression is that Martin is just suggesting the we split evals >> > >>> into non-collective and collective 'groups', and have as few >> > >>> versions of eval in each group. With fewer evals, we have few >> > >>> code paths to think about. >> > >>> The interface for collective evals is orthogonal to this, and >> > >>> something we need to look closely at. >> > >>> Garth >> > >>> Exactly. The collective variants of eval will always end up >> > >>> calling non-collective variants of eval but never the other way >> > >>> around. Also the collective eval variants won't be virtual >> > >>> functions. The evaluation of a GenericFunction on cells that are >> > >>> known to be on the local process differs between Function and >> > >>> Expression, but can be handled with a single signature "virtual >> > >>> void eval(v,x,cell) const = 0;" >> > >>> in GenericFunction, implemented by Function and by the users >> > >>> subclass of Expression. The only thing we need to clean up there >> > >>> in the ufc::function >> > >>> is the ufc::cell vs dolfin::cell. >> > >>> For a collective eval over multiple points, I currently see only >> > >>> two variants: >> > >>> a) each process pass a different set of point(s) and wants the >> > >>> corresponding results >> > >>> b) each process actually pass the same point(s) and wants all >> > >>> results In both cases there is first communication of 'who owns >> > >>> which points', then each process makes calls to the >> > >>> non-collective evals for its own points, >> > >>> then the results are communicated back to the right place. >> > >>> I believe a) is what Chris was talking about as the most common >> > >>> operation. >> > >> >> > >> OK, makes sense. I suppose I am thinking ahead to an >> > >> implementation of 'interpolate', where fitting the collective >> > >> calls in might be tricky. >> > >> >> > >> Chris >> > >> >> > >> -- >> > >> Chris Richardson >> > >> BP Institute >> > >> Madingley Road >> > >> Cambridge CB3 0EZ >> > >> _______________________________________________ >> > >> fenics mailing list >> > >> fenics@fenicsproject.org <mailto:fenics@fenicsproject.org> >> > >> http://fenicsproject.org/mailman/listinfo/fenics >> > >> <http://fenicsproject.org/mailman/listinfo/fenics> >> > >> >> > -- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916
_______________________________________________ fenics mailing list fenics@fenicsproject.org http://fenicsproject.org/mailman/listinfo/fenics
