librelist archives

« back to archive

Proper way to determine if a quadrant in contained in a tree.

Proper way to determine if a quadrant in contained in a tree.

From:
Ethan Hereth
Date:
2014-10-21 @ 19:37
Good day everybody,

I have a simple question, I have a situation where I have a p4est and a
mesh derived from the p4est. I need to perform a 'flood fill' type
algorithm on all the quadrants in the p4est. The naive and brute force
implementation that I have so far looks something like this:

 // Loop over local trees...
 for (p4est_topidx_t iTree = p4est->first_local_tree; iTree <=
p4est->last_local_tree; ++iTree) {
          p4est_tree_t * tree = p4est_tree_array_index(p4est->trees, iTree);

          // Loop over all quadrants in tree...
          for (size_t iQuad = 0; iQuad < tree->quadrants.elem_count;
++iQuad) {

            // If any of the tree-local quadrants match my target quadrant,
I've found the tree.
            if (nbrOrig == p4est_quadrant_array_index(&(tree->quadrants),
iQuad)) {
              whichTree = iTree;
            }
            if (whichTree != -1)
              break;
          }
          if (whichTree != -1)
            break;
        }

Now that I've located the correct tree, I loop over the mesh as usual,
using the neighbor connectivity to perform the flood fill...

 p4est_mesh_face_neighbor_init(&mfn, p4est, p4estGhost, p4estMesh,
whichTree, nbrOrig);

 while ((nq = p4est_mesh_face_neighbor_next(&mfn, &whichTree, &whichQuad,
&nface, &nrank)) != NULL) {
...
}

This has always worked in the past until I ran into a really big case where
it simply takes forever. And, understandably so I suppose; this is very
brute force and inefficient. Profiling the code shows that most of the time
seems to be spent in p4est_quadrant_array_index and p4est_tree_array_index;
no surprise. I began to look into better ways to do this. I came across the
p4est_quadrant_is_inside_tree function; it would seem that this is exactly
what I am looking for. However, I cannot make it work; the way I am using
it anyway. If I execute the following loop over local trees it seems that
the target quadrant is in every tree.

// Loop over local trees...
 for (p4est_topidx_t iTree = p4est->first_local_tree; iTree <=
p4est->last_local_tree; ++iTree) {
          p4est_tree_t * tree = p4est_tree_array_index(p4est->trees, iTree);
          if (p4est_quadrant_is_inside_tree(tree, nbrOrig))
            std::cout << nbrOrig << " is in tree " << tree << "[" << iTree
<< "]" << std::endl;
          else
            std::cout << nbrOrig << " is NOT in tree " << tree << "[" <<
iTree << "]" << std::endl;
 }

I.e. the second statement never gets printed. Am I misunderstanding
something here? I checked out the source a bit but I'm not familiar enough
with the inner workings of p4est to know for sure what's going on. There
seems to be something a little strange in there though: in the develop
branch, in p4est_bits.c, in p4est_quadrant_is_inside_tree(...) there is the
following bit of code:

 /* check if the end of q is not after the last tree quadrant */
  p4est_quadrant_last_descendant (q, &desc, P4EST_QMAXLEVEL);
  if (p4est_quadrant_compare (&tree->last_desc, q) < 0)
    return 0;

I wonder if there is a problem as the last descendant is calculated but not
used? Maybe a copy and paste error from p4est_quadrant_overlaps_tree(...)
directly above? I tried changing it and the behavior did not change so I'm
probably way off. Is this a bug or do I misunderstand the functionality of
the function?

Ultimately, my question is: how should I be doing this? what's the optimal
way? I need to remove this bottleneck from my code.

Thank you all for your patience,

Ethan Alan

Re: [p4est] Proper way to determine if a quadrant in contained in a tree.

From:
Tobin Isaac
Date:
2014-10-21 @ 21:09
Hi Ethan,

On Tue, Oct 21, 2014 at 03:37:17PM -0400, Ethan Hereth wrote:
> Good day everybody,
> 
> I have a simple question, I have a situation where I have a p4est and a
> mesh derived from the p4est. I need to perform a 'flood fill' type
> algorithm on all the quadrants in the p4est. The naive and brute force
> implementation that I have so far looks something like this:
> 
>  // Loop over local trees...
>  for (p4est_topidx_t iTree = p4est->first_local_tree; iTree <=
> p4est->last_local_tree; ++iTree) {
>           p4est_tree_t * tree = p4est_tree_array_index(p4est->trees, iTree);
> 
>           // Loop over all quadrants in tree...
>           for (size_t iQuad = 0; iQuad < tree->quadrants.elem_count;
> ++iQuad) {
> 
>             // If any of the tree-local quadrants match my target quadrant,
> I've found the tree.
>             if (nbrOrig == p4est_quadrant_array_index(&(tree->quadrants),
> iQuad)) {
>               whichTree = iTree;
>             }
>             if (whichTree != -1)
>               break;
>           }
>           if (whichTree != -1)
>             break;
>         }


