Subject: [UFC-dev] Facet tensor

From: Anders Logg <[email-masked]>
Date: Tue, 5 Dec 2006 19:03:01 +0100

How should we incorporate our newly found (and still expanding)
knowledge of assembly for DG methods into the UFC?

Now, we have a class element_tensor that does two things. It
represents contributions from cells ("interior") as well as from
boundary facets ("boundary").

For DG, we need a third type, namely contributions from interior
facets. We could add a third tabulate function (and has_contribution)
to element_tensor, but perhaps instead we should have three different
classes:

    cell_tensor
    exterior_facet_tensor
    interior_facet_tensor

Possibly, these three could inherit from element_tensor.

Opinions?

/Anders

---

From: [email-masked]
Date: Tue, 5 Dec 2006 21:20:10 +0100 (CET)

> Possibly, these three could inherit from element_tensor.

You mean we should have the name cell_tensor instead of
element_tensor?  I kind of like the word element_tensor. But I
do not have a strong opinion.

Instead of *_facet_tensor, we could also call it macro_element_tensor, since
a facet tensor typically consist of two (or several elements?).

Do you have H(curl) like DG methods where the discontinuity is along the
lines of the tetrahedra (instead of faces), leading macro elements
containing all elements with a common line ?

Kent

---

From: Anders Logg <[email-masked]>
Date: Wed, 6 Dec 2006 10:33:35 +0100

> do not have a strong opinion.

I also like element tensor. We could call it

    element_tensor
    exterior_facet_tensor
    interior_facet_tensor

> Instead of *_facet_tensor, we could also call it
> macro_element_tensor, since a facet tensor typically consist of two
> (or several elements?).

Maybe...

> containing all elements with a common line ?

I don't know. I'm just learning DG as I go along.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Wed, 6 Dec 2006 11:29:18 +0100

> > macro_element_tensor, since a facet tensor typically consist of two

Actually, in a certain way I don't like "element_tensor". I was going
to ignore that, but since we're discussing names...

Assume
  ufc::element_tensor et;

A tensor is something with a rank and shape, which "et" has, that can
be indexed, which "et" can not. So technically, "et" is not a tensor.

However, "et" is an object that can _create_ the element tensor, it is
a sort of operator from the cell (or macro cell, and coefficients) to
the element tensor. The element tensor is really the "double *A". But
of coarse, one could say that "A" is an explicit representation of the
element tensor while "et" is another representation of the same thing.

martin

---

From: Anders Logg <[email-masked]>
Date: Wed, 6 Dec 2006 23:11:44 +0100

> > > > knowledge of assembly for DG methods into the UFC?

I've thought about this some more and we should definitely try to
(in some way) collect the three different cases into one interface.

There are (at least) two options:

1. Put three different cases as members in one class

2. Create three different classes which each has a tabulate() function
and then a forth class (ufc::form) that holds each of the others as a
member function.

Since it is somewhat problematic to keep member data in the interfaces
(it makes the interface more complicated and they would have to be
pointers) I suggest option 1. It also minimizes the number of classes
we have to specify.

So I would suggest something like this:

class form
{

  ...

  virtual void tabulate_element_tensor(double* A,
                                       const double * const * w,
                                       const cell& c) const = 0;

  virtual void tabulate_exterior_facet_tensor(double* A,
                                              const double * const * w,
                                              const cell& c,
                                              unsigned int facet) const = 0;

  virtual void tabulate_interior_facet_tensor(double* A, const double * const * w,
                                              const cell& c0,
                                              const cell& c1,
                                              unsigned int facet0,
                                              unsigned int facet1) const = 0;

  ...

};

The obvious problem here is the very long names. We could reduce this
to

    tabulate_element
    tabulate_exterior_facet
    tabulate_interior_facet

or overload tabulate_facet on the arguments:

    tabulate_element
    tabulate_facet      
    tabulate_facet

or be consequent about the naming:

    tabulate_cell
    tabulate_facet
    tabulate_facet

