# Lecture 27: Finite Elements in 2D (part 2)

Flash and JavaScript are required for this feature.

Instructor: Prof. Gilbert Strang

The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high-quality educational resources for free. To make a donation, or to view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

PROFESSOR STRANG: So, let's see, you probably guessed on that quiz problem three, it wasn't what I meant. I get a zero for that problem. But you'll get probably good numbers, so that's an election gift if it comes out that way. So they're all in the hands of the TAs to be graded. We have a holiday Monday, I think. We come back to Fourier. Now, so we just have a concentrated shot at Fourier, just about eight or nine lectures in November. So stay with it and that'll be of course the subject of the third quiz. Which will have no mistakes. It'll be solved by the TAs in advance and we'll spot things. So, and if we have the quizzes to return to you by Wednesday that will be great. I hope so, but they have a big job. A little bit, looking far ahead, the end of Fourier, the quiz is December 4th, I think that's a Thursday. And that's the end of the course. So December 4th. So we'll be ending the course a little bit early. Because I'll be in Hong Kong, to tell the truth. And and we've done a lot, and with the review sessions we're really doing well. So, that's the future, Fourier. Today is an important day too, finite elements in 2-D, that's a major part of computational science and engineering. The finite element idea, the idea of using polynomials, you can find in some early papers by Courant, a mathematician in New York, and by a guy in China, neat guy named Feng Kang. But those papers were sort of, you could do it this way if you wanted. It was really the structural engineers in Berkeley and elsewhere who made it happen ten years later. And the whole idea has just blossomed. Continues to grow. So I had an early book, in the '70s, actually, about the mathematical underpinnings, the math basis for the finite element method. And many other finite element books have come. Professor Bathe you know, teaches a full of course on that.

But, I think we can get the idea of finite elements here. We did them in 1-D, and now there's a MATLAB problem and I'd like to just describe that particular problem if I can, as an example. And, of course, you would use the code that's printed in the book, and that's available on the website just to download. But the problem is not on a square domain. It starts on a circle, so that the first lines of the code, calling the MATLAB command squaregrid, are not applicable. So you have to create, then, a mesh. Well, I have a suggested mesh so I'll draw that, and then from that you want to make a list of all the node points. A list, P, of-- So what the code needs is two lists. Well, let me draw a picture of, well, it's a circle. And I'm going to be solving Poisson's equation. The equation will be -u_xx-u_yy=4, in the circle. So it's Poisson but with a constant right-hand side. That will mean that all the integrals of f times v, the right-hand side of our discrete equation will be, the integrals are all easy because we just have a constant there times the trial function. OK, and then on the boundary it's going to be u=0. On the boundary. So it's a classic problem. And we can say what the solution is. So it's one with a known solution. I think it would be x squared-- No, I guess one, one minus x squared minus y squared. This should all be on the .086 site. I just didn't have a chance to look this morning to be sure it got up. So you can watch, here. So that, I hope, does solve the problem. Two x derivatives give us a two, two y derivatives another two, so we get four.

So we know the answer; the question is, and I'm interested in this question, for research reasons too, is what's the error when you go to a polygon? You go to a-- These curved boundaries don't get correctly saved. You approximate them by straight lines. That would be the first idea. And with this, all this symmetry, let's keep the problem nice and use a regular polygon. So maybe I'll try to draw one with about eight sides, but-- OK. So we impose u=0 at these nodes. So u is zero at those nodes, and then we have a mesh. So we want to create a mesh. OK, so with all the symmetry here, the natural idea would be to start with eight pieces, or M pieces if I have, this is a regular M side, let's say, and I'll take M to be eight in this picture. And I think we can work on just one triangle. By rotational symmetry, all those triangles are going to be the same. So I think our domain is really this one triangle here. That's where we're working. And in that triangle, I think we have zero boundary conditions. And across this edge I think we have natural boundary conditions. Slope zero, if I see the picture correctly. The rotational symmetry would mean that things are not changing. That every triangle is the same. So I think on these boundaries it's the Neumann condition, dU/dn=0. And I'm frankly not sure what to do at the origin, so I'll maybe just try both ways and see.

And so this point would be (cos(pi/8), sin(pi/8)). We know where they are. But that's going to be what we have to list. We have to list where are the coordinates of all the mesh points. So let me describe the rest of the mesh points and then see what that list would look like. And once you've created the list, the code will take over. And then you plot the results and see what's going on. So let me suggest a mesh. It's pretty straightforward. I just divided this piece into N, this center thing into N-- So I called that distance h. So Nh gets me out to here. Whatever that is, that's cos(pi/8), I guess, that's the x coordinate of that line. So those are mesh points and then let me keep drawing these. So these will be mesh points too, and these will be mesh points too. So at this point, I've got probably 12 or 13 mesh points. But I've got quads, right? Well, I've got a couple of triangles here that I'm not going to touch, those are fine. But these are quads and they could be used with the Q_1 element, but I'm thinking let's stay with triangles. So I just suggested to put in triangles, put in these diagonals, keeping symmetry. And so there's the mesh. There's the mesh.

And then what does the code ask for? So I've got 13 mesh points. And the code, first of all it wants a list of the coordinates of all those mesh points. So that the code will-- And I better number them, of course. So let me number them one, shall I number them the center guys first, one, two, three, four, five, and then up here six, seven, eight, nine, now I don't know if this is a good numbering. 10, 11, 12, 13. Why don't we, just to have some consistency within plans, why don't you, you don't have to take N=4. I hope you'll take, what did I take N as four or five? Yeah, right. N is 4. But you'll want to try different N's. N=4 would be a good crude start to see what's going on, but then I hope you'll go higher and get better accuracy. And you can see how the accuracy improves, how you get closer to that, as N gets bigger. OK, we've got the nodes numbered. Oh, I better number the triangles. OK, how shall we number the triangles? Shall we do along the top, or this? I don't know. What do you want to do for the numbering of the triangles? Maybe run along the top, and then run along the bottom, because then it'll practically be a copy. So I didn't leave myself much space, but one, two, three, four, five, six, seven. Seven triangles along the top and seven along the bottom, so I have 14 triangles in this mesh. So it's a mesh with 14 triangles and 13 nodes. And I know the positions of everyone, right? I know the x,y coordinates of every one. So what the code will want is a list of those coordinates. So a list p, p will be a list of coordinates. The first guy will be-- Of the nodes so 13 rows, 3 columns. So it's a little 13 by 3 matrix that tells you where all the nodes are.

So the first one on that list would be (0,0). That's for node one. And the second one would be whatever the coordinates of that are, something, zero. (h,0), I guess it is. The third one will be (2h,0), and so on, and then complete the list of 13 positions. So you have then told the code where all the nodes are. What else do you have to tell it? Not much. You now have to tell it about the triangles. So now for every, why do I say three columns? Maybe only two, is it? You see the point already. I've forgotten, maybe, yeah. I don't see why. Well for triangles, I'm going to need three, for nodes maybe it's only got two. Maybe it's 13 by 2. I don't see why I need three. But, anyway. Then the other list is triangles. So this will be the list t, and so it takes triangle number one. Which is right here. That very first triangle. And what does it have to tell us about triangle one? The three nodes. If it tells us the three node numbers, and this list, p, gave their positions, we've got it. So how many triangles did that I have? 14? So t will be 14 by 3, and so the first guy will be just one, node number one, node number two, and node number six. That will tell us which is the first triangle. And the second triangle, I guess I've drawn from two to seven to six, right? That's a very skinny triangle up there but it's the one that started at node two, went up to seven and back to six.

