Implementation of Nedelec interpolation on tetrahedral and hexahedral elements

[gijswl]

First attempt at implementing a VectorInterpolation on a three-dimensional reference element, based on the discussion in #1159 .

@KnutAM I am still unsure of get_direction. Right now it passes the continuity test, but I am not completely convinced that this is right. Perhaps I'm misunderstanding what that function is supposed to return.

[gijswl]

What is this expected to return?

[KnutAM]

It should return +1 if the edge associated with shape function j is positive, and -1 if negative. A positive edge is defined by the global node number of the first vertex being smaller than that of the second.

[gijswl]

What about shape functions associated with the cell faces (such as would be necessary for a Raviart-Thomas interpolation)?

[KnutAM]

I don't think we have defined what a positive face is in Ferrite, but maybe the definition that @fredrikekre suggested is already in place in some parts:
A face is positive if, starting from the lowest node number, the next vertex' node number is smaller than the previous vertex number.

Here, the vertices are defined by Ferrite.faces(cell), e.g. for cell = Tetrahedron((3,5,4,7)), we have
Ferrite.faces(Tetrahedron((3,5,4,2)) = ((3, 4, 5), (3, 5, 2), (5, 4, 2), (3, 2, 4))

Edit (was tired):

• face 1: (3, 4, 5) - positive: 5-3-4 (5 > 4)
• face 2: (3, 5, 2) - positive: 5-2-3 (5 > 3)
• face 3: (5, 4, 2) - negative: 4-2-5 (4 < 5)
• face 4: (3, 2, 4) - negative: 3-2-4 (3 < 4)

(Please correct me if I'm wrong with the convention @fredrikekre ! (Edit: Tag correct Fredrik :) )

So far, that was not required to use since we didn't have faces as facets for H(div) interpolations.

[gijswl]

That appears to be right. Implementing this reasonsing (for RaviartThomas{RefTetrahedron,1} and RaviartThomas{RefHexahedron,1}), I get passing continuity tests.
I initially thought we should be able to default all the faces to positive, since the description of Ferrite.faces says the tuples it returns are oriented (normals pointing outward). That also appears to be the approach taken for the 2D elements.

[gijswl]

Is this a valid generalization, or do we need to implement additional tests for interpolation functionals with rdim > 2

[KnutAM]

For $H(\mathrm{curl})$ (edge elements) I think this should be sufficient, but haven't validated - need to check what the actual functionals defined on DefElement says if the same can be used - but I don't see why not.

For $H(\text{div})$, we need two parameters in general and have to integrate a face, so havent' implemented how that should be done.

[gijswl]

I've made an attempt at validating Raviart-Thomas interpolation on tetrahedral and hexahedral elements, by computing the functional through integration over the faces. quadgk is not very well suited for that, so improvement is probably possible, but I couldn't find another method which allows the integration boundaries to depend on an integration variable (for the triangular faces).

[KnutAM]

We could also use FacetValues I guess, so far I decided to avoid since to make the test purer and avoid using these before they have been tested.

Will look at the code later this week and give feedback!

Really nice work so far!

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 74 lines in your changes missing coverage. Please review.

Project coverage is 0.00%. Comparing base (6eead25) to head (8f47022).
Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
src/interpolations.jl	0.00%	74 Missing ⚠️

❗ There is a different number of reports uploaded between BASE (6eead25) and HEAD (8f47022). Click for more details.

HEAD has 5 uploads less than BASE

Flag	BASE (6eead25)	HEAD (8f47022)
	6	1

Additional details and impacted files

@@            Coverage Diff             @@
##           master   #1162       +/-   ##
==========================================
- Coverage   93.74%   0.00%   -93.75%     
==========================================
  Files          39      39               
  Lines        6238    6234        -4     
==========================================
- Hits         5848       0     -5848     
- Misses        390    6234     +5844     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[KnutAM]

It surprised me that it looks like it passed some tests for the wrong definitions here. Would need to investigate a bit more to ensure that our tests cover sufficiently well.

[KnutAM]

It should return +1 if the edge associated with shape function j is positive, and -1 if negative. A positive edge is defined by the global node number of the first vertex being smaller than that of the second.

[gijswl]

Thanks for reviewing. I agree that test coverage should be investigated.
One thing I ran into while trying to define further vector interpolations (attempted Nedelec{RefTetrahedron,2} and RaviartThomas{RefTetrahedron,1}) is that the unit tests where the functionals are verified do not automatically handle the different combinations of dofs associated with edges and faces. For example a 2nd order Nedelec interpolation has 12 dofs associated with the edges (which are handled), but also 8 dofs associated with the faces (for which there is no implemented method to verify the functional).

I would expect there is a way to implement this in a general manner, but I have to think how that might be done.

[termi-official]

8 dofs associated with the faces

I am not sure if dof distribution in 3D is ready for this.

[gijswl]

For now I have no need of higher-order elements, but good to know that there may be a limitation in the dof distribution. Maybe that should also define the scope of this PR:

• 1st order Nedelec interpolation on RefTetrahedron and RefHexahedron.
• 1st order Raviart-Thomas interpolation on RefTetrahedron and RefHexahedron.
• Careful review of the interpolation unit test coverage for vector interpolations in 2 and 3 spatial dimensions.

[KnutAM]

8 dofs associated with the faces

I am not sure if dof distribution in 3D is ready for this.

Yes, Ferrite is currently limited to one interior face dof, since the ordering of interior facedofs is not properly defined. This is ongoing work to define even for standard nodal interpolations.

[gijswl]

This is not very elegant

[KnutAM]

I think that is quite reasonable, but we could possibly abstract the logic to be more generally applicable.

Just realized that we have

Ferrite.jl/src/Grid/grid.jl

Lines 763 to 779 in bd65884

	function OrientationInfo(path::NTuple{2, Int})
	flipped = first(path) < last(path)
	return OrientationInfo(flipped, 0)
	end
	function OrientationInfo(surface::NTuple{N, Int}) where {N}
	min_idx = argmin(surface)
	shift_index = min_idx - 1
	if min_idx == 1
	flipped = surface[2] < surface[end]
	elseif min_idx == length(surface)
	flipped = surface[1] < surface[end - 1]
	else
	flipped = surface[min_idx + 1] < surface[min_idx - 1]
	end
	return OrientationInfo(flipped, shift_index)
	end

But that gives that "flipped" is positive according to the definition as I gave above - @termi-official, @AbdAlazezAhmed, @fredrikekre - is this intentional or did this use a different definition (essentially reversed)?

IMO we shouldn't use this function directly, but perhaps it would make more sense to do something like,

get_face_direction(cell, facenr) = get_face_direction(faces(cell)[facenr])

get_face_direction(facenodes::Tuple)
    min_idx = argmin(facenodes)
    if min_idx == 1
        positive = surface[2] < surface[end]
    elseif min_idx == length(surface)
        positive = surface[1] < surface[end - 1]
    else
        positive = surface[min_idx + 1] < surface[min_idx - 1]
    end
    return positive ? 1 : -1
end

or your code, which is equally good (except that I think argmin is cleaner (and potentially faster?) than findmin).

Then we should call get_face_direction in OrientationInfo.

We could do the same get_edge_direction abstraction as well, even if the logic is easier, maybe it makes the code easier to read to have this logic in the same place.