# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

import theme
theme.load_style()

# # The Euler-Bernoulli Beam Element
#
# <a rel="license" href="http://creativecommons.org/licenses/by/4.0/"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by/4.0/80x15.png" /></a>
#
# This lecture by <NAME> is licensed under the
# <a rel="license" href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution 4.0 International License</a>.  All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT).

# # Topics
# - [Introduction](#intro)
# - [Beam Element](#beam-el)
#   - [$Discretization of the Domain](#disc)
#   - [Derivation of Element Equations](#derive)

# # <a id='intro'></a> Introduction
#
# In the last two lectures we discussed the finite element formulation of one dimensional second order BVPs and we completed examples in thermal, fluid, and solid mechanics.  Today we discuss the finite element formulation of the one dimensional fourth order BVP associated with the Euler-Bernoulli beam theory.  The formulation of the finite element equations is obtained by following the same steps as in the second order BVP finite element formulation.

# # <a id='beam-el'></a> Euler-Bernoulli Beam Element
#
# In Euler-Bernoulli beam theory, it is assumed that plane sections remain
# plane. In this theory, the transverse deflection $u$ of the beam is governed
# by
#
# $$
# \frac{d^2}{dx^2}\left( a\frac{d^2u}{dx^2}\right) +q(x)=0
# $$
#
# In the figure below, a typical beam, with sign conventions is shown
#
# <img src='./Images/beam.d/beam-detail.png' style='width:80%'/>
#
# Here $a=a(x)$ and $q=q(x)$ are the given problem data and $u$ is the dependent
# variable. The function $a=EI$ is the flexural rigidity of the beam, $q$ is the
# transversely distributed load, and $u$ is the transverse deflection. In
# addition to satisfying the beam bending equation, $u$ must also satisfy the
# boundary conditions. Since beam bending is governed by a fourth order DE, four
# boundary conditions are required to solve it. As in the second order DEs we
# previously solved, the appropriate boundary conditions arise from the weak
# form of the governing equation. We now develop the finite element equations
# for a beam.
#
# ## <a id='disc'></a> Discretization of the Domain
#
# The domain (length) of the beam is divided into $N-1$ line elements, each having at least two end nodes, as shown below.  Note that the element geometry is the same as with the bar element, however, the number and form of the primary and secondary unknowns at each node is different and are dictated by the variational (weak) formulation of the beam equation.
#
# <img src='./Images/beam.d/beam-disc.png'/>
#
# ## <a id='derive'></a> Derivation of Element Equations
#
# We now isolate a typical element $\Omega^{(e)}=(x_A,x_B)$ and construct the
# weak form of the beam bending equation over that element. Recall that the weak
# form provides the primary and secondary variables of the problem. We will then
# choose approximation for the primary variable and derive the element
# equations.

