dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Mar 26, 2007 at 09:04:37AM +0200, Garth N. Wells wrote:
> 
> 
> Anders Logg wrote:
> > 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.)
> > 
> 
> The vast majority of analyses use Lagrange basis functions, so I think 
> that we should support it, and for more complex basis functions the most 
> flexible approach would be weak enforcement of Dirichlet boundary 
> conditions.

Yes, definitely. We should support everything we need to do with
Lagrange elements, but I still think we will be able to do everything
in a uniform way for Lagrange elements and (not all but many) other
types of elements.

> > 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.) 
> 
> Sounds like something as simple as a slip boundary condition would be 
> become complex (or impossible) to impose.

I realize I forgot to point out something important. FFC now treats
vector-valued Lagrange elements as mixed elements, so each component
is a sub system. This means we can set Dirichlet conditions for
individual components.

Some vector-valued elements like Lagrange will be tensor-products of
scalar elements (what we call MixedElement in FFC), while others like
BDM elements will not be orthogonal. In the former case, we can
specify things for individual components, while in the latter case we
can't (but we can specify normal component = 0 as a Dirichlet
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

Correction. Should be something like

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

The last argument specifies the sub system we set the Dirichlet
condition for. For Stokes, we have a "nested mixed system" so if we
want to set a Dirichlet condition for u[0], we would need to do

  set_bc(A, x, u0, [0, 0]);

I'm not sure what type of argument the last one should be. We could
use a vararg (...) and do

  set_bc(A, x, u0, 0, 0);

So an argument list i, j, k, l would mean sub system l of sub system k
of sub system j of sub system i.

In Python, it would be simple to just send in a list [i, j, k, l] of
"component indices".

> > 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).
> >
> 
> I would like to see the current functionality as well as the possibility 
> to set boundary conditions on subdomains (very useful when the geometry 
> is not exact).

Yes!

/Anders



> > 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
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/mailman/listinfo/dolfin-dev
> > 
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev