\begin{equation} S_{l}=\begin{bmatrix}\frac{1}{E_{r}}&-\frac{v_{tr}}{E_{t}}&-\frac{v_{lr}}{E_{l}}&0&0&0\\ -\frac{v_{rt}}{E_{r}}&\frac{1}{E_{t}}&-\frac{v_{lt}}{E_{l}}&0&0&0\\ -\frac{v_{rl}}{E_{r}}&-\frac{v_{tl}}{E_{t}}&\frac{1}{E_{l}}&0&0&0\\ 0&0&0&\frac{1}{2G_{tl}}&0&0\\ 0&0&0&0&\frac{1}{2G_{lr}}&0\\ 0&0&0&0&0&\frac{1}{2G_{rt}}\end{bmatrix}\\ \end{equation} \begin{equation} \label{eq:stiffness_mat} \label{eq:stiffness_mat}C=S^{-1}\\ \end{equation} \begin{equation} \label{eq:trans_mat} \label{eq:trans_mat}G=\begin{bmatrix}a_{r}^{x}a_{r}^{x}&a_{r}^{y}a_{r}^{y}&a_{r}^{z}a_{r}^{z}&a_{r}^{y}a_{r}^{z}&a_{r}^{x}a_{r}^{z}&a_{r}^{x}a_{r}^{y}\\ a_{t}^{x}a_{t}^{x}&a_{t}^{y}a_{t}^{y}&a_{t}^{z}a_{t}^{z}&a_{t}^{y}a_{t}^{z}&a_{t}^{x}a_{t}^{z}&a_{t}^{x}a_{t}^{y}\\ a_{l}^{x}a_{l}^{x}&a_{l}^{y}a_{l}^{y}&a_{l}^{z}a_{l}^{z}&a_{l}^{y}a_{l}^{z}&a_{l}^{x}a_{l}^{z}&a_{l}^{x}a_{l}^{y}\\ 2a_{t}^{x}a_{l}^{x}&2a_{t}^{y}a_{l}^{y}&2a_{t}^{z}a_{l}^{z}&a_{t}^{y}a_{l}^{z}+a_{t}^{z}a_{l}^{y}&a_{t}^{x}a_{l}^{z}+a_{t}^{z}a_{l}^{x}&a_{t}^{x}a_{l}^{y}+a_{t}^{y}a_{l}^{x}\\ 2a_{r}^{x}a_{l}^{x}&2_{r}^{y}a_{l}^{y}&2a_{r}^{z}a_{l}^{z}&a_{r}^{y}a_{l}^{z}+a_{r}^{z}a_{l}^{y}&a_{r}^{x}a_{l}^{z}+a_{r}^{z}a_{l}^{x}&a_{r}^{x}a_{l}^{y}+a_{r}^{y}a_{l}^{x}\\ 2a_{r}^{x}a_{t}^{x}&2a_{r}^{y}a_{t}^{y}&2a_{r}^{z}a_{t}^{z}&a_{r}^{y}a_{t}^{z}+a_{r}^{z}a_{t}^{y}&a_{r}^{x}a_{t}^{z}+a_{r}^{z}a_{t}^{x}&a_{r}^{x}a_{t}^{y}+a_{r}^{y}a_{t}^{x}\end{bmatrix}\\ \end{equation}

Where \(a^{j}_{i}\) is the directional cosine from \(j\) to \(i\) (global to local), \(j=x,y,z\) and \(i=r,t,l\).

\begin{equation} \label{eq:cm_trans} \label{eq:cm_trans}\boldsymbol{C=G^{T}C_{l}G}\\ \end{equation} \begin{equation} \boldsymbol{G^{T}=D^{-1}}\\ \end{equation} \begin{equation} \label{eq:stress_tans} \label{eq:stress_tans}\boldsymbol{\sigma_{l}=D\sigma}\\ \end{equation}

The stiffness matrix \(\boldsymbol{C}\), stress vector \(\boldsymbol{\sigma}\) and strain vector \(\boldsymbol{\epsilon}\) are in the global coordinate system while \(\boldsymbol{C_{l}}\),\(\boldsymbol{\sigma_{l}}\),\(\boldsymbol{\epsilon_{l}}\) are in the local system.

The forces being applied in the global system (through the transformation of the elastic constants) results in stresses being calculated in the global system. As the experimental work was conducted in the local system the proportionality limit stresses are in the local system. Equation \ref{eq:stress_tans} was used in order to transform the stresses in the global system (calculated by the model) back into the local system in order to evaluate the state of failure.

The forces due to gravity from self weight of the stem and canopy were applied to appropriate domains. The force due to gravity from self weight of the stem was applied as a body force to the whole domain \(\omega\) (i.e. the whole stem) in Equation \ref{eq:Bs}. Green density \(\rho\) is calculated at each point using Equation \ref{eq:constants} and g is gravity (\(9.81m/s\)). The canopy is assumed to take the geometry of the upper half of an ellipsoid described in Equation \ref{eq:Vc}, where \(z_{0}\) is the height at the bottom of the current element and \(z_{1}\) is the height at the top of the current element, with the top of the ellipsoid at the top of the stem, height \(h\). The start height \(S_{c}\) of the canopy was estimated using Equation \ref{eq:Sc} which was derived from data presented by \citet{waghorn_influence_2007}. In order to estimate crown radius, \(r_{c}\) , Equation \ref{eq:rc} is used, which was derived from the assumption that the maximum radius which can be achieved by the crown is half the distance between two trees evenly spaced in the stand. The weight of the canopy on the stem is applied as a body force described in Equation \ref{eq:Bc1}.