On Mon, Jan 13, 2014 at 12:45:11PM +0000, Garth N. Wells wrote: > I've just pushed some changes to master, which for > > from dolfin import * > mesh = UnitCubeMesh(128, 128, 128) > mesh.init(2) > mesh.init(2, 3) g> > give a factor 2 reduction in memory usage and a factor 2 speedup. > Change is at > https://bitbucket.org/fenics-project/dolfin/commits/8265008. > > The improvement is primarily due to d--d connectivity no being > computed. I'll submit a pull request to throw an error when d--d > connectivity is requested. The only remaining place where d--d is > (implicitly) computed is in the code for finding constrained mesh > entities (e.g., for periodic bcs). The code in question is > > for (MeshEntityIterator(facet, dim) e; . . . .) > > when dim is the same as the topological dimension if the facet. As > in other places, it would be consistent (and the original intention > of the programmer) if this just iterated over the facet itself > rather than all facets attached to it via a vertex.
I don't see why an error message is needed. Could we not just add the possibility to specify what d -- d means? It might be useful for other algorithms to be able to compute that data. One option could be to throw an error message when d -- d connectivity has not been specified. -- Anders > Garth > > > On 2014-01-08 17:00, Garth N. Wells wrote: > >On 2014-01-08 16:48, Anders Logg wrote: > >>On Wed, Jan 08, 2014 at 04:41:17PM +0000, Garth N. Wells wrote: > >>>On 2014-01-08 16:17, Anders Logg wrote: > >>>>On Wed, Jan 08, 2014 at 04:07:35PM +0000, Garth N. Wells wrote: > >>>>>On 2014-01-08 15:21, Anders Logg wrote: > >>>>>>I would claim that the connectivity is stored very efficiently. The > >>>>>>simplest example is D--0 connectivity for a mesh which is for each > >>>>>>cell which vertices belong to it. That data is stored as one big array > >>>>>>of integers, something like > >>>>>> > >>>>>> 0 1 2 3 0 1 2 4 ... > >>>>>> > >>>>>>which means that tetrahedron 0 has vertices 0 1 2 3, tetrahedron 1 has > >>>>>>vertices 0 1 2 4 etc. > >>>>>> > >>>>>>There are D + 1 different kinds of entities in a mesh. For a > >>>>>>tetrahedron we have vertex, edge, face, cell and so the total number > >>>>>>of connectivities that may be computed is (D + 1)*(D + 1). Not all > >>>>>>these are needed, but they can all be computed. > >>>>>> > >>>>>>For solving Poisson's equation, we first need D--0. To see which > >>>>>>connectivities are computed, one can do info(mesh, True) to print the > >>>>>>following little table: > >>>>>> > >>>>>> 0 1 2 3 > >>>>>> 0 - - - - > >>>>>> 1 - - - - > >>>>>> 2 - - - - > >>>>>> 3 x - - - > >>>>>> > >>>>>>To set boundary conditions, one needs to create the faces of the mesh > >>>>>>(which are not stored a priori). One will then get the following > >>>>>>connectivities: > >>>>>> > >>>>>> 0 1 2 3 > >>>>>> 0 - - - x > >>>>>> 1 - - - - > >>>>>> 2 x - - - > >>>>>> 3 x - x x > >>>>>> > >>>>>>Here we see the following connectivities: > >>>>>> > >>>>>> 3--0 = the vertices of all cells = the cells > >>>>>> 0--3 = the cells connected to all vertices > >>>>>> 3--3 = the cells connected to all cells (via vertices) > >>>>>> 2--0 = the vertices of all faces = the faces > >>>>>> 3--2 = the faces of all cells > >>>>>> > >>>>>>These get computed in order: > >>>>>> > >>>>>> 0--3 is computed from 3--0 (by "inversion") > >>>>>> 3--3 is computed from 3--0 and 0--3 > >>>>>> 2--0 and 3--2 are computed from 3--3 by a local search > >>>>>> > >>>>>>3--3 is needed in order for a cell to communicate with its neighboring > >>>>>>cells to decide on how to number the common face (if any). > >>>>>> > >>>>> > >>>>>3--3 connectivity is *undefined*. > >>>> > >>>>Yes, in general, but we have *chosen* to define it as 3--0--3. > >>>> > >>> > >>>Which I think we should attempt to remove > >>> > >>>>The reason for this choice is that we start with 3--0, from which we > >>>>can easily compute 3--0--3. > >>>> > >>>>In other words, we start with cell-vertex data and hence the only way > >>>>for two cells to figure out if they share common entities (for example > >>>>faces) is to go via the vertices. Unless you have some other brilliant > >>>>idea. > >>>> > >>> > >>>The function GraphBuilder::compute_local_dual_graph does this. It > >>>loops over cells and builds a hashed map (boost::unordered_map) from > >>>a (sorted) vector of the facet vertices to the cell index. If the > >>>key (the facet vertices) is already in the map, then we know a cell > >>>that the current cell shares a facet with. If the key is not already > >>>in the map, then we add it. With appropriate hashing, the > >>>unordered_map look-ups are O(1). > >> > >>Now that's a brilliant idea. Why is this not already the default? > >> > > > >Reluctance to touch the beautifully crafted mesh library? :) I'll add > >it alongside the current code so that we can test it for speed. > > > >It's used in GraphBuilder to construct a dual graph (cell-facet-cell > >connections) to feed to a graph partitioner because in GraphBuilder we > >don't yet have a Mesh object - just cell connectivity data. > > > >>Application of boundary conditions can bypass TopologyComputation.cpp > >>and build the facets using the dual graph and it would not be in > >>conflict with TopologyComputation.cpp (which may still be useful for > >>computing other connectitivies). > >> > > > >The approach in GraphBuilder::compute_local_dual_graph could be used > >for other connection types. > > > >Garth > > > > > >_______________________________________________ > >fenics mailing list > >[email protected] > >http://fenicsproject.org/mailman/listinfo/fenics _______________________________________________ fenics mailing list [email protected] http://fenicsproject.org/mailman/listinfo/fenics
