Subject: [FFC-dev] dof locations

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 04 Dec 2006 14:14:51 +0100

Could we add something to FFC to describe where the various degrees of 
freedom live (on vertices, edges, internal)?

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 4 Dec 2006 14:21:01 +0100

> freedom live (on vertices, edges, internal)?

Yes we could, but I'd rather not. Why do we need it? I'd prefer if
DOLFIN did not know anything about dofs, other than how to reorder
them.

I see two options:

 1. Reorder the dofs computed by FFC, not knowing why and how they
    are numbered.

 2. Compute a different mapping in FFC.

I can be convinced otherwise of course, but I don't understand why
this is needed.

/Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 04 Dec 2006 14:35:21 +0100

> DOLFIN did not know anything about dofs, other than how to reorder

Not sure that it's this simple if you want to assemble some terms 
block-wise. Also, for parallel assembly we might need to know where dofs 
lie. I'm still thinking about this so I can't be concrete in what's 
needed just yet.

>  1. Reorder the dofs computed by FFC, not knowing why and how they

This is one option I'm considering. Just take the FFC numbering, and run 
a renumbering algorithm over it to minimise bandwidth.

>  2. Compute a different mapping in FFC.

This could be difficult as it depends on properties of the mesh, 
partition, etc. I think that this should remain simple.

> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev
>

---

From: "Johan Hoffman" <[email-masked]>
Date: Mon, 4 Dec 2006 15:26:56 +0100 (MET)

In connection with the dof-discussion; could we add geometry info to FFC?
Like diameter, normals, tangents etc. From what I understand this info is
already in FIAT, so it would amount to adding a mapping to FFC?

/Johan

---

From: Anders Logg <[email-masked]>
Date: Mon, 4 Dec 2006 15:30:08 +0100

These are functions that are best treated as Functions (coefficients).
We already have MeshSize as a default function in DOLFIN and we will
soon have (probably this week) a function FacetNormal that we will use
for DG.

The mapping we need to make is from FFC to DOLFIN so that one can
write h = MeshSize() in FFC without having to declare the
corresponding function in DOLFIN.

/Anders

On Mon, Dec 04, 2006 at 03:26:56PM +0100, Johan Hoffman wrote:
> DOLFIN-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/dolfin-dev

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 04 Dec 2006 15:32:15 +0100

> already in FIAT, so it would amount to adding a mapping to FFC?

I think that this belongs more to the mesh. Diameter is already in 
DOLFIN, and normals should be added shortly (needed for DG which is 
currently being implemented).

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 4 Dec 2006 15:33:40 +0100

On Mon, Dec 04, 2006 at 08:00:28AM -0600, Kirby, Robert C wrote:
> Note that all of this information is conveniently stored
> in FIAT in the DualBasis.entity_ids dictionary.

Yes, this is what is used by FFC to generate the code that computes
the mapping. If needed, FFC could generate code that propagates this
info to DOLFIN.

/Anders

> As a reference point on this discussion, only the evaluation
> engine in Sundance seems to need to know about locations
> (to get mappings and to figure out who should be nonzero on
> boundaries).  On the other hand, as a more integrated environment
> than DOLFIN/FFC/FIAT, it's much easier to get the information
> whenever you need it -- same runtime library used in defining
> elements,forms, operations, etc.
> 
> Robert C. Kirby
> Associate Professor
> Department of Mathematics and Statistics
> Texas Tech University
> 
> "Mathematical software should be mathematical"
> 
> 
> 
> -----Original Message-----
> From: [email-masked] on behalf of Anders Logg
> Sent: Mon 12/4/2006 7:21 AM
> To: [email-masked]
> Subject: Re: [FFC-dev] dof locations
>

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 04 Dec 2006 16:19:18 +0100

Martin Sandve Alnæs wrote:
> On 12/4/06, Garth N. Wells <[email-masked]> wrote:
>> Anders Logg wrote:
> I've just done some work on parallell assembly for PyCC, using the UFC
> interface. This is basically how I make the node partition:
> 
>   // fill index_set with all global dofs used by elements in the local
> grid partition
>   std::set<int> index_set;
> 
>   // iterate over all mesh cells in local mesh partition
>   pycc::UfcCellIterator *iter = mesh_iter_fac.create_cell_iterator();
>   for(; !iter->end(); iter->next())
>   {
>     const ufc::cell & ufc_cell = iter->get_cell();
> 
>     // get loc2glob from nm
>     ufc_node_map.tabulate_nodes(rows, ufc_mesh, ufc_cell);
> 
>     // insert loc2glob entries into index_set
>     for(int i=0; i<row_size; i++)
>     {
>       index_set.insert(rows[i]);
>     }
>   }
>   delete iter;
> 
>

