apply nodal force for unstruct mesh

[Aminofa70]

Dear Ferrite Developer,

In the following code, when I want to apply nodal force the results is incorrect.
For traction it works but for nodal force is incorrect.
is the function apply_nodal_force_3d! correct?
Could you kindly provide your guidance on this?
Thank you.

using Ferrite
using Comodo
function apply_nodal_force_3d!(grid, node_set_name, load_vector, f_ext)
node_set = getnodeset(grid, node_set_name)
for node_id in node_set
f_ext[3node_id-2] += load_vector[1]
f_ext[3node_id-1] += load_vector[2]
f_ext[3node_id] += load_vector[3]
end
end
function get_material_matrix_3d(E, ν)
C_voigt = E / ((1 + ν) * (1 - 2 * ν)) * [
1-ν ν ν 0 0 0;
ν 1-ν ν 0 0 0;
ν ν 1-ν 0 0 0;
0 0 0 (1-2ν)/2 0 0;
0 0 0 0 (1-2ν)/2 0;
0 0 0 0 0 (1-2ν)/2
]
return fromvoigt(SymmetricTensor{4,3}, C_voigt)
end
function assemble_cell_3d!(ke, cell_values, E, ν)
C = get_material_matrix_3d(E, ν)
for qp in 1:getnquadpoints(cell_values)
dΩ = getdetJdV(cell_values, qp)
for i in 1:getnbasefunctions(cell_values)
∇Ni = shape_gradient(cell_values, qp, i)
for j in 1:getnbasefunctions(cell_values)
∇δNj = shape_symmetric_gradient(cell_values, qp, j)
ke[i, j] += (∇Ni ⊡ C ⊡ ∇δNj) * dΩ
end
end
end
return ke
end
function assemble_global_3d!(K, dh, cell_values, E, ν)
n_basefuncs = getnbasefunctions(cell_values)
ke = zeros(n_basefuncs, n_basefuncs)
assembler = start_assemble(K)
for (cell_index, cell) in enumerate(CellIterator(dh))
reinit!(cell_values, cell)
fill!(ke, 0.0)
material_E = typeof(E) <: AbstractVector ? E[cell_index] : E
assemble_cell_3d!(ke, cell_values, material_E, ν)
assemble!(assembler, celldofs(cell), ke)
end
return K
end

Control parameters

Lx = 1.0
Ly = 1.0
Lz = 1.0
pointSpacing = 0.08

Dimensions for the box and spacing

boxDim = [Lx, Ly, Lz]

Generate tetrahedral mesh

E, V, Fb, CFb_type = tetbox(boxDim, pointSpacing)

Find extreme coordinates

min_x = minimum([v[1] for v in V])
max_x = maximum([v[1] for v in V])
min_y = minimum([v[2] for v in V])
max_y = maximum([v[2] for v in V])
min_z = minimum([v[3] for v in V])
max_z = maximum([v[3] for v in V])

Adjust vertex positions to start from the origin

V = [(v[1] - min_x, v[2] - min_y, v[3] - min_z) for v in V]

Recalculate extreme coordinates after adjustment

min_x = minimum([v[1] for v in V])
max_x = maximum([v[1] for v in V])
min_y = minimum([v[2] for v in V])
max_y = maximum([v[2] for v in V])
min_z = minimum([v[3] for v in V])
max_z = maximum([v[3] for v in V])
cells = [Ferrite.Tetrahedron((e[1], e[2], e[3], e[4])) for e in E]
nodes = [Ferrite.Node((v[1], v[2], v[3])) for v in V]
grid = Grid(cells, nodes)

addnodeset!(grid, "top", x -> abs(x[3] - (max_z)) < 1e-6)
addfacetset!(grid, "bottom", x -> abs(x[3] - (min_z)) < 1e-6)

