TOP Webinar Tutorials: Level set methods 1

Added: 2026-05-15

[00:00:24]running.

[00:01:24]I'm not going Okay, it is time. But I'll just give it a bit more time. People are still coming in. So Okay, I think it's time to start. So, um, welcome and thank you so much for being here. Um so this is the kickoff of the uh the second lecture series uh as part of the top webinar and the um the topic this time is on uh the level set method or level set methods for topology optimization and um the recording is currently live on YouTube and it will be available there later as well. And um as always, you can follow along uh on our newest uh plans on the website and register for the mailing list there as well. And um I'm currently moving trying to move to a more automized mailing list system. Um I see a hand up.

[00:04:42]Is the sound okay and everything?

[00:04:47]What is the hand?

[00:04:50]Oh, okay. More hands.

[00:04:55]Yes. Yes. Very clear. Okay. So, I don't know what the hands mean. Um, so I'll just continue. So, today we have the esteemed professor Kra and uh also senior researcher Charles Taponyi. Um so I think Grego needs no introduction and if I were to do it justice uh just like Ola we would be here all day. Um so I'll give a short one. Um so Gregoire has been one of the influential mathematical voices in topology optimization for for several decades and uh his work began in the late 80s with homogenization theory and he was also among the sort of pioneers who brought topology optimization into the mathematical sort of mainstream you could say and uh he is perhaps best known for his contributions to the level set method for structural optimization and uh he has also uh continued working on bridging homogenization and optimal design and uh Charles was a PhD student of Kagua and uh has been contributing to bringing level sets into the the newest generation with remshing which will be the the topic of the third lecture. Um so the lecture will be around 4550 minutes long. So 50 minutes long and then we will have questions afterwards.

[00:06:22]Please use the Q&A uh button. So at the bottom of your screen you should see a Q&A like a question mark. You can press on that and enter your your question during the the presentation and then we will take it at the end. And um if we have a lot of questions we may go a little bit over time and uh but I will try to keep the time. Okay. And uh with that uh Gregoire please if I stop sharing then you should be able to take over.

[00:07:05]Can you see my screen on the slide full screen?

[00:07:10]>> Yes, I see it now.

[00:07:11]>> Yes.

[00:07:12]>> Good. Excellent. So thank you very much Joey for the introduction and for the invitation. It's a great pleasure to uh present this uh series of three lectures on the level set method for structural topology and it's a pleasure to do it with uh my former student colleague and friend Shiao. All right so um so I I don't think I need to say a lot about shape on topology optimization and its importance in the design of a physical system. uh you can see here a picture which has been done by some of our colleague at NCIS and u u level set method is one of the method of course it's a bit later it appears a bit later after the density function method like omization or simp but it has some advantages and to introduce a topic let's say that it it is based on two ingredients uh a first ingredient which is the Adam method for geometric shape variation that I will very briefly discuss today but which will be mostly the topic of the next webinar and then the second ingredient is the oar city algorithm for advicting level level set functions which was existing uh some years before it was introduced in the context of shape and topology optimization.

[00:08:38]Uh okay maybe uh just a teaser you know this movie on a classical cant program uh problem sorry where you can see the evolution of this kind of sweet Swiss cheese initialization uh that moves so topology change as you can see you have also some slight oscillation that I will explain maybe later or in the next lecture good so the content of the webinar today is to introduce uce the level set algorithm for interface capturing then explain how it can be implemented in the context of shape and topology optimization and then show some numerical experiment. The next second webinar will be about how you compute shape derivative by the so-called Adamar method and then all the uh advanced tricks uh when you implement a level set algorithm uh the question of extending the velocity the choice of the optimization algorithm the fact that you can use more general mesh simplicial meshes all those very important details will be discussed in the next webinar and then the final webinar the third one will be about the last variant of the level set method which is based on remshing and then uh meaning that you always do uh computation on a body fitted domain. So the mesh is always perfectly adapted to the level set and so to the geometry of your of your shape.

[00:10:11]Good. Um so uh level set method were introduced uh around the uh end of the 80s so probably in 88 something like this by both Osher and Citian for interface capturing and you can see on the top some pictures of bubbles uh that you know after three instant change shapes and so it was really used first for combusion model and fluid mechanics model but of course it has been introduced a little bit later early 2000 in shape on topology optimization. Uh and today we are going to discuss some first presentation and first numerical experiment about this.

[00:10:55]Uh let me start with a basic present principle and let me first present the setting of shape on topology optimization. So a shape for me is a domain omega in R2 or R3. And a shape on topology optimization problem is a minimization of some objective function, some cost function J of omega among various shapes omega with constraint.

[00:11:19]And for simplicity, let me only consider equality constraint. Of course, inequality constraint are completely similar. Also sometimes a bit harder to treat from a numerical point of view.

[00:11:30]But um uh today I will not discuss that much constraint that will be the topic of the second webinar. Uh and I will discuss mostly unconstrained problem where I minimize uh an objective function. Uh there will always be a volume constraint but further constraint will be discussed uh next week. Of course, the difficulty of the business is that usually the objective function as well as a constraint do not depend directly on the shape, but they depend on the shape through a state. For example, the solution of the elasticity equation in your domain omega. And so that's a very indirect implicit nonlinear nonconvex dependence that makes uh things a bit hard.

[00:12:18]Okay, so here is the first example.

