Finite Element

Methods of Weighted Residual to Obtain Governing Weak Form

Consider: $-(k\hat{u}')'+b\hat{u}'+c\hat{u}=f$
Residual: $r(\hat{u})=-(k\hat{u}')'+b\hat{u}'+c\hat{u}-f$
Multiply by test function v such that v and u are in the same space (thus v also satisfies EBC) and integrate over domain: $$\int_0^L r(\hat{u})vdx = 0 \rightarrow \int_0^L [-(k\hat{u}')'+b\hat{u}'+c\hat{u}-f]vdx=0 $$ Suppose v is sufficiently smooth so that we can integrate the 1st term by parts: $$\int_0^L -(k\hat{u}')'vdx=-(k\hat{u}')v|^L_0 + \int_0^Lk\hat{u}'v'dx$$ Suppose the initial boundary is zero, then u(0)=0, it follows that that v at x=0 must also be 0 and we have $$\int_0^L k\hat{u}'v'dx + \int_0^Lc\hat{u}vdx - \int_0^Lfvdx - T_vv(L)=0$$

Numerical Approximation: Galerkin Method in Weak Form

Let $\hat{u} \sim u_N = \sum_{k=1}^N \alpha_k \phi_k(x)$, then $\hat{u} \sim u_N = \sum_{i=1}^N \beta_i \phi_i(x)$ $$\int_0^L k \left( \sum_{k=1}^N \alpha_k \phi_k'(x) \right)\left( \sum_{i=1}^N \beta_i \phi_i'(x) \right)dx + \int_0^L c \left( \sum_{k=1}^N \alpha_k \phi_k(x) \right)\left( \sum_{i=1}^N \beta_i \phi_i(x) \right)dx - \int_0^L f \left( \sum_{i=1}^N \beta_i \phi_i(x) \right)dx-T_v \left( \sum_{i=1}^N \beta_i \phi_i(L) \right)dx$$ Let $$K_{kj}=\int_0^L k(x)\phi_j'(x)\phi_k'(x)dx=K_{jk}$$ $$C_{kj}=\int_0^Lc(x)\phi_j(x)\phi_k(x)dx$$ $$f_k=\int_0^L f\phi_k(x)dx+\phi_k(L)$$ The final equation to solve for is $$\sum_i B_i \left\{ \sum_k (K_{jk}+C_{jk})a_k -f_i \right\}=0 \rightarrow \left\{ \sum_k (K_{jk}+C_{jk})a_k -f_i \right\}=0 \forall B_i, \quad i=1, \dotsc ,n$$ The goal is to solve for $\alpha_k$ in $D\alpha=f$ where $D_{jk}=K_{jk}+C_{jk}$

Selecting Appropriate Basic Functions

Consider the following linear basis functions

Trial functions: $U_n(x)=\sum_{i=1}^n \alpha_i \phi_i(x)$
Test functions: $V_n(x)=\sum_{i=1}^n \beta_i X_i(x)$
How do we choose a $\phi(x)$ such that $\hat{u}$ satisfies the boundary conditions?

If we consider using quadratic basis

Local Functions: $$\left\{ \begin{array}{c} \phi_A^e(\xi)=-\frac{1}{2}(\xi)(1-\xi)\\ \phi_B^e(\xi)=\frac{1}{2}(\xi)(1+\xi)\\ \phi_C^e(\xi)=(1-xi)(1+\xi)\\ \end{array} \right.$$ If we notice that $x=x_c+\xi \frac{h}{2} \rightarrow \xi=\frac{2(x-x_c)}{2}$ $$\left\{ \begin{array}{c} \phi_A^e(\xi)=\frac{-x(x_C)}{h} \frac{h-2x+2x_C)}{h}\\ \phi_B^e(\xi)=\frac{x(x_C)}{h} \frac{h+2x-2x_C)}{h}\\ \phi_C^e(\xi)=\frac{h-2x+2x_C}{h} \frac{h+2x+2x_C)}{h}\\ \end{array} \right. \to \left\{ \begin{array}{c} \phi_A^e(\xi)=\frac{(x_C-x)(x_B-x)}{(x_C-x_A)(x_B-x_A)}\\ \phi_B^e(\xi)=\frac{(x_A-x)(x_B-x)}{(x_A-x_C)(x_B-x-C)}\\ \phi_C^e(\xi)=\frac{(x_A-x)(x_C-x)}{(x_A-x_B)(x_C-x_B)}\\ \end{array} \right.$$ Procedure: Compute $K_{kj}=\int_0^1 \bar{x}^{1/2}\phi_j'\bar{x}\phi_j'\bar{x}d\bar{x}$ after normalization $$\left[ \begin{array}{ccc} \int_0^h \sqrt{\bar{x}}\phi_1'(\bar{x})\phi_1'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_1'(\bar{x})\phi_2'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_1'(\bar{x})\phi_3'(\bar{x})d\bar{x} \\ \int_0^h \sqrt{\bar{x}}\phi_2'(\bar{x})\phi_1'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_2'(\bar{x})\phi_2'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_2'(\bar{x})\phi_3'(\bar{x})d\bar{x} \\ \int_0^h \sqrt{\bar{x}}\phi_3'(\bar{x})\phi_1'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_3'(\bar{x})\phi_2'(\bar{x})d\bar{x} & \int_0^h \sqrt{\bar{x}}\phi_3'(\bar{x})\phi_3'(\bar{x})d\bar{x} \\ \end{array} \right]$$ Compute $C_{kj}=2\int_0^1\phi_j(\bar{x})\phi_k(\bar{x})d\bar{x}$ $$ \left [ \begin{array}{ccc} \int_0^h \phi_1(\bar{x})\phi_1(\bar{x})d\bar{x} &\int_0^h \phi_1(\bar{x})\phi_2(\bar{x})d\bar{x} &\int_0^h \phi_1(\bar{x})\phi_3(\bar{x})d\bar{x} \\ \int_0^h \phi_2(\bar{x})\phi_1(\bar{x})d\bar{x} &\int_0^h \phi_2(\bar{x})\phi_2(\bar{x})d\bar{x} &\int_0^h \phi_2(\bar{x})\phi_3(\bar{x})d\bar{x} \\ \int_0^h \phi_3(\bar{x})\phi_3(\bar{x})d\bar{x} &\int_0^h \phi_3(\bar{x})\phi_2(\bar{x})d\bar{x} &\int_0^h \phi_3(\bar{x})\phi_3(\bar{x})d\bar{x} \\ \end{array} \right] $$ Solve for $\alpha_k$ in $(K_{kj}+C_{kj})\alpha_k = f_k \to \alpha_k=(K_{kj}+C_{kj})^{-1}f_k$ Approximated Solution: $\hat{u} \sim u_N=\sum_{k=1}^N \alpha_k(\bar{x})$

Discussion:

• Approximated solution using quadratic basic function converges to the analytical solution quicker than linear functions.
• We generally work with the weak form to minimize the residual
• The smaller the mesh, the more accurate $u(x,\phi)$ will be
• In the discussion of error analysis, quadratic algorithm yields a smaller error than linear due to a larger big O.

Suppose we have the schematic diagram below and the assoicated nodes.

We first note the symmetry of the membrane and simplify the model by considering only the first quadrant. Equation of Ellipse is $$(\frac{x}{a})^2+(\frac{y}{b})^2=1$$ We segment each of the quadrant into two individual meshes: Biquadratic Quadralateral and Bilinear Triangle
Rectangular Element: $Nodes [1 \;2 \;3 \;5 \;6 \;7 \;8 \;9 \;0]$
Triangular Element: $Nodes [3 \;4 \;12 \;11 \;10 \;7]$

Numerical Integration

For the integration, we use numerical method to determine the weights of the Gaussian Quadrature.
$$\int_{-1}^1 f(x)dx=\sum_{i=1}^nw_if(t_i)$$ $$(n+1)L_{n+1}(t)-(2n+1)tL_n(t)+nL_{n-1}(t)=0, L_0(t)=1, L_1(t)=t$$ $$ L_{n+1}(t)=\frac{(2n+1)tL_n(t)-nL_{n-1}(t)}{n+1} $$ For all n values between 2 and 7, we have the following table

n	$\mathbf{L_{n+1}}$
2	$L_2(t)=\frac{3tL_1(t)-L_0(t)}{2}= \frac{3t^2-1}{2}$
3	$L_3(t)=\frac{5tL_1(t)-2L_1(t)}{3}= \frac{5t^3-3t}{2}$
4	$L_4(t)=\frac{7tL_1(t)-3L_2(t)}{4}= \frac{35t^4}{8}-\frac{15t^2}{4}+\frac{3}{8}$
5	$L_5(t)=\frac{63t^5}{8}-\frac{35t^3}{4}+\frac{15t}{8}$
6	$L_6(t)=\frac{231t^6}{16}-\frac{315t^4}{16}+\frac{105t^2}{16}-\frac{5}{16}$
7	$L_7(t)=\frac{429t^7}{16}-\frac{693t^5}{16}+\frac{315t^3}{16}-\frac{35t}{16}$

Torsional Properties

Physical quantities of interest, such as shear stresses and the relationship between the twisting moment, or torque, T, and the angle of twist $\theta$, per unit length of the shaft, can be determined as follows, once the potential function u is known. $$T =2G\theta \int_\Omega ud\Omega$$ The shear stresses on the elliptical cross-section are given by the following expressions: $$\sigma_{xz}=2G\theta \frac{\partial u}{\partial y}, \; \sigma_{yz}=-2G\theta \frac{\partial u}{\partial x}$$ The stress function may be written as $$\phi=B\lbrace(\frac{x}{a})^2+(\frac{y}{b})^2-1)\rbrace$$ But since $$\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}=-2G\theta$$ We get that $$B=-\frac{a^2b^2G\theta}{a^2+b^2}$$ $$ \sigma_{xz}=\frac{\partial \phi}{\partial y}=\frac{2By}{b^2}\quad \sigma_{yz}=-\frac{\partial \phi}{\partial x}=-\frac{2Bx}{a^2} $$ $$ T =2G\theta \int_\Omega ud\Omega=-\pi B ab $$ We can also verfiy that $$\phi=B\lbrace(\frac{x}{a})^2+(\frac{y}{b})^2-1)\rbrace=\phi=-\frac{a^2b^2G\theta}{a^2+b^2}\lbrace(\frac{x}{a})^2+(\frac{y}{b})^2-1)\rbrace$$

