dolfin team mailing list archive

[Andre Massing <massing@xxxxxxxxx>]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 08/15/2011 08:05 PM, Johan Hake wrote:
> On Monday August 15 2011 06:51:09 Andre Massing wrote:
>> Hi!
>>
>> I am jumping a bit late into discussion, having been lost the Norwegian
>> mountains :)
>>
>> On 08/12/2011 08:52 AM, Johan Hake wrote:
>>> Any objections to me merging this into trunk?
>>
>> There  exists a blueprint
>>
>> https://blueprints.launchpad.net/dolfin/+spec/bbtree-all-codim
> 
> Now I see that this blueprint is relevant for me. I somewhat overlooked it 
> before, as the description is not that extensive ;)

I admit the description is rather "condensed" ;)

>  
>> which implementation would be include probably most of your
>> functionality. So it would be cool if we could discuss a
>> interface/implementation which meets your specific wished and
>> implemented functionality and the goals of the bluprint without
>> reproducing to much code. Very nice that somebody else is using and
>> working on similar functionality!
> 
> I agree. The only thing that is not covered is to look at a BoundaryMesh 
> instead of the whole mesh.
> 
>>> Will add unittest for all methods soonish.
>>>
>>> The code now resides in:
>>>   lp:~dolfin-core/dolfin/hake
>>
>> I looked into the code, looks nice! I guess the wrapper you wrote could
>> be  simplified, generalized and merged into the IntersectionOperator
>> class by slightly generalizing the conctructor of the
>> IntersectionOperater class.
> 
> Yes. I started out doing this. But then I realized that my "mesh" was a 
> BoundaryMesh and that I needed to do some changes to include MeshFunctions. 
> This would have added BoundaryMesh specific constructors _and_ variables. The 
> latter drove me to implement a specified Class for BoundaryDistances.
> 
>> Right now a mesh has a (default) IntersectionOperator object, which has
>> a tree for all cells of a mesh.
>> We now want to provide additional IntersectionOperator objects
>> each of them having a tree for a subset of entities for some codim, as
>> Johan for instance has implemented, a subset of the boundary mesh.
> 
> If a general subset/whole mesh intersection operator is implemented, can my 
> usercase be included in the intersection operator of a BoundaryMesh. Then my 
> use case will not need any special treatment. We could however add a method to 
> BoundaryMesh to turn a FacetFunction from the original mesh into a 
> CellFunction at the BoundaryMesh. This could be added for VertexFunctions too, 
> if needed.

Yes, that sounds like a good idea.

>  
>> All what the IntersectionOperator(Implementation) class constructor
>> actually does in the end is to pass an  iterator range to the bounding
>> box tree to be constructed. Therefor we just have
>> 1) generalized the constructor that the IntersectionOperator is building
>> a tree for a *subset* of entities (*any* codim) of the mesh instead of
>> all cells.
> 
> Sounds fine. I tried using SubsetIterator but failed miserably with template 
> errors over the whole screen. You might succeed here ;). 
> 
> I fallbacked on using SubMesh, which in it self is a Mesh and I could then 
> just use MeshEntityIterator, which worked fine. But another Mesh is generated 
> with SubMesh and it is suboptimal as it involves copying.

Yes, I stumbled about the same issue (being forced to generate a proper
submesh instead of having just a view of the original mesh).
I guess in the general case where you want to have the proper
connectivity information etc. in the submesh is harder to realize via a
view concept. But in our simple use cases a simple subsetiterator should
be okay.

> 
>> 2) to expose all search/distance related functions from the underlying
>> CGAL bbtree (simply by adding member function delegating it to the CGAL
>> tree)
> 
> Sounds good!
> 
>> 3) think about what concrete instances we want to provide, since we are
>> using templates. The  question is also whether these still should be a
>> part of the mesh and finding a nice and not clumsy interface to keep
>> track of the different IntersectionOperator objects.
> 
> What IntersectionOperators are generated in the present implementation? Only 
> intersection of any entity with the cells of a mesh? 

Yes, exactly.

> 
> Because CGAL does not support traversing a tree with only a subset (right?) we 
> need to generate a CGAL tree for each subset we want to use. 