[00:12:20]That's the first example of thermal conduction. And so instead of being a shape optimization problem, it's really an interface optimization problem. You can see here uh square domain which is filled by two different material, two phases. Um so typically on on part of the external boundary you have a durish layer boundary condition. So the temperature is fixed to be zero. the remaining part of the boundary is isolated. So it's a n man boundary condition, homogeneous nan boundary condition. Uh and then you have a let's say a source term f uh inside your domain and then you want to find the best arrangement of the two-phase alpha and beta. So for example, you minimize the average temperature. So if you minimize the average temperature in some sense if the the phase beta is a good conductor compared to the phase alpha which is a bad conductor. So the phase beta is yellow while the bad conductor alpha is red. You want to have a kind of a geometry of the good conductor that will extract it which is produced everywhere uh up to the durish layer boundary where it can be evacuated because you have the fixed temperature equal to equal to zero. Okay. So in this type of problem it's not the external shape of the square that you want to optimize but it's the interface between the two the two phases.

[00:13:48]Now a more uh classical example in solid mechanics is where you indeed have only one phase and you know outside this phase this is void and so you are really doing shape optimization. So for example uh you have some part of your structure here gamma d which is clumped and on some other part gamma n you have some applied loads and then all the other part are free boundary. So nman boundary condition and what you want to do if you solve the linearized elasticity equation in your domain omega the solution u omega is a displacement e of u which is a symmetriized gradient is a strain tensor for simplicity I assume no bulk forces and I simp simply have a surface load on the boundary gamma n and a classical objective function is a compliance which is either the stored elastic energy or the work done by the load. Uh and you want to minimize this with a volume constraint and uh on on on on top you have a sketch of a classical 2D cany level but uh below you see a more complicated pylon that we will discuss a bit later in this in this talk.

[00:15:02]Good. But there are many other possible application. You could also do shape on top optimization in in fluid mechanics.

[00:15:10]uh on on the left you have an example of heat exchanger. So the white parts are a kind of an obstacle where you can extract it which is carried by fluid outside those white obstacles. on the left uh on the right sorry you have an obstacle inside a fluid flow and of course nowadays other application like let's say for electrical machine on on on the bottom left or for nanopotonic converter on the bottom right so there are many other possible application of structural optimization uh using level sets okay so one of the ingredients of the level set method is a method of differentiation with respect to the shape and let me explain uh the basis of this. Of course you know that when you uh optimize a function you need to compute derivative. So here you have a sketch of how you compute the derivative frime of x of a function f defined in r2 with value in r and so you do a tailor expansion. So by by you know computing f of x plus h and you do a tailor expansion like this. So that's very classical. But of course when the optimization variable is a shape, what is a small variation of a shape? You know, you cannot do addition of shapes.

[00:16:32]There is nothing like a vector structure on shapes. So what is a convenient way of defining increment and computing a derivative for shape functional? Uh this was provided by Adamar at the beginning of the 20th century. And the idea of Adama was to look at a variation of a shape with respect to a reference shape. So on this picture here you see a reference shape with this black boundary omega. And then you look at very specific variation of this shape omega which are called omega index theta and which are defined as identity plus theta applied to omega.

[00:17:14]What does that mean? That mean that any point x in omega is displaced to a new position which is x plus theta of x. So sitta is in some sense an increment a displacement for the shape omega. And so you obtain the new shape with a boundary which is a red dotted line. And every point x which was on the boundary of the initial domain omega moves to a new position x plus theta of x that define the boundary of the new shape omega.

[00:17:46]Okay. So in some sense we are looking at a very specific class of variation which is kind of a restricted class of variation and those variation are parameterized by this vector field which we assume to be smooth let's say C1 uh in the SQL from RD into RD and then um oops sorry excuse me uh I'm going too fast and then um if SA is small u a very classical result is that omega and omega sitta have the same topology. It's well known that's one definition of topology that if you make a continuous deformation of a given domain then this continuous transformation has the same topology. So topology in 2D means the same number of holes in 3D it's a bit more complicated. It means the same number of holes and uh okay so I will not go too much into details here but there is this very uh well-known book by and Pierre about ship variation and optimization that you can look at if you want more details on those aspects good and then the definition of what is a derivative is it's going to be a derivative with respect to this variation uh to this parameterization of the variation so derivative with respect to this vector field theta So the definition of a derivative is you make a tailor expansion for your objective function J evaluated at your pur domain omega index theta. It will start like J of omega and then you have the derivative which is a linear form with respect to this variable parameterization of your deformation plus a small remainder term that mean that must be small with respect to the norm in the space of c1 function. Okay.

[00:19:37]Uh and then the uh main interest of computing this derivation is that of course if you want to decrease because you minimize your objective function you are going to choose a vector field so that the directional derivative is negative and then if you take a small increment small decent step toe which is positive uh if toe is sufficiently small uh this will allow you to neglect the rem term in the tar expansion and then in the direction of Sitha for a small decent step toe you are sure to decrease the objective function and that's the beginning of a gradient algorithm okay let me give you u in this first lecture only two example of shape derivative very simple one uh the first one is when you integrate a fixed smooth function f ofx on the domain omega so the domain omega appears only in the fact that it is the domain of integration of your objective function.

