Distributed Force Element Method

• Distributed Force Element Method is a force representation approach that replaces point loads with spatially resolved force elements integrated to yield global results.
• It applies across diverse fields such as tactile sensing, linear elasticity, incompressible flow, immersed boundary coupling, and beam theory, ensuring enhanced physical modeling.
• The methodology leverages discretization techniques like quadrature, finite elements, and spectral analysis to achieve controlled convergence, improved stability, and precise force balance.

Distributed Force Element Method denotes a family of force-centric formulations in which local tractions, interface loads, pressure stresses, or distributed source terms are represented explicitly on surface bins, boundary segments, patches, element edges, or embedded quadrature sets, and only then integrated to obtain resultant quantities such as total force, torque, pressure drag, or structural response. In tactile sensing, the method appears as mesh-independent aggregation of finite-element contact forces into spatial force bins; in linear elasticity, as quadrature-based discretization of singular interface forces; in incompressible flow, as a Green’s-function decomposition of wall pressure into distributed volume and surface sources; in immersed boundary coupling, as a distributed Lagrange multiplier; in equilibrium finite elements, as integrated traction degrees of freedom; in beam theory, as element-level integration of distributed forces and moments; and in surface force apparatus inversion, as recovery of local pressure from a measured global force curve (Sferrazza et al., 2019, Asghar et al., 2024, Zhou et al., 9 Sep 2025).

1. Conceptual scope and recurring definition

Taken together, these works suggest that DFEM is not a single universally standardized algorithm, but a recurring methodological pattern: replace point loads, purely nodal forces, or only integrated resultants by a spatially resolved force-element representation whose algebra respects the governing continuum mechanics.

Context	Force element	Primary objective
Vision-based tactile sensing	Equal-area surface bins carrying 3D force vectors	Learn distributed surface force from images
Linear elasticity with singular loading	Interface segments or patches with quadrature weights	Approximate $g\,n\,\delta_\Gamma$ stably
Incompressible flow	Convection volume elements and surface source elements	Attribute local wall pressure to flow structures
Immersed boundary coupling	Distributed Lagrange multiplier on the structure	Enforce fluid–structure kinematics variationally
Equilibrium FEM	Integrated traction DOFs on edges or faces	Enforce traction continuity and force balance
Geometrically exact beams	Cumulative distributed loads along the element axis	Solve element ODEs without equivalent nodal loads
Surface force apparatus inversion	Surface pressure elements linked by a kernel matrix	Convert global force to local pressure and density

The unifying operation is distribution followed by integration. In the tactile formulation, nodal contact forces are summed into bins and then to global force and torque. In elasticity with singular interface loading, a line or surface distribution is discretized by midpoint quadrature. In flow decomposition, local wall pressure is expressed as the sum of direct and scattered contributions from volumetric and boundary sources. In equilibrium finite elements and immersed-boundary formulations, the elemental or Lagrangian force field is the primary unknown, not a by-product of displacement recovery (Sferrazza et al., 2019, Asghar et al., 2024, Boffi et al., 2014, Olesen et al., 2016).

2. Core mathematical structures

A first canonical form is the surface-traction representation. For the tactile sensor, the contact traction on the top surface $\Gamma_c$ is

$t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$

and Abaqus returns discrete nodal contact forces $\{f_j\}$ on the surface mesh. These are aggregated into equal-area bins,

$F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$

with total contact force and torque obtained by summation:

$$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$

This creates a compact force-element label space that is independent of the FE mesh (Sferrazza et al., 2019).

A second form is the distributional interface-force model in linear elasticity. For a smooth closed interface $\Gamma\subset\Omega$, the force density is written

$-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$

with the weak-form load contribution

$$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$

DFEM then assembles the load vector by interface quadrature rather than by singular nodal forcing:

$$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$

Because the displacement due to Dirac-supported loading is singular on $\Gamma_c$0, the exact solution is not in $\Gamma_c$1 for $\Gamma_c$2; the rigorous convergence statement is therefore formulated away from the support, with

$\Gamma_c$3

The quadratic order comes from midpoint-rule accuracy on the interface and stability of the elasticity operator (Asghar et al., 2024).

A third form is pressure-source decomposition in incompressible flow. The pressure field is split as $\Gamma_c$4, where $\Gamma_c$5 solves a Poisson equation driven by convective acceleration and $\Gamma_c$6 is harmonic in the fluid with Neumann data induced by boundary vorticity and boundary acceleration. The distributed wall pressure is then written as the sum of three source classes,

$\Gamma_c$7

with $\Gamma_c$8 further decomposed into direct radiation and boundary scatter. The wall traction is

$\Gamma_c$9

Here the force elements are not structural DOFs but source-attribution elements: convection volume elements, vorticity surface elements, and acceleration surface elements (Zhou et al., 9 Sep 2025).

A fourth form is the variational interaction-force density used in immersed boundary formulations. The distributed Lagrange multiplier $t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$0 lives on the structure in Lagrangian coordinates and appears in distributional form as

$t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$1

This replaces pointwise interpolation and spreading by FE pairings $t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$2 and $t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$3, so that the coupling force is distributed through the same variational operators that enforce the kinematic constraint (Boffi et al., 2014).

A fifth form appears in force-based finite elements. In the higher-order equilibrium method, the DOFs are integrated traction components along edges,

$t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$4

and discrete equilibrium is expressed topologically as

$t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$5

In the geometrically exact beam element, the distributed loads $t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$6, $t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$7, and $t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$8 enter through cumulative resultants

$t(x)=\sigma(x)n(x), \qquad x\in\Gamma_c,$9

so internal forces and moments are explicit functions of integrated distributed loading rather than equivalent nodal loads (Olesen et al., 2016, Jirasek et al., 2024).

3. Discrete realizations and numerical machinery

The tactile formulation uses a high-fidelity forward FE model of a multilayer elastomer sensor with incompressible hyperelastic Ogden materials of order $\{f_j\}$0, large deformation, frictional contact, and mixed-element discretization. Abaqus/Standard with NLgeom solves the nonlinear residual by Newton–Raphson, using C3D4H tetrahedra and C3D8RH hexahedra with local refinement of characteristic size $\{f_j\}$1 mm, tie constraints at the layer interface, hard normal contact, and isotropic Coulomb friction with constant coefficient $\{f_j\}$2. The resulting nodal contact forces are binned into $\{f_j\}$3 equal-area force elements over the $\{f_j\}$4 surface (Sferrazza et al., 2019).

In the linear-elastic interface-force setting, the central discretization is the midpoint rule on a geometric partition of $\{f_j\}$5: line segments in 2D and triangular patches in 3D. The bulk mesh and the interface discretization can be independent, because one only needs the evaluation of shape functions $\{f_j\}$6 at interface quadrature points. The analysis recommends measuring errors on $\{f_j\}$7 and notes that higher-order quadrature could raise the convergence order away from $\{f_j\}$8 if $\{f_j\}$9 and $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$0 are sufficiently smooth (Asghar et al., 2024).

The flow-decomposition version uses the Method of Fundamental Solutions. Virtual sources are placed inside the solid body, boundary collocation points are used on $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$1, and the Green’s functions $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$2 and $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$3 are represented as linear combinations of free-space fundamental solutions. The coefficient system

$F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$4

is solved by truncated SVD to control ill-conditioning. Complexity is dominated by assembling $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$5 in $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$6 and by SVD in $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$7 or $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$8, depending on aspect ratio (Zhou et al., 9 Sep 2025).

In the immersed-boundary DLM formulation, discretization yields a monolithic saddle-point system coupling fluid velocity and pressure, structure position, and multiplier. The semi-implicit time step evaluates the kinematic constraint on the previous configuration $F_b=\sum_{j\in N_b} f_j\in\mathbb{R}^3,$9, treats viscosity implicitly, and uses the same configuration in both coupling forms, which is essential for the discrete energy argument (Boffi et al., 2014).

In the higher-order equilibrium finite element method, the numerical machinery is spectral. A primal grid of Gauss–Lobatto–Legendre points supports stress and traction variables via edge polynomials, and a dual grid of Gauss–Legendre points supports displacement and rotation multipliers. Inter-element traction continuity is strong because adjacent elements share integrated traction DOFs. Force balance depends on the incidence matrix rather than on metric interpolation, which is why equilibrium can be satisfied pointwise in the discrete polynomial space (Olesen et al., 2016).

The beam formulation uses finite differences and the shooting method. Because $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$0 and $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$1 depend only on midpoint values of $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$2 while $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$3 depends only on grid-point values of $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$4 and $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$5, central differences produce an explicit staggered update. Unknown end forces are recovered by Newton iteration on the shooting map from left-end generalized forces to right-end generalized coordinates (Jirasek et al., 2024).

The surface-force-apparatus transform is algebraically distinct but conceptually aligned. After variable changes, the measured force satisfies

$$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$6

where $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$7 is the local pressure on the cylindrical surface elements and $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$8, $$F_{\mathrm{tot}}\approx \sum_j f_j \approx \sum_b F_b,\qquad \tau \approx \sum_j (x_j-r_0)\times f_j.$$9 are matrices obtained from analytically integrated cell kernels. Because the integral equation is of the first kind, regularized inversion or SVD truncation is recommended (Hashimoto et al., 2015).

4. Representative domains and results

In tactile sensing, the FE-generated distributed force elements serve as supervision for a vision-based inverse map from optical flow to three-dimensional surface forces. The sensor uses a compact RGB camera imaging particles embedded in Ecoflex GEL; optical flow is computed by the Dense Inverse Search algorithm, averaged over $\Gamma\subset\Omega$0 equal-area image regions, and mapped by a feedforward DNN from a $\Gamma\subset\Omega$1-dimensional feature vector to a $\Gamma\subset\Omega$2-dimensional force label. The network has fully connected hidden layers of sizes $\Gamma\subset\Omega$3 with sigmoid activations and dropout $\Gamma\subset\Omega$4, and is trained with Adam, learning rate $\Gamma\subset\Omega$5, batch size $\Gamma\subset\Omega$6, on $\Gamma\subset\Omega$7 indentations with an $\Gamma\subset\Omega$8 train/test split. Validation against a six-axis ATI Mini 27 Titanium sensor yields RMSE between FEM total forces and F/T measurements of $\Gamma\subset\Omega$9 N in $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$0, $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$1 N in $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$2, and $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$3 N in $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$4, matching sensor resolution. Test-set force-distribution errors are RMSE $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$5 N, $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$6 N, and $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$7 N per bin in $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$8, $-\nabla\cdot\sigma(u)=f+g\,n\,\delta_\Gamma,$9, and $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$0, and end-to-end inference runs at $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$1 Hz on a dual-core $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$2 GHz laptop CPU without GPU acceleration (Sferrazza et al., 2019).

In singularly forced elasticity, the key result is analytical rather than empirical: for a superposition of midpoint-applied point forces converging to an integral force density on a cell boundary or interface, the displacement error away from the support converges in $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$3 with the same quadratic order as the midpoint rule. The paper verifies this in $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$4 and $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$5 for circular and spherical configurations and also records deterioration when evaluation domains intersect the singular support, exactly where the assumptions of the theorem no longer hold (Asghar et al., 2024).

In incompressible flow, DFEM is used to attribute instantaneous wall-pressure distributions to distinct physical source families. For laminar flow over a stationary circular cylinder at $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$6, the method reconstructs wall-pressure distributions in close agreement with direct Navier–Stokes solutions; negative pressure near the rear sides is primarily due to the direct radiation of convection volume sources, and pressure fluctuations and lift fluctuations are dominated by volume sources. For an oscillating circular cylinder at $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$7, a pronounced destructive cancellation between acceleration surface elements and convection volume elements suppresses lift fluctuation. For a sphere at $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$8, strong negative volume source in the boundary layer creates low-pressure regions and the separated shear layer exhibits near-zero pressure because the rotation and deformation contributions nearly cancel. For a turbulent sphere at $$\langle F_\Gamma,\phi\rangle=\int_\Gamma g\,n\cdot\phi\,dS.$$9, mean pressure coefficients and skin-friction profiles are reproduced consistently with prior DNS and experiments (Zhou et al., 9 Sep 2025).

In immersed boundary coupling, the principal result is unconditional stability of the semi-implicit fully discrete scheme with respect to the time step size. The discrete energy

$$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$0

satisfies a time-step-independent dissipation inequality. The same work also reports typical observed space-convergence rates of $$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$1 for velocity and $$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$2 for pressure when pressure is discontinuous across the immersed interface, and first-order convergence in $$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$3 for the Euler scheme (Boffi et al., 2014).

In planar linear elasticity, the higher-order equilibrium method reports force-equilibrium residuals of order $$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$4 to $$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$5 on orthogonal and curvilinear meshes, optimal $$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$6 rates for smooth solutions, strong traction continuity across interfaces by shared DOFs, and complementary energy convergence from above. On an L-shaped domain with a stress singularity, the force-based method maintains equilibrium to machine precision while a traditional displacement FEM exhibits large equilibrium residuals near the singular point (Olesen et al., 2016).

In beam mechanics, the finite-difference force-based element is tested on midspan point loading, uniform load, distributed moment, rigid-joint offsets, and buckling and post-buckling problems. The paper emphasizes quadratic convergence under refinement of the internal finite-difference grid, accurate handling of variable stiffness and rigid segments without increasing global DOFs, and a comparison between Reissner and Ziegler sectional models in relation to Haringx and Engesser stability theories (Jirasek et al., 2024).

In colloidal and interfacial measurement theory, the crossed-cylinder transform reconstructs a pressure profile $$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$7 from the measured force curve and then converts it to a normalized density profile through the rigid-wall relation

$$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$8

Here DFEM is an inverse force-distribution procedure: a global force measured by a surface force apparatus is distributed back onto local cylindrical surface elements through a geometry-specific kernel (Hashimoto et al., 2015).

5. Relation to adjacent methods and common misconceptions

A recurring misconception is that DFEM is synonymous with any inverse reconstruction of force from displacement or image data. The tactile work draws a sharper distinction: prior inverse FEM approaches, including GelSlim-related formulations, estimated surface forces by inverting a stiffness operator under linear elasticity, whereas the reported method uses forward FEM under large-deformation hyperelasticity with contact and friction to generate ground-truth distributed labels, and only then learns a fast inverse map at runtime (Sferrazza et al., 2019).

A second misconception is that force-element methods necessarily operate only at the level of resultant forces. The fluid-mechanics formulation is explicit on this point: the traditional force-element method decomposes the resultant aerodynamic force, but it cannot identify how specific flow structures map to local distributed pressure stress, which is the stated reason for introducing DFEM (Zhou et al., 9 Sep 2025).

A third misconception is that a distributed force representation removes singular behavior. The linear-elastic interface-force analysis proves controlled $$b_i=\int_\Gamma N_i\,g\,dS \approx \sum_{k=1}^N w_k\,g(s_k)\,N_i(s_k).$$9 convergence away from $\Gamma_c$00, but it also stresses that the exact solution is not in $\Gamma_c$01 for $\Gamma_c$02 and that evaluation on the support itself still requires regularization or weighted spaces (Asghar et al., 2024).

A fourth misconception is that distributed force elements are merely equivalent nodal loads under a different name. In the immersed-boundary DLM formulation, the multiplier is a variationally transferred interaction force density that enforces the kinematic constraint exactly in the chosen spaces; in the equilibrium finite element method, integrated tractions are primary DOFs shared across interfaces; and in the beam formulation, distributed loads enter the element ODEs continuously through $\Gamma_c$03, $\Gamma_c$04, and $\Gamma_c$05 rather than being condensed into nodal load vectors before the element solution (Boffi et al., 2014, Olesen et al., 2016, Jirasek et al., 2024).

These contrasts suggest that the most defensible general definition of DFEM is structural rather than nominal: a method belongs to the DFEM family when local force distributions are explicit unknowns, explicit labels, or explicit source terms in the governing formulation, and when global quantities are obtained by integrating those distributed objects.

6. Assumptions, limitations, and significance

The family resemblance across DFEM formulations does not erase substantial differences in modeling regime. The tactile method assumes incompressible hyperelastic Ogden materials, neglects strain-rate dependence, aging, and Mullins effect, uses a single constant friction coefficient, and limits spatial detail by a bin size of about $\Gamma_c$06 mm; extension to complex shapes and multiple contacts requires broader data collection (Sferrazza et al., 2019). The singular-elasticity formulation assumes small strains, homogeneous isotropic linear elasticity, smooth closed interfaces, and homogeneous Dirichlet conditions in the analysis; its convergence guarantee explicitly excludes a tubular neighborhood of the singular support (Asghar et al., 2024).

The fluid-pressure formulation assumes incompressible flow and a Poisson pressure description in the low-Mach “inner region,” and is numerically sensitive to velocity-gradient accuracy, source placement, source spacing, and SVD truncation. For truncated domains, a pressure correction may be needed, especially for the direct-radiation component (Zhou et al., 9 Sep 2025). The immersed-boundary DLM method introduces an additional field and coupling blocks, and its well-posedness depends on compatibility among $\Gamma_c$07, $\Gamma_c$08, $\Gamma_c$09, and $\Gamma_c$10; the full inf–sup analysis depends on the particular discrete pairing (Boffi et al., 2014).

The higher-order equilibrium method satisfies moment equilibrium only weakly when rotation is placed on the Gauss–Legendre grid, while the stronger GLL placement can introduce singular kinematic modes on deformed grids (Olesen et al., 2016). The beam formulation is currently two-dimensional, uses linear sectional laws $\Gamma_c$11, $\Gamma_c$12, and $\Gamma_c$13 despite fully nonlinear kinematics, and relies on sufficiently smooth distributed loading for accurate finite-difference integration (Jirasek et al., 2024). The crossed-cylinder inversion is an ill-posed first-kind integral equation, so regularization is essential and the rigid-wall mapping is best aligned with hard-sphere-like interactions and long-cylinder geometry (Hashimoto et al., 2015).

Despite these differences, the literature consistently treats DFEM as a way to preserve physically meaningful locality in force representation. Whether the object being distributed is a contact traction field, an interface Dirac load, a wall-pressure source attribution, a Lagrange-multiplier interaction force, an edge traction DOF, or a beam-load resultant, the methodological payoff is the same: local force information becomes directly computable, learnable, or invertible, while resultant force and moment remain accessible by exact or controlled integration.