Visualization of Ellipsoids

Visualization of Ellipsoids

Bachelor’s Thesis in Physics

Presented by:

Holger G¨ otz May 14, 2018

Friedrich-Alexander-Universit¨ at Erlangen-N¨ urnberg

Supervisor: Prof. Dr. Klaus Mecke

Carried out at: Department of Physics, Ume˚ a University

Supervisor in Ume˚ a: Peter Olsson

Contents

Contents

1. Introduction 3

2. Information about the simulations 3

2.1. Shape of used particles . . . . 4

2.2. A simulation cycle . . . . 4

3. Quaternions 5 3.1. Definition and properties . . . . 5

3.2. Describing rotations . . . . 6

3.2.1. Examples . . . . 7

3.3. Quaternions to rotation matrix . . . . 8

4. Blender 10 4.1. The Python API . . . . 10

4.2. Animations . . . . 11

4.3. Colouring and texture . . . . 12

5. The visualization script 12 5.1. The basic procedure . . . . 12

5.2. The input file . . . . 13

5.3. Applying of one particle configuration . . . . 13

5.3.1. Setting the location . . . . 13

5.3.2. Contact point highlighting . . . . 14

5.4. The render process . . . . 15

5.5. How to use this script . . . . 15

A. The Code 17 A.1. Python files . . . . 17

A.2. Cpython files . . . . 35

References 38

1. Introduction

1. Introduction

The intention of this project is to create movies of particle simulations with Blender [1].

Where in this specific case the simulations are contact only simulations of frictionless ellipsoids within the context of the so called jamming transition.

Within this document we first give more information about the simulations that are the base for the project. After that we introduce quaternions that are used by the simulations to describe the rotations of the particles. And finally we write about Blender, the program we used, and the actual script we created.

2. Information about the simulations

In order to better understand why specific things are done, this section is intended to provide more information about the simulations that shall be visualized. Though we encourage looking at e.g. [2] for further information.

The simulations are molecular dynamics ones of frictionless spheres or ellipsoids with contact only interactions. One special thing about them is that they are carried out at zero temperature. This means that the system, if not disturbed, will get stuck in a steady state with balanced forces after some time. The purpose of the simulations is to investigate the so called jamming transition. The jamming transition is a transition between the solid and the liquid state where both are equally disordered.

To prevent getting a crystal the simulations are done with two different sizes of particles such that there are equally as many small as big ones. To further prevent the system of getting into a steady state and to introduce the jamming transition the simulation is disturbed via shearing. This means that the simulation box is tilted with the time dependent tilt being described via the tilt variable γ = t ˙γ. ˙γ here is the shear-rate that determines the change of γ in time. For the purpose of the simulations the tilt variable is split into an integer part (γ

₁

) and a floating point part (γ

₀

) where |γ

₀

| < 0.5 so that one then jumps directly from γ

₀

= 0.5 to γ

₀

= −0.5.

Lees-Edwards boundary conditions are used to introduce this tilt. This means that a

particle at the position (x, y, z) is identical to the particle at the position (x ± γL

_y

,

y ± L

_y

, z), where L

_y

is the size of the simulation box in the y direction. For the x and

z coordinates periodic boundary conditions apply. With this setup the shearing indeed

can go on indefinitely, as the simulation has to only account for the γ

0

part. For this

purpose the red marked part in figure 1 is shifted to the left side when γ

₀

jumps from

+0.5 to −0.5. With this in mind the visualization will use only γ

₀

, too.

2. Information about the simulations

Figure 1: Top view of the simulation box with a tilt of γ

₀

= +0.5 on the left and γ

₀

= −0.5 on the right side. The unsheared box is indicated with green dots.

Identical parts are marked with the same colour.

2.1. Shape of used particles

As the simulated particles are ellipsoids, they are constructed via their axes. These are calculated by using two parameters α and d.

d

_x

= dα

^2/3

, d

_y

= dα

^−1/3

, d

_z

= dα

^−1/3

,

Note that d

_x

d

_y

d

_z

= d

³

so that the volumes of the particles do not depend on α. To account for the two different sizes of the particles d is set to be d

_b

= 1.4 or d

_s

= 1.0 for the big and small ellipsoids respectively. For the description of the rotation we use quaternions. An example of one big and one small particle is shown in figure 2.

2.2. A simulation cycle

Before the simulations start the particles are placed in system with arbitrary positions and rotations. To initialize the starting positions and rotations of the particles for the actual simulation run the simulations are executed for a few iterations without recording anything. Alternatively the final configuration of a previously done simulation may be used.

In the actual simulation run then the positions and rotations are saved for each step so that one effectively obtains them as a function of time.

Typically the simulations are done with around 1024 particles and a shear-rate ˙γ between

10

⁻¹¹

and 10

⁻¹

. Mostly there are about 1000 to 10000 time steps with a length of about

0.05. Though to be able to visualize the simulations they are executed with a reduced

number of only 64 particles.

3. Quaternions

Figure 2: Two particles with α = 1.2 in two different sizes and colours.

3. Quaternions

As the simulations use quaternions to specify the rotations of the particles we here give a short introduction to them. This section is based on [3] and [4].

3.1. Definition and properties

Quaternions are an extension of the complex numbers by the additional imaginary units j and k. These fulfil the following properties

i

²

= j

²

= k

²

= ijk = −1.

With this definition the imaginary units are anti commutative so that e.g. ij = −ji.

A general quaternion can then be written as

q = c

₀

+ c

_i

i + c

_j

j + c

_k

k ≡ q

₀

+ q

₁

+ q

₂

+ q

₃

. Alternatively a quaternion may be represented in the form

q = q

0

+ ~ q = [~ q,q

0

],

where ~ q = c

_i

i + c

_j

j + c

_k

k is a vector within the associated three dimensional vector space spanned by the basis {i,j,k}.

For further descriptions we define two quaternions q = [~ q,q

₀

] and q

⁰

= [~ q

⁰

,q

⁰₀

]. The

following properties can then be found:

3. Quaternions

Multiplication: qq

⁰

:= [~ q × ~ q

⁰

+ q

₀

~ q

⁰

+ q

₀⁰

~ q,q

₀

q

₀⁰

− ~ q · ~ q

⁰

].

Addition: q + q

⁰

:= [~ q + ~ q

⁰

,q

₀

+ q

₀⁰

].

Conjugate: q

^∗

:= [−~ q,q

₀

], (qq

⁰

)

^∗

= q

^0∗

q

^∗

, (q

^∗

)

^∗

= q.

Norm: N (q) := qq

^∗

=

3

P

i=0

q

_i²

with N (qq

⁰

) = N (q)N (q

⁰

), N (q

^∗

) = N (q).

Inverse: q

⁻¹

= q

^∗

/N (q), if N (q) 6= 0.

Where × is the cross product and · is the scalar product between three dimensional vectors.

3.2. Describing rotations

Quaternions are a useful way to describe rotation as e.g. they only require four numbers to be saved whereas e.g. a rotation matrix requires 9 (or at least 6) numbers [5]. Also the axis and angle of a rotation can easily be obtained by or converted into a quaternion.

