Curl error of interpolation is too large for non-conforming Nedelec elements

[wei3li]

Hi @amartinhuertas,

When the mesh is non-conforming and Nedelec elements are used, the error of the FEFunction produced by the interpolate function is too far away from the FEM solution based on the same mesh with the same order. I have the following MWE:

using Gridap, GridapDistributed, GridapP4est
using PartitionedArrays
import MPI
import Printf: @printf

function refine_model(octree_model, eh, refine_frac)
  ∇xeh, dΩ = ∇ × eh, Measure(Triangulation(octree_model), 10)
  err_indicator = map(dc -> sqrt.(get_array(dc)), local_views(∫(eh ⋅ eh + ∇xeh ⋅ ∇xeh)dΩ))
  strategy = FixedFractionAdaptiveFlagsMarkingStrategy(refine_frac, 1e-2)
  ref_coarse_flags = map(partition(get_cell_gids(octree_model.dmodel))) do indices
    flags = zeros(Cint, length(indices))
    flags .= nothing_flag
  end

  update_adaptivity_flags!(
    ref_coarse_flags,
    strategy,
    partition(get_cell_gids(octree_model)),
    err_indicator;
    verbose=false)

  refined_model, _ = GridapP4est.adapt(octree_model, ref_coarse_flags)
  return refined_model
end

function compute_l2_and_hcurl_norms(uh, dΩ)
  ∇xuh = ∇ × uh
  l2normsq, hcurlseminormsq = ∑(∫(uh ⋅ uh)dΩ), ∑(∫(∇xuh ⊙ ∇xuh)dΩ)
  sqrt(l2normsq), sqrt(l2normsq + hcurlseminormsq)
end

function curr_item(items)
  local_views(items).item_ref[]
end

MPI.Init()

u(x) = VectorValue(cospi(x[1])cospi(x[2]), sinpi(x[1])sinpi(x[2]))
f(x) = ∇(∇ ⋅ u)(x) - Δ(u)(x) + u(x)
ranks = distribute_with_mpi(LinearIndices((prod(MPI.Comm_size(MPI.COMM_WORLD)),)))
coarse_model = CartesianDiscreteModel((0, 1, 0, 1), (1, 1))
model = OctreeDistributedDiscreteModel(ranks, coarse_model, 4)
reffe = ReferenceFE(nedelec, Float64, 1)

for i in 0:3
  V = FESpace(model, reffe; conformity=:Hcurl, dirichlet_tags="boundary")
  U = TrialFESpace(V, u)
  Ω = Triangulation(model)
  dΩ, dΩ⁺ = Measure(Ω, 5), Measure(Ω, 10)
  a(u, v) = ∫((∇ × u) ⋅ (∇ × v) + (u ⋅ v))dΩ
  l(v) = ∫(f ⋅ v)dΩ⁺
  uh_fem = solve(AffineFEOperator(a, l, U, V))
  uh_ipl = interpolate(u, U)
  eh_fem, eh_ipl = uh_fem - u, uh_ipl - u
  fem_l2err, fem_hcurlerr = compute_l2_and_hcurl_norms(eh_fem, dΩ⁺)
  ipl_l2err, ipl_hcurlerr = compute_l2_and_hcurl_norms(eh_ipl, dΩ⁺)
  println("\n----------------------- Iteration $(i) -----------------------")
  @printf "FEM           errors: L2 = %.6f, H(curl) = %.6f.\n" fem_l2err fem_hcurlerr
  @printf "Interpolation errors: L2 = %.6f, H(curl) = %.6f.\n" ipl_l2err ipl_hcurlerr
  @printf "Interpolation  / FEM: L2 = %.6f, H(curl) = %.6f.\n" ipl_l2err/fem_l2err ipl_hcurlerr/fem_hcurlerr
  println("-----------------------------------------------------------\n")
  vtk_fields = ["FEM error" => curr_item(eh_fem), "FEM curl error" => curr_item(abs ∘ (∇ × eh_fem)), 
  "Interpolation error" => curr_item(eh_ipl), "Interpolation curl error " => curr_item(abs ∘ (∇ × eh_ipl))]
  writevtk(Triangulation(curr_item(model)), "fem_vs_interpolation", cellfields=vtk_fields)
  model = refine_model(model, uh_fem - u, 0.1)
end

The output of the above script is:

----------------------- Iteration 0 -----------------------
FEM           errors: L2 = 0.001016, H(curl) = 0.006460.
Interpolation errors: L2 = 0.001016, H(curl) = 0.006460.
Interpolation  / FEM: L2 = 1.000139, H(curl) = 1.000023.
-----------------------------------------------------------

