Development for nodal force in Ferrite

[Aminofa70]

Dear @KnutAM,
@KianAghani, @Kevin-Mattheus-Moerman and I developed the following code for nodal force, for now just a simple case.

The code works well for quad4 element for different boundary conditions, but for tri3 when we have shared nodes with Dirichlet boundary condition only works when
Lx, Ly = 1.0, 1.0
nx, ny = 1, 1

When we don't have share node it works for tri3 as well. Also node_to_dof function works well for assigning dof for nodes and we can add it to ferrite if we could solve the problem for tri3 that I think comes from apply!(K, f_ext, ch),

What do you think? then in two examples we could do a pull request to add the functionality if you like it.
We think, node_to_dof(dh) function does the job.

using Ferrite, SparseArrays

# Material
E, ν = 1.0, 0.3
C_voigt = E * [1.0 ν 0.0; ν 1.0 0.0; 0.0 0.0 (1-ν)/2] / (1-ν^2)
C = fromvoigt(SymmetricTensor{4,2}, C_voigt)

# Mesh
Lx, Ly = 1.0, 1.0
nx, ny = 10, 10
corners = [Vec((0.0,0.0)), Vec((Lx,0.0)), Vec((Lx,Ly)), Vec((0.0,Ly))]
grid = generate_grid(Triangle, (nx,ny), corners)

# Dirichlet BCs
addfacetset!(grid, "support_x", x->x[1] ≈ 0.0)  # fix x
addnodeset!(grid, "support_y", x->x[1] ≈ Lx && x[2] ≈ 0.0) # fix y at bottom-left corner

# # Nodal force
addnodeset!(grid, "nodal_force", x->x[1] ≈ 0.0 && x[2] ≈ Ly)


# DofHandler
dh = DofHandler(grid)
add!(dh, :u, Lagrange{RefTriangle,1}()^2)
close!(dh)

# CellValues
qr = QuadratureRule{RefTriangle}(2)
cv = CellValues(qr, Lagrange{RefTriangle,1}()^2)

# Allocate stiffness
K = allocate_matrix(dh)

# Assemble element stiffness
function assemble_cell!(ke, cellvalues, C)
    for q in 1:getnquadpoints(cellvalues)
        dΩ = getdetJdV(cellvalues, q)
        for i in 1:getnbasefunctions(cellvalues)
            ∇Nᵢ = shape_gradient(cellvalues, q, i)
            for j in 1:getnbasefunctions(cellvalues)
                ∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q, j)
                ke[i,j] += (∇Nᵢ ⊡ C ⊡ ∇ˢʸᵐNⱼ) * dΩ
            end
        end
    end
end

function assemble_global!(K, dh, cellvalues, C)
    ke = zeros(getnbasefunctions(cellvalues), getnbasefunctions(cellvalues))
    assembler = start_assemble(K)
    for cell in CellIterator(dh)
        reinit!(cellvalues, cell)
        fill!(ke, 0.0)
        assemble_cell!(ke, cellvalues, C)
        assemble!(assembler, celldofs(cell), ke)
    end
end

assemble_global!(K, dh, cv, C)

# ConstraintHandler
ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfacetset(grid,"support_x"), (x,t)->0.0, 1))
add!(ch, Dirichlet(:u, getnodeset(grid,"support_y"), (x,t)->0.0, 2))
close!(ch)

# Map node -> DOFs
function node_to_dof(dh)
    map = Dict{Int,Vector{Int}}()
    for cellid in 1:getncells(dh.grid)
        cell = getcells(dh.grid, cellid)
        dofs = celldofs(dh, cellid)
        n_nodes = length(cell.nodes)
        per_node = length(dofs) ÷ n_nodes
        for (i,node) in enumerate(cell.nodes)
            range_dofs = ((i-1)*per_node+1):(i*per_node)
            map[node] = get(map,node,Int[]); append!(map[node], dofs[range_dofs])
        end
    end
    for node in keys(map)
        sort!(map[node]); unique!(map[node])
    end
    map
