base points associated to dof

[edljk]

Thank you for developing Ferrite.jl!
I am using Lagrange finite elements.
In that case, is it possible for every dof to get the coordinates of a point representing the associated basis element. I am thinking to a function analogous to the following in getfem++

DOFpts = gf_mesh_fem_get(mesh_fem MF, 'basic dof nodes'[, mat DOFids])
Get location of basic degrees of freedom.

Thanks

[fredrikekre]

If I understand you correctly you want to find the coordinate x for a given input dof y? I don't think there is an easy way to do that right now, but probably a hacky way. Can I ask what you need it for?

[edljk]

You are welcome! I am working with polygons / polytops with strange periodicity conditions. To enforce the global regularity of a fem function on such domains, I need to be able to merge dofs which are equivalent up to these periodic conditions. I understand that this is a non standard use of fem..
Is the motivation clear enough?

[edljk]

Does the following hack (inspired by test functions) looks reasonable to you?

function dof_points(grid, qr, ip, ip_geom)
    dim = length(grid.nodes[1].x)
    fx = x -> x[1]
    fy = x -> x[2]
    valsx = _analytical(grid, qr, ip, ip_geom, fx)
    valsy = _analytical(grid, qr, ip, ip_geom, fy)
    L2p = L2Projector(ip, grid, geom_ip = ip_geom)
    xdof = project(L2p, valsx, qr, project_to_nodes = false)
    ydof = project(L2p, valsy, qr, project_to_nodes = false)
    if dim == 2
        return hcat(xdof, ydof)
    else
        fz = x -> x[3]
        valsz = _analytical(grid, qr, ip, ip_geom, fz)
        zdof = project(L2p, valsz, qr, project_to_nodes = false)
        return hcat(xdof, ydof, zdof)
    end
    nothing
end
function _analytical(grid, qr, ip, ip_geom, f)
    qp_values = []
    cv = CellScalarValues(qr, ip, ip_geom)
    for cell in CellIterator(grid)
        reinit!(cv, cell)
        r = [f(spatial_coordinate(cv, qp, getcoordinates(cell))) for qp in 
             1:getnquadpoints(cv)]
        push!(qp_values, r)
    end
    return identity.(qp_values)
end

[edljk]

Previous approach seems to be ok using P2 elements on linear (standard) triangular mesh

But for curvelinear mesh of order 2, even if it seems almost correct,

there is a small difference of coordinates between the initial vertices of triangles and the corresponding dof points. This should come from the L2 projection. Something should not be correct in my implementation, I expected a perfect match at these points. See the zoom below

[edljk]

I identified my mistake. I was using an inexact quadrature order: Switching from 3 to 4 for fem of degree 2 gives the good dof points also with curved mesh (of order 2). Thanks.

[fredrikekre]

There is probably a more efficient way of doing it, but great that you got it working!

Note that there is some support for periodic boundary conditions, where you don't have to do this manipulation manually like this: https://ferrite-fem.github.io/Ferrite.jl/stable/manual/boundary_conditions/#Periodic-boundary-conditions but currently it only works for hypercubes. I am currently generalizing though, so should be much simpler in the near future!