On Sun, Apr 08, 2007 at 11:27:47AM +0200, Garth N. Wells wrote:
Anders Logg wrote:
On Sat, Apr 07, 2007 at 02:35:07PM +0200, Garth N. Wells wrote:
Any further thoughts on how to apply Dirichlet boundary conditions? I've
been playing around with it, and with the new UFC interface I don't see
how we can retain the simplicity of the old approach.
Garth
How about the following:
- A boundary condition can be applied in one of two ways: either by
overloading eval() in the class BoundaryCondition, or by suppliying a
Function that defines the boundary condition.
- In the first case, the eval() function looks the same as the eval()
method in Function, but instead of real we have BoundaryValue:
void eval(BoundaryValue* values, const real* x);
When using this method, one sets the values of dofs direcly, so for a
BDM element, one would specify the value of the normal component, not
the function value.
Don't we need a function like pointmap to do this? tabulate_coordinates
will give the position, but we are still missing something that
translates the component of a vector function (0, 1 or 2) to the
appropriate dof (like pointmap did).
Ah, yes. I forgot about this already... Since it's not always possible
to create such a list (only for Lagrange), I would like to avoid it
(and preferrably also tabulate_coordinates) and just use evaluate_dof.
If you look in the generated code for PoissonSystem.form (one of the
FFC demos), we have the following code in evaluate_dof:
virtual double evaluate_dof(unsigned int i,
const ufc::function& f,
const ufc::cell& c) const
{
...
// Nodal coordinates on reference cell
static double X[6][2] = {{0, 0}, {1, 0}, {0, 1}, {0, 0}, {1, 0}, {0, 1}};
// Components for each dof
static unsigned int components[6] = {0, 0, 0, 1, 1, 1};
...
// Evaluate function at coordinates
f.evaluate(values, coordinates, c);
// Pick component for evaluation
return values[components[i]];
}
FFC knows that the element is a tensor-product Lagrange element so it
generates the components array, evaluates the function and picks the
dof from the appropriate component. It would be good if we didn't need
to worry about the components in DOLFIN.
So then my suggestion would be to skip the first option and require
that one always sets a boundary condition by a Function. This can be
very flexible. The Function can either be something user-defined (with
an overloaded eval) or a DiscreteFunction obtained as the solution
from a previous problem.
The problem is then that we need to be able to set boundary conditions
only on a certain sub domain of the mesh, but this should be right way
to do it anyway. The current coordinate-based approach (if x == 0
then...) has only been a temporary solution until we figured out how
to label sub domains (by using a MeshFunction<uint>).
Perhaps we could provide some functionality for defining sub domains
based on coordinates?
/Anders
- The class BoundaryValue remains the same as today. So if one does
values[0] = 3.5;
then the value is set. Otherwise it remains Neumann.
- In the second case, one specifies that the boundary condition should
be u = g on the boundary for g a given Function. For a BDM element,
one would then say that the normal component of u should be the same
as the normal component of g on the boundary. Or do we want to say
that g should be the Dirichlet values of u, so that g would be a
scalar that returns the normal components in the case of BDM?
For more complicated elements, we could just impose Dirichlet boundary
conditions weakly.
Garth
Agree.
/Anders
_______________________________________________
DOLFIN-dev mailing list
DOLFIN-dev@xxxxxxxxxx
http://www.fenics.org/mailman/listinfo/dolfin-dev
