Incorrect facet normal with variable layer extrusion

[stephankramer]

Following on from #1835, I get nans in any surface integral that includes the facet normal over the left boundary (this is an extruded upwards mesh with an increasing number of cells going left to right: 1 on the left and 21 on the right) :

from firedrake import *
from math import ceil

L = 100
H1 = 2.
H2 = 42.
dy = 5.0
ny = round(L/dy)
dz = 2.0

# create mesh
mesh1d = IntervalMesh(ny, L)
layers = []
cell = 0
yr = 0
for i in range(ny):
    yr += dy  # y of rhs of column (assumed to be the higher one)
    height = H1 + yr/L * (H2-H1)
    ncells = ceil(height/dz)
    layers.append([0, ncells])
    cell += ncells

mesh = ExtrudedMesh(mesh1d, layers, layer_height=dz)
y = mesh.coordinates.dat.data_ro[:,0]
z = mesh.coordinates.dat.data_ro[:,1]
mesh.coordinates.dat.data[:,1] = np.minimum(H1 + y/L * (H2-H1), z)
n = FacetNormal(mesh)
print(assemble(Constant(1.0)*ds_v(1, domain=mesh)))
print(assemble(dot(n, as_vector((1,0)))*ds_v(1, domain=mesh)))
print(assemble(n[0]*ds_v(1, domain=mesh)))
print(assemble(n[1]*ds_v(1, domain=mesh)))
y, z = SpatialCoordinate(mesh)
print(assemble(conditional(z<2.0, n[0], 0)*ds_v(1, domain=mesh)))

produces:

2.0
nan
nan
nan
-2.0000000000000004

So I imagine this has something to do with the squashed facets on the left boundary?

[wence-]

I can't immediately see what the setup is. Do you have some geometrically measure-zero facets?

[stephankramer]

Yes, but that's by design, and will be hard to avoid (if I understand your question right): the left boundary has two facets, one of which has zero length:

Presumably a similar issue would occur with dS_v as well - assuming it's to do with measure-zero facets? - but that seems to behave:

print(assemble(avg(dot(n,n))*dS_v(domain=mesh)))

gives a believable number.

[wence-]

I think that my original implementation claims that it integrates over all topologically "exposed" exterior faces (so none of them are allowed to be zero-sized).

For in the interior faces, we integrate over the intersection of the exposed faces from both sides (so when there's a step, we don't integrate over that face)

[stephankramer]

Ah ok, I see. So I suppose as a workaround I should just fiddle with the top boundary mesh nodes to ensure some minimal vertical edge length. A little awkward in 3d and steep bathymetries, but doable.

Any thoughts on whether this could be improved? Would it be possible to just leave these facets out of the set of exterior facets? I guess the underlying issue is that we don't specify what we want when the mesh is created.

[wence-]

I think it's difficult to cull them a priori because I set this up when I only have topology, not geometry. We could add an option for you to say "these bits should be squashed out", but it would probably be a little fiddly. Is there a reason that in the real cases you would end up in a situation like this?

[stephankramer]

In what I would describe as a real case with z-layers, you would actually have a bathymetry (or height) that's located at the vertices of the horizontal mesh, i.e. it's described as a continuous piecewise linear function (of the horizontal coordinates). Whenever that bathymetry intersects with any of the fixed z-layers in a (horizontal) cell adjacent to the boundary this issue would occur; in particular it's guaranteed to happen when the bathymetry is steep dx>dz, where you'd typically end up with more than one zero-measure facet. What I was planning to do was to write a kernel that takes that bathymetry (and/or top height), establishes the required number of layers on top of each horizontal vertex, and then take the maximum of that in each horizontal cell over its adjacent vertices to establish the number of required vertical cells in that column. In the interior I can then just move the "extra" nodes, i.e. the nodes that are outside the range of the specified bathymetry and top height, to within the range.

Oh hm writing that out I thought of another case where this goes wrong in 3D: if two adjacent horizontal triangles ask for more layers in the two opposite vertices they don't share than in the two vertices they do share, than above the shared edge it will have zero-area facets that it still wants to integrate over (because the corresponding layer exists in both adjacent columns). Ugh...

[wence-]

I think a cell column has to have a number of layers that is the min rather than the max of the number of layers of the vertices it touches, that might fix things?

I guess that would mean you wouldn't be able to represent

         +  
        /|\ 
       / | \
+--+--+--+--+--+--+
|  |  |  |  |  |  |
|  |  |  |  |  |  |
+--+--+--+--+--+--+
|  |  |  |  |  |  |
|  |  |  |  |  |  |
+--+--+--+--+--+--+
|  |  |  |  |  |  |
|  |  |  |  |  |  |
+--+--+--+--+--+--+

?

[wence-]

Oh no, that would be OK I think.

[stephankramer]

Do you mean

         +  
        /|\ 
       / | \