Could you tell us a little more about where the nbrOrig pointer comes
from?

> 
> Now that I've located the correct tree, I loop over the mesh as usual,
> using the neighbor connectivity to perform the flood fill...
> 
>  p4est_mesh_face_neighbor_init(&mfn, p4est, p4estGhost, p4estMesh,
> whichTree, nbrOrig);
> 
>  while ((nq = p4est_mesh_face_neighbor_next(&mfn, &whichTree, &whichQuad,
> &nface, &nrank)) != NULL) {
> ...
> }
> 
> This has always worked in the past until I ran into a really big case where
> it simply takes forever. And, understandably so I suppose; this is very
> brute force and inefficient. Profiling the code shows that most of the time
> seems to be spent in p4est_quadrant_array_index and p4est_tree_array_index;
> no surprise. I began to look into better ways to do this. I came across the
> p4est_quadrant_is_inside_tree function; it would seem that this is exactly
> what I am looking for. However, I cannot make it work; the way I am using
> it anyway. If I execute the following loop over local trees it seems that
> the target quadrant is in every tree.
> 
> // Loop over local trees...
>  for (p4est_topidx_t iTree = p4est->first_local_tree; iTree <=
> p4est->last_local_tree; ++iTree) {
>           p4est_tree_t * tree = p4est_tree_array_index(p4est->trees, iTree);
>           if (p4est_quadrant_is_inside_tree(tree, nbrOrig))
>             std::cout << nbrOrig << " is in tree " << tree << "[" << iTree
> << "]" << std::endl;
>           else
>             std::cout << nbrOrig << " is NOT in tree " << tree << "[" <<
> iTree << "]" << std::endl;
>  }
> 
> I.e. the second statement never gets printed. Am I misunderstanding
> something here? I checked out the source a bit but I'm not familiar enough
> with the inner workings of p4est to know for sure what's going on. There
> seems to be something a little strange in there though: in the develop
> branch, in p4est_bits.c, in p4est_quadrant_is_inside_tree(...) there is the
> following bit of code:
> 
>  /* check if the end of q is not after the last tree quadrant */
>   p4est_quadrant_last_descendant (q, &desc, P4EST_QMAXLEVEL);
>   if (p4est_quadrant_compare (&tree->last_desc, q) < 0)
>     return 0;

This might be an error (thanks again for helping us find these),
but it's unrelated to why this isn't working for you, see below.

> 
> I wonder if there is a problem as the last descendant is calculated but not
> used? Maybe a copy and paste error from p4est_quadrant_overlaps_tree(...)
> directly above? I tried changing it and the behavior did not change so I'm
> probably way off. Is this a bug or do I misunderstand the functionality of
> the function?

p4est_quadrant_is_inside_tree() only looks at the coordinate/level
data of the quadrant: if it comes from a tree with a different index,
then its return value is invalid.

> 
> Ultimately, my question is: how should I be doing this? what's the optimal
> way? I need to remove this bottleneck from my code.

I'm assuming that your flood-fill is based on a stack of pointers to
quadrants that need to be processed, correct? The problem with this
approach is that, because the tree index isn't stored with the
quadrant, the "whichTree" is lost and has to be found: in your loop
above, you try to find by testing pointer equality.  I would change
the type of your stack from (p4est_quadrant_t *) to a (p4est_topidx_t,
p4est_locidx_t) tuple, indicating the tree and the index in
tree->quadrants of the quadrant that needs to be processed.  You would
form one of these tuples from the (whichTree, whichQuad) values
returned by p4est_mesh_face_neighbor_next().

Hope that's helpful,
  Toby

> 
> Thank you all for your patience,
> 
> Ethan Alan

Re: [p4est] Proper way to determine if a quadrant in contained in a tree.

From:
Ethan Hereth
Date:
2014-10-22 @ 17:17
On Tue, Oct 21, 2014 at 5:09 PM, Tobin Isaac <tisaac@ices.utexas.edu> wrote:

