If this is the wrong community for this question, please let me know.
Say I have a 4 noded quadrilateral element: Quad element in NX
(I can't add an embedded image yet...)
I know the global coordinates of each node G1 through G4.
My idea was the following:
Find the intersection of the two diagonals = D
Find the angle alpha
Find the intersection point between Axis X and edge G2-G3 = G23
Global direction of local Axis X = G23 - D
The problem is that for a quad element, the nodes are not necessarily on the same plane. Therefor, the global coordinates of the location of the element's local coordinate system is not exactly at the intersection of the two diagonals.
Assuming all the nodes are on the same plane, the intersection of the diagonals would be trivial.
A general rule in FEM is that quad elements should be as flat as possible. But how does a FEM solver usually flatten a non-flat element, or how does it define the global position of the element's local coordinate system?
It seems that I made a small mistake in the comment claiming that for a bilinear basis the resulting axes can be curved. It seems that this is not the case although it initially did not feel very intuitive to me. Let me try to explain this and how isoparametric mappings work.
Suppose we have the reference quadrilateral $D = \{(X,Y) : -1 < X < -1, -1 < Y < 1\}$. The bilinear basis functions on this domain are $$\varphi_1 = \frac14(1-X)(1-Y),\\ \varphi_2 = \frac14(1-X)(1+Y),\\ \varphi_3 = \frac14(1+X)(1+Y),\\\varphi_4 = \frac14(1+X)(1-Y).$$ The basis functions are visualized in the following figure.
The isoparametric mapping from $D$ to arbitrary quadrilateral embedded in 3D space is $$\begin{pmatrix} x\\y\\z \end{pmatrix} = F(X,Y) = \sum_{j=1}^4 \begin{pmatrix} x_j\\y_j\\z_j \end{pmatrix} \varphi_j(X,Y)$$ where $(x_j,y_j,z_j),~j\in\{1,2,3,4\},$ are the locations of the corner nodes in the global coordinate system. We can now visualize the surface $F(D)$ that results from mapping the domain $D$ through $F$ with the example corner locations $$(x_1,y_1,z_1) = (-1,0,0),\\(x_2,y_2,z_2) = (0,1,1),\\(x_3,y_3,z_3) = (1,1,0),\\(x_4,y_4,z_4) = (3/2,-1,1).$$
As you can see, the resulting quadrilateral can have a quite wonky shape. Not flat at all! This is why I assumed that if we mapped the axes through $F$ (i.e. compute the sets $F(\{(x,y) : x=0\})$ and $F(\{(x,y) : y=0)$) they could also look curved.
Then I tried doing exactly that and it was revealed that what we get are always straight lines indeed (the blue lines correspond to mapped axes):
I should have figured this out since the bilinear functions are linear exactly in the direction of the local axes.
In conclusion, I think that the best way to find these local axes is to map points through the mapping $F$.