dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Tuesday September 13 2011 11:31:38 Anders Logg wrote:
> On Tue, Sep 13, 2011 at 11:25:00AM -0700, Johan Hake wrote:
> > On Tuesday September 13 2011 11:15:06 Anders Logg wrote:
> > > On Tue, Sep 13, 2011 at 10:51:28AM -0700, Johan Hake wrote:
> > > > 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.
> > > > 
> > > > Pushed it to trunk. You need to enable the unit test and fix some
> > > > errors which is throwned from the C++ code. :) Also see a README in
> > > > MeshValueCollection.h.
> > > 
> > > Thanks! That was fast! I will look at it tomorrow.
> > 
> > Ctrl-w Ctrl-y + Meta-% MeshFunction MeshValueCollection
> > 
> > Is bleeding fast!
> 
> I thought about it but I didn't dare... :-)

May the SWIG force be with you!

Johan