end

node_dof_map = node_to_dof(dh)

# Apply nodal force safely
f_ext = zeros(ndofs(dh))
for node_id in getnodeset(grid,"nodal_force")
    dof_y = node_dof_map[node_id][2]  # y-direction
    if dof_y ∉ ch.prescribed_dofs
        f_ext[dof_y] = -1.0
    end
end

# Apply BCs and solve
apply!(K, f_ext, ch)
u = K \ f_ext

println("Max Displacements: ", maximum(abs.(u)))

[KnutAM]

Nice work! There are some potential improvements to the above code that the following should fix:

1. Making it work with multiple fields in the dofhandler
2. Throwing an error for non-isoparameteric elements
3. Except when the field is only on part on the domain, we always have the full number of nodes. So instead of using Dict{Int, Vector{Int}}, we can just index a Matrix{Int} directly as [dimnr, nodenr] (this reduces allocations and gc load, and will also make lookup faster (although probably not important))

Ferrite.get_n_copies(::Interpolation) = 1 # Extend internal function to work for non-vectorized interpolations

function node_to_dof(dh, fieldname)
    grid = Ferrite.get_grid(dh)
    ip = Ferrite.getfieldinterpolation(first(dh.subdofhandlers), fieldname);
    dofmap = zeros(Int, Ferrite.get_n_copies(ip)::Int, getnnodes(grid));
    for sdh in dh.subdofhandlers
        node_to_dof!(dofmap, sdh, fieldname)
    end
    return dofmap
end

function node_to_dof!(dofmap::Matrix{Int}, sdh, fieldname)
    if geometric_interpolation(getcelltype(sdh)) != Ferrite.get_base_interpolation(Ferrite.getfieldinterpolation(sdh, fieldname))
        throw(ArgumentError("A node to dof map is only possible for isoparametric elements (same function and geometric interpolation)"))
    end
    grid = Ferrite.get_grid(dh)
    ndofs_per_node = size(dofmap, 1)
    for cell in CellIterator(sdh)
        for (i, nodeidx) in enumerate(getnodes(cell))
            for j in 1:ndofs_per_node
                dofidx = dof_range(sdh, fieldname)[ndofs_per_node * (i - 1) + j]
                dofmap[j, nodeidx] = celldofs(cell)[dofidx]
            end
        end
    end
    return dofmap
end

In the case you've tested above, the result is the same so I don't think this will actually make it work differently. Could you elaborate what you mean by "doesn't work" for Triangle cells?

[Aminofa70]

@KnutAM brilliant, thanks.
I tested for tri3, and still get wrong result.
When "nodal_force" shares nodes with "support_x", the result is not accurate.
For example, the following code for MMB beam
MATLAB: max(abs(U)): 12.8208
Ferrite Julia: Max displacement: 13.044649865617167

When nodal force does not share nodes with "support_x" and "support_y", the result is accurate.
For quad4, we don't have this problem.

using Ferrite, SparseArrays
##
# Material
E, ν = 1.0, 0.3
C_voigt = E * [1.0 ν 0.0; ν 1.0 0.0; 0.0 0.0 (1-ν)/2] / (1-ν^2)
C = fromvoigt(SymmetricTensor{4,2}, C_voigt)

# Mesh
Lx, Ly = 1.0, 1.0
nx, ny = 10, 10
corners = [Vec((0.0,0.0)), Vec((Lx,0.0)), Vec((Lx,Ly)), Vec((0.0,Ly))]
grid = generate_grid(Triangle, (nx,ny), corners)

# Dirichlet BCs
addfacetset!(grid, "support_x", x->x[1] ≈ 0.0)  # 
addnodeset!(grid, "support_y", x->x[1] ≈ Lx && x[2] ≈ 0.0) 