[00:20:39]Then a classical result tell you that the directional derivative of your function J in the direction of SA is simply the integral on the boundary of the domain omega of the same function F multiply by the normal component of the vector field SA. So to understand uh this you just do this small sketch. So you see the domain omega with a boundary in black and then the variation omega sitta as a boundary which is a dotted line in red and you can see that if sitta is small then the difference between the two domain is small and the difference is concentrated along the boundary of omega and the sickness of this difference is exactly sa n. Okay.

[00:21:24]So uh if you make the difference between the two integrals of f on omega and omega then it's something which is carried by the boundary of omega with sickness sa n and then of course you still have the integon f. So that's a kind of intuitive proof of this formula here when you make a tailor expansion of your objective function.

[00:21:50]uh of course you can make a rigorous proof of this and then this rigorous proof is really u just an simple consequence of the so-called Reynolds or transport theorem which is well known for people who do fluid mechanics which is that if you compute the time derivative of an integral of a function that depends on time on a domain that depends on time then uh by the chain rule lema it's a combination of first the partial derivative of the function f that depends on time on the domain omega t. But then you also have a contribution because your domain is moving and this contribution is exactly uh of the same type of the objective function that we of the derivative of the objective function that we have just considered.

[00:22:35]In particular, if you apply this result to the function f which is one, then your objective function is just the volume. And then you recover that the derivative of the volume is simply the integral on the boundary of the normal component of your vector field theta.

[00:22:49]And by integration by part it's the integral on the domain omega of the divergence of this vector field. And in particular if you have chosen a direction of derivation which is divergence free it's well known that it does not change the volume. So the derivative is going to be zero.

[00:23:08]Okay. A second example which is very uh practical is when the objective function is the integral again of a smooth fixed function J but now it's an integral only on the boundary of of the domain omega D omega and then the shape derivative in the direction SA is again an integral on the boundary of omega of the normal derivative of your function G plus the mean curvature kapa is the mean curvature of the domain times the function G again multiplied by the normal component of your vector field theta. In particular, if you apply this result to the function G, which is a constant equal to one, then the objective function is a perimeter. Okay.

[00:23:53]And then it's well known that the derivative of the perimeter is proportional to the mean curvature.

[00:24:00]Uh actually I have again u small sketch here that shows you a domain omega with a dotted black line and then omega sa is the domain with the red boundary. uh and then uh it's well known that if you want to minimize a perimeter you should choose a descent direction SA which is proportional to minus the mean curvature times the normal and this specific vector field has a tendency to smear the bumps of omega. So you see here you smear the bump on on on on the red parts here but on the blue part you will seal the hole. So that means that you will try to have the more possible flat or convex shape that will in some sense minimize the perimeter.

[00:24:59]Good. uh more generally there is a theorem which is called the adamar structure theorem that tell you that for any objective function j of omega when you compute the derivative in a direction theta it's always written as an integral on the boundary of some integrant v depending on omega times the normal component of your vector field theta okay so this of course this integrant v index omega that depends on the specific specific choice of your objective function and in particular that will depend on the state uh u omega which define your objective function and of course on some adjint state but we will review this theory uh in the next webinar.

[00:25:46]Uh okay this is an integral that may hold in the sense of distribution. So what does that mean? That means that it may also depends on tangential derivative of your vector field theta on the normal component of your vector field theta.

[00:26:02]Okay. One one uh key consequence of this Adamar structure theorem is that this derivative in the direction SA depends only on the normal component of SA on the boundary. So in particular if your vector field uh has a normal component that vanish meaning that what you do is you simply u you know convect the boundary on it itself. So it's only kind of a rotation then it does not change your objective function because you do not change the geometry of your of your domain. Okay. So a tangential vector field it's simply kind of a convection or you just slide the boundary on itself. So that does not change the the the domain. So of course the the shape derivative should vanish. Nevertheless, uh if you choose a different vector field, in particular, if you choose a vector field which is proportional to the normal with a scalar with a scalar coefficient in front of it, which is minus V of omega, you ensure that the directional derivative is strictly negative. And so uh by this uh change of the geometry induced by SATA you will decrease your objective function.

[00:27:21]Good. Now let's discuss how we uh uh implement those principle. And so the idea is to do of course a gradient algorithm that we will initialize with some initial guess omega 0 and then for every iteration n we will compute the solution of your state equation u index omega n. Then possibly an adjint p index omega n. And then from those two state u and adjint p we can compute the shape derivative. And from this shape derivative we can induce a descent direction sa n. And then we are going to move the uh previous shape omega n to a new shape omega n + one by simply moving all point x in the domain to a new position which is x plus descent step to n times theta n of x where sa n is this descent direction inferred from the shape derivative. And then you iterate like this until you uh uh converge to an optimizer shape.

[00:28:27]Okay. So uh this algorithm so far is kind of theoretical and in particular we have to decide how we update you know omega n to omega n plus one. How do we apply this vector field to your to your shape.

[00:28:44]So what was classical was a lagron approach which is more or less the approach which is used in geometric shape optimization. Very classical geometric shape optimization. So what does that mean? Each shape omega n has a triangular mesh tn. Okay. And this triangular mesh of course is used to compute by the finite element method the state and the adjint. Right? And then the update from omega n to omega n + one is simply realized by pushing the nodes of the mesh tn by this vector field with a distance step to n and times theta n.

