I'm working with static beams with Euler–Bernoulli model which ODE is
$$ \dfrac{d^2}{dx^2} \left(EI \cdot \dfrac{d^2w}{dx^2}\right) = q(x). $$
With a beam along the $x$ axis, the solution consists of finding $\vec{U} = (u_x, \ u_y, \ u_z, \ \theta_x, \ \theta_y, \ \theta_z)$ for every point on the beam, which I can transform into forces and momentums $\vec{F} = (F_x, \ F_y, \ F_z, \ M_x, \ M_y, \ M_z)$.
Question: Once I find the forces and momentums, how can I get the stress tensor $\left(\overline{\overline{\sigma}}\right)$ from the vector $\vec{F}$? I know that it must satisfy:
\begin{gather*} \left(F_x, \ F_y, \ F_z\right) = \int_{D} \overline{\overline{\sigma}} \cdot \vec{e}_x \ \ dy \ dz \\ \left(M_x, \ M_y, \ M_z\right) = \int_{D} \left(0, \ y, \ z\right) \times \left(\overline{\overline{\sigma}} \cdot \vec{e}_x\right) \ dy \ dz. \end{gather*}
I meant to get a PDE in terms of $y$ and $z$:
$$ f\left(\overline{\overline{\sigma}}, \ \dfrac{\partial \overline{\overline{\sigma}}}{\partial y}, \ \dfrac{\partial \overline{\overline{\sigma}}}{\partial z}, \ \dfrac{\partial^2\overline{\overline{\sigma}}}{\partial y^2}, \ \dfrac{\partial^2 \overline{\overline{\sigma}}}{\partial y\partial z}, \ \dfrac{\partial^2 \overline{\overline{\sigma}}}{\partial z^2}, \cdots \right) = 0. $$
But I have no idea how to find this PDE. I only have the intuition it should be an elliptic PDE because it has a unique and stable solution that minimizes the total strain energy.
The boundary is that there's no stress on the external surface:
$$ \overline{\overline{\sigma}} \cdot \vec{n} = \vec{0} \ \ \ \ \text{on} \ \partial D. $$
Motivation: I'm making software to compute the stress of each beam using finite elements. Solving the ODE, I can find $\vec{U}$ and then the effort $\vec{F} = [K] \cdot \vec{U}$, which I already have. Then I would like to compute the stress on $D$.