>
> Hi Ethan,
>
> On Tue, Oct 21, 2014 at 03:37:17PM -0400, Ethan Hereth wrote:
> > Good day everybody,
> >
> > I have a simple question, I have a situation where I have a p4est and a
> > mesh derived from the p4est. I need to perform a 'flood fill' type
> > algorithm on all the quadrants in the p4est. The naive and brute force
> > implementation that I have so far looks something like this:
> >
> >  // Loop over local trees...
> >  for (p4est_topidx_t iTree = p4est->first_local_tree; iTree <=
> > p4est->last_local_tree; ++iTree) {
> >           p4est_tree_t * tree = p4est_tree_array_index(p4est->trees,
> iTree);
> >
> >           // Loop over all quadrants in tree...
> >           for (size_t iQuad = 0; iQuad < tree->quadrants.elem_count;
> > ++iQuad) {
> >
> >             // If any of the tree-local quadrants match my target
> quadrant,
> > I've found the tree.
> >             if (nbrOrig == p4est_quadrant_array_index(&(tree->quadrants),
> > iQuad)) {
> >               whichTree = iTree;
> >             }
> >             if (whichTree != -1)
> >               break;
> >           }
> >           if (whichTree != -1)
> >             break;
> >         }
>
>
> Could you tell us a little more about where the nbrOrig pointer comes
> from?
>

nbrOrig comes from a std::queue of 'seed' p4est_quadrant_t *s that I select
from the p4est. I am not sure that answers your question or if it's even
relevant.

>
> >
> > Now that I've located the correct tree, I loop over the mesh as usual,
> > using the neighbor connectivity to perform the flood fill...
> >
> >  p4est_mesh_face_neighbor_init(&mfn, p4est, p4estGhost, p4estMesh,
> > whichTree, nbrOrig);
> >
> >  while ((nq = p4est_mesh_face_neighbor_next(&mfn, &whichTree, &whichQuad,
> > &nface, &nrank)) != NULL) {
> > ...
> > }
> >
> > This has always worked in the past until I ran into a really big case
> where
> > it simply takes forever. And, understandably so I suppose; this is very
> > brute force and inefficient. Profiling the code shows that most of the
> time
> > seems to be spent in p4est_quadrant_array_index and
> p4est_tree_array_index;
> > no surprise. I began to look into better ways to do this. I came across
> the
> > p4est_quadrant_is_inside_tree function; it would seem that this is
> exactly
> > what I am looking for. However, I cannot make it work; the way I am using
> > it anyway. If I execute the following loop over local trees it seems that
> > the target quadrant is in every tree.
> >
> > // Loop over local trees...
> >  for (p4est_topidx_t iTree = p4est->first_local_tree; iTree <=
> > p4est->last_local_tree; ++iTree) {
> >           p4est_tree_t * tree = p4est_tree_array_index(p4est->trees,
> iTree);
> >           if (p4est_quadrant_is_inside_tree(tree, nbrOrig))
> >             std::cout << nbrOrig << " is in tree " << tree << "[" <<
> iTree
> > << "]" << std::endl;
> >           else
> >             std::cout << nbrOrig << " is NOT in tree " << tree << "[" <<
> > iTree << "]" << std::endl;
> >  }
> >
> > I.e. the second statement never gets printed. Am I misunderstanding
> > something here? I checked out the source a bit but I'm not familiar
> enough
> > with the inner workings of p4est to know for sure what's going on. There
> > seems to be something a little strange in there though: in the develop
> > branch, in p4est_bits.c, in p4est_quadrant_is_inside_tree(...) there is
> the
> > following bit of code:
> >
> >  /* check if the end of q is not after the last tree quadrant */
> >   p4est_quadrant_last_descendant (q, &desc, P4EST_QMAXLEVEL);
> >   if (p4est_quadrant_compare (&tree->last_desc, q) < 0)
> >     return 0;
>
> This might be an error (thanks again for helping us find these),
> but it's unrelated to why this isn't working for you, see below.
>
> >
> > I wonder if there is a problem as the last descendant is calculated but
> not
> > used? Maybe a copy and paste error from p4est_quadrant_overlaps_tree(...)
> > directly above? I tried changing it and the behavior did not change so
> I'm
> > probably way off. Is this a bug or do I misunderstand the functionality
> of
> > the function?
>
> p4est_quadrant_is_inside_tree() only looks at the coordinate/level
> data of the quadrant: if it comes from a tree with a different index,
> then its return value is invalid.
>

What is the purpose of this function then? You're saying that it is not
correct to use this function with a quadrant that comes from a tree with a
different index than the tree input to the function? I must be missing
something here.

