Subject: [UFC-dev] About changes to UFC suggested in Delft, an alternative approach

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 16 Nov 2006 15:03:04 +0100

The ufc.h version presented in Delft is attached.

It was suggested (outside the sessions) that we add a class "form",
which can have one element_tensor object for each integral of
different dimension. Then we would skip the
"has_boundary_contribution" etc. and just have one tabulate(...) in
class element_tensor. An alternative could be to extend tabulate(...)
with two new arguments to represent several integrals of different
dimension at once:

    /// Tabulate the element tensor contribution for all mesh entities
in this cell
    virtual void tabulate(double* A, const double * const * w, const
cell& c, unsigned int ** entities, unsigned int *num_entities) const =
0;

The new part here is
  unsigned int ** entities, unsigned int *num_entities
which would mean something like this:
      entities[0][i] = which local vertices to integrate over, i<num_entities[0]
      entities[1][i] = which local edges to integrate over, i<num_entities[1]
      entities[2][i] = which local faces to integrate over, i<num_entities[2]
      entities[3][i] = which local volumes to integrate over, i<num_entities[3]

The generated element_tensor can then iterate over entities[2][i] if
it has a boundary contribution in a 3D problem. The contents of this
array is filled by the assembly algorithm, based on which local mesh
entities are on the relevant subset of the mesh (most commonly, a
particular part of the boundary, but could f.ex. be a line inside a 3D
problem).

The alternative of introducing the form class as we mentioned above,
would be something like this:

  class element_tensor
  {
  public:

    /// Return the topological dimension of the integral
    virtual unsigned int integral_dimension() const = 0;

    /// Tabulate the element tensor contribution for this mesh entity
    virtual void tabulate(double* A, const double * const * w, const
cell& c, unsigned int entity) const = 0;

  };

  class form
  {
  public:
    virtual unsigned int rank() const = 0;
    virtual unsigned int num_coefficients() const = 0;

    virtual unsigned int num_tensors() const = 0;

    virtual element_tensor*  create_element_tensor(unsigned int i) const = 0;

    virtual node_map*        create_node_map(unsigned int i) const = 0;
    virtual finite_element*  create_finite_element(unsigned int i) const = 0;
  };

Here, a line integral inside a 3D problem could mean
form::num_tensors() would return 3, and create_element_tensor(i) could
then create three different element_tensor objects with
integral_dimension() equal to 1 (for the line), 2 (for the boundary)
and 3 (for the volume). This is probably a clearer approach, but
introduces one more class. The form class separates the "component
factory" role from the computations, since it has no other role.

Another question that arises now is, should we allow a form to create
element_tensor objects of different ranks? That is, a form could
contain both the matrix and the right hand side contributions. So far
we have said that this should be dealt with as separate element_tensor
components, but if we go for using the form class as a "component
factory" I think perhaps this is a good idea. The backside is that the
interface grows and is now not so minimal anymore.

Discuss!

martin

---

From: Anders Logg <[email-masked]>
Date: Thu, 16 Nov 2006 23:26:16 +0100

> with two new arguments to represent several integrals of different

I have been thinking about this and I suggest we extend the tabulate()
function and remove the stuff about "boundary" and "interior". The
assembly is simply one iteration over all cells and then we compute
all integrals on the cell, including any integrals over facets of the
cell.

One may have a couple of objections to this:

1. It is inefficient if you just need to integrate over the interior
of each cell.

This is not true since in that case the generated function would just
compute integrals over the interior. It would simply skip putting in
any code for the facet contributions.

2. It is inefficient if you just need to integrate over the boundary,
which you may want to do repeatedly for time-dependent boundary
conditions.

This is a good point, but if someone wants to make that optimization,
knowing that some terms are time-dependent and others not, then I
suggest compiling that contribution separately and adding to the
assembly library a function that just integrates over cells that are
on the boundary, something like assembleBoundary().

> particular part of the boundary, but could f.ex. be a line inside a 3D

I suggest a slightly different solution:

  void tabulate(double* A, const double * const * w,
                const cell& c, int boundary_facet);

The integer boundary_facet specifies which facet of the cell is on the
boundary if any. If no facet is on the boundary, then boundary_facet
is set to -1.

> interface grows and is now not so minimal anymore.