# # Nodal force
addnodeset!(grid, "nodal_force", x->x[1] ≈ 0.0 && x[2] ≈ Ly)


# DofHandler
dh = DofHandler(grid)
add!(dh, :u, Lagrange{RefTriangle,1}()^2)
close!(dh)

# CellValues
qr = QuadratureRule{RefTriangle}(2)
cv = CellValues(qr, Lagrange{RefTriangle,1}()^2)

# Allocate stiffness
K = allocate_matrix(dh)

# Assemble element stiffness
function assemble_cell!(ke, cellvalues, C)
    for q in 1:getnquadpoints(cellvalues)
        dΩ = getdetJdV(cellvalues, q)
        for i in 1:getnbasefunctions(cellvalues)
            ∇Nᵢ = shape_gradient(cellvalues, q, i)
            for j in 1:getnbasefunctions(cellvalues)
                ∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q, j)
                ke[i,j] += (∇Nᵢ ⊡ C ⊡ ∇ˢʸᵐNⱼ) * dΩ
            end
        end
    end
end

function assemble_global!(K, dh, cellvalues, C)
    ke = zeros(getnbasefunctions(cellvalues), getnbasefunctions(cellvalues))
    assembler = start_assemble(K)
    for cell in CellIterator(dh)
        reinit!(cellvalues, cell)
        fill!(ke, 0.0)
        assemble_cell!(ke, cellvalues, C)
        assemble!(assembler, celldofs(cell), ke)
    end
end

assemble_global!(K, dh, cv, C)

# ConstraintHandler
ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfacetset(grid,"support_x"), (x,t)->0.0, 1))
add!(ch, Dirichlet(:u, getnodeset(grid,"support_y"), (x,t)->0.0, 2))
close!(ch)

# Map node -> DOFs
Ferrite.get_n_copies(::Interpolation) = 1 # Extend internal function to work for non-vectorized interpolations
function node_to_dof!(dofmap::Matrix{Int}, sdh, fieldname)
    if geometric_interpolation(getcelltype(sdh)) != Ferrite.get_base_interpolation(Ferrite.getfieldinterpolation(sdh, fieldname))
        throw(ArgumentError("A node to dof map is only possible for isoparametric elements (same function and geometric interpolation)"))
    end
    grid = Ferrite.get_grid(dh)
    ndofs_per_node = size(dofmap, 1)
    for cell in CellIterator(sdh)
        for (i, nodeidx) in enumerate(getnodes(cell))
            for j in 1:ndofs_per_node
                dofidx = dof_range(sdh, fieldname)[ndofs_per_node * (i - 1) + j]
                dofmap[j, nodeidx] = celldofs(cell)[dofidx]
            end
        end
    end
    return dofmap
end

function node_to_dof(dh, fieldname)
    grid = Ferrite.get_grid(dh)
    ip = Ferrite.getfieldinterpolation(first(dh.subdofhandlers), fieldname);
    dofmap = zeros(Int, Ferrite.get_n_copies(ip)::Int, getnnodes(grid));
    for sdh in dh.subdofhandlers
        node_to_dof!(dofmap, sdh, fieldname)
    end
    return dofmap
end

node_dof_map = node_to_dof(dh, :u)   # size: (ncomp, n_nodes)
n_nodes = getnnodes(grid)
ncomp = size(node_dof_map, 1)        # e.g. 2 for ux, uy

f_ext = zeros(ndofs(dh))   

for node_id in getnodeset(grid, "nodal_force")
    for comp in 1:ncomp
        dof = node_dof_map[comp, node_id]   #
        if dof ∉ ch.prescribed_dofs
            if comp == 2   
                f_ext[dof] = -1.0                  
            end
        end
    end
end

apply!(K, f_ext, ch)   # 
u = K \ f_ext

println("Max displacement: ", maximum(abs.(u)))

[KnutAM]

It is a bit hard to know if/why it is different to your? matlab implementation.

