behavior of vertices or mesh.init()

[Former user]

Hello,

I run the following code

from dolfin import *

mesh = UnitSquare(16,16)

mesh.init()

for v1 in vertices(mesh):
    print "vertex", v1.index()
    for v2 in vertices(v1):
          print "neighbors", v2.index()

with fenics-1.0.x, the last version of fenics got from ppa:fenics-packages/fenics and the binaries

In fenics-1.0.x I get the following result
vertex 0
neighbors 1
neighbors 18
neighbors 17
vertex 1
neighbors 0
neighbors 18
neighbors 2
neighbors 19
vertex 2
neighbors 1
neighbors 19
neighbors 3
neighbors 20
vertex 3
neighbors 2
neighbors 20
neighbors 4
neighbors 21

etc

in the last version of fenics for ubuntu and in mac I get

vertex 0
neighbors 0
vertex 1
neighbors 1
vertex 2
neighbors 2
vertex 3
neighbors 3
vertex 4
neighbors 4

[Johan Hake]

This is caused by changed behavior for computation of connectivity between vertices and vertices. I think it solved some weird behavior for 1D computation. See commit:

0098c79218d8e83e92fe0b03c7c2189e7ac4d608

However I agree that computing the connectivity between vertices to work like you expects. Maybe we need to make the if test at line 237 in Topologycomputation.cpp also check for topological dimension of mesh == 1?

Johan

[Fotini Karakatsani]

Hi Johan,

the previous comment was mine. I did not understand. Are you going to correct the connectivity between the vertices? I need to implement something and I need that connectivity. Or should I install a previous version of fenics?

thank you in advance Fotini

[Johan Hake]

Sorry for not responding. I do not get email updates on the issued I do not watch. I am not used to the new system here on Bitbucket...

I will "fix" this. The change was done to fix a previous issue. I am not sure if we need to live with the side effect demonstrated by your issue, or if we could fix this corner case for top dim != 1.

[Prof Garth Wells]

I think the 1.0 behaviour is inconsistent and that the new behaviour is correct.

[Johan Hake]

So you think one should not be able to get vertex neighbors to vertices? I think one should.

I realize that there were some ambiguities for mesh with top dim 1, but I wonder if that introduced a bug for mesh with higher top dim.

[Prof Garth Wells]

And how to you define the neighbour of a vertex?

It is perfectly possible to get vertices that are 'connected' via an entity of a selected dimension via the present interface.

[Johan Hake]

That is the precise functionality that is broken for 0-0 connectivities.

from dolfin import * mesh = UnitSquareMesh(16,16) mesh.init() for v1 in vertices(mesh): print v1.index(), "is connected with", v1.entities(0)

[Anders Logg (Chalmers)]

On Tue, Apr 30, 2013 at 09:34:34AM -0000, Garth Wells wrote:

vertex-vertex relations can be arbitrarily defined in terms of other entities. I'm happy with the current definition where one needs to explicitly work via the connecting entities.

-- Anders

[Prof Garth Wells]

Sample code from Johan does as I would expect.

d-d connectivities are computed on the basis of entities of dimension d that share an entity of dimension 0. 1.0 had an inconsistency that for 0-0, it considered neighbours to share a cell. Removing this inconsistency fixed some bugs for 1D problems.

I think that a more consistent d-d behaviour would be for the entity to just return itself. E.g., the the below code you report only '1' connection for each edge:

from dolfin import * 
mesh = UnitSquareMesh(16, 16) 
mesh.init() 
for v1 in edges(mesh): 
   print v1.index(), "is connected with", len(v1.entities(1))

I don't expect that this would break any library code.

[Johan Hake]

What is the point of just returning itself for d-d connectivites? There are plenty of times one would like to know what other entities of the same topological dimension that is connected to itself. This is in my mind broken for 0-0 connectivities.

Also isn't the inconsistency experiences in 1.0 only valid for 1D meshes where the facet topological dimension is the same as the vertex topological dimension? One could easily fix this by adding an extra check for the topological dimension of the mesh.

[Prof Garth Wells]

You still have not answered "what constitutes a connection?". A shared vertex, shared edge, shared facet. . . ? We can't discuss this properly unless you assert what you think constitutes a connection.

[Johan Hake]

I see what you mean. I bet if we implemented the "sieve" this would just work out. Not having that luxury available i would go for an intuitive and useful approach:

