Unexpected behaviour, if evaluating function on grid

[LukasDomenig]

Hi,

I am encountering an unexpected behaviour, if I execute the following code:

using  GridapGmsh
using  Gridap
using  Gridap.Geometry
using  LineSearches: BackTracking
using  LinearAlgebra: norm
import Random
using Plots
using DataFrames
using CSV

model = GmshDiscreteModel("../../mesh/for_gridap2d/mlm2d.msh")
writevtk(model ,"mesh_info/mlm2d")

Ω = Triangulation(model)

function Hs(x)
    J = 1e6
    H_sx = 0.0
    if 0.025 <= x[1] && x[1] <= 0.0375 && 0.075 <= x[2] && x[2] <= 0.125
        H_sy = J*(x[1]-0.025)
    elseif 0.0375 < x[1] && x[1] <= 0.0625 && 0.075 <= x[2] && x[2] <= 0.125
        H_sy = J*0.0125 
    elseif 0.0625 < x[1] && x[1] <= 0.075 && 0.075 <= x[2] && x[2] <= 0.125
        H_sy = J*(0.075-x[1])
    else
        H_sy = 0.0
    end

    return VectorValue(H_sx,H_sy)
end

writevtk(Ω,"Hs_field",cellfields=["Hs"=>Hs])

This is the field magnitude I am obtaining over the grid, if this code is executed: Zoomed in version:

On the top and bottom edge the function should go discontinuously to zero, however somehow this is not happening, but there is I think a continuous linearly decrease. Am I doing something wrong, or does gridap behave strange in this case?

I tried to project the Hs field into several function spaces (H1,Hcurl,L2), however this did not help. I have also tried to interpolate the Hs field similar to the tutorial 15 (section: Interpolating vector-valued functions), but this also did not fix this wrong behaviour.

I also attached the .msh file as .txt (The path to generate the model has to be changed, to try the code) mlm2d.txt

I would appreciate any insights or assistance to resolve this issue.

[ericneiva]

Hi, @LukasDomenig, this is the expected result of running that code.

ParaView linearly interpolates the grid nodal values to draw the values on the whole domain. If you plot the nodes or elements on top of the plot you show, you'll immediately see this.

If you want to have a sharp discontinuity on those edges, I suggest you to work with L2 function spaces and modify Hs such that it discriminates whether the node you are evaluating belongs to an element inside the region of interest or not (maybe, e.g., tagging the region in Gmsh).

[LukasDomenig]

I'll try it out. Thank you for the quick reply!