Specifically if we want to describe a rotation by the angle θ around an unit vector ~ u we define an unit quaternion q as

q = q

₀

+ ~ q = cos( θ

2 ) + sin( θ 2 )~ u.

Then the vector

~

w

⁰

= q ~ wq

^∗

is the vector obtained from the rotation of the vector ~ w by the angle θ around ~ u.

Proof

We define the quaternions v

₀

= [~ v

₀

,0] and v

₁

= [~ v

₁

,0] with ~ v

₀

,~ v

₁

∈ R

³

, k~ v

₀

k = k~v

₁

k = 1,

~ v

0

×~v

1

6= 0. We further define q = v

1

v

^∗₀

= [~ v

0

×~v

1

,~ v

0

·~v

1

] = [~ v cos

^θ₂

, sin

^θ₂

] with ~ v =

_k~^~v_v^0×~v1

0×~v1k

and θ = 2∠(~v

0

,~ v

₁

), where ∠(~v

0

,~ v

₁

) is the angle between the vectors ~ v

₀

and ~ v

₁

. We also remember from the definition in Sec. 3.1, that the multiplication with an unit quaternion does not change the norm.

Now let’s look at v

⁰₀

v

1

:= (qv

0

q

^∗

)v

^∗₁

,

(qv

₀

q

^∗

)v

₁^∗

= (qv

₀

(v

₁

v

₀^∗

)

^∗

)v

₁^∗

= q v

₀

|{z}

=−v^∗₀

v

₀

v

^∗₁

|{z}

=−v1

v

₁^∗

= q(−1)(−1)

= v

₁

v

₀^∗

,

3. Quaternions

which means that

~ v

₁

· ~v

₀

= k~ v

₁

kk~v

₀

k · cos θ

2 = ~ v

⁰₀

· ~v

₁

, and

~

v

₁

× ~v

⁰₀

= −~ v

₁

× ~v

₀

.

So that ~ v

⁰₀

,~ v

₀

and ~ v

₁

are in the same plane and ∠(~v

0

,~ v

⁰₀

) = θ. So ~ v

⁰₀

is obtained by rotating ~ v

₀

by the angle θ around ~ v.

Now we look at v

₁⁰

v

₀^0∗

:= (qv

₀

q

^∗

)v

₁^∗

v

^0∗₀

,

(qv

1

q

^∗

)v

₀^0∗

= (q(qv

0

)q

^∗

(qv

0

q

^∗

)

^∗

)

= q = v

₁

v

₀^∗

,

which means that

~

v

⁰₁

· ~v

⁰₀

= ~ v

₁

· ~v

₀

, and

~ v

⁰₀

× ~v

⁰₁

= ~ v

₀

× ~v

₁

.

So that ~ v

0

,~ v

⁰₀

,~ v

1

,~ v

⁰₁

are in the same plane and ∠(~v

⁰

,~ v

1

) = ∠(~v

⁰0

,~ v

⁰₁

). With that we can conclude, that again ~ v

⁰₁

is obtained from a rotation of ~ v

₁

by θ around ~ v. Now we have covered every angle that is perpendicular to ~ v. So we only need to look at v

⁰

:= qvq

^∗

,

qvq

^∗

= [ ~ q × ~ v

| {z }

=cos^θ₂~v×~v=~v×~q

,~ q · ~ v]q

^∗

= [~ v × ~ q + ~ v · ~ q]q

^∗

= vqq

^∗

= v.

Now we have shown, that every vector ~ u

⁰

defined with [~ u

⁰

,0] := q[~ u,0]q

^∗

is ~ u rotated by the angle θ around ~ v.

3.2.1. Examples Trivial Quaternions

The quaternion q = i = ~ e

₁

then can be used to apply a rotation of π around ~ e

₁

. The

quaternions j and k can be interpreted similarly. Whereas q = 1 or q = −1 would

represent a rotation of the angle 0, 2π or rather no rotation at all.

3. Quaternions

(a) θ = 0 (b) θ = π/4 (c) θ = π/2

Figure 3: A cuboid with different rotations around the x-axis (displayed with a red arrow).

Rotation around the x axis

For a rotation of the angle θ around the x-axis we accordingly define the quaternion to be

q = [~ e

_x

sin( θ

2 ), cos( θ

2 )]. (1)

E.g. for θ =

^π₂

this would be [~ e

_x

√ 2 2

,

√ 2 2

] =

√ 2

2

[~ e

_x

,1] and for θ =

^π₄

this would be

≈ [~e

_x

0.38,0.92]. For further illustration these rotations are applied to a cuboid shown in figure 3.

3.3. Quaternions to rotation matrix

For the purpose of creating a rotation matrix we look at the quaternion as if it were a 4-vector

[~ q,q

₀

] ≡







 q

_x

q

y

q

_z

q

₀









.

3. Quaternions

We first look at pv with v = [~ v,v

₀

],

qv =[~ q × ~ v + q

₀

~ v + v

₀

~ q,q

₀

v

₀

− ~ q · ~ v]

=









q

_y

v

_z

− q

_z

v

_y

+ q

₀

v

_x

+ q

_x

v

₀

q

_z

v

_x

− q

_x

v

_z

+ q

₀

v

_y

+ q

_y

v

₀

q

x

v

y

− q

y

v

x

+ q

0

v

z

+ q

z

v

0

q

₀

v

₀

− q

_x

v

_x

− q

_y

v

_y

− q

_z

v

_z









=









q

₀

−q

_z

q

_y

q

_x

q

_z

q

₀

−q

_x

q

_y

−q

y

q

x

q

0

q

z

−q

_x

−q

_y

−q

_z

q

₀















 v

_x

v

_y

v

z

v

₀









=:R

_q

v.

So we have the rotation matrix for the left multiplication. Now we have to look at vp

^∗

(remember: q

^∗

= [−~ q,q

₀

]),

vq

^∗

=[~ v × ~ q

^∗

+ ~ vq

₀^∗

+ v

₀

~ q

^∗

,v

₀

q

₀^∗

− ~v · ~ q

^∗

]

=









q

y

v

z

− q

z

v

y

+ q

0

v

x

− q

x

v

0

q

_z

v

_x

− q

_x

v

_z

+ q

₀

v

_y

− q

_y

v

₀

q

_x

v

_y

− q

_y

v

_x

+ q

₀

v

_z

− q

_z

v

₀

q

₀

v

₀

+ q

_x

v

_x

+ q

_y

v

_y

+ q

_z

v

_z









=









q

0

−q

z

q

y

−q

x

q

_z

q

₀

−q

_x

−q

_y

−q

_y

q

_x

q

₀

−q

_z

q

_x

q

_y

q

_z

q

₀















 v

x

v

_y

v

_z

v

₀









=:R

q^∗

v.

4. Blender

To get the total rotation matrix for the rotation represented by q we have to multiply R

_q

and R

_q^∗

.

R

_q

R

_q^∗

=









q

₀

−q

_z

q

_y

q

_x

q

_z

q

₀

−q

_x

q

_y

−q

_y

q

_x

q

₀

q

_z

−q

_x

−q

_y

−q

_z

q

₀

















q

₀

−q

_z

q

_y

−q

_x

q

_z

q

₀

−q

_x

−q

_y

−q

_y

q

_x

q

₀

−q

_z

q

_x

q

_y

q

_z

q

₀









=









q

²₀

− q

_z²

− q

²_y

+ q

²_x

−2q

₀

q

_z

+ 2q

_x

q

_y

2q

₀

q

_y

+ 2q

_x

q

_z

0 2q

₀

q

_z

+ 2q

_x

q

_y

−q

_z²

+ q

₀²

− q

²_x

+ q

_y²

2q

_y

q

_z

− 2q

₀

q

_x

0

−2q

₀

q

_y

+ 2q

_x

q

_z

2q

_y

q

_z

+ 2q

_x

q

₀

−q

²_y

− q

_x²

+ q

₀²

+ q

_z²

0

0 0 0 q

_x²

+ q

_y²

+ q

_z²

+ q

²₀









=









N (q) − 2(q

_y²

+ q

_z²

) 2(q

_x

q

_y

− q

₀

q

_z

) 2(q

_x

q

_z

+ q

₀

q

_y

) 0 2(q

_x

q

_y

+ q

₀

q

_z

) N (q) − 2(q

_x²

+ q

_z²

) 2(q

_y

q

_z

− q

₀

q

_z

) 0 2(q

_x

q

_z

− q

₀

q

_y

) 2(q

_y

q

_z

+ q

₀

q

_x

) N (q) − 2(q

_x²

+ q

²_y

) 0

0 0 0 N (q)









We now have the total rotation matrix R = R

q

R

q^∗

with Rv = qvq

^∗

. For three dimen- sional rotation of a vector we can of course only use the 3 × 3 part of R, as v

₀

= 0 for a vector.

4. Blender

Blender is a 3D creation suit. This means that it provides a wide range of possibilities to create images, animations or movies or even the graphics for a video game. As such it includes a variety of tools and methods for modeling, rendering, animation, video editing and even more [1]. It also supports scripting and extending via Python and is released under GNU General Public License (GPL). It can run on multiple platforms such as Linux, macOS and Windows.

4.1. The Python API

For us the most important part is Blender’s Python API that may be used for every action available in Blender. Python scripts in Blender can be used as a plugin to extend the functionality of Blender as well as to represent libraries for other Python scripts, to set various Blender preferences or even to be executed on startup to define some GUI elements or operators. Additionally there is the possibility to execute the script once when it is needed. As we are not really interested in extending Blender but in creating movies with it we chose to make a script that can be executed as needed. For that purpose we have three possibilities to execute them:

The Python console: With the console you get a very interactive and direct access to

the object’s data. This makes it a good way to get to know the Python API or test

4. Blender

single pieces of code. For the execution of larger code-snippets we do not think of it as a very good possibility.

The text editor area: The editor offers the possibility to open, edit and execute scripts directly from the Blender GUI. Especially the possibility to see the result of a script in the 3D-view without the need of rendering anything is very useful in some cases. Though it might be worth mentioning that the Blender GUI is frozen during the execution of the script.

The command line interface: This interface gives a very convenient way to execute a Python script via Blender without the need to open any graphical interface. This has the big advantage, that it can be executed without any monitor present. Also the render process might be faster as none of the rendered frames has to be shown directly. Additionally it is possible to specify any options via custom command-line options rather than editing the script itself.

The script discussed in this document is designed to work with the command line inter- face, but it still is possible to execute it from inside Blender via the Text Editor area.

A Python documentation of all the available methods and functions can be accessed through the Blender website [1].

4.2. Animations

As we want to create movies rather than pictures of the simulations we have to make animations. For this purpose we use so called keyframes. A keyframe is defined as a marker which stores values of specific properties for a specific point in time. For us these properties will mainly be the location and rotation of the particles. Keyframes also allow automatic interpolations to be done by Blender if they are not defined in two consecutive frames. Interpolation here means that Blender creates estimated states of properties specified in the keyframes for the frames in between so that a change might appear smoother. If the use of these interpolations is not desired, just as in our case, Blender has be told not to do them. Though we normally define a keyframe every frame, so no interpolation would occur anyway.

For our purpose the data that is contained in the keyframes is limited to the location and quaternion of a particle, as the further properties do not change. Only for the objects that are used to highlight the contact points we also insert keyframes for the size and visibility.

In the latter case we use keyframes for the visibility in order not to have to remove and

recreate the objects used for highlighting. This increases the speed of the movie creation

process significantly as the creation of Blender objects takes considerably longer than

defining keyframes.

5. The visualization script

4.3. Colouring and texture

We create the particles in two colours, red and blue, in order to make it easier to observe a specific particle. This is achieved by the use of materials. Blender materials can represent various properties of an object like shininess, reflectance or transparency.

We however use it for the colouring purpose only.

Additionally, we add a cloud texture in grey to the materials in order to make it easier to follow the rotations of the ellipsoids as this can be quite hard if there is no pattern on the surface. The standard cloud texture we use is created by Blender using Perlin noise which is a gradient noise named after the developer Ken Perlin. For illustration purpose three particles, one without material and texture, one with a material only and one with both, are displayed in figure 4.

Figure 4: Three particles, one with no material (left), one with a red material (middle) and one with a material and a cloud texture (right).

5. The visualization script

The visualization script is intended to be used for the visualization of moving particles (as described in section 2.1). As these simulations are done with contact only interactions the contact points between the particles might be of interest and are therefore also made visible.

In order to function the script has to be executed via Blender or somehow else be supplied with the necessary libraries. The script was tested with Blender version 2.69.

5.1. The basic procedure

After the input file is read the script creates the visible objects in Blender that will

represent the particles. A keyframe is then inserted for the rotations and locations of

the particles for every frame. Also the contact points are determined and highlighted

5. The visualization script

via the use of additional graphical objects. When everything is specified the frames are rendered and saved by Blender.

5.2. The input file

The input file should look like the configuration file defined in the instructions. This means it should contain the following for every time step, and thereby frame, to be displayed:

• N (int): the number of particles.

• n

_var

(int): the number of variables used to describe the configuration (in this specific case that would be 13N + 2).

• L

_x

, L

_y

, L

_z

(3 double): the size of the simulation box.

• x, y, z, v

_x

, v

_y

, v

_z

, ω

_a

, ω

_b

, ω

_c

, q

₀

, q

₁

, q

₂

, q

₃

(13N double): the locations, velocities, rotational velocities and the quaternion describing one particle. Please note that the angle part of the quaternion should be specified as q

₃

, whereas q

₀

, q

₁

, q

₂

are the parts for the x, y and z axis respectively. Also note that the values of the (rotational) velocities are not relevant for the visualization but still have to be supplied.

• γ

0

, γ

1

(2 double): The tilt-variable γ split up in two parts, where γ

1

is the integer part and |γ

₀

| ≤ 0.5.

Every even particle is small and every odd one is big (starting numbering with 0).

5.3. Applying of one particle configuration

After reading the configuration file the animation data has to be applied. Specifically this means that for every frame the locations of the particles have to be set and the contact points have to be highlighted. So in this section we describe how this is done for one frame.

5.3.1. Setting the location

The setting of the locations may be very easy. If the particles shall be displayed as

specified in the input file the properties can simply be transferred to the displayed

objects. But as the simulations may specify the particles to be outside of the cuboid-

shaped (not tilted) simulation box due to the shearing, it might still be wanted to display

the particles inside the box. In this case some further calculations have to be done in

order to relocate some ellipsoids back into the simulation box. Additionally this then

prevents jumps in the created movies that are introduced when the tilt variable suddenly

5. The visualization script

Figure 5: Two consecutive frames with a jump from γ

₀

= 0.5 to γ

₀

= −0.5

changes from γ

0

= 0.5 to γ

0

= −0.5. Figure 5 shows two consecutive frames where the described jump happens.

To implement this relocation we check if a particle is actually outside of the simulation box. If this is the case it should be relocated by ±L

_x,y,z

, where L

_x,y,z

is the size of the simulation box in the respective direction. As this may lead to unwanted back and forth jumping of a particle in two or more consecutive frames we introduce a small tolerance that a particle may deviate from the simulation box. To further also detect such jumps if they are introduced by the input file we observe the position of a particle in relation to the origin.

So before we relocate a particle we check if x,y and z, the coordinates of the particle, are greater than zero in the current and in last frame. If then e.g. x > 0 in both frames or x < 0 in both frames the particle is moved by ±L

_x

along the x axis if |x| > 0.5 · L

_x

+ τ

_x

, with τ

x

the tolerance for the x-axis. The sign (plus or minus) and thereby the direction of the shift along the axis is determined by the sign of the coordinate. This means that if e.g. |x| > 0.5 · L

_x

+ τ

_x

and x > 0 the particle is shifted by −L

_x

along the x-axis. If x < 0 the particle would be shifted by +L

_x

along the x-axis.

But if x > 0 in one frame and x < 0 in the other frame the particle is relocated by ±L

x

if |x| > 0.5 · L

_x

− τ

_x

along the x-axis. Again the direction is determined by the sign of the coordinate. Of course x is only an example here and the same is also done for the y and z coordinate of the particle.

If none of the described options apply the particle will not be relocated. In any case we then insert a keyframe for the object with the new location and quaternion so that the visible object is also changing its position.

5.3.2. Contact point highlighting

After the locations are set for one frame, the contact points are highlighted. Obviously the first thing that has to be done is to check if and where two objects are having contact.

For this purpose we use a provided C-function that is wrapped with Cython and returns

a struct (or rather a Python dictionary) with information about the contact point if it

exists. Cython [6] is a compiler that makes it possible to write C-extensions for Python.

5. The visualization script

So then if two ellipsoids are in contact we have to check if the returned points are actually valid for the visible objects. This is done in order not to highlight points where the contacting objects are not displayed next to each other. The check is necessary because the C-function that returns the contact point between two ellipsoids is aware of the boundary conditions. If indeed the visible objects have contact as well we set a red cylinder to appear at the intersection points so that its top and bottom surface is tangent to the surface of the ellipsoids. The scale of the cylinders is calculated, as if there were two spheres of radius r

_max

touching. Where r

_max

is the maximum of the semi-axes of the touching ellipsoids. Particularly the scale is then calculated with:

1.1 · p

r

_max²

− d

²

where d is the distance between the intersecting surfaces calculated via the provided contact points on the surfaces of each of the cylinders. This may seem as a rather arbitrary choice but it yields the desired appearance.

It might be also noteworthy how we handle the cylinders that highlight the intersections.

Every cylinder created once can and will be reused for any other contact point that may occur. For this purpose we save the object of every created cylinder in a list. For each frame then the contacting points are calculated and highlighted with an element from that list. The properties such as the location, rotation and even the size of this object will then be adjusted via keyframes. If the list does not contain enough objects to mark all the points new ones are created as needed. If on the other hand there are too many cylinders in the list the rest will be set to invisible for that frame, so that they don’t have to be removed and can be reused in later frames. We reuse the cylinders as creating and deleting new or old cylinders would make no sense and would take considerably longer.

5.4. The render process

Before now finally the render process starts, the camera is set to view top down onto the particles and focus them all. Also the lamp is set to light from the same direction as the camera. Shadows are turned off, as this leads to a better view of the particles that are not in the top plane.

When all is set the render process is finally started and Blender creates a movie that is saved to the desired location. For illustration of the result, figure 6 shows one frame of a created movie.

5.5. How to use this script

This script is intended to be used via the command-line-interface of Blender. A basic command would look like

bash$ blender --background --python ./script/render_particles.py \

-- Lat/0064_r7000_GDf100_KDk500_Ml100_EL120_idt20

5. The visualization script

Figure 6: Sample frame created by the described script.

where everything after ” - - ” is passed as a command line argument to the Python script. The following (and more) command line parameters can be used:

-o ofilename With this option an output file can be specified. Please note, that Blender might add some numbers describing the rendered frames. If -o is not defined the script uses the input filename as a base.

-of ofolder With this option it is possible to define an output folder to save the generated movie in, if -o is not provided. If neither -of nor -o are specified the script will use the directory of the input file.

-nf number Tells the script to apply only every number th frame (default is 1).

-f Specifies that the particles should be relocated in order to always be displayed in the simulation box.

-rp number Sets the resolution to 0.01·number ·(1920,1080) (default is 100).

–end number Tells the script to stop after reaching frame number (default is 1000).

-h, –help Shows a help message and displays all available options.

If the script is executed directly e.g. via ./render particles.py the script will only read the configurations and then exit. Also a command how to use the script via Blender is printed out.

A second possibility is to execute the script from inside Blender via the Text Editor

area. With this method however it is required to directly specify the command line

A. The Code

arguments within the script. Also Blender has to be started from the directory with the custom libraries, as otherwise there might be some import problems. But as this is not the intended usage, no further explanation is provided at this point.

A. The Code

A.1. Python files

The file render particle.py, displayed in listing 1, is the main script. It imports methods from utilities.py (shown in listing 2) as well as from the library compiled with Cython (see Appendix A.2) and uses them to create the movie. A short description of what each function does is included as comments in the listnings.

Listing 1: render particles.py

1 # !/ usr / bin / env p y t h o n 3 2

3 i m p o r t sys , os 4 i m p o r t m a t h 5 i m p o r t r a n d o m 6 i m p o r t a r r a y

7 f r o m s t a t i s t i c s i m p o r t m e a n 8

9 # Add the c u r r e n t d i r e c t o r y to the python - p a t h if it is not a l r e a d y 10 # i n c l u d e d

11 if os . p a t h . d i r n a m e ( os . p a t h . a b s p a t h ( _ _ f i l e _ _ )) not in sys . p a t h : 12 sys . p a t h . i n s e r t (0 , os . p a t h . d i r n a m e ( os . p a t h . a b s p a t h ( _ _ f i l e _ _ ))) 13 p r i n t(" A d d e d file - d i r e c t o r y to p y t h o n p a t h in o r d e r to " \ 14 " f i n d the c u s t o m l i b r a r y ")

15

16 # Try s o m e s t e p s to i m p o r t f u n c t i o n s f r o m the c u s t o m C - l i b r a r y 17 try:

18 f r o m m o d u l e i m p o r t check , c h e c k _ g e t _ o v e r l a p 19 e x c e p t:

20 sys . p a t h . i n s e r t (0 , os . p a t h . a b s p a t h (" . "))

21 try:

22 f r o m m o d u l e i m p o r t check , c h e c k _ g e t _ o v e r l a p 23 e x c e p t:

24 r a i s e I m p o r t E r r o r (" C o u l d not i m p o r t c u s t o m l i b r a r y . \ n " \

25 " P l e a s e c o n s i d e r s t a r t i n g B l e n d e r f r o m " \

26 " the d i r e c t o r y t h a t c o n t a i n s it ")

27

28 # i m p o r t f u n c t i o n s f r o m u t i l i t i e s . py

29 f r o m u t i l i t i e s i m p o r t r e a d _ c o n f , g e t _ a l p h a , g e t _ s e m i _ a x e s 30

31 # S p e c i f y s o m e g l o b a l v a r i a b l e s .

32 # S o m e of t h e m m i g h t be o v e r w r i t t e n l a t e r . 33 m a i n = 0

34 V E R B O S E = 0 35 T O L E R A N C E = N o n e

A. The Code

36 r e n d e r = 1 37

38 # Try to i m p o r t l i b r a r i e s p r o v i d e d by B l e n d e r 39 try:

40 i m p o r t bpy

41 f r o m m a t h u t i l s i m p o r t Vector , Q u a t e r n i o n 42

43 m a i n = 0

44 a r g s = sys . a r g v [ sys . a r g v . i n d e x (" - - ") + 1:]

45 e x c e p t:

46 m a i n = 1

47 a r g s = sys . a r g v [ 1 : ]

48 p r i n t(" E R R O R : Not p o s s i b l e to i m p o r t the m o d u l e \" bpy \ " . \ n "\ 49 " In o r d e r to r e n d e r p l e a s e e x e c u t e \" b l e n d e r " \

50 " - - b a c k g r o u n d - - p y t h o n {} - - { } \ " . \ n ".f o r m a t( \

51 sys . a r g v [0] , " ". j o i n ( sys . a r g v [1:])) , f i l e= sys . s t d e r r ) 52

53

54 # # P a r s e A r g u m e n t s 55 if len( a r g s ) <1:

56 # D e f i n e A r g u m e n t s h e r e in o r d e r to 57 # e x e c u t e d via the T e x t E d i t o r A r e a 58 p r i n t(" not m a i n ")

59 a r g s = [ ’ ../ Lat / 0 0 6 4 _ r 7 0 0 0 _ G D f 1 0 0 _ K D k 5 0 0 _ M l 1 0 0 _ E L 1 2 0 _ i d t 2 0 ’, \

60 ’ - r ’, ’ 0 ’, ’ - - s t a r t ’, ’ 420 ’, ’ - - end ’, ’ 920 ’, ’ - f ’]

61

62 m a i n = 0

63 r e n d e r = 0

64

65 # D e f i n e the a r g u m e n t s 66 i m p o r t a r g p a r s e

67 p a r s e r = a r g p a r s e . A r g u m e n t P a r s e r ( d e s c r i p t i o n =’ C r e a t e a n i m a t i o n s . ’) 68

69 p a r s e r . a d d _ a r g u m e n t (’ f i l e 1 ’, m e t a v a r =’ f i l e n a m e ’, t y p e=str, \

70 h e l p=’ the i n p u t f i l e n a m e ’)

71 p a r s e r . a d d _ a r g u m e n t (’ - o ’, d e s t =’ o f i l e ’, t y p e=str, \

72 h e l p=’ o u t p u t f i l e n a m e ’)

73 p a r s e r . a d d _ a r g u m e n t (’ - of ’, d e s t =’ o f o l d e r ’, t y p e=str,\

74 h e l p=’ o u t p u t f o l d e r - N o t e t h a t t h i s w i l l be ’ \

75 ’ n e g l e c t e d if \ ’ - o \ ’ is u s e d ’)

76 p a r s e r . a d d _ a r g u m e n t (’ - nf ’, d e s t =’ n f r a m e ’, \

77 t y p e=int, d e f a u l t =1 , \

78 h e l p=’ use o n l y any \ ’ n f r a m e \ ’ th ’ \

79 ’ p a r t i c l e - c o n f i g u r a t i o n ’)

80 p a r s e r . a d d _ a r g u m e n t (’ - vf ’, ’ - - video - f r a m e ’, d e s t =’ v f r a m e ’, \

81 t y p e=int, d e f a u l t =1 , \

82 h e l p=’ a p p l y one p a r t i c l e c o n f i g u r a t i o n e v e r y ’ \

83 ’ \ ’ v f r a m e \ ’ video - f r a m e . O b v i o u s l y t h i s ’ \

84 ’ may l e a d to non - s m o o t h m o v e m e n t ’)

85 p a r s e r . a d d _ a r g u m e n t (’ - f ’, ’ - - f o r c e ’, a c t i o n =’ s t o r e _ t r u e ’, \

86 h e l p=’ f o r c e o b j e c t s to be in the box ’)

87 p a r s e r . a d d _ a r g u m e n t (’ - r ’, ’ - - r e n d e r ’, d e s t =’ r e n d e r ’, \

A. The Code

88 t y p e=int, d e f a u l t =1 , \

89 h e l p=’ d i s a b l e r e n d e r w i t h \ ’ - r 0\ ’ ’)

90 p a r s e r . a d d _ a r g u m e n t (’ - rp ’, ’ - - r e s o l u t i o n - p e r c e n t a g e ’,\

91 d e s t =’ r e s o l u t i o n _ p e r c e n t a g e ’, \

92 t y p e=int, d e f a u l t =100 , \

93 h e l p=’ set r e s o l u t i o n p e r c e n t a g e , at 100 ’ \

94 ’ the r e s o l u z i o n is 1 9 2 0 , 1 0 8 0 ’)

95

96 p a r s e r . a d d _ a r g u m e n t (’ - - s t a r t ’, d e s t =’ s t a r t _ f r a m e ’, \

97 t y p e=int, d e f a u l t =0 , \

98 h e l p=’ s t a r t w i t h the \ ’ s t a r t _ f r a m e \ ’ th f r a m e ’) 99 p a r s e r . a d d _ a r g u m e n t (’ - - end ’, d e s t =’ e n d _ f r a m e ’, \

100 t y p e=int, d e f a u l t = -1 , \

101 h e l p=’ end w i t h the \ ’ e n d _ f r a m e \ ’ th f r a m e ’)

102 p a r s e r . a d d _ a r g u m e n t (’ - - measure - t i m e ’, a c t i o n =’ s t o r e _ t r u e ’, \

103 d e s t =’ m e a s u r e _ t i m e ’, \

104 h e l p=’ M e a s u r e s the t i m e ’)

