Subject: [UFC-dev] create_sub_element (mixed functions)

From: Anders Logg <[email-masked]>
Date: Fri, 2 Feb 2007 21:45:47 +0100

We had a meeting (at Simula) a few days ago about the inclusion of the
function create_sub_element() in the UFC specification.

What create_sub_element(i) does is that it returns a finite element
for sub element i, if the finite element on which the function is
called is a mixed element, that is, generated as a tensor-product of a
set of finite element spaces.

This can be used for example when someone writes down a mixed
formulation of the Stokes equations and solves for w = (u, p).

Then one may (in DOLFIN) do

    w = pde.solve();
    u = w[0];
    p = w[1];

There are two options:

1. We include create_sub_element() in the specification, but then we
also need to have tabulate_sub_dofs(i) and and create_sub_map(). The
function tabulate_sub_dofs() can be used to extract the dofs for u (or
p) from the vector of dofs for w. Then one may use create_sub_map() to
create a dof_map that can be used to address the extracted dofs. One
will also need num_sub_elements() and num_sub_dofs().

We also need to decide how to do with vector-valued functions that can
be written as tensor-products (like vector Lagrange elements). The
natural would be to treat them in the same way so one may do

    u = w[0];
    p = w[1];
    u0 = u[0];
    u1 = u[1];

2. Remove create_sub_element and related functions entirely. Then to
extract a component or sub function, one needs to work through a form
that defines the extraction. For example:

    a = v*u0*dx
    L = v*u[0]*dx

So one computes the projection. The disadvantage is that a linear
solve is needed (which should be fast since it's just a mass matrix).

The advantage would be that this generalizes to all kinds of elements
and all kinds of extractions:

    w = (p, u[1])
    w = (u[0], p, u[1])
    w = (du[0]/dx, p + du[1]/dy)

Note that simple things like extracting components can be
automatically generated by the compiler, so one would still be able to
do

    u = w[0];
    p = w[1];
    u0 = u[0];
    u1 = u[1];

in DOLFIN, but behind the scenes a mass matrix would be created and a
linear system solved.

I prefer option 2 since it minimizes the complexity of the UFC
interface (removing I think 5 functions) and takes advantage of what
is already in the interface to accomplish something new.

Opinions?

/Anders

---

From: [email-masked]
Date: Sun, 4 Feb 2007 11:32:02 +0100 (CET)

I half-way agree here.

I don't think solving a mass matrix is cheap compared to a slice operation.
I guess you would use about 5 iterations, and each iteration
contains the matrix-vector product which is > 4*n. Hence,
at least 20*n. The slice operation is  < n.

One the other hand I do not think the sub_element abstraction is the
right one. The restriction might be from a sub element, but it might
also be from a sub grid such as the boundary, a coarser gridx, or
a part of the grid in a parallel setting. When we in PyCC have done this
earlier, we
have simply created functions which returns index sets which
are sub sets of the global index set. That way we can use
'take' in numpy.

Martin, can you give us some input on this. You have implemented
several such functions already in the PyCC/UFC framework.

Kent

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

---

From: Anders Logg <[email-masked]>
Date: Sun, 4 Feb 2007 12:58:36 +0100

> at least 20*n. The slice operation is  < n.

Yes, not cheap compared to a slice operation, but cheap compared to
computing the solution that you take the slice from and probably also
the operations you then make on the slice.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 5 Feb 2007 11:04:56 +0100

In one way Anders is definitely right. Projections can do so much
more, it's a solution for a whole range of problems using existing
tools. I would still like to do efficient slice operations, without
computing a projection, but then again I'd also like to do this for
the things Kent mentioned, not only for sub elements. Most importantly
in the UFC setting, I'd like to do this for sub components of a vector
element, even if it's not mixed. Vector elements and mixed elements
are really the same, just with equal or different spaces for the
components.

I also think we need to know these dofs for mixed boundary conditions.
Say I want to set the x-component of the displacement to 0 at one
boundary, but keep the other two components tension free. I can't do
that with a projection, can I? I would need to know the dof numbers.

I'll throw in another suggestion here, which I think should work for
both mixed elements and vector/tensor elements in the same fashion:

- Remove the sub_element functions

- Introduce functionality to get the dofs of a component, f.ex.:

class dof_map
{
    /// Tabulate the local-to-global mapping of dofs for one component
in the value space
    virtual void tabulate_component_dofs(
                                     unsigned int* dofs,
                                     const mesh& m,
                                     const cell& c) const = 0;

    /// Return the number of dofs for component i in the value space
    virtual unsigned int num_component_dofs(unsigned int i) const = 0;
}

Maybe we need tabulate_component_facet_dofs(...) as well, but then we
need num_component_facet_dofs(i) too. I'm sure a better choice of
functions is possible.

martin

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Mon, 5 Feb 2007 11:06:41 +0100

Then again, I guess not all elements have separate dofs for each
vector component? Then this would fail.

martin

---

From: Anders Logg <[email-masked]>
Date: Mon, 5 Feb 2007 12:19:23 +0100

> vector component? Then this would fail.

True.

About the Dirichlet boundary conditions for sub components, FFC
currently has a function

    void pointmap(Point points[], unsigned int components[], ...)

that tabulates the points where the dofs are located on an element and
at the same time tabulates the components associated with these
dofs. DOLFIN then checks which of these points are on the boundary and
for a given component, and may thereby set Dirichlet boundary
conditions for vector valued functions.

This also fails when the vector-valued element is not a simple
tensor-product of scalar elements.

Perhaps we need to do the right thing and set the boundary condition
as a constraint that generally involves more than one dof on each
element?

/Anders