So a list like that. And then the code will do the rest. I hope. Almost all the rest. The code will create the matrix K. It'll create the matrix K with, it'll be singular. Boundary conditions won't yet be in there. And then a final step after that, K, or maybe we could call it K_0, is created, a final step will be to fix u. At least at these three points. So these three will be boundary nodes. And as I say, I'm not too sure about that one, I apologize. At those boundary nodes I'm going to take the values to be zero. So this is going to be zero along the whole edge, because if it's zero there, zero there, zero there and zero there, and if it's linear, it's zero. So the final sort of subroutine in the code, the final group of commands you want to impose, zeroes here, that should then make the matrix K invertible, and then you've got KU=F to solve. So what the code is doing is creating K and F. You see the overall picture? I jumped right into this particular mesh, particular problem, but now I really should back up to where it starts. This is going to be the weak form of Laplace, Poisson, maybe I'll make a little space to put in Poisson's name too. You have a picture already of what this weak form is about, so now I'm really backing up to the start. I take the equation and I get its weak form. And remember that's in the continuous case, as it was in the quiz problem. The first step is the continuous weak form, and then the second step is choose test functions and trial-- I'm sorry, gosh, I'm in bad shape here. Because they're the same. This isn't my worst error, today. But those are trial functions, and these are test functions. OK, questions at this point, because I've-- Yeah, thank you. Good. Let's look at this picture.

AUDIENCE: [INAUDIBLE]

But did I answer that question OK? Yeah. So this is the problem I would like to know the answer to. If I have this equation, zero boundary conditions on a regular polygon with M sides, what's the answer? And it's going to be close to this, but it won't be the same. Because this does not vanish on the polygon edges. And I would like to compare the slopes, too. So the homework problem asked you not only to compare u_circle with u_(P_1), but also the slopes. The slopes here are easy, slopes here are easy because it's a bunch of flat functions. So the slopes are just constant in each triangle. OK, I'm guessing that the error gets smaller as you go in. I think that if you plot the error, it'll be largest out here and get small there. But remains to be seen. So I hope you enjoy-- yeah, good.

AUDIENCE: [INAUDIBLE]

PROFESSOR STRANG: Rather than seven.

AUDIENCE: [INAUDIBLE]

PROFESSOR STRANG: No, the middle. There's nothing magic about any particular mesh. I just chose this mesh as pretty good, and actually, I'm imagining M could get pretty big. That would be interesting. M=8 would be interesting, M=16, M=1,024, now then I'd really get interested. OK, but so if M is 1,024, then this side would be very small. Right? And I just wanted more triangles. Actually, I would like more than I've got. I'd like, if M was really big, then probably N should be at least that big. So I should have a thousand this way, if this is just a tiny bit. I just want little tiny h, and then, yeah. Actually, that might not be too bad. If M and N were roughly comparable, then that length would be roughly comparable to these lengths. And the triangles would be pretty good shapes. And that's what you're looking for. I think there's a lot of experiments to be done here. So, I'm thinking then of M and N. Here I took N to be just four. When M was eight. That's fine. But if you keep M and N roughly the same size, then you've got triangles that are not too long and skinny. I'll tell you when you might want. So generally you want nice shaped triangles. You don't want angles very small or very large, usually. But there would be, anybody in Course 16 can imagine that if I am computing the flow field past a wing, that long, thin triangles in the direction of the wing are natural. I mean, somehow a problem like true aerodynamics is by no means isotropic. I mean, the direction of the wing is kind of critical to whether the plane flies, right? So don't make the wing vertical. And if you want accuracy, then you have long, thin triangles in the direction of the flow.

But here we're not doing a flow problem. We haven't got shocks, or trailing edges, and other horrible stuff that makes planes fly. We just got Poisson's equation. OK. Thanks for those good questions, another one. AUDIENCE: [INAUDIBLE]