# ### The Weak Form
#
# We derive the weak form in the same way as with the second order BVP, namely,
# we multiply the governing equation by a weight function $w(x)$, integrate by
# parts to reduce the continuity requirements on $u$, and evaluate the resulting
# boundary conditions.
#
# $$
# \begin{gather}
# \int_{x_A}^{x_B}\left(\frac{d^2}{dx^2}\left(a\frac{d^2u}{dx^2} \right) + q \right) w dx = 0 \\
# \int_{x_A}^{x_B}\left( -\frac{dw}{dx}\frac{d}{dx}\left( a \frac{d^2u}{dx^2}\right) + wq \right) dx + \left[ w\frac{d}{dx}\left( a\frac{d^2u}{dx^2} \right) \right]_{x_A}^{x_B} =0 \\
# \int_{x_A}^{x_B}\left( a\frac{d^2w}{dx^2}\frac{d^2u}{dx^2}+wq\right) dx + \left[ w\frac{d}{dx}\left( a\frac{d^2u}{dx^2} \right)-\frac{dw}{dx}a\frac{d^2u}{dx^2} \right]_{x_A}^{x_B} =0
# \end{gather}
# $$
#
# Note that, because of the two integrations by parts, $u,w\in C^1$. We have
# reduced the continuity requirements on $u$ from $C^3$ to $C^1$. Also, because
# of the integration by parts, there are now two boundary terms to be evaluated
# at $x=x_A$ and $x=x_B$. Let's look closely at the two boundary conditions:
#
# $$
# \left[ w\frac{d}{dx}\left( a\frac{d^2u}{dx^2} \right)-\frac{dw}{dx}a\frac{d^2u}{dx^2} \right]_{x_A}^{x_B}
# $$
#
# From the boundary conditions, it follows that the specification of $u$ and
# $du/dx$ at the endpoints of the elements constitute the essential boundary
# conditions and the specification of
#
# $$
# \frac{d}{dx}\left( a\frac{d^2u}{dx^2} \right)=V \quad (\text{shear force})
# $$
#
# and
#
# $$
# a\frac{d^2u}{dx^2} =M \quad (\text{bending moment})
# $$
#
# at the endpoints constitute the natural boundary conditions. Thus, their are
# two essential and two natural boundary conditions. Therefore, we identify $u$
# and $du/dx$ as the primary variables at each node. The natural boundary
# conditions always remain in the weak form and end up as the source vector of
# the matrix equations. To simplify notation, we introduce the the notation
# $\phi=du/dx$ and
#
# $$\begin{gather}
# Q_1^{(e)}=V_1^{(e)}=\left[ \frac{d}{dx}\left( a\frac{d^2u}{dx^2} \right) \right]\Bigg|_{{x_A}}, \quad  Q_2^{(e)}=M_1^{(e)}=\left( a\frac{d^2u}{dx^2} \right) \Bigg|_{x_A}\\
# Q_3^{(e)}=V_2^{(e)}=-\left[ \frac{d}{dx}\left( a\frac{d^2u}{dx^2} \right) \right]\Bigg|_{x_B}, \quad  Q_4^{(e)}=M_2^{(e)}=-\left( a\frac{d^2u}{dx^2} \right) \Bigg|_{x_B}\\
# \end{gather}
# $$
#
# where $Q_1^{(e)}$ and $Q_3^{(e)}$ denote the shear forces, and $Q_2^{(e)}$ and
# $Q_4^{(e)}$ denote the bending moments shown in \rfig{bd0}. Since the set
# $Q_i^{(e)}$ contain both point and ''bending'' forces, we will call the set
#
# $$
# \begin{Bmatrix} Q_1^{(e)} \\ Q_2^{(e)}\\ Q_3^{(e)}\\ Q_4^{(e)} \end{Bmatrix}
# $$
#
# the *generalized forces*. We call the corresponding set of displacements and
# rotations the *generalized displacements*:
#
# $$
# \begin{Bmatrix} u(x_A) \\ u'(x_A) \\ u(x_B) \\ u'(x_B) \end{Bmatrix} =
# \begin{Bmatrix} u_1^{(e)} \\ u_2^{(e)} \\ u_3^{(e)} \\ u_4^{(e)}\end{Bmatrix} =
# \begin{Bmatrix} u_1 \\ \phi_1 \\ u_2 \\ \phi_2\end{Bmatrix}
# $$
#
# With this notation, we can now write
#
# $$
# \int_{x_A}^{x_B}\left( b\frac{d^2w}{dx^2}\frac{d^2u}{dx^2}+qw\right) dx + \left[ wV-\frac{dw}{dx}M\right]\Bigg|_{x_A}^{x_B} =0
# $$