[00:29:24]So what does that mean? You take the position x here of every node and in the direction of the arrow which is the direction of your vector field theta n you move the points. Of course, the difficulty is that if your vector field or your distance step is too big when you move the points of your mesh, maybe the new resulting mesh is not valid in the sense that you have inverted triangles. Okay, you may have very flat triangles. So the difficulty when you do this is that you need to check uh that uh the new position of the nodes of your new mesh uh are still valid and maybe you need to remesh because uh the the the the elements of your mesh are not of a good quality uh inverted too flat or whatever. So that's the classical difficulty when you do a geometric optimization that when you move the notes maybe you need to adapt or to remesh at some point. Nevertheless, it works and let me show you an example for a so-called cany level problem. So you see the geometry the dash part is where the structure is fixed clamped on gamma n on the right part you have this vertical load and I start an initial configuration with seven holes and then I minimize the compliance with a constraint on the volume.

[00:30:49]So you should see uh the movie where you see the mesh which is deformed of course each iteration is computing a derivative on moving the mesh. Uh so it seems to work pretty well but but you see it stopped at some point and actually it's not even a local minima.

[00:31:10]It's not even a local minima because the the the leftmost bar vertical bar here becomes thinner and thinner and at some point uh it wants to become so thin that the meshing algorithm cannot mesh the the new mesh or you invert triangle. So actually it does not converge um because it has reached a local minima. It stopped at some point because the meshing algorithm tell you well you I cannot mesh I have two small elements in this simbar. So the algorithm really wants to break the bar. Okay. So that's the limit of classical geometric optimization and the interest of the level set algorithm is that we are going to uh go beyond this limitation of geometric shape optimization.

[00:32:00]So a few words about the level set method. So as I told you before this was devised by Osher and Sitan in 88 uh to represent to to help them to solve a combustion problem. And here is an example. Maybe I should uh oh well I think it will cycle. So you will see it again. This is an example in fluid mechanics. That's a collapse of a water column. Blue is water. Red is air. And you see how water under the action of gravity is falling down and you have this interface between water and air which is described by a level set function. So that's a classical application of the level set algorithm. Uh there is a picture below on image segmentation by the method of so-called active contour which is again using level set. So this were previous application of the level set but of course today we are going to discuss application of the level set method for shape and topology optimization of mechanical structures.

[00:33:04]Good. So first uh let me come back to the lagon point of view when we have a a family of shape omega of t which are evolving with time. So there is no time in optimization but in some sense you could think about time being absodo time corresponding to the decent step. So we have a family of uh uh shapes which are better and better for optimization problem indexed by the decent steps or the sodo sim t. So those shapes evolve have a boundary gamma t and they evolve with some given velocity field v of t and x and this v correspond to the previous vector field theta which was a descent direction. Okay. Now let me define the concept of characteristic curves. So for a point initial point x0 k depending on the initial position of x0 and t is a solution of the od uh where the right hand side of the od is this velocity field v evaluated at the same point kai. So that's exactly what we call characteristic. So if you see uh on the picture below I have a plot of my vector field v that's the blue arrows and starting from some point x0 I have a line which is uh the solution kai of this od. So that mean that at every point on this line the tangent vector field is exactly my vector field v. Of course for some different init initial point x0 x1 and x2 you have different uh characteristic curves.

[00:34:38]And so an intuitive definition of an evolving domain is that omega t is just precisely the set of points k which are on the characteristic curves that started at some point x0 which was in my initial uh domain omega 0. So that's a a precise definition a lonjian definition of an evolving domain omega t. So you have a sketch of uh of this definition below. So if you start from some uh domain omega t0 you take a point x0 let's say on the boundary of this domain. So you evolve along your velocity field v and you arrive at some point t on the uh characteristic point k starting from x0 at point t. Uh and that defines a new point on your domain omega t. Okay. So that's a definition of evolving domain based on solving OD's uh computing characteristic curves.

[00:35:40]So this langen viewpoint is very nice but it has a big drawback which is exactly uh exemplified on this picture.

[00:35:49]Imagine that the velocity field is such that the boundary of your domain you know the velocity field here is in red at some point they will cross. So now uh the the boundary of your new domain has has crossed and then there is a very tough question because this intersection uh zone here is both inside view from the right but it is uh also from the uh uh outside. So you have to decide you know am I inside or am I outside? So where is my where is my domain? In particular, if you have to compute your velocity field by knowing where is the inside or the outside or by knowing what is the normal vector of the curvature vector at the intersection point you have two normal vectors. So which one is the right one.

[00:36:39]Okay. So you have a difficulty which is the classical difficulty for lonjan method and then the main idea of the level set method is we are going to completely forget about this lronen setting. In this langen setting by the way uh displacement is reversible because you could solve the OD for the characteristic of for positive time or for negative time. So if you just revert the sign of the velocity you can go uh back and forward uh it's completely reversible. The level set method is going to forget about the lagon setting forget about reversibility and it will use now aan setting which will be irreversible and that will give you a definition of what happened after crossing uh boundaries.

[00:37:31]So uh the idea as probably most of you knows is to introduce a function fi which is called a level set function which is going to be by definition zero on the boundary negative inside and positive outside. Okay. So for example if you have a shape here of this cartoon character uh you define a function here which is represented in 3D. This picture here is in 2D. The function file is now in in a higher dimension. It's negative.

[00:38:01]So blue inside the yellow region of your domain omega and it becomes increasingly green and red outside. So positive.