function create_values()
order = 1
dim = 3
ip = Lagrange{RefTetrahedron,order}()^dim
qr = QuadratureRule{RefTetrahedron}(2)
qr_face = FacetQuadratureRule{RefTetrahedron}(1)
cell_values = CellValues(qr, ip)
facet_values = FacetValues(qr_face, ip)
return cell_values, facet_values
end

Function to create a DOF handler

function create_dofhandler(grid)
dh = Ferrite.DofHandler(grid)
Ferrite.add!(dh, :u, Ferrite.Lagrange{Ferrite.RefTetrahedron,1}()^3)
Ferrite.close!(dh)
return dh
end

Function to create boundary conditions

Function to create boundary conditions

function create_bc(dh, grid)
ch = Ferrite.ConstraintHandler(dh)
dbc = Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> [0.0, 0.0, 0.0], [1, 2, 3])
add!(ch, dbc)
Ferrite.close!(ch)
return ch
end

dh = create_dofhandler(grid)
ch = create_bc(dh, grid)

cell_values, facet_values = create_values()

loads = [0.0, 0.0, 1e4]
E = 210e6
ν = 0.3 # Poisson's ratio
Neumann_bc = "top"

Allocate and assemble the global stiffness matrix

K = allocate_matrix(dh)
assemble_global_3d!(K, dh, cell_values, E, ν)

Initialize external force vector

f_ext = zeros(ndofs(dh))
apply_nodal_force_3d!(grid, Neumann_bc, loads, f_ext)
apply!(K, f_ext, ch) # Apply constraints
u = K \ f_ext # Solve the system Ku = f_ext

DD = u;
display(maximum(DD))
VTKGridFile("my_solution", grid) do vtk
write_solution(vtk, dh, DD)
end;

VTKGridFile("boundary-conditions", dh) do vtk
Ferrite.write_constraints(vtk,ch)
end

[KnutAM]

is the function apply_nodal_force_3d! correct?

It seems to assume that dof number is 3(n-1) + d where n is the node number and d the component of the vector. This is not true. I assume you want to apply a traction vector, and in that case you have to integrate the contribution on the surface, see e.g. https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/linear_elasticity/#Boundary-conditions

[Aminofa70]

Dear @KnutAM

For the traction vector it works perfectly. Actually I want to apply a force to node (nodal force) not traction.
The function that I have provided directly apply force to force vector f_ext corresponding to degree of freedom.
I don't know if we can do it in Ferrite or not.

[KnutAM]

You can do it, but it is not straight-forward. See the discussion here for a script provided by @lijas (see the version for master, but not sure if this still works as there has been some breaking changes since then).

[Aminofa70]

You can do it, but it is not straight-forward. See the discussion here for a script provided by @lijas (see the version for master, but not sure if this still works as there has been some breaking changes since then).

Thank you so much @KnutAM I try to implement it.

[Aminofa70]

Dear @KnutAM
I have implemented this:
you comments and ideas are appreciated.
using Ferrite
function get_material_matrix(E, ν)
C_voigt = E * [1.0 ν 0.0; ν 1.0 0.0; 0.0 0.0 (1-ν)/2] / (1 - ν^2)
return fromvoigt(SymmetricTensor{4,2}, C_voigt)
end
function assemble_cell!(ke, cell_values, E, ν)
C = get_material_matrix(E, ν)
for qp in 1:getnquadpoints(cell_values)
dΩ = getdetJdV(cell_values, qp)
for i in 1:getnbasefunctions(cell_values)
∇Ni = shape_gradient(cell_values, qp, i)
for j in 1:getnbasefunctions(cell_values)
∇δNj = shape_symmetric_gradient(cell_values, qp, j)
ke[i, j] += (∇Ni ⊡ C ⊡ ∇δNj) * dΩ
end
end
end
return ke
end

function assemble_global!(K, dh, cell_values, E, ν)
n_basefuncs = getnbasefunctions(cell_values)
ke = zeros(n_basefuncs, n_basefuncs)
assembler = start_assemble(K)