Procedure in Calculating Torsional Rigidity

$$u_h^e=\sum_{j=1}^{N_e}u_j^e \Phi_j^e$$ $$T =2G\theta \int_\Omega u_h^e d\Omega$$

1. With $u_h^e$ from part (i), multiply this with $\Phi_j$ for quadrilateral and triangular elements from page 198 and 204, respectively. You could ignore the nodes on boundary since it will be zero. For simplicity, consider working with one quadrant.
2. Use numerical integration to compute the integral similar to step 8 from part (i)
3. Repeat Steps 1 and 2 for the other element.
4. Sum everything from both elements.
5. Assume G and $\theta$ to be 1, multiply result from step (4) by 8 (2 from the constant term of the equation, and 4 since we have four quadrants.
6. Make sure to add node 1 three times since it's shared by all quadrants.

Discussion:

• Maximum shear stress occurs at the extreme values, namely at a and b (as it approaches the boundary). Furthermore, from Figure 9, we see that greater max stress occurs at the end of the minor axis of the ellipse
• Stress along the two axes of the centerlines is symmetrical due to the geometric symmetry of the membrane. This is expected since the heaviest concentration is the middle and distributed evenly across the membrane.

Derivation of Matrices

Evaluations of $K_{ij}$
$$K_{ij}^e=\int_{\Omega_e} k \left[\frac{\partial \phi_i}{\partial x} \frac{\partial \phi_y}{\partial x}+\frac{\partial \phi_i}{\partial y} \frac{\partial \phi_y}{\partial y} \right]$$ For rectangular elements $$\phi_1(x,y)=\frac{(x-a)(y-b)}{ab}$$ $$\phi_2(x,y)=\frac{-x(y-b)}{ab}$$ $$\phi_3(x,y)=\frac{xy}{ab}$$ $$\phi_4(x,y)=\frac{-(x-a)y}{ab}$$ Evaluations of $f_i$
$$f_i^e=\int_{\Omega_e} f\phi_idxdy=\frac{1}{2A_e}\int_{\Omega_e}f(x,y)(\alpha_i +\beta_ix+\gamma_iy)dxdy$$ From the Law of Conservation Energy, the strain and potential energy must be equal and opposite from each other
$S=-W=\frac{1}{2}u^{h^T} ku^h$

Sample Output

For a fine mesh, students should see the greater yield of higher stress distribution from the centroid of the membrane.
Hint: Play with the mesh size on the element slider of the GUI. What would happen if there was a slit (cut) on the membrane?

Discussion:

• The 4 by 4 matrix elementary stiffness matrix is the same for any given n.
• Stress deflection along the two axes of the centerlines is symmetrical due to the geometric symmetry of the membrane
• Introduction of a slit on any plane causes an uneven symmetry about both axes
• Strain and potential energy derivation is consistent with the law of conservation of energy. Membrane is able to �storemore energy as we refine better meshes
• Suppose a slit cut is on the left of the membrane, because the distribution of energy originates from the center of the membrane, one can presume that the stress on the right is always greater than the left.

Linear Triangular Elements

$$u^h(x,y)=a+bx+ct$$ $$\phi_i(x,y)=\frac{1}{2A_e}[\alpha_i +\beta_ix+\gamma_iy]$$ where $\alpha_i=x_jy_k-x_ky_j$
$\beta_i=y_i-y_k$
$\gamma_i=x_k-x_j$
$A_e=\frac{\alpha}{2}$

Derivation of Matrices

Evaluations of $K_{ij}$
$$K_{ij}^e=\int_{\Omega_e} k \left[\frac{\partial \phi_i}{\partial x} \frac{\partial \phi_y}{\partial x}+\frac{\partial \phi_i}{\partial y} \frac{\partial \phi_y}{\partial y} \right]$$ since

$\frac{\partial \phi_i}{\partial x}=\frac{\beta_i}{2A_e}$

$\frac{\partial \phi_j}{\partial x}=\frac{\beta_j}{2A_e}$

$\frac{\partial \phi_i}{\partial x}=\frac{\gamma_i}{2A_e}$

$\frac{\partial \phi_j}{\partial x}=\frac{\gamma_j}{2A_e}$

We could rewrite $K_{ij}$ as
$$K_{ij}^e=\int_{\Omega_e} k \left[\frac{\beta_i}{2A_e} \frac{\beta_j}{2A_e} + \frac{\gamma_i}{2A_e} \frac{\gamma_j}{2A_e} \right]$$ Suppose k was 1 for simplicity, we then have
$$K_{ij}^e=\frac{1}{2A_e \times 2A_e}(\beta_i \beta_j + \gamma_i \gamma_j)\int_{\Omega_e} dxdy= \frac{1}{4A_e}(\beta_i \beta_j + \gamma_i \gamma_j)$$ Evaluations of $f_i$
$$f_i^e=\int_{\Omega_e} f\phi_idxdy=\frac{1}{2A_e}\int_{\Omega_e}f(x,y)(\alpha_i +\beta_ix+\gamma_iy)dxdy$$ Again, to simplify the concepts for the GUI implementation, we set f=1 to get
$f_i^e=\frac{A_e}{3}$

Potential Strain Energy and Error Analysis

From the Law of Conservation Energy, the strain and potential energy must be equal and opposite from each other
$S=-W=\frac{1}{2}u^{h^T} ku^h$ The error computation is derived using the triangle inequality $$||u^h-u^{h/2}||=||u-u^{h/2}||+||-u^h+u^h|| \leq ||u-u^{h/2}||+||u^h-u^h||$$ The energy norm is $||u^h-u^{h/2}\leq ch^{k+1}||$ in the order of 2

Discussion:

• The 3 by 3 matrix elementary stiffness matrix is the same for both upright and inverted triangle (normal and upside down triangle) for any given n.
• Thermal heat from a vertex to the middle of the opposite end looks the same for any vertex due to the geometric symmetry of the membrane.
• The contour plot of the membrane should look like a dome instead of a volcano. As we refine the meshes even more, we are more likely to detect the center node that all the lines from the vertex to the opposite end will intersect, and that node will have most concentrated thermal.
• Membrane is able to �storemore energy as we refine better meshes.

The 27-node hexahedron is the analog of the 8-node �serendipityquadrilateral
For example, the general formulas for the midside nodes are $$N_j=\frac{1}{4}(1-\xi^2)(1+\eta_j\eta)(1+\varsigma_j\varsigma)$$ $$N_j=\frac{1}{4}(1+\xi_j\xi)(1-\eta^2)(1+\varsigma_j\varsigma)$$ $$N_j=\frac{1}{4}(1+\xi_j\xi)(1+\eta_j\eta)(1-\varsigma^2)$$ Hexahedral element with tri-quadratic approximation functions
$N_1=\frac{1}{8}(1-\xi)(1-\eta)(1-\varsigma)$
$N_2=\frac{1}{8}(1+\xi)(1-\eta)(1-\varsigma)$
$N_3=\frac{1}{8}(1+\xi)(1+\eta)(1-\varsigma)$
$N_4=\frac{1}{8}(1-\xi)(1+\eta)(1-\varsigma$
$N_5=\frac{1}{8}(1-\xi)(1-\eta)(1+\varsigma)$
$N_6=\frac{1}{8}(1+\xi)(1-\eta)(1+\varsigma)$
$N_7=\frac{1}{8}(1+\xi)(1+\eta)(1+\varsigma)$
$N_8=\frac{1}{8}(1-\xi)(1+\eta)(1+\varsigma)$
$N_9=\frac{1}{4}(1-\xi^2)(1-\eta)(1-\varsigma)$
$N_{10}=\frac{1}{4}(1+\xi)(1-\eta^2)(1-\varsigma)$
$N_{11}=\frac{1}{4}(1-\xi^2)(1+\eta)(1-\varsigma)$
$N_{12}=\frac{1}{4}(1-\xi)(1-\eta^2)(1-\varsigma)$
$N_{13}=\frac{1}{4}(1+\xi)(1-\eta)(1+\varsigma^2)$
$N_{14}=\frac{1}{4}(1+\xi)(1-\eta^2)(1+\varsigma^2)$
$N_{15}=\frac{1}{4}(1-\xi^2)(1+\eta)(1+\varsigma)$
$N_{16}=\frac{1}{4}(1-\xi)(1-\eta^2)(1+\varsigma)$
$N_{17}=\frac{1}{4}(1-\xi)(1-\eta)(1-\varsigma^2)$
$N_{18}=\frac{1}{4}(1+\xi)(1-\eta^2)(1-\varsigma^2)$
$N_{19}=\frac{1}{4}(1+\xi)(1+\eta)(1-\varsigma^2)$
$N_{20}=\frac{1}{4}(1-\xi)(1+\eta^2)(1-\varsigma^2)$
$N_{21}=\frac{1}{2}(1-\xi^2)(1-\eta^2)(1-\varsigma)$
$N_{22}=\frac{1}{2}(1-\xi^2)(1-\eta^2)(1+\varsigma)$
$N_{23}=\frac{1}{2}(1-\xi^2)(1-\eta)(1-\varsigma^2)$
$N_{24}=\frac{1}{2}(1+\xi^2)(1-\eta^2)(1-\varsigma^2)$
$N_{25}=\frac{1}{2}(1-\xi^2)(1+\eta)(1-\varsigma^2)$
$N_{26}=\frac{1}{2}(1-\xi)(1-\eta^2)(1-\varsigma^2)$
$N_{27}=(1-\xi^2)(1-\eta^2)(1-\varsigma^2)$

Derivatives of Shape Functions \begin{eqnarray} \frac{\partial N_i^{(e)}}{\partial x}=\frac{\partial N_i^{(e)}}{\partial \xi}\frac{\partial \xi}{\partial x}+\frac{\partial N_i^{(e)}}{\partial \eta}\frac{\partial \eta}{\partial x}+\frac{\partial N_i^{(e)}}{\partial \varsigma}\frac{\partial \varsigma}{\partial x}\\ \frac{\partial N_i^{(e)}}{\partial y}=\frac{\partial N_i^{(e)}}{\partial \xi}\frac{\partial \xi}{\partial y}+\frac{\partial N_i^{(e)}}{\partial \eta}\frac{\partial \eta}{\partial y}+\frac{\partial N_i^{(e)}}{\partial \varsigma}\frac{\partial \varsigma}{\partial y}\\ \frac{\partial N_i^{(e)}}{\partial z}=\frac{\partial N_i^{(e)}}{\partial \xi}\frac{\partial \xi}{\partial z}+\frac{\partial N_i^{(e)}}{\partial \eta}\frac{\partial \eta}{\partial z}+\frac{\partial N_i^{(e)}}{\partial \varsigma}\frac{\partial \varsigma}{\partial z} \end{eqnarray} The infinitesimals $d\xi$, $d\eta$ and $d\varsigma$ transform into $dx$, $dy$, and $dz$ can be written in matrix form as $$ \left[ \begin{array}{c} dx \\ dy \\ dz \end{array} \right]= \left[ \begin{array}{ccc} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} & \frac{\partial x}{\partial \varsigma}\\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} & \frac{\partial y}{\partial \varsigma}\\ \frac{\partial z}{\partial \xi} & \frac{\partial z}{\partial \eta} & \frac{\partial z}{\partial \varsigma} \end{array} \right] \left[ \begin{array}{c} d\xi \\ d\eta \\ d\varsigma \end{array} \right] $$ Jacobian Matrix $$ \mathbf{J}=\frac{\partial (x,y,z)}{\partial (\xi,\eta,\varsigma)}= \left[ \begin{array}{ccc} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} & \frac{\partial z}{\partial \xi}\\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} & \frac{\partial z}{\partial \eta}\\ \frac{\partial x}{\partial \varsigma} & \frac{\partial y}{\partial \varsigma} & \frac{\partial z}{\partial \varsigma} \end{array} \right] $$ The isoparametric definition of hexahedron element geometry is $$ x=\sum_{i=1}^{20}x_iN_i^{(e)} \qquad y=\sum_{i=1}^{20}y_iN_i^{(e)} \qquad z=\sum_{i=1}^{20}z_iN_i^{(e)} $$ Whenever $|J|\neq 0$, we can write $$ \left[ \begin{array}{ccc} \frac{\partial \xi}{\partial x} & \frac{\partial \xi}{\partial y} & \frac{\partial \xi}{\partial z}\\ \frac{\partial \eta}{\partial x} & \frac{\partial \eta}{\partial y} & \frac{\partial \eta}{\partial z}\\ \frac{\partial \varsigma}{\partial x} & \frac{\partial \varsigma}{\partial y} & \frac{\partial \varsigma}{\partial z}\\ \end{array} \right] = \mathbf{J^{-1}} $$ The derivatives of $\Phi_j^e$ are obtained by the chain rule: \begin{eqnarray} \frac{\partial \Phi_j^e}{\partial x}=\frac{\Phi_j}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\Phi_j}{\partial \eta} \frac{\partial \eta}{\partial x}+\frac{\Phi_j}{\partial \varsigma} \frac{\partial \varsigma}{\partial x}\\ \frac{\partial \Phi_j^e}{\partial y}=\frac{\Phi_j}{\partial \xi} \frac{\partial \xi}{\partial y}+\frac{\Phi_j}{\partial \eta} \frac{\partial \eta}{\partial y}+\frac{\Phi_j}{\partial \varsigma} \frac{\partial \varsigma}{\partial y}\\ \frac{\partial \Phi_j^e}{\partial z}=\frac{\Phi_j}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\Phi_j}{\partial \eta} \frac{\partial \eta}{\partial z}+\frac{\Phi_j}{\partial \varsigma} \frac{\partial \varsigma}{\partial z}\\ \end{eqnarray} Numerical Integration Over Hexahedral, the total number of Gauss points is $p^3$ $$ \int_{-1}^{1}\int_{-1}^{1}\int_{-1}^{1}F(\xi,\eta,\varsigma)d\xi d\eta d\varsigma \approx \sum_{i=1}^{p_1}\sum_{j=1}^{p_2}\sum_{k=1}^{p_3} w_iw_jw_k F(\xi_i,\eta_j,\varsigma_k) $$