+--+--+  |  +--+--+
|  |  |  |  |  |  |
|  |  |  |  |  |  |
+--+--+--+--+--+--+
|  |  |  |  |  |  |
|  |  |  |  |  |  |
+--+--+--+--+--+--+
|  |  |  |  |  |  |
|  |  |  |  |  |  |
+--+--+--+--+--+--+

I think that only works because all nodal columns have the same number of layers What if I want:

+  
|\ 
| \
+--+
|  |\
|  | \
+--+--+
|  |  |\
|  |  | \
+--+--+--+
|  |  |  |\
|  |  |  | \
+--+--+--+--+
|  |  |  |  |

By taking the minimum, as you suggest, I would get (before moving any nodes):

+--+
|  |
|  |
+--+--+
|  |  |
|  |  |
+--+--+--+
|  |  |  |
|  |  |  |
+--+--+--+--+
|  |  |  |  |

To make that follow the desired geometry, in a continuous sense, I can move up the top left noded, but in the rest of the domain I have to move up interior nodes:

+ 
|\
| \
|  +
| /|\
|/ | \
+  |  +
|  | /|\
|  |/ |  \
+--+  |  +
|  |  | /|\
|  |  |/ |  \
+--+--+  |  +
|  |  |  |  |

thus creating zero-measure facets again. And I think you would still have a similar problem in 3D: in any nodal column you need to move up any interior node that is initially at a level equal to the minimum top level in any adjacent nodal column, up to the position of the top exterior node in its column, So if you have two triangles that agree on the number of levels in their prismatic column, the shared vertices (nodal columns) might still be connected to nodal columns with fewer nodes/levels, and thus the corresponding node in the shared nodal columns need moving up, creating facets that face both prismatic columns and have zero area.

[wence-]

Hmm. The interior measure zero facets are fine I think. It's just if you want your domain to have a wedge-shaped (in the vertical direction) element at a domain boundary.

But I'm bad at picturing all these cases.

We could add some facility to turn facets off and on in a data-dependent fashion: the topological computations remain the same, we just need some extra data telling us where to compute

[stephankramer]

Yeah, it does my head in as well a little.

I'm not sure I understand why you say interior measure zero facets are fine. Previously you said they are fine, because "we integrate over the intersection of the exposed faces from both sides" - but as I argued above that logic does not work in 3D

And indeed:

from firedrake import *
from math import ceil

# create mesh
mesh2d = UnitSquareMesh(1,1)
layers = [[0,2], [0,2]]
mesh = ExtrudedMesh(mesh2d, layers, layer_height=1.0)
x = mesh.coordinates.dat.data_ro[:,0]
y = mesh.coordinates.dat.data_ro[:,1]
z = mesh.coordinates.dat.data_ro[:,2]
mesh.coordinates.dat.data[:,2] = np.minimum(1 + abs(x+y-1), z)

n = FacetNormal(mesh)
print(assemble(avg(dot(n,n))*dS_v(domain=mesh)))
File('u.pvd').write(Function(FunctionSpace(mesh, "CG", 1)))

gives me a nan.

If you say I should set the number of cells in each cell column to be the minimum over the number of layers required in each adjacent nodal column (as opposed to maximum that was used in this figure) - then in the above I would loose the the two heighest vertices, and have to move up the vertices below, from z=1 to z=2, to get the same geometry again. That works if there are just 2 triangles, but if one of those two higher nodal columns of height 2, say nodal column i, is connected to a triangle/prismatic column of which the other two nodal columns also have height 2, then the number of layers in that prismatic column would have to be two (the minimum over the original number of layers/edges (n/o nodes+1) in each adjacent nodal column) which would create a zero length vertical edge at the top of nodal column i. If this is true for two adjacent triangles/prismatic columns c and d, both adjacent to nodal column i and a second nodal column j - then it could easily be true that the same logic applies to nodal column j, and it also would have a zero-length edge on top, in which case we end up again with an interior facet of zero area that faces both column c and column d.

If you say, well let's reduce the number of layers in c and d after we have established that we want fewer layers in column i, then you just end up reducing the number of layers in more and more columns.

[wence-]

Ah crap. OK, it seems like we can't get away with setting up the topology right. So I suppose one needs some data-dependent way of masking these things.

[stephankramer]

In terms of a user interface it would be nice to have a firedrake function that generates the correct extruded mesh from a provided P1 bathymetry (and/or top height). I think that would be generic functionality that would be useful to various users. Internally, I don't how hard it is to add this masking of zero-measure facets - is it just a matter of restricting the subsets associated with ds_v and dS_v - does it require pyop2 changes? For now I'm happy to just go with the workaround of having very short rather than zero vertical length edges. Shall I just mark this as a wishlist item then?

[wence-]

there are minor pyop2 changes (he says hopefully) required, I think.