I'm trying to compute the shadow of 3D objects on the ground. To do this I'm assuming parallel rays from Sun origin.
Let's assume the Sun direction is given by azimuth $\varphi$ and elevation $\theta$ angles. So the direction unit vector is retrieved as:
$$ \hat{\mathbf{s}} = \mathbf{R}_y\left(\theta\right)\mathbf{R}_z\left(\varphi\right)\hat{\boldsymbol{\imath}}$$
where $\mathbf{R}_{x/y/z}$ are the rotation matrices along the $x, y, z$ axes given by, for example, Euler-Rodrigues equation, and $\hat{\boldsymbol{\imath}}$ is the unit vector along $x$-axis.
A point on a ray-line has the equation: \begin{equation} \mathbf{r} = \mathbf{r}_p + \lambda\hat{\mathbf{s}} \tag{1} \end{equation}
being $\mathbf{r}_p$ a point on that line.
A point on a general plane is given by the equation:
\begin{equation} \left(\mathbf{r} - \mathbf{r_0}\right)\cdot \hat{\mathbf{n}} = 0 \tag{2} \end{equation}
being $\mathbf{r}_0$ a point on the plane and $\hat{\mathbf{n}}$ the unit vector normal to that plane.
Ensuring the validity of both equations (1) and (2) allows to compute the $\lambda$ multiplier and retrieve the position vector of the projected point (writing $\tilde{\mathbf{r}} = \mathbf{r_p} - \mathbf{r_o}$)
\begin{gather} \lambda = - \left(\frac{\tilde{\mathbf{r}}\cdot \hat{\mathbf{n}}} {\hat{\mathbf{s}}\cdot\hat{\mathbf{n}}}\right) \\ \mathbf{r}_\text{proj} = \mathbf{r}_p - \left(\frac{\tilde{\mathbf{r}}\cdot \hat{\mathbf{n}}} {\hat{\mathbf{s}}\cdot\hat{\mathbf{n}}}\right)\hat{\mathbf{s}} \tag{3} \end{gather}
Applying eq. (3) to all the vertices of a mesh $\left\{\mathbf{r}_{p_i}\right\}$ gives the shadow.
A problem can be noted: $\hat{\mathbf{s}}\cdot \hat{\mathbf{n}}$ can be zero (e.g. Sun is horizontal, i.e. the ray is ortogonal to the normal of plane).
Are there some conceptual mistakes in my reasoning? Is there a way to rewrite the operation as a linear operator applied to the set of vertices position vectors?
$$ \left[\mathbf{A}\right]\left\{\mathbf{r}_{p_i}\right\} = \mathbf{r}_\text{proj}$$
I address here your second question, but I think that it also answers a little the first one.
Consider oblique projection where unit vertical vector $\vec{k}$ has for its shadow $\left(\begin{array}{c}a\\b\end{array}\right)$ on the $xOy$ plane.
Its matrix is:
$$\left(\begin{array}{ccc}1 & 0 & a\\0 & 1 & b\\0 & 0 & 0\end{array}\right)$$
This is how I have realized the following figure representing the projection (in black) of an helix (in red):
The explanation relies plainly on the way a transformation matrix is built: its successive columns are the images of vectors $\vec{i},\vec{j},\vec{k}$ (directing $x$,$y$,$z$ axis resp.)
$\vec{i},\vec{j}$, being already on $x0y$ plane are projected on themselves.
$\vec{k}$ is projected on its shadow $\left(\begin{array}{c}a\\b\\0 \end{array}\right)=\left(\begin{array}{c}\cos(\theta)\sin(\lambda)\\ \sin(\theta)\sin(\lambda)\\ 0 \end{array}\right)$
($\theta$ being the azimuth and $\lambda$ the latitude of the sun's direction).
Edit: Another way to define the third vector using the sun's direction $\vec{s}$ is $\left(\begin{array}{c}a\\b\\0 \end{array}\right)=\left(\begin{array}{c}\left(\dfrac{\vec{s}.\vec{i}}{\vec{s}.\vec{k}}\right)\\\left(\dfrac{\vec{s}.\vec{j}}{\vec{s}.\vec{k}}\right)\\ 0 \end{array}\right).$