Direction of a rectangle, part 2

Numerical orientation of a polygon continues

Hi again! Last time, I wrote about how I was given some a very basic FEM mesh, and had to determine the magnetization direction of a permanent magnet from it. In simpler words, I had a bunch of triangles determined by their corner coordinates, I knew their union formed a rectangle, and I had to determine the long and short axes of it. The direction of a rectangle, you could say.

I concluded that the axes are the eigenvectors of the matrix $\hat{\mathbf{M}}$

$\hat{\mathbf{M}} = \begin{bmatrix} \int\limits_{\Omega} x^2 \mathrm{d}x\mathrm{d}y & \int\limits_{\Omega} xy \mathrm{d}x\mathrm{d}y \\ \int\limits_{\Omega} xy \mathrm{d}x\mathrm{d}y & \int\limits_{\Omega} y^2 \mathrm{d}x\mathrm{d}y \end{bmatrix}$ ,

where $\Omega$ denotes the area that the magnet occupies.

Matrix evaluation

Now, the remaining problem is how the evaluate the matrix entries. Luckily for me, it can be easily done with some basic finite element codes that I already have. You see, in finite element analysis, we approximate the solution $A(x,y)$ as a weighted sum

$A(x,y) \approx \sum\limits_{k=1}^N a_k \shapefunction_k(x,y)$ .

Here, $a_k$ are the unknown coefficients to solve – let’s us not concern ourselves with them now.

More important for today are the shape functions $\shapefunction_k(x,y)$ . They are some scalar-valued (in our case at least) functions of the position $(x,y)$ , each associated with a certain node $k$ . The main point about them today, however, is that I know how to manipulate them inside and out. That’s the side-effect of having a functioning FEM library.

You see, each shape function $\shapefunction_k$ gets a value of 1 at the point $(x,y)$ corresponding to the node $k$ , and a zero value at the remaining nodes. Thus, we can do the following stupid-looking substitution

$x = \sum\limits_{k=1}^N x_k \shapefunction_k(x,y)$

$y = \sum\limits_{k=1}^N y_k \shapefunction_k(x,y)$ .

In other words, we can express the $x$ – and $y$ -coordinates with the shape functions and the coordinates of the nodes $(x_k, y_k)$ .

That dummy-relation can then be utilized the calculate the integrals for the $\hat{\mathbf{M}}$ matrix. For example, we can write

$\int\limits_{\Omega} x^2 \mathrm{d}x\mathrm{d}y$

$= \int\limits_{\Omega} \left( \sum\limits_{k=1}^N x_k \shapefunction_k(x,y) \right)^2 \mathrm{d}x\mathrm{d}y$ .

This, in turn, can be written as

$= \mathbf{x}^\mathrm{T} \mathbf{M} \mathbf{x}$ ,

where $\mathbf{x}$ contains the x-coordinates of the nodes

$\mathbf{x}^\mathrm{T} = \begin{bmatrix} x_1 & x_2 & \ldots x_n \end{bmatrix}$ ,

and the newly-introduced matrix $\mathbf{M}$ has the entries

$\mathbf{M}_{rc} = \int\limits_{\Omega} \shapefunction_r(x,y) \shapefunction_c(x,y) \mathrm{d}x\mathrm{d}y$ .

This is the so-called mass matrix, one of the most important matrices used in FEM. In other words, given the mesh data – which we have – assembling it is no problem at all. With the mass matrix $\mathbf{M}$ , the other matrix $\hat{\mathbf{M}}$ used in the axis calculation can be compactly written as

$\hat{\mathbf{M}} = \begin{bmatrix} \mathbf{x} & \mathbf{y} \end{bmatrix}^\mathrm{T} \mathbf{M} \begin{bmatrix} \mathbf{x} & \mathbf{y} \end{bmatrix}$ ,

the vector $\mathbf{y}$ containing the y-coordinates of the nodes.

Centering

Bear in mind one thing, however. Last time I assumed that the permanent magnet – the domain $\Omega$ – is “centered” around zero in some way. An exact definition would be that the center of mass of the magnet is in the origin. In other words, we want

$x_0 = \frac{ \int\limits_{\Omega} x \mathrm{d}x\mathrm{d}y }{ \int\limits_{\Omega} 1 \mathrm{d}x\mathrm{d}y } = 0$ ,

and equivalently for $y_0$ . Using a similar reasoning as above, this can be written as

$x_o = \frac{ \mathbf{x}^\mathrm{T} \mathbf{f} }{ \mathbf{1}^\mathrm{T} \mathbf{f} }$ .

Here, $\mathbf{1}$ is a vector of all-ones, and $\mathbf{f}$ is the prototype finite element load vector

$\mathbf{f}_r = \int\limits_{\Omega} \shapefunction_r(x,y) \mathrm{d}x\mathrm{d}y$ .

Once the center of mass has been calculated, the centered coordinates can be obtained by substituting

$x \leftarrow x - x_0$ and

$y \leftarrow y - y_0$ .

This has been illustrated in the figure you already saw last time. Indeed, the green circle denotes the center of mass of the magnet.

Conclusion

So all in all, a problem stemming from FEM – the direction of a rectangle, or PM – could easily be solved with FEM codes. Kinda boring.

However, on a deeper level the problem itself has more diverse, more multi-disciplinary aspects than this.

Firstly, the problem of determining the orientation of a rectangle is nothing else than computational geometry. Secondly, finding the previously-mentioned vector-matrix-vector product maximum was pure linear algebra.

And finally – and this is the most interesting one for me – the problem could also be approached from a purely mechanical angle. The rectangle could be interpreted as a solid object, and its long and short axes as its principal axes. Indeed, this principal axis interpretation ties our weird $\hat{\mathbf{M}}$ matrix to the inertia matrix often used in mechanics.

Mechanical engineers – let’s say we’re even now, shall we?

UPDATE: Check out this post for another possible application example of the algorithm!

-Antti

Check out EMDtool - Electric Motor Design toolbox for Matlab.

Need help with electric motor design or design software? Let's get in touch - satisfaction guaranteed!

Direction of a rectangle, part 2

Antti Lehikoinen

Direction of a rectangle, part 2

Matrix evaluation

Centering

Conclusion

Leave a Reply Cancel reply