Your boundary conditions seem to be not quite correct according to the mechanical problem. Sorry, I don't have the time to go through your code today, but I got a version running, although this will take some time and might be an overkill, since it is based on the full 3D theory. I have to go home now, I will try to take a look at your code again tomorrow, if nobody else finds the error.

I treated the problem as follows (length in `x` direction `L1`, in virtual `y` direction `L2` and in `z` direction `L3`)

[![enter image description here][1]][1]

First, let's get a reference solution of the 1D theory:

    (*Geometry - in m*)
    L1 = 2;
    L2 = 0.1;
    L3 = 0.2;
    Iy = L2*L3^3/12;
    (*Force and densities - in N*)
    F = 10;
    qA = F/(L1*L2); (*area density - for 3D FEM*)
    ql = qA*L2; (*line density - for 1D theory*)
    (*Material parameters*)
    Em = 2.1*10^9; (*Young's modulus*)
    nu = 0.3;(*Poisson's ration*)
    (*1D theory*)
    wsol1D = DSolveValue[{
        Em*Iy*D[w[x], {x, 4}] == ql
        , (w[0]) == 0, (w[L1]) == 0
        , (w''[0]) == 0, (w''[L1]) == 0
        }, w, x];
    My = -Em*Iy*wsol1D''[x];(*Moment*)
    sig = My/Iy*z;(*normal stress*)
    GraphicsRow[{
      Plot[wsol1D[x], {x, 0, L1}, AxesLabel -> {"x", "w(x)"}]
      , Plot[sig /. x -> L1/2, {z, -L3/2, L3/2}, 
       AxesLabel -> {"z", "\[Sigma](x=L1/2,z)"}]
      }
     , ImageSize -> Large
     ]

[![enter image description here][2]][2]

Now, let's get the full 3D FEM solution (takes 1.4 seconds for me) with a area force density

    (*FEM solution*)
    Needs["NDSolve`FEM`"]
    (******************************)
    (*Region definition*)
    reg = Cuboid[{0, -L2/2, -L3/2}, {L1, L2/2, L3/2}];
    (******************************)
    (*Isotropic material stiffness - fourth-order tensor*)
    (*Identities*)
    I2 = IdentityMatrix@3;
    IdI = TensorProduct[I2, I2];
    I4 = TensorTranspose[IdI, {1, 3, 2, 4}];
    IS = (I4 + TensorTranspose[I4, {1, 2, 4, 3}])/2;
    (*Isotropic projectors*)
    P1 = 1/3*IdI;
    P2 = IS - P1;
    (*Isotropic stiffness*)
    Ciso = l1*P1 + l2*P2;
    l1 = 3*Km;
    l2 = 2*Gm;
    Km = 1/3*Em/(1 - 2*nu);
    Gm = 1/2*Em/(1 + nu);
    (******************************)
    (*Equations*)
    eq = Table[
       Inactive[Div][
         Ciso[[i, ;; , 1, ;;]].Inactive[Grad][u[x, y, z], {x, y, z}], {x, 
          y, z}]
        + Inactive[Div][
         Ciso[[i, ;; , 2, ;;]].Inactive[Grad][v[x, y, z], {x, y, z}], {x, 
          y, z}]
        + Inactive[Div][
         Ciso[[i, ;; , 3, ;;]].Inactive[Grad][w[x, y, z], {x, y, z}], {x, 
          y, z}]
       , {i, 3}
       ];
    (******************************)
    (*BCs*)
    (*Dirichlet*)
    bcD = {
       DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0, 
         w[x, y, z] == 0}, x == 0 && z == 0]
       , DirichletCondition[{v[x, y, z] == 0, w[x, y, z] == 0}, 
        x == L1 && z == 0]
       };
    (*Neumann*)
    bcN = {0, 0, NeumannValue[-qA, z == -L3/2]};
    (******************************)
    (*Solution*)
    {usol, vsol, wsol} = 
       NDSolveValue[{eq == bcN, bcD}, {u, v, w}, Element[{x, y, z}, reg], 
        Method -> {"PDEDiscretization" -> {"FiniteElement", 
            "MeshOptions" -> {"MaxCellMeasure" -> 0.0001, 
              "MeshOrder" -> 2}}}
        ]; // AbsoluteTiming

> {1.4435, Null}

You can take a look at the deformed mesh if you want

    mesh = usol["ElementMesh"];
    Show[{
      mesh["Wireframe"]
      , ElementMeshDeformation[mesh, {usol, vsol, wsol}, 
        "ScalingFactor" -> 10^4][
       "Wireframe"[
        "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]
      }, Axes -> True, AxesLabel -> {x, y, z}]

[![enter image description here][3]][3]

The FEM solution (FEM solution as red points) is in good accordance with the analytical 1D theory

    Plot[wsol1D[x], {x, 0, L1}, 
     Epilog -> {PointSize -> Medium, Red, 
       Point[Table[{x, wsol[x, 0, 0]}, {x, 0, L1, L1/10}]]}]

[![enter image description here][4]][4]

You can get the stress distribution at any point with the full 3D Hooke's law $\sigma_{ij} = C_{ijkl} u_{k,l}$ (FEM solution as red points)

    uv[x_, y_, z_] := {usol[x, y, z], vsol[x, y, z], wsol[x, y, z]}
    eps[xs_, ys_, zs_] := 
     D[uv[x, y, z], {{x, y, z}, 1}] /. {x -> xs, y -> ys, z -> zs}
    (*linear map of second order tensor B over fourth-order tensor A*)
    lm[A_, B_] := TensorContract[TensorProduct[A, B], {{3, 5}, {4, 6}}]
    (*Get Cauchy stress at x==L1/2 depending on z with 3D Hooke's law*)
    sigloc = lm[Ciso, eps[L1/2, 0, z]][[1, 1]];
    siglocdata = Table[{zi, sigloc /. z -> zi}, {zi, -L3/2, L3/2, L3/10}];
    Plot[sig /. x -> L1/2, {z, -L3/2, L3/2}, 
     AxesLabel -> {"z", "\[Sigma](x=L1/2,z)"}, 
     Epilog -> {PointSize -> Medium, Red, Point@siglocdata}]

[![enter image description here][5]][5]

 


  [1]: https://i.sstatic.net/doRfJ.png
  [2]: https://i.sstatic.net/eKWuI.png
  [3]: https://i.sstatic.net/mxGeG.png
  [4]: https://i.sstatic.net/LEmvI.png
  [5]: https://i.sstatic.net/E9ugP.png