----------------------- Iteration 1 -----------------------
FEM           errors: L2 = 0.001006, H(curl) = 0.005423.
Interpolation errors: L2 = 0.001004, H(curl) = 0.010116.
Interpolation  / FEM: L2 = 0.998606, H(curl) = 1.865236.
-----------------------------------------------------------

----------------------- Iteration 2 -----------------------
FEM           errors: L2 = 0.000971, H(curl) = 0.004370.
Interpolation errors: L2 = 0.000967, H(curl) = 0.014372.
Interpolation  / FEM: L2 = 0.995770, H(curl) = 3.288731.
-----------------------------------------------------------

----------------------- Iteration 3 -----------------------
FEM           errors: L2 = 0.000881, H(curl) = 0.003094.
Interpolation errors: L2 = 0.000873, H(curl) = 0.018460.
Interpolation  / FEM: L2 = 0.990350, H(curl) = 5.966791.
-----------------------------------------------------------

At iteration 0, the mesh is conforming, the curl error of the interpolation is nearly the same as that of the FEM solution. However, at iterations 1-3, the mesh becomes non-conforming, and the curl error of the interpolation is increasingly different from that of the FEM solution. I also have the mesh and the curl error comparison at iteration 3.

Notice that the curl error is really large around the jump of different level of refinement elements.

The interpolation has discontinuity in its curl even though my analytic function is very smooth.

[amartinhuertas]

Can you try to interpolate a FE function in the FE space to see if the interpolation error is zero?

[wei3li]

Can you try to interpolate a FE function in the FE space to see if the interpolation error is zero?

Yes, I interpolated the FEM solution into the FE space, and the interpolation error was zero.

[amartinhuertas]

Can you try to do the same test with lagrangian elements?

[amartinhuertas]

Yes, I interpolated the FEM solution into the FE space, and the interpolation error was zero.

Can you try to interpolate an analytical function which you know is going to be in the FE space into the FE space and see whether the interpolation error is zero?

[wei3li]

Can you try to do the same test with lagrangian elements?

Do you mean solving the same Maxwell's equation using Lagrangian elements?

Can you try to interpolate an analytical function which you know is going to be in the FE space into the FE space and see whether the interpolation error is zero?

Yes, the error is zero up to machine precision.

[amartinhuertas]

Do you mean solving the same Maxwell's equation using Lagrangian elements?

No. To solve a poisson problem with lagrangian elements, and see whether you have a similar behaviour as with Nedelec (i.e.. higher interpolation errors at the regions where there are jumps in refinement level).

If the interpolation error is zero with functions in FE space, then I am positive the hanging DoF constraints coefficients are computed correctly, so I am not sure if there is anything wrong in the results you are showing.

In these experiments, you are only using a single quadtree in the forest, correct?

[wei3li]

Here are the results of a Poisson problem test.

----------------------- Iteration 0 -----------------------
FEM           errors: L2 = 0.000246, H1 = 0.012765.
Interpolation errors: L2 = 0.000246, H1 = 0.012766.
Interpolation  / FEM: L2 = 1.000962, H1 = 1.000077.
-----------------------------------------------------------

----------------------- Iteration 1 -----------------------
FEM           errors: L2 = 0.000218, H1 = 0.011541.
Interpolation errors: L2 = 0.000219, H1 = 0.011575.
Interpolation  / FEM: L2 = 1.004337, H1 = 1.002911.
-----------------------------------------------------------

----------------------- Iteration 2 -----------------------
FEM           errors: L2 = 0.000166, H1 = 0.009167.
Interpolation errors: L2 = 0.000168, H1 = 0.009228.
Interpolation  / FEM: L2 = 1.011550, H1 = 1.006596.
-----------------------------------------------------------

----------------------- Iteration 3 -----------------------
FEM           errors: L2 = 0.000111, H1 = 0.006566.
Interpolation errors: L2 = 0.000114, H1 = 0.006671.
Interpolation  / FEM: L2 = 1.029868, H1 = 1.016066.
-----------------------------------------------------------

The Lagrangian interpolation seems correct.

If the interpolation error is zero with functions in FE space, then I am positive the hanging DoF constraints coefficients are computed correctly, so I am not sure if there is anything wrong in the results you are showing.

Is Nedelec interpolation on non-conforming mesh something new? It is for sure an issue that the Nedelec interpolation is not curl-conforming.

In these experiments, you are only using a single quadtree in the forest, correct?

Yes. I also tried 4 quadtrees, the results are similar, i.e. the interpolation error is much larger than the FEM error.

[amartinhuertas]

ok, I will take a look when I have some time.