PROFESSOR STRANG: What would the dimension look like? Ah, would you like me to show you something about quadratics? Yeah. Shall I jump into quadratics, it's kind of fun. Quadratics, so let me just do-- So I'll come back to the weak form, right it's totally-- Oh, I'll do it now. It's so simple I don't want to forget it. The weak form, so I write the equation down, -u_xx-u_yy=f. This is the continuous weak form, equal f(x,y), OK? So that's the strong form. And I've made it the Laplacian here to keep it simple, and any right-hand side. OK, how do I get to the weak form? Just remind me, I multiply both sides by any test function v(x,y). Multiply by v(x,y), and then what do I do? I integrate over the whole region. So that's the weak form, dxdy, this is for all v, all v(x,y), all, I'll say all admissible v(x,y). So that's the weak form. If this holds for this great family of v's, the idea behind it is, that if this holds for all these trial functions, test functions, v(x,y), the only way that can happen is for this to actually equal that. That's a fundamental lemma in this part of math, and of course it has to be spelled out more than I'm doing in words. But the idea is that if these hold for such a large class of v(x,y), then the only way that can happen is for the strong form to hold. For this to actually match this. OK, so that's the start. But then what's the next step in the weak form? I like the right-hand side but I'm not so crazy about the left-hand side. I'm not crazy about it because this has second derivatives of u, and my little roof functions, pyramid functions, haven't got second derivatives. So I would be dead in the water without doing the natural step that makes everything beautiful, which is? Integration by parts.

Integrate by parts. Move derivatives off of u, onto v. One derivative onto v, off of u, so then u and v each have one derivative. I can use my piecewise linear, piecewise quadratic, all my finite elements are going to go fine. So I integrate by parts. So integrate by parts, and what is that mean in 2-D? Of course I have a double integral here. So integrate by parts, that mean you use the Green's formula. That was the key point of this Green, or Gauss-Green formula. Can I do it first in, this is -div(grad u), times v dxdy, we can write out all the terms. We can use vector notation. I could use that nabla, that upside down triangle notation, or whatever. But maybe good to see it a few different ways. So what's the point? When I integrate by parts, that minus disappears to a plus, I have a double integral then, and these derivatives move off of-- I'm taking one derivative off of here, the divergence moves over there, but when the divergence moves onto v it becomes? The transpose. It becomes gradient. And so this is gradient of u, gradient of v, dxdy. Plus boundary terms. The integral of, what is, let's see. What do I have in this integral, I have grad u dot n, times v around the boundary. And with my boundary conditions that's going to be gone, so I can come back to that.

Now, you all looked a little uncertain when I wrote Green's formula this way. For this problem I can write it more easily. This is my left side. I want to write the answer, I just want to write this weak form in a much simpler form. So let me say, what have I got here. Well, all I've got is one derivative is moving off of u and onto v. And the minus sign is disappearing, so I have du/dx times dv/dx. Right? One off of u, onto v. The other term, one y derivative, moving off of this and onto v. Minus sign again going to a plus. du/dy, dv/dy. That's the integral. That's it, that's cool. Easy to do. And on the right-hand side of course I have no change. The integral of f(x,y) v(x,y) dy. Now, that's the weak form. dxdy. Here it is, weak form. That's pretty nice. Beautifully symmetric, so the matrix that comes up when we plug in finite specific trial functions and test functions is going to be a symmetric matrix K. And the integrals are of first derivatives, so as long as our functions, our trial functions and test functions are continuous, that is, they shouldn't jump. If the trial functions or test functions jump, then if I have a jump, then the derivative would be a delta. I'd have another delta here, I'd have an integral delta, a delta times delta, and I don't want that. That's infinite. Those discontinuous elements would not be conforming, and that's a whole new world of discontinuous Galerkin. I'd have to impose penalty stuff, and Professor Peraire I mentioned, and others, Professor Darmofal in aero, are experts on this. We're doing continuous Galerkin. CG. Our piecewise linear, piecewise quadratic, they'll be continuous. All I have to do is these derivatives. Integrate those things and that's what the code will do.

