dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Wed, Nov 23, 2005 at 11:55:58AM +0100, hetzel.devel@xxxxxx wrote:

> Hi, I like to set/read the array element of a function by hand. If I
> have quadratic order (or higher) the sequence of the elements is not
> the sequence I defined by my mesh. This is because of adding a point
> in the middle of each line (for a quadratic element I need 3 points)
> 
> This is clear so far to me, but how is this sequence organzied? Is
> there an Iterator or something, that helps me to identify the n-th
> value of the function array with the connected coordinates?
> 
> /Haiko

First, you shouldn't need to get/set the degrees of freedom by hand
but I can understand there are situations when you need to do this.

The right way to do it is to iterate over the cells in a mesh and for
each cell ask the finite element of the cell how the degrees of
freedom are organized on the current element.


// Mapping from local to global degrees of freedom
int dofs[element.spacedim()];

// Iterate over cells in mesh
for (CellIterator cell(mesh); !cell.end(); ++cell)
{
    // Compute maps from local to global degrees of freedom
    element.dofmap(dofs, *cell, mesh);

    // Then use the list dofs to get/set the values
    for (unsigned int i = 0; i < element.spacedim(); i++)
    {
        // dofs[i] gives you the global number of the local number i


    }
}


Also take a look at the other member functions of FiniteElement, such
as pointmap() and vertexeval().

If you really want to do things manually, here's a some quick info on
how the degrees of freedom are stored for a Lagrange element: first
the values at all vertices, then all values associated with edges,
then all interior degrees of freedom. This holds both locally and
globally, so in the Vector of degrees of freedom you will first have a
list of mesh.noNodes() values associated with vertices (we should
change name from Node to Vertex at some point) then for quadratics
there will be a list of mesh.noEdges() values at edges. For cubics
the list will be of length 2*mesh.noEdges() followed by a list of
mesh.noCells() values associated with the interior.

For higher order elements, you need to look at the code generated by
FFC for the function FiniteElement::dofmap(), but hopefully you can
do everything you need by calling this function so you shouldn't need
to know how it's implemented.

/Anders