Subject: [UFC-dev] New attempt

From: Anders Logg <[email-masked]>
Date: Fri, 8 Dec 2006 00:22:29 +0100

ok, so let me make a new attempt. How about this:

class cell_tensor
{

  virtual bool is_nonzero() const = 0;
  virtual void tabulate(...) const = 0;

};

class exterior_facet_tensor
{

  virtual bool is_nonzero() const = 0;
  virtual void tabulate(...) const = 0;

};

class interior_facet_tensor
{

  virtual bool is_nonzero() const = 0;
  virtual void tabulate(...) const = 0;

};

class form
{

  ...

  virtual cell_tensor* create_cell_tensor() const = 0;

  virtual exterior_facet_tensor* create_exterior_facet_tensor() const = 0;

  virtual interior_facet_tensor* create_interior_facet_tensor() const = 0;

  ...

};

This way, we can call everything tabulate, there are no overloaded
functions and we follow the convention we've used for the other
classes of having factory functions instead of cluttering the
interface with too many functions.

/Anders

---

From: Ola Skavhaug <[email-masked]>
Date: Fri, 8 Dec 2006 00:27:04 +0100

> interface with too many functions.

I like that design. Similar interfaces for similar classes is good
for understanding the code.

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

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 8 Dec 2006 09:02:34 +0100

On 12/8/06, Ola Skavhaug <[email-masked]> wrote:
> Anders Logg skrev den 08/12-2006 f�lgende:

Three unneeded functions!

> >   virtual interior_facet_tensor* create_interior_facet_tensor() const = 0;

In the assembly:

cell_tensor * ct = form->create_cell_tensor();
if(ct)
{
  // iterate over cells
}

etc.

martin

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 8 Dec 2006 09:07:14 +0100

> cell_tensor * ct = form->create_cell_tensor();

What about:

class cell_integral
{

etc.

class form
{
  virtual cell_integral* create_cell_integral() const = 0;
};

so we get in the assembly:

cell_integral * cintegral = form->create_cell_integral();
if(cintegral)
{
  for( all cells c )
  {
     cintegral->tabulate(..., c)
  }
}

The concepts "cell tensor" and "facet tensor" are a bit strange I
think. Like I argued before, these objects are not really tensors.

martin

---

From: Anders Logg <[email-masked]>
Date: Fri, 8 Dec 2006 10:32:44 +0100

On Fri, Dec 08, 2006 at 10:12:32AM +0100, Martin Sandve Alnæs wrote:

> >Maybe it is strange to think about tensors, but this is difficult to
> >avoid since what actually gets computed by the tabulate() function is
> >a tensor.
> >
> >We have a given form a(v0, v1, ...) corresponding to a global tensor
> >A and we generate code that can be used to compute the global tensor
> >using a specific algorithm, namely assembly of local contributions
> >from the entities of the mesh. Each such local contribution is a
> >tensor and this is what the code computes. The person writing the code
> >for assembly and responsible for calling the tabulate() function must
> >know that what the tabulate() function computes is a tensor.
> 
> Exactly: What the tabulate() function ___computes___ is a tensor.
> The tensor is computed by tabulating the values of a particular integral
> with various input, over a given cell.

True.

> >Also, if you do
> >
> >    integral->tabulate()
> >
> >then what gets tabulated? How do you tabulate an integral?
> 
> It's not the integral that _gets_ tabulated, the integral object
> tabulates the tensor.

ok.

> >If you have
> >
> >    cell_tensor->tabulate()
> >
> >then you tabulate a cell tensor. Granted "cell tensor" is not a simple
> >concept but it is what gets tabulated.
> 
> 
> Given a cell integral, you tabulate the integral values over a given
> cell into the tensor. The tensor is "what gets tabulated", yes,
>    integral->tabulate(tensor, cell, ...)
> 
>    cell_tensor->tabulate(tensor, cell, ...)
> 
> 
> What about:
>    cell_integral->tabulate_tensor(tensor, cell, ...)
> 
> Similar to:
>    dof_map->tabulate_nodes(nodes, cell, ...)

Looks ok to me (but "nodes" --> "dofs").

So do we all agree on this new interface?

> One question: is the tensor the same for each integral? It can be for
> cell and exterior facets, but the tensor for interior facets is maybe
> larger?

Yes, it is larger. It is exactly twice the size in each dimension. So
one needs to call tabulate_nodes() on each of the neighboring cells
and create the larger dof map. Here's how we do this currently in
DOLFIN (still with the old FFC interface which is slightly
different from UFC, with eval() instead of tabulate()):

  // Initialize local data structures (FIXME: reuse)
  const uint m = a.test().spacedim();
  const uint n = a.trial().spacedim();
  real* block = new real[4*m*n];
  int* test_nodes = new int[2*m];
  int* trial_nodes = new int[2*n];
  
  // Compute dof map for cell 0 (+ side)
  a.test().nodemap(a.test_nodes, cell0, mesh);
  a.trial().nodemap(a.trial_nodes, cell0, mesh);

  // Copy dofs to common array
  for (uint i = 0; i < m; i++)
    test_nodes[i] = a.test_nodes[i];
  for (uint i = 0; i < n; i++)
    trial_nodes[i] = a.trial_nodes[i];

  // Compute dof map for cell 1 (- side)
  a.test().nodemap(a.test_nodes, cell1, mesh);
  a.trial().nodemap(a.trial_nodes, cell1, mesh);

  // Copy dofs to common array
  for (uint i = 0; i < m; i++)
    test_nodes[m+i] = a.test_nodes[i];
  for (uint i = 0; i < n; i++)
    trial_nodes[n+i] = a.trial_nodes[i];

  // Compute interior facet tensor
  a.eval(block, map0, map1, det, facet0, facet1, alignment);

  // Add exterior facet tensor to global tensor
  A.add(block, test_nodes, 2*m, trial_nodes, 2*n);

  // Delete local data structures
  delete block;
  delete test_nodes;
  delete trial_nodes;

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Fri, 8 Dec 2006 10:49:13 +0100

> So do we all agree on this new interface?

I added these classes to ufc.h in the hg repos.
You should update the tabulate_tensor signatures properly,
and I haven't changed the comments.

martin

---

From: Anders Logg <[email-masked]>
Date: Fri, 8 Dec 2006 10:54:47 +0100

> and I haven't changed the comments.

ok, I'll have a look at this. I'll be offline for a few hours (going
to the dentist) but will be back later.

/Anders

---

From: Anders Logg <[email-masked]>
Date: Fri, 8 Dec 2006 21:24:12 +0100

> and I haven't changed the comments.

Fixed now.

/Anders