Find displacement of vertices of a prism (polyhedron)

I am an engineering student trying to write a simple visualization program for a rigid block model subject to earthquake loading (using discrete element method). Now my DEM code is capable of outputting the displacement for each of the vertices in each block (each block being a 12-sided prism). However, I would like to minimize the data that I write to the disc and was trying to see how many vertices per block are needed to completely recover the displacements at all vertices and calculate their position (thus finding the entire block position and orientation).

After thinking about this problem I realized that I will need at least three points. It is embarrassing that it took me a while to figure this out since it seems rather obvious and this is stated everywhere (you need three non-collinear points).

Now I am trying to understand how I go from the three-components of displacements at three vertices to calculate displacements at the remaining vertices. So if I imagine a cube, this means that I will need to calculate the displacement of the 5 remaining vertices and then move those points to the appropriate coordinate to get the final block position.

As you can see I am not the brightest bulb around, so I would appreciate a very low-level explanation with as much real examples so I can understand what is going on. I had taken a linear algebra class many years ago and I was pretty bad it even then. Hopefully I can start over again.

Regards, Al

PS: I asked this question on mathOverflow and it was recommended that this is a better forum for this type of question. Hopefully that is okay.

-
Your blocks are 12-sided prisms? What shapes are the cross sections? Are the displacements associated just with the block itself, or is there a unique displacement for each vertex? –  Victor Liu Jun 21 '12 at 5:29
Well most of the blocks are 12 sided. The cross-sections are irregular convex polygons. The displacements are associated with each vertex but since all blocks are rigid, they really describe the motion of the block itself. –  Cannibal Flea Jun 21 '12 at 7:02

I assume what you mean is that your model consists of many prisms, with each vertex of the prisms having a displacement vector, but these vertices are not shared among prisms, so we need only look at one prism at a time. Furthermore, I assume that the prisms can only undergo rigid body motions, meaning that their initial and final states are related by a Euclidean transformation.

Overview of rigid transformations