I see that this generates a set of degrees of freedom for a given 
process for a given cell partitioning, but where do you renumber? Do you 
renumber somewhere so that nodes 0 to m-1 are on processs 1, m to 2m-1 
on process 2, etc?

Garth

> Here the mesh iterator knows the local partition of mesh cells, and
> the ufc_node_map knows the local to global mapping of nodes for a
> given cell, so I simply insert the nodes for each cell into a sorted
> set. Then I extract the node indices of the local node partition from
> index_set afterwards.
> 
> This set of nodes will be overlapping on the inner boundaries, but if
> that isn't what you want you could iterate over the inner boundary and
> call ufc_node_map.tabulate_nodes(...) on the cell at each side of the
> boundary facet to find common nodes.
> 
> This is completely independent on the type of element (except perhaps
> the rare elements that Rob mentioned in Delft, but I think that's a
> small extension).
> 
> martin

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 04 Dec 2006 16:35:38 +0100

Martin Sandve Alnæs wrote:
> I don't, really. I use Epetra in this case, and Epetra_Map handles
> this renumbering. An Epetra_Vector is constructed with a particular
> Epetra_Map, and the vector class has accessor functions for global
> value i or local value i. Epetra_FEVector and Epetra_FECrsMatrix have
> a function SumIntoGlobalValues which I use for the assembly. So in
> this case the local assembly routine actually uses global indices for
> mesh entities.
>

Is this efficient? On a given process, it's likely that you're 
assembling entries that "belong" to a sub-matrix that is residing on 
another process. Epetra probably takes care of the communication (PETSc 
does), but you'll be communicating a lot of values back and forth which 
will hurt performance severely. With appropriate renumbering, only terms 
on the boundaries will be communicated between processes.

Garth

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 04 Dec 2006 17:19:37 +0100

Martin Sandve Alnæs wrote:
> I don't see where values other than on the boundaries would be
> communicated? Only the nodes that are shared between cells in
> different mesh partitions (== nodes on the boundary) will be in
> different node partitions and thus lead to communication. Or am I
> missing something?
> 
> The performance is currently horrible because of a horrible mesh
> partition, so I can't test this right now. So I don't really know.
>

The partition is related to this. Entries in a matrix associated with a 
particular partition are formed on the process to which the mesh 
partition belongs - no problems there. But you need to make sure 
(through the dof numbering) that (nearly) all of the terms computed by a 
process are also stored by that process.

The FFC mapping for vector-valued equations is particularly unsuited to 
this as two unknowns corresponding to a single node (say x and y 
components) are located far from each other (the distance is 1/2 of the 
matrix size).

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 4 Dec 2006 15:32:10 +0100

> lie. I'm still thinking about this so I can't be concrete in what's 

ok. I'm certainly open to generating everything we need to do
efficient assembly, but my feeling is that it won't be needed.

> a renumbering algorithm over it to minimise bandwidth.

ok.

> partition, etc. I think that this should remain simple.

Agree.

/Anders

---

From: "Martin Sandve Alnæs" <[email-masked]>
Date: Tue, 5 Dec 2006 09:29:28 +0100

> process are also stored by that process.

In the index_set built in the beforementioned algorithm, all nodes
related to a cell partition are included. Nodes shared between two
processes (partitions) are stored on both processes. This is
independent of the global dof numbering.

During the assembly, I'm not sure matrix->SumIntoGlobalValues(...)
does any communication at all. At the end of assembly there is a call
I didn't mention, matrix->GlobalAssemble(), which does the
communication of shared values which are then added together.

> components) are located far from each other (the distance is 1/2 of the

I don't see this as being the same issue. Which nodes reside on a
particular process is defined by the cell partition, independent of
the global dof numbering. This is what the algorithm does; given a
global dof numbering and a cell partition, select all the global dofs
touched by cells in this partition. (In the local renumbered dof
vector, the distance between these nodes will be 1/2 of the local
vector size.)

Note that I haven't done much parallell stuff at all, so my arguments
here are pretty much colored by the design of Epetra as I've
understood it so far. But this seems very flexible.

martin