 Hello everyone. Welcome to Dynamic Deformables, Implementation and Production Practicalities. So just some brief speaker introductions. I'm Ted Kim. I'm currently a professor at Yale, but from 2015 to 2019, I was a senior research scientist at Pixar where I designed and implemented a bunch of the algorithms that we'll be seeing today. My co-presenter, David Everly, is the simulation tools R&D lead at Pixar. So in addition to also implementing and designing some of the algorithms that you'll see today, he has been in charge of FISTEA off and on for the last six years. So FISTEA, what is FISTEA? FISTEA is Pixar's internal physics simulator. So FIS, as in physics, T is a tool. So this was originally David Baroff's production implementation of the classic 1998 paper, Large Steps and Cloth Simulation. So it's been in production use at Pixar for about 20 years now, and it's not exactly the same as it was back then. Of course, it has evolved quite a bit, and that's a lot of what we're going to describe today. All right. So let's see a few shots in which FISTEA was used. Okay. So here's some of the newer stuff. Here's a shot from Bow. This is the short that ran in front of Incredibles 2. All right. So the mom here, so you can see the clothing on her. So that was simulated. We actually ran a body sim on her as well. So her flesh was actually simulated as a pre-pass to make it easier for the cloth simulation to give useful results. All right. So here's also a shot from Toy Story 4. We did a skin simulation on the cat named Dragon here so that you can see the fur, the way it slides over the cat's body. That was actually a simulation as well. All right. And here's a shot from Toy Story 4. Woody's pull string, which you're going to see, he's going to use to last so forky. That was actually a volume simulation as well. All right. So these are just some of the improvements. And we've presented a lot of these improvements piecemeal at SIGGRAPH over the last few years. So the big thing that we've added here is the handling of solids. And that's one of the things that I'm going to go over. There's been lots of improvements also in collision handling and in fast matrix assembly. And that's the sort of stuff that Dave is going to go over. Okay. So it's not just me and Dave that worked on this. There were actually lots of collaborators and co-authors that we worked at and we'd like to acknowledge them all here. So if you've done any collaborative work, then you know that good work doesn't happen unless you have good collaborators. And we had a set of really great collaborators. So we'd like to call them out here. All right. Now before we get too far into this, I am going to warn you this is not an introductory class. So I'm going to give you a refresher on deformation mechanics at the very front. But if this is the first time you've ever seen this stuff, don't feel bad if you don't get it right off the bat. So there's quite a few good resources out there if you want to fill in the background on some of this more introductory material. So there's the classic notes on this. So there's the Wiccan-Bearaf notes from 2001. There's more recent stuff too. So Sofakis and Barbic had a very nice introduction to solid mechanics back in 2012. And Bartel and Shanar have been doing this introduction to physics-based animation course at SIGGRAPH for the last few years. And you can take a glance at those as well. So they're also very good. So you can get the basics on how this stuff works from these notes. And once you're comfortable with that, you can always come back here and see how it all fits together in production. Okay. So I'd like to mention there's also a very large document that accompanies this course. So sometimes you go to these courses at SIGGRAPH and they say we have course notes and you go look at them. And it's just a big PDF of the slides from the PowerPoint presentation. So that is not what these notes are. All right. So I'm going to show you the table of contents for these notes. That's this. So it actually is a very large text. It's 183 pages of text that goes into the concepts that we're going to go into in lots of detail. So actually, we're only going to scratch the surface of the contents of the course notes here. If you want more detail, there's lots and lots of detail in the course notes themselves. Okay. So with that out of the way, let's get to the main show. So if you've ever implemented your own simulator, you know that one big piece of that simulator is how to compute forces. And one of the big things that we added to FISD is support for volumes. So how do you compute the forces in a volume? So we're specifically going to look at internal forces. So how does a piece of material resist as you squash or you stretch it? All right. So this happens for triangles, if you have cloth. So if you try to stretch out a piece of cloth, it will resist you, right? The same thing happens with a volume, though. So we have a tetrahedron here, and we're going to squish it, right? Squish it. And then what happens? Well, it's going to exert some forces because it doesn't want to stay squished, right? And it's going to push back on the fingers that you're using to squish it. So how much does it push back? There's actually no one true way to do this. There's multiple ways to do it, to compute these forces. And we're going to go over our particular answer that we have in FISD for answering this question. Okay. So there's more than one kind of force that can appear in a simulation, right? So again, there's no one true way to handle this. And we're going to go over how FISD does this. But for example, if we have a collision, so if something collides with the tetrahedron, that is a different kind of force, right? So what forces appear in that case? And again, Dave will go over the FISD way to do this. All right. So again, we're dealing with a baref and wittgen-style solver. In fact, we're dealing with the actual original baref and wittgen-style solver. So one of the big pieces of that kind of solver is that it is an implicit integrator. So we do need these forces, which we saw before. There's another piece, though. So you actually need the derivative of those forces, the spatial derivative. So these are also known as the force known as the force Jacobians. All right. So if as a many student who's writing their first simulator has discovered, deriving the forces actually is not that easy the first time you try to do it. But the pain and the difficulty actually goes way up when you have to then compute the derivative of these forces. So this is the really painful part, figuring out what the force Jacobians should be. And if that wasn't bad enough, the first time you try this, it probably will explode pretty horribly on you because there's a few hidden gotchas here that aren't always that well disclosed. All right. So let's go over what we need in a little bit more detail. Okay. So we need the forces, right? This is going to be a 12 vector. So we have four nodes on our tetrahedron, and each of them have three degrees of freedom, right? So four times three is 12. All right. And we also need the Jacobian of those forces. So this is a 12 by 12 matrix. So it is a derivative of a vector with respect to a vector. And we get a matrix out of it. All right. So there's actually one other really hard piece, and this is why the simulator can sometimes explode even if you did do the math correctly. All right. You actually need the eigen decomposition of this matrix as well, the force Jacobian. Now, why do we need to do this? Oh, again, we're dealing with a barefoot Wiccan style solver. So you're going to use precondition and conjugate gradients, PCG as your numerical solver, right, in order to invert your matrix. Now, for PCG to work, your underlying system has to be semi-positive definite. So the eigenvalues of the system matrix all have to be greater than or equal to zero. Otherwise, that solver, it will explode. And that's what you will see oftentimes. Now, unfortunately, any elastic energy can produce negative eigenvalues. So actually, even the ones in the barefoot Wiccan paper. So they don't really mention this. It sort of pretends like maybe this won't happen, but it does. So any energy that's not just some zero length spring, it actually can produce negative eigenvalues. It can produce an indefinite system. So you do have to do some sort of filtering. And in fact, if you want to look at exactly when this happens in barefoot Wiccan cloth, I am going to present, I presented paper earlier this week at SCA that describes exactly how this happens. So you can go take a look at that paper if you're really curious about barefoot Wiccan cloth. All right. So that is why we need the eigen decomposition of each of the tetrahedron. So once you have the eigenvalues in hand, you have to filter them. So if any of them are negative, you need to make sure that those get snapped to zero instead. All right. So there's good news here. And that is that you can actually do this analytically. So I'm going to hand you some closed form expressions for the eigenvalues and the eigen matrices. So this is one of the main results from our paper last year. And if I can just pause for one second here, I think this is pretty amazing, right? So usually, if you want to get eigenvalues or eigenvectors, what have you, you sort of push it through MATLAB or you push it through eigen. You just get this more ass, this like primordial C of numbers and you can't really make heads or tail of it. You just chug ahead with the actual numerical values. And that's not what's happening here. We actually have closed form expressions here. You can actually see the expressions for what these are. I think it's pretty neat. Okay. So if you actually just want to skip to the end and you just want to see the results of all of our analysis, you can go to chapter seven and look at figure 7.7 in that chapter. There, I actually give you the MATLAB code that will do everything. So this code will actually directly spit out the analytic eigen composition of the forced Jacobian for any isotropic energy that you can think of. So it doesn't matter if you came up with your own or you used one out of a textbook, this one will give you the analytic eigenvalues. So if that's all you want from me, you can go ahead and just scrub this video and skip to Dave's part. But if you actually do want to understand what's going on, let's get started on that. All right. So as promised, we're going to begin with a brief refresher on deformation mechanics, just so that we're on the same page with the style of analysis that we're going to be using. All right. So we want to detect when an object has been deformed, i.e., it's been squashed or stretched, like the tet that we squashed right here. Now, based on how much it is squashed or stretched, we can then compute how much force that the tet will then resist with, right? How much it will try to go back to its original shape. So there's basic questions of how do we actually measure how much something is deformed? So what we could do is we could look at just the positions of the nodes on the tet and see how those have changed. Now, the problem here is that if you use a position-based measure, it's actually really easy to trick it into thinking that something that is not a deformation actually does look like a deformation. So, for example, if we're just subtracting these node positions from each other, if you just went and rotated a tetrahedron, it would suspiciously look like it has deformed. But we know it hasn't deformed, right? You haven't actually squished the tet at all. So we need a little bit better way to measure whether something has deformed or not. Okay, so because of that, we're going to use something called the deformation gradient. So this is a 3x3 matrix, and this matrix is derived from the node positions of the tet, but it is derived in a way that translation, so if you translate the tet, it can't be tricked into thinking that a translation is actually a deformation. So it is a better way to measure deformation. All right, so a big neon sign here. This is going to be our primary variable. So everything that I do from here on out, it's going to be based on the deformation gradient F. So forget about positions. F is the important thing. Now, I'm not going to tell you how to derive and compute F here, because that would take a little bit of time. Appendix D in the course notes actually goes over it in quite a bit of detail, and at the end of it, I actually even hand you C++ code to help you out. So if you're really curious about how to compute this, you can pause the video and go read that section or just go read it afterwards, and we will tell you. But for the time being, let's just assume that we have this matrix F in hand. All right, so F can tell us how much that this tet has deformed. So once we know that, then the next question is, well, what forces should appear at the vertices to try to restore this tet to its original shape? Well, so we can do this a barefoot can away, which uses positions. So we'll go over this first, then we'll look at how to do it the deformation gradient way. All right. So if we have our positions, we can stack them all into just one big vector here. It's a 12 vector. And the barefoot can away is then to define an elastic potential here. So this is a function that takes an X that that stacked vector of positions, and it spits out a single scalar. All right. So this seems a little odd. What does a scalar have to do with forces? So it turns out that if you take the negative gradient of this energy, then it will actually spit out the forces that you care about. Okay. And once you have those forces, you can just go ahead and apply them to the notes, like you see here. All right. So this is exactly what we wanted. Now, the thing to remember here is from vector calculus, if you have a scalar function, and you take the derivative of it with respect to a vector, you'll get a vector out of it. All right. So that was the position based way, which is not the way that we're going to do it. And we're going to do it based on deformation gradients. We're going to call this the finite element way. So let's take a look at this instead. So instead of node positions for using deformation gradient, we're still going to define an elastic potential here. But instead of taking X as an input, we're going to take F as the input. It still spits out a scalar though. Okay. All right. So then based on this, how do we get forces out of this function? So we can take the derivative with respect to F, but the derivative of a scalar with respect to a matrix is actually a matrix. So last time what happened was we got this nice list of forces out, right? And we're not getting that in this case. This is not a nice list of forces. What are we supposed to use this on our nodes? Okay. So we can actually convert it back to the list of forces using the chain rule. And when we do, then we get the vector back that we expect. And then again, we can just go ahead and apply these forces to each of the nodes. Okay. There's actually an evil professor trick that I'm doing here. And I'll come back to it later. But just for the time being, let's take for granted that this works. So let's take for granted you can write some elastic potential in terms of F. And with a little work, you can get some forces out of it. All right. So we've written down this function psi a few times now, these elastic energies. What do these actually look like? Can you give me an example of one of these elastic energies? Yes, I will. All right. I will in a second. So the traditional way to write down elastic energies is using what are called the Cauchy Green invariance. So these are three invariants that describe every kind of deformation that F could be trying to communicate to us. Now, I'm not going to go into detail over what, over how these are derived or how they're computed in too much detail, especially because they involve stuff like trace and determinant. So Section 5.2 in the notes, this goes into a lot of detail on this. So if you really want to know lots of detail about these invariants, you can look in that section. Instead, let's look at what these invariants intuitively describe instead of what they have their computed. Okay. What are these things describe? So what we want to know is what's going really going on inside this deformation gradient. So if you ever really want to know what's really going on inside a matrix, you take its SVD, its singular value decomposition. So let's do that here. We took the SVD of F, and let's see how this relates to the invariants. So what we really care about is the singular values here. So the diagonal matrix in the middle of the SVD. So what we're going to see is that these describe how much things have been stretched or squashed in different directions. Okay. Say that we have a cube of material here, and we can actually very intuitively interpret what the singular values are in terms of this. All right. So if we take each singular value, it actually corresponds to the length of the edges of this cube that we're squashing and squishing. All right. So we have the sigma zero here, and sigma one here, and then the last direction is sigma two. All right. So it's just the edges of the cube. So if we look at the first Cauchy-Green invariant, it's actually just the squared sum of the lengths of the cube. So it is asking how much did the lengths of this cube change after you squished and squashed it? All right. If you look at the second invariant, it looks at the areas of the faces of these cubes. All right. So you can see that this is the area of this top face here, the area of this front face, and the area of this side face. All right. So it is asking how much did the area of the faces of this cube change? And then we have one invariant left, and that is this one. And this third invariant, this measures the volume of the entire cube. All right. How much did the volume of the cube change? So we have these three Cauchy-Green invariants, and they tell you whether you've experienced any change in length, area, or volume. So that's a pretty good set of things to measure, right? If something has deformed, it's going to show up as a change in one of these things. That's for sure. So if you're going to try to figure out how your tet is deformed, then these three variables should tell you everything that you need to know. All right. So let's look at some actual elastic energies. So there's a bunch of ones here that you can see. Now, these aren't that pretty looking, right? They look pretty complicated. But if we drill down into them, actually, it's all just a bunch of mixing and matching of the three invariants. So it's a big balancing act. So it's asking, what do you care about the most? Do you care about area changing the most, about volume, about length changing? And based on what you actually care about, you're just mixing and matching those three terms, those three terms in some combination. All right. So now I'm going to toot the horn of our paper from last year a little bit. So I just showed you a bunch of elastic energies, and I said, yeah, it's just all big mixing and matching game of these three Cauchy-Green invariants. There's actually a whole bunch of other elastic potentials that people use in practice that you can't write in terms of these invariants. Now, the real bummer here is that arguably the most popular energy out here, ARAP, as rigid as possible, which is one that we use in FISTI in some form, it actually it's one of the ones that can't be written in terms of these invariants. This is a real bummer for lots of reasons. Okay. So why does this happen? How can we can't write it in terms of those invariants, which actually looks like they were a very reasonable set of measures. So we're going to do an expansion of the ARAP energy. This actually is the same as the equation from the previous slide, but we did a little bit of rewriting. And we look at this trace S term right here. So this is S from the polar decomposition of F. And if we look closely at what trace S actually means, it actually corresponds to the sum of the singular values. So the sum of the lengths of the edges of the cube. So what's the trouble here? Didn't we just see this? We just saw an invariant that actually does measure the change in lengths of the size of the cube, right? Actually we didn't. So this is the Cauchy-Green invariant that measured length. And it's actually the sum of the squared sum of the lengths of the sides of the cube. So it doesn't actually express the same thing. It is actually not good enough if it actually wants to express the as rigid as possible energy. So I'm going to advocate and we actually advocated last year in our paper. So don't use the Cauchy-Green invariants. It's a little scary to not use them because they've been around for at least 80 years. But I advocate that you should use ours instead. So from Smith et al 2019. So they can do everything that the Cauchy-Green invariants do. We actually do show that in the paper and they can do more. They can actually express a whole bunch of other energies that the Cauchy-Green invariants couldn't. So these energies from before that I said couldn't be expressed using the old set of invariants. If we use our new set of invariants we can express them just fine. So we can get invariant forms of those energies just fine. All right. And then you can do cool stuff like simulate inflating tires and cars three. All right. So this guy's tires are going to inflate. These were a volume simulation using a wrap in FISD. And you can get a nice looking gag. Okay. So now let's start getting to some other really ugly parts of these simulations. So the force Jacobians. All right. So we have an elastic energy in hand. We need to take its first derivative to convert it to a force. And because we're doing a bare of what can style solve, we also need its second derivative to convert it to a force Jacobian. All right. So again, let's just review what we're looking at here. We have forces on a 10. All right. And we need the forces at each node of this squished element. And the way we got that is we used an F based energy. So when we take the derivative of that, we get a matrix out of it. But we can convert it back to our list of forces using this just this change of variables. Right. Okay. And now we want the other thing. We want the force Jacobian. And this is usually what stands since the students screaming in the room because it is not that easy to derive an implement. Okay. How do we compute this one? So we're going to take the second derivative of our elastic energy. So this is the energy Hessian. So the second derivative with respect to F. Right. So let's notice there's a two here. So we're going to apply that same trick that we did before. We're going to do the change of basis, except we're going to both left and right multiply the change of basis matrix with respect to the energy Hessian. And when we do that, then the matrix that we want, this actually will pop out. Okay. Now, I said before there was an evil professor trick lurking in here. And actually, I just apply the evil professor trick yet again. And again, I'm not going to dig into it at this exact moment. But let's remember it for later. Okay. All right. So I said before that all the energies that you see, they're just a big mix and match game between length, area and volume. So we're just mixing and matching a few invariants here. Right. Now we have this trouble that we have to compute the energy Hessian. So we have to compute the second derivative of this mix and match game that we built for ourselves. It turns out it's just a big mix and match game between the second derivatives of these pieces as well. Right. So let's write out what the second derivative of the energy Hessian of an energy is. So here's a generic form. It's not pretty, right. This is pretty ugly looking. But we can actually break it down. So its individual pieces are a little bit more friendly. All right. So these are the second derivatives of the original invariants. And we can think of these as just preformed Lego pieces. These are some things that we've already worked out. I can just hand these to you. So you don't actually have to work these out yourself. These terms over here, these are gradients. So these are just the first derivatives of each of these invariants. And again, these are Lego pieces. So you don't have to derive these yourselves. We've already done it for you. And it has a generic form for every single energy out there. All right. I don't want to downplay some of the easy parts and the hard parts. So this actually was a bit of a hard part. So how long did it take to actually find these Lego pieces? It actually took quite a bit of effort. Right. So we actually had to find each of these Lego pieces across two papers that were presented at SIGGRAPH. But we have figured these out for you. So I can hand them to you and you should use them. All right. So those were the Lego pieces that I can hand you. That means that what's left over is these six expressions right here. So actually, these are all scalar derivatives. This is just like introductory calculus. So you're just taking the derivative of a scalar function with respect to another scalar function. So it's actually even better, you're in a better situation than initial calculus. You can go ahead and just use some sort of automatic differentiation program. So you can just go to Wolfram Alpha and have it take all the derivatives for you. And once you have those in hand, you actually have the force Jacobian of your energy. So don't feel bad if you have to use automatic differentiation like this. I do it all the time. It's not cheating. This isn't a test. All you want is a robust simulator. And this is totally in balance to do. Okay. All right. So I could give you the code that does this. So everything I've covered up to this point, I could give you MATLAB code that will spit all this out for you. But I'm not going to do it. So what's going to happen is it'll spit out some force Jacobian code for you. And it'll go indefinite at some point. And this will be very frustrating and difficult to debug. What will happen is your simulation will just chug along merely until some mysterious condition is triggered. You're not going to know what. And then the simulation just explodes. All right. So let's go a little further along before we start generating code because we'd like to fix that. All right. So remember from before, we have to keep the entire system positive definite. So how do we guarantee that all the eigenvalues that come out of these force Jacobians are all greater than or equal to zero? All right. So that is the last piece here. Okay. So we have an elastic energy in hand. And we have its forces. We have its force Jacobians. And we want to compute the eigen decomposition of the force Jacobian. So once again, we can mix and match a few generic length area and volume Hessians to get whatever Hessian you want. Now it turns out that you can get analytical expressions for the eigen decomposition of these Hessians as well. And you can mix and match these two. Now you can't do this generally in linear algebra. So if you know the eigen decompositions of two matrices and you just like sum them together, you actually can't say anything about the sum. This is a special case though. So as far as the invariance go, there are some rules that we can exploit to actually get the sum of the eigen decompositions. All right. So let's go through what I mean here. Okay. I said that we can get the eigen decompositions of the invariance. What are these eigen decompositions actually look like? All right. So to see what the structure of the eigen decomposition looks like, finally, finally, we're going to get to this evil professor trick that I've been alluding to this entire time. All right. So what's so evil about this thing that I was talking about? All right. This expression in the middle, this is not actually a matrix. What do I mean? What do you mean it's not a matrix? Okay. So here it pays to recall that the original function was a scalar. And remember we took its first derivative and that first derivative was already a matrix, right? Because f is a matrix, you take the derivative of a scalar function with respect to a matrix, you get a matrix out. All right. Now we actually want the second derivative of this with respect to f. What does that mean? All right. Well, we already took the derivative with respect to a matrix, and if we take this derivative again, what does that even look like? Well, to answer this, let's try to take an easier looking derivative instead. All right. So instead of the second derivative with respect to the entire matrix f, let's take it with respect to one of its scalar entries. All right. So f is zero, just one of its scalar entries. And if we do that, the derivative of a matrix with respect to a scalar, still a matrix. All right. So we'll get a matrix out of it. So we can do the same thing for entry f1. And we can do the same thing for the entry f2. All right. We'll get a matrix over every single time. So we can do this for all nine entries. And in the end, you have this big pile of nine matrices. So what are you supposed to do with them? I have a pile of matrices. Do I arrange them into another matrix? What's going on? All right. So this is a higher order tensor. So this is a fourth order tensor. It is not a matrix anymore. It is a higher order tensor. Now, there's lots of different ways to arrange this. There's no one way to arrange these. Actually, there's no one standard scheme that everyone seems to have converged on. So we're going to adopt this following approach. We're going to arrange them as a matrix of matrices. So it's not like we're stacking this all into one big matrix. Don't do that. The arithmetic actually is a bit different. This is a matrix of matrices. The operations are actually quite a bit different. All right. So let's get rid of all these nice shapes that we had here. And you can see that each of these entries is actually a matrix within itself. And this is actually, this is pretty hard to look at, right? There's lots of scalars here, and it just becomes just even visually confusing. So I said that the operations here are not the same as the traditional matrix operations, right? We can actually convert this back to a second order tensor, the usual matrices that we're used to, using what's called a vectorization operation. So that's what we do here. This is vectorization. So this will flatten things out into a matrix, the thing that you're used to seeing. This is still pretty difficult to look at, right? This is just a lot of information. So what you usually do here is, once you have this matrix in hand, you can just push it to eigen or MATLAB or something like that. And if you want the eigen decomposition, you call whatever eigen decomposition routine it has. Now it will just give you a C of numbers, but at least you'll have something in hand that you can then use in your code somewhere. Okay, that's what you usually do, convert it to a matrix, and then throw your usual matrix tools at it. I'm going to advocate. Don't do this, actually. If you keep it in tensor form, when you look at its eigen vectors, you'll actually see something that's extremely simple and has a very beautiful structure. Now the problem is, is that if you do this, you're not going to get an eigen vector, you're going to get an eigen matrix, right? So what is an eigen matrix? Alright, so this is not a standard form that usually shows up in your usual linear algebra textbooks. In your usual linear algebra textbooks, this is what you usually see, right? So this is your usual eigen value problem. So you have a matrix and you multiply by a vector, and then you get a scaled version of that vector, right? So matrix vector vector, that is the form of your traditional eigenvalue problem. Now if you have a fourth order tensor, this is what the eigen value problem looks like, okay? Instead, you have a tensor on the left, you get an eigen matrix out, and then a scaled version of that eigen matrix shows up on the right-hand side, right? So instead of matrix vector vector, it is tensor matrix matrix. Alright, so why is this important? So if you actually look at this eigen matrix, it has a very simple structure, right? So let's just look at one of these, one of these eigen matrices, so Q0, the first eigen matrix. So in order to see the structure, we're going to have to look back at the SVD of f that we showed you before, alright? So remember this, we looked at that middle matrix sigma before, right? Because if we use that, we could interpret the invariance and what the invariance meant. I never said anything about the other two pieces of this, right? There's a U matrix and there's a V matrix. These are the left and right singular vectors, alright? But if we actually look at the eigen matrix that comes out of the eigen decomposition of the force Jacobian, it actually takes this form, right? So the U matrix and the V matrix from f, it shows up as the left and right matrices here, and then it just sandwiches a constant matrix in the middle. So actually the eigen matrix is just a constant with respect to the deformation gradient f, alright? Now, this actually applies to every single isotropic energy out there. So all those energies I showed you before, this is an eigen vector of all of them. So this is an, again, eigen matrix, eigen vector, these actually are interchangeable using vectorization operation, but it's actually an eigen vector of all of them, right? So that was one eigen vector, one eigen matrix. And here's the first six, alright? So there's actually nine eigen matrices in all. I'm showing you just the first six. But again, these eigen matrices, these are the eigen matrices of every single isotropic energy out there. They all deep down inside have the exact same eigen vectors. So we interpret these a little bit more in chapter seven of the notes if you're interested. Now I told you there's nine of these. If you're lucky, then this is what the last three look like, alright? And actually, you get lucky a lot. So because it's really nice when this form appears, we said, okay, so you hit the eigen value jackpot if it turns out that all nine of these work out to this. So sometimes you get slightly less lucky. So they don't separate as cleanly. And instead, you get a linear combination of these. So it's still a pretty good closed form solution to this. But it's not as nice as those three matrices you see on the bottom here. Okay, so that is the eigen matrices, which you can actually vectorize into the eigen vectors. What about the eigen values? These are the ones that we actually cared about the most, right? Because we want to know when they're less than zero. Okay, so there are expressions to get the closed form eigen values. Now they're not that pretty, right? Now the good news here is I did type these all into a MATLAB script for you. So you don't have to transcribe these yourselves and figure out, you know, did you transcribe them correctly or are you actually getting the wrong values out of them? And you can go ahead and take any energy that you want and just get the analytic eigen values for them. It's pretty fast actually. I'm just going to show you a run that I did myself. All right, it's actually loading up SMPI. So it would actually be slightly faster than if you did it yourself. It loaded up the hyperlastic energy. And here's the eigen values. This is the ARAP energy, right? So you hit the jackpot actually. All the eigen matrices cleanly separate. And this is it. These are the eigen values. All right, so we actually have the closed form eigen values for the ARAP energy, so the eigen vectors and the eigen values. All right, so the eigen vectors themselves, these are listed in chapter 7.3 of the notes. God, here they are again. Okay, so you can actually run this for whatever energy your heart desires. So in section 7.3.3, I give you the generic script to do this. And I actually go through and crank a bunch of energies through it. So St. Kanoff, Kirchhoff, Stable Neohokian, the corrotational model, symmetric Dirichlet and symmetric ARAP. We give you the analytic eigen systems for all of these. So if you ever wanted these energies, you can just look them up right now. All right, so let's put it all together. So let's say that you have your own energy. So as I said before, all these energies, it's just a big mix and match game between the length, the area, and the volume. So if you decided, well, I care about length, area, and volume in a slightly different way than everyone else has in the past, so I'm going to design my own energy. Okay, well, that's great. So how do I go ahead and compute the force radiance? Well, one thing you can do is you can plug your amazing new energy into line 11 of figure 7.7 from the notes, and that'll crank through. It won't take that long as we just saw, and it'll give you the eigenvalues for the eigen decomposition of the energy hessian. All right, so these other matrices here, the eigenvectors, these are just Lego blocks. So again, I give you these, give these to you, and then you can just combine everything together to get the full energy hessian. All right, so this term here, got it now. All right, so this still is not a force Jacobian within itself, so we're going to left and right multiply this by the DFTX matrices here. These are Lego blocks too, actually, and I give you the code for this. This is in figure E.1 of the notes. This is even in C++, so you can actually use this in fast code, so not just in MATLAB code. All right, and once you have that, you're all done. You've got your force gradients now, so that thing that was super painful for students taking their first simulation class, hopefully it's slightly less painful to do now. All right, so in the actual Pixar-FISTI implementation, we use this version of ARAP, the one written in terms of our invariance, and we invoke this when Poisson's ratio has been set to zero. So when the Poisson's ratio is non-zero, then we use stable neo-hookian energy, which is the one that we presented two years ago. All right, so there's a bunch of left out, but again, there are these course notes out there, so if you're interested to fill in some more gaps, you can go ahead and read those notes. All right, and if you're really keen, even after having read those notes to go even further, come and apply to YLCS, come work with me, and we'll work on problems like this, and it'll be tons of fun. All right, and with that, that wraps up my part, and I'm going to hand you off to Dave. Welcome to the second part of Dynamic Deformables. My name is Dave Everly, and I don't smile for photos. I'm also the tool sim lead at Pixar Animation Studios. For the next hour, I'll be discussing some of the implementation details behind our in-house simulator FISTI, which is used to simulate cloth, skin, and volumes. FISTI has evolved over 20 years at Pixar, and has been used in over 19 films. It was first used on Monster's Inc., and is still largely based on the paper Large Steps in Cloth Simulation. Here is a brief summary of what I'll attempt to cover. I'll need to gloss over some details in the talk, but on most of the topics, there are more details in the course notes. Let's get started. If you've not computed forced Jacobians before, this section won't make much sense. Skip it and come back after you've had some experience. It has some useful information that is not covered in our course notes. So I'm going to try something a little different to keep your attention. If it's too much for you to handle, skip ahead and you'll be in more familiar territory. In the first part of this talk, Ted showed us how to compute the internal forces from the deformation gradient which he denoted uppercase F. He chose to denote the forces as psi, but I'm going to break with that convention just in this section and denote the forces acting on node i as f sub i. Let's consider one of our constitutive models acting on a 3D tetrahedral element. Here we have our 3 by 1 forces stacked into a 12 by 1 column vector. We need to compute the Jacobian of F with respect to our position state x. This yields our symmetric 12 by 12 Jacobian matrix of 3 by 3 blocks. Hungry yet? To go further, let's rewind and revisit a basic property of the forces. Our internal forces must sum to zero. This means that our block columns must also sum to zero. We know the Jacobian is symmetric, so this means that our block rows must also sum to zero. Thanks to Tim Haynes, the smart monkey first introduced me to this property at Disney's secret lab. He used it with respect to the global system Jacobians, which I'll talk about later. Terms from boundary conditions, which are also discussed in the notes, make it less applicable to our global system blocks, but there is still quite a lot we can do with this knowledge at the individual element level. If we are computing these terms in the blocks by hand, we only bother computing the upper triangular blocks. With our new found insights, we can write some of these difficult blocks in terms of others. Let's say the partial of f2 with respect to x3 was something we wished to avoid, but we had the second block from the top row of f1 and had the partial of f2 with respect to x2 and the partial of f2 with respect to x4. We can compute the partial of f2 with respect to x3 as follows. We can ultimately play this game for other block terms in the rows further down, saving us from unnecessary work. Let's examine another use of our new found knowledge. Under ideal conditions, our forces sum to zero, but numerical roundoff prevents us from actually being the case. We can easily enforce this property, though. The blocks present us with a similar accuracy problem that must be handled differently. Even when forming the full 12 by 12 Jacobian by tensor operations, our Jacobians do not come out to be 100% symmetric due to numerical roundoff issues. We first average the values of the off diagonal blocks and assign them to the upper diagonal block. We then assign the transpose of the result to the corresponding lower diagonal block. We then assign all of the diagonal blocks to be the negative sum of the row blocks. The row sums and column sums won't be exactly zero, but this tends to make the errors smaller. I have more hard earned advice to give on computing Jacobians, so don't let the cheesy graphics bore you. Just listen, you're not going to regret it. The payoff is huge. Let's say you compute the terms by hand, verify them with an analytical math package, and then verify that you've made absolutely no mistakes trying and scribing the expressions into code. I bet you'd be pretty proud of yourself. I know why I was, but it's a trap! The right math is wrong. During a simulation, there are times when the element Jacobian will not be semi-definite. This means that the left hand side of your global sparse matrix for methods like backward whaler may become non-invertible. Even the initial baref wick and cloth implementation in FISTE suffered from this and forced them to use an adaptive time stepping scheme. The symptoms can be mysterious, ugly, and a nightmare to debug, unless you take my advice. You can go through the math of analyzing the eigenstructure of your tensors and take the minimal appropriate action, or you'll need to adjust the math in some ad hoc way. In practice, this usually means discarding certain terms entirely or discarding them under the conditions that might cause the indefiniteness. My advice to anyone who isn't absolutely confident their Jacobians are always well behaved is to have an alternate running mode in their simulator that does an eigen decomposition of the local element matrices to uncover these issues. Sometimes they are there but small enough that you won't notice, but other times they can wreak havoc. If you're saying, there's no way I'm taking any advice from a guy who illustrates matrix math with cheese, I reluctantly respect that, but you have been warned. Relax, the rest of this presentation will be more typical of SIGGRAPH and not involve through products. We will now move on to discussing the driving purpose of computing their forced Jacobians, implicit integration with backward Euler. Large steps in cloth simulation ushered in a new era in deformable simulation to computer graphics in 1998. It clearly demonstrated the advantages of implicit integration. I'll start by giving an overview of the most simple first order accurate method we have at our disposal, backward Euler. In backward Euler, we want to update the velocity with delta v and our delta x depends on this updated velocity. Delta v however is determined by Newton's second law with the forces evaluated at our unknown state at the end of the time step. How can we do this? To approximate the force at the future state, we use a first order Taylor series approximation. This requires us to have the Jacobian of our force with respect to both the position and velocity state. The Jacobian with respect to velocity state is needed to account for any damping forces in our system, although we did not explicitly discuss damping forces in the previous sections for simulating volumes we use really damping and for cloth our damping model follows the one described by the original paper. Using the Taylor series and some algebra we arrive at the linear system we must assemble and solve at each time step denoted H. In practice the linear system is very sparse and an iterative solver like conjugate gradient converges to a solution accurate enough and fast enough for our purposes. It is possible to implement backward Euler without assembling the global matrix. Instead one can iterate over the elements and accumulate their effect on the fly in the solver, but in practice this is much slower than traversing the entries of a global matrix during the sparse matrix vector multiply. Using a global sparse system representation allows for a more coherent traversal of the memory. Given the current particle state we can compute the desired delta v necessary to achieve a target velocity or position. We are not limited to constraining all three degrees of freedom or doffs however. Keeping with the notation of the original paper we denote the vector of the desired delta v projected onto the constrained degrees of freedom as z. To maintain the z values during the solution of our linear system the original paper introduced the filtered precondition conjugate gradient method. The solution is initialized to z and the values are enforced by applying per particle constraint filters s to the variables that can change the solution during the conjugate gradient iterations. The definition of s allows us to constrain the motion in any combination of the three degrees of freedom which can be described by an arbitrary frame. Here p hat and q hat are unit vectors that are orthogonal to each other. While this gets the job done the application of these filters at each iteration impacts the rate of solver convergence and impedes performance. In 2015 Jones and Tamstorff presented the pre-filtered precondition conjugate gradient technique as a part of the Disney paper smooth aggregation multigrid for cloth simulation. ppcg is a way to pre-filter the system once before passing it to the conjugate gradient solver. They transform the original system by first defining the vectors c and y from the original system variables. The constrained system can then be expressed in the following form. Once they are done solving for y they transform the solution back to x. The transformed system leads to a substantial reduction in iterations required by the solver and also outperform the modified pcg method of Asher and Boxerman. The ppcg method led to an average solver increase of 50% over a set of examples depicted in the course notes. Let's examine the linear system assembly in terms of this simple 5x5 cloth mesh. There are per triangle stretch and shear forces, there are dihedral edge bending forces, and there are per particle forces. Example of per particle forces might include gold springs, viscous drag, aerodynamic drag, etc. Some per particle forces like gravity are considered constant and have a zero Jacobian. In FISTI, the connectivity of our cloth and volume meshes remains fixed, and we'll be discussing how we exploit that. Each of the force or mesh component domains can be viewed as fixed elements in our system. Together these fixed elements will contribute their three by three blocks into a fixed global sparse matrix structure, which we will refer to as the FGSM. We don't assemble the entries of A directly in the FGSM, but rather just the global DFDV and DFDX terms needed to construct A. Each of our fixed elements has storage for its particle forces and the related three by three blocks of the Jacobian's DFDV and DFDX. Let's examine the block storage for a dihedral element that acts on four particles. The blue blocks here represent the storage for both DFDV and DFDX terms. Since the Jacobian is symmetric, we only store the upper triangular blocks. The numbers in each box represent the block entries with respect to the element's local particle indices. Triangle and per particle elements would naturally have smaller storage requirements. Since the FGSM is fixed, it can be flattened into a simple contiguous vector once at the beginning of our simulation. Each global block in the FGSM has pointers to all the fixed element blocks that contribute to it. There is a similar global accumulation structure in pointers for the forces. We perform the assembly of the global DFDV and DFDX terms with a two-pass strategy. First, we do a parallel pass over all the fixed elements in our system to compute the forces and their Jacobians, which again are locally stored. Then we do a second pass over the vector blocks of the FGSM to accumulate the contributions from the elements. The fixed forces are also accumulated from the elements in a straightforward fashion as a part of this two-pass algorithm. While both steps are parallel, a lot of random memory access is occurring, which impedes the performance and scalability. But before we attempt to address that, we have an important detail of the assembly to discuss. As I mentioned earlier, the elements have a local particle ordering. Every particle in our system also has a unique global index that indicates its row in the global system. As the block entries from the elements are accumulated into the FGSM blocks, the global indices must be examined. If the global block pair indices of the FGSM block have a different ordering than the global indices of the local element block, we must first transpose the block before adding it. Some forced Jacobians have symmetric blocks, and this will not matter, but others do not. Making a mistake here can lead to bizarre simulation behavior and hours of nightmarish debugging. You have been warned. Instead of having a variety of fixed elements to support particle configurations over different domains, we define a single fixed force element that has storage for forces and Jacobians of up to four particles. To understand why, let's examine the number of elements in associated blocks that need to be dealt with if there is a separate element type for each domain. In our 5x5 example, we would have 32 stretched shear elements with six blocks, 40 dihedral elements with 10 blocks, and there would be 25 per particle blocks. There's a multiple of two here because we have blocks for both DFDV and DFDX. This means we would have to traverse 97 elements in our first pass and gather from 1,266 blocks in the second. To improve upon this, I assign the fixed force element some unique set of the domain elements. Naturally, some elements will have more work than others, but let's compare the numbers. The number of elements traversed in our first pass is cut by half, and the number of blocks accessed in the second pass is cut by almost a third. Why does this matter? Modern CPUs and GPUs can do amazing amounts of number crunching. Accessing the data and feeding the process efficiently is what matters. This is the general theme of the majority of optimization efforts we will discuss. I don't recommend this aggregated element approach for a first implementation. If you do implement it, remember the bookkeeping necessary to perform the local to global transposition needs to also partially be handled at the element level and not completely during the accumulation pass. Up to now, I've told you that our linear system structure is fixed. I am altering the deal. FIST's linear system does not have a completely fixed structure. The majority of the block entries are fixed, but there are a number of forces that come and go during the course of the simulation which change the sparsity pattern. Let's examine one case of these transient elements and discuss how we handle them. Contact response between deformable objects and self-collision contact demand a transient representation in our system. I'll discuss the details further in a later section, but for now, it is only important to know that our response generates implicit forces and block entries in our linear system. Other sources of transient elements are user forces used to control the simulation. Here, the pink particles represent the set from a hypothetical vertex face contact. Transient elements do not have local storage for their forces or Jacobian blocks. Instead, they directly put their contributions into global structures that have storage per thread to avoid locks. For the matrix entries, we creatively name this structure the transient global sparse matrix or TGSM. Our storage per thread is a flap map per matrix row. After all the transient elements have been traversed, a reduction step is performed on each row. This reduction is performed in parallel over the rows. If we overlay the block sparsity pattern of the FGSM and TGSM, we can see that there are many blocks that exist in both. In parallel row fashion, these overlapping blocks are removed from the TGSM and have their contribution added to the FGSM. This leaves us with two global structures containing the accumulated DFDV and DFDX blocks needed to construct our A matrix from before. Note, I've added some more pink blocks to represent other entries from hypothetical transient forces. At this point, we have the global DFDV and DFDX matrices represented by the FGSM and TGSM structures, along with the vector F of all the accumulated forces. With this, in our current velocity state, we can construct the right hand side B of our original linear system. Any desired constraint values reside in our vector Z. With our per-particle masses and constraint filters, we are able to construct both sides of the pre-filtered linear system simultaneously in row parallel fashion. The combined block terms on the left hand side for both the FGSM and TGSM must be transformed by any of the necessary S filters. We refer to our pre-filtered left hand side as A prime and the associated right hand side as B prime. It is important to clarify that we have not constructed a single matrix representation of our left hand side, but rather a split representation of it. When we are constructing the filtered blocks of A prime, we store the values in two separate block compressor storage, or BCRS structures. BCRS is a common and simple sparse matrix format, and for brevity, I'll omit discussing it here. These BCRS structures can be populated very efficiently and almost entirely in parallel. When we encounter a diagonal block, which will always be in our fixed pattern portion, we invert it and store it for our block diagonal Jacobi preconditioner. We have found this simple preconditioner works best in practice. Details on populating the BCRS structures and the results of some experience with larger block preconditioners are in the course notes. For a first implementation, I would not suggest keeping the left hand side split and writing a custom solver that accepts the split matrix. You can easily combine the two structures and put them in whatever storage format, a numerical library of your choice requires. We are now done talking about assembly and can focus on improving solver performance. So far, our computation has been in double precision, but when we populate the BCRS structures, we convert to single precision to reduce memory bandwidth during our sparse matrix vector multiplies. Here is another common technique to improve solver performance. This is the sparsity pattern of a tetrahedral mesh output from some geometry library. An iterative solver will perform many sparse matrix vector multiplies. While the block values in the BCRS structures are in continuous memory, the access to the vector during multiplication is not. We can gain some performance if this random access is reduced. At the beginning of our simulation, we examine the structure of the FGSM and perform bandwidth reduction, which reorders the global indices of our particles. I found that the boost grafts libraries reversed cutthill mickey ordering works well, but other orderings are also available in the library. Applying this reordering yielded a 5 to 10 percent improvement on our solver. Other solver optimizations are discussed in the course notes. You're probably wondering, what was the payoff of all of this effort? This is a comparison of our performance in early 2015 versus 2016. Initial times per frame were before the addition of continuous collision detection, and later times include its additional overhead. The original bare-footed implementation had locking in the system assembly. We observed around a 3x gain on our system assembly. By using PPCG, bandwidth reduction, converting our solver to single precision, and a handful of optimizations in the notes, we achieved a gain of around 4x on the solver performance. Our largest cloth sims on cocoa were the dresses of these dancers, each with 108,000 simulated vertices. We ran these at 20 steps per frame for better energy conservation with a second order accurate integration scheme. Each had an average simulation time of 20 seconds per frame with 8 threads. Even though great gains were made with respect to the initial bare-footed implementation, we know a lot of optimization and performance gain potential still remains in FISTE. Let's examine for a moment where we are. You've got some constitutive models. You've got system assembly and constraints covered. You just need to write or use a linear system solver and you've got a full simulator in hand, right? Without collision detection and response, a simulator would not be of any use to production or impress any of your simulation friends. Yes, simulation programmers have friends. Anyway, there are a wide variety of collision detection techniques, but here we will only focus on the approaches used in FISTE. FISTE's collision surfaces are represented by manifold triangle meshes. We use bounding volume hierarchies of axis-aligned bounding boxes. Again, we refer to the non-simulated geometry as kinematic. The bowl here is an example of a static kinematic object. If you've tuned out from looking at the 20-meter-long hypno-ribbon, wake up! Pixar's Universal Scene Description has a much broader utility and graphics, but I highly recommend it to debug simulations. Create an option for your simulator to export subframe collision data and geometry. Examining the subframe data in USDVU can reveal a lot about the behavior that would otherwise go unnoticed at whole frames or just with unit tests. In fact, I wrote some additional visualization code for this course, which led to the discovery and fixing of some ancient bugs. In this example, the colored lines are proximity contacts, which we will discuss next. Proximity detection is one of many categories of collision detection techniques. Given the geometry state at a discrete snapshot in time, its purpose is to examine which feature pairs are close enough to be considered in contact. The contact criteria is typically governed by user-defined thresholds. If the mesh features are deemed sufficiently close, reapply a response to attempt to prevent them from actually colliding over the time step. FIST's kinematic meshes are oriented in mostly closed surfaces. We only perform particle-face proximity queries between dynamic and kinematic objects. Let's examine a profile of a kinematic surface here shown in blue. First, we define edge bisector planes for each edge on the mesh. The details of the construction are in the notes. We have a user-defined inset and offset parameters to define our collision envelope. Using this information, we create clip cells that partition the space above and below the surface. Belief-level bounding box of our BVH encloses the cell region for each face. This concept of partitioning is very useful for ensuring that we don't get unwanted multiple contacts for a single particle. Multiple contacts per particle are certainly still possible, though, but I'll get into that case later. Our first step in the proximity query is to determine whether a particle is in the clipped cell of a given face. This can be done by simple distance queries to the plane of our face and its edge bisector planes as described in the notes. If it is outside our cell, we ignore it. If it is inside the cell, our work is not done. When our surface becomes highly creased as shown here, it is possible for a particle to be inside the cell, but outside our collision envelope. Once the point is established to be in the cell, we perform a closest point to triangle query, which gives us the barycentric coordinates of the closest point and the normal on the surface. Let's discuss our response for the dynamic particle kinematic face contact. We first compute the face relative coordinates of the particle with respect to the face at the current time. Using the kinematic vertex positions from the face at the end of the time step, we can compute the target position of our particle. Using the current position and velocity, we trivially compute the necessary delta V required to reach this future relative position. We then guarantee that the delta V in the collision normal direction at the end of the time step is achieved by adding the constraint filter S A to our PPCG system assembly. So far, we've just enforced that the particle will match the velocity of the kinematic face in the n-normal direction. This allows the particle to slide over the face in the other two under-constrained degrees of freedom. We can be sure that with this hard constraint mechanism that the vertex will not collide with the face during the time step. If we want to achieve any desired offset distance at the end of the time step, we use the position alteration technique discussed in the notes and the original barefooted paper. The constraint release criteria is also discussed there as well. Let's now talk about the scenario that can lead to multiple contacts per particle. Finding multiple faces per particle is most commonly the result of pinching situations between kinematic surfaces. The fly-papering algorithm was introduced in the paper untangling cloth to handle this case. We still use a variant of that technique described in the notes, but have adopted a different approach on recent films to reduce its need. We use FIS-T's volume simulation on a tetrahedral mesh to perform a deep penetration pass on our characters. The subdivision surface that feeds the collision geometry for our cloth simulation is first deformed by the surface of this volume sim. What you see here is 100% out of the box. The only obvious artifact is the temporary snag on barley slum towards the end. These volume simulations typically take 2-3 seconds per frame with 6-8 cores. This process saves our simulation TDs from tedious hours of manual sculpting. Further details on the technique can be found in our 2017 SIGGRAPH talk. Besides a robust and efficient constitutive model, this process heavily relies on FIS-T's collision handling between dynamic surfaces. FIS-T performs vertex face and edge-edge proximity queries between dynamic surfaces, including self-collision. Because meshes like cloth can be opened, there isn't a notion of an inset and offset, but instead a self-offset and offset, which define envelopes that are both equal and opposite on both sides of the face. Our previous particle and cell check is not a sufficient rejection criteria for dynamic surfaces. For example, this particle is in our collision envelope but outside our definition of a face cell because of the boundary edge. The necessary changes are discussed in the notes, along with a contact visibility criteria between edges within the collision envelope. Before discussing our response, let's cover some conventions for a dynamic-to-dynamic vertex face proximity contact. We denote the two points of contact as PA and PB, and we have the convention that the collision normal points from B to A. The positions on each side of the contact are computed using the particle positions and their associated barycentric coordinates. We also define sine barycentric weights for the particles involved in our contact, according to which side they reside on. For the vertex face case, our sine weights are defined by the vector W. This makes computing the relative position and relative velocity of our contact configuration a simple weighted sum. Our collision normal is just the unit vector of the relative position. Using the length of the relative position, a desired offset, the relative velocity, and our collision normal, along with some spring coefficients, we can compute a damped penalty spring force that acts on PA. The equal and opposite force acts on our contact point PB. Penalty forces are frowned upon in the literature. Stability is not a concern, though, because we're using implicit integration. Although penalty forces cannot make guarantees to prevent an impending collision, they are computationally cheap to the alternatives mentioned in the notes and are quite effective in practice. A simple application of our sine weights allows us to distribute the response force to the particles. This comes from the chain rule and the fact that the response forces at PA and PB are equal and opposite. Our spatial Jacobian block is given by the following equation, but do not blindly insert this into your code. Do you remember in cheesy and important tips on Jacobians I mentioned that you should always check that your Jacobians remain well behaved? Well, because we only apply a contact penalty force when the length of delta x is less than the desired offset, we have a spring that is in compression. The second term of the equation can make our linear system become indefinite. Besides causing trouble, this term filters motion orthogonal to the direction of our response, so we discard it without much concern. The damping Jacobian has a similar form. Everything we've discussed so far easily generalizes to response on edge-edge proximity contacts. But we still have not discussed how to distribute the Jacobian block from our point of contact PA to the blocks relating to the particle pairs in our contact configuration. This turns out to be a product of our sine particle weights. The same relationship holds for the damping Jacobian, but it generalizes much further. It is now possible to define a virtual particle anywhere on a mesh via barycentric weighting and have an arbitrary force applied to that point. You have all the necessary techniques to distribute the forces in Jacobians to the particles of your system. Think of the cool things that you can construct with this power. Fisty's performance demands, unfortunately, rule out the proper handling of friction. I discuss one of the tangential forces used to fake friction between dynamic surfaces in the course notes, but I won't cover it here. While such ad hoc approaches have some drawbacks, they can be more intuitive for TDs to use rather than relying on an accurate friction model to address a director's notes. What do we do to account for interaction of arbitrary force connections between simulated and kinematic particles? One could add rows to our linear system and consider these kinematic particles as fully constrained, but that would be an expensive way to handle things. For kinematic meshes, we can compute the velocity and delta v from the known position using finite differences. We label these with the subscript k. The forces on the simulated particles and Jacobian blocks representing the interactions between them go into our global system just as they did before. The off-diagonal Jacobian blocks for dynamic kinematic interaction are accounted for by adding the appropriate terms to the right hand side. First, we add the delta v term from the left hand side. Then, we add the velocity term. This accounts for the interaction without increasing the size of our linear system. With what we covered up to now, you can have arbitrary forces acting between any desired points on the meshes in your simulator. In FISTI, we don't use boundary conditions explicitly to resolve dynamic kinematic collisions, but when large parts of our dynamic meshes are fully constrained, they are essentially kinematic. We remove these particles from the system to reduce the linear system size. By using boundary conditions, we can have the dynamic surfaces collide with these removed regions and use the same penalty force model we described earlier and get the correct behavior. With that said, we are still very far from done discussing collisions. Proximity queries and their response will only get you so far in handling collisions. They can only anticipate an attempt to prevent potential collisions from a discrete snapshot of the geometry state. Fast moving and thin geometry will result in collisions being missed due to temporal aliasing, unless a ridiculous number of steps per frame are taken. This would slow down our simulation. One of the primary goals and advantages of using implicit integration is to take larger and fewer steps. Continuous collision detection, or CCD, is what we need to handle these scenarios efficiently. CCD was added to FISTI in 2015 to deal with the fast moving thin geometry of our skeleton characters on cocoa. In FISTI, CCD is performed between all dynamic mesh pairs, including self-collision and all dynamic kinematic mesh pairs. CCD is performed after we solve our linear system and compute a candidate position state using the updated velocity. Each of our mesh vertex positions, X sub i has a current start position and a candidate end position. These are used to define a linear trajectory over the time step. Let's examine the vertex face continuous collision detection test. We can define the normal of the triangle as a quadratic function of t from its linear vertex trajectories. To examine whether a particle ever enters the plane of the triangle in the interval, we want to know if the projection of its trajectory minus that of a point on the face will be zero. Said another way, we want to determine if the volume of the tetrahedra between the four points ever goes to zero in the interval. This gives us a cubic equation in T for which we must find all roots in the interval zero to one. We traverse the roots in order. At each root, we must perform a geometric test to determine whether the point at the root is inside the triangle at that time. Looks like this one missed. If a collision had been found, we'd stop. Otherwise, we examine the rest of the roots until they are all exhausted. The same process can be performed for edge-edge collisions and is covered in the notes. Easy, right? Cubics have an analytical solution, so you might as well consider this as good as done. Nope, not so fast. Numerical roundoff issues make the straightforward process difficult to implement robustly in practice. Robust and efficient CCD has been the topic of many publications. These papers represent a small sample. In Fishti, we take a less rigorous approach than what the literature offers, partly for performance and also because we don't believe that the assumptions on which the guarantees these papers offer holds up in a production environment. I'll have more to say about that later. First, we solve for our roots with analytical methods. I recommend reading Jim Blinn's series on solving cubic equations. We then polish the roots using an implementation of Schroder's method in the Boost Math library. Our routines pass millions of procedurally generated trajectories at a variety of scales within typical use. With CCD and our response, we are then able to handle academic scenarios such as this funnel test. We see in the next movie that it remains intersection-free, and if we disable CCD, we show the resulting intersections for our comparison. We still need to discuss CCD response. Our formulation closely follows the notation of the paper, robust treatment of simultaneous collisions. Their paper focuses on simultaneous contact, but their notation also provides a compact way to discuss the response for a single collision. The paper derives things in terms of velocities, momentum, and impulses. I deviate from this and operate on displacement directly. After our entire process has finished with all collisions, we make the velocities consistent with the displacements over the time step. Let's again examine the vertex face collision case. Each vertex has a displacement defined by its endpoints. We then define the configuration state vector q to be the vector of all displacements involved in the collision we wish to resolve. We then define the scalar function c of q that determines the relative normal displacement between the points of contact on each side of the collision. NC hat is our normal at the time of impact, or TUI, found by the detection. We orient the collision normal such that c of q is less than zero because we know a collision has occurred in the interval. We then take the gradient of c with respect to q. Continuing their derivation and skipping over the details, we arrive at this equation. The paper's formulation seeks the global minimum change in the vertex momentum that results in the relative normal displacement becoming zero. Again, I'm using displacements instead of velocity, but the desired effect will be the same. After solving for lambda, we update the end positions of our particles. If any particles were constrained and we wish to have that constrained motion respected by our response, the inverse mass can be filtered by the appropriate per-particle constrained filters as from the previous section. Particles of our kinematic geometry are always considered to be fully constrained. One needs to be careful including these constraints into the response, however. There are configurations that will lead to undesirable large deformations or division by zero. Handle these situations however you see fit. You have been warned. As I mentioned before, the response will stop all relative motion in the contact normal direction. If the points of the contact started below our offset threshold, this is exactly what we want. However, if the initial relative distance is greater than our offset, we'd like to halt the motion at the offset. The following modification to the right-hand side of our equation achieves this and is discussed in the notes. I should also mention that CCD response is completely decoupled from the system dynamics and can lead to surface stretching artifacts in the extreme cases. Any deformation required to keep the surfaces intersection free will be encountered at the next time step. So are we done with CCD yet? Almost, but not quite. Let's consider a short stacking ribbon such that the bottom particles are colliding against the ground plane and not moving over the time step. In a single time step, a lot of continuous collisions can occur. Also, each of the vertices in our mesh can be a part of multiple collisions. How do we sort out this mess? Well, sorting is definitely part of the answer and a justification of why we want accurate time of impacts for our collisions. We must process the responses in chronological order, but unfortunately it's more tricky than just doing that. First, some collisions might not occur once earlier ones are resolved, so before resolving a collision we need to see if it still exists between the feature pairs. Second, by responding to a collision, we might cause a new collision to occur that did not previously exist, so we need to run CCD detection between applications of response. In the notes, I describe in pseudocode how to handle this and how to make the collision detection passes more efficient after the initial one. In FISD, the response iterations are fixed to a maximum number and then we just give up. We don't employ the fail-safe methods of rigid impact zones or inelastic projection from the literature to guarantee against failures. Instead of incurring the additional cost and driving ourselves mad trying to enforce an intersection-free guarantee under assumptions that don't hold, we accept that failure is often necessary in production work. Let's see some examples. Funnel and cloth ribbon tests are just the beginning to making a production-capable simulator. They exercise things when the inputs are clean and everything should go according to plan. Production is messy. Here, Miguel's legs are intersecting with the crypt and the guitar intersects with his body. Different parts of a character's costume may also initially intersect. The CCD literature assumes that we start in a clean state, but the production reality is that we often don't. It would be a pain for our TDs to guarantee a pristine state, and even then, a shot might require that things be allowed to initially intersect. How can we allow for pinching between kinematic objects and recovery for cloth intersections and still have robust CCD? Why doesn't our CCD in response prevent the cloth from entangling here? The answer is global intersection analysis, which was introduced in the paper, Untangling Cloth. We have two intersecting surfaces depicted here. The algorithm first finds all the edge-face intersections between the meshes and constructs closed-contour lines between the meshes. The algorithm first finds all the edge-face intersections between the meshes and constructs closed-contours from the data. Then, a flood fill is performed on both sides of the contour. It is assumed that the flood fill with the minimum area is the region of intersection, which is generally true in practice. In FISTI, we perform global intersection analysis, or GIA, between all kinematic mesh pairs, all dynamic mesh pairs, and all dynamic kinematic mesh pairs. Why will become clear soon? We are now going to discuss how incorporating information from global intersection analysis into our proximity and CCD response allows us to handle a number of common but challenging scenarios. The information from GIA allows FISTI to fail when it is necessary and later recover gracefully when it is possible. Here, we have two intersecting dynamic surfaces with the GIA intersection regions colored in blue. We will view the recovery at subframe steps in USD view. The GIA information serves two purposes in this example. During our vertex-face penalty response to the proximity query, we examine whether the vertex-face pair was labeled as intersecting and belonging to the same intersection. If they are, we modify the response to pull the particle to the opposite side of the face instead of pushing it away. Second, when continuous collision examines the mesh feature pair, we look to see if it was labeled as being an intersection prior to taking the step. If it was, we can early out from the CCD test. Once the GIA region disappears, proximity response and CCD can work robustly to prevent the intersection from recurring without the need for special heuristics. I initially tried many, but all of them could lead to collision failures as well. Using the GIA information provides us with the desired one-way pass-through behavior. In FISTI, GIA is also computed between kinematic and dynamic mesh pairs to help improve their proximity response. If a particle is too deep inside our kinematic object, as depicted on the left, our proximity response will push it out to the wrong side of the object. Instead, we only allow a particle beneath the surface to escape if the face has been colored by GIA. Instead, particles near the intersection are pushed out first and then pull out the ones deeper inside. This was needed for the thin skeleton characters on Coco. In red, we see an initial deep intersection in the red box. The two areas you see in green are actually modeled holes in the simulated garment. FISTI is able to properly resolve the situation in the first frame of the simulation. We also perform GIA on and between kinematic meshes. In this example, we have omitted our volumes and passed entirely. We perform GIA on the geometry at the end of the time step and cache the information from the previous one. This way, we can know locally if an intersection is persisting, occurring, or vanishing over the time step. Before our body simulation, this was crucial to make sure our CCD response didn't always max out the iteration limit from particles ricocheting between pinched kinematic regions. Even with our body simulation, this is still important today. Characters collide with a variety of geometry in the scene that we do not simulate with volume meshes. Also, there can sometimes be a discrepancy between the TEP mesh and collision subdivision surface resolutions of our body sims. This causes small pinched regions in our cloth simulations to occur sometimes, but the GIA information allows us to handle these situations. Okay, I just want to make this point one last time. With global intersection analysis information, you can make superior response decisions for proximity queries and continuous collisions. Don't spend months of development writing algorithms that can theoretically never fail. Collision failure is often necessary in production work. In my opinion, one should accept that and instead focus on making a system that handles these failures gracefully. Let's recap what we've covered in part two of our course. We discussed the common mistake people make of assuming that the correct math will provide the correct behavior in the simulator when implementing forced Jacobians for implicit methods. We also saw some tricks to reduce the hand crank derivations needed and how to make the Jacobian numerically more consistent with its properties. We covered backward Euler integration and the efficient PPCG method for constraining any of our particle's degrees of freedom. Efficient linear system assembly and tips on speeding up the solver were also covered. We then looked at proximity and response for dynamic kinematic object pairs and dynamic object pairs. The discussion of two-way coupled response generalized to allow for arbitrary connection points between simulated meshes. We then extended this generalization with boundary conditions to properly handle interactions at arbitrary connections between dynamic and kinematic meshes. We then solved our temporal aliasing issues with continuous collision detection and response. We then discussed global intersection analysis and finally we illustrated how GIA can be used to enhance our response to proximity and continuous collisions. That was a lot to cover so please refer to the course notes for more details. We hope you've enjoyed the course or at the very least have found a reliable cure for insomnia. I have anticipated some questions and have prepared answers for them if they are asked otherwise please feel free to ask something else. An updated version of the course notes is available at these locations. Since the upload to SIGGRAPH we've added some content and have been actively making corrections. I also want to thank the artist at Pixar for making this work fun and interesting. In particular though I'd like to thank Chris Campbell for helping me with an issue in the making of the cheese drop. Enjoy the rest of SIGGRAPH and we'll see you at the Q&A.