Bending of a Uniform, Homogeneous Elastic Beam (Euler-Bernoulli Theory)

Strong Form: $$\int_0^L\lbrace (EIw'')''-q\rbrace vdx-\lbrace -(EIw'')(L)-M_L\rbrace v'(L)+\lbrace (EIw'')(L)-v_L\rbrace v(L)=0 \quad \forall \;v \quad s.t \; v(0)=0, \; v'(0)=0$$ Weak Form: $$\int_0^L EIw''v''dx-\int_0^L qvdx+M_Lv'(L)-v_Lv(L)=0 \quad \forall \; v \quad s.t \; v(0)=0, \; v'(0)=0$$ Principal of Virtual Work $$\int_0^L EIw''v''dx=\int_0^L qvdx-M_Lv'(L)+v_Lv(L)$$ But since we don't have any applied moment of load at the supported end,our governing equation is actually $$\int_0^L EIw''v''dx=\int_0^L qvdx$$ The simplest Bernoulli-Euler plane beam element with two end nodes has four degrees of freedom $\mathbf{u^e}=[v_1 \theta_1 v_2 \theta_2]^T$.
Shape functions for this problem are conveniently expressed in terms of the dimensionless coordinate $$\xi=\frac{2x}{h}-1 \qquad \frac{dx}{d\xi}=\frac{1}{2}h \qquad \frac{d\xi}{dx}=\frac{2}{h}$$ \begin{eqnarray} N_1(\xi)=\frac{1}{4}(1-\xi)^2(2+\xi)\\ N_2(\xi)=\frac{h}{8}(1-\xi)^2(1+\xi)\\ N_3(\xi)=\frac{1}{4}(1+\xi)^2(2+\xi)\\ N_4(\xi)=\frac{h}{8}(1+\xi)^2(\xi-1)\\ \end{eqnarray} Local Stiffness Matrix $$K_{ij}=\int_eEI\frac{d^2N_i(x)}{dx^2} \frac{d^2N_j(x)}{dx^2}dx=\int_{-1}^1 EI\frac{d^2N_i(x)}{dx^2} \frac{d^2N_j(x)}{dx^2}\frac{1}{2}h d\xi$$ For a generic case with both free nodes, the matrix becomes \[ K^e = \frac{EI}{h^3}\left[ \begin{array}{cccc} 12 & 6h & -12 & 6h \\ 6h & 4h^2 & -6h & 2h^2 \\ -12 & -6h & 12 & -6h \\ 6h & 2h^2 & -6h & 4h^2 \end{array} \right]\] Local Force Vector $$f_i^{(e)}=\int_e q(x)N_i(x)dx=\int_{-1}^1 q(x)N_i(x)\frac{1}{2}h d\xi$$ For uniform load $q_0$ $$f^{(e)}=\frac{q_0h_e}{12} \left\{ \begin{array}{c} 6\\ h_e \\ 6 \\ -h_e \end{array} \right\} $$ System array is defined with dimensions $n\;x\;m$ where $n$ is 4 (the number of degree of freedom per element) and $m$ is 3 (the number of elements).

Distribution of Bending Moment Along the Beam

Since moment is related to the distributed load by its second derivative, the approximate solution is $$M(x)=-EI\frac{d^2w}{dx^2}=-EI\sum_{j=1}^4 u_j^e \frac{d^2 \Phi_j^e}{dx^2}$$

Reactions at the Beam Ends

$$V(x)=\frac{dM}{dx}=-\frac{d}{dx}\left(EI\frac{d^2w}{dx^2}\right)=-EI\sum_{i=1}^4 u_j^e \frac{d^3 \Phi_j^e}{dx^3}$$