dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Oct 17, 2006 at 02:35:16PM +0200, Garth N. Wells wrote:
> OK. It's taking me a while to get into the new structure.

It should be simple and intuitive, so if it's not I might have
overlooked something. Any suggestions for improvements are welcome.

Maybe it's just that it's (a little) different and more general than
the old mesh library but maybe there is room for improvement.

> So, for a 2D mesh, mesh.init(1, 2) computes which cells share an
> edge

Yes.

> and mesh.init(2, 1) computes which edges belong to a cell?

Yes.

The connectivity 1 - 2 is actually computed as the transpose of 2 - 1
so what happens when 1 - 2 is computed is that 2 - 1 is first computed
and 1 - 2 follows from that.

This computation happens in the class TopologyComputation. For a mesh
of topological dimension d, there are (d+1)^2 different things to
compute and there is some logic that takes care of computing one from
the other. By doing things in the right order, everything is O(N) and
should be quite fast. (There are some hidden factors for local
iterations in there as well.)

> Inside
> mesh.init(2, 1) would therefore be the place to compute local edge ID's?

I don't think you need to modify the computation (if there is no
bug). Everything should be available. I'm not sure what you mean by
local edge IDs but I assume you mean the local numbering of edges
within a cell. Then you just need to do mesh.init(2, 1) for a 2D mesh
and you have it. For any cell in that mesh, you may then do

    uint* edges = cell.connections(1);

Which gives you the global edge numbers for the three edges of the
cell. The local numbers are 0, 1, 2.

Is that what you mean?

/Anders



> Garth
> 
> 
> Anders Logg wrote:
> > On Tue, Oct 17, 2006 at 12:16:40PM +0200, Garth N. Wells wrote:
> >> I'm running into a few problems trying to understand the new mesh
> >> structure. If I do
> >>
> >>     UnitSquare mesh(1,1);
> >>
> >>     for (FacetIterator facet(mesh); !facet.end(); ++facet)
> >>       cout << "Number of connected cell "
> >>           << (*facet).numConnections(mesh.topology().dim()) << endl;
> >>
> >> I get the error
> >>
> >>     *** Assertion (entity < num_entities) failed
> > 
> > This is correct behaviour, but perhaps the error message should be
> > more informative.
> > 
> > When you read a mesh from file or create a simple mesh like a mesh of
> > the unit square or unit cube, the only connectivity that exists is
> > 
> >     mesh -> vertices
> >     mesh -> cells
> > 
> >     cell -> vertices
> > 
> > In the above example, you use the connectivity
> > 
> >     cell -> facets (edges in this case)
> > 
> > This is possible since the iterators automatically generate the
> > requested connectivity the first time they are used. So when you
> > create a FacetIterator over the mesh, the facets of the mesh are
> > created. 
> 
> The next time you iterate over facets, the facets don't need
> > to be created.
> > 
> > But then the problem comes when you want to use the connectivity
> > 
> >     facet -> cells
> > 
> > This does not exist and you don't create an iterator that
> > automatically computes it. So in this case you need to build the
> > connectivity. This can be done by
> > 
> >     mesh.init(1, 2); // facet --> cell
> > 
> > If you call this first, then the above example should work.
> > 
> > The reason this worked for the old mesh is that then all the
> > connectivity was computed at startup. Now, only the connectivity that
> > is actually used is computed.
> > 
> > Any ideas for how to make this more user-friendly are welcome. My idea
> > is that as long as you work with the iterators, everything gets
> > generated as needed and otherwise you need to explicitly compute the
> > data you need.
> > 
> >> but it I create a boundary mesh first,
> >>
> >>     UnitSquare mesh(1,1);
> >>     MeshFunction<dolfin::uint> vertex_map;
> >>     MeshFunction<dolfin::uint> cell_map;
> >>     BoundaryMesh boundary(mesh, vertex_map, cell_map);
> >>
> >>     for (FacetIterator facet(mesh); !facet.end(); ++facet)
> >>       cout << "Number of connected cell "
> >>           << (*facet).numConnections(mesh.topology().dim()) << endl;
> >>
> >> everything works as expected. I guess the mesh is not building the list
> >> of facets, and this is being done when creating a boundary mesh. Is this
> >> a bug?
> >>
> >> Garth
> > 
> > This works because to compute the boundary, the connectivity facet -> cells
> > is needed. So during the computation of the boundary, this gets
> > computed and then exists when you ask for it.
> > 
> > The above example should work also with just a simple
> > 
> >   BoundaryMesh boundary(mesh);
> > 
> > The vertex and cell maps are don't needed for this.
> > 
> > /Anders
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/mailman/listinfo/dolfin-dev
> > 
> 
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev