Generalized assembly over any entity type

Just sketching an idea here and looking for comments.

Is it feasible to merge the integral types and matching assembly routines into a single do-it-all routine if we generalize the ufc integral interface? Will this result in a major performance hit?

I'm pretty sure this is technically possible and it might even simplify code some places.

Several features people are working on are limited by the UFC abstractions of integral types (cell, exterior facet, interior facet) in various ways. Overlapping meshes, functions living on different meshes, functions living in lower dimensional manifolds, possibly other things. Maintenance of the assemblers in dolfin seems to require quite a bit of work and will only become harder with more features added. Adding integral types to UFC will result in adding new assembly routines to the assemblers in dolfin. The combination of all this seems to be a bottleneck in the fenics development.

Here are some UFC ideas for a more generalized interface:

// We may need to reintroduce the ufc cell class:
struct ufc::cell
{
    vector<double> vertices; // flattened vertices array
    ufc::cell_shape shape;
    int entity; // local facet number, local vertex number, cell number in macro cell
    int orientation;
};

struct ufc::generic_integral
{
    // We can easily throw in a bunch of utility functions in the
    // generated code to make the integral object queryable:

    vector<size_t> A_dims(); // all virtual const=0
    vector<size_t> A_strides();
    size_t A_size();

    size_t num_cells();
    size_t num_coefficients();

    size_t domain_gdim();
    size_t domain_tdim();
    size_t domain_id();
    size_t domain_cellshape();

    // And now for the generic tabulate_tensor interface:
    virtual void tabulate_tensor2(
        vector<double> & A,
        const vector<ufc::cell> & cells,
        const vector< vector< vector<double> > > & w, // w[cell][coeff][dof]
        ) = 0;
};

// Adaptor functions for specialized integrals is easy
struct interior_facet_integral: public generic_integral
{
    // Current interface:
    virtual void tabulate_tensor(double* A,
                             const double * const * w, // w[coeff][dof]
                             const double* vertex_coordinates_0,
                             const double* vertex_coordinates_1,
                             std::size_t facet_0,
                             std::size_t facet_1) const = 0;

    // I think this represents the additional runtime overhead we get:
    virtual void tabulate_tensor2(
        vector<double> & A,
        const vector<ufc::cell> & cells,
        const vector< vector< vector<double> > > & w, // w[cell][coeff][dof]
        )
    {
        tabulate_tensor(A,
                    w[0],
                    cells[0].vertices, cells[1].vertices,
                    cells[0].entity, cells[1].entity);
    }
};

struct ufc::form
{
    // We could even generalize create_foo_integral somehow.
    // I don't know what integraldomainrepr is, but something more
    // generic than (cell,ext_facet,int_facet) + domain_id:
    map<integraldomainrepr, generic_integral *> create_integrals();
};

Adapting the code generation to this interface is not the main hurdle. The generalized dolfin assembler will be rather complex. However it might replace multiple routines with one, so the question is if that's desirable.

Some notes:

We'll probably still generate code where all dimensions are fixed
Integration is always performed over cell 0
The integration cell and coefficient cells are not necessarily the same
In particular the generated code already knows how many cells there are
The generated code also knows which coefficient lives on which cell(*)
This interface can handle coefficients on different meshes, or on facets
To implement macro cell integrals, the size of A will vary
The interface can of course be introduced without having to implement all possible features before they're needed, but having the flexibility in ufc will make it easier for e.g. Andre Massing to experiment with cut cell configurations etc. :)

(*) This is maybe the point I'm the most unsure about how hard will be to handle in dolfin.

There is a major assumption here that the input "makes sense". In particular the integration cell must be contained by every coefficient cell for coefficients to be defined everywhere.

Comments (8)