>
> >
> > Ultimately, my question is: how should I be doing this? what's the
> optimal
> > way? I need to remove this bottleneck from my code.
>
> I'm assuming that your flood-fill is based on a stack of pointers to
> quadrants that need to be processed, correct? The problem with this
> approach is that, because the tree index isn't stored with the
> quadrant, the "whichTree" is lost and has to be found: in your loop
> above, you try to find by testing pointer equality.  I would change
> the type of your stack from (p4est_quadrant_t *) to a (p4est_topidx_t,
> p4est_locidx_t) tuple, indicating the tree and the index in
> tree->quadrants of the quadrant that needs to be processed.  You would
> form one of these tuples from the (whichTree, whichQuad) values
> returned by p4est_mesh_face_neighbor_next().
>
> Thanks a bunch Toby, that's what I needed and I just didn't think about it
that way. This solved my problem!

Thanks for the advice,

Ethan Alan



> Hope that's helpful,
>   Toby
>
> >
> > Thank you all for your patience,
> >
> > Ethan Alan
>

Re: [p4est] Proper way to determine if a quadrant in contained in a tree.

From:
Tobin Isaac
Date:
2014-10-22 @ 17:30

On October 22, 2014 12:17:39 PM CDT, Ethan Hereth 
<advocateddrummer@gmail.com> wrote:
>On Tue, Oct 21, 2014 at 5:09 PM, Tobin Isaac <tisaac@ices.utexas.edu>
>wrote:
>
>>
>> Hi Ethan,
>>
>> On Tue, Oct 21, 2014 at 03:37:17PM -0400, Ethan Hereth wrote:
>> > Good day everybody,
>> >
>> > I have a simple question, I have a situation where I have a p4est
>and a
>> > mesh derived from the p4est. I need to perform a 'flood fill' type
>> > algorithm on all the quadrants in the p4est. The naive and brute
>force
>> > implementation that I have so far looks something like this:
>> >
>> >  // Loop over local trees...
>> >  for (p4est_topidx_t iTree = p4est->first_local_tree; iTree <=
>> > p4est->last_local_tree; ++iTree) {
>> >           p4est_tree_t * tree =
>p4est_tree_array_index(p4est->trees,
>> iTree);
>> >
>> >           // Loop over all quadrants in tree...
>> >           for (size_t iQuad = 0; iQuad <
>tree->quadrants.elem_count;
>> > ++iQuad) {
>> >
>> >             // If any of the tree-local quadrants match my target
>> quadrant,
>> > I've found the tree.
>> >             if (nbrOrig ==
>p4est_quadrant_array_index(&(tree->quadrants),
>> > iQuad)) {
>> >               whichTree = iTree;
>> >             }
>> >             if (whichTree != -1)
>> >               break;
>> >           }
>> >           if (whichTree != -1)
>> >             break;
>> >         }
>>
>>
>> Could you tell us a little more about where the nbrOrig pointer comes
>> from?
>>
>
>nbrOrig comes from a std::queue of 'seed' p4est_quadrant_t *s that I
>select
>from the p4est. I am not sure that answers your question or if it's
>even
>relevant.
>
>>
>> >
>> > Now that I've located the correct tree, I loop over the mesh as
>usual,
>> > using the neighbor connectivity to perform the flood fill...
>> >
>> >  p4est_mesh_face_neighbor_init(&mfn, p4est, p4estGhost, p4estMesh,
>> > whichTree, nbrOrig);
>> >
>> >  while ((nq = p4est_mesh_face_neighbor_next(&mfn, &whichTree,
>&whichQuad,
>> > &nface, &nrank)) != NULL) {
>> > ...
>> > }
>> >
>> > This has always worked in the past until I ran into a really big
>case
>> where
>> > it simply takes forever. And, understandably so I suppose; this is
>very
>> > brute force and inefficient. Profiling the code shows that most of
>the
>> time
>> > seems to be spent in p4est_quadrant_array_index and
>> p4est_tree_array_index;
>> > no surprise. I began to look into better ways to do this. I came
>across
>> the
>> > p4est_quadrant_is_inside_tree function; it would seem that this is
>> exactly
>> > what I am looking for. However, I cannot make it work; the way I am
>using
>> > it anyway. If I execute the following loop over local trees it
>seems that
>> > the target quadrant is in every tree.
>> >
>> > // Loop over local trees...
>> >  for (p4est_topidx_t iTree = p4est->first_local_tree; iTree <=
>> > p4est->last_local_tree; ++iTree) {
>> >           p4est_tree_t * tree =
>p4est_tree_array_index(p4est->trees,
>> iTree);
>> >           if (p4est_quadrant_is_inside_tree(tree, nbrOrig))
>> >             std::cout << nbrOrig << " is in tree " << tree << "["
><<
>> iTree
>> > << "]" << std::endl;
>> >           else
>> >             std::cout << nbrOrig << " is NOT in tree " << tree <<
>"[" <<
>> > iTree << "]" << std::endl;
>> >  }
>> >
>> > I.e. the second statement never gets printed. Am I misunderstanding
>> > something here? I checked out the source a bit but I'm not familiar
>> enough
>> > with the inner workings of p4est to know for sure what's going on.
>There
>> > seems to be something a little strange in there though: in the
>develop
>> > branch, in p4est_bits.c, in p4est_quadrant_is_inside_tree(...)
>there is
>> the
>> > following bit of code:
>> >
>> >  /* check if the end of q is not after the last tree quadrant */
>> >   p4est_quadrant_last_descendant (q, &desc, P4EST_QMAXLEVEL);
>> >   if (p4est_quadrant_compare (&tree->last_desc, q) < 0)
>> >     return 0;
>>
>> This might be an error (thanks again for helping us find these),
>> but it's unrelated to why this isn't working for you, see below.
>>
>> >
>> > I wonder if there is a problem as the last descendant is calculated
>but
>> not
>> > used? Maybe a copy and paste error from
>p4est_quadrant_overlaps_tree(...)
>> > directly above? I tried changing it and the behavior did not change
>so
>> I'm
>> > probably way off. Is this a bug or do I misunderstand the
>functionality
>> of
>> > the function?
>>
>> p4est_quadrant_is_inside_tree() only looks at the coordinate/level
>> data of the quadrant: if it comes from a tree with a different index,
>> then its return value is invalid.
>>
>
>What is the purpose of this function then? You're saying that it is not
>correct to use this function with a quadrant that comes from a tree
>with a
>different index than the tree input to the function? I must be missing
>something here.
>