[00:38:11]Okay. So that's one way of representing a shape. Of course there are many such function and usually we prefer to take the one which is a sign distance but any other function will will do the job. Uh the first interest of this is that it defines implicitly the geometry and in particular it is the computation of a lot of geometrical quantities like the normal vector. So if you look at the normal vector one definition if you have a level set function is that it's also equal to the gradient of your level set function divided by its modulus at least where the gradient is not zero and then that give you also an extension of the normal vector not only on the zero level set but on on any value of the level set function ph of x being equal to a constant. So you have an extension of the normal vector everywhere. Uh similarly you can compute the mean curvature which by definition is a divergence of the normal vector. So you can compute the mean curvature everywhere again as an extension. And of course you can use this level set function uh by means of ani function to compute integrals of another function defined on the domain omega or on this boundary gamma.

[00:39:29]Good. Now what is the equation for the evolution of such a level set function?

[00:39:33]So let me uh come back to my lrangian definition of evolving domain. So starting from a domain omega 0 if it evolves according to some velocity field I have a lronian definition but now I want to translate this langian definition into an ularan definition of a moving domain.

[00:39:53]So remember that again u u only the normal component of the vector field v matters if you have a tangential component right uh like it's plotted here on the left picture it's only uh a kind of sliding the do the domain on itself so the green part will simply move on a new position the red part will move to a new position but you do not change the general form of your of your domain so we can really ignore the tangential component of V that does not matter in the transformation of your domain.

[00:40:29]Okay. So uh when you have a domain evolving according to a velocity field V uh there are three different cases. The first case is when the velocity field V does not depend at all on omega T. So V is defined on transport world. So that's what is can be called a passive transport. That's the easy case. Uh a second more complicated case is when v depends on local features of omega. So what do I mean by local features?

[00:40:56]Typically it depends on the normal vector n or on the mean curvature. So that's quantities that are computed locally at every point x on the boundary. And then the third cases which is a case of interest for us is when V depends globally on the domain omega because for example it depends on the solution of a model of a PTE both on the domain omega. So that means that you need all information on the domain to compute the velocity. It's not just given local formula.

[00:41:31]Okay. Let me give you some example of the easy case which was actually the original example of Osher on CT which is the representation of omega t being a burnt region. So that's the region where the gas has already burned and the combustion is represented by the fact that the velocity which evolve the domain omega the burnt region is proportional to the normal vector with a constant normal speed. So C is a constant positive and then you see the evolution of the region of the burned gas. It's simply an evolution along the normal with a constant velocity. So you start with a small bubble and then the bubble increase from left to right.

[00:42:16]Uh something which is a little bit uh more complicated because here in this case the velocity depends on the property of the domain is the evolution by the mean curvature. So that's also a classical problem. So the velocity field is minus the mean curvature times the normal. And you see the evolution from the left to the right. And you see that this has uh the the the the influence of this type of murvature flow is that it will uh fill the bump feel the crease and uh flatten the bumps. So that makes this non-convex initial shape become more and more convex uh like it is. And actually it is a theorem proved by Grayson that in the end any shape even if they are nonconvex at the initialization will end up in something which is convex.

[00:43:10]Okay. In a more general setting um let's look at what will be uh the u the evolution equation for for the level set five that represent an evolving domain omega. Okay. So if by definition phi is always zero on the boundary, okay, you can replace a point x in a boundary by the characteristic uh the solution of the characteristic equation k starting from a point x0 which was initially on the boundary of the initial domain. Okay. So if x0 was initially on the boundary of the domain, k of x0 at time t is on the boundary of your moving domain. Okay. So by definition of the level set function, this should be zero. So now I differentiate this definition uh of the level set function being zero on the boundary. I differentiate it with respect to time. Time appears at two uh position in this equation. it appears first as the dependency on the level set function fi but it also appears in the dependency of the characteristic curve.

[00:44:23]So that's why I have two when I differentiate respect to time I have two terms first the partial derivative of ph with respect to t and second by the general rule lema the gradient of phi multiply by the time derivative of the characteristic point and this time derivative this is by the od equal to my velocity capital v right so I get the first level set equation which is just an equation of adection for phi dy over dt plus v which is a vector field nabli equal to zero.

[00:45:00]So that's the level set advction equation. But remember remember and that's very important uh only the normal component of the vector field matter.

[00:45:11]Everything which is tangential does not change uh the the the the shape. It it change maybe the parameterization of the shape but it does not change the shape.

[00:45:20]So actually only the normal component is important. So let me uh replace capital V by only its normal component small V.

[00:45:29]And then the normal field in the level setting can be represented by Nablafi divided by its modulus. So if I replace capital V by this representation where only the normal velocity matters uh I end up with a new version of the level set equation which is an Hamilton Jacobi equation which is a nonlinear equation because you see that you have the modulus of the gradient so that the norm of the gradient so that makes the equation nonlinear but the interest of this equation is that it is guaranteed that you always evolve in the direction of the normal vector at every time. So that mean that you are always very efficient. It's a in some sense a more efficient way of writing the evolution that the previous one because in the previous one maybe you do some unnecessary adction in the tangential direction and that is not good because that may from a numerical point add numerical diffusion.

[00:46:29]Good. um I did this computation for the zero level set but the same computation work if I was looking at the C uh uh level set. So if I was looking at five equals C the same computation all true.

[00:46:44]So that mean that uh my equation for the level set at vection or the Hamilton Jacco they all true everywhere because everywhere you will sit on the level set of five equal to some constant. Okay, good.

