← Back to team overview

dolfin team mailing list archive

Re: pointmap

 

On Tue, Apr 10, 2007 at 10:26:38AM +0200, Garth N. Wells wrote:
> 
> 
> Anders Logg wrote:
> > 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]];
> >   }
> > 
> 
> What about making evaluate_dof a templated function?
> 
>    template <class T>
>    virtual T evaluate_dof(unsigned int i,
>                                 const ufc::function& f,
>                                 const ufc::cell& c) const { . . . }
> 
> User-supplied functions could then supply anything (real, integers, 
> boolean). For applying boundary conditions, a user-supplied function 
> could return
> 
>    std::pair<double, bool> bc
> 
> where bc.second == true implies Dirichlet bc, otherwise no Dirichlet bc.
> 
> Garth

I think the evaluate_dof must be a real. What comes out of
evaluate_dof is what goes into the vector x. It is a real-valued
functional.

Here's my "new" suggestion. We have the following boundary conditions:

    u|_{\Gamma_i} = g_i  for i = 1, 2, ...

where each \Gamma_i is a subset of the boundary.

It looks to me like what we need to specify these boundary conditions
are the following items:

    a matrix and vector (the linear system we want to modify)
    
    a set of functions (the functions g_i)
        
    a mesh (defining the domain, including the boundary)

    a set of markers for sub domains on the mesh

So we need something like the following:


    apply_bc(Matrix& A,
             Vector& b,
             Function& g,
             Mesh& mesh,
             MeshFunction<uint> sub_domains,
             uint sub_domain);

One then calls apply_bc() for each sub domain:

    apply_bc(A, b, g0, mesh, sub_domains, 0);
    apply_bc(A, b, g1, mesh, sub_domains, 1);

(The order of arguments and names can be discussed.)

The MeshFunction is defined on the *facets* of the mesh, meaning the
entities of dimension one less than the cells, and it would be ok to
for example set a Dirichlet condition somewhere inside the domain.

If we implement the above, I think we can cover most cases that we
need to handle. Then, we can provide simple wrappers/utilities for
simple cases. For example, not needing to supply a MeshFunction.
We could put the following in the class BoundaryCondition:

class BoundaryCondition
{
    virtual void eval(real* values, const real* x) = 0;

    bool mark(const real* x) = 0;
};

So a BoundaryCondition would essentially be a UserFunction with an
extra function that returns true in mark() for the part of the
boundary it wants to set Dirichlet conditions on.

/Anders


References