Yes, right. The search structure is a in principle a binary tree based
upon geometric information like bounding boxes etc. (see below), so only
considering a subset makes the tree highly unbalanced and render it
rather inefficient (since  for instance bounding boxes including several
objects are largely oversized for only a subset of them)

> This can 
> potentially generate several trees. I think the optimal thing is to keep these 
> in the Mesh object, but I am not sure how the interface should look like. 

Yes, that could be a possibility. If we for instance generalize the
IntersectionOperator constructor to take a meshfunction and label (as
the SubsetIterator) maybe we can somehow identify the
IntersectionOperator via the latter. On the other side you probably
don't want to have delegate member functions in the mesh class
except from "canonical"/default intersection operations as there are now
in Mesh, as it will clutter the Mesh interface. But what would be then
the point in keeping them inside the Mesh? All calls would be something like
mesh.intersection_operator(SOME_IDENTIFIER).fmember_function(), so there
little benefit from keeping it inside the mesh, except i.e.
Hmmh....? What do you think?

> 
> Also as I wrote above, my BondaryMesh usecase can be handled by a general 
> subset IntersectionOperator by just using the operator on a BoundaryMesh.

Yes, I guess most of it boils down to have an IntersectionOperator
constructor taking meshfunction and a label as you  did in your

DistanceToBoundaryComputation(const MeshFunction<unsigned int>& labels,
			      uint label,
			      const std::string&
			       kernel_type="SimpleCartesian");
I will have a look at that tomorrow.
	

> 
> 
>> A short note about parallelization.  It is true that the common approach
>> for using the bbtrees is rather serially designed. I thought about this
>> a bit and at least for the search algorithms an idea would be to have a
>> bbtree for each mesh partition (as we have right now I guess) and to
>> merge these forrest than to a tree. 
> 
> Does CGAL support merging of trees with only the knowledge of each tree?

No unfortunately not. CGAL AABB tree is extremely restrictive when is
comes to exposing data. Actually you just have member functions
computing stuff but almost nothing gives you access to the data.
And adding functions like traversing of two tree is not possibly either
with their interface right now. You have to copy paste their code,
remove private attributes/or add member member functions... :(


> 
>> The merged upper part of the tree
>> could be on one (or a few) processor(s) then. 
> 
> What do you mean with upper part?

The part in the tree closer to the root. A bounding box tree is usually
build like this:

0) List of bboxes of all entities in question (the leaves of the tree)
1) Find common bounding for that list of bboxes
2) partition the list into to list (binary treee) according to some
criteria.
3) Start from 1) with each of the lists.

The idea would to use the bbtree for each mesh partition,
which gives you a list of bboxes, namely the one root bbox for each mesh
partition and then build a tree from these ones. These bbtree the is
then the "upper" part, as in higher level/closer to the root in the
overall tree. This upper part could be stored at one processor.
In the usual depth first approach you descend recursively into the next
level at each positive intersection test for pairs of boxes. You don't
want to do this in a parallel implementation, since you then will
constantly switching processors while you are traversing up and down
through the levels. So one idea would be first check all nodes at a
current level (breadth first) for the tree having the root bboxes of the
mesh partition as leaves. That will be at one processor only, no
communication needed. Only in case you an intersection of the root boxes
of the mesh partition, you send the data once to the processor and then
traverse locally on the processors the bbtree of the mesh partition
(usual depth-first search)
As I said, that is just a rough idea to how get tree working in parallel
without to much communication. But literature on bbtrees in parallel is
really sparse and rather uncommon.

> 
>> To reduce communication I
>> would suggest a hybrid breadth first / depth first search. That means a
>> breadth search on the one processor to find out which of the merged
>> trees might be involved. And then communicate the needed data to the
>> corresponding processors / mesh and to do the usual depth search first
>> search. That avoid the back and forth communication in the usual
>> recursive depth first search.
>> Of course that is just a rough thought and will need definitely some
>> work and some thoughts about scalability. But I am quite sure that we
>> could it get running in parallel.
> 
> I am not into the terminology here. How would it relate to the functions in 
> the CGAL tree?

As said above, CGAL is very restrictive and I don't see a way without
copy pasting their code and rewriting some essential parts or  starting
from scratch. :/

- --
Andre

