The proper way of setting concentrated loads on the right-hand side vector (node to DoF index)

[suurj]

This might be a simple issue, but is there a handy way to assign arbitrary values on the right-hand side vector that correspond to the degrees of freedom of a specific node and field? That is, I want a function like map_node_to_dof(nodeset["s"], :field) => DoFindices

Let us say that I want to set a concentrated load in a linear elasticity problem on one node in the mesh. In that case, I would not use the typical Neumann boundary condition assembly procedure (i.e. a scenario of distributed load) which iterates through the Neumann faces of the mesh and computes the surface integrals along them before assembling the surface integrals onto the load vector. Instead, I would need the DoF indices of the variable of interest at the given node(s) and modify directly the right-hand side vector.

[lijas]

There has never been a straight forward way of extracting the dofs on nodes in Ferrite. The reason for this is because we distribute the dofs on cells, which means we never store the node-to-dof information. The reason we dont have utility functions for this (yet), is because the weak form is always expressed in terms of surface integrals (and not points). So typically the external forces are applied as neumann integrals.

That being said, it is of course a common thing to want to apply forces to nodes, so we have talked about adding utility functions for this 🙂 . But in the meanwhile, maybe this example code can help you?

using Ferrite

function vertexdofs(dh::DofHandler, vertexid::VertexIndex)

    cellid, lvidx = vertexid
    nfields = length(dh.field_names)

    local_vertex_dofs = Int[]
    for ifield in 1:nfields
        field_name = dh.field_names[ifield]
        offset = Ferrite.field_offset(dh, field_name)
        field_dim = Ferrite.getfielddim(dh, field_name)
        field_ip = dh.field_interpolations[ifield]
        vert = Ferrite.vertexdof_indices(field_ip)[lvidx]
        
        for vdof in vert, d in 1:field_dim
            push!(local_vertex_dofs, (vdof-1)*field_dim + d + offset)
        end
    end

    dofs = zeros(Int, ndofs_per_cell(dh, cellid))
    celldofs!(dofs, dh, cellid)
    return dofs[local_vertex_dofs]
end

function nodeid_to_vertexindex(grid::Grid, nodeid::Int)
    for (cellid, cell) in enumerate(grid.cells)
        for (i, nodeid2) in enumerate(cell.nodes)
            if nodeid == nodeid2
                return VertexIndex(cellid,i)
            end    
        end
    end
    error("Node $(nodeid) does not belong to any cell")
end

grid = generate_grid(Quadrilateral, (10,10))
dh = DofHandler(grid)
add!(dh, :u, 2, Lagrange{2,RefCube,2}()) #displacement dofs
add!(dh, :p, 1, Lagrange{2,RefCube,1}()) #pressure dofs
close!(dh)

#Node id to add force on
nodeid = 11

#convert nodeid to VertexIndex (cellid, local-node-id) since we need the 
# the cell id to extract the dofs
vertexid = nodeid_to_vertexindex(grid, nodeid)

dof = vertexdofs(dh, vertexid)
udofs = dof[1:3]
#f[udof] = [0.0, 0.0, 1.0] #add force

[suurj]

Thanks for the reply! There might be a version issue since I get this error while trying to run the code:
ERROR: MethodError: no method matching field_offset(::DofHandler{2, Grid{2, Quadrilateral, Float64}}, ::Symbol)

I am using the master branch version of Ferrite.

[lijas]

This was on version 0.3.14. I will see if I can modify the script tomorrow.

[lijas]

Hello again. I think this should work on master

using Ferrite

function vertexdofs(dh::DofHandler, vertexid::VertexIndex)

    cellid, lvidx = vertexid
    nfields = length(dh.field_names)
    sdh = dh.subdofhandlers[dh.cell_to_subdofhandler[cellid]]
    local_vertex_dofs = Int[]
    for ifield in 1:length(sdh.field_names)
        offset = Ferrite.field_offset(sdh, ifield)
        field_dim = Ferrite.getfielddim(sdh, ifield)
        _field_ip = sdh.field_interpolations[ifield]
        if _field_ip isa Ferrite.VectorizedInterpolation
            field_ip = _field_ip.ip
        else
            field_ip = _field_ip
        end
        vert = Ferrite.vertexdof_indices(field_ip)[lvidx]
        
        for vdof in vert, d in 1:field_dim 
            push!(local_vertex_dofs, (vdof-1)*field_dim + d + offset)
        end
    end

    dofs = zeros(Int, ndofs_per_cell(dh, cellid))
    celldofs!(dofs, dh, cellid)

    return dofs[local_vertex_dofs]
end

function nodeid_to_vertexindex(grid::Grid, nodeid::Int)
    for (cellid, cell) in enumerate(grid.cells)
        for (i, nodeid2) in enumerate(cell.nodes)
            if nodeid == nodeid2
                return VertexIndex(cellid,i)
            end    
        end
    end
    error("Node $(nodeid) does not belong to any cell")
end

grid = generate_grid(Quadrilateral, (10,10))
dh = DofHandler(grid)
add!(dh, :u, Lagrange{2,RefCube,2}()^2) #displacement dofs
add!(dh, :p, Lagrange{2,RefCube,1}()) #pressure dofs
close!(dh)

#Node id to add force on
nodeid = 11

#convert nodeid to VertexIndex (cellid, local-node-id) since we need the 
# the cell id to extract the dofs
vertexid = nodeid_to_vertexindex(grid, nodeid)

dof = vertexdofs(dh, vertexid)
udofs = dof[1:3]
#f[udof] = [0.0, 0.0, 1.0] #add force

[suurj]

Many thanks, it works! As you pointed out, something like this should be really incorporated into the official codebase.