Fracture Formation and FEM

We’ll be looking at the phenomena of fracture formation, and modeling it using the finite element method based on this thesis by P. Federl. The general idea is that we have a material which consists of a surface with two layers that are connected.  The bottom layer grows, causing the top layer to deform.  When the amount of deformation passes a threshold, a fracture is formed.  You can assume most if not all of the material presented is from the thesis, while any mistakes are my own.

Tensors and Growth

The growth is modeled using a growth tensor of rank 2.  Think of it as a matrix, or an extension of a vector.  A vector has magnitude and one direction.  A tensor of rank 2 has magnitude and 2 directions associated with it.  If we have a vector and we want to obtain another vector, we can multiply it by a scalar or a matrix.  The matrix can change the direction, while the scalar cannot (except for transforming it in the opposite direction).  Here is more info on tensors if you’re curious, they’re pretty wild.  We can solve for the velocity of growth at any point from this matrix by solving for the velocity components of the following set of linear equations, where T is the growth tensor.  The velocity components can then be calculated by a given position.

So for any point on the surface, we can get a growth vector, and this gives us a field of velocities.

The above is from the thesis, and shows you different types of growth that can occur on a surface. Isogonic growth is growth where the rate of growth is the same in all directions at all locations.  If I apply this to a circle, it will still be a circle.  Uniform anisotropic growth is that the growth is the same at all locations, but changes given a direction.  If I apply this to a circle, it will become deform into another shape.  Note that the growth operates on the shape such that the shape is at the origin.

In this example we’ll simply use growth along the x axis, where R is the rate of expansion.

And the corresponding growth tensor is simply

Creating the Surface

We begin by creating a triangulated surface with semi-random triangles.  Because of the nature of FEM, uniform triangulation and random triangles can have undesirable effects.  With random triangles, some points can be too close.  Place the nodes on the surface randomly, and repel them.  I used the separation function we made in the Boids algorithm, adding a drag force proportional to velocity.  Run this until all of the nodes are in equilibrium.  You can implement this as running the loop until all of the velocities are zero, or you can choose to run the loop for a certain amount of time as the particles separate (but this does not guarantee the particles will fully separate evenly). Now triangulate the nodes using a Delauney triangulation.

This gives you a triangular mesh on a plane.  Project this plane in the normal direction to get a second layer.  Each triangle now becomes a wedge element.  Each wedge element has six nodes.

We are going to focus on this element for now. Then we will see how to combine all of them to calculate the stress on the nodes.

Finite Element Method – Qualitative

The finite element method can provide approximate numerical solutions to a wide range of engineering problems.  This is most useful when encountering problems where an analytic solution is difficult to find, or simplifying assumptions cannot be made. We begin first with something familiar: finite differences.  These are also used to find approximate numerical solutions.  If we have a shape, we can approximate it using a grid to cover the shape.  More points (a finer grid) produces a more accurate covering.  The finite element method is similar, only it uses triangles instead and can much better approximate the region.

The above figure is from here, which provides a good introduction to the ideas of FEM. When we created smooth particle hydrodynamic simulations, we approximated continuous functions and properties over a region using kernel functions and interpolation.  Similarly with FEM we divide the region into nodal points and use interpolation functions to approximate the property of interest over the region.

The FEM allows one to divide a region into elements, calculate the solutions for individual elements, and then combine the element solutions to form a soltuion for the region.  In our case we find the stiffness of each element, and then combine them to obtain the stiffness of the entire material.

There are a variety of ways to calculate the properties of the individual elements.  We choose the direct approach.  Other approaches include the variational approach, and the weighted residuals approach.  We divide the region into elements, calculate the element properties with an interpolation function with the element nodal points, then combine each of these element matrices into a giant matrix to represent the entire system. We impose the boundary conditions (any values we know about the system and the nodes) and then solve the matrix equation for the system.

Finite Element Method – Quantitative

So from the perspective of an element, the amount of force generated is proportional to the displacement of the nodes.  Sound familiar?  In the matrix form is looks similar to F=kx, Hooke’s law to determine force for a spring.

Where F is the force of the element, K is its stiffness, and q is its nodal displacements.  For each element we have six nodes.  Each node has three degrees of freedom, in the x direction, y direction, and z direction.  This gives us a total of 18 degrees of freedom.  [F] is an 18 x 1 vector with F_x, F_y, and F_z for each node.  Similarly [q] is an 18×1 vector with the corresponding q_x, q_y, q_z for each node.

If we want to determine the value of some function over the entire element, it is helpful to define a coordinate system within the element.  These coordinates would be called isoparametric coordinates. We can define three variables to integrate over in the wedge.

We can now use interpolation with the six nodes to determine the value of any point within the element, as a summation of the value at the node times a weighting function N.  The weight is called a shape function, and is some value that depends on the isoparametric coordinates, and is inbetween 0 and 1.  Mathematically this looks like

And if we have a matrix of N_i * I and concatanate them for i = 1, …, 6, where ‘I’ is the identity matrix, we have an 3×18 matrix for N.  Given V is the 1×18 vector of concatanated values for each of the nodes, we can represent the above as a matrix equation [v] = [N][V].  If Q is the 1×18 vector of position points of all the nodes, we have [x,y,z]=[N][Q], converting the coordinates in the element to cartesian coordinates. Given a displacement [u,v,w] in an element, we can interpolate it over the nodes in an element using the shape function N.

Where u, v, w, and N are dependent on the isoparametric coordinates.

Before contintuing onward to calculate the element stiffness matrix, it is helpful to have a better understanding of the shape function.  It should be pretty familiar in concept, and fortunately it is much simpler than the analog we used in SPH.

For any of the shape functions, you can vary the coordinates and see how it changes within the element. We can integrate over these coordinates to find [K_e], the element stiffness matrix.  We can approximate the integral using the Gaussian integration technique.  Instead of computing the continuous integral, we compute the discrete sum over select points that are weighted.

Now that we know how to go about integrating, we need to determine what exactly to integrate over.  To find [K] we integrate over two other values, the material stiffness matrix [D], and the matrix to relate displacement with strain, [B].  The material is uniform throughout, so the stiffness matrix is I * constant.  We can derive strain in each direction as a summation over the shape fuctions, since we can relate the shape functions and the displacement. There is also the Jacobian matrix [J] that transforms between coordinate systems (cartesian and isoparametric). The final equation for [K_e] is

where we use n Gaussian integration points.

Note we can rewrite the matrix equation for stiffness as [F] = [K][P]+[R] if we notice [q] = [P] – [Q] , where [Q] is the undeformed coordinates.  So to obtain the global stiffness equations we must find

For a system of n nodes, F_g is a vector of length 3n, since each node has three force components.  Similarly the coordinate vector [P_g] has length 3n. The global stiffness matrix [K_g] is simply adding the appropriate values at the correct indices of the element stiffness matrices.

In our simulation the bottom layer is fixed and the top layer is free to move.  We set some initial strain (the growth of the bottom layer), insert the boundary conditions (force on all nodes is zero), and solve the equations (calculate the nodal coordinates of the free nodes.  We now solve the familiar equation

Our set of linear equations is pretty sparse, since we created the global stress matrix by simply combining the element matrices.  Because of this, we can use the conjugate gradient method to solve the set of linear equations.

We’ll take a brief moment to discuss a few concepts.

And here is an algorithm for it, from Wikipedia.

\mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0  \,
\mathbf{p}_0 := \mathbf{r}_0 \,
k := 0 \,
repeat
\alpha_k := \frac{\mathbf{r}_k^\mathrm{T}  \mathbf{r}_k}{\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k}  \,
\mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k  \mathbf{p}_k \,
\mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k  \mathbf{A p}_k \,
if rk+1 is sufficiently small then exit loop end if
\beta_k := \frac{\mathbf{r}_{k+1}^\mathrm{T}  \mathbf{r}_{k+1}}{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}  \,
\mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k  \mathbf{p}_k \,
k := k + 1 \,
end repeat
The result is xk+1

We calculate the eigenvalues and eigenvectors of the stress tensor on the mesh to determine the principle stresses at given points.  If the largest of these stresses exceeds the material’s stress, a fracture is formed.

There are a few ways to represent the fracture.  One way is to delete the element.  This can produce fractures, though it is not realistic. Another method is to subdivide the elements along the principle stress plane.

To optimize we can use dynamic subdivision. If a given element is a candidate for fracture, we subdivide and recalculate the stresses.  Once an element is small enough, a fracture can form on element relaxation.

More info will be posted soon!

Emergence: Fireflies

One of the most intriguing phenomena within complexity science is that of emergence.  It is a salient feature of complex systems, so much so that definitions of complex systems include the quality of exhibiting emergence.  Many exciting solutions to the world’s bigger problems can be understood through collective dynamics and emergence.

Google’s search algorithm, to a heavy extent, depends on the nature of linking between websites.  On the scale of the website and what links to it, Google is able to provide an answer the global problem of findability on the internet.

Muhammad Yunus developed a micro loan system that re-frames loans to take into account collective dynamics of a group.  Through challenging the assumption of who one can give a successful loan to, and by tapping into group trust, he has created an effective method of alleviating poverty.

What is Emergence?

Simply put, emergence is when behavior at a smaller scale of a system produces global behavior that is not entirely intuitive given the behavior at the smaller scale.  The global behavior is not directly programmed in the behavior at the local level. Something surprising happens at the larger scale.  From simple local interactions, we can arrive at global behavior that is immensely complex, sometimes seemingly random. Think of the classic example of Conway’s Game of Life. The rules to define the simulation are simple, but the behavior can be pretty unpredictable. Or on the other hand, we can think of systems which have a lot of randomness at the smaller scale, but order emerges at the larger scale. This underlies the concept of self-organization, a process which is widespread in the development of biological and societal form.

Another way to think of emergence involves the notion that the whole is greater than the sum of its parts. The behavior of the system is a result of not only its parts, but their interactions. The actions at the level of the individual do not imply that of the behavior of the system. We need to understand not only the system at the individual level, but at the level of multiple individuals to observer the effect of their interactions. One of the most exciting examples of emergence involves evolution. At a smaller scale we have natural selection, species in a population surviving or dying over time based on selection pressures from the environment, with possible species variations through reproduction and sexual selection and the random chance of mutation. At a larger scale we have the development of biodiversity, a wide range of organisms from bacteria to whales, and an environment with so many different species able to thrive within different niches. Global biodiversity is not programmed in the local interactions; it is an emergent property.

Why is it important?

Understanding emergence is essential when working with complex systems. For the most part, complex systems exist because they are grown, i.e. they have evolved to be what they are. Think about what makes up a complex system. They are composed of many parts that interact across scales, from individuals to groups to the whole. When you think of a city, it grows through the interaction of all of the individuals and organizations which comprise the city. The city is not formed by a single individual, it is formed by its population. Emergent behavior can be difficult to understand, since by its very nature it is unintuitive and unpredictable (the very reason why it has received so much attention). In understanding global behaviors of a system, it helps to understand how individual interactions led to the global behaviors. In this way, one can attempt to reproduce the global behaviors through recreating the individual interactions.

Example of Emergent Behavior: Fireflies

A popular example which illustrates emergent phenomena is that of fireflies which can spontaneously synchronize. There is no head firefly with the ability to communicate with all of the others, so how is it possible that they can flash in synchrony, when no firefly can see all of the rest? In their paper ‘Synchronization of pulse-coupled biological oscillators’, R. Mirollo and S. Strogatz seek to answer this question by developing a model that introduces a mechanism of simple interactions between the fireflies. Through this local interaction, the global behavior is that of synchrony. The key concept with respect to emergence is that the local interactions do not mention anything about synchronizing. Each firefly can be thought of as an oscillator which fires at time T, and resets to zero only to fire again at time T. All of the fireflies have the same period T, but start at different times in their period. At each step, all of the fireflies increment in their period. A firefly has the ability to see neighbors within a certain detection radius. When one of the neighbors of a firefly flashes, the firefly increments in its period, leading it closer to flashing. The nature of this firing response is important. Say that a firefly is at time t’ in its period. The firing response is defined as follows: the amount of increment increases as t’ increases. So if a firefly is later in its period, the increment will be higher. If its just beginning its period, the increment will be lower. The next is not so important for the concepts, but mathematically it can be more accurately stated as such:

Here, t” is the new internal time if a firefly experiences a flash at time t’. As defined in the model, f is the firing function and epsilon is a small constant < 1 .
From these specifications, the fireflies will synchronize. The synchronization can be viewed in the following applet. I approximate the increment function as described in the paper ‘Firefly-Inspired Sensor Network Synchronicity with Realistic Radio Effects‘ Werner-Allen, et. al. as

Where epsilon is 0.1 .

Click on the image to see the simulation

The simulation is done in Processing. I encourage you to open it up and play with the number of fireflies and the period.
It is difficult to really see where in the period each firefly is, so I created another model which simulates the fireflies as if they are standing in a single line. They bounce up and down periodically, flashing when they hit the ground. If a neighbor flashes, they increment in their period, evident in a sharp jump. This one has an interface to allow you to change the number of fireflies, the period, and the detection radius (how far a firefly will look to its left and right to detect a neighbor.

Click on the image to see the simulation. This one includes an interface to play with!

Its pretty exciting to see the simulation converge to synchrony. From a simple rule given a neighbor flashing, organization emerges as the fireflies flash as a group. The group is able to act as a whole through local interactions.

Cool! … now what?

Think about environments that your are in or have seen, and reframe them in the systems perspective. How do the parts combine to form the whole. Is the emergent behavior obvious given local interactions?

What if there is a particular global quality you want to achieve. What are ways in which local interactions can achieve those global qualities? Because the global behavior can be unpredictable, it can become a difficult problem to determine the local interactions. But other times, those interactions can be found through the ingenuity of a group and global behaviors can be achieved. The question then also becomes, how do you create and environment that will promote such creativity? This is a question of emergence, in designing local interactions to promote global creativity.

Design and Dissemination

An example of the utility of thinking in terms of systems and emergence is in the world of creativity and design. If a designer and the consumer have similar environments, it is very likely the designer can find successful solutions for the consumer. In cases where the designer and the user are in completely different environments (including different experiences and customs), it will be immensely difficult for the designer to make something useful. In some cases, products will be developed, only to suffer completely different uses in the field (which is not always bad, but it is a gamble). So if instead the designer works with the user, together they can design something successful. From the complex process of improving experiences through design, the collaboration leads to the emergence of successful, effective solutions. This may appear intuitive, but it has not always been the case that those who believe they have solutions have taken the time to fully understand to conditions those who they are trying to help.

With respect to dissemination, the goal is to have emergent behavior that spreads information, or a product. The question becomes: what sort of local interactions will result in this global goal? There are many ways to approach this problem, but a few methods include understanding the topology of the society in disseminate within. Certain people are more likely to spread information or a product, so if one targets those people, the product will spread faster and further. Another approach is to focus on the role of the individual. If one is empowered to be included in the product design and construction, they will not only help in characterizing the problem, but will have a deeper connection with the product. It is more likely they will, in their enthusiasm and success, spread the information and product to others. The task becomes including those who use the product to become part of the design process, be it through business, construction, or design, and they will in their own interests spread the product.

An apt analogy is to think of the process of gardening and growth. The designer is a seed, and the soil are those who will use the product in whatever environment they are in. The soil contains the nutrients, that which the seed needs to form into a healthy plant. From the interaction of the two, a product is formed, the plant. Though the seed may be the same, in two different environments, the plant will be different. From the environment of collaboration, many solutions are possible. Such solutions may spawn more seeds that in turn develop more solutions.

This process is not too difficult to grasp. What becomes harder to fathom is the type of solutions developed in the aggregate. The variety of solutions and even ideas for possible solutions of the group far surpass that of the individual. It becomes difficult to predict the creative capabilities that emerge in a group.

Collectively we are capable of much more than any of us alone can imagine.

The above article describing emergence was created as a supplement for Eric Reynold’s presentation on creative capacity building and systems science at the 2010 International Development and Design Summit in Fort Collins, Colorado.

A Generative Method for Infrastructure Emergence

Social systems are becoming more complex from technological advancements and increased connectivity. Individuals are further empowered with the capability to augment their memory and communication through computers, the internet, and cell phones. Every society has structures which influence collective behavior, and with all of the possible configurations of people in a population, the question emerges for designers of how to implement a method to use the collective information and create a successful design solution [1].

Cities have been shown to have fractal geometry. In this paper we show how the fractal shape can emerge from a generative process that takes information on the scale of individuals or groups, and uses it to design a permanent infrastructure on the scale of a city. In this sense, we grow cities consisting of individuals and roads, starting from just individuals. [from introduction of paper, see PDF for rest]

A Generative Method for Infrastructure Emergence from Kawandeep Virdee on Vimeo.

Complexity and Creativity

A talk I gave at the December 2009 Research Club brunch at Tribute Gallery on applying ideas of complexity to creative projects.  Applications of complex systems ideas in the projects HEXAGON and PDX I Love You are described in the presentation.  The talks were limited to five minutes each, and a transcript is included below.

Brunch #1 / Lecture #3 / Kawandeep Virdee Talks About Complexity from Research Club on Vimeo.