Then we also need three corresponding functions that return a bool
whether or not the tabulation returns nonzero.

/Anders

---

From: Ola Skavhaug <[email-masked]>
Date: Thu, 7 Dec 2006 08:54:44 +0100

> whether or not the tabulation returns nonzero.

Being consequent about the naming, we loose the word element, which I think is
a bit unfortunate. So, I opt for either your last suggestion, or a something
a-la:

     tabulate_element
     tabulate_<something>
     tabulate_<something>

The facet is really an element, right? I have no idea how to name it though.
Maybe Ridg could comment?

Regarding the three corresponding functions returning a bool; would it be
possible to change the signature of the tabulate_* methods to return a bool?
Or should the methods be called prior to calling the tabulate functions? We
could possibly save three methods here...

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

---

From: Kent-Andre Mardal <[email-masked]>
Date: Thu, 07 Dec 2006 09:38:23 +0100

>                                               unsigned int facet1) const = 0;

In the case with more elements meeting should also be included we could
have 
int numCells, cell* c_set, 

Are the facet0 and facet1 needed ? The meeting edge can be computed from
c0 and c1. I guess it is more efficient, but it does not cost many flops
to compute the meeting edge I guess. 

Kent

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Dec 2006 10:52:31 +0100

On Thu, Dec 07, 2006 at 08:54:44AM +0100, Ola Skavhaug wrote:
> Anders Logg skrev den 06/12-2006 f?lgende:

Didn't you want to loose the word "element"?

> Maybe Ridg could comment?

Maybe you could think of it as an element, but I don't think it is
natural to do it. The basis functions are still defined on each of the
neighboring cells. The difference is that we compute the integral on a
subset of the cell. So in DG, we have two elements which each is
defined on a cell, but then we compute the integral on the common facet.

> could possibly save three methods here...

Do you mean something like this:

The case when there is a contribution:

    bool tabulate_element_tensor(double* A,
                                 const double * const * w,
                                 const cell& c) const
    {
      if ( !A )
        return true;

      A[0] = ...

      return true;
    }

The case when there is no contribution:

    bool tabulate_element_tensor(double* A,
                                 const double * const * w,
                                 const cell& c) const
    {
        return false;
    }

The problem with this is that one needs to call this with arguments

    if ( form.tabulate_element_tensor(0, 0, c) )
       ...

So we can pass 0 as argument for the data, but we still need to
provide a cell and if we don't want to assemble over the cells, then
one still needs to create a dummy cell just to pass as argument.

/Anders

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Dec 2006 10:55:58 +0100

> int numCells, cell* c_set, 

Perhaps, but then we need to specify in some way how they meet and
support all different cases. So there would need to be extra arguments
stating how they meet. I suggest limit ourselves to the cases that we
know that we need.

> to compute the meeting edge I guess. 

I agree that they can be computed, but it seems natural to have
arguments (c0, c1, f0, f1) in this case when we have arguments (c, f)
in the case where we integrate over a facet on the boundary.

/Anders

---

From: Kent-Andre Mardal <[email-masked]>
Date: Thu, 07 Dec 2006 10:58:58 +0100

> stating how they meet. I suggest limit ourselves to the cases that we

Agree! 

But then I think the above functions are fine. 

Kent

---

From: Ola Skavhaug <[email-masked]>
Date: Thu, 7 Dec 2006 11:00:23 +0100

> one still needs to create a dummy cell just to pass as argument.

I'm clearly missing a central point here. I thought that the assembler didn't
know anything about what actions should be taken, and stupidly calls all the
tabulate-methods of the ucf through the interface. To clear up my mind, could
anybody scetch a typical loop over a mesh for me?

Ola

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Dec 2006 11:00:33 +0100

> But then I think the above functions are fine. 

But what do we name them? tabulate_interior_facet_tensor is a little
long.

/Anders

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Dec 2006 11:06:29 +0100

> anybody scetch a typical loop over a mesh for me?

