dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Sun, Mar 25, 2007 at 04:46:49PM +0200, Garth N. Wells wrote:
> 
> 
> Anders Logg wrote:
> > On Sun, Mar 25, 2007 at 10:44:58AM +0200, Garth N. Wells wrote:
> >>
> >> Anders Logg wrote:
> >>> On Sat, Mar 24, 2007 at 06:52:04PM +0100, Garth N. Wells wrote:
> >>>> In the old FFC output format, the function
> >>>>
> >>>>    void pointmap(Point points[], unsigned int components[],
> >>>>                  const AffineMap& map) const
> >>>>
> >>>> returns in components[] a degree of freedom identifier (e.g. 0 for u, 1 
> >>>> for v, 2 for p, etc) for each entry in the element tensor. How can we 
> >>>> get this information with the new UFC format? (or how can we avoid 
> >>>> requiring it?)
> >>>>
> >>>> Garth
> >>> Do you mean for evaluating dofs on user-defined functions (to get the
> >>> expansion coefficients in the nodal basis to put in the array w) or
> >>> for setting boundary conditions?
> >>>
> >>> I think that in both cases it should be enough to evaluate the dofs on
> >>> the ufc::function, but we might have missed something. The function
> >>> evaluate_dof() takes a function f that may or may not be vector-valued
> >>> and computes the scalar value of dof i. So for a 2D vector-valued
> >>> Lagrange element of degree 1, dof 0 will be f_0(v0), dof 1 will be
> >>> f_0(v1), dof 2 will be f_0(v2), dof 3 will be f_1(v0) etc. The
> >>> function evaluate() in ufc::function needs to compute all values of
> >>> the possibly tensor-valued function.
> >>>
> >> What I don't see is how to make the link between a user-defined function 
> >> in terms of x and j (j being u0, u1 and p for Stokes), and ufc::function 
> >> which is in terms of x and i (i = 0 --> space_dimension).
> >>
> >> Garth
> > 
> > The ufc::function needs to evaluate all values at once, not one at a
> > time. So for 2D Stokes, the array values would be of length 3.
> > 
> > We could keep the current interface eval(p, i) for Function and then
> > call this for each i to fill in the array values, but perhaps we
> > should change the interface of the class Function in DOLFIN so that
> > for a vector-valued function one needs to set all the values at once.
> > An advantage of this would be that one would not need the
> > 
> >    if ( i == 0 )
> >      return 3.0;
> >    else if ( i == 1 )
> >      return 5.0;
> >    else
> >      return 6.0;
> > 
> > Instead, one would write
> > 
> >    values[0] = 3.0;
> >    values[1] = 5.0;
> >    values[2] = 6.0;
> > 
> > The array values is a contiguous array that may represent something
> > tensor-valued as well. If the function takes values in R^{3x3} then
> > the array has length 9.
> >
> 
> For a vector and scalar function this would be be nicer, but we have to 
> make sure that it works nicely with mixed elements. When applying 
> boundary conditions, not all vector components will always be supplied 
> and something like an boolean array will be required to denote Dirichlet 
> bc's.

This is more difficult than I thought... The problem seems to be that
if we use evaluate_dof from the UFC interface to get the value of dof
number i located on the boundary when evaluated on a given function f
(the Dirichlet boundary condition), we want to have the option of not
saying anything about the value for certain *components* of the
function (for which we want Neumann or do nothing boundary
conditions). So we would need a mapping from component to dof in order
to do this, but we don't have this and I don't think we can have it in
general. (If the element is not a simple tensor product of scalar
elements like vector Lagrange.)

Perhaps we could solve the problem by requiring that a boundary
condition must always be given by a Function? And that the boundary
condition applies to all components that the Function represents. (No
option of not setting the boundary condition.) In the case when one
has a mixed system, then one needs to supply one Function for each sub
system that one wants to set Dirichlet conditions for, and don't set
conditions for sub systems that shouldn't have a Dirichlet
condition. So one could do something like

  set_bc(A, x, u); // u is a Function
  set_bc(A, x, p); // p is a Function

This still does not solve the problem when one wants to set a
Dirichlet condition for u in one part of the domain but not in
another, but then we should support setting boundary conditions for
sub domains which are specified by a MeshFunction (over facets).

On a related note, I have mentioned earlier that Kalkulo (kalkulo.no)
is developing a mesh tool for DOLFIN (paid for by Simula). They are
still working out a few bugs, but it should soon be ready (free and
available from fenics.org as part of FEniCS). The mesh tool lets you
mark (graphically) subsets of the boundary and writes the results as a
DOLFIN XML MeshFunction file.

/Anders