105 p a r s e r . a d d _ a r g u m e n t (’ - - t r a n s p a r e n t ’, a c t i o n =’ s t o r e _ t r u e ’, \

106 d e s t =’ t r a n s p a r e n t ’, \

107 h e l p=’ Set the b a c k g r o u n d to t r a n s p a r e n t ’)

108 p a r s e r . a d d _ a r g u m e n t (’ - - i m a g e ’, a c t i o n =’ s t o r e _ t r u e ’, \

109 d e s t =’ i m a g e ’, \

110 h e l p=’ R e n d e r i m a g e i n s t e a d of movie , - - s t a r t ’ \

111 ’ and - - end s h o u l d be the s a m e for t h i s ’ \

112 ’ o p t i o n ’)

113

114 # P a r s e the a r g u m e n t s

115 a r g s = p a r s e r . p a r s e _ a r g s ( a r g s ) 116

117 if a r g s . m e a s u r e _ t i m e : 118 i m p o r t t i m e

119 s t a r t _ t i m e = t i m e . t i m e () 120

121 if a r g s . o f o l d e r and not a r g s . o f i l e :

122 a r g s . o f i l e = os . p a t h . j o i n ( a r g s . ofolder , \

123 os . p a t h . b a s e n a m e ( a r g s . f i l e 1 ))

124 e l s e:

125 a r g s . o f i l e = a r g s . o f i l e if a r g s . o f i l e e l s e a r g s . f i l e 1 126

127 a r g s . f i l e 1 = os . p a t h . a b s p a t h ( os . p a t h . e x p a n d u s e r ( a r g s . f i l e 1 )) 128 a r g s . o f i l e = os . p a t h . a b s p a t h ( os . p a t h . e x p a n d u s e r ( a r g s . o f i l e )) 129

130 def g e t _ s u n ():

131 """

132 T h i s m e t h o d r e t u r n s a L a m p o b j e c t if f o u n d 133 """

134 if(len( bpy . d a t a . c a m e r a s ) == 1):

135 r e t u r n bpy . d a t a . o b j e c t s [’ L a m p ’]

136 r a i s e L o o k u p E r r o r (" L a m p not f o u n d ") 137

138 def g e t _ c a m e r a ():

139 """

A. The Code

140 T h i s m e t h o d r e t u r n s the c a m e r a o b j e c t if f o u n d 141 """

142 if(len( bpy . d a t a . c a m e r a s ) == 1):

143 r e t u r n bpy . d a t a . o b j e c t s [’ C a m e r a ’]

144 r a i s e L o o k u p E r r o r (" C a m e r a not f o u n d ") 145

146 def s e t _ c a m e r a _ p o s ( pos =" n o r m a l ", l i m i t s =(0 ,0 ,0) , a d j u s t _ s u n =T r u e):

147 """

148 T h i s m e t h o d can be u s e d to set the p o s i t i o n of the camera ,

149 as w e l l as the l i g h t s o u r c e

150 """

151 c a m e r a = g e t _ c a m e r a () 152 if pos ==" n o r m a l ":

153 p a s s

154 e l i f pos == " n o r t h ": 155 c a m e r a . l o c a t i o n . x = 0 156 c a m e r a . l o c a t i o n . y = 0

157 c a m e r a . l o c a t i o n . z = 4* l i m i t s [2]

158

159 c a m e r a . r o t a t i o n _ e u l e r . x = 0 160 c a m e r a . r o t a t i o n _ e u l e r . y = 0 161 c a m e r a . r o t a t i o n _ e u l e r . z = 0 162 if a d j u s t _ s u n :

163 sun = g e t _ s u n ()

164 sun . l o c a t i o n . x = c a m e r a . l o c a t i o n . x 165 sun . l o c a t i o n . y = c a m e r a . l o c a t i o n . y 166 sun . l o c a t i o n . z = c a m e r a . l o c a t i o n . z 167 f o c u s _ a l l ()

168

169 def f o c u s _ a l l ():

170 """

171 S e l e c t s e v e r y o b j e c t in o r d e r to f o c u s the c a m e r a on t h e m . 172 """

173 # S e l e c t o b j e c t s t h a t w i l l be r e n d e r e d

174 bpy . ops .o b j e c t. s e l e c t _ a l l ( a c t i o n =’ D E S E L E C T ’) 175 for obj in bpy . c o n t e x t . v i s i b l e _ o b j e c t s : 176 if not ( obj . h i d e or obj . h i d e _ r e n d e r ):

177 obj . s e l e c t = T r u e

178

179 bpy . ops . v i e w 3 d . c a m e r a _ t o _ v i e w _ s e l e c t e d () 180 bpy . ops .o b j e c t. s e l e c t _ a l l ( a c t i o n =’ D E S E L E C T ’) 181

182 def r e n d e r _ a n d _ s a v e ( fname , f r a m e _ e n d , f r a m e _ s t a r t =0 , \

183 r e s o l u t i o n _ p e r c e n t a g e =50 , s h a d o w s =False, \ 184 a r g s =None, t r a n s p a r e n t =0 , i m a g e =F a l s e):

185 """

186 S t a r t s the r e n d e r p r o c e s s and s a v e s the m o v i e

187 to a f i l e of f o r m a t AVI w i t h the H 2 6 4 c o d e c

188 """

189

190 # S e t t i n g r e n d e r o p t i o n s

191 p r i n t(" S e t t i n g r e n d e r o p t i o n s ... ")

A. The Code

192 s c e n e = bpy . c o n t e x t . s c e n e 193

194 s c e n e . f r a m e _ s t a r t = f r a m e _ s t a r t 195 s c e n e . f r a m e _ e n d = f r a m e _ e n d 196

197 r e n d e r = s c e n e . r e n d e r 198

199 if not i m a g e :

200 r e n d e r . i m a g e _ s e t t i n g s . f i l e _ f o r m a t =" A V I _ J P E G "

201 r e n d e r . f f m p e g .f o r m a t = " AVI "

202 r e n d e r . f f m p e g . c o d e c = " H 2 6 4 "

203

204 # s e t t i n g s h a d o w s

205 r e n d e r . u s e _ s h a d o w s = s h a d o w s 206 if t r a n s p a r e n t :

207 r e n d e r . a l p h a _ m o d e =" T R A N S P A R E N T "

208

209 # s e t t i n g o u t p u t

210 r e n d e r . u s e _ p l a c e h o l d e r = T r u e 211 r e n d e r . u s e _ f i l e _ e x t e n s i o n = F a l s e 212 if not " . " in f n a m e :

213 if a r g s :

214 f n a m e += " _nF { : 0 2 d } ".f o r m a t( a r g s . n f r a m e )

215 if not i m a g e :

216 r e n d e r . f i l e p a t h = f n a m e + " . avi "

217 e l s e:

218 r e n d e r . f i l e p a t h = f n a m e + " _ {}. png ".f o r m a t( f r a m e _ e n d -1) 219 e l s e:

220 r e n d e r . f i l e p a t h = f n a m e 221