(There should be no need to do the filtering: checking if dof in ch.prescribed_dofs, but changing this should make no difference.

To debug the differences you can check the free and constrained parts of your matrices, e.g.

cdofs = ch.prescribed_dofs;
fdofs = ch.free_dofs;
Kff = K[fdofs, fdofs]; 
af = Kff \ f_ext[fdofs]; 

Although the numbering will be different, you can check that you have the same results and the same norm.

Also, please make sure you use the exact same mesh. If e.g. your triangle in the upper left corner is flipped you will get a different local stiffness (this is one reason why applying loads to a single node is normally discouraged).

[Aminofa70]

@KnutAM, I found the problem. The code now works for both tri and quad.
So, I believe that we can now add the functions to the src and provide some examples for it, no?

I agree with you that the traction is more realistic. For example, if you take a look at professional codes, such as FEbio or commercial codes like Abaqus, they all support nodal forces.
Why don't we think "one-stop shop" for Ferrite? no?
What do you think?

[KristofferC]

I found the problem. The code now works for both tri and quad.

It is mandatory to also say what the problem was in situations like this 😄

[Aminofa70]

@KristofferC
The problem was the mesh.
The mesh in MATLAB and Ferrite was different and I used the mesh in Ferrite for MATLAB.
Then result is the same. Now nodal force is applicable for both quad4 and tri3 and most probably others, Tet4 and Hex8.
It is gonna be practical to add it to Ferrite, isn't it?

[KnutAM]

What I'm still not sure about is

1. Should we provide the node_to_dof function that builds the complete map, or does it make more sense to have a utility that builds this for a given list of nodes (no need to materialize the entire matrix if you only need a single dof)
2. There has been discussions on doing something similar for e.g. facets (in that case with the coords, but similar would be relevant for the dofs), how should that interface look? (we have to think carefully about this, since we are at v1, we cannot implement something now and remove it later) (ref Face coordinates #702, cc: @lijas)
3. How do we avoid users jumping to this solution for e.g. applying Neumann boundary conditions, instead of doing it the proper way with FacetValues. (I foresee questions coming with why fext[boundarydofs_x] += force / length(boundarydofs) doesn't give the right result)

[termi-official]

I am also not sure if this is a good idea to do in first place. Thinking about this problem analytically, what you technically do is the following.

Your node $i$ has several surface elements attached to them. If we denote this patch with $S$ then the expression in your force vector and assuming you write 1 as your entry, then this expression would be $\int t_n \cdot \delta u_i dS = 1$, where $t_n$ is your traction vector and $\delta u_i$ the test function of node $i$. The issue I have here is that with mesh refinement we essentially make the surface area of this patch smaller, such that we obtain an infinitesimal surface in the limit (of e.g. infinite uniform mesh refinements). This implies that you apply an infinite force as $t_n$ goes to infinity. Is it still guaranteed that your underlying analytical problem does have a well-defined solution in that case and do you have some literature showing this?

[KristofferC]

Even though you, in the mesh refinement limit, apply what is effectively a Dirac delta function for the force I think you should still be able to get ok global behavior (for example displacements).

[KnutAM]

Even though you, in the mesh refinement limit, apply what is effectively a Dirac delta function for the force I think you should still be able to get ok global behavior (for example displacements).

In the case of linear problems, I think you are right that this should be "safe" when not looking too close around the load application point.

I am also not sure if this is a good idea to do in first place.

While, in principle, this is "a bad idea" from the variational form, there are some practical use-cases in e.g. topology optimization when the mesh is "fixed" (one could argue also here that it should be a field, but for practical purposes, many apply to a node). So I think supporting this is valid, but I just want to add it in a way that doesn't encourage users to jump to this solution in more cases than needed. It may appear easier when starting out (especially if coming from a structural mechanics background, when loads are indeed applied on nodes for beams/bars), but can come back to bite when doing e.g. a mesh convergence study as @termi-official mentions.