Overlap Computation
4.3 Face and contact area determination
4.3.1 Face determination
FACE AND CONTACT AREA DETERMINATION 85
vertex and Algorithm 4.3 together with Algorithm 4.4 for generated vertices), we obtain all the vertices of the overlap polyhedron, e.g. in Fig. 4.9.
P1
P2
Fig. 4.9 The vertices of the overlap polyhedron (left) obtained after the computation of inherited vertices (circles) and of generated vertices (stars) for the two intersecting polyhedra (right, also in Fig. 4.6). The vertices obtained are points scattering in space and we need to find their topological relations, namely the faces, to determine the overlap polyhedron.
86 OVERLAP COMPUTATION
V
(a) (b)
P1 P2
F P1
P2 V g
V g V i
F i V i
V i F g
V g
V i V g V g
g
Fig. 4.10 Example of an inherited face Fi and a generated face Fg: In (a), for the two polyhedra P1 and P2, since the three vertices (Vi) of the triangular face (filled by dark gray) of P1 are all inside P2, the face is an inherited face for the overlap polyhedron;
In (b), for P1 and P2, there is only one vertex (Vi) of P1 inside P2, thus it is not an inherited face. The gray triangle on the face F of P2, which consists of three generated vertices Vg, is a generated face of the overlap polyhedron by the face intersections of F with the faces where Viis located. Faces which consist of both generated verticesVg and inherited vertices Vi are also generated faces.
parts from all intersecting faces of P1 and P2 which become the generated faces of the overlap polyhedron.
With the list of the intersecting face pairs (contact_f ace_pair from Algorithm 4.3) from computing the generated vertices, what remains to be determined for those faces are the indices of the generated vertices which they are located on (since one face can not intersect only one face). Instead of finding the inherited faces directly by checking their vertices, we make use of the VERTEX_FACE_TABLE array (in which for each vertex all the faces it is located on are stored, as explained in section 3.1.1) for the inherited vertices.
For an inherited vertex, we check all its faces: if the face has already been registered as a face of Po, we register it on the face ofPo; if not, we register the face as a new entry in the face list of Po and register the inherited vertex for this newly registered face. By doing so, we can not only register the inherited vertices on those generated faces it may belong to, but also find the inherited faces. The algorithm to find the faces of the overlap polyhedron Po of P1 and P2 is summarized in Algorithm 4.5.
FACE AND CONTACT AREA DETERMINATION 87
Algorithm 4.5 Determine the face of the overlap polyhedronPo
1: {generated faces fromP1}
2: for allthe faces ofP1: Fi(P1) do
3: if Fi(P1) is intersecting withP2 then
4: register Fi(P1) as a face ofPo: Fk(Po)
5: find the entries of Fi(P1) incontact_f ace_pair
6: find the vertex indices for all the intersection points of Fi(P1) from the corre-sponding entries in intersect_point_pair_idx
7: register all the vertex indices for Fk(Po)
8: end if
9: end for
10: {inherited faces fromP1}
11: for allthe inherited vertices fromP1: Vini (P1) do
12: for all the facesFi inVERTEX_FACE_TABLEof Vini (P1) do
13: if Fi is in the list of faces ofPo (Fi =Fk(Po ) then
14: register theVini (P1) for the faceFk(Po )
15: else
16: registerFi as a new entry in the face list of Po
17: registerVini (P1) for the face entry corresponding toFi
18: end if
19: end for
20: end for
21: {Repeat the same process for intersecting and inherited faces fromP2}
Though the faces of the original polyhedraP1andP2are triangles, the generated faces which consist of two inherited verticesViand two generated verticesVgare not necessarily triangular, as can be seen in Fig. 4.10 (a). Since our formulae (and the corresponding subroutines in the DEM code) for computing the physical properties of a polyhedron are based on triangular faces (section 3.2), to obtain the volume and the center of mass of the overlap polyhedron, we need to triangulate those generated faces with more than three vertices. In this case, we devised two algorithms, one using the centroid (Fig. 4.11) the other using an edge (Fig. 4.12) to determine the relative orientations of the vertices and order them counterclockwise. For the method which uses the centroid, we first need to set up a reference system for determining the orientations of the vertices to be ordered:
we choose the origin as the centroid C(Cx, Cy, Cz) of a set of k vertices Vi(Vix, Viy, Viz) which is given as the algebraic average of the vertex-coordinates,
C=
∑k i=1
Vi/k; (4.20)
the w-axis as the normal of the generated face nf, which can either be found in the
88 OVERLAP COMPUTATION
(v) f
θ 2
θ 5 θ 6
C
V1 (V’1) V3 (V’5)
V5 (V’4)
V2 (V’2) V6 (V’3) V4 (V’6)
Vb (u) Va
n
Fig. 4.11 Order the vertices on a generated face with respect to the centroid: Vi is the entry in the vertex indices list of the generated face before ordering and V′i is the entry in the list after ordering; C is the average of the vertex coordinates; nf is the normal of the face; a unit vector from C to the first vertex V1 is selected as the base unit vector Vb; an auxiliary unit vectorVa is chosen asnf×Vb. The angle between theVb and the vector CVi is θi from Eq. (4.23). The vertices Vi are ordered according to the values of θi. For the ordered list V′i, we can obtained a triangulation with triangles (C,V′i,V′i+1) fori= 1, . . . ,5 and (C,V′6,V′1).
FACE_EQUATION array (section 3.1.1) or computed from the three non-collinear vertices (Eq. (3.2)). The u-axis is chosen as the unit base vector Vb from Cto the first vertex in the list V1
Vb= V1−C
||V1−C||; (4.21)
and an auxiliary vector Va as the v-axis from
Va =nf×Vb. (4.22)
The sin(θi)and cos(θi)for theVi(i= 2, . . . k) can be computed from vector inner products sin(θi) = (Vi−C)·Va,
cos(θi) = (Vi−C)·Vb.
Using the atan2(y, x) function from FORTRAN5, we can obtain the angle θi forVi with
5The atan2(y,x) function gives the angle in radians for a point (x,y) to the positive x-axis for all four quadrants. The results are positive for counter-clockwise angles when y>0 and negative for clock-wise
FACE AND CONTACT AREA DETERMINATION 89
respect to the unit base vectorVb as
θi = atan2(sin(θi),cos(θi)), (4.23a) θi = 2π+θi, if θi <0. (4.23b) Sort Vi according to θi for the vertex list of a generated plane leads to a list of vertices ordered counter-clock-wisely, as shown in Fig. 4.11.
Vb 6
θ 5 θ 3 θ 4
V1 (V’1) V3 (V’5)
V5 (V’4)
V2 (V’2) V6 (V’3) V4 (V’6)
θ
Fig. 4.12 Order the vertices for a generated face with respect to an edge: the first two entries in the vertex list for a generated face in our algorithm always one of its edges.
We use the edgeV1V2 as unit base vectorVb and compute the cos(θi) for the remaining vertices Vi. The larger the cos(θi) is, the closer the Vi are to the edge V1V2. A trian-gulation is obtained automatically for the ordered list with triangles (V′1), V′i,V′i+1) for i= 2, . . . ,5.
This algorithm is self-sufficient, which means no additional information is necessary except the coordinates of the vertices to be ordered. Thus it can be used also in other appli-cations as a general approach for ordering a set of disordered points on a plane. In contract, in our overlap computation, we obtain convenient information which can facilitates the ordering process, namely the information of intersection segments (the intersect_point_-pair_idxarray from Algorithm 4.4), which are the edges of the overlap polyhedron. When we register the vertex indices of the intersection points for a generated face, the first two vertices are always recorded the list as a pair, which means for a generated face, we know at least one edge from the first two vertices of its vertex indices list, as shown in Fig. 4.12.
Since the angles of all the other vertices to the edgeV1V2 cannot be larger thanπ, instead
angles when y<0.
90 OVERLAP COMPUTATION
of computing the angle it is sufficient to compute only the cosine using the norms with Vb = V2−V1
||V2−V1||, cos(θi) = (Vi−V1)·Vb
||Vi−V1|| ,
for Vi (i= 3, . . . k) and then sort the Vi according to the cos(θi). The result from this method can either be counter-clock-wise or clock-wise, since the normal of the face has not been taken into account. To triangulate the generated face with a ordered vertex list is then trivial: the two different ordering methods themselves have already served as two different triangulation methods as shown in Fig. 4.11 and Fig. 4.12, usingCorV1 and the two vertices which form the opposite side of the angles θi′ (corresponds to the vertex V′i in the ordered list) become the triangles for the triangulation. The latter method which makes use of the edge information is implemented in our DEM simulation.
Once the coordinates of the vertices and triangular faces are both determined, the volume and the center of mass of the overlap polyhedron can be computed with Eq. (3.13) and Eq. (3.18) based the polyhedron decomposition method introduced in section 3.2.