222 r e n d e r . r e s o l u t i o n _ p e r c e n t a g e = r e s o l u t i o n _ p e r c e n t a g e 223 r e n d e r . r e s o l u t i o n _ x = int( r e n d e r . r e s o l u t i o n _ x )

224 225

226 p r i n t(" S t a r t i n g r e n d e r p r o c e s s ")

227 if i m a g e :

228 bpy . ops . r e n d e r . r e n d e r ( w r i t e _ s t i l l =T r u e) 229 e l s e:

230 bpy . ops . r e n d e r . r e n d e r ( a n i m a t i o n =T r u e) 231

232

233 def c l e a n u p ():

234 """

235 R e m o v e s all c u r r e n t l y e x i s t i n g spheres , c u b e s and c y l i n d e r s . 236 V e r y h e l p f u l if t h i s s c r i p t is e x e c u t e d f r o m i n s i d e b l e n d e r . 237 """

238 bpy . ops .o b j e c t. s e l e c t _ a l l ( a c t i o n =’ D E S E L E C T ’) 239 for i in bpy . d a t a . o b j e c t s :

240 if " S p h e r e " in i . n a m e or " C u b e " in i . n a m e or \ 241 " C y l i n d e r " in i . n a m e :

242 bpy . d a t a . o b j e c t s [ i . n a m e ]. s e l e c t = 1 243 bpy . ops .o b j e c t. d e l e t e ()

A. The Code

244

245 def c r e a t e _ c y l i n d e r _ m a t ():

246 """

247 C r e a t e s ( if not a l r e a d y e x i s t a n t ) and r e t u r n s a m a t e r i a l for 248 the c y l i n d e r s , t h a t are u s e d to d i s p l a y the c o n t a c t p o i n t s . 249 """

250 try:

251 r e t u r n g e t _ c y l i n d e r _ m a t () 252 e x c e p t:

253 p a s s

254

255 c y l i n d e r _ m a t = bpy . d a t a . m a t e r i a l s . new ( n a m e =" C y l i n d e r ") 256 c y l i n d e r _ m a t .t y p e=" S U R F A C E "

257

258 c y l i n d e r _ m a t . d i f f u s e _ c o l o r = (255 ,0 ,0) 259 r e t u r n c y l i n d e r _ m a t

260

261 def g e t _ c y l i n d e r _ m a t ():

262 """

263 R e t u r n s the c y l i n d e r m a t e r i a l c r e a t e d by

264 ’ c r e a t e _ c y l i n d e r _ m a t ’ if e x i s t i n g , o t h e r w i s e an a p p r o p r i a t e

265 E R R O R w i l l o c c u r .

266 """

267 r e t u r n bpy . d a t a . m a t e r i a l s [" C y l i n d e r "]

268

269 def r e m o v e _ a l l _ c y l i n d e r s ():

270 """

271 R e m o v e s all o b j e c t s t h a t c o n t a i n " C y l i n d e r " in t h e i r n a m e . 272 H e l p f u l if s c r i p t is u s e d and t e s t e d via the B l e n d e r GUI 273 """

274 bpy . ops .o b j e c t. s e l e c t _ a l l ( a c t i o n =’ D E S E L E C T ’) 275 for i in bpy . d a t a . o b j e c t s :

276 if " C y l i n d e r " in i . n a m e :

277 bpy . d a t a . o b j e c t s [ i . n a m e ]. s e l e c t = 1 278 bpy . ops .o b j e c t. d e l e t e ()

279

280 def c h e c k _ i n _ e l l i p s o i d ( vec , obj , o b j C o o r d = 1 ) : 281 """

282 C h e c k s if a v e c t o r vec is c o n t a i n e d in the e l l i p s o i d obj . 283 """

284 vec = vec . c o p y ()

285 q u a t = Q u a t e r n i o n ( obj . r o t a t i o n _ q u a t e r n i o n ). c o p y () 286 m a t r i x = q u a t . t o _ m a t r i x (). c o p y ()

287 m a t r i x . i n v e r t () 288 loc = obj . l o c a t i o n 289 s c a l e = obj . s c a l e 290

291 # If vec is in the c o o r d i n a t e s y s t e m of the e l l i p s o i d no 292 # t r a n s f o r m a t i o n has to be d o n e

293 if o b j C o o r d :

294 r e t u r n g e t _ n u m b e r ( vec , s c a l e ) 295

A. The Code

296 vec = vec - loc

297 vec = m a t r i x * vec 298

299 r e t u r n g e t _ n u m b e r ( vec , s c a l e ) 300

301 def g e t _ n u m b e r ( vec , s c a l e ):

302 """

303 C a l c u l a t e s the e l l i p s o i d f u n c t i o n and r e t u r n s the v a l u e 304 """

305 ret = 0

306 for i in r a n g e( 3 ) :

307 ret += ( vec [ i ]/( s c a l e [ i ] ) ) * * 2 308 r e t u r n ret

309

310 def h i g h l i g h t _ a l l _ i n t e r s e c t i o n s _ e l l i p s o i d ( p a r t i c l e s , objects , \

311 alpha , L , gamma , e l e m e n t s = [] , \

312 mat =None, f r a m e =1 , l i m i t s =N o n e):

313 """

314 T h i s f u n c t i o n u s e s the c u s t o m P y t h o n m o d u l e in o r d e r to 315 h i g h l i g h t the i n t e r s e c t i o n b e t w e e n the E l l i p s o i d s w i t h a

316 red c y l i n d e r .

317 """

318 l e n g t h = len( p a r t i c l e s ) 319

320

321 n _ e l e m e n t s = len( e l e m e n t s )

322 n _ u s e d = 0

323

324 # I t e r a t e t h r o u g h all p a r t i c l e s 325 for i in r a n g e( l e n g t h ):

326 for j in r a n g e( i +1 , l e n g t h ):

327 # d e f i n e the a c t u a l p o s i t i o n of the p a r t i c l e , as t h i s 328 # m i g h t h a v e b e e n c h a n g e d by the u p d a t e _ p o s m e t h o d .

329 p a r t _ i = []

330 p a r t _ i . e x t e n d ( o b j e c t s [ i ]. l o c a t i o n ) 331 p a r t _ i . e x t e n d ( p a r t i c l e s [ i ] [ 3 : ] ) 332

333 p a r t _ j = []

334 p a r t _ j . e x t e n d ( o b j e c t s [ j ]. l o c a t i o n ) 335 p a r t _ j . e x t e n d ( p a r t i c l e s [ j ] [ 3 : ] ) 336

337 # c h e c k if the two p a r t i c l e s are t o u c h i n g

338 o v e r l a p = g e t _ i n t e r s e c t i o n _ e l l i p s o i d _ o v e r l a p ( part_i , i , \

339 part_j , j , alpha , L , g a m m a )

340 if not o v e r l a p : c o n t i n u e 341

342 e l e m e n t = N o n e

343 if n_used < n _ e l e m e n t s : e l e m e n t = e l e m e n t s [ n _ u s e d ] 344

345 # r e s h o w the c y l i n d e r s t h a t w i l l be r e u s e d t h i s f r a m e s 346 e l e m e n t = h i g h l i g h t _ i n t e r s e c t i o n _ e l l i p s o i d ( o b j e c t s [ i ] , \

