Basic explanation of shape function

well this is my first time posting here. I just started studying FEM in a more structured basis compared to what I used to do during my undergraduate courses. I'm doing this because, despite the fact that I can "use" the "FEM tool" in commercial (and other not-so-commercial) software, I would like to really-really understand the underground techniques that support the method. That's why I'm coming here with such, at least for the experienced user of the technique, stupid and basic question.

Now I'm reading a quite popular (I think) and "engineer-friendly" book called "Finite element method- The basics" from Zienkwicz. I've been reading this book from the first page but I yet can't understand the concept of shape function in the way Zienkwicz explains it.

What I know about from the things that I'd read is that a "Stifness" matrix, the one that relates the unknowns with the result (A in: A*k=b), has his components from the "relationships between the nodes", and if that "relationship" changes, (i.e. if we change it to a Higher order interpolant), that stifness matrix changes, because the relationship between the nodes does.

But in this book, the definition is quite fuzzy for me, because in some point it says that you can arbitrarily chose the function as, i.e. the identity matrix:

The only explanation I found is in this blog , but it is still not so clear for me. So, somebody can give me a simple plain explanation of what is a Shape functon and how it is done to "put it" in the stifness matrix?

-
It is basically some approximation of how the solution varies within the element. E.g., in a single, small 1D linear element (think of a bar/rod) it is reasonable to assume that as you move from one end to the other, the solution (e.g., displacement or temperature) would vary linearly. Off course you choose higher order polynomials (instead of linear) as well. Hope that helps. –  stali Feb 2 at 16:45
Im not going to try to answer this directly, but the best explanation I have seen of FEM is (ironically) in a book on mesh-free methods if you can get ahold of it. Meshfree methods : moving beyond the finite element method / G.R. Liu. –  Nick Feb 3 at 2:02

I've always found the approach to describing finite element methods that focuses on the discrete linear system and works backward unnecessarily confusing. It is much clearer to go the other way, even if that involves a bit of mathematical notation in the beginning (which I'll try to keep to a minimum).

Assume that you are trying to solve an equation $A u = f$ for given $f$ and unknown $u$, where $A$ is a linear operator that maps functions (e.g., describing the displacement at every point $(x,y)$ in a domain) in a space $V$ to functions in another space (e.g., describing the applied forces). Since the function space $V$ is usually infinite-dimensional, this system cannot be solved numerically. The standard approach is therefore to replace $V$ by a finite-dimensional subspace $V_h$ and look for $u_h \in V_h$ satisfying $Au_h = f$. This is still infinite-dimensional due to the range space, so we just ask for the residual $Au_h-f$ to be orthogonal to $V_h$ - or equivalently, that $v_h^T(Au_h-f) = 0$ for every basis vector $v_h$ in $V_h$. If we now write the $u_h$ as a linear combination of these basis vectors, we are left with a linear system for the unknown coefficients in this combination. (The terms $v_i^TAu_j$ are exactly the entries of the stiffness matrix $K_{ij}$, and $v_j^Tf$ are the entries of the load vector. If $A$ is a differential operator, one usually performs integration by parts at some point, but this is not important here.)

None of this so far is specific to finite element methods, but applies to any so-called Galerkin method or method of weighted residuals. The finite element method is characterized by a special choice of $V_h$: The computational domain is decomposed into a number of elements of the same basic shape (e.g., triangles; the process is often called triangulation), and the space $V_h$ is chosen such that restricted to each element, functions in $V_h$ are polynomials (e.g, linear in $x$ and $y$). Furthermore, the basis functions are chosen such that they are non-zero only in (the neighborhood of) one of the elements. The point of this choice is that you can build such a basis of $V_h$ fairly easily by finding a basis $\{\psi_j\}$ of the polynomial space on a single reference element (such as the triangle with vertices $(0,0)$, $(0,1)$ and $(1,0)$) and then using an affine transformation to map these basis functions to basis functions on each element in the triangulation. These $\psi_j$ are the shape functions. Usually, one requires that the local basis functions take the value $1$ at only one of the vertices and $0$ at the others (called a nodal basis), which is what the page you linked is talking about.

(Other choices of $V_h$ lead to other methods; in fact, there are spectral methods where the basis functions are chosen such that the stiffness matrix is the identity. Of course, there's no free lunch, so other parts of the procedure become more difficult with this basis.)

-

My take is in lecture 4 at http://www.math.tamu.edu/~bangerth/videos.html . In particular, it gives you an idea of why we choose the hat functions we usually use when we use the finite element method -- namely, because they lead to the important concept of sparsity, even though many other choices of basis functions would have been equally valid.

-

The most important thing to know about "shape functions" is that they describe how the dependent variable(s) you want to calculate (e.g. displacement) vary as a function of the spatial coordinates of the element (e.g. x and y) in terms of some unknown scalar parameters.

Often the shape functions are simple polynomials and the scalar parameters are the values of the dependent variables at the element nodes.

Forming the finite element equations using these shape functions requires a few other fundamental concepts such as establishing a "weak form" of the partial differential equation you are trying to solve.

There is a lot of unnecessary "mysticism" associated with the finite element method so I encourage your approach of trying to get a thorough understanding of the fundamentals.

-

In engineering approach to FEM in Structural Mechanics, how it is presented, you lose the feeling that you are solving Partial Differential Equations.

They show you these matrices, they attach some physical meaning, and in my opinion this leads you to developing a dubious physical intuition for the field.

It may be helpful to think about the subject it terms of geometry. The solution to a boundary value problem for PDE is some shape. V.I. Arnol'd once said praising Newton's accomplishments in the field, to paraphrase - he did a marvelous thing by creating the field of differential equations by allowing us to reformulate the problems of natural sciences to geometrical problems of curves in plane and surfaces in space.

In FEM you approximate the solution (in FD and FVM you approximate the governing equation).

Enter Boris Gligorievich Galerkin. What did B.G. Galerkin say?

He said: “I want you, not to be able to make residual with the same basis functions, you used to create the solution.

(P.S. This story is completely not true, and I urge my readers to find a better explanation of (Bubnov-) Galerkin method, if it exists.)

Basis functions, or trial functions are those that you use to build the solution. You use them to approximate the shape of the solution.

Galerkin wants you to think of functions in abstract way analogous to vectors in ordinary Euclidean space, and using the notion of orthogonality to create enough algebraic conditions to find the necessary parameters to build the solution. These algebraic conditions when all put together create $Ku=f$.

FEM is when you approximate shape of the solution function peace by peace. In each peace - the element-you have some basic shapes (shape functions) that have some flexibility and can approximate various solutions - but only one is the solution to your problem. There are shape functions as much as there are nodes. I-th shape function is equal to one in i-th node, in other nodes it is zero. Solving $Ku=f$ we find parameters that define shape functions - their values at nodes, the mix of shape functions in an element gives the solution over that element, collecting the shapes peace by peace over elements we get global shape-which is a sought solution to our BVP.

-