# ### Interpolation Functions
# The weak form requires that the interpolation functions of an element be continuous with nonzero derivatives up to order two.  The approximation of the primary variables over a finite element should be such that it satisfies the essential boundary conditions of the element:
#
# $$
# u(x_A)=u_1, \quad u(x_B)=u_2, \quad \phi(x_A)=\phi_1, \quad \phi(x_B)=\phi_2
# $$
#
# Since there is a total of four conditions in an element (two per node), a four parameter polynomial must be selected for $u$:
#
# $$
# u(x)=u_0 + u_1x + u_2x^2 + u_3x^3
# $$
#
# We now express the $c_i$ in terms of the general displacements and require that
#
# $$
# u(x_A)=u_1^{(e)}=u_0 + u_1x_A+u_2x_A^2+u_3x_A^3\\
# \phi(x_A)= u_2^{(e)}=u_1 + 2u_2x_A+3u_3x_A^2\\
# u(x_B)=u_3^{(e)}=u_0 + u_1x_B+u_2x_B^2+u_3x_B^3\\
# \phi(x_B)= u_4^{(e)}=+u_1 + 2u_2x_B+3u_3x_B^2
# $$
#
# or
#
# $$
# \begin{Bmatrix} u_1^{(e)} \\ u_2^{(e)} \\ u_3^{(e)} \\ u_4^{(e)} \end{Bmatrix} =
# \begin{bmatrix} 1 & x_A & x_A^2 & x_A^3 \\ 0 & 1 & 2x_A & 3x_A^2 \\ 1 & x_B & x_B^2 & x_B^3 \\ 0 & 1 & 2x_B & 3x_B^2\end{bmatrix}
# \begin{Bmatrix} u_0 \\ u_1 \\ u_2 \\ u_3 \end{Bmatrix}
# $$
#
# Inverting this expression to express the $u_i$ in terms of the $u_i^{(e)}$ and substituting the result gives
#
# $$
# u^{(e)}(x)=u_1^{(e)}N_1^{(e)}+u_2^{(e)}N_2^{(e)}+u_3^{(e)}N_3^{(e)}+u_4^{(e)}N_4^{(e)}=\sum_{j=1}^4u_j^{(e)}N_j^{(e)}
# $$
#
# where
#
# $$
# N_1^{(e)}=1-3\left(\frac{x-x_A}{h_e}\right)^2+2\left(\frac{x-x_A}{h_e}\right)^3, \ \ \ N_2^{(e)}=\left( x-x_A\right)\left(  \frac{x-x_A}{h_e}-1\right)^2\\
# N_3^{(e)}=3\left(\frac{x-x_A}{h_e}\right)^2-2\left(\frac{x-x_A}{h_e}\right)^3, \ \ \ N_4^{(e)}=\left( x-x_A\right)\left(\left(\frac{x-x_A}{h_e}\right)^2-\frac{x-x_A}{h_e}\right)
# $$
#
# Note that the cubic interpolation functions were derived by interpolating $u$
# and its derivatives at the nodes. Interpolation polynomials derived in this
# way are known as Hermite interpolation functions.
#
# The shape functions $N_i^{(e)}$ can be expressed in terms of the local coordinates $\hat{x}$:
#
# $$
# N_1^{(e)}=1-3\left(\frac{\hat{x}}{h_e}\right)^2+2\left(\frac{\hat{x}}{h_e}\right)^3, \ \ \ N_2^{(e)}=\hat{x}\left( \frac{\hat{x}}{h_e}-1\right)^2\\
# N_3^{(e)}=3\left(\frac{\hat{x}}{h_e}\right)^2-2\left(\frac{\hat{x}}{h_e}\right)^3, \ \ \ N_4^{(e)}=\hat{x}\left(\left(\frac{\hat{x}}{h_e}\right)^2-\frac{\hat{x}}{h_e}\right)
# $$
#
# The Hermite cubic interpolation functions satisfy the following interpolation properties ($i,j=1,2$)
#
# $$
# N_{2i-1}^{(e)}(\hat{x}_j)=\delta_{ij}, \ \ \ \ N_{2i}^{(e)}(\hat{x}_j)=0, \ \ \ \ \sum_{i=1}^2N_{2i-1}^{(e)}=1\\
# \frac{dN_{2i-1}^{e}}{dx}\Bigg|_{\hat{x}_j}=0, \ \ \ \ \frac{dN_{2i}^{(e)}}{dx}\Bigg|_{\hat{x}_j}=\delta_{ij}
# $$
#
# where $\hat{x}_1=0$ and $\hat{x}_2=h_e$ are the local coordinates of nodes 1 and 2 of the element $\Omega^{(e)}=(x_A,x_B)$.\\