[00:47:01]Uh [clears throat] now um of course uh the evolution depends uh uh on this definition of this velocity field and that's something which is very clear if your domain omega and your velocity field are smooth. But what happen if uh singularity appears?

[00:47:19]And it turns out that even if you start from very smooth initial shapes and you have a very smooth vector field uh at some point singularity will appear. Let me explain that uh in a in a very specific setting.

[00:47:33]That's my initial domain here in in yellow. And I'm looking at this part of the boundary which is a kind of a curved part. And I'm looking at a very simple velocity field V which is simply equal to the normal vector. Okay. So that means that I'm at vecting with a constant velocity one in the direction of the normal and I plot here different position of the boundary of the domain uh after some increasing time. So you move in the direction of the normal. So the normal is towards the exterior. So you see that you move and at some point you arrive at this red curve where the blue point is a singularity. It's a corner. Let me explain why you uh end up in a corner. Because if you end up after the corner the corner here is here this red curve then you have exactly this process of crossing domain. Okay. So uh if you if you do a lrangian evolution of all those point with velocity one uh you end up after this uh first point of singularity uh with boundaries which are cross because if you move the points independently well you know you you you end up with a multiffolded boundary. So what is the meaning of this? And the idea which was proposed by Osher and Cition was to impose a kind of entropic critarium saying what is burnt. So this yellow part is burnt will stay burned.

[00:49:01]So that means that all this triangular part is completely irrelevant that's inside the domain. So you should eliminate this crossing and so what should be retained by this entropy criterion is exactly uh only the part which is not intersected. So eliminate this triangular part here and that's what is the uh correct boundary in the sense of the level set equation.

[00:49:27]>> Sorry to cut in quickly. I just wanted to let you know we have about eight minutes left.

[00:49:33]>> Okay good. So I should I should go a little bit quickly. So uh so there is a a notion of viscosity solution that explains this in mathematical details uh and I'm not going into more details but that's exactly a big picture here showing the same uh the same idea that's the way it evolves and that's the way uh it allows you to change topology also.

[00:49:59]So for example uh you see on this cany lever example the evolution of the level set equation uh and the change of topology that correspond exactly to this notion of viscosity solution.

[00:50:12]Good. Uh so I pass um u the the floor to SH.

[00:50:20]>> Yes. So maybe I can share my slides on my on my own.

[00:50:26]Okay. And can you see my screen?

[00:50:30]>> Yes.

[00:50:31]>> Okay.

[00:50:33]All right. So, thanks for passing me the hands. Uh, all right. So, let's now deal about how we are going like to actually solve level set equations on the on the computer. So, as Gagar said, we have a shape the setting is this. We have a shape omega t that describes an evolving shape and it is described a level set function meaning that omega t at any instant t is the set of points where phi is negative and it evolves through a velocity vtx a normal velocity which in our case of shape optimization describes a shape gradient. So the Hamilton Jacobe equation reads like this exactly like what Gregor presented and the velocity field is actually very complicated at any point t the shape gradient is going to be some kind of taste of the shape derivative of the objective function under minimization. So the only thing that we can do possibly is to freeze to to to cut the time interval into a series of sub intervals and on each such sub interval to freeze the velocity fields. Meaning for all time t in the period tn tn plus one vtx will be frozen to its actual value at tn and then we are back to solving a steady state equation with a veloc with a normal velocity v of x like this. So how we going to solve this equation which is not at all trivial because it's a nonlinear equation. It makes appear here the the the absolute value the the norm of the of the gradient of the function of the solution function phi. So we are looking for first order amid jacobi equations and there is a well established theory to solve such equations uh which appeals to the theory of viscous solutions that ka mentioned.

[00:52:13]So it proposes in particular adequate numerical schemes to catch the good solution the good discosity solution uh to this equation the one that has physical sense in the sense that it imposes this kind of entropic criterion in the cases where it is relevant and the key point in the device of such schemes is the word upwinding so that's what I'm going to explain in the simple case where we are in 2D and the the computational support is a cationian great so just to make to try to to give a hint of Why upwinding is needed? Let's consider the very particular very simple case where the normal component of the velocity is constant. It's C. It's a positive constant. Then it's some after some math rolling up the math. You can see that the viscosity solution to this guy is this function here. It's a distance shifted by the coefficient C of T. Meaning that when C is positive, then the lever function zero isol function.

[00:53:07]We do that with rounding corners that go away and away from the from the square.

[00:53:12]And when the velocity normality is negative then the shape will shrink into itself like this and the corners here they will stay corners actually. So what you see is that there is some kind of irreversible uh uh pattern that is happening here because if you shrink and then go back in the in the other sense then you will see that the corners that appear here will get rounded meaning you don't obtain the same thing by applying the forward solution of this jacob equation then backwards or the other way around. there's some kind of irreversibility here that only upweening uh can catch and that there has to be some kind of translation of that in the numerical scheme. So there is some something that will appeal to people that speak to people that know a little bit about hyperbolic equations. If you consider this the equation that I mentioned in 1D and formally that you take the derivative with respect to the special variable in this guy you obtain this one-dimensional hyperbolic equation then it's well known in some sense in this in this word that the theory of shock waves predict that if you consider very simple data uh here with a value ul on the left part and the value ul on the x on the right part then depending on the values of ul and you are you may have a shock wave or a reflection cone.

