dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #24424
Re: Exposing MeshMarkers in Python
On Tuesday September 13 2011 10:21:24 Anders Logg wrote:
> On Tue, Sep 13, 2011 at 10:18:58AM -0700, Johan Hake wrote:
> > On Tuesday September 6 2011 14:01:10 Anders Logg wrote:
> > > On Tue, Sep 06, 2011 at 01:57:20PM -0700, Johan Hake wrote:
> > > > On Tuesday September 6 2011 13:47:02 Anders Logg wrote:
> > > > > On Tue, Sep 06, 2011 at 12:21:03PM -0700, Johan Hake wrote:
> > > > > > On Tuesday September 6 2011 12:18:33 Anders Logg wrote:
> > > > > > > On Tue, Sep 06, 2011 at 11:31:03AM -0700, Johan Hake wrote:
> > > > > > > > On Tuesday September 6 2011 11:04:23 Anders Logg wrote:
> > > > > > > > > On Tue, Sep 06, 2011 at 05:45:33PM +0100, Garth N. Wells
wrote:
> > > > > > > > > > On 6 September 2011 17:31, Johan Hake
> > > > > > > > > > <johan.hake@xxxxxxxxx>
> > > >
> > > > wrote:
> > > > > > > > > > > On Monday September 5 2011 00:09:58 Anders Logg wrote:
> > > > > > > > > > >> On Sun, Sep 04, 2011 at 11:23:04PM -0700, Johan Hake
wrote:
> > > > > > > > > > >> > On Friday September 2 2011 23:19:22 Anders Logg
wrote:
> > > > > > > > > > >> > > On Fri, Sep 02, 2011 at 02:35:57PM -0700, Johan
> > > > > > > > > > >> > > Hake
> >
> > wrote:
> > > > > > > > > > >> > > > What is the different between a MeshMarker and a
> > > > > > > > > > >> > > > MeshFunction? Is MeshMarker a MeshFunction but
> > > > > > > > > > >> > > > instead of storing the values in line with its
> > > > > > > > > > >> > > > global entity index it stores it wrt the global
> > > > > > > > > > >> > > > cell entity index together with its local entity
> > > > > > > > > > >> > > > index?
> > > > > > > > > > >> > >
> > > > > > > > > > >> > > Yes, that and values don't need to be stored on
> > > > > > > > > > >> > > the entire mesh, only for a subset, so you can
> > > > > > > > > > >> > > mark just 3 facets without needing to store
> > > > > > > > > > >> > > markers for a million facets.
> > > > > > > > > > >> >
> > > > > > > > > > >> > ok, I will see what I can do.
> > > > > > > > > > >>
> > > > > > > > > > >> Thanks!
> > > > > > > > > > >>
> > > > > > > > > > >> > > Copy paste from the MeshMarker docstring:
> > > > > > > > > > >> > > /// The MeshMarkers class can be used to store
> > > > > > > > > > >> > > data associated
> > > > > > > > > > >> > >
> > > > > > > > > > >> > > with /// a subset of the entities of a mesh of a
> > > > > > > > > > >> > > given topological /// dimension. It differs from
> > > > > > > > > > >> > > the MeshFunction class in two ways. /// First,
> > > > > > > > > > >> > > data does not need to be associated with all
> > > > > > > > > > >> > > entities /// (only a subset). Second, data is
> > > > > > > > > > >> > > associated with entities /// through the
> > > > > > > > > > >> > > corresponding cell index and local entity number
> > > > > > > > > > >> > > /// (relative to the cell), not by global entity
> > > > > > > > > > >> > > index, which means /// that data may be stored
> > > > > > > > > > >> > > robustly to file.
> > > > > > > > > > >> > >
> > > > > > > > > > >> > > > Also, will this take over for the way we use
> > > > > > > > > > >> > > > MeshFunctions in the assembler, or will a
> > > > > > > > > > >> > > > MeshFunction be generated by a MeshMarker before
> > > > > > > > > > >> > > > assemble gets called?
> > > > > > > > > > >> > >
> > > > > > > > > > >> > > I think we will do that as a first step (convert
> > > > > > > > > > >> > > from MeshMarker to MeshFunction) since then we
> > > > > > > > > > >> > > don't need to touch the assembler. Then later we
> > > > > > > > > > >> > > can think about using MeshMarkers directly.
> > > > > > > > > > >> >
> > > > > > > > > > >> > Ok.
> > > > > > > > > > >> >
> > > > > > > > > > >> > > > I think I also get confused with the naming
> > > > > > > > > > >> > > > here. If my explaination of what MeshMarker is
> > > > > > > > > > >> > > > doing is correct, a MeshMarker and a
> > > > > > > > > > >> > > > MeshFunction are essentially doing the same
> > > > > > > > > > >> > > > thing. What differs is the way the data is
> > > > > > > > > > >> > > > stored. This is not reflected in the naming of
> > > > > > > > > > >> > > > the classes
> > > > > > > > > > >> > >
> > > > > > > > > > >> > > It was the best I could come up with. Feel free to
> > > > > > > > > > >> > > suggest something else. SubsetMeshFunction would
> > > > > > > > > > >> > > also be confusing since it's not really a
> > > > > > > > > > >> > > MeshFunction.
> > > > > > > > > > >> > >
> > > > > > > > > > >> > > Either way, I expect the MeshMarkers class to be
> > > > > > > > > > >> > > used mostly internally by the MeshDomains class.
> > > > > > > > > > >> >
> > > > > > > > > > >> > Ok.
> > > > > > > > > > >> >
> > > > > > > > > > >> > Not sure these are better, but they might reflect
> > > > > > > > > > >> > the difference between this guy and a MeshFunction
> > > > > > > > > > >> > in a slightly more intuitive way.
> > > > > > > > > > >> >
> > > > > > > > > > >> > MeshEntityFunction, LocalMeshEntityFunction,
> > > > > > > > > > >> > LocalMeshFunction, SubMeshFunction
> > > > > > > > > > >>
> > > > > > > > > > >> I'm not sure those are much better, and I don't think
> > > > > > > > > > >> it would be correct to call them something containing
> > > > > > > > > > >> "Function" since they are not really functions. With
> > > > > > > > > > >> a MeshFunction, one can input x (a mesh entity) and
> > > > > > > > > > >> get y = f(x) (the value of the MeshFunction at that
> > > > > > > > > > >> entity). That's not possible with MeshMarkers; they
> > > > > > > > > > >> are just a collection of markers, not really a
> > > > > > > > > > >> function since the value is only defined on a subset
> > > > > > > > > > >> and one would need to loop through the list of values
> > > > > > > > > > >> to get the value at any entity where the value is
> > > > > > > > > > >> actually defined.
> > > > > > > > > > >
> > > > > > > > > > > What with MeshValueCollection? As it is a templated
> > > > > > > > > > > class I do not think Marker is an appropriated name.
> > > > > > > > > >
> > > > > > > > > > Agree.
> > > > > > > > > >
> > > > > > > > > > > 'Collection' says that the class is not
> > > > > > > > > > > defined over the whole Mesh.
> > > > > > > > >
> > > > > > > > > I don't see what the templating has to do with the name
> > > > > > > > > "markers" but MeshValueCollection sounds good.
> > > > > > > > >
> > > > > > > > > > > Two questions:
> > > > > > > > > > >
> > > > > > > > > > > How can the following code work:
> > > > > > > > > > > // Get marker data
> > > > > > > > > > > const std::vector<uint>& marker = _markers[i];
> > > > > > > > > > > const uint cell_index = marker[0];
> > > > > > > > > > > const uint local_entity = marker[1];
> > > > > > > > > > > const T marker_value = marker[2];
> > > > > > > > > > >
> > > > > > > > > > > when _markers is declared as:
> > > > > > > > > > > // The markers
> > > > > > > > > > > std::vector<std::pair<std::pair<uint, uint>, T> >
> > > > > > > > > > > _markers;
> > > > > > > > >
> > > > > > > > > The above code doesn't work. I suspect the code hasn't yet
> > > > > > > > > been instantiated so it wasn't detected by the compiler.
> > > > > > > > >
> > > > > > > > > The markers need to be accessed as follows (from
> >
> > XMLMeshMarkers.h):
> > > > > > > > > for (uint i = 0; i < mesh_markers.size(); ++i)
> > > > > > > > > {
> > > > > > > > >
> > > > > > > > > pugi::xml_node entity_node =
> > > > > > > > > mf_node.append_child("marker"); const
> > > > > > > > > std::pair<std::pair<uint, uint>, T>& marker =
> > > > > > > > >
> > > > > > > > > mesh_markers.get_marker(i);
> > > > > > > > >
> > > > > > > > > entity_node.append_attribute("cell_index") =
> > > > > > > > > marker.first.first;
> > > > > > > > > entity_node.append_attribute("local_entity") =
> > > > > > > > > marker.first.second;
> > > > > > > > > entity_node.append_attribute("marker_value") =
> > > > > > > > > marker.second;
> > > > > > > > >
> > > > > > > > > }
> > > > > > > > >
> > > > > > > > > > The above also permits multiple entries. Perhaps we want
> > > > > > > > > >
> > > > > > > > > > boost::unordered_map<std::pair<std::pair<uint, uint>,
> > > > > > > > > > T>
> > > > > > > > > >
> > > > > > > > > > > _markers;
> > > > > > > > >
> > > > > > > > > Yes, maybe but I'm not sure what the cost would be for the
> > > > > > > > > lookup on each cell during assembly.
> > > > > > > > >
> > > > > > > > > > > What is the logic behind:
> > > > > > > > > > > // Set all value of mesh function to maximum value
> > > > > > > > > > > (not all will // be set) by markers below
> > > > > > > > > > > mesh_function.set_all(maxval);
> > > > > > > > > > >
> > > > > > > > > > > Isn't it more natural to initiate the values to zero?
> > > > > > > > > > > Also it makes no sense in conjunction with boundary
> > > > > > > > > > > markers. Then all boundary faces gets marked with the
> > > > > > > > > > > largest marker value. I cannot see how that could be
> > > > > > > > > > > correct.
> > > > > > > > > >
> > > > > > > > > > I don't get ' mesh_function.set_all(maxval);' or the code
> > > > > > > > > > comment.
> > > > > > > > >
> > > > > > > > > The point is that one should be able to define a form with
> > > > > > > > > domains say dx(0), dx(1) and dx(2) and then have a mesh
> > > > > > > > > file that marks a subset of the cells with '0', '1' and
> > > > > > > > > '2'.
> > > > > > > > >
> > > > > > > > > Then the conversion to MeshFunction inserts '3' for all
> > > > > > > > > other (unmarked) cells. This allows a user to specify only
> > > > > > > > > the interesting cells and no need to mark the rest with -1
> > > > > > > > > or None or similar.
> > > > > > > >
> > > > > > > > That would make sense if the code would be:
> > > > > > > > mesh_function.set_all(maxval+1);
> > > > > > >
> > > > > > > Yes, that is the intention! Thanks for proofreading my code
> > > > > > > before I've even had a chance to test it. :-)
> > > > > > :
> > > > > > :)
> > > > > >
> > > > > > You gave the impression that it was test when you asked me to
> > > > > > wrap it to Python. Give me a ping when it is ready and I will
> > > > > > have a bite at the SWIG code.
> > > > >
> > > > > I wanted the Python wrappers so that I could write the unit tests
> > > > > for it (in Python) and then find the bugs... :-)
> > > >
> > > > Ahhh!
> > > >
> > > > You fix the backward compatability of the file format of the
> > > > MeshFunction and I start on the SWIG code. Deal?
> > >
> > > Deal! But not until after tomorrow. I have a paper I need to revise.
> >
> > Where should I push my fix? To
> >
> > lp:~dolfin-core/dolfin/logg
> >
> > ?
>
> I've merged with trunk so you can just push there.
Where should I push. That was just slightly ambigous :)
Johan
References