Formally, let's assume a definite order of all the vertices, so that we can arrange the displacement vectors into a matrix, each column of which is one displacement vector of a particular vertex. Call it $V$. $V$ has 3 rows and (typically) 12 columns. Let's also assume you can do something similar for the prism vertex position vectors, forming a matrix $P$ (we won't be interested in $P$ itself, but it's useful for the discussion. The initial position of the prism is specified by $P$, and the final position is $Q=P+V$. The fact that you only allow rigid transformations means that $P+V=RP+\mathbf{t}\mathbf{1}^T$, where $\mathbf{1}$ is the vector of all 1's (in this case typically of length 12), $R$ is a rotation matrix with determinant 1, and $\mathbf{t}$ is a translation vector. So you can see that really, the entire set of displacement vectors is specified by just $R$ and $\mathbf{t}$, which is 12 real numbers, with the implicit assumption that you already know $P$ the prism position.

Your problem basically boils down to computing a rigid transformation given only 3 points. See, for example this paper.

Using only 3 vertices

Suppose you only consider the first three columns of $P$ and $V$ (call them $P_3$ and $V_3$, correspondingly $Q_3$), which is like only looking at 3 of the vertices. Then we have $$Q_3=RP_3+\mathbf{t}\mathbf{1}^T$$ The translation is very annoying, so let's get rid of it first. The centroid of the vertices is unaffected by the rotation, so let's work relative to the centroids. Call the centroids $p_c=\frac{1}{3}P_3\mathbf{1}$ and $q_c=\frac{1}{3}Q_3\mathbf{1}$. Define new translated matrices $$Q_c = Q_3 - q_c\mathbf{1}^T$$ $$P_c = P_3 - p_c\mathbf{1}^T$$ We now have $Q_c=RP_c$. You can't just solve for $R = Q_cP_c^{-1}$ because the matrices $Q_c$ and $P_c$ will only have rank 2. You need to compute the SVD of $P_cQ_c^T=U\Sigma V^T$, and set $R = VU^T$. To calculate any $i$-th displacement vector, $$v_i = R(p_i - p_c) + q_c$$

More robust scheme

If your displacements are not entirely consistent, then the rotation you compute may not be accurate (it may not be an exactly orthogonal matrix). It is instead best to store the $R$ matrix and the displaced centroid $q_c$ (redefined below) rather than just $V_3$. You want to compute the "best" $R$ given all the displacements you have. Following the linked paper, similarly compute the $3\times n$ centroid-relative matrices $$Q_c = Q - q_c\mathbf{1}^T$$ $$P_c = P - p_c\mathbf{1}^T$$ for analogous definitions of the centroids $q_c$ and $p_c$. You once again have $Q_c=RP_c$, but now $Q_c$ and $P_c$ are wider than they are tall. The best $R$ in the least squares sense is $R=VU^T$ where $V$ and $U$ are derived from the singular value decomposition of $P_cQ_c^T=U\Sigma V^T$. If your simulations are self-consistent, $\Sigma$ should be the $3\times 3$ identity matrix, but due to rounding errors, it may be slightly off.

Example

Suppose you have three points $(2,-1,4)$, $(-3,0,1)$, and $(5,1,-3)$, with respective displacements $(-0.75,-0.3,-0.58)$, $(0.65,-0.24,-0.72)$, $(0.07,-0.92,0.38)$. We would then form the following matrices: $$P_3 = \begin{bmatrix} 2&-3&5\\-1&0&1\\4&1&-3\end{bmatrix}$$ $$Q_3 = \begin{bmatrix} -0.75& 0.65&0.07 \\ -0.3& -0.24 &-0.92\\-0.58&-0.72&0.38\end{bmatrix}$$ The centroids are $p_c=(4/3,0,2/3)$ and $q_c = (1.323,-0.487,0.36)$. The recentered matrices are $$P_c = \begin{bmatrix} 2/3 & -13/3 & 11/3 \\ -1 &0& 1\\10/3 & 1/3 & -11/3 \end{bmatrix}$$ $$Q_c = \begin{bmatrix} -0.073 & -3.673 & 3.747 \\ -0.813 &0.246& 0.566\\3.06 & -0.08 & -2.98 \end{bmatrix}$$ Then compute the SVD of $P_cQ_c^T$, letting $\{U,\Sigma,V\}=svd(P_cQ_c^T)$ (yes, this $V$ is a different $V$ now). Finally, $R = VU^T$ results in $$R = \begin{bmatrix} 0.990 & 0.075 & -0.122 \\ -0.071 &0.997& 0.036\\0.125 & -0.027 & 0.991 \end{bmatrix}$$ Note that $R$ is almost the identity matrix (diagonal elements are almost 1, off diagonal are small) since our displacements were small. Note also that $RP_c \approx Q_c$. It's not exactly equal because the displacements I chose are only approximately corresponding to a rigid transformation (I rounded the numbers, so this introduces some stretch and shear).

Now, if we had a new point $p_4=(3,3,3)$, it's displacement is $$v_4 = q_c + R(p_4-p_c) - p_4 = (-0.087,-0.531,-0.199)$$

-
I am a bit lost when you are talking about the centroids (p.c and q.c). You said P.3 would be a 3x3 matrix with the original position of each of the three points in each of the columns (x y z). Is p.c the location of the centroid? I'm having problem doing this for a simple example of a rectangular prism. If I have three points on the face of the prism, then p.c = 1/3 P.3 {1 1 1} won't give you the centroid of the prism! I think I am missing something here, I just don't know what! –  Cannibal Flea Jun 21 '12 at 13:57
Also, do you think it would be possible to provide a simple example with your answer, maybe looking at a rectangular prism. This way through the application the concept becomes easier. –  Cannibal Flea Jun 21 '12 at 14:48
The centroids $q_c$ and $p_c$ are simply the averages of the set of points you're considering. So if you're using just 3 points, it's going to be the average of those three points, not the centroid of the prism itself. –  Victor Liu Jun 22 '12 at 1:49
Victor, Thank you for the example and extra explanation. However, there is something I still don't get. In the first section you say that the displacements of the i-th point would be $v_i = R(p_i - p_c) + q_c$ . However in the example you add an extra term ($-p_i$) . I follow the example up to this point. I was trying to see if I get the displaced matrix $Q_c$ but I get $$RP_c = \begin{bmatrix} 0.176&-4.329&4.153\\-0.924&0.322&0.602\\3.416&-0.21&-3.206\end{bmatrix}$$ which is very different from $Q_c$. Could you help me out a bit more? –  Cannibal Flea Jun 24 '12 at 19:09