347 o b j e c t s [ j ] , overlap , element , frame , i , j , l i m i t s )

A. The Code

348 if e l e m e n t :

349 if n_used >= n _ e l e m e n t s :

350 e l e m e n t s . a p p e n d ( e l e m e n t )

351 n _ u s e d +=1

352

353 # Set all c y l i n d e r s t h a t are not u s e d to h i g h l i g h t an 354 # i n t e r s e c t i o n t h i s f r a m e to h i d e .

355 for i in r a n g e( n_used , n _ e l e m e n t s ):

356 e l e m e n t s [ i ]. h i d e _ r e n d e r =T r u e 357 e l e m e n t s [ i ]. h i d e =T r u e

358 e l e m e n t s [ i ]. k e y f r a m e _ i n s e r t ( d a t a _ p a t h =" h i d e ", i n d e x = -1) 359 e l e m e n t s [ i ]. k e y f r a m e _ i n s e r t ( d a t a _ p a t h =" h i d e _ r e n d e r ", \

360 i n d e x = -1)

361

362 r e t u r n e l e m e n t s 363

364 def h i g h l i g h t _ i n t e r s e c t i o n _ e l l i p s o i d ( obj1 , obj2 , overlap , element , \

365 f r a m e =1 , i =None, j =None, l i m i t s =N o n e):

366 """

367 H i g h l i g h t s one s p e c i f i c i n t e r s e c t i o n if the v i s i b l e o b j e c t s

368 are i n t e r s e c t i n g too .

369 N o t e : T h i s is c h e c k e d in o r d e r to p r e v e n t m a r k i n g s t h a t 370 w o u l d a p p e a r on ’ b o r d e r ’ - p a r t i c l e s .

371 """

372 g l o b a l V E R B O S E 373

374 e = c h e c k _ i n _ e l l i p s o i d ( V e c t o r ( o v e r l a p [’ r _ c 1 ’]) , obj1 , 0) 375 f = c h e c k _ i n _ e l l i p s o i d ( V e c t o r ( o v e r l a p [’ r _ c 2 ’]) , obj2 , 0) 376

377 # To c h e c k if the v i s i b l e o b j e c t s are t o u c h i n g , too 378 if ( e -1 > 1 e -6 or f -1 > 1 e - 6 ) :

379 if VERBOSE >=1 and j !=N o n e and i !=N o n e:

380 p r i n t(" Not c o n t a c t i n g b e t w e e n o b j e c t {0} and {1} " \ 381 " ( at l e a s t in b l e n d e r ) ".f o r m a t( i , j ))

382 r e t u r n N o n e

383 384

385 loc = t u p l e( o v e r l a p [’ r_c ’]) 386

387 # C a l c u l a t e the r o t a t i o n of the c y l i n d e r , so t h a t it is n i c e l y 388 # a d j u s t e d to the s u r f a c e of the c o n t a c t i n g p o i n t s .

389 d21 = V e c t o r (( o v e r l a p [’ r_c ’][0] - o v e r l a p [’ r _ c 1 ’][0] , \ 390 o v e r l a p [’ r_c ’][1] - o v e r l a p [’ r _ c 1 ’][1] , \ 391 o v e r l a p [’ r_c ’][2] - o v e r l a p [’ r _ c 1 ’] [ 2 ] ) ) 392 q u a t = d21 . t o _ t r a c k _ q u a t (" Z ", " X ")

393

394 d i s t = d21 . l e n g t h

395 r _ m a x = max(max( o b j 1 . s c a l e ) , max( o b j 2 . s c a l e )) 396

397 # T h i s s c a l e is a a p p r o x i m a t e of the r a d i u s r e q u i r e d for the 398 # c y l i n d e r s . T h i s is c a l c u l a t e d , as if t h e r e w e r e two S p h e r e s 399 # of r _ m a x t o u c h i n g

A. The Code

400 s c a l e = 1 . 3 * m a t h . s q r t ( ( r _ m a x *2 - d i s t / 2 ) * d i s t /2) 401 if scale < 0 . 0 8 * r _ m a x :

402 s c a l e = 0 . 0 8 * r _ m a x 403

404 s c a l e = [ scale , scale , r _ m a x * 0 . 0 5 ] 405

406 mat = c r e a t e _ c y l i n d e r _ m a t () 407

408

409 r e t u r n c r e a t e _ o r _ u p d a t e _ e l e m e n t ( loc , scale , quat , mat , frame , \

410 element , l i m i t s )

411

412 def g e t _ i n t e r s e c t i o n _ e l l i p s o i d ( ia , ii , ja , jj , alpha , L , g a m m a ):

413 """

414 R e t u r n s the c o n t a c t p o i n t b e t w e e n ia and ja or 0 if

415 t h e r e is no c o n t a c t .

416 """

417 r c _ c o n t a c t = a r r a y . a r r a y (’ d ’, [0 ,0 ,0])

418 if c h e c k ( r c _ c o n t a c t , [0] , ia , ii , ja , jj , alpha , L , g a m m a ):

419 r e t u r n r c _ c o n t a c t 420 r e t u r n 0

421

422 def g e t _ i n t e r s e c t i o n _ e l l i p s o i d _ o v e r l a p ( ia , ii , ja , jj , \

423 alpha , L , g a m m a ):

424 """

425 R e t u r n s a d i c t i o n a r y w i t h m o r e i n f o r m a t i o n a b o u t the 426 c o n t a c t b e t w e e n the p a r t i c l e s .

427 """

428 r c _ c o n t a c t = a r r a y . a r r a y (’ d ’, [0 ,0 ,0])

429 o v e r l a p = c h e c k _ g e t _ o v e r l a p ( r c _ c o n t a c t , [0] , ia , ii , ja , jj , \

430 alpha , L , g a m m a )

431 r e t u r n o v e r l a p 432

433 def c r e a t e _ o r _ u p d a t e _ e l e m e n t ( loc , scale , quat , mat =None, f r a m e =1 , \

434 c y l i n d e r = None, l i m i t s = N o n e):

435 """

436 U p d a t e s the o b j e c t p r o p e r t i e s for the g i v e n f r a m e . If no 437 o b j e c t is p a s s e d a new c y l i n d e r w i l l be c r e a t e d .

438 """

439 if not c y l i n d e r :

440 # C r e a t e a new c y l i n d e r if non p r o v i d e d 441 bpy . ops . m e s h . p r i m i t i v e _ c y l i n d e r _ a d d () 442 bpy . ops .o b j e c t. s h a d e _ s m o o t h ()

443 c y l i n d e r = bpy . c o n t e x t . s c e n e . o b j e c t s . a c t i v e 444

445 c y l i n d e r . r o t a t i o n _ m o d e = " Q U A T E R N I O N "

446

447 if mat : c y l i n d e r . d a t a . m a t e r i a l s . a p p e n d ( mat ) 448

449 for i in r a n g e( f r a m e ):

450 c y l i n d e r . h i d e =T r u e

451 c y l i n d e r . h i d e _ r e n d e r = T r u e