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