dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Monday August 15 2011 16:18:22 Andre Massing wrote:
> 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.

That was what I also thought. But something funny happens in the 
MeshPrimitive::getEntity. For some reason it gets called using an index 
outside the number of Cells in the domain. So the search tree is not properly 
generated using the SubsetIterator. I suspect it might have something to do 
with the initialization of the search tree is beeing done using size of the 
whole mesh instead of the size of the Subset. 

> >> 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.

Ok.

> > 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)

Ok.

> > 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?

Probably no point of keeping Subset based search trees in the Mesh object. But 
maybe as you point out we can use the canonical search tree to be used for the 
Mesh memberfunctions.

> > 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.

Yes. I just realized that I use the knowledge of all vertices in a Mesh to 
compute the distances between the vertices and the Boundary (the reason I 
added this class in the first place...). This can, however be solved by 
passing a Mesh to the method call.

> >> 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... :(

Sounds complicated...

> >> 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.

Maybe for a reason ;)
 
> >> 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. :/

Sounds complicated...

Johan

> > 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
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp