A comparison between the finite element method and a kinematic model derived from robot swarms for first and second gradient continua

In this paper, we consider a deformable continuous medium and its discrete representation realized by a lattice of points. The former is solved using the classical variational formulation with the finite element method. The latter, a 2D discrete “kinematic” model, instead is conceived to determine the displacements of the lattice points depending on interaction rules among them and thus provides the final configuration of the system. The kinematic model assigns the displacements of some points, so-called leaders, by solving Newton’s law; the other points, namely followers, are left to rearrange themselves according to the lattice structure and the flocking rules. These rules are derived from the effort to describe the behaviour of a robot swarm as a single whole organism. The advantage of the kinematic model lies in reducing computational cost and the easiness of managing complicated structures and fracture phenomena. In addition, generalizing the discrete model to non-local interactions, such as for second gradient materials, is easier than solving partial differential equations. This paper aims to compare and discuss the deformed configurations obtained by these two approaches. The comparison between FEM and the kinematic model shows a reasonable agreement even in the case of large deformations for the standard case of the first gradient continuum.

The origin of this idea lies in the Computer Graphics industry [2]; in video games, it is necessary to calculate the deformation and breakage of objects very quickly. The body is discretized, and PBD tries to obtain a plausible description of its deformation, considering only the relative positions of the points to describe the action of the internal forces. It is imperative in the visual effects industry to provide movies with simulated versions of complex physical phenomena that are too expensive or impossible to solve precisely. In such a commercial environment, these simulations must be as efficient and fast as possible while remaining realistic to the viewer's eye.
Initially, approaches to simulating dynamic systems in computer graphics were often based on force [3] and Newton's equations. Internal forces (i.e. elasticity) and external forces (i.e. gravity) are computed to obtain the displacements; they are derived from Newton's law using the discrete time integration method to update the velocities, from the accelerations, and finally, the positions of the discretized object.
This classical approach presents many problems which are typical for numerical integration systems. Moreover, collision constraints and penetrations are difficult to manage. Also, while external forces are relatively easy to apply, problems occur in calculating internal forces. Sometimes these forces are represented by a massless linear spring (spring model); this model constitutes an important class of models that can be characterized in different ways by changing the geometry of the springs and their type. Some difficulties may arise when one has to separate bending stiffness from torsional stiffness. It is possible to establish micro-macro relationships to link spring constants with the constitutive parameters of bulk materials (i.e. Young's modulus and Poisson's ratio). Some other models to describe internal forces may be unstable or require high computation time.
Another possibility involves the use of energy minimization methods and the actions principles [14]. These are guided by variational principles, which state that a constrained surface will take the form that minimizes the total strain energy. Unfortunately, differential equations are still needed to calculate the force by expressing the energy of a surface in terms of its local deformation and differentiating the energy with respect to position. The problem is even more difficult in the case of the second gradient, where the energy is written in a more complex way, and Lamé's two constants increase by five. This because in the case of strain gradient theory of a centro-symmetric and isotropic material, there are 5 additional constants (so they are 7 instead of 2 in the simplest case).
Another possibility emerged in the 1970s with smoothed particle hydrodynamics (SPH) [15]. But in this case one always starts with the Navier-Stokes equations to be discretized and thus the computational cost is high.
We emphasize how in computational sciences the main focus is on accuracy. In contrast, the main issues in the Computer Graphics industry are stability, robustness, and computational speed, keeping the results plausible. So far, a new method has been developed in which it is preferable to have direct control over the positions of objects, often represented by the vertices of a mesh, avoiding collision between them and managing their displacements instead of dialoguing with forces. The method is known as Position-Based Dynamic (PBD) and was introduced by Jakobsen in 2001 [16]. It works directly on positions, which makes manipulations and deformations easy to manage. Moreover, with the PBD approach, overshooting and energy problems, related to explicit integration can be avoided because the integration can be directly controlled. The advantage of this technique, even over a mass-spring model, is its stability, even when taking large time steps, which is not always possible in the mass-spring model. The new position can be expressed as an algebraic function of the old position and the time step. The algorithm is iterative and can be stopped at any time if a certain inaccuracy is accepted: If the time step is larger the algorithm can go faster and thus a time/accuracy trade-off must be chosen. Problems with internal forces are avoided by replacing them with a set of geometric constraints. In the simple case, the constraints are holonomic and bilateral. To simulate elastic behaviour, for example, the massspring model can be converted to a constraint system, which requires the distance between connected nodes to be equal to the initial distance. This system can be solved sequentially and iteratively, moving nodes to satisfy each constraint until sufficient accuracy is achieved. One of the main limitations may be the dependence of the results on the number of iterations of the solver. Still, one of the significant advantages of a position-based approach is its controllability. Constitutive parameters of the material are hidden in the rules of interaction between first neighbours, hence in the constraints.
In some previous papers [17][18][19][20], we have shown how a 2D kinetic model is able to handle different types of phenomena, such as elasticity, plasticity, and fracture, using the same framework. The ultimate goal is not to obtain an exact description of the deformation, but a plausible visible description. This work aims to obtain an algebraic equation that can calculate the position of a point using the placement of other points as input.
This kinematic approach has the advantage of working with simple algebraic equations and the complexity of the system increases only linearly, rather than quadratically, as the number of points increases. We want to emphasize now that, in this paper, we use an approach that is only partially similar to PBD because, so far, we only use position data.
In the context of underwater swarm robotics, the problem related to flocking rules to determine what the behaviour of a single robot should be so that the swarm reaches a certain assigned configuration is often faced. One option to accomplish this task is to consider what a given number of its neighbours are doing. Using this information, a pure kinematic model can be developed with which we try to reproduce some behaviours of two-dimensional deformable bodies according to the standard Cauchy model as well as to the second gradient theory. All the constitutive properties of the material are within the rules determining the kinematic displacements of the particles, each with respect to the others. This tool has an advantage in terms of computational cost and is very flexible to be adapted to samples of complex geometry and different behaviours. The results are encouraging and even fracture can be considered easily and in many different ways.
The other PBD methods we found in the literature require knowledge of the particle velocity, are related to Newton's equation, and use mass explicitly. We avoid, for now, the explicit use of these concepts by using a purely positional kinematic approach. Therefore, in this work, we refrain from considering basic aspects such as mass, velocity, and conservation of momentum but only evaluate kinematic characteristics and transformation rules. The most crucial problem, then, is how to reconnect the classical constitutive equations, which determine the characteristics of the model and the material, to the parameters of the kinematic model under study.
The proposed model is based on swarm intelligence theory [21]. The algorithm that exploits is parallelizable and runs with the CUDA programming language, so we can obtain numerical simulations in a few seconds, even if they are very complex, using a modest desktop PC. Remark that the model does not yet have any physical meaning. It concerns only a set of variables, which are points in a Cartesian plane, that change coordinates in successive steps by iterations of an algorithm. Therefore, physical concepts such as "motion", "time" or "velocity" are not considered at this stage of the model development. It will be the subject of future work to connect the model parameters with the constitutive equations of the material.
This tool could be helpful in complex microstructures not easily analysed by the Cauchy continuum [22]. Classical Cauchy continua cannot accurately predict highly inhomogeneous microstructures [23][24][25][26][27] and [28]. However, generalizations must be introduced by considering additional degrees of freedom to account for kinematics at the microstructure level or by including higher displacement gradients than the first one in the strain energy density [29][30][31][32][33][34][35][36]. The latter is particularly relevant when considering the technological interest in developing exotic mechanical metamaterials capable of performing targeted tasks; therefore, investigating new and efficient algorithms, such as our tool, is of great interest.
In the second gradient theory, the Lagrangian function also depends on the second derivatives of the displacement field and not only on the first ones. The origin of this model lies in the fact that higher-order derivatives have sometimes been included to obtain more general models; Lagrangian functions dependent on higher-order derivatives are often called non-local in the literature [37][38][39][40][41][42][43][44][45][46][47][48][49][50]. This terminology can be explained by considering that, while for a standard material, in a discrete framework with a lattice of points, a simple interaction between a given point and the immediately neighbouring points, in the continuum limit involves an energy function based on the first derivative of the placement field according to the Piola's ansatz. In the case of the nth gradient, the interaction among points involves all the nth frame neighbouring points in the lattice. Therefore, in the continuum limit, the energy depends on the nth derivative of the placement [51][52][53][54][55]. Indeed, in practice, one must measure finite differences of higher order in a discrete framework to estimate the values of the higher derivatives. Furthermore, by introducing new characteristic lengths, higher-order models also introduce new parameters that require appropriate methods for their determination [56][57][58][59][60][61][62].
Sometimes remote actions are excluded, which has led to the introduction of local Lagrangians; however, in some situations, non-local Lagrangians must be considered. For example, in composite materials where microfibers are present, these can introduce non-local interaction between points that are not close in space due to the connection between distant points provided by the fibres.
We are aware that many mathematical results are challenging to satisfy, but we believe that the proposed discrete algorithm can be helpful in a wide variety of mechanical applications.
The paper is organized as follows. The first section contains a description of the swarm robot model. The second section shows numerical results comparing a 2D continuum with a first and second gradient energy and the proposed kinematic model. Future developments on this topic and open questions are presented in the conclusions.

The tool
The use of robot swarms in underwater work is becoming increasingly important; flocking rules have been applied to determine the displacement of the individual robot in order to reach the assigned final configuration of the entire swarm [63]. The new position of the robot is determined by the position of its neighbours; we applied the same algorithm to describe the deformation of a lattice of material particles. This is a purely kinematic approach to the problem, without the use of the concepts of mass or force.
In practice, we proposed to consider a transformation operator, with constraints, between matrices representing the particle configuration, C t , for a discrete set of time steps t 1 t 2 ,...t n. This operator can be represented by a fourth-order tensor that transforms a two-dimensional matrix into another two-dimensional matrix. The continuous two-dimensional body is discretized into a finite number of particles lying, in the undeformed configuration, in the lattice node. A complete description of the algorithm can be found in [17,18,20] where we obtained a plausible description of success even in the case of material fracture. In the previous articles, we have defined four types of particles that make up the lattice but, if necessary, new types can be introduced, to describe other properties, due to the modular structure of the algorithm. They were the leaders whose movement is assigned and determines the movement of the other particles. The followers whose movement is determined by the flocking rule with other particles; we also had the frame, introduced to avoid edge effects and the fourth type concerns fracture These last two will not be used in this paper, where. a model without the frame avoiding edge effect with different flocking rules is considered.
Some hints from graph theory help us to note that time and velocity are concepts that do not belong to this model; therefore, we prefer to talk about pseudo-time as an update step of the algorithm.
We then consider a graph as a mathematical structure used to model pairwise relationships between objects. We denote a graph G = (V, L) as an ordered pair of a set P of points, a set L of links. We will consider only undirected graphs, so there is no direction associated with the links. Considering n points V i , with i ∈ {1, 2, ..., n}. We define the adjacency matrix as: the matrix A has dimension n × n. This is not the only possible definition of adjacency matrix; to satisfy some properties of our swarm we can change it, i.e. change the definition of neighbour.
A path on the graph is defined as a finite or infinite sequence of links joining a sequence of points; it will be useful to define the concept of neighbourhood order, related to gradient and non-local order.
To describe the time evolution of our swarm we consider a set of p-time (pseudo-time) steps, i.e. a set of ordered discrete values of the time variable t that we will denote with T = {t 1 t 2 , ..., t n }. We talk about pseudotime because tit is only an index parameter of the loop. A swarm, denoted as S,is composed by n material points P i (i from 1 to n). At pseudo-time t m the i element of the swarm has position P i (t m ) = (x i , y i )(t m ) where x, y, are the coordinates in a reference system. So far P i is the material point and P i (t m ) is its position represented in Cartesian coordinate. At time t 0 , we denote by C 0 as the reference configuration of S.Therefore, S is the swarm as collection of points P i while the configurations C tm are the set of coordinates of them at pseudo-time t m . We assume that C 0 is the representation of a bidimensional crystalline lattice to be specified. For each elementi in Swe define a fixed set of neighbour's points, calculated in the reference configuration C 0 . Their number, typically, is nc, the coordination number of the lattice. We call this set of point the first neighbours. The process is iterative so we can define the neighbours of each neighbours and call this set the second neighbours.
The first neighbours of a point P i so are all the points of S at time t 0 adjacent to the point P i in the meaning of the graphs (so it is related to the lattice coordination number). The second neighbours of a point P i are all the points of S whose minimum path without repetitions consists of only two links. Therefore, defining l(i, j) as the minimum path without repetitions between points i and j, the k-th neighbours of a point i is the set So, Ne k (i)is the set of point neighbours of i order k. We call it the shell of order k. We can also consider pairs of opposite neighbours of chosen point; this will be useful in defining the rule to compute an alignment term. The set of pairs k-order opposite neighbours of a point P i is the set Therefore, we have defined, assigned a point P i and its kneighbours a couple of point composed by the kneighbours and its opposite neighbours h. So, the j-th pair of the point P i will be indicated with ρ i j and its average point of the pair as < ρ i j > for each t m . Note that Ne k (i)is the i-th row of the adjacency matrix: this means that our adjacency matrix A is not nxn but nxm where m is the number of the total neighbours considered until the k-th order.
We select a set of leaders L; in analogy to biological swarms they are a subset, L ,of C 0 whose trajectory is determined. These elements drive the evolution and the deformation of the swarm and they have a "priori" imposed motion and are not influenced by other elements in the swarm.
All the remaining elements of C 0 form the set of the Followers, F. The movement of the follower is guided by the flocking rules. A possible choice of the rules is: Where Int i (t m ) is an interaction term between the point i and its neighbours; one of the possibilities is to choose the interaction term so Where k is the number of neighbours of the point i.
Where the quantities are explained as follow.
is the new position of the point P i . Int i (t m )is the Interaction term among the point P i and all its neighbours in the h-th order shell until the k-th order. N i is the number of the neighbours of the point P i . r , j (i,t m ) is the new position of the j-th neighbours of the point P i at the p-time step t m . α j is the coefficient that modulates the intensity of the "elastic" interaction (W i, j ) and the convergence velocity. It depends by the j-th neighbour's shell order. β l is the coefficient that modulates the intensity of the "alignment" interaction (Y ,i,l ). It depends by the j-th neighbour's shell order; the lis the label of the couple ρ i,l . <> has the meaning of simple average value of the distances, The term W i, j, can simulate an elastic behaviour of the material and the rule contains a dependence on the distances between the point and its neighbours trying to emulate a similarity with Hooke's law. We remark that the proposed algorithm is still different from a particle system with spring interactions and some constrained, because we do not give a constitutive equation for elastic potential. We implement a recursive set of rules in order to define the evolution. It should be noted how in Eqs. 5 and 6 although the interaction is between two particles it occurs for all neighbours, so this would be a 'multi-body' Hooke's law.
The Y i,l term was introduced to handle two problems. The first is to obtain a way to modulate the bending response of the swarm. Indeed, the Y i,l term pushes the point i to be located in the centroid of the considered pair of opposite neighbours. It is an alignment term but we will not always consider the Y i,l term in this paper using other methodologies. The second problem is to avoid the overlap of points during displacement.
In previous articles we have supposed the existence of a frame, but not here. This is because, unlike [17,18], to balance the interaction on the edges of the swarm we do not need all shells to be complete (i.e. each shell has the maximum number of neighbours expected for the considered order). When used, Frame points have their own rule of motion and are moved in a different pseudo-time than the other points in the swarm. In their absence therefore, here all non-Leader points are moved simultaneously and this makes the algorithm leaner, smoother and faster. However, we noticed that the edges of the swarm suffer from local overlaps and oscillations; they cause exotic behaviours that propagate within the whole lattice, especially increasing the term α j Thus, our idea about the incompleteness of neighbouring shells was wrong, and to correct it, we correct the term 1/N i with 1/(2Z k − N k (i)), in Eq. 6, where Z k is the number of points in the shell of order k when it is complete.
The boundary elements may have a larger displacement than the inner followers, due to the lack of points at the edge of the swarm not balancing the interaction. This because sometimes a border point can have a neighbour's number less than k. This reduces the weight of the edge points and coregulates the imbalance. This modified rule will be used thereafter.
Several methods have been used to avoid particle overlap and to overcome the alignment parameter Y i,l ; one of them is the use of intermediate propagation stages where the leaders are frozen and the lattice resettles.
In previous work [19] it has been observed how the deformation information in the swarm propagates through its elements as the pseudo-temporal steps move forward. The strain information will propagate through the swarm proceeding in shells that are proportional to the k-order of the neighbours considered in the interaction.
To avoid point's overlap, we introduce the parameter γ ; It consists in a certain number of settling steps while the leaders are stopped. In the previous model it was necessary to pay attention to the pseudo-velocity of the leaders, assign them an extremely low value (resulting in longer computation times) to avoid the overlap problem (especially in compression simulations) and give enough time to the followers to adapt according to the deformations. We introduced algorithm cycles γ between two pseudo-time steps; in these intermediate cycles the leaders are stopped while the followers continue to resettle. This has the significance of changing the information about the rate of propagation of deformation within the swarm. What changes in the visible motion is the stiffness of the deformation but the final configuration obtained is the same.
With the parameter γ we can modulate how many shells are involved in a pseudo-temporal step; it is then related to the rate of propagation of the movie deformation in the swarm. This concept is analogous to the concept of correlation length in lattice models (such as the Ising model) and the concept of persistence length in polymer models.
It needs to find an optimum between the pseudo-speed (which can be higher) and the increase in machine time due to the increase in computational steps.
Practically, we start with a swarm occupying the nodes of a crystal lattice, move the leaders to an assigned position, and calculate the positions of the followers by means of the above equations.
The code was initially developed in Mathematica by Wolfram research. It was later converted into the Python language to take advantage of the CUDA libraries to parallelize the calculation. The simulation involves particles (or agents) moving in two-dimensional space and interacting with each other according to certain rules or behaviour. The code imports several libraries, including matplotlib, NumPy, math, itertools, SciPy, timeit, numba and tqdm. Another section of the code includes several functions for creating the swarm. These functions include ADJMatrix, which creates an unordered adjacency matrix and is called by the lattice function, and NBORD, which reorders the elements of the adjacency matrix and is also called by lattice. The purpose of these functions is to create a network of connections between the agents in the swarm, allowing them to interact with each other. This is done by calculating the distance between each agent and its neighbours and creating a matrix representing these connections. The remaining code includes functions to simulate the behaviour of the agents in the swarm, updating their positions based on interactions with each other and the environment. The code can also include functions to visualize the swarm as it evolves over time using CUDA. It takes several input parameters, such as step size, start and stop times, number of particles, their initial positions and velocities, and various interaction parameters. The function then prepares the data to invoke a CUDA kernel, which is responsible for the calculation of the swarm's evolution. The kernel uses the input parameters to calculate the new positions and velocities of the particles, based on their interactions with neighbouring particles and the swarm leader. The function then returns the updated positions and velocities of the particles, as well as the time step at which the function finished its calculation. The kernel is structured in such a way that it performs the particle interactions in parallel on several GPU threads.

Results
In this section, we present some comparisons between our tool and the solution obtained by COMSOL software using finite elements analysis on a square lattice.
Let us consider a swarm S constituted by a finite number of elements, and indicate by C 0 the reference configuration of S a time t = t 0 . In C 0 the elements of S occupy the nodes of the lattice. We consider a set of time steps, i.e. a set of ordered discrete values for the time variable t that we will denote with T m = {t 1 , t 2 , t m }. Consider an orthonormal reference system with axes that are parallel to the lattice, whose unit length is equal to the side of the lattice and whose origin coincides with the left bottom vertex of the swarm. To identify the elements of the swarm, we use the coordinate (x i j , y i j )(t m ) of the element occupying the node (i, j) at the time t m . For each element (i, j), we associate a set of neighbours kth order as the -th coordination number of the lattice cell. By this, we define, in case of k = 1 and k = 2, the definition of first and second neighbours, including the neighbours of the first neighbours. Actually, we use a definition of k-th neighbours that is Lagrangian, i.e. the definition has memory of the identity of the neighbours of a given element of S. The reason is that we are interested in describing matter in solid state, where the behaviour of constitutive particles depends on local molecular interactions which are preserved during the evolution of the system. In the next papers, we shall relax this hypothesis.
So, we select a leader (or more than one) that is an element L placed at a node {i, j} of the lattice with an assigned displacement M, i.e.
That means a selection of points of S whose initial coordinates are a subset of C 0 .
In particular, the deformations we have compared are related to a rectangular specimen of size 1.2 m × 0.3 m that undergoes simple tensile stress with and without a slot (dimensions of the slot are 0.9 m × 0.075 m), posed in the centre of the sample. Points on the left side are clamped while the points on the right side of the sample are moved up to the final position which is assigned. Each time unit is composed by two steps. The movement of the leaders and the movement of the followers. No fracture is considered in this case. The simple rule, governing follower's motion is expressed by Eq. 4. This rule can be defined as "barycentre of its neighbours" in case of complete shell that means when all neighbours are considered. The neighbours are determined by the coordination number of the lattice; therefore, the leader's motion implies a displacement of the first layer that propagates in successive "time steps" to the other particles, according to the γ frames used. This means that the displacements, at each time step, involve a larger shell of points until it regards all the lattice points. The quotation marks are to remind that in our model the time is only fictitious because it is related to the iterative process. The model solved with COMSOL treats stationary case; therefore, the time is not important at all. This led us to two problems: how to link our model's parameters (α i , β i , γ ) with the constitutive parameters (e.g. ν, E) and how to simulate higher-order gradients continua.
In this section, we will show how we managed an empirical link between the Poisson coefficient ν with the only parameter α i in order to reproduce with tiny errors a first-order gradient continuum; meanwhile, for second-order gradient continua also β i plays an important role.
After that, we will show a quantitative comparison between our simulations and FEM solutions, an important step forward with respect to the previous papers [64] and [18].
We want to underline that in this paper we are not considering the fracture. Moreover, the dynamic evolution of the swarm is only a mean to obtain the final equilibrium configuration; in the future we shall discuss as to link this evolution with a real motion. Remember as the value of γ does not influence the final configuration as shown in [18], but only the intermediate steps.

First order gradient
In the previous papers [64] and [18], we stated that the interaction term W i,j is sufficient to describe qualitatively a first-order gradient deformation. In this subsection, we want to validate this statement showing that our model can describe first-order gradient deformation also quantitatively, comparing our results with those obtained by FEM analysis.
In [64], we noted that using all the α i with the same value α (Eq. 1) the Poisson effect does not exist In order to obtain Poisson's effect we can assign different values of α for each neighbour in all the shells of the swarm; by this way, we can vary the Poisson effect. In particular, in the case of a square lattice, we assign a value α for neighbours positioned in vertical and horizontal directions with respect to the centre point of the shell, and a different value, α d , for neighbours positioned in the diagonal direction.
Then we analysed the Poisson effect varying α d and keeping α fixed. We observed that for α d = 0 there is no Poisson effect; indeed, the displacement of each particle depends only on the interaction with vertical and horizontal neighbours. Therefore, the vertical and horizontal displacement are uncoupled. Increasing α d , keeping α fixed, the Poisson effect become more and more significant. Furthermore, we noted that varying both α d and α keeping their ratio fixed does not alter the Poisson effect (for values of the parameters in standard range as specified in [64]).
Thus, we hypothesize that only the ratio (α d / α) affects the intensity of the Poisson effect. Particularly, we expect that there is an increasing monotone function ν = ν (α d /α) because we are increasing the horizontal and vertical links. To test this hypothesis, we perform several simulations varying the ratio (α d /α) and evaluating "geometrically" the Poisson coefficient ν for each simulation. This means that we compute as Poisson coefficient the ratio of the transverse strain to the longitudinal strain of the specimen; the transverse strain is considered in the middle of the beam. We obtained the relation ν(α d /α) interpolating the collected data, see Fig. 1 (Consider the value of β = 0. Inverting this relation, we can evaluate the ratio (α d /α) that reproduces the desired Poisson coefficient. Remember as this relationship is purely empirical and there is no theoretical link between our parameter and Poisson coefficient. Forcing the coupling between the ratio (α d / α) to higher values, we are introducing a sort of anisotropy in our model showing Poisson's ratio larger than 0.5.
We underline that this analysis was performed only for swarms without the slot. This is because we try to simulate the same material with different shapes.

Second-order gradient
In previous papers [18,64], we supposed that second-order gradient phenomena could be obtained extending the W i,j interaction to second neighbours. We have observed that the Y i,l interaction is just enough to give good results, and we believe that it is a better candidate to describe second-order gradient phenomena with respect to the use of the second neighbours. Therefore, in the following we will consider only the W i,j and Y i,l interactions extended to the first neighbours in order to simulate the second gradient continua. No second neighbours will be considered in this paper. Since βis linked to three alignment points it is able to introduce non-local interaction albeit in a simple way. For these reasons, we explore the effect of β in reproducing second gradient phenomena.
As in the first gradient case, we performed several simulations varying the parameters (α, α d , β) considering the same βfor all the couples of first neighbours. For each simulation we valued geometrically the Poisson coefficient. In Fig. 1 is shown how β influences ν: fixed the ratio (α d /α), the Poisson coefficient decreases if β increases. A possible explanation could be that increasing β the system becomes stiffer, so Poisson effect is decreasing.
In this case, the relation ν = ν (α d /α, β) is not trivial: different triplets (α, α d , β) lead to the same Poisson coefficient but with different global deformations of the swarm. So, being unable to select the right triplets for now, we proceeded in heuristic way to find the triplets that better fit COMSOL simulations.
Also, in this case the magnitude of the Poisson effect does not depend directly on α, but the ratio (α d / α) affects the intensity of the Poisson effect, as shown in Fig. 1 for different β.

Comparison
To obtain a quantitative measure of the agreement between the results we consider the coordinates of a subset of points, N g , in the deformed configuration and we measure the distances between the corresponding points using the two methods, FEM and our tool. A colour plot with the differences shown as a very good agreement can be reached. Three cases are considered with three different nominal Poisson coefficients (0.2, 0.3 and 0.4). The swarm is composed of about 2500 elements; the pictures show a sub sampling of about 300 elements. The value of β is 0 in the case of the first gradient, while 1.3 in the second gradient model; we found that this value has the better match with the FEM solutions. The value of α/α d is given in Table 2. The motion is characterized by the features presented in Table 1.
We then modulated the parameters of our model so that the results coincided as closely as possible with those obtained by the FEM. We have considered two kinds of sample (with and without slot), three Poisson coefficients and first and second gradient cases. Finding the match was not easy, especially in the case of the second gradient, where the parameters involved are more than two owing to the presence of β. Having obtained very similar results, we can make the following considerations.
The differences between our tool solution and FEM are very small. In Table 2, there are max and mean errors, computed between the coordinates of the lattice, for all the cases studied. Now we will see where the differences are located within the sample. In the case of the first gradient, without slot, these are located in the lateral area (coloured yellow, see Fig. 2). In the centre of the sample, the differences are almost zero, while near the leaders, where the deformations are greater, the area is more critical. However, the location of the differences is not quite the same; increasing Poisson coefficient we have different areas where the mismatch is greater, also if they always are distributed close to the corners as can be seen from Fig. 2. Far from the leaders, where are concentrated deformation owing to edge effect and at the centre of the specimen, we note a uniform deformation, as it is expected in this kind of tests. In contrast, in the case of the second gradient, the main differences are always located in the same position, i.e. about one-quarter and three-quarters of the sample, and in the centre with respect to the vertical axis. See Fig. 3. This area corresponds to the location of maximum second gradient deformation. Also, in the case of the sample with the slot there are differences. It would seem that for the second gradient the difference between our method and the FEM solution is more pronounced close to the slot towards its centre; the slot is a critical area.
Generally, the agreement between models is greater in the first gradient case. But there is also a qualitative fact. Take for example the case of the sample with the slot in the second gradient condition. For high values of Focusing our attention on a particular area of our sample, in the worst case with the highest Poisson's modulus, we can see the quantitative and qualitative differences described above, see Fig. 6. Moreover, look now at this interesting phenomenon, in Table 2. In the FEM analysis, we assigned the Poisson coefficient and performed the simulation; if we measure "a posteriori" the Poisson coefficient from its geometrical meaning, as a variation of the sample width, we get a rather different value from the nominal one used in the calculation. This is because the deformation is not purely uniaxial. Again, the agreement with our model (without a nominally assigned Poisson coefficient) is very good as can be seen in Table 2. In the case of the second gradient, the differences between the nominal and measured values are even greater due to the higher stiffness of the sample.
In Table 2, ν N is the nominal Poisson coefficient that was used in the constitutive continuum equation for FEM analysis, ν S is measured by FEM simulation and ν M is measured by our tool. From Table 2, we note that the trend of these values closes to 0.4 is different going towards its limit corresponding to the incompressible behaviour. Note as β =1.3 in the second gradient and 0 in the first gradient to approximate FEM solutions.

Limit of the model: an example
In this section, we intend to give some insights about the limitations of our model, particularly on the second gradient effects. For this purpose, we consider the same specimen as in the previous cases subjected to a shear test, with Poisson coefficient 0.4. We can see that the deformed external shape of the specimen, solved by means of the COMSOL system, is in very good agreement with our model in both the first and second gradient cases (Fig. 7). However, the internal strain distribution is quite different, especially in the case of the second gradient. In particular, our model shows excessive rotation of the vertical sections, compared with the solution found with the COMSOL system. Probably, the linearity of the rule (Eq. 4) fails to capture some aspects of the strain. In this paper, we stop here but these results deserve more attention in a future paper.

Conclusion
This paper deals with the possibility of describing a deformable body employing the rules governing the behaviour of a swarm of robots.
We have presented a simple algorithm easily adaptable to different problems to describe complex physical effects plausibly. This discrete model is purely kinematic, so only the positions of the particles describing the body appear. Deformation is applied via a displacement assigned to the leading particles, while the follower Fig. 7 Final configuration of the sample undergone a shear test. First gradient case (up) and second gradient case (down). Units on the axes are metres particles move according to "position rules" until a final equilibrium configuration is reached. These displacements are determined by the relative position of the particle's neighbours, as in a swarm of birds. Therefore, the deformed configuration is not computed from Newton's law, but only from the relative positions between the particles in the system, saving machine time for computation since this approach requires only the solution of an algebraic equation. The model is also easily extensible in the case of second gradient continua.
In this work, we compared the results obtained by solving a simple tensile problem by FEM and using our kinematic tool. Different values of Poisson's modulus were used and applied to a compact sample and another with a crack in the centre to compare the solutions at critical points. The results obtained are in good agreement, allowing us to hope that an improvement in the model, together with an understanding of the physics related to its parameters, may lead us to satisfactory quantitative predictions. Having demonstrated that we are able to obtain the same results, it is now necessary to find the physical meaning of the parameters of our model to make predictions of behaviour starting from defined initial conditions. Therefore, the results presented in this paper could be a promising theoretical tool for discrete approximation of continuous deformable bodies with generalized strain energy, encompassing non-Cauchy continua. In our opinion, what is attractive in the proposed model is its remarkable simplicity and flexibility in the choice of the type of interaction considered and the contact actions imposed. It allows us to describe non-local interaction and the use of second neighbours' particles can be interpreted as a second gradient theory, even though, in this case, the comparison results with the FEM approach are not as good as in the case of the first gradient model.
The proposed model can be enriched and generalized in many ways; more complex interactions can be considered, for example, cases where the interaction can be weighted by local density, distance, directionality, or other characteristics. A micromorphic approach is the natural generalization of the tool, where a material body is described as a continuous collection of deformable particles of finite size endowed with internal structure.
In our paper, we showed how flocking rules, used to determine the position in a robot swarm, can describe the deformation of a bidimensional continuum.
Our tool has many similarities to PBD, but unlike the PBD methods used in computer graphics, we do not require knowledge of velocity and we do not introduce any force to account for external actions. We remain in the framework of a fully kinematic system.
In this article only, static final configurations are shown. Here, we are not interested in evolution, but the model is able to show evolution in pseudo-time and specimen failure (see [18]).
We want to emphasize again that our model does not have a physical meaning yet, but we are working on a physical analogy. This is because we aim to create a complex swarm-based algorithm that can describe material deformations and not obtain another algorithm based on solving physical equations. Clearly, to assay the goodness of the model, we have to compare it with already available proven models and keep in mind that it is still developing.
The results presented are interesting but are still at a preliminary stage; however, they show good agreement with the deformations calculated by traditional methods. It is clear that in order to have a predictive theory of behaviour, we should know the displacement of the leaders instead of imposing it a priori, as done in this paper. The most important future development is to understand how the constitutive parameters of the materials are related to the parameter choices we make in our tool so that we can connect with the usual methods of continuum mechanics. However, the tool has shown enough flexibility to suggest that many other behaviours can be described once the association with the constituent parameters is established. Generalization to 3D should present no difficulty, but there is a need for some optimization in the code to keep computation time in the order of seconds using a standard desktop PC.
The model, based on the behaviour of robotic systems introduced here, will need further investigation and generalization in both theoretical and numerical aspects. Still, it promises to be of great interest.
Funding Open access funding provided by Ente per le Nuove Tecnologie, l'Energia e l'Ambiente within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.