# ### Finite Element Model
#
# We now substitute
# $$
# u^{(e)}(x)=\sum_{j=1}^4 u_j^{(e)}N_j^{(e)}
# $$
# and
# $$
# w=\sum_{i=1}^4 w_i^{(e)}N_i^{(e)}
# $$
#
# Substituting and simplifying gives
#
# $$
# \sum_{j=1}^4\left(\int_{x_A}^{x_B}a\frac{d^2N_i^{(e)}}{dx^2}\frac{d^2N_{j}^{(e)}}{dx^2}dx\right) u_j^{(e)} + \int_{x_A}^{x_B}N_i^{(e)}qdx-Q_i^{(e)}=0
# $$
#
# or
#
# $$
# \sum_{j=1}^4k_{ij}^{(e)}u_j^{(e)}+F_i^{(e)}=0
# $$
# where
#
# $$
# k_{ij}^{(e)}=\int_{x_A}^{x_B}a\frac{d^2N_i^{(e)}}{dx^2}\frac{d^2N_j^{(e)}}{dx^2}dx, \ \ \ \ F_i^{(e)} = \int_{x_A}^{x_B}N_i^{(e)}qdx-Q_i^{(e)}
# $$
#
# In matrix notation this can be written
#
# $$
# \begin{bmatrix} 
# k_{11}^{(e)} & k_{12}^{(e)} & k_{13}^3 & k_{14}^{(e)} \\ 
# k_{21}^{(e)} & k_{22}^{(e)} & k_{23}^3 & k_{24}^{(e)} \\ 
# k_{31}^{(e)} & k_{32}^{(e)} & k_{33}^3 & k_{34}^{(e)} \\ 
# k_{41}^{(e)} & k_{42}^{(e)} & k_{43}^3 & k_{44}^{(e)} 
# \end{bmatrix} 
# \begin{Bmatrix} 
# u_1^{(e)} \\  u_2^{(e)} \\ u_3^{(e)} \\ u_4^{(e)} 
# \end{Bmatrix} = 
# -\begin{Bmatrix} q_1^{(e)} \\  q_2^{(e)} \\ q_3^{(e)} \\ q_4^{(e)} \end{Bmatrix} + \begin{Bmatrix} Q_1^{(e)} \\  Q_2^{(e)} \\ Q_3^{(e)} \\ Q_4^{(e)} \end{Bmatrix}
# $$
#
# For the case in which $a=EI$ and $w$ are constant over an element, the coefficients of each element equation are found as with the bar element:
#
# $$
# k_{11}^{(e)}=EI\int_0^{h_e}\frac{d^2N_1^{(e)}}{d\hat{x}^2}\frac{d^2N_1^{(e)}}{d\hat{x}^2} \ d\hat x=EI\int_0^{h_e}\left(-\frac{6}{h_e^{2}}+12\,{\frac {\hat x}{h_e^{3}}}\right)^2 \ d\hat x=\frac{12EI}{h_e^3}
# $$
#
# $$
# F_1^{(e)}=q\int_0^{h_e}N_1^{(e)} \ d\hat x - Q_1^{(e)} =q\int_0^{h_e}\left(-\frac{6}{h_e^{2}}+12\,{\frac {\hat x}{h_e^{3}}}\right) d\hat x - V_1^{(e)} = \frac{qh_e}{2} - V_1^{(e)}
# $$
#
# All other coefficients can be found in a similar manner and combine to give the beam equations for a single two node element:
#
# $$
# \frac{2EI}{h_e^3}\begin{bmatrix}
# 6 & 3h_e & -6 & 3h_e \\
# 3h_e & 2h_e^2 & -3h_e & h_e^2 \\
# -6 & -3h_e & 6 & -3h_e \\
# 3h_e & h_e^2 & -3h_e & 2h_e^2 \end{bmatrix}\begin{Bmatrix} u_1^{(e)} \\ \phi_1^{(e)} \\ u_2^{(e)} \\ \phi_2^{(e)} \end{Bmatrix}=-\frac{wh_e}{12}\begin{Bmatrix}6\\h_e\\6\\-h_e\end{Bmatrix}+\begin{Bmatrix}V_1^{(e)}\\M_1^{(e)}\\ V_2^{(e)}\\M_2^{(e)}\end{Bmatrix}
# $$

