360.246 - Assignment 1 - Solved

General information

•   You may use any compiler of your liking, however we recommend additional warning output. When using GCC, at least the following flags should be used: -std=c++11 -Wall -Wextra -pedantic

•   Make sure you structure your programs in a readable fashion. It will make you more efficient and help us identify possible errors easily, i.e.: save you marks if there are problems.

1           Dense Grid Signed Distance Function
 The signed distance function Φ(~x) is a method by which a surface can be defined implicitly. It gives the shortest distance of the point ~x to the surface it describes. The sign of the distance relates to the point being outside or inside the surface. Although there is no strict convention, most simulation frameworks use negative values to describe points inside the surface, which is the convention you should follow.

1.1         Task
Implement a dense rectangular grid on which to store the signed distance function of a surface. The grid should fulfil the following:

•   The grid extent in each dimension, as well as the grid spacing(∆x), should be variable • Both reflective and periodic boundary conditions should be implemented (see Fig. 2)

 
Figure 1: Dense grid with
• Block allocated memory should be used for your data structures
highlighted boundary points.
Represent the following surfaces defined by parameters shown in parentheses, by filling your grid with the respective signed distance function:

•   Sphere (centre, radius)

•   Axis-parallel rectangle (minimum Corner, maximum Corner)

In your report, describe how you generated each surface and how the grid resolution of the implicit representation as a level set influences the accuracy of the surface description. Your program should be callable using:

./grid x-size y-size spacing [Sphere | Rectangle] [parameters]

Hint: Use helper functions to translate neighbour points outside of the grid to indices inside your grid according to the chosen boundary condition. If the boundary condition is reflective, a point at xmax +d should be mapped to xmax −(d−spacing). If it is periodic, the same point should be mapped to xmin +(d−spacing).

Hint: In order to analyse your results, you may use the provided vtkOutput.hpp file to generate a VTK Rectilinear Grid. The data can then be visualised using ParaView, which also provides functionality for analysing the data, e.g.: generating the explicit surface (zero level set) under Filters− > Alphabetical− > Contour using the contour value 0.

                                                Reflective                                                                                      Periodic

 

Figure 2: The centre rectangle with white background is the simulation domain, where a signed distance function is stored with the surface (the zero level set) shown in red. The grey rectangles are the neighbourhood of the simulation domain for reflective and periodic boundary conditions. The neighbourhood is here only shown to visualise the boundary conditions and will not be represented in your grid, which only consists of the centre rectangle.

2           Geometric Variables of a Signed Distance Function
During physical simulations, the geometric parameters of a surface, such as the surface normal or the curvature, are often required. In implicit surface representations, these need to be generated from the signed distance function.

2.1         Numerical Derivatives
 Due to the properties of signed distance functions, geometric variables can be approximated using its derivatives. Numerically, they are calculated using the forward  and backward  differences

                                             ,                                 (1)

where Φ(~x) is the signed distance function at ~x, ~ei the vector to the next

grid point in direction i and ∆x the grid spacing. The numerical derivative is then given by the central difference, which is the average of the two:
Figure 3: Next neighbours of point at ~x in a dense grid.
                                                    .                                        (2)

2.2         Surface Normal
Since the signed distance function increases linearly away from the surface, the direction of maximum change(i.e. the derivative) gives the normal vector to the level set:

                                                                                                .                                                                                     (3)

Therefore, the i-th component of the normal vector can be approximated numerically using central differences:

                                                                                    ,                                                                        (4)

where D is the number of dimensions used in the simulation, in our case D = 2.

2.3         Curvature
Since the curvature is essentially the second derivative, it is given by

                                                                                .                                                                    (5)

Therefore, the simplest approximation to the curvature is the sum of central differences of the normal vectors around the current point:

                                                                                                                (6)

2.4         Task
Implement functions to calculate the normal vector and curvature for any given point in the dense signed distance grid described earlier. The program should output the normal vector and the curvature at point (x,y) when invoked as follows:

./grid x-size y-size spacing [Sphere | Rectangle] [parameters] x y

3           Moving the Zero Level Set to Advance a Surface
Since the explicit surface is defined by the zero level set, the surface can be moved by changing the signed distance function. In microelectronic process simulations this would correspond to the passage of time during a fabrication step which etches or deposits on a surface. The simplest way to do this is by adding a constant value to the signed distance function, which will simply shift the interface along the surface normals. Therefore, the time evolution is given by

                                                                                               ,                                                                                   (7)

where V (~x) is a scalar velocity field describing the movement of the surface at any point in the domain.

However, for curved surfaces this would lead to the signed distance function losing its normalisation, which is why the scalar field has to be normalised using the gradient of the signed distance function:

                                                                                      .                                                                         (8)

Therefore, the change of the values in the signed distance function in one time step is
 
                                                                     ∆Φ(~x) = −V (~x) | ∇Φ(~x) | ∆t     .
(9)
Since the right-hand side of Eq. (8) is a Hamiltonian, many numerical schemes exist to solve it, the simplest being the Engquist-Osher scheme

                                    ,                 (10)

which will result in a properly normalised signed distance function in the whole domain, as long as the time steps taken are small enough to satisfy the Courant-Friedrichs-Lewy (CFL) condition

                                                                                          .                                                                           (11)

Therefore, in order to advance the simulation time by t, equation Eq. (10) has to be solved t/∆t times to guarantee a stable movement of the surface.

3.1         Task
Investigate how advancing the surface changes the signed distance function. Compare the following ways of advancing the surface:

•   Simply subtracting a velocity value from every grid point as in Eq. (7)

•   Several steps of the Engquist-Osher scheme Eq. (8) - Eq. (11)

Compare the results for V (~x) = const for a positive velocity greater than several grid spacings for the following shapes:

•   Sphere (radius of 10) for V = 10 at the resolutions ∆x = (1,0.25) for the times t = (0.1,1)

•   Rectangle (with side lengths 5, 20) for V = 10 at the resolutions ∆x = (1,0.25) for the times t = (0.1,1) Describe what effect different advection schemes have on the result and how the grid spacing influences the results of each scheme.

Use the normal vectors at each point to introduce a vector velocity function  ), which can move the surface directionally. The scalar velocity value used in all algorithms, V (~x), can be generated easily by taking the dot product of the vector velocity and the normal vector at that point:

                                                                                     V (~x) = V~ (~x) · ~n(~x)                                                                             (12)

Use the normal vectors and curvatures generated earlier to investigate the following, using the Engquist-Osher scheme:

•   How does the rectangle behave over many time steps when it is moved by the vector velocity function V~ (~x) = (1,0) ?

•   Investigate what happens when the curvature is used as the velocity,

i.e. V (~x) = −κ(~x). What could this be used for?