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
