cell diameter computed incorrectly; hmax gets bigger with refinement

[Douglas Arnold]

The following code uses mshr to create a mesh, creates a second mesh by refining the first, and reports hmax for both of them. Note that it reports that hmax gets larger with refinement, which, unless I am badly understanding, is certainly not true.

from dolfin import *
from mshr import *
cylinder = Cylinder(Point(0., 0., -1.), Point(0., 0., 1.), 1., 1.)
mesh1 = generate_mesh(cylinder, 1)
mesh2 =  refine(mesh1)
print("mesh1.hmax() = {}  mesh2.hmax() = {}".format(mesh1.hmax(), mesh2.hmax()))

The output under Dolfin version 1.6.0 is

mesh1.hmax() = 1.99438988373  mesh2.hmax() = 3.49259294361

I did some experimenting, and discovered that if c is a cell of the mesh, then c.diameter(), which should be the diameter of the cell, is calculated incorrectly. In the following code I calculate the diameter of each element directly and compute my own hmax. This decreases as it should.

from dolfin import *
from mshr import *
import numpy as np

def diam(c):
    # reimplementation of c.diameter(): diameter of a cell c
    verts = c.mesh().coordinates()[c.entities(0)] # list of vertex coordinates
    d = 0.
    for i in range(4):
        for j in range(i):
            d = max(d, np.linalg.norm(verts[i] - verts[j]))
    return d

def hmx(mesh):
    # reimplementation of mesh.hmax(): largest cell diameter in mesh
    h = [diam(c) for c in cells(mesh)]
    return max(h)

cylinder = Cylinder(Point(0., 0., -1.), Point(0., 0., 1.), 1., 1.)
mesh1 = generate_mesh(cylinder, 1)
print("hmax = {}  hmx = {}".format(mesh1.hmax(), hmx(mesh1)))
mesh2 =  refine(mesh1)
print("hmax = {}  hmx = {}".format(mesh2.hmax(), hmx(mesh2)))

The output is:

hmax = 1.99438988373  hmx = 1.99259014276
hmax = 3.49259294361  hmx = 1.22433824285

[Douglas Arnold]

On further study I see that I am running into this previously reported bug:

https://bugs.launchpad.net/dolfin/+bug/1132490

The diameter() method is returning twice the circumradius of the simplex, not what is commonly called the diameter. E.g., according to Wikipedia, "the diameter of a subset of a metric space is the least upper bound of the set of all distances between pairs of points in the subset," which is the standard usage in math and finite elements. For a simplex, this diameter is the maximum edge length and it can easily happen that it tends to zero for a sequence of simplices, while the circumradius tends to infinity. Similarly hmax() is returning the maximum circumradius of a mesh, not what is commonly denoted hmax in finite element theory, which is the maximum diameter in the sense of maximum edge length.

I think this usage should be changed, or at least carefully documented in the doc strings. Is there a Dolfin method to compute the diameter and hmax in their usual senses?

[Prof Garth Wells]

A 'workaround' is to use Cell::inradius. With

from dolfin import *
from mshr import *
cylinder = Cylinder(Point(0., 0., -1.), Point(0., 0., 1.), 1., 1.)
mesh1 = generate_mesh(cylinder, 1)
mesh2 =  refine(mesh1)
print("mesh1.hmax() = {}  mesh2.hmax() = {}".format(mesh1.hmax(), mesh2.hmax()))

# Test inradius
h1 = []
for cell in cells(mesh1):
    h1.append(cell.inradius())
h2 = []
for cell in cells(mesh2):
    h2.append(cell.inradius())
print("mesh1 hmax* = {}  mesh2 hmax* = {}".format(max(h1), max(h2)))

I get

mesh1.hmax() = 1.99438988373  mesh2.hmax() = 3.49259294361
mesh1 hmax* = 0.261217910523  mesh2 hmax* = 0.133399237347

We need to fix this. Changing what diameter returns could silently change the result of some user code. We could change it's behaviour, or we could remove diameter (I don't like papering over something misleading with a docstring). My preference would be to change the behaviour of diameter and add a new function Cell::circumradius.

