dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Jul 07, 2005 at 04:30:35PM +0200, Karin Kraft wrote:
> Hello!
> 
> We have made some progress. We have made a simple implementation of the
> general boundary conditions. As a start we have assumed homogeneous
> dirichlet with a constant penalty factor and this works reasonably
> well. 

Excellent!

> Some issues we have discovered:
> 1. We need to compute the local edge numbers. At the moment we do this
> by comparing to all the edges in the cell like so:
> 
>   // Iterate over all edges in the boundary
>   for (EdgeIterator edge(boundary); !edge.end(); ++edge)
>   {
>       ...
>       for(int i = 0; i < cell.noEdges(); i++)
>       {
> 	  if(cell.edgeID(i) == edge->id())
> 	  {
> 	      segment = i;
> 	  }
>       }
> 
> 
> However, this is inefficient. Perhaps the mesh should be extended to
> store this information?

Matt Knepley is working on a new parallel mesh component for FEniCS
and the plan is that we will wrap this in DOLFIN (same as with PETSc)
once it is ready, so maybe we should not invest too much time in
extending the mesh data structures with new data.

I suggest that we use your temporary implementation (looping over
edges locally to find the local edge number), but create it as a
function in Cell (put it close to the functions for computing edge and
face alignments). Something like

    uint Cell::edgeNumber(const Edge& edge) const;

We will also need

    uint Cell::faceNumber(const Face& edge) const;

By putting this as a member function in Cell, you can use the local
data structures directly instead of iterators (the array ce) which
should be a little faster.

> 2. Our eval looks something like this and we want to avoid the
> conditional evaluation. Any ideas?
>  
>  void eval(real block[], const AffineMap& map, uint boundary) const
>   {
>     // Compute geometry tensors
>     real G0_ = map.scale*1e6;
> 
>     // Compute element tensor
> 
>     if( boundary == 0 ) {
> 	block[0] = 0.0;
> 	block[1] = 0.0;
> 	block[2] = 0.0;
> 	
> 	block[3] = 0.0;
> 	block[4] = 3.333333333333318e-01*G0_;
> 	block[5] = 3.333333333333318e-01*G0_;
> 
> 	block[6] = 0.0;
> 	block[7] = 3.333333333333318e-01*G0_;
> 	block[8] = 3.333333333333318e-01*G0_;
>     }
> 
> Regards,
> Karin and Johan

I managed to avoid conditionals when building the dofmap() function by
creating arrays. Maybe there is some solution where the segment
(boundary) argument is used to index an array to get the values, but
it doesn't look like that would be easy to do here.

One solution would be to create three different eval functions
(eval_segment_0, eval_segment_1. eval_segment_2) and then have an
array of function pointers to these functions and use the argument
segment (boundary) to index into this array, get the correct function
and then call it. I'm not sure how this compares to having a switch or
a series of if statements inside eval.

It bugs me that we first have to search for the edge number, and then
when we have the number go through a switch statement. Are we doing
something wrong? Maybe there is a shortcut that avoids both the
searching and the switch?

/Anders