> 
> Johan
> 
>>> Also, what would the best way to make it work in parallel. The distance
>>> from all vertices in a mesh to the closest boundary might not be easy to
>>> compute in Parallel as some vertices residing in one mesh might have the
>>> closest distance in the mesh on another processor.
>>>
>>> I am inclined to think that this is a bad side effect a user need to be
>>> aware of when using this function in parallel. But then I know someone
>>> who would disagree ;)
>>>
>>> Another thing is that the present implementation takes a GenericVector
>>>
>>> representing the return values of the distances at each vertex. Somthing 
> like:
>>>   distances = Function(V)
>>>   
>>>   # Compute distance to Boundary for each vertex
>>>   distance_computation = DistanceToBoundaryComputation(mesh)
>>>   distance_computation.vertex_distances(distances.vector())
>>>
>>> In vertex_distances() I check that the local size of the passed vector
>>> has the same size as the mesh.num_vertices() this gives an error when
>>> running in
>>>
>>> parallel:
>>>  Expected a vector with the same local size as the number of vertices
>>>  (1449)
>>>
>>> in the mesh, got 1427
>>>
>>>  Expected a vector with the same local size as the number of vertices
>>>  (1457)
>>>
>>> in the mesh, got 1441
>>>
>>> I suspect that it has something to do with shared vertices. How do I
>>> access the "correct" number of vertices of a mesh partition and how do I
>>> know which one is only owned by local mesh?
>>>
>>> I figure I have to look into ParallelData, which btw is not wrapped to
>>> Python. We need to add it to dolfin_mesh.h. Will do later...
>>>
>>> Johan
>>>
>>> On Wednesday August 10 2011 14:02:34 Johan Hake wrote:
>>>> Hello!
>>>>
>>>> I have created a class, DistanceToBoundaryComputation, similar to
>>>> IntersectionOperator, which takes a Mesh or a FacetFunction and compute
>>>> distances between any point and the closest point to the Boundary or a
>>>> subset of the boundary given by the FacetFunction.
>>>>
>>>> I have published it together with two demos, cpp and python to
>>>> illustrate some of its functions.
>>>>
>>>>    lp:~johan-hake/dolfin/distance-to-boundary
>>>>
>>>> If the distances from each vertex is computed, will the result be
>>>> similar (not always equal) to the signed distance function, or the
>>>> eikonal equation, but it computes faster.
>>>>
>>>> I am not sure how this best can be integrated into the present dolfin.
>>>> It generates some data, like a BoundaryMesh, which it need to store, so
>>>> it might not be something we want to put into the Mesh class? If I use
>>>> the same lazy initialization that Andre used it might be possible.
>>>>
>>>> Please feel free to have alook at one of the demos to see it in action.
>>>>
>>>> Johan
>>>>
>>>> _______________________________________________
>>>> Mailing list: https://launchpad.net/~dolfin
>>>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>>> Unsubscribe : https://launchpad.net/~dolfin
>>>> More help   : https://help.launchpad.net/ListHelp
>>>
>>> _______________________________________________
>>> Mailing list: https://launchpad.net/~dolfin
>>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>>> Unsubscribe : https://launchpad.net/~dolfin
>>> More help   : https://help.launchpad.net/ListHelp
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~dolfin
>> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~dolfin
>> More help   : https://help.launchpad.net/ListHelp

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.17 (GNU/Linux)
Comment: Using GnuPG with SUSE - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJOSak+AAoJEA79ggnbq9dmdtsIALT88hMuz5ke4Yp7fIblbS6k
1GPlCWP8VEbDBucCKF8jSardcjzSBp4/UHKqd2/Bdw3K0XMWhIStpljk4G0I21V4
bwW8q4LkOlwlXWmUGNYJYdopo6DzGacQWnkzLaP1G/V3QprvKSipRqpPaX8IwI5X
SR/0O5i/7QeFlmHTRbT+66m4EDZJsok4H6mxKkoqWnOsWqsnrpOyDdlynmqaTfN0
H7A9VisRXOj2o3+E4fDGN/K1SEvQoYbi7gB8P4laBY4j9YeNWoVMcWW70w/9GRj1
PN+KwQNt+8VGIn46soq4/TMVKK6JI4U/PMiIvrxOgHK4Y/Ye99ddigeqdlWJOzI=
=wYvH
-----END PGP SIGNATURE-----