Nedelec element on boundaries

[HelgeGehring]

I'm trying to use a Nédélec element on a boundary, but the creation of the FESpace fails.
Is this an error in Gridap? Or am I just creating it wrong?

Thanks a lot for your help and the nice library!

Here a minimal example producing the error

using Gridap
using Gridap.Geometry

domain = (0, 2, 0, 2, 0, 2)
partition = (10, 10, 10)
oldmodel = CartesianDiscreteModel(domain, partition)

labels = get_face_labeling(oldmodel)
face_mask = get_face_mask(labels, ["boundary"], 2)
model = DiscreteModelPortion(DiscreteModel(Polytope{2}, oldmodel), face_mask)

reffe = ReferenceFE(lagrangian,Float64,1)
V = FESpace(model,reffe) # Works

reffe = ReferenceFE(nedelec,1)
V = FESpace(model,reffe) # Doesn't work

[HelgeGehring]

ups, there was a small error in the example, I've updated it.

[simbilod]

+1 this would be very useful!

[HelgeGehring]

Here a simplified version which generates the same error:

I generated two meshes using (python code)

import gmsh

for z in [0,1]:
    gmsh.initialize()

    gmsh.model.add("t1")

    lc = 5
    gmsh.model.geo.addPoint(0, 0, z, lc, 1)
    gmsh.model.geo.addPoint(1, 0, z, lc, 2)
    gmsh.model.geo.addPoint(1, 1, z, lc, 3)
    gmsh.model.geo.addPoint(0, 1, z, lc, 4)

    gmsh.model.geo.addLine(1, 2, 1)
    gmsh.model.geo.addLine(3, 2, 2)
    gmsh.model.geo.addLine(3, 4, 3)
    gmsh.model.geo.addLine(4, 1, 4)

    gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)
    gmsh.model.geo.addPlaneSurface([1], 1)

    gmsh.model.geo.synchronize()
    gmsh.model.addPhysicalGroup(1, [1, 2, 4], 5)
    gmsh.model.addPhysicalGroup(2, [1], name = "My surface")

    gmsh.model.mesh.generate(2)
    gmsh.write(f"t1-{z}.msh")

    gmsh.finalize()

Both meshes are the same, just the z-coordinate for all points differs in t1-1.msh to let Gridap handle it with 3D point dims
Here the meshes: t1.zip

The error from the previous example can be reproduced with:

using Gridap
using GridapGmsh

model = GmshDiscreteModel("t1-1.msh")

println("num_cell_dims ", num_cell_dims(model))
println("num_point_dims ", num_point_dims(model))

reffe = ReferenceFE(lagrangian, Float64, 0)
V = FESpace(model, reffe)

reffe = ReferenceFE(nedelec, 1)
V = FESpace(model, reffe)

While the generation of the lagrangian FESpace works just fine, the nedelec FESpace fails with the same error.

Would be great if someone could have a look, I don't really know how to go on debugging :/

Thanks a lot!

[fverdugo]

Hi @HelgeGehring

It seems that Nedelec RefFEs are not working on manifolds.

Fixing this will provably require to allocate VectorValues here

Gridap.jl/src/FESpaces/CurlConformingFESpaces.jl

Line 29 in ffe4667

face_moments_out = basis.face_moments[face]

with length equal to the dimension of the ambient space (instead of the reference space) and use the pseudoinverse here:

Gridap.jl/src/ReferenceFEs/NedelecRefFEs.jl

Line 318 in ffe4667

v⋅transpose(inv(Jt)) # we multiply by the right side to compute the gradient correctly

[HelgeGehring]

Thanks @fverdugo !

I just tried to fix this, but it seems to me, that I cannot just allocate VectorValues of a different dimension in that line.
The basis-struct defines the dimension and I cannot just swap it out :/

Would we need to fix this maybe within the Basis?

[amartinhuertas]

Hi @HelgeGehring !

FYI ... in the following PR:

#577

I introduced a set of misc changes to support RT FEs on general manifolds.

I do not have time now to take a careful look at nedelec elements, but I would guess the changes required would be similar ...

Hope this helps.
Best regards,
Alberto.