Explain the flow of introducing Shape Function and FEM with Linear Matrix

About

I don’t know much about how foreign books explain finite element methods, weak forms, and structural analysis using them, but frankly, Japanese textbooks are not very easy to understand. I’m not sure why, but I’ve come to think that it might be because many of them simply list the formulas that directly express structural analysis, which is an application, using examples.

For example, here’s something I wrote while reading someone else’s article. I feel like there must be a simpler explanation, but I’m not sure how to go about it.

Why is that? I’ve started to think that maybe I should write it using linear numbers a little more effectively.
What is the principle of virtual work in the weak form?
Don’t you think it’s necessary to explain the intention of the weak form as well?

The heart of weak-form

PDEs, or differential equations, require that the equations hold at certain points, but this is not suitable for numerical calculation. This is where the weak form comes in. This is where the weak form comes in, and it is integrated. In other words, it is enough to satisfy the equation on average.

This converts the point equation into the overall energy term. This is the principle of virtual work.

Shape Function with Linear Form

From this point, shape function will bring meanings.

A field u will be represented with shape function, there is

u(x)=\sum_i N_i(x) u_i

This is an expression that limits the number of degrees of freedom of the field. In other words, the field is expressed by a finite number of nodes, their coordinates, and the values ​​at those nodes.

What \nabla N represents for

\nabla u
=
\sum_i (\nabla N_i) u_i

Shape Matrix and Jacobian

Introduce M matrix, which is constructed with XY 2D coordinates, and then compute shape matrix N and the Jacobian (Assuming triangle elements).

M =
\begin{pmatrix}
1 & x_1 & y_1 \\
1 & x_2 & y_2 \\
1 & x_3 & y_3
\end{pmatrix}\\

\\

M^{-1} =
\begin{pmatrix}
a_1 & a_2 & a_3 \\
b_1 & b_2 & b_3 \\
c_1 & c_2 & c_3
\end{pmatrix}

A shape function is introduced to describe field and its derivatives

N_i(x,y) = a_i + b_i x + c_i y 

This N is represented with M.

N(x,y)
=
M^{-1}
\begin{pmatrix}
1 \\ x \\ y
\end{pmatrix}

By taking derivatives of shape function like this, PDE will be solved.

\nabla N\\

\nabla N_i =
\begin{pmatrix}
\partial N_i/\partial x \\
\partial N_i/\partial y
\end{pmatrix}
=

\begin{pmatrix}
b_i\\
c_i
\end{pmatrix}

\nabla N =
\begin{pmatrix}
b_1 & b_2 & b_3 \\
c_1 & c_2 & c_3
\end{pmatrix}

where

J =
\begin{pmatrix}
x_2-x_1 & x_3-x_1 \\
y_2-y_1 & y_3-y_1
\end{pmatrix}

Therefore,

\begin{align}
\nabla N
=
J^{-1}
\begin{pmatrix}
-1 & 1 & 0 \\
-1 & 0 & 1
\end{pmatrix}
\end{align}

But practically the M is NOT used for actual implementation. Instead, reference element is used.

In triangle case, (0,0), (1,0), (0,1) are reference nodes.

N_1(\xi,\eta) = 1 - \xi - \eta\\
N_2(\xi,\eta) = \xi \\
N_3(\xi,\eta) = \eta

what is important here is N_i(\xi_j,\eta_j) = \delta_{ij}, that is the value on an node equals 1, otherwise 0.

u(\xi,\eta) =
N_1 u_1 + N_2 u_2 + N_3 u_3

The projection from reference nodes to physicall coordinates is

x(\xi,\eta) =
N_1 x_1 + N_2 x_2 + N_3 x_3\\
y(\xi,\eta) =
N_1 y_1 + N_2 y_2 + N_3 y_3

So, the derivative of shape function is

J =
\begin{pmatrix}
\frac{\partial x}{\partial \xi} &
\frac{\partial x}{\partial \eta} \\
\frac{\partial y}{\partial \xi} &
\frac{\partial y}{\partial \eta}
\end{pmatrix}\\

\nabla_x N =
J^{-1}
\nabla_\xi N

Here,

\nabla_{\xi}\mathbf{N}
=
\begin{pmatrix}
\dfrac{\partial N_1}{\partial \xi} &
\dfrac{\partial N_2}{\partial \xi} &
\dfrac{\partial N_3}{\partial \xi}
\\
\dfrac{\partial N_1}{\partial \eta} &
\dfrac{\partial N_2}{\partial \eta} &
\dfrac{\partial N_3}{\partial \eta}
\end{pmatrix}
=
\begin{pmatrix}
-1 & 1 & 0 \\
-1 & 0 & 1
\end{pmatrix}

This is as same as (1)

Example distortion

Distortion is determined by the gradient of the shape function,

\varepsilon=\nabla u \\

\nabla u
=
\sum_i (\nabla N_i) u_i

and shape function is constant when linear shape function like this.

\nabla N_i=
\begin{pmatrix}
\partial N_i/\partial x \\
\partial N_i/\partial y
\end{pmatrix}
=
\begin{pmatrix}
b_i \\
c_i
\end{pmatrix}

therefore, stress is constant in an element.

Quadrature Loop

The integration points do not affect the shape function itself.

The only thing that matters is where you evaluate it.

• Shape function N_i → Determined by element and order (fixed)
• Integration point x_q→ Evaluation point for numerical integration

And this is how loop runs for assembling.

for element:
    X = node_coords[element]

    for qp:
        J = X.T @ gradN_ref[qp]
        measure = w[qp] * det(J)
        gradN = gradN_ref[qp] @ inv(J)
        contribution += integrand(gradN) * measure

But, there are 2 for loop.

phi[e, d, d, q]

basis[q, i]
grad_basis[q, i, d]