I suggest the above solution since (i) I think it solves the problem
and (ii) it actually removes 3 functions from the interface at the
cost of only one extra int argument.

/Anders

> // This is UFC (Unified Form-assembly Code) v. 0.0.
> // This code is released into the public domain.
> //
> // The FEniCS Project (http://www.fenics.org/) 2006.
> 
> #ifndef __UFC_H
> #define __UFC_H
> 
> namespace ufc
> {
> 
>   /// This class defines the data structure for a finite element mesh.
> 
>   class mesh
>   {
>   public:
> 
>     /// Constructor
>     mesh(): topological_dimension(0), geometric_dimension(0), num_entities(0) {}
> 
>     /// Destructor
>     virtual ~mesh() {}
> 
>     /// Topological dimension of the mesh
>     unsigned int topological_dimension;
>     
>     /// Geometric dimension of the mesh
>     unsigned int geometric_dimension;
> 
>     /// Array of the global number of entities of each topological dimension
>     unsigned int* num_entities;
> 
>   };
>   
>   /// This class defines the data structure for a cell in a mesh.
>   
>   class cell
>   {
>   public:
> 
>     /// Constructor
>     cell(): entity_indices(0), coordinates(0) {}
> 
>     /// Destructor
>     virtual ~cell() {}
> 
>     /// Array of global indices for the mesh entities of the cell
>     unsigned int** entity_indices;
> 
>     /// Array of coordinates for the vertices of the cell
>     double** coordinates;
>     
>   };
> 
>   /// This class defines the interface for a general tensor-valued function.
> 
>   class function
>   {
>   public:
> 
>     /// Constructor
>     function() {}
> 
>     /// Destructor
>     virtual ~function() {}
> 
>     /// Evaluate the function at the point x = (x[0], x[1], ...) in the cell
>     virtual void evaluate(double* values, const double* x, const cell& c) const = 0;
> 
>   };
>     
>   /// This class defines the interface for a local-to-global mapping of nodes.
> 
>   class node_map
>   {
>   public:
> 
>     /// Constructor
>     node_map() {}
> 
>     /// Destructor
>     virtual ~node_map() {}
>     
>     /// Return a string identifying the node map
>     virtual const char* description() const = 0;
>     
>     /// Return true iff mesh entities of topological dimension d are needed
>     virtual bool needs_mesh_entities(unsigned int d) const = 0;
> 
>     /// Initialize node map for a new mesh and return true iff
>     /// the finite element should be initialized for each cell
>     virtual bool init_mesh(const mesh& mesh) = 0;
> 
>     /// Initialize node map for given cell
>     virtual void init_cell(const mesh& mesh, const cell& cell) = 0;
> 
>     /// Return the dimension of the global finite element function space.
>     /// This is only valid after calling initialization.
>     virtual unsigned int global_dimension() const = 0;
>     
>     /// Return the dimension of the local finite element function space
>     virtual unsigned int local_dimension() const = 0;
>     
>     /// Tabulate the local-to-global mapping of nodes (degrees of freedom)
>     virtual void tabulate_nodes(unsigned int *nodes, const mesh& m, const cell& c) const = 0;
>   };
> 
>   /// This class defines the interface for a finite element.
> 
>   class finite_element
>   {
>   public:
> 
>     /// Constructor
>     finite_element() {}
> 
>     /// Destructor
>     virtual ~finite_element() {}
> 
>     /// Return a string identifying the finite element
>     virtual const char* description() const = 0;
> 
>     /// Return the dimension of the local finite element function space
>     virtual unsigned int local_dimension() const = 0;
> 
>     /// Return the rank of the value space
>     virtual unsigned int value_rank() const = 0;
> 
>     /// Return the dimension of the value space for axis i
>     virtual unsigned int value_dimension(unsigned int i) const = 0;
> 
>     /// Evaluate basis function i at the point x = (x[0], x[1], ...) in the cell
>     virtual void evaluate_basis(double* values, const double* x, unsigned int i, const cell& c) const = 0;
>     
>     /// Evaluate node i on the function f
>     virtual double evaluate_node(unsigned int n, const function& f, const cell& c) const = 0;
>     
>     /// Interpolate vertex values from nodal values
>     virtual void interpolate_vertex_values(double* vertex_values, const double* nodal_values) const = 0;
>   
>     /// Return the number of sub elements (for a mixed finite element)
>     virtual unsigned int num_sub_elements(unsigned int i) const = 0;
> 
>     /// Return sub element i (for a mixed finite element)
>     virtual const finite_element& sub_element(unsigned int i) const = 0;
>     
>   };
> 
>   /// This class defines the interface for the computation of the
>   /// element tensor corresponding to a form with r + n arguments,
>   /// that is, a mapping
>   ///
>   ///     a : V1 x V2 x ... Vr x W1 x W2 x ... x Wn -> R
>   ///
>   /// with arguments v1, v2, ..., vr, w1, w2, ..., wn. The rank r
>   /// element tensor AK is the local contribution from an element K
>   /// to the rank r global tensor A defined by
>   ///
>   ///     Ai = a(V1, V2, ..., Vr, w1, w2, ..., wn),
>   ///
>   /// where each argument Vj represents the application to the sequence
>   /// of nodal basis functions of Vj and w1, w2, ..., wn are given fixed
>   /// functions (coefficients).
>   
>   class element_tensor
>   {
>   public:
> 
>     /// Constructor
>     element_tensor() {}
> 
>     /// Destructor
>     virtual ~element_tensor() {}
> 
>     /// Return a string identifying the form
>     virtual const char* description() const = 0;
> 
>     /// Return the rank of the element tensor (r)
>     /// Return the number of coefficients (n)
>     /// Return true iff the form has a contribution from the interior
>     virtual bool has_interior_contribution() const = 0;
> 
>     /// Return true iff the form has a contribution from the boundary
>     virtual bool has_boundary_contribution() const = 0;
> 
>     /// Tabulate the interior contribution to the element tensor
>     virtual void tabulate_interior(double* A, const double * const * w, const cell& c) const = 0;
>     
>     /// Tabulate the boundary contribution to the element tensor
>     virtual void tabulate_boundary(double* A, const double * const * w, const cell& c, unsigned int facet) const = 0;
>  
>     /// Create a node map for argument function i (caller responsible for deletion)
>     /// Create a finite element for argument function i (caller responsible for deletion)
> }
> #endif

> _______________________________________________
> UFC-dev mailing list
> [email-masked]
> http://fenics.org/mailman/listinfo/ufc-dev

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 17 Nov 2006 10:35:21 +0100

> boundary if any. If no facet is on the boundary, then boundary_facet

I don't get this. If this is the signature of tabulate, how can it do

> all integrals on the cell, including any integrals over facets of the

I agree that my alternative with two new arguments to tabulate is ugly.

But in your suggestion here, do you mean we will still have _one_
element_tensor object like before, in which the tabulate() function
should compute all kinds of contributions? Then this function must be
called once for each boundary facet (there may be more than one!), and
it still doesn't support f.ex. a line integral in 3D, like the problem
Ridge mentioned. Unless the integer is "unpacked" to something like
this:
-1 : cell
[0,num_facet-1]: face number
[num_facet,num_facets+num_edges-1]: edge number
etc., which I think is really ugly. Or am I wrong?

Note that with the introduction of the "form" class, tabulate() will
look just like you said. With this approach, your counterarguments (1)
and (2) above does not apply, and we can still do one iteration over
the mesh. The number of virtual calls is  exactly the same, the
difference being that you would pass a different integer argument to
"switch" between different behaviour (interior/boundary) and with the
form class the calls are made to different functions instead.

martin

---

From: Anders Logg <[email-masked]>
Date: Fri, 17 Nov 2006 11:16:25 +0100

> > all integrals on the cell, including any integrals over facets of the

The generated code for a form with both an interior contribution and a
boundary contribution typically looks something like

    void tabulate_interior(...)
    { 
      ...
      A[0] = ...
      A[1] = ...
      ...
    }

    void tabulate_boundary(..., facet)
    {
      ...
      switch ( facet )
      {
      case 0:
        A[0] = ...
        A[1] = ...
        ...
      case 1:
      ...  
    }

Now, instead it would be

    void tabulate(..., boundary_facet)
    { 
      // Interior contribution
      ...
      A[0] = ...
      A[1] = ...
      ...
      // Boundary contribution
      if ( boundary_facet == -1 )
        return; // Not on the boundary

      switch ( boundary_facet )
      {
      case 0:
        A[0] += ...
        A[1] += ...
        ...
      case 1:
      ...
    }

This also makes it possible to generate code for discontinuous
Galerkin methods where we also have a contribution from each interior
facet.

> should compute all kinds of contributions?

Yes.

> Then this function must be

I don't understand. The tabulate function gets called once on each
cell meaning that the tabulate function can add contributions from all
parts of the cell it wants, including facets, edges etc.

The integer boundary_facet tells the tabulate function if any of the
cell's facets is on the boundary which the tabulate function may or
may not care about. (It may want to integrate over all facets including
internal in a DG method.)

Also note that the assembly function is not forced to iterate over all
cells in the mesh. It can for example iterate over all cells that
touch the boundary, which is essentially what we do now with
tabulate_interior anyway: iterate over boundary facets and for each
facet pick out the interior cell it is part of.

> form class the calls are made to different functions instead.

I don't think this is necessary.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 17 Nov 2006 13:10:04 +0100

>     void tabulate_boundary(..., facet)

Ok.

>       switch ( boundary_facet )

If on one cell, facet 0 and 1 are both on the boundary, the tabulate()
here must be called twice, once for each of the facets:

tabulate(..., 0)
tabulate(..., 1)

Now, with the above code, the interior is assembled twice for the same
cell. This is of course easily fixed by instead writing:

But how does the line integral fit into this? We can't just speak of
cells and facets if we are to support the example Ridge came up with.
This approach is fine for cells and facets, but it doesn't scale well
if we include more than these two different integrals.

> Galerkin methods where we also have a contribution from each interior

I don't get this. A contribution from _each_ interior facet would
never be a problem, since then tabulate() would know to integrate over
all facets. Then it doesn't need the boundary_facet at all. The
problem is when _some_ of the facets should be integrated over, and
which facets that is would be defined by the mesh. Also, speaking of
only cells and facets was not general enough for Ridges problem.

> parts of the cell it wants, including facets, edges etc.

As said above, this doesn't add up with your example code. The single
integer cannot possibly carry enough information about which faces,
edges, etc should be handled.

> facet pick out the interior cell it is part of.

Of course. The assembly function is always free to iterate over any
subset of the mesh.

> I don't think this is necessary.

Maybe not, we'll see. But I still like the concept. It also allows for
extending a form with a right hand side, which would be a (set of)
element_tensor(s) of different rank. That would be very flexible.

martin

---

From: Anders Logg <[email-masked]>
Date: Fri, 17 Nov 2006 13:37:17 +0100

> here must be called twice, once for each of the facets:

Good point. I forgot that more than one facet may be on the boundary...

>        return; // Not on the boundary

How about putting in markers for facets in the cell data structure?
Something that specifies for each facet if it is on the boundary and
perhaps also which boundary it is part of if there are many?

> if we include more than these two different integrals.

I'm not sure I understand where this comes up. If you have a 3D
domain, where and how would you have a line integral in the
formulation? And is this line integral along a part of the boundary or
is it a line integral throughout the domain?

> which facets that is would be defined by the mesh. 

Exactly. We may want to integrate over interior facets and facets on
the boundary, either both or just one of them. From the argument
boundary_facet the tabulate function will know if any of the facets is
on the boundary (modulo the case where two facets are on the boundary
that I had missed). If it needs to do something on all facets, it just
ignores the boundary_facet argument, if it needs to do something only
on boundary facets, it checks for boundary_facet != -1 and if it needs
to do something on interior facets only it checks for boundary_facet
== -1.

> Also, speaking of

> edges, etc should be handled.

I think it does, if we temporarily ignore that more than one facet in
any given cell may be on the boundary. So we just need to add this
information.

> element_tensor(s) of different rank. That would be very flexible.

Perhaps, but it adds more complexity. It is no problem for the code
generator to generate one single file containing the element tensor
for a bilinear form as well as the element tensor for a linear form:

class my_left_hand_side : public ufc::element_tensor
{

};

class my_right_hand_side : public ufc::element_tensor
{

};

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 17 Nov 2006 15:15:37 +0100

Ok.

> Something that specifies for each facet if it is on the boundary and

Could be. Or a simplified version of the extended tabulate I suggested
earlier, only supporting cells and facets this time:

  tabulate(..., int num_boundary_facets, int *boundary_facets)

Putting these arguments in the cell is basically equivalent. This
function should be capable of assembling the contribution from the
cell as well as all boundary contributions in one call.

But I have another objection to this. What if I only want the
contribution from the boundary? This version doesn't support doing
_only_ the interior or _only_ the boundary contribution. This could be
solved by something like this:

  tabulate(..., bool do_cell, int num_boundary_facets, int *boundary_facets)

which lets the assembly function turn off the integration over a cell.
Then we have:

void element_tensor::tabulate(...)
{
  if(do_cell)
  {
    // assemble from interior
  }
  for(int i=0; i<num_boundary_facets; i++)
  {
     switch(i)
     {
     case 0:
         // assemble from boundary 0
      etc
  }
}

This should work. But I'm still not sure about this, we're just
pushing the complexity to another place.

> perhaps also which boundary it is part of if there are many?

I think the ufc shouldn't bother if the user specifies different data
sources on different parts of the boundary. So "which boundary" isn't
relevant for the ufc, only what coefficients it gets and if it should
integrate or not.

> is it a line integral throughout the domain?

He had a surface with a line boundary somewhere in the 3D domain. But
I'm don't know the details. Could you ask Ridge about this?

<about class form>
> class my_right_hand_side : public ufc::element_tensor

It does add a few lines (but only a few) and a new class to the ufc
specs, but it will at the same time simplify the class element_tensor
and allow a collection of f.ex. my_left_hand_side and
my_right_hand_side be constructed and passed around in a unified way.

martin

---

From: Anders Logg <[email-masked]>
Date: Fri, 17 Nov 2006 15:25:12 +0100

> cell as well as all boundary contributions in one call.

Looks ok to me.

Another option is a single argument bool* facet_on_boundary that
specifies for each facet if it is on the boundary.

> pushing the complexity to another place.

I don't think we need to support this. The user that writes the form
decides what should be assembled. If someone writes

  a = dot(grad(v), grad(u))*dx + v*u*ds

then we should integrate over the interior and the boundary.

If someone wants to integrate only over the boundary, then he may just
write

  a = v*u*ds

and put the rest in a separate form.

It is not natural that if someone writes A = B + C we should generate
code that enables someone to compute B and C separately. We should
generate code for computing A.

Then if someone for some reason (optimization) wants to compute B and
C separately, then compile them separately.

> relevant for the ufc, only what coefficients it gets and if it should

ok.

> I'm don't know the details. Could you ask Ridge about this?

ok, I'll ask him.

> my_right_hand_side be constructed and passed around in a unified way.

I think this should be done at the application level (like in DOLFIN).
DOLFIN has a class named LinearPDE that holds a bilinear and a linear
form. This can be passed around, but the UFC does not need to cover
this abstraction.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 17 Nov 2006 15:48:51 +0100

> specifies for each facet if it is on the boundary.

I have yet another objection. If a cell has two (or more) boundary
facets on different boundaries, where the data sources are different,
setting the coefficients (at the assembly side) might be difficult
when all facets are handled in one tabulate() call?

> then we should integrate over the interior and the boundary.

In one way you're right, we should do both. But doing both isn't the
same as doing both simultaneously. I think this is making a big
assumption: it forces the assembly function to do the boundary and
interior contributions in the same loop. What if the assembly function
for some reason really wants to iterate over the interior mesh first
and then the boundary mesh? (Or maybe I'm constructing a nonexisting
problem here?)

> > I'm don't know the details. Could you ask Ridge about this?

It may be that his problem is better solved by having two forms, one
3D and one 2D (embedded in 3D). But if they should share unknowns,
that may become complicated... The class form could make this easier
by having shared nodemaps between element_tensors. But this is also
possible to do at the application level by comparing
node_map::description().

<about class form>

Ok.

martin

---

From: Anders Logg <[email-masked]>
Date: Fri, 17 Nov 2006 17:26:05 +0100

> when all facets are handled in one tabulate() call?

This on the other hand is a very good point. I need to think about
this one.

/Anders