dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Mon, Apr 24, 2006 at 12:24:39PM -0500, Robert C. Kirby wrote:
> 
> On Apr 24, 2006, at 2:07 AM, Johan Hoffman wrote:
> 
> >It appears that the most straightforward (and general?) way to  
> >implement
> >surface integrals would be to just integrate the cell basis  
> >functions over
> >a boundary/interior facet. Possibly by first throwing away dofs  
> >that are
> >zero on the facet. Using the "trace_tabulate" of FIAT this should  
> >be too
> >much work?
> >
> 
> You have to be a little careful here -- yes you can throw away zero  
> dofs on a facet *if* you are integrating the function.  Sometimes you  
> want to integrate derivatives.  That makes things more complicated.

This is no problem. FFC will know which entries of A^K are zero and
could generate code that only inserts the nonzeros into A but I have a
feeling it's not worth the effort.

> >On the DOLFIN side we would get element matrices corresponding to  
> >all dofs
> >of the cell also when looping over (boundary) edges, typically with  
> >a lot
> >of zero entries. There may be some optimizations here to consider,  
> >but the
> >general algorithm seems rather clear, or?
> >
> >On the FFC side, the reference tensor needs to be modified (using FIAT
> >trace_tabulate?) to reflect integration over an edge/face instead  
> >of over
> >the whole cell.
> >
> 
> Yes, you also need to use the lower-dimensional quadrature rules,  
> etc.  This should just require being careful and not necessarily  
> addressing any major conceptual problems. Except maybe getting the  
> orientation correct -- you might need different routines for doing  
> boundary integrals if for each possibly local boundary number/ 
> orientation.

The code that needs to be generated will look something like

void eval(real block[], AffineMap& map, unsigned int facet)
{
  // Compute geometry tensor
  const real G00 = ...
  const real G01 = ...
  ...

  // Compute reference tensor
  switch ( facet )
  {
  case 0:
    block[0] = ...
    block[1] = ...
    ...
    break;
  case 1:
    ...
  }
}

The geometry tensor looks the same for each facet (orientation) but
the reference tensor differs.

/Anders