## 3 Mixed FE-GDQ Numerical Formulation
\tag{26}\]For a unified treatment of the problem, the degrees of freedom (DOFs) can be divided into the domain- and boundary-type DOFs, as follows:  
\[\mathbf{U}_{d}=\left\{ \begin{array}{l}U_{12}\\ U_{13}\\ \vdots\\ U_{N(N_{2}-1)}\end{array} \right\},\ \ \mathbf{V}_{d}=\left\{ \begin{array}{l}V_{12}\\ V_{13}\\ \vdots\\ V_{N(N_{2}-1)}\end{array} \right\},\ \ \mathbf{W}_{d}=\left\{ \begin{array}{l}W_{12}\\ W_{13}\\ \vdots\\ W_{N(N_{2}-1)}\end{array} \right\}, \tag{27}\]  
\[\mathbf{U}_{b}=\left\{ \begin{array}{l}U_{11}\\ U_{21}\\ \vdots\\ U_{N{N_{2}}}\end{array} \right\},\ \ \ \ \mathbf{V}_{b}=\left\{ \begin{array}{l}V_{11}\\ V_{21}\\ \vdots\\ V_{N{N_{2}}}\end{array} \right\}^{\prime},\ \ \ \ \mathbf{W}_{b}=\left\{ \begin{array}{l}W_{11}\\ W_{21}\\ \vdots\\ W_{N{N_{2}}}\end{array} \right\}^{\prime},\]  
In which \(W_{mn}=W_{m}(z_{n},t),V_{mn}=V_{m}(z_{n},t),U_{mn}=U_{m}(z_{n},t)\). This means that Equations (21)-(23) can be rearranged in matrix form as:  
\[\left[ \begin{array}{cc}\mathbf{k}_{bb}&\mathbf{k}_{bd}\\ \mathbf{k}_{db}&\mathbf{k}_{dd}\end{array} \right]\left[ \begin{array}{cc} \mathbf{d}_{b}\\ \mathbf{d}_{d}\end{array} \right]=P\left[ \begin{array}{cc}\mathbf{0}&\mathbf{ 0}\\ \mathbf{0}&\mathbf{G}\end{array} \right]\left[ \begin{array}{cc}\mathbf{d}_{b} \\ \mathbf{d}_{d}\end{array} \right] \tag{28}\]  
where the stiffness quantities \(\mathbf{k}_{bb}\), \(\mathbf{k}_{db}\), \(\mathbf{k}_{bd}\), \(\mathbf{k}_{dd}\) refer to the boundary, \(b\), and domain, \(d\), weighting coefficients of the plate respectively, \([\mathbf{d}_{b}\ \mathbf{d}_{d}]^{T}\) is the displacement vector, \(P\) refers to the buckling load and \(\mathbf{G}\) is the stability matrix due to the in-plane stresses. At the same time, from Equations (23)-(26), the boundary weight coefficients can be replaced by the domain weight coefficients as follows:  
\[\mathbf{d}_{b}=\mathbf{k}_{bb}{}^{-1}\mathbf{k}_{bd}\mathbf{d}_{d} \tag{29}\]  
By substitution of \(\mathbf{d}_{b}\) from Equation (29) into Equation (28), and considering the harmonic solution \(\left[ \begin{array}{c}\mathbf{U}_{d}\\ \mathbf{V}_{d}\\ \mathbf{W}_{d}\end{array} \right]=\left[ \begin{array}{c}\overline{\mathbf{U} }_{d}\\ \overline{\mathbf{W}}_{d}\\ \overline{\mathbf{W}}_{d}\end{array} \right]e^{iwt}\), the governing equations of the problem can be rewritten in terms of the domain unknowns, as follows:  
\[\mathbf{K}\left[ \begin{array}{c}\overline{\mathbf{U}}_{d}\\ \overline{\mathbf{V}}_{d}\\ \overline{\mathbf{W}}_{d}\end{array} \right]-P\mathbf{G}\left[ \begin{array}[] {c}\overline{\mathbf{U}}_{d}\\ \overline{\mathbf{W}}_{d}\end{array} \right]=0 \tag{30}\]  
where \(\mathbf{K}\) is the stiffness matrix, defined as:  
\[\mathbf{K}=\mathbf{K}_{dd}-\mathbf{K}_{db}\mathbf{K}_{bb}{}^{-1}\mathbf{K}_{bd} \tag{31}\]  
The solution of Equation (30) corresponds to the critical buckling load of the structure under in-plane conditions, labeled hereafter as \(\lambda\). In what follows, three different BCs for the buckling analysis of the plate structure are considered:  
(i) Simply supported BCs at all edges (SSSS):  
\[ \begin{array}{l}w(0,y,z)=w(a,y,z)=w(x,0,z)=w(x,b,z)=0\\ u(0,y,z)=u(a,y,z)=v(x,0,z)=v(x,b,z)=0\end{array} \] (32) (ii)  
Clamped BCs at edges parallel to the \(y\)-axis (i.e., at \(x=0,\ a\)) and free BCs at edges parallel to the \(x\)-axis (i.e., at \(y=0,\ b\)) (CFCF):  
\[u,v,w(0,y,z)=u,v,w(a,y,z)=0 \tag{33}\]  
decrease or increase, under the same assumptions of porosity coefficient and distribution, depending on the selected BC. Differently from the FCFC case, a clear reduction of the buckling load is noticed for a SSSS and CFCF porous plate with an increased aspect ratio. This confirms the strict dependence of the stiffness on the geometrical dimensions and BCs of the structural member.