Something like this:

if ( a.cell_tensor_nonzero() )
{
   for (CellIterator cell(mesh); !cell.end(); ++cell)
   {
      a.tabulate_cell_tensor(block, ...);
      A.add(block);
   }
}

if ( a.facet_tensor_nonzero() )
{
   BoundaryMesh boundary(mesh);
   for (FacetIterator cell(boundary); !facet.end(); ++facet)
   {
      a.tabulate_facet_tensor(block, ...);
      A.add(block);
   }
}

So if there is no contribution from the boundary, we don't need to
iterate over the boundary, and in particular we don't need to compute
the boundary.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 7 Dec 2006 11:16:30 +0100

I vote for this:

class form
{
  bool has_cell_contribution()
  bool has_interior_facet_contribution()
  bool has_exterior_facet_contribution()

  void tabulate_cell_contribution(...)
  void tabulate_interior_facet_contribution(...)
  void tabulate_exterior_facet_contribution(...)

  ... rank(), create_dof_map(...) etc
}

But if better names for "interior_facet" and "exterior_facet" come up,
that's great.

martin

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Dec 2006 11:20:11 +0100

> But if better names for "interior_facet" and "exterior_facet" come up,

How about dropping "_contribution" from the last three?

class form
{
 
  void tabulate_cell(...)
  void tabulate_interior_facet(...)
  void tabulate_exterior_facet(...)

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 7 Dec 2006 11:20:48 +0100

> But if better names for "interior_facet" and "exterior_facet" come up,

Try reading this as "English":

form.has_cell_contribution()
meaning "a(u,v) has a contribution from a cell integral"

form.tabulate_cell_contribution(Ae, ...)
meaning "tabulate the contribution from the cell integral of a(u,v)
into the element tensor Ae"

What I'm worried about now is that "exterior_facet" and
"interior_facet" is meaningless for everyone with no DG experience. It
would be great if the regular boundary integral didn't have a name
that is "complicated" by the DG issues.

martin

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Thu, 7 Dec 2006 11:23:52 +0100

>    ... rank(), create_dof_map(...) etc

I don't think it reads correctly.

form.tabulate_cell(Ae, ...)

"let a(u,v) tabulate its cell into Ae"

"let a(u,v) tabulate its cell contribution into Ae"

It's not the cell that is tabulated.
Alternative:

form.tabulate_cell_integral(Ae, ...)

"let a(u,v) tabulate its cell integral into Ae"

martin

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Dec 2006 13:03:21 +0100

> that is "complicated" by the DG issues.

What we have is a generalized concept of "element stiffness matrix".

We drop "stiffness" since it is not specialized to the Laplacian.

    --> element matrix

We change "matrix" to "tensor" since it is not specialized to
matrices. (It also includes vectors, scalars and perhaps higher order
tensors.)

    --> element tensor

An element tensor is something that we compute on each element or cell
and normally one may use "cell" or "element" interchangeably since
each element is defined on a cell and we integrate over the cell.

But now we want to cover two more cases:

   2. Integrate over a facet of a given cell/element
   3. Integrate over a common facet for two given cells/elements

So we have three cases:

   1. One cell, one element
   2. One facet, one cell, one element
   3. One facet, two cells, two elements

So we need a name for each of three cases and my suggestion is

   cell tensor
   exterior facet tensor (a facet that has only one cell neighbor)
   interior facet tensor (a facet that has two cell neighbors)

Perhaps there are other options?

/Anders

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Dec 2006 13:20:18 +0100

> Perhaps there are other options?

Other options would be to use the notation we have in FFC for these
tensors:
  
    long name               LaTeX       Code
    ----------------------------------------
    cell tensor             $A^K$       AK
    exterior facet tensor   $A^{S_e}$   ASe
    interior facet tensor   $A^{S_i}$   ASi

We could then have

    tabulate_ak
    tabulate_ase
    tabulate_asi

but this looks a little cryptic.

/Anders