[Prof Garth Wells]

Also, Mesh::hmax() should probably take an enum to pick which cell size measure is used.

[Douglas Arnold]

Neither circumradius or inradius is an acceptable substitute for cell diameter. For example, if you consider an isoceles triangle with base 1 and height tending to zero, the inradius tends to zero, the diameter remains 1, and the circumradius tends to infinity. The diameter of a simplex is almost universally used as the measure of its size in finite elements (typically denoted by h), so the fact that it is missing from FEniCS is unfortunate. Inradius and, to a lesser extent, circumradius, are used in shape measures. The current Cell::diameter should really be called Cell::circumdiameter or, better yet, replaced by Cell::circumradius, which would be 1/2 the current value.

I am sympathetic with the issue of backward compatibility, but it might be worth living with the incompatibility in this case in order to restore compatibility with common usage. Alternatively, one could introduce Cell::h for the true cell diameter, and retire Cell::diameter with a warning that its old meaning is now available as twice Cell::circumradius while the true cell diameter is now available as Cell::h. The name Cell::h would be consistent with Mesh::hmax which would be redefined as the maximum of the Cell::h values over all cells in the mesh.

[Prof Garth Wells]

I'd like to hear opinions on two options:

1. Deprecate Cell::diameter, and add Cell::h and Cell::circumradius
2. Change the behaviour of Cell::diameter, and add Cell::circumradius

My preference is for 2.

[Anders Logg (Chalmers)]

My preference is for 1.

[Douglas Arnold]

I'll vote for a variant of 1:

1A. Add Cell::h and Cell::circumradius, leaving Cell::diameter as is for backward compatibility. Mesh::hmax would change to the maximum of the Cell::h.

However, I don't feel strongly as I don't have any code that would be affected by the disappearance or changing of Cell::diameter.

[Martin Sandve Alnæs]

You can express the maximum (and minimum) edge length of each cell in variational forms using h = MaxCellEdgeLength(mesh) and some other similar types:

    CellVolume, Circumradius, FacetArea, 
    MinCellEdgeLength, MaxCellEdgeLength,
    MinFacetEdgeLength, MaxFacetEdgeLength,

Here MaxFacetEdgeLength is only available in facet integrals, and MaxCellEdgeLength must be restricted to + or - side in an interior facet integral.

You can use these in assemble and project, and to avoid the expensive projection you can assemble the cell versions to a DG0 space (untested code draft):

from fenics import *
parameters["form_compiler"]["representation"] = "uflacs"
mesh = UnitSquareMesh(3, 8)
hmax = MaxCellEdgeLength(mesh)
hmin = MinCellEdgeLength(mesh)
ratio = hmax / hmin

V0 = FunctionSpace(mesh, "DG", 0)
v = TestFunction(V0)
H = CellVolume(mesh)
r = Function(V0)
assemble(ratio/H*dx, tensor=r.vector())
plot(r)

Not sure if the uflacs representation is necessary.

[Martin Sandve Alnæs]

Note that Circumradius is also in ufl. Its definition shouldn't be changed because that will silently change the numerical results of users codes. If necessary it's better to introduce replacements with other names and deprecate it.

[Prof Garth Wells]

Add function Cell::h(), which computes greatest distance between any two vertices of a cell. For simplices, this is the greatest edge length.

• Deprecate Cell::diameter()
• Add Cell::circumradius()
• Change Mesh::hmin/max to use Cell::h()

Fixes Issue #664.

→ <<cset a565417d05dd>>

[Prof Garth Wells]

I implemented Cell::h() to be the greatest distance between any two vertices such that it will work for quads and hexes. Is this acceptable for quads and hexes (we don't support quads and hexes in the form compiler yet, but DOLFIN does now support quad and hex meshes)?

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)

[Jan Blechta]

Is this acceptable for quads and hexes

Yes. That's exact for P1 or Q1 cells. For higher-order cells it can give arbitrarily poor approximation.