On Mon, 2 Sep 2013 13:20:46 +0200 Martin Sandve Alnæs <[email protected]> wrote:
> Some ideas then: > > File formats for MeshFunction should be used for storage of > mesh-associated data. Conversion tools relate to the file format, not > to arbitrary reordering algorithms in dolfin. In principle you're right. Problem is that 1. there is an issue https://bugs.launchpad.net/dolfin/+bug/1081635. 2. mechanism of reordering user's data not supplied within mesh (but still dependent on mesh ordering) is not available > > Dolfin needs functionality for conversion between MeshFunction<float> > and Function, the most typical cases are cell and vertex data: Yeah and this is what `vertex_to_dof_map` is for. There should be also `cell_to_dof_map` for DG spaces. > > V = CG1 space > u = Function(V) > m = MeshFunction("double", mesh, 0) > u.assign(m) # Only valid if dofs are vertex associated > m.assign(u) # Ditto Do you think that you can handle all the combination of function spaces and entity functions robustly? Maybe you finally end up doing projection or interpolation which is what one wanted to avoid so introduced `vertex_to_dof_map`. Jan > > V = DG0 space > u = Function(V) > m = MeshFunction("double", mesh, dim) > u.assign(m) # Only valid if dofs are cell associated > m.assign(u) # Ditto > > Although "double" is not sufficient for vector/tensor valued > functions. > > Martin > > > > On 2 September 2013 13:09, Johan Hake <[email protected]> wrote: > > > On Monday September 2 2013 12:07:15 Jan Blechta wrote: > > > On Mon, 02 Sep 2013 11:52:55 +0200 > > > > > > Johan Hake <[email protected]> wrote: > > > > On Monday September 2 2013 11:45:21 Jan Blechta wrote: > > > > > On Mon, 02 Sep 2013 11:36:28 +0200 > > > > > > > > > > Johan Hake <[email protected]> wrote: > > > > > > On Monday September 2 2013 10:09:56 Garth N. Wells wrote: > > > > > > > On 2 September 2013 09:58, Johan Hake > > > > > > > <[email protected]> > > > > > > > > > > > > > > wrote: > > > > > > > > On Monday September 2 2013 09:42:49 Garth N. Wells > > > > > > > > wrote: > > > > > > > >> On 2 September 2013 09:30, Johan Hake > > > > > > > >> <[email protected]> > > > > > > > >> > > > > > > > >> wrote: > > > > > > > >> > Seems like updating the doc string wont help as > > > > > > > >> > enough people have tried > > > > > > > >> > to > > > > > > > >> > use the vertex_to_dof_map and failed. > > > > > > > >> > > > > > > > > >> > I agree that the left to right reading does not > > > > > > > >> > apply to the example Garth > > > > > > > >> > presented. If that is the expected behavior, and I > > > > > > > >> > guess it is given the > > > > > > > >> > comments in this treahd, we should just rename the > > > > > > > >> > methods. That would generalize the methods and > > > > > > > >> > probably fit better into the general interface > > > > > > > >> > of DofMap. > > > > > > > >> > > > > > > > > >> > However that would limit the scope of the map and > > > > > > > >> > remove one important motivation for adding the map > > > > > > > >> > in the first place, namely to turn general > > > > > > > >> > > > > > > > >> > vector function values ordered as: > > > > > > > >> How can renaming limit scope? The functionality > > > > > > > >> remains the same. > > > > > > > > > > > > > > > > Because it is not enough to just rename it. We also > > > > > > > > need to remove the functionality for vector function > > > > > > > > spaces. Your example does not make sense > > > > > > > > if> > > > > > > > > > > > > > > > > you change: > > > > > > > > V = FunctionSpace(mesh, "Lagrange", 1) > > > > > > > > > > > > > > > > to > > > > > > > > > > > > > > > > V = VectorFunctionSpace(mesh, "Lagrange", 1) > > > > > > > > > > > > > > My example was deliberately simple. The functionality can > > > > > > > be retained. > > > > > > > > > > > > How? With a map of maps, where the second is the local dof > > > > > > index? > > > > > > > > > > > > > >> > vertex_index*dofs_per_vertex+local_dof > > > > > > > >> > > > > > > > > >> > to an array which could be feed directly into a > > > > > > > >> > vector of a Function in a > > > > > > > >> > VectorFunctionSpace (or similar mixed CG1 function > > > > > > > >> > spaces). The present functionality also works for > > > > > > > >> > parallel runs, as seen by the following example: > > > > > > > >> > > > > > > > > >> > mpirunt -np 2 python vertex_to_dofs.py > > > > > > > >> > > > > > > > > >> > # vertex_to_dofs.py > > > > > > > >> > > > > > > > > >> > from dolfin import * > > > > > > > >> > import numpy as np > > > > > > > >> > > > > > > > > >> > mesh = UnitSquareMesh(20,20) > > > > > > > >> > V = VectorFunctionSpace(mesh, "CG", 1) > > > > > > > >> > u = Function(V) > > > > > > > >> > vertex_to_dof_map = > > > > > > > >> > V.dofmap().vertex_to_dof_map(mesh) > > > > > > > >> > > > > > > > > >> > data = np.reshape(mesh.coordinates()[:], > > > > > > > >> > (mesh.num_vertices()*2)) > > > > > > > >> > > > > > > > >> This is problematic - it makes an assumption of the > > > > > > > >> ordering in mesh.coordinates(). > > > > > > > > > > > > > > > > The only assumption is that you have some data (possible > > > > > > > > vector or tensor data) which are ordered based on the > > > > > > > > mesh (vertices). > > > > > > > > > > > > > > > >> I have seen that a good re-ordering of mesh data > > > > > > > >> can give up to a 50% speed up for assembly, and which > > > > > > > >> will be added in the future. We should not be exposing > > > > > > > >> low-level storage. > > > > > > > > > > > > > > > > Not sure what you mean. This has nothing to do with > > > > > > > > assemble. Only transferring vertex based data into a > > > > > > > > Function. > > > > > > > > > > > > > > Exposing low level storage (e.g. (mesh.coordinates()), > > > > > > > violates data hiding, which then can affect all parts of > > > > > > > a code. If the mesh data ordering is changed, say to make > > > > > > > assembler faster, your example code will likely break. > > > > > > > > > > > > Using mesh.coordinates() was just an example on some data > > > > > > which comes together with the mesh. I just high-jacked > > > > > > coordinates to represent some vector field aligned with the > > > > > > mesh. Real data often comes aligned with the mesh and we > > > > > > need a simple and intuitive way to read such data into a > > > > > > Function. This is basic functionality alot of users need > > > > > > for their applications. > > > > > > > > > > > > If you intend to include some mesh reordering, I suggest > > > > > > that you also include some mapping that brings mesh data to > > > > > > reordered mesh data, and then we need a map to bring > > > > > > reordered mesh data to dof ordering. > > > > > > > > > > Why? This is just changing vertex indices. Currently they're > > > > > already ordered "somehow irregularly". So they will ordered > > > > > in the other way... > > > > > > > > Because the average Joe have a mesh from somewhere together > > > > with some data. If the mesh is reordered by some algorithm we > > > > need to provide a way for Joe to map his data to the new mesh > > > > ordering, so he can call > > > > > > > > u.assign_vertex_data(reordered_vretex_based_data) > > > > > > > > or what not. > > > > > > If I get it correctly, mesh will be reordered once when loaded by > > > DOLFIN. This is already being done and the method how to reorder > > > some data accordingly is not available. > > > > AFAIK the present reordering is limited to shuffling local vertices > > around in > > a cell. This will screw up any facet based mesh values, but all > > vertex and cell based values are preserved. > > > > > For example, there is an open bug on launchpad, that when > > > unordered DOLFIN XML mesh with mesh function is loaded and > > > reordered, mesh function breaks. > > > > > > Loading of mesh functions is currently supported and safe to the > > > extent supported in `meshconvert.py`. > > > > > > Note also that changing vertex indices will require to renumber > > > cell-local vertices and facets according to UFC numbering scheme > > > (=mesh.order()). Hence in addition you will have to provide a way > > > how to keep up-to-date facet based data (FacetFunction), as this > > > is not currently available as I said above. > > > > True. But the ufc ordering is well documented, and implemented in > > several mesh > > conversion methods in meshconvert.py. > > > > I am fine with any solutions that makes it possible for a user to > > read in mesh > > based values and assigning that easily to a Function. This should > > of course work in parallel too. > > > > Johan > > > > > > > > > > > Jan > > > > > > > Johan > > > > > > > > > Jan > > > > > > > > > > > > >> > u.vector().set_local(data[vertex_to_dof_map]) > > > > > > > >> > plot(u, interactive=True) > > > > > > > >> > > > > > > > >> Why not just use Function::compute_vertex_values(...) > > > > > > > >> (plus any necessary generalisation)? > > > > > > > > > > > > > > > > The comparison with compute_vertex_values is > > > > > > > > appropriate. It was raised when we discussed the > > > > > > > > inclusion of the map in the first place. However the > > > > > > > > (present) vertex_to_dof_map give the mapping from > > > > > > > > vertex based data to a Function, where > > > > > > > > compute_vertex_values does the opposite. > > > > > > > > > > > > > > Yes, but two functions were added to GenericDofMap. One > > > > > > > seems to duplicate existing functionality. > > > > > > > > > > > > True. But the two maps, can only be used on data defined on > > > > > > vertices (CG1). compute_vertex_values are more general as it > > > > > > works for Functions on a lot more FunctionSpaces (CG2, DG0, > > > > > > aso)? > > > > > > > > > > > > > > The map is also just > > > > > > > > computed once and can therefore be reused by the user > > > > > > > > if that is needed. > > > > > > > > > > > > > > I don't see the benefit if one can use > > > > > > > Function::compute_vertex_values. > > > > > > > > > > > > See Martin's answer. > > > > > > > > > > > > Johan > > > > > > > > > > > > > Garth > > > > > > > > > > > > > > > Johan > > > > > > > > > > > > > > > >> Garth > > > > > > > >> > > > > > > > >> > Johan > > > > > > > >> > > > > > > > > >> > On Saturday August 31 2013 10:20:21 Simone Pezzuto > > > > > > > >> > wrote: > > > > > > > >> >> Hi, > > > > > > > >> >> > > > > > > > >> >> I'm familiar with these two maps since I > > > > > > > >> >> use > > > > > > > >> >> > > > > > > > >> >> them > > > > > > > >> >> > > > > > > > >> >> for a gradient > > > > > > > >> >> > > > > > > > >> >> recovery technique. > > > > > > > >> >> > > > > > > > >> >> I can assure you that first time I used > > > > > > > >> >> vertex_to_dof_map I was a bit confused, > > > > > > > >> >> since the convention should be left to right (as > > > > > > > >> >> Garth pointed out). > > > > > > > >> >> > > > > > > > >> >> Example: eps2pdf fig.eps ---> fig.pdf > > > > > > > >> >> > > > > > > > >> >> vertex2dof vertex_id --> dof_id > > > > > > > >> >> dof2vertex dof_id --> vertex_id > > > > > > > >> >> > > > > > > > >> >> So at the moment is really confusing. Maybe we can > > > > > > > >> >> introduce new functions > > > > > > > >> >> {vertex2dof,dof2vertex}_map > > > > > > > >> >> (no name collision) and deprecate the old one, so > > > > > > > >> >> the user is aware of the > > > > > > > >> >> change but its code doesn't brake. > > > > > > > >> >> > > > > > > > >> >> Simone > > > > > > > >> >> > > > > > > > >> >> 2013/8/31 Jan Blechta <[email protected]> > > > > > > > >> >> > > > > > > > >> >> > On Fri, 30 Aug 2013 23:47:35 +0100 > > > > > > > >> >> > > > > > > > > >> >> > "Garth N. Wells" <[email protected]> wrote: > > > > > > > >> >> > > On 30 August 2013 23:37, Johan Hake > > > > > > > >> >> > > > > > > > > > >> >> > > <[email protected]> wrote: > > > > > > > >> >> > > > On Friday August 30 2013 23:19:09 Garth N. > > > > > > > >> >> > > > Wells > > > > > > > >> >> > > > > > > > > > > >> >> > > > wrote: > > > > > > > >> >> > > >> On 30 August 2013 22:50, Johan Hake > > > > > > > >> >> > > >> <[email protected]> > > > > > > > > > > > > wrote: > > > > > > > >> >> > > >> > On Friday August 30 2013 15:47:28 Garth N. > > > > > > > >> >> > > >> > Wells > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > wrote: > > > > > > > >> >> > > >> >> The functions > > > > > > > >> >> > > >> >> GenericDofmap::vertex_to_dof_map and > > > > > > > >> >> > > >> >> GenericDofMap::dof_to_vertex_map are not > > > > > > > >> >> > > >> >> properly documented (the doc string is > > > > > > > >> >> > > >> >> the same for both), and I think that they > > > > > > > >> >> > > >> >> are back to front. The docstring in > > > > > > > >> >> > > >> >> DofMap has inconsistencies. I would > > > > > > > >> >> > > >> >> expect that > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> map0 = > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> GenericDofmap::vertex_to_dof_map(...) > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> would mean a map from vertex to dof, i.e. > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> map0[vertex_index] -> dof index > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> and that > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> map1 = > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> GenericDofmap::dof_to_vertex_map(...) > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> would mean a map from dof index to > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> map1[dof_index] -> vertex index > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> Tests (see below code) and the return > > > > > > > >> >> > > >> >> types also indicate that > > > > > > > >> >> > > >> >> things are back to front. Can someone > > > > > > > >> >> > > >> >> clarify the situation? > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > The map was introduced to help a user map > > > > > > > >> >> > > >> > vertex based data onto > > > > > > > >> >> > > >> > a Function.> > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > from dolfin import * > > > > > > > >> >> > > >> > import numpy as np > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > mesh = UnitSquareMesh(20,20) > > > > > > > >> >> > > >> > V = VectorFunctionSpace(mesh, "CG", 1) > > > > > > > >> >> > > >> > u = Function(V) > > > > > > > >> >> > > >> > vertex_to_dof_map = > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > V.dofmap().vertex_to_dof_map(mesh) > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > data = np.reshape(mesh.coordinates()[:], > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > (mesh.num_vertices()*2)) u.vector()[:] = > > > > > > > >> >> > > >> > data[vertex_to_dof_map] > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > plot(u, interactive=True) > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > The size of the data array should be: > > > > > > > >> >> > > >> > mesh.num_vertices()*u.value_size() > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > The documentation should be improved, and > > > > > > > >> >> > > >> > not least properly mapped from C++ to > > > > > > > >> >> > > >> > Python. > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > The name refer to the mapping that turn > > > > > > > >> >> > > >> > vertex based data to dof > > > > > > > >> >> > > >> > based and reads quite well when used as > > > > > > > >> >> > > >> > above. I can see that the word map can be > > > > > > > >> >> > > >> > missleading. It is not a "map" data > > > > > > > >> >> > > >> > structure. It is an index set that "maps > > > > > > > >> >> > > >> > values". > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> > Still confused? > > > > > > > >> >> > > >> > > > > > > > >> >> > > >> I'm not confused. It's clear that the > > > > > > > >> >> > > >> function names are back-to-front. It doesn't > > > > > > > >> >> > > >> matter what they were included for - they > > > > > > > >> >> > > >> are members of GenericDofMap and must make > > > > > > > >> >> > > >> sense in that context. > > > > > > > >> >> > > >> > > > > > > > >> >> > > >> Since reading from left to right is a well > > > > > > > >> >> > > >> established convention, > > > > > > > >> >> > > >> I propose that (a) the function names be > > > > > > > >> >> > > >> fixed by reversing them; > > > > > > > >> >> > > >> and (b) the doc strings be fixed. > > > > > > > >> >> > > > > > > > > > > >> >> > > > Agree on (b). I am not fully convinced by (a). > > > > > > > >> >> > > > > > > > > > > >> >> > > > I am not sure what your example tries to > > > > > > > >> >> > > > show. You are not using the mapping the > > > > > > > >> >> > > > intended way and I am therefore confused > > > > > > > >> >> > > > about the whole back-to-front, front-to-back > > > > > > > >> >> > > > discussion. > > > > > > > >> >> > > > > > > > > > >> >> > > Just read the function names aloud from left to > > > > > > > >> >> > > right > > > > > > > >> >> > > - 'vertex_to_dof_map' should be a 'vertex to dof > > > > > > > >> >> > > map', i.e. a map from > > > > > > > >> >> > > a > > > > > > > >> >> > > vertex *to* a dof. > > > > > > > >> >> > > > > > > > > >> >> > Just read from left to right - 'vertex_to_dof_map' > > > > > > > >> >> > stands for a map which turns a vertex map into a > > > > > > > >> >> > dof map (when used as a right composition). > > > > > > > >> >> > > > > > > > > >> >> > Yes, I was confused at first when I saw this and > > > > > > > >> >> > agree with Garth it should be 'left to right'. > > > > > > > >> >> > But does it worth switching it? Is the whole > > > > > > > >> >> > concept of indexing by > > > > > > > >> >> > > > > > > > > >> >> > vertex_index*dofs_per_vertex+local_dof > > > > > > > >> >> > > > > > > > > >> >> > sustainable? Or should it be replaced by some more > > > > > > > >> >> > robust types which > > > > > > > >> >> > would handle non-injective map (and its > > > > > > > >> >> > inversion)? > > > > > > > >> >> > > > > > > > > >> >> > There were some user codes using these functions > > > > > > > >> >> > as seen in discussions. > > > > > > > >> >> > > > > > > > > >> >> > Jan > > > > > > > >> >> > > > > > > > > >> >> > > Garth > > > > > > > >> >> > > > > > > > > > >> >> > > > Johan > > > > > > > >> >> > > > > > > > > > > >> >> > > >> Garth > > > > > > > >> >> > > >> > > > > > > > >> >> > > >> > Johan > > > > > > > >> >> > > >> > > > > > > > > >> >> > > >> >> Garth > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> from dolfin import * > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> mesh = UnitSquareMesh(4, 4) > > > > > > > >> >> > > >> >> V = FunctionSpace(mesh, "Lagrange", 1) > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> dof_to_vertex = > > > > > > > >> >> > > >> >> V.dofmap().dof_to_vertex_map(mesh) > > > > > > > >> >> > > >> >> vertex_to_dof = > > > > > > > >> >> > > >> >> V.dofmap().vertex_to_dof_map(mesh) > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> for c in cells(mesh): > > > > > > > >> >> > > >> >> print "Cell index:", c.index() > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> # Get cell dofs > > > > > > > >> >> > > >> >> dofs = V.dofmap().cell_dofs(c.index()) > > > > > > > >> >> > > >> >> print " Cell dofs:", dofs > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> # Get vertices from cell > > > > > > > >> >> > > >> >> cell_vertices0 = sorted([v.index() > > > > > > > >> >> > > >> >> for v in vertices(c)]) > > > > > > > >> >> > > >> >> print " Cell vertex indices (from > > > > > > > >> >> > > >> >> cell):", cell_vertices0 > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> # Get vertices from dof_to_vertex > > > > > > > >> >> > > >> >> cell_vertices1 = > > > > > > > >> >> > > >> >> sorted([dof_to_vertex[dof] > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> for > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> dof in > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> dofs]) print " Cell vertex indices (from > > > > > > > >> >> > > >> >> dof_to_vertex_map):", > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> cell_vertices1 > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> # Get vertices from vertex_to_dof_map > > > > > > > >> >> > > >> >> cell_vertices2 = > > > > > > > >> >> > > >> >> sorted([vertex_to_dof[dof] > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> for > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> dof in > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> dofs]) print " Cell vertex indices (from > > > > > > > >> >> > > >> >> vertex_to_dof_map):", > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> cell_vertices2 > > > > > > > >> >> > > >> >> > > > > > > > >> >> > > >> >> _______________________________________________ > > > > > > > >> >> > > >> >> 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 > > > > > > > > _______________________________________________ > > > > 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
