> On 08 Sep 2015, at 10: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 :)
Ok. Not something to be expected soon then. Meanwhile, with current design of evaluate_dofs, it needs to be done more or less in a similar fashion to what I have done. The algorithm works, I've got my implementation working for any element I have tested, even higher order Nedelec that does not implement tabulate_coordinates. The implementation is completely generic and does not assume anything on the mesh nor the distribution. The routine would be necessary for example for parallel simulations where one stores the solution with a certain mesh decomposition and then wants to restart using a different number of processors (often the case for turbulence simulations). M > > Martin > > > On 8 September 2015 at 10:23, Jan Blechta <blec...@karlin.mff.cuni.cz > <mailto:blec...@karlin.mff.cuni.cz>> wrote: > On Tue, 8 Sep 2015 08:58:54 +0200 > Mikael Mortensen <mikael.morten...@gmail.com > <mailto:mikael.morten...@gmail.com>> wrote: > > > > > > On 07 Sep 2015, at 15:35, Martin Sandve Alnæs <marti...@simula.no > > > <mailto: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> > > > <mailto: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> > > > <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> > > >> <mailto: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> > > >>> <mailto: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> > > >>>> <mailto: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> > > >> <mailto:fenics@fenicsproject.org <mailto:fenics@fenicsproject.org>> > > >> http://fenicsproject.org/mailman/listinfo/fenics > > >> <http://fenicsproject.org/mailman/listinfo/fenics> > > >> <http://fenicsproject.org/mailman/listinfo/fenics > > >> <http://fenicsproject.org/mailman/listinfo/fenics>> > > > >
_______________________________________________ fenics mailing list fenics@fenicsproject.org http://fenicsproject.org/mailman/listinfo/fenics