• 0-0 are connected if they share a 1 entity (potential exception for topdim = 1)
• 1-1 are connected if they share a 0 entity
• 2-2 are connected if they share a 1 entity
• 3-3 are connected if they share a 2 entity

[Prof Garth Wells]

I would say that it 'works out' now. The assumption for d-d is that connections go via a common vertex.

Your proposal is for connections d-d via a common facet rather than a common vertex. This is in fact what the 0-0 case does in the present implementation - it returns vertices that 'share' a facet (which is the vertex itself).

I don't mind how we have it now, but do think that d-d connections being a connection to the entity itself would be unambiguous across all dimensions. It would be up to the users to decide what constitutes a connection between d - d entities.

[Johan Hake]

Well then I guess we disagree. I would rather have the connectivity decided by a common facet with the exception for d=0, where it is the edge (not the facet) that determines the connectivity.

• It creates smaller data structures.
• I would also argue that it is easier to get the vertex based d-d connectivities from a facet based connectivity.

[Prof Garth Wells]

So you're arguing for dimension-dependent code, and I'm arguing for dimension-independent code? Striving for dimension-independent code is a FEniCS virtue.

Getting entities connected by a facet is presently very easy, and the presently implemented d-d connectivity makes this possible with using common code for any dimension. It's a key part of the parallel code. It would be crazy to unnecessary write special cases.

[Johan Hake]

The only dimension dependent part of my argument is for 0-0 connectivity which is calculated based on sharing a common edge for all topological mesh dimensions. Would that be crazy?

[Prof Garth Wells]

Crazy - yes. You would need to program special cases to run in parallel in 1D.

[Johan Hake]

That do sounds crazy. However I am not sure my proposal suggest that.

I propose to use common d - 1 entities to connect all d - d entities, also for d < mesh.topology().dim(). However, for d=0 we use a common 1 entity. I would argue that this is consistent for meshes of all topological dimensions. I have no clue how general it would be to implement but I guess it would not be any worse than we have now.

[Prof Garth Wells]

Using d-1 is bad. It makes a special case for d = 0.

Review the code for building a BoundaryMesh and the changes that were made more all dimensions to run in parallel without special cases.

This report should be closed because it's not a bug. It is a discussion on the design of the Mesh library.

[Johan Hake]

Ok I will have a look and I rest my case for the moment ;)

FYI, we already have a special case for the 0-0 connectivity. See TopologyComputation.cpp line 237.

[Anders Logg (Chalmers)]

Just some notes on this issue:

• An entity is always defined in terms of its vertices
• Two entities are adjacent if they share a vertex

It follows that a vertex is only connected to itself

Because it is sometimes useful to loop over vertices connected to a vertex (a connectivity which then needs to be defined), we previously defined this as sharing a common cell. This is reasonable since it covers the maximum case (it includes connections through edges and faces).

[Johan Hake]

Ok, it make sense. However, my intuition started out somewhere else and it is still, apparently, caught in a local minimum.

Johan

[Fotini Karakatsani]

Hello,

if I understood well the example in fenics book which gives the neighbors of a vertex is not valid anymore. I would like to compute the local extremum of a finite element function. Can you please tell me how I implement that in new version of fenics?

thank you in advance Fotini

[Johan Hake]

Yes, the text about MeshIterators and how vertices are connected are not longer valid.

To collect a neighborhood of vertices you can do something like:

for v1 in vertices(mesh):
    print "vertex", v1.index()
    neighbors = set();
    for e in edges(v1):
        neighbors.update(e.entities(0))
    neighbors.remove(v1.index())
    print "connections", neighbors

[Fotini Karakatsani]

OK Johan, thanks a lot! I just want to make a general comment as a user. It would be nice if the different version of FEnics are compatible. I had used another FE library and I stopped using it because of the incompatibility of the versions. It was better to write the code from the beginning than to adapt your code to the new version. I hope that there aren't many changes between the different versions of FEnics.

Thank you again Fotini

[Prof Garth Wells]

For reference, here's the bug report (from Launchpad, https://bugs.launchpad.net/dolfin/+bug/1063013) that showed up the inconsistency in older versions of DOLFIN.

BoundaryMesh gets it wrong in 1D (which is breaking 1D in parallel). Instead of picking up the end vertices, it picks up vertices one position away from the end. The code

    from dolfin import *
    mesh = UnitInterval(3)
    bmesh = BoundaryMesh(mesh)
    info(bmesh, True)

prints

    0: 0.333333
    1: 0.666667

when it should print

    0: 0.0
    1: 1.0