[00:54:31]The shock wave would correspond to a sharp corner of psy and the reflection cone will correspond to a rounded corner of. So this is really something that tells you that you need to do upwind schemes to find this very particular behavior of u of the solution. So how do you do in practice? So let's consider the very uh the very concrete situation where we have a custo grid of space. We discretize the spa the time by a series of time steps like this. And we discretize the space by this cartian grides here with a step size delta x and delta y in the x and y directions. So we introduce a few difference quantities dj x and dj minus which is here. And this quantity here is just taking the value here minus the value here and to divide by delta x.

[00:55:21]And when I consider the minus x then I'm basically doing it backwards like this vi j minus vi minus1 j etc etc. So the level set of engine equation can be solved by this particular scheme here.

[00:55:36]So this is an explicit first order scheme where the update is this discretization this seemingly weird discretization here of the gradient term which is over here which is over here sorry and this is exactly what we need because it features discretization nabla I plus and I minus of five which are upwind and downwind discretizations of the gradient nor which is here. Why are there upwind discretizations?

[00:56:05]because they only use values among the neighbor values of fij that are larger or smaller depending on whether it is downwind or upwind respectively. So meaning that the discretization of the hamiltonian which is here. If you just look at this formula in di then you can see that it is upwind in the sense that if vig is positive then this value is used for the discretization of the gradient norm which only uses value of vij at a point i j that are less than fj that are smaller than fj and if vig is negative meaning the information goes in the other sense it only uses values of among the neighbors of j that are larger than vi and this scheme can be proved to be converted to the unique viscosity solution to the am jacobi solutions. So this seems a little bit weird and a little bit difficult but really the formula is like this and if you program it this like this is like 10 lines in Python and it works exactly what it does exactly what you need to solve the hon Jacobi equation. So you also have high order variance of this strategy and all these things they are you have very nice numerical packages that allow to solve this equation.

[00:57:17]So another point that is very important in the numerical practice of the level set methods is the need to initialize or redistiate send distance level set functions. So as they were mentioned there is not one level set function associated to a particular shape. There are a lot of them but there is certainly one which is more beautiful than all the other ones. This is a sign distance function to the shape which is minus the distance to the boundary when you are inside zero when you're on the boundary and plus the distance of the boundary when you're outside. So this is a 1D visualization of any an arbitrary level set function of the interval 01 and this is a visualization of the sign distance function to this interval and this is the solution is often preferred because it has unique gradient no almost everywhere.

[00:58:03]>> So actually >> sorry to cut in just saying the 50 minutes up but we have five more minutes of the sort of allowed hour. So I suggest we we spend that five minutes on you finishing up.

[00:58:15]>> Okay.

[00:58:15]>> We just have a few questions. So it's okay we use that time.

[00:58:19]>> I'm I'm finishing one minute and then leave the floor for to to conclude for the session. Uh okay. Uh and so here it's very important in numerical practice to have levels and function that are close to send function. This is not something that the theory dictates.

[00:58:35]This is something which is dictated by the need to for for accuracy. This is it creates numerical artifacts when you don't have this feature. So in practice it's important to initialize distance functions as sign distance functions and it's important also to periodically reinitialize the current fun the current level function as a send distance function. So you have adequate schemes based on the force merging method for instance to do that and I'm not entering into any more details uh for lack of time but these are classical issues in the numerical practice of the send distance of the level set method. So here are some final clues about the level set method. It's usually very cheap. It's easy to implement at least on cartian grids. You have very nice uh open source libraries on the web. uh if you don't want to implement them they are very very cheap uh they are say some kind of linear or n login like algorithms and in the next webinar we'll see how to implement these type of techniques on computational support that are not casian grids but are triangular grids for instance and so gagwa I think the floor is yours I think I can stop sharing Okay.

[01:00:09]Okay. Sorry. Very quickly, let's go to some example before um Yeah. Uh before we we [clears throat] answer some questions. So uh here is an example for compliance minimization under volume constraint. So it's a so-cal electric pylon where you see some loads which are applied to the blue points. Um oops.

[01:00:37]Yeah. And here is the movie.

[01:00:40]uh so the initialization was a full domain and uh in 3D because I'm sure I will have questions in 3D change of topology are very easy because uh most of the topology that we are looking at areles and not internal holes while of course in 2D uh uh change of topology is a bit more complicated there are some obstruction and so that's why we very often initialize with a so-called Swiss cheese initialization with many but again uh is another example of a so-called dome which is optimized under hydrocatic pressure with five anchor points which are red at the bottom and again uh this one you see that the initialization is completely trivial but it's able to change topology very easily. So uh change of topology is much easier in 3D than than in 2D where I must admit it's a bit more uh tricky.

[01:01:40]All right. Here is an example of a compliant mechanism. So uh we act on the blue part and we want the jaws which are red to close and uh you you find this optimal gripping mechanism which is in the undeformed configuration on the left and deformed. You see that jaws indeed are closed on the right. Again a very complicated topology has been obtained with the level set method. It works also for other objective functions. So here is an example with uh LP norm of the stress which is concentrated in this corner here. Uh so the optimal shape for compliance is this one. And if you do an LP norm for P equal to two five and 10 you see that for 10 indeed you get a rounding of the corner as should be expected which makes a big difference with a compliance solution.

