dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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