Maybe the name is misleading: the function indicates whether a quadrant is
in the *local* section of a tree. A quadrant may be a descendant of the 
tree root, but it may lie outside the range assigned to the current 
process.

>>
>> >
>> > Ultimately, my question is: how should I be doing this? what's the
>> optimal
>> > way? I need to remove this bottleneck from my code.
>>
>> I'm assuming that your flood-fill is based on a stack of pointers to
>> quadrants that need to be processed, correct? The problem with this
>> approach is that, because the tree index isn't stored with the
>> quadrant, the "whichTree" is lost and has to be found: in your loop
>> above, you try to find by testing pointer equality.  I would change
>> the type of your stack from (p4est_quadrant_t *) to a
>(p4est_topidx_t,
>> p4est_locidx_t) tuple, indicating the tree and the index in
>> tree->quadrants of the quadrant that needs to be processed.  You
>would
>> form one of these tuples from the (whichTree, whichQuad) values
>> returned by p4est_mesh_face_neighbor_next().
>>
>> Thanks a bunch Toby, that's what I needed and I just didn't think
>about it
>that way. This solved my problem!

Glad to hear it!
  Toby

>
>Thanks for the advice,
>
>Ethan Alan
>
>
>
>> Hope that's helpful,
>>   Toby
>>
>> >
>> > Thank you all for your patience,
>> >
>> > Ethan Alan
>>

Re: [p4est] Proper way to determine if a quadrant in contained in a tree.

From:
Carsten Burstedde
Date:
2014-10-24 @ 10:50
> > >  /* check if the end of q is not after the last tree quadrant */
> > >   p4est_quadrant_last_descendant (q, &desc, P4EST_QMAXLEVEL);
> > >   if (p4est_quadrant_compare (&tree->last_desc, q) < 0)
> > >     return 0;
> >
> > This might be an error (thanks again for helping us find these),
> > but it's unrelated to why this isn't working for you, see below.
> >
> > >
> > > I wonder if there is a problem as the last descendant is calculated but
> > not
> > > used? Maybe a copy and paste error from p4est_quadrant_overlaps_tree(...)
> > > directly above? I tried changing it and the behavior did not change so
> > I'm
> > > probably way off. Is this a bug or do I misunderstand the functionality
> > of
> > > the function?

That last_descendent call is indeed unnecessary.  Fixed this.

Thanks a lot,

Carsten