OK, I've got to the weak form. That's the weak form. Now comes the finite element idea. So there is our weak form, now ready for the finite element idea. OK, so what was that idea? That's the continuous problem. Now, the finite element idea is, plug in U as a combination-- Let me write out the terms. You know what's coming here. If I'm using finite elements, I'm going to choose nice polynomials phi, say N of them. That would be like, one for every node, so I would have 13 functions here. I'm going to choose the V's to be the same as the phis. And then, I'm working then in 13 dimensions instead of infinite dimensions. So what do I do? For this limited subspace, this finite element subspace, this piecewise polynomial, piecewise linear subspace, I plug that into the weak form and I test it against 13 V's, which are phis. So I plug that in, so now what is K? Now let me just say, so I now have the integral of, yeah I guess I'd better plug it in. K_ij would then be the integral -- I'm just copying the weak form in -- Of dU/dx d-- no, sorry, I'd better just plug it in first. dU/dx*dV/dx plus dU/dy*dV/dy, those are the integrals I have to do. And on the right hand side I have to do the f-integral. fV dxdy. OK, plug that in. That's the integral over the whole domain. When I plug it in this U is a combination of known functions and the V's will be the same guys. So what am I going to get here? It's just as in 1-D. So no new ideas entering here. The new idea's going to enter when I construct these phis.

Let me just say, though, one thing. In 1-D, we pretty much had a choice of-- When it was one dimension. Just remember that. In one dimension, when I had these hat functions, when I had these guys, integrated against these guys, I pretty much had a choice of did I want to think about integrating that hat function against that one. Or actually it was their derivatives. It was the integral of U, yeah. Of phi, what I needed was all the integrals of phi_i', phi_j'. Those are what I needed, these go into K. Into the matrix K. In fact, that's what equals K_ij, the integral of phi prime. In 1-D. OK. Now, what I was going to say, I could do it this way if I wanted. But you remember the other way to do it? Was element at a time. So this was one method here. That found the entries of K separately, one by one. The other way was take the elements, one by one. So the other way was take an element like this element. It's got two functions, two trial functions are involved there. There's a little two by two, so this is four. Two by two element matrices K_e. And the quiz recalled that part. That approach. So what I want to say is that's the right way to do it in two dimensions. A triangle at a time. That's the way the code will do it. It creates these little element matrices, and then it stamps them into the big matrix K. Alright. So I want to do this integral one triangle at a time. Is the good way. OK, and that's what the code will do. Actually, I think that the best way to learn these steps is just to read the lines of the code. You can read them in the book, page 303 or something. And you'll see it just doing all the steps that need to be done. One triangle at a time.

So, now. Now comes the fun. I get to answer, what do these piecewise linear elements look like? What do the quadratic elements look like? What do the Q_1, quad elements look like? This was the golden age of finite elements, when people invented these ways to create piecewise polynomials. And it continues. People are still inventing, I had a email this week, somebody says I've got spectral elements. People are going higher and higher degrees. You, know sixth degree, eighth degree. In order to get more accuracy. OK, let's start with P_1. How do I describe a P_1 element inside a triangle? So in a triangle, the unknowns will be the value-- This has a height U_1, this has a height U_2, and a height U_3 at those nodes. Inside the triangle, the function U is linear. a+bx+cy. Then, you see that if I know these three values, then I know these three numbers. And vice versa. There's a three by three matrix, right? There has to be a three by-- Any time you see pictures like this, like the good part of 18.085 is to realize that if I have three numbers here, three values and I've got three coefficients, that there's some three by three matrix that connects them that you're going to need. That's like a meta-message of this course. Is, you've got to translate between the node values and the coefficients. Because the node values are the unknowns, right? These are the guys that are multiplying the pyramid function, this is multiplying a pyramid function with height one at that point, going down to zero. So this one will be a pyramid function of height U_2 times one, going down to zero. And U_3. So we've got a flat function in here. And it looks exactly like that. OK?

