dolfin team mailing list archive

[Johan Jansson <jjan@xxxxxxxxxx>]

Bartosz Sawicki wrote:

> On 09/01/09 09:38 AM, Johan Jansson wrote:
> > Hi!
> >
> > This looks interesting, could you please post references to the 
> > papers you mentioned about the algorithm?
>
> Basic concept is from Rivara, but I'd like to avoid recursive 
> algorithms, which has usually efficiency problems.


I don't think recursive algorithms have efficiency problems, what do you 
mean by that? Quicksort is a recursive algorithm for example, and it's 
one of the fastest sorting algorithms.


> I've found in some papers [1], where they use similar algorithm but 
> with  iteration procedure. Propably something like this is also 
> implemented in FEtk.org (GAMer package). I write "probably" because 
> the code is still unpublished, and I can't find any documentation 
> about it.
>
> [1] N. A. Baker, D. Sept, M. J. Holst, J. A. McCammon: The adaptive 
> multilevel finite element solution of the Poisson–Boltzmann equation 
> on massively parallel computers,IBM J. RES. & DEV. VOL. 45 NO. 3/4 
> MAY/JULY 2001


Thanks for the reference, I'll look this paper up.


> Algorithm with I've implemented. It wasn't directly taken from any paper:
>
>   while( there are any marked cells ):
>      forech cell which is marked:
>         find the longest edge
>         if the edge hasn't been bisected:
>            bisect longest edge e
>            remember bisected edge and new vertex
>         else :
>            use vertex which was remembered
>         bisect cell
>      foreach bisected edge:
>         mark all cells which incidence the edge
>
> The weak point of the algorithm is the storing of bisected edges, and 
> serching in them. I use double index STL mapping for that purpose - 
> which works pretty well, but that's the slowest part of the algorithm.
>
> I have some ideas, which can make it faster, but it hasn't been tested 
> yet. Stay tuned ...
>
> The mesh quality is fine, because _only the longest edge_ is bisected.


Ok, good! I see that the algorithm you've implemented is indeed an 
iterative variant of the Rivara algorithm. I think this could be very 
useful to know. But as you mention your implementation has some 
performance issues, could it also be because of the mesh 
re-initialization (connectivity) in each iteration?

Here is a bundle of the Rivara refinement:

http://www.csc.kth.se/~jjan/dolfin-transfer/rivara.bundle

Which exists as:

dolfin/mesh/RivaraRefinement.cpp

together with a performance test in:

sandbox/refinement-perf

You can run the test as:

./demo rivara

for the recursive Rivara refinement and:

./demo iterative-Rivara

for the iterative Rivara refinement.

I realize the iterative version is an early implementation, but these 
are the times I get on my computer:

recursive Rivara: 2.3s
iterative Rivara: 416s

So the recursive implementation is over 100 times faster...

Is it really only the map which makes the iterative version so slow, or 
could it be something simple? Perhaps a re-implementation with a dynamic mesh as in my implementation would speed it up?

Best,
 Johan

PS. Re-sending with link instead of attachment since the list didn't post the original message.