# ### Assembly of Global Equations
#
# We assemble the global equations in the same way as the bar element equation, except we take into account the increased degrees of freedom at element nodes.  The assembly of the global equations is based on:
#
#
# - Interelement continuity of primary variables (deflection and slope)
# - Interelement equilibrium of the secondary variables (shear force and bending moment) at the nodes of common elements
#
# We demonstrate the assembly of the global equations by examining the following two element model
#
# <img src='./Images/beam.d/beam-example.png' style='width:70%'/>
#
# The continuity of the primary variables implies the following relation between the element degrees of freedom $u_i^{(e)}$ and the global degrees of freedom $u_i$ and $\phi_i$:
#
# $$ \begin{matrix}
# u_1^1=u_1 & u_2^1=\phi_1 & u_3^1=u_1^2=u_2 \\
# u_4^1=u_2^2=\phi_2 &  u_3^2=u_3 & u_4^2=\phi_3 \end{matrix}
# $$
#
# The equilibrium of the forces and moments at the connecting node requires that
# $$
# Q_3^1+Q_1^2={\rm applied \ external \ point \ force}=0
# $$
#
# $$
# Q_4^1+Q_2^2={\rm applied \ external \ bending \ moment}=0
# $$
#
# If no external applied forces are given, the sum should be equated to zero.  In equating the sums to the applied generalized forces (i.e. the force or moment) the sign convention for the element force degrees of freedom should be followed.  Forces are taken positive acting upwards and moments positive taken clockwise.\\
#
# We now expand each element equation to the size of the global equation and combine them to form the global set of equations:
#
# ### Element 1
# $$
# \begin{bmatrix} k_{11}^1 & k_{12}^1 & k_{13}^1 & k_{14}^1 & 0 & 0\\
#  		          k_{21}^1 & k_{22}^1 & k_{23}^1 & k_{24}^1 & 0 & 0\\
# 			k_{31}^1 & k_{32}^1 & k_{33}^1 & k_{34}^1 & 0 & 0\\
# 			k_{41}^1 & k_{42}^1 & k_{43}^1 & k_{44}^1  & 0 & 0\\
# 			0 & 0 & 0 & 0  & 0 & 0\\
# 			0 & 0 & 0 &0  & 0 &0\\
# \end{bmatrix}
# \begin{Bmatrix}u_1 \\ \phi_1 \\ u_2 \\ \phi_2 \\ u_3 \\ \phi_3 \end{Bmatrix}
# =-\frac{wh_e}{12}\begin{Bmatrix}6\\h_e\\6\\-h_e\\0\\0\end{Bmatrix}+\begin{Bmatrix}Q_1^1\\Q_2^1\\Q_3^1\\Q_4^1\\0\\0\end{Bmatrix}
# $$
#
# ### Element 2
# $$
# \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0\\
#  		           0 & 0 & 0 & 0 & 0 & 0\\
# 			 0 & 0 & k_{11}^2 &k_{12}^2  & k_{13}^2 & k_{14}^2\\
# 			 0 & 0& k_{21}^2 &k_{22}^2  & k_{23}^2 & k_{24}^2\\
# 			0 & 0 & k_{31}^2 & k_{32}^2  & k_{33}^2 & k_{34}^2\\
# 			0 & 0 & k_{41}^2 & k_{42}^2  & k_{43}^2 & k_{44}^2\\
# \end{bmatrix}\begin{Bmatrix}u_1 \\ \phi_1 \\ u_2 \\ \phi_2 \\ u_3 \\ \phi_3 \end{Bmatrix}
# =-\frac{wh_e}{12}\begin{Bmatrix}0\\0\\6\\h_e\\6\\-h_e \end{Bmatrix}+\begin{Bmatrix}0\\0\\Q_1^2\\Q_2^2\\Q_3^2\\Q_4^2\end{Bmatrix}
# $$
#
# Combining the equations for elements one and two give
#
# $$
# \begin{bmatrix} 
# k_{11}^1 & k_{12}^1 & k_{13}^1 & k_{14}^1 & 0 & 0\\
# k_{21}^1 & k_{22}^1 & k_{23}^1 & k_{24}^1 & 0 & 0\\
# k_{31}^1 & k_{32}^1 & k_{33}^1+k_{11}^2 & k_{34}^1+k_{12}^2  & k_{13}^2 & k_{14}^2\\
# k_{41}^1 & k_{42}^1 & k_{43}^1+k_{21}^2 & k_{44}^1+k_{22}^2  & k_{23}^2 & k_{24}^2\\
# 0 & 0 & k_{31}^2 & k_{32}^2  & k_{33}^2 & k_{34}^2\\
# 0 & 0 & k_{41}^2 & k_{42}^2  & k_{43}^2 & k_{44}^2\\
# \end{bmatrix}
# \begin{Bmatrix}
# u_1 \\ \phi_1 \\ u_2 \\ \phi_2 \\ u_3 \\ \phi_3 
# \end{Bmatrix}
# =-\frac{wh_e}{12}\begin{Bmatrix}6\\h_e\\12\\0\\6\\-h_e 
# \end{Bmatrix}+
# \begin{Bmatrix}
# Q_1^1\\Q_2^1\\Q_3^1+Q_1^2\\Q_4^1+Q_2^2\\Q_3^2\\Q_4^2
# \end{Bmatrix}
# $$
#
# We now apply the boundary conditions and partition the equations to account for the homogeneous boundary conditions:
#
# $$
# \frac{2EI}{h_e^3}\left[\begin{array}{cc|cccc} 
# 6 & 3h_e & -6 & 3h_e & 0 & 0\\
# 3h_e & 2h_e^2 & -3h_e & h_e^2 & 0 & 0  \\ \hline
# -6 & -3h_e & 12 & 0  &-6 & 3h_e\\
# 3h_e & h_e^2 & 0 & 4h_e^2  & -3h_e & h_e^2\\
# 0 & 0 & -6 & -3h_e^2  & 6 & -3h_e\\
# 0 & 0 & 3h_e & h_e^2  & -3h_e & 2h_e^2\\
# \end{array}\right]
# \begin{Bmatrix}0 \\ 0 \\\hline u_2 \\ \phi_2 \\ u_3 \\ \phi_3 
# \end{Bmatrix}
# =-\frac{wh_e}{12}\begin{Bmatrix}
# 6\\h_e\\\hline 12\\0\\6\\-h_e \end{Bmatrix}
# +\begin{Bmatrix}
# Q_1^1\\Q_2^1\\\hline 0\\0\\-F_0\\M_0
# \end{Bmatrix}
# $$
#
# We can now write the system as
#
# $$
# \begin{bmatrix} \left[ K^{11}\right] & \left[ K^{12}\right] \\ \left[ K^{21}\right] & \left[ K^{22}\right] \end{bmatrix}
#  \begin{Bmatrix} \left\{ u^{1}\right\} \\ \left\{ u^{2}\right\} \end{Bmatrix}=
#   \begin{Bmatrix} \left\{ F^{1}\right\} \\ \left\{ F^{2}\right\} \end{Bmatrix}
# $$
#
# and
#
# $$
# K^{11}u^{1} + K^{12}u^{2}=F^{1}
# $$
# $$
# K^{21}u^{1} + K^{22}u^{2}=F^{2}
# $$
#
# From which
#
# $$
# F^{1} =K^{11}u^{1} + K^{12}u^{2}
# $$
# $$
# K^{22}u^{2}=F^{2} -K^{21}u^{1}
# $$
#
# We can now solve for $u^2$ and then find the reaction forces $F^1$:
#
# $$
# u^2=\begin{Bmatrix} u_2 \\ \phi_2 \\ u_3 \\ \phi_3\end{Bmatrix}
#        =\frac{h_e^3}{2EI}\begin{bmatrix} 12 & 0  &-6 & 3h_e \\
# 			                 0 & 4h_e^2  & -3h_e & h_e^2 \\
# 			                 -6 & -3h_e^2  & 6 & -3h_e \\
# 			                 3h_e & h_e^2  & -3h_e & 2h_e^2 \\
# \end{bmatrix}^{-1}
# \begin{Bmatrix}  -wh_e \\ 0 \\ -\frac{wh_e}{2}-F_0 \\ \frac{wh_e^2}{12}+M_0\end{Bmatrix}
# $$
#
# Which gives
#
#
# $$
# \begin{Bmatrix}u_2\\\phi_2\\u_3\\\phi_3\end{Bmatrix}=
# \frac{h_e}{6EI} \begin{bmatrix}
# 3M_0h_e-4.25wh_e^3-5F_0h_e^2 \\
# 6M_0 - 7wh_e^2-9h_eF_0\\
# 12M_0h_e -12wh_e^3-16F_0h_e^2\\
# 12M_0-8wh_e^2-12h_eF_0
# \end{bmatrix}
# $$
#
# The reactions $Q_1^1$ and $Q_2^1$ can be found to be
#
# $$
# \begin{Bmatrix} Q_1^1  \\  Q_2^1 \end{Bmatrix} =
# \frac{h_e^3}{16EI}\begin{bmatrix} -6 & 3h_e & 0 & 0 \\ -3h_e & h_e^2 & 0 & 0 \end{bmatrix} \begin{Bmatrix} u_2 \\ \phi_2 \\ u_3 \\ \phi_3 \end{Bmatrix} +\frac{wh_e}{12}\begin{Bmatrix} 6 \\ h_e \end{Bmatrix} =
#  \begin{Bmatrix}
#  2wh_e+F_0\\
#  2h_e\left( wh_e+F_0\right)-M_0
#  \end{Bmatrix}
# $$



StarcoderdataJupyterScriptsDedupFiltered