[01:02:35]Okay to conclude there are some assets of the method uh and some drawbacks. Uh clearly one drawback is that the shape is not mesh exactly. That's the same drawback as for density method and in particular the choice of theat material could be tricky for the design of compliant method. um it's quite sensitive to the initial design but mostly in 2D uh and then um uh admittedly it may be a bit complicated but actually there are many open source libraries as Sh said and so that should made the implementation of the level set method easy for newcomers. Thank you for your attention.

[01:03:15]>> Yeah, thank you for your attention.

[01:03:17]>> Thank you so much. Thank you and sorry for pressuring you on life.

[01:03:23]>> We understand.

[01:03:25]>> So, uh we have some questions. So, I will just ask them uh for the people who have submitted them. So, um one question is why does the shape derivative only depend on the normal component of the vector field on the boundary in practical simulations? How do you compute the normal vector accurately on complex measures?

[01:03:48]>> Okay, so why it depends only on the normal uh component is exactly explained uh by this kind of a picture because uh if you move tangentially it's only a kind of a sliding of the boundary on itself. So you don't see any change globally. Okay. So uh if you if you have a a vector field that move tang tangentially you don't change anything.

[01:04:11]So the shape stays the same. So that's why it has a zero derivative. So only the normal component matters. Now how do you compute the normal u vector? So u uh the starting point is to say that um this was maybe before sorry yes that's here uh if you have a level set function nabli divided by it modulus is an extension of the normal vector. So that's one way uh kind of naive way of computing the normal vector. We'll discuss maybe some uh details in the next webinar for more general measures.

[01:04:48]>> You you have actually a lot of methods to do that. You could do that by exactly this formula that use guar. You could use this formula with a list square taste. You could use some kind of larger stencil reconstruction. Some people also when you work on cagen grids use essentially nonatory different differences. So I mean you have like a lot a lot a lot a lot of numerical methods to do that but when you are in a level set context it's probably better to use to do something based on the formula that pointed out.

[01:05:19]>> Okay thank you. Um so we have a question that came quite early on so potentially it's been answered on the the way but um so what is the simplest way to connect the shape derivative of the objective function? Oh wait sorry what is the simplest way to connect the derivative of the objective function to the shape derivative? Is it as simple as using the Hamilton Jacobi equation to connect the normal velocity to the time derivative of the level set function?

[01:05:52]>> We'll we'll discuss this uh in details in the next webinar. Uh so that's what we call the inbur extension.

[01:06:00]>> Okay.

[01:06:00]>> So yeah, sorry not to say more but that will be typically the topic of the next webinar. One of the topic of the next webinar.

[01:06:08]>> Great. So then FAD you will need to come back next time. Yes, I'm sorry about that. [laughter] >> And uh we have another question. So in the iconal equation, can C be zero? What would happen to the shape then?

[01:06:23]>> Well, nothing happens if the velocity zero the shape stays the same.

[01:06:29]>> Yeah. Okay. Yeah.

[01:06:33]And um do we need re reinitialization even if we do small displacements and then extract our contour of interest?

[01:06:44]>> Actually redistantiation it's more it's really something like a bit empirical. I mean you need to redistiate not necessarily at every step of the practice of method. You need to resist redistiate often enough to have a level set function which is neither too steep nor too flat but possibly not too often to avoid that it be large too large computational burden or to avoid puring too much level set function. Uh so it's a bit usually you do resistation every five four five iterations but it's really something like uh depending on the on the sense of the winds.

[01:07:26]>> Yeah.

[01:07:28]Um then there's one question I'm not quite sure I understand what they mean but it's is the level set solved on immersed cartesian grits? [snorts] M >> um for the moment um in this presentation here we have presented very briefly a discussion of the algorithm for the Milton Jacobi equation here which was on cartisian grid. So admittedly for the moment we discuss only the case of cartisian grid. So that means that you know for simplicity you could take the same mesh for your finite elements but we will discuss uh uh the same uh solution of algorithm but on simplicial mesh next time next week next next webinar sorry. So um >> okay >> for the moment we discuss only the case of cartisian grid but we will discuss the case of simplic measures next week next we know sorry >> well it is actually next week >> that's next week okay okay I've never I I don't remember [laughter] >> that's okay and then um yeah that's the questions we have for now um if you have any more questions come next time and ask them uh and Um, also I know uh you have posted the slides on uh Charles's website.

[01:08:58]>> So maybe you have the link here or otherwise it's >> just a minute. Yes, I have it here. It was before people got connected. So I I'll put the link in the chat. Yes. Can everyone does everyone?

[01:09:11]>> Yeah, that was only for us. I've put it for everyone here. So here's the link if wants to see.

[01:09:18]>> Um yeah, thanks a lot. So I will just >> steal the share again and then I will just show that the next two lectures in the series are next week. So that's next Wednesday, the 20th and then the week after on the Thursday. And um you can go to our website there you can find the zoom registration links. So just Google top webinar um and and you should find it.

[01:09:53]>> Okay. So very sorry for being too long.

[01:09:55]Uh we shall uh we shall try to be shorter next time. [laughter] >> That's that's perfectly fine. Thank you very much again for uh for your presentations and thank you uh thank you all for joining and uh see you next time. Thank you for this invitation and thank you for attending.

[01:10:15]>> Thank you.

[01:10:15]>> Thank you.

[01:10:16]>> Bye-bye. Bye-bye.