So what do I want to say? We know the positions of these three nodes from our list p, right? These were the crucial things we needed. The positions of all the nodes, we know where they are. Then there has to be a three by three matrix that will now connect to the coefficients. Why do we want the coefficients? Because those are what we do when we integrate. The coefficients are what we need, we need to integrate dU/dx, dU/dy, dU/dz. Sorry, dU/dx, dU/dy. Are you visualizing this overall solution capital U, yeah. So the overall solution capital U, you should visualize with a combination of all the little U's, is zero here and then it's going to go up in these triangles and bend around and, I don't know, maybe down again. Or maybe, no, maybe it keeps going up. This is probably the largest value, because it's the largest value in the correct solution. So is this is probably going to be the highest point of this, what's-- The Forbidden City, right? In China, in Beijing is like, or a single-- Do pagodas have flat--? No. We we will meet pagoda functions. But this would be just an ordinary western roof, I guess. Just flat pieces. Yeah. OK, see, you've got to see the whole thing and then you look at each piece. Each piece looks like that, and the integrals are doable.

OK, so while I'm going here, I want to do quadratics. You'll get the idea right away. So, same triangle, now I'm going to have quadratics. So I'm now going to have, so this won't be the arrow, this arrow will now go this way. I'm going to have d x squared, exy, and f y squared. So now how many coefficients have I got to determine a quadratic? Six, right? a, b, c, d, e, f. How many nodes do I need? Six. Where are they? Well, the natural positions are those guys in the mid-points. So now, those are all nodes now. Some nodes are at vertices of triangles, some nodes are at midpoints. But remember, we've got other triangles hooking on here, many other triangles, all with their own six nodes. Well, not their own, because they share. That's a big point. So there's a grid of triangles, with nodes for quadratic. And we've got one, two, three, four, five, six, seven, eight, nine, ten, 11, 12, 13, 14, 15, 16 nodes, I think. And within each triangle, this is what we've got. So there's a six by six matrix for each triangle. A six by six matrix which will connect the values U_1, U_2, U_3, U_4, U_5, U_6 for this triangle. Connect those six heights with these six numbers. And what will the roof look like within that triangle? Well, sort of curved. A parabola, right? A parabola somehow in 2-D, it'll look like this, yeah. Yeah. And here's the key question. Will that roof, that curvy roof, fit the one over there? Because if it didn't fit, we're in trouble. This derivative would have a delta function, and we've got delta functions, and integrals squaring them would give infinite.

So here's the question. Why does this roof, using these six points, fit onto the roof that uses U_7, U_8, U_9, and U_3, U_4 and U_5? Why do those two roofs fit together? This is the fun of piecewise polynomials? Of course, the slope will change. But the roof won't have a gap. Water won't go through it. Why's that? Do you see why? Because what do they share, what do those two curvy roofs share? They share a side. They share the same values along the side. And are those three values that are shared along the side sufficient to make it match all along the side? Yes. That's the important question. Finite elements lives or dies on that question. The answer is yes, because along that side, if I just focus on that side, where these three values are shared on both sides, by the triangle on both sides. Along that edge, what kind of a function have I got? It's second degree. This is whatever, when I restrict this to just run along a line, it's a parabola. And the parabola is determined by those three values. So having it right at three points means I have it right the whole way. Yeah. So there you see what quadratic elements would look like, and you could extend the code in the book and on the CSE site to work for quadratic elements. And you want to just guess what cubic elements could look like? I'm sorry, we've run five minutes over, but maybe finite elements is worth it. So if I had cubic elements, any idea how many? So I'm now going up to, I'm adding g x cubed, h, i, j, any idea how many coefficients I now have? Four new ones plus these six is ten. I need ten nodes. Where I am I going to put ten nodes in this triangle? I want to put them, I'd like to have some on the edges. Because the edges help me make triangles match each other. They'll just be like bowling balls. So here's six, oops, that wouldn't be believable. Is that right? Four, three, two, and one. Yeah. Yeah. OK.

So, now I've got a bubble node inside and I've got four nodes of vertices and two points, at two 1/3 points, and that will then match the triangle next to it. Because four points determine a cubic. There you go, I hope you have fun, I hope you have a great holiday. I'll see you Wednesday for Fourier and always open for questions on the MATLAB.