On 23 June 2015 at 15:31, Martin Sandve Alnæs <[email protected]> wrote: > Both assemble and interpolate have loops over all cells in the mesh, > so at least for single-mesh situations the cell will be known. > > Looking at this from another angle: > what exactly is the problem with a Mesh owning its related BBTree? > Is it only a matter of making updates safe or is there something else? > If we get rid of unsafe ways to modify the Mesh (e.g. MeshGeometry::x()), > and add safe ways to modify the Mesh that mark the BBTree as invalidated, > does that make all issues go away? >
Yes, some issues would go away. In any case, we should avoid user-level functions that expose the underlying storage layout. The design question persists as whether or not it's good for objects to transparently cache data structures that might be quite large. Garth > Martin > > > On 23 June 2015 at 16:05, Anders Logg <[email protected]> wrote: >> >> Who owns the bbtree necessary for Function::interpolate, and for assembly? >> >> For interpolate, we could have an optional tree argument (if not supplied >> it will then be built and thrown away). >> >> For assembly, I can't think if a good solution. Sending in a bbtree to >> assemble() seems far-fetched? >> >> -- >> Anders >> >> >> tis 23 juni 2015 kl 15:43 skrev Garth N. Wells <[email protected]>: >>> >>> On 23 June 2015 at 14:02, Martin Sandve Alnæs <[email protected]> wrote: >>> > On 23 June 2015 at 14:13, Jan Blechta <[email protected]> >>> > wrote: >>> >> >>> >> On Tue, 23 Jun 2015 12:18:05 +0100 >>> >> "Garth N. Wells" <[email protected]> wrote: >>> >> >>> >> > On 23 June 2015 at 10:43, Martin Sandve Alnæs <[email protected]> >>> >> > wrote: >>> >> > > On 23 June 2015 at 11:00, Jan Blechta <[email protected]> >>> >> > > wrote: >>> >> > >> >>> >> > >> On Tue, 23 Jun 2015 08:54:03 +0000 >>> >> > >> Anders Logg <[email protected]> wrote: >>> >> > >> >>> >> > >> > Yes exactly, that's the point, but it would get rid of a few >>> >> > >> > lines of code and look simpler. Users who want to write more >>> >> > >> > lines of code are free to do this: >>> >> > >> > >>> >> > >> > BoundingBoxTree tree; >>> >> > >> > tree.build(mesh); >>> >> > >> > const std::size_t cell_index = >>> >> > >> > tree.compute_first_entity_collision(x); // this may fail and >>> >> > >> > give a non-context related error message Cell cell(mesh, >>> >> > >> > cell_index); function.eval(values, x, cell); >>> >> > >> > >>> >> > >> > Compare with: >>> >> > >> > >>> >> > >> > AutoCell cell; // we don't even need to supply the mesh here >>> >> > >> > function.eval(values, x, cell); >>> >> > >> >>> >> > >> Ok, I'm just saying that this is equivalent to >>> >> > >> >>> >> > >> BoundingBoxTree tree(mesh); >>> >> > >> function.eval(values, x, tree); >>> >> > >> >>> >> > >> while not introducing new class. But if it creates some benefit >>> >> > >> that AutoCell inherits from Cell, then why not... >>> >> > >> >>> >> > >> Jan >>> >> > > >>> >> > > >>> >> > > Because it's a strange abstraction and doesn't solve any problems >>> >> > > that needs solving, but introduces new hidden pitfalls: >>> >> > > >>> >> > > AutoCell cell; // we don't even need to supply the mesh here >>> >> > > function.eval(values, x, cell); // but the mesh is secretly >>> >> > > attached here >>> >> > > >>> >> > > So now cell contains a big data structure (the tree) and a >>> >> > > reference to the mesh, >>> >> > > and the user still needs to carry the AutoCell object around to >>> >> > > avoid to gain >>> >> > > any benefits related to reuse of the tree. Why not just carry the >>> >> > > tree around? >>> >> > > >>> >> > > With Jan's example, it's much clearer how the data flows: >>> >> > > BoundingBoxTree tree(mesh); >>> >> > > function.eval(values, x, tree); >>> >> > > and now the user has a reusable tree in the same amount of code. >>> >> > > >>> >> > >>> >> > I like having the user manage the tree. What about >>> >> > >>> >> > tree = BoundingBoxTree(mesh) >>> >> > cell = tree.find_cell(x) >>> >> > function.eval(values, x, cell) >>> >> > >>> >> > This way Function::eval doesn't need to be parallel-aware. If x is >>> >> > not >>> >> > inside cell, Function::eval can throw and error. >>> >> >>> >> Yes, but what would be the interface for other situations which may >>> >> trigger non-matching eval internally, like interpolate, assemble, ...? >>> >> >>> >> Jan >>> > >>> > >>> > >>> > If the inside check can be kept outside Function::eval, I think it's >>> > best if Function::eval can assume that x is contained in the given >>> > cell to keep things efficient: >>> > >>> >>> I think we should have an (optional) check that x is contained in the >>> cell. The check could be in debug mode only, or turned on/off via a >>> function argument. >>> >>> >>> > tree = BoundingBoxTree(mesh) >>> > cell = tree.find_cell(x) >>> > status = not cell.empty() >>> > if status: >>> > function.eval(values, x, cell) >>> > >>> > In addition there can be a wrapper: >>> > status = function.eval(values, x, tree) >>> > >>> > which does: >>> > cell = tree.find_cell(x) >>> > status = not cell.empty() >>> > if status: >>> > function.eval(values, x, cell) >>> > return status >>> > >>> > So if you have a tree, you can avoid recreating it, >>> > and if you have a cell, you can avoid searching for it. >>> > >>> >>> I like the first, and I'm not bothered either way with the second. >>> >>> Garth >>> >>> > Martin >>> > >>> > >>> >> >>> >> > >>> >> > Garth >>> >> > >>> >> > >>> >> > > Martin >>> >> > > >>> >> > > >>> >> > >> > >>> >> > >> > -- >>> >> > >> > Anders >>> >> > >> > >>> >> > >> > >>> >> > >> > tis 23 juni 2015 kl 10:48 skrev Jan Blechta >>> >> > >> > <[email protected]>: >>> >> > >> > >>> >> > >> > > On Tue, 23 Jun 2015 08:42:49 +0000 >>> >> > >> > > Anders Logg <[email protected]> wrote: >>> >> > >> > > >>> >> > >> > > > How about something like this: >>> >> > >> > > > >>> >> > >> > > > 1. Require additional Cell& cell argument to Function::eval >>> >> > >> > > > >>> >> > >> > > > 2. Add new class AutoCell handling this for users who don't >>> >> > >> > > > want to explicitly work with a BoundingBoxTree >>> >> > >> > > > >>> >> > >> > > > class AutoCell : public Cell >>> >> > >> > > > { >>> >> > >> > > > public: >>> >> > >> > > > AutoCell(Mesh &mesh); >>> >> > >> > > > BoundingBoxTree tree; >>> >> > >> > > > } >>> >> > >> > > > >>> >> > >> > > > AutoCell cell(mesh); >>> >> > >> > > > function.eval(values, x, cell); >>> >> > >> > > > >>> >> > >> > > > 3. If the mesh moves, a user can do something like: >>> >> > >> > > > >>> >> > >> > > > cell.invalidate() >>> >> > >> > > >>> >> > >> > > This is semantically equivalent to passing bounding box. >>> >> > >> > > AutoCell would be just a wrapper for BoundingBoxTree, not >>> >> > >> > > doing anything new. >>> >> > >> > > >>> >> > >> > > Jan >>> >> > >> > > >>> >> > >> > > > >>> >> > >> > > > -- >>> >> > >> > > > Anders >>> >> > >> > > > >>> >> > >> > > > >>> >> > >> > > > >>> >> > >> > > > tis 23 juni 2015 kl 10:30 skrev Garth N. Wells >>> >> > >> > > > <[email protected]>: >>> >> > >> > > > >>> >> > >> > > > > On 22 June 2015 at 18:00, Anders Logg >>> >> > >> > > > > <[email protected]> wrote: >>> >> > >> > > > > > The challenge in moving the bounding box tree outside >>> >> > >> > > > > > of >>> >> > >> > > > > > the mesh is >>> >> > >> > > > > that it >>> >> > >> > > > > > has always been part of the mesh (it replaced an >>> >> > >> > > > > > earlier >>> >> > >> > > > > > data structure >>> >> > >> > > > > that >>> >> > >> > > > > > was there) so a user expects to be able to do >>> >> > >> > > > > > >>> >> > >> > > > > > v = Function(V) >>> >> > >> > > > > > print v(x) >>> >> > >> > > > > > >>> >> > >> > > > > >>> >> > >> > > > > This is fine for Expressions, but for a Function I don't >>> >> > >> > > > > think it's bad for the interface to make obvious to the >>> >> > >> > > > > user that they are performing a potentially expensive >>> >> > >> > > > > operation. If the user was required to pass the cell, all >>> >> > >> > > > > would be fine. It would also fix the issues with Function >>> >> > >> > > > > evaluations in parallel. >>> >> > >> > > > > >>> >> > >> > > > > > without needing to instantiate some cryptic >>> >> > >> > > > > > BoundingBoxTree data >>> >> > >> > > > > structure. >>> >> > >> > > > > > Furthermore, a user expects that on subsequent calls >>> >> > >> > > > > > v(x) is fast since >>> >> > >> > > > > the >>> >> > >> > > > > > tree has already been built. >>> >> > >> > > > > > >>> >> > >> > > > > > I don't see a way around automatic handling of building >>> >> > >> > > > > > the search tree. >>> >> > >> > > > > Are >>> >> > >> > > > > > there some clever suggestions? >>> >> > >> > > > > > >>> >> > >> > > > > >>> >> > >> > > > > We have a fundamental problem/flaw that MeshGeometery is >>> >> > >> > > > > mutable and a Mesh owns a bounding box object. One of the >>> >> > >> > > > > two needs to give. >>> >> > >> > > > > >>> >> > >> > > > > Garth >>> >> > >> > > > > >>> >> > >> > > > > >>> >> > >> > > > > > Handling this for PointSource and MultiMesh is >>> >> > >> > > > > > unproblematic. (But for MultiMesh, I would rather want >>> >> > >> > > > > > to move it myself.) >>> >> > >> > > > > > >>> >> > >> > > > > > -- >>> >> > >> > > > > > Anders >>> >> > >> > > > > > >>> >> > >> > > > > > >>> >> > >> > > > > > mån 22 juni 2015 kl 17:54 skrev Marco Morandini < >>> >> > >> > > > > [email protected]>: >>> >> > >> > > > > >> >>> >> > >> > > > > >> >>> Besides, the mesh bounding_box_tree used to find >>> >> > >> > > > > >> >>> the >>> >> > >> > > > > >> >>> colliding mesh entity is cached. I fear this could >>> >> > >> > > > > >> >>> be a source of "strange" results, because its use >>> >> > >> > > > > >> >>> is here completely transparent to the user, who >>> >> > >> > > > > >> >>> may >>> >> > >> > > > > >> >>> be unaware of the need to update it. >>> >> > >> > > > > >> >>> >>> >> > >> > > > > >> >> >>> >> > >> > > > > >> >> I really don't like magical caching. How about >>> >> > >> > > > > >> >> having >>> >> > >> > > > > >> >> >>> >> > >> > > > > >> >> PointSource::apply(GenericVector& b, Cell& c, >>> >> > >> > > > > >> >> double >>> >> > >> > > > > >> >> magnitude); >>> >> > >> > > > > >> >> >>> >> > >> > > > > >> >> The user is responsible for finding the cell, and >>> >> > >> > > > > >> >> thereby also responsible for handling meshes that >>> >> > >> > > > > >> >> move, etc. >>> >> > >> > > > > >> >> >>> >> > >> > > > > >> >> PointSource::apply(...) presently uses >>> >> > >> > > > > >> >> Mesh::bounding_box_tree(), which I would like to >>> >> > >> > > > > >> >> get >>> >> > >> > > > > >> >> rid of from Mesh since mesh geometry is mutable. If >>> >> > >> > > > > >> >> the search tools are not cached, the user takes >>> >> > >> > > > > >> >> responsibility for managing the bounding boxes. >>> >> > >> > > > > >> > >>> >> > >> > > > > >> > For the record, this issue is filed as >>> >> > >> > > > > >> > https://bitbucket.org/fenics-project/dolfin/issue/89 >>> >> > >> > > > > >> > >>> >> > >> > > > > >> >>> >> > >> > > > > >> Right now Mesh::bounding_box_tree() is used by >>> >> > >> > > > > >> >>> >> > >> > > > > >> void PointSource::apply(GenericVector& b) >>> >> > >> > > > > >> void Function::eval(Array<double>& values, const >>> >> > >> > > > > >> Array<double>& x) const >>> >> > >> > > > > >> >>> >> > >> > > > > >> and >>> >> > >> > > > > >> >>> >> > >> > > > > >> MultiMesh ( + MultiMeshDirichletBC ) >>> >> > >> > > > > >> >>> >> > >> > > > > >> It would be pretty easy to change PointSource and >>> >> > >> > > > > >> Function. But I think that, for MultiMesh, I should >>> >> > >> > > > > >> move the bboxes there, and leave MultiMesh::build() >>> >> > >> > > > > >> unchanged. >>> >> > >> > > > > >> >>> >> > >> > > > > >> I can go this route if there is consensus. >>> >> > >> > > > > >> >>> >> > >> > > > > >> Marco >>> >> > >> > > > > >> >>> >> > >> > > > > >> >>> >> > >> > > > > >> >>> >> > >> > > > > >> _______________________________________________ >>> >> > >> > > > > >> fenics mailing list >>> >> > >> > > > > >> [email protected] >>> >> > >> > > > > >> http://fenicsproject.org/mailman/listinfo/fenics >>> >> > >> > > > > _______________________________________________ >>> >> > >> > > > > fenics mailing list >>> >> > >> > > > > [email protected] >>> >> > >> > > > > http://fenicsproject.org/mailman/listinfo/fenics >>> >> > >> > > > > >>> >> > >> > > >>> >> > >> > > _______________________________________________ >>> >> > >> > > fenics mailing list >>> >> > >> > > [email protected] >>> >> > >> > > http://fenicsproject.org/mailman/listinfo/fenics >>> >> > >> > > >>> >> > >> >>> >> > >> _______________________________________________ >>> >> > >> fenics mailing list >>> >> > >> [email protected] >>> >> > >> http://fenicsproject.org/mailman/listinfo/fenics >>> >> > > >>> >> > > >>> >> > > >>> >> > > _______________________________________________ >>> >> > > fenics mailing list >>> >> > > [email protected] >>> >> > > http://fenicsproject.org/mailman/listinfo/fenics >>> >> > > >>> >> > _______________________________________________ >>> >> > fenics mailing list >>> >> > [email protected] >>> >> > http://fenicsproject.org/mailman/listinfo/fenics >>> >> >>> > >>> _______________________________________________ >>> fenics mailing list >>> [email protected] >>> http://fenicsproject.org/mailman/listinfo/fenics > > _______________________________________________ fenics mailing list [email protected] http://fenicsproject.org/mailman/listinfo/fenics