for (cell_index, cell) in enumerate(CellIterator(dh))
    reinit!(cell_values, cell)
    fill!(ke, 0.0)

    # Check if E is scalar or vector
    local_E = isa(E, AbstractVector) ? E[cell_index] : E

    assemble_cell!(ke, cell_values, local_E, ν)
    assemble!(assembler, celldofs(cell), ke)
end

return K

end
function vertexdofs(dh::DofHandler, vertexid::VertexIndex)
cellid, lvidx = vertexid
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.n_components(sdh, ifield)
    field_ip = isa(sdh.field_interpolations[ifield], Ferrite.VectorizedInterpolation) ?
               sdh.field_interpolations[ifield].ip :
               sdh.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
function create_grid(Lx, Ly, nx, ny)
corners = [
Ferrite.Vec{2}((0.0, 0.0)), Ferrite.Vec{2}((Lx, 0.0)),
Ferrite.Vec{2}((Lx, Ly)), Ferrite.Vec{2}((0.0, Ly))
]
grid = Ferrite.generate_grid(Ferrite.Quadrilateral, (nx, ny), corners)
addnodeset!(grid, "support_1", x -> x[1] ≈ 0.0) # left side; Dirichlet BC
addnodeset!(grid, "support_2", x -> x[2] ≈ 0.0) # bottom side; Dirichlet BC
addnodeset!(grid, "force", x -> x[1] ≈ Lx) # right side; Neumann BC
return grid
end

Function to create CellValues and FacetValues

function create_values()
dim, order = 2, 1
ip = Ferrite.Lagrange{Ferrite.RefQuadrilateral,order}()^dim
qr = Ferrite.QuadratureRule{Ferrite.RefQuadrilateral}(2)
qr_face = Ferrite.FacetQuadratureRule{Ferrite.RefQuadrilateral}(1)
cell_values = Ferrite.CellValues(qr, ip)
facet_values = Ferrite.FacetValues(qr_face, ip)
return cell_values, facet_values
end
function create_dofhandler(grid)
dh = Ferrite.DofHandler(grid)
Ferrite.add!(dh, :u, Ferrite.Lagrange{Ferrite.RefQuadrilateral,1}()^2)
Ferrite.close!(dh)
return dh
end
function create_bc(dh)
ch = Ferrite.ConstraintHandler(dh)
Ferrite.add!(ch, Ferrite.Dirichlet(:u, Ferrite.getnodeset(dh.grid, "support_1"), (x, t) -> [0.0], [1]))
Ferrite.add!(ch, Ferrite.Dirichlet(:u, Ferrite.getnodeset(dh.grid, "support_2"), (x, t) -> [0.0], [2]))
Ferrite.close!(ch)
return ch
end
Lx, Ly = 1.0, 1.0 # Plate dimensions
nx, ny = 8,8 # Number of elements along x and y
grid = create_grid(Lx, Ly, nx, ny) # Generate the grid

Create DOF handler and constraints

dh = create_dofhandler(grid)
ch = create_bc(dh)

Create CellValues and FacetValues

cell_values, facet_values = create_values()

Find nodes on the right edge (x = Lx)

nodeid = Ferrite.getnodeset(grid, "force") # Nodes on the edge

Get the coordinates of the nodes

coords = [Ferrite.get_node_coordinate(grid, id) for id in nodeid]

Sort nodes based on their primary edge direction (assume edge in x or y)

Automatically detect if the edge is vertical (x constant) or horizontal (y constant)

is_vertical = all(abs(coords[1][1] - coord[1]) < 1e-8 for coord in coords)

if is_vertical
primary_direction = 2 # y-direction
else
primary_direction = 1 # x-direction
end

Convert OrderedSet to Vector for indexing

nodeid_vec = collect(nodeid)

Sort nodes based on the primary direction

sorted_indices = sortperm([coord[primary_direction] for coord in coords])
sorted_nodeid = nodeid_vec[sorted_indices]

Compute segment lengths (dy or dx) along the primary direction

sorted_coords = [coords[idx] for idx in sorted_indices]
dy = diff([coord[primary_direction] for coord in sorted_coords])

Define the total force to be applied along the edge

total_force = [16.0, 0.0] # Example: Fx = 5.0, Fy = 0.0

Initialize the global force vector

f = zeros(Float64, ndofs(dh))

Map forces to the DOFs of nodes

for i in eachindex(dy)
# Get the node IDs for the current segment
node1 = sorted_nodeid[i]
node2 = sorted_nodeid[i + 1]

# Calculate the segment's contribution to the total force
segment_force = total_force * dy[i] / sum(dy)  # Proportional to segment length

# Get the vertex indices for these nodes
vertex1 = nodeid_to_vertexindex(grid, node1)
vertex2 = nodeid_to_vertexindex(grid, node2)

# Get the DOFs associated with these vertices
dofs1 = vertexdofs(dh, vertex1)
dofs2 = vertexdofs(dh, vertex2)
# Distribute force equally to the DOFs of node1 and node2
f[dofs1[1:2]] .+= segment_force / 2  # Half the force to node1
f[dofs2[1:2]] .+= segment_force / 2  # Half the force to node2

end
@info sum(f)

E = 1.0
ν = 0.3 # Poisson's ratio

Allocate and assemble the global stiffness matrix

K = allocate_matrix(dh)
assemble_global!(K, dh, cell_values, E, ν)

Apply constraints and solve the system

Ferrite.apply!(K, f, ch)
u = K \ f # Solve the linear system
DD = u;
display(maximum(DD))
display(minimum(DD))
VTKGridFile("my_solution", grid) do vtk
write_solution(vtk, dh, DD)
end;
u_nodes = evaluate_at_grid_nodes(dh, u, :u)

Extract ux and uy components

ux = [vec[1] for vec in u_nodes] # First component of each vector
uy = [vec[2] for vec in u_nodes] # Second component of each vector

nothing

[Aminofa70]

and please refer to this in 3d
Thank you so much for your attention and please let me know you comments
I believe that I you could add this kind of load that is common to Ferrite, it will be interesting.
I will try to verify the code with books or some commercial codes

Convert node IDs (OrderedSet) to a Vector for indexing

node_id_vec = collect(node_ids)

Sort nodes based on the primary and secondary directions for consistent ordering

sorted_indices = sortperm([(coord[primary_direction], coord[secondary_direction]) for coord in node_coords])
sorted_node_ids = node_id_vec[sorted_indices]

Get the sorted coordinates for the nodes

sorted_coords = [node_coords[idx] for idx in sorted_indices]

Compute segment lengths in full 3D space

segment_lengths = [
norm([
sorted_coords[i + 1][j] - sorted_coords[i][j] for j in 1:3
]) for i in 1:(length(sorted_coords) - 1)
]

Define the total force to be applied along the edge (Fx, Fy, Fz)

total_force = [0.0, 0.0, 1.0e4] # Example: forces in x, y, z directions

Initialize the global force vector

f = zeros(Float64, ndofs(dh))

Map forces to the DOFs of nodes

for i in eachindex(segment_lengths)
# Get the node IDs for the current segment
node1 = sorted_node_ids[i]
node2 = sorted_node_ids[i + 1]

# Calculate the segment's contribution to the total force
segment_force = total_force .* (segment_lengths[i] / sum(segment_lengths))  # Proportional to segment length

# Get the vertex indices for these nodes
vertex1 = nodeid_to_vertexindex(grid, node1)
vertex2 = nodeid_to_vertexindex(grid, node2)

# Get the DOFs associated with these vertices
dofs1 = vertexdofs(dh, vertex1)
dofs2 = vertexdofs(dh, vertex2)

# Distribute forces equally to the DOFs of node1 and node2
for j in 1:3  # Iterate over x, y, z directions
    f[dofs1[j]] += segment_force[j] / 2  # Half the force to node1
    f[dofs2[j]] += segment_force[j] / 2  # Half the force to node2
end

end