(*finite element using Ritz method for axial loaded beam
by Nasser M. Abbasi
version 6/2/2012*)
Using the Ritz method,the solution to axial deformation is found.As more elements are used,the solution is \
seen to converge to the exact solution.*)

Manipulate[process[totalLength, numberOfElements, area , youngModulus, traction, endLoad],
 Grid[
  {
   
   {
    Control[{{totalLength, 1, Column[{Style["Total", 10], Style["Length", 10]}]}, .1, 10, .1, Appearance -> "Labeled", 
      ImageSize -> Tiny }],
    Control[{{numberOfElements, 3, Column[{Style["number", 10], Style["of elements", 10]}]}, 2, 14, 1, Appearance -> "Labeled", 
      ImageSize -> Tiny}],
    Control[{{area , 1, Column[{Style["A, cross", 10], Style["section", 10], Style["area", 10]}]}, 1, 10, 1, 
      Appearance -> "Labeled", ImageSize -> Tiny}]
    },
   
   {
    Control[{{youngModulus , 1, Column[{Style["E, Young's", 10], Style["modulus", 10]}]}, .01, 5, .01, Appearance -> "Labeled", 
      ImageSize -> Tiny}],
    Control[{{traction , 1, Column[{Style["c, traction", 10]}]}, .01, 5, .01, Appearance -> "Labeled", ImageSize -> Tiny}],
    Control[{{endLoad , 0, Column[{Style["P, force at", 10], Style["end", 10]}]}, 0, 10, .01, Appearance -> "Labeled", 
      ImageSize -> Tiny}]
    }
   },
  Frame -> All, FrameStyle -> Directive[AbsoluteThickness[.1], Gray], Spacings -> {2, 0},
  ItemSize -> Automatic
  ],
 
 {{plotWidth, 280}, ControlType -> None},  
 {{plotHeight, 200}, ControlType -> None},
 {{plotImagePadding, 38}, ControlType -> None},
 (*ContentSize\[Rule]Automatic,*)
 FrameMargins -> 0,
 ImageMargins -> 0,
 ContinuousAction -> False,
 SynchronousUpdating -> True,
 ControlPlacement -> Top,
 FrameMargins -> 0,
 AutorunSequencing -> {{2, 40}},
 ImageMargins -> 0
 , Initialization :>
  (
   process[totalLength_, numberOfElements_, area_ , youngModulus_, traction_, endLoad_] := 
     Module[{numberOfNodes, \[CapitalNu], L, \[CapitalEpsilon], c, \[CapitalAlpha], p, y, z, \[CapitalBeta], \[CapitalPi], 
       strainEnergy, externalLoad, pVector, gamma, d, x, eqs, b, kmat, output, sol, data, p2, p1, p3, p4, exactSolution, strain, 
       strainData, stressData, actualStrain, actualStress, k, pVectorString},
      numberOfNodes = numberOfElements + 1;
      
      (\[CapitalNu] = {{(L - s)/L, s/L}}) // MatrixForm;
      (\[CapitalBeta] = D[\[CapitalNu], s]) // MatrixForm;
      \[CapitalPi] = 0;
      strainEnergy = Table[0, {numberOfElements}];
      externalLoad = Table[0, {numberOfElements}];
      pVector = Table[0, {numberOfNodes}];
      
      gamma = Table[0, {numberOfElements}];
      Do[
       d = Array[{Subscript[u, #]} &, 2, k];
       strainEnergy[[k]] = (Simplify[(\[CapitalAlpha] \[CapitalEpsilon])/2 \!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(L\)]\(\((Transpose[d] . Transpose[\[CapitalBeta]] . \[CapitalBeta] . 
                d)\) \[DifferentialD]s\)\)])[[1, 1]];
       x = (k - 1)*L;
       gamma[[k]] = Simplify[c Transpose[d].\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(L\)]\(\((Transpose[\[CapitalNu]] \((x + s)\))\)\ \[DifferentialD]s\)\)];
       
       If[k == numberOfElements, externalLoad[[k]] = d[[2]]*p, externalLoad[[k]] = 0];
       
       \[CapitalPi] = \[CapitalPi] + (strainEnergy[[k]] - gamma[[k]] - externalLoad[[k]])[[1, 1]];
       , {k, 1, numberOfElements}
       ];
      \[CapitalPi] = Simplify[\[CapitalPi]];
      
      (*Print["Strain energy initially=",strainEnergy];*)
      
      d = Array[Subscript[u, #] &, numberOfNodes, 1];
      eqs = Table[0, {numberOfNodes}];
      Do[
       eqs[[k]] = D[\[CapitalPi], d[[k]]] == 0;
       , {k, 1, numberOfNodes}
       ];
      
      {b, kmat} = Normal[CoefficientArrays[eqs, d]];
      
      pVector[[-1]] = p;
      pVectorString = pVector;
      pVectorString[[-1]] = "p";
      output = Text[
        ToString[Style["\!\(\*FractionBox[\(\[CapitalEpsilon]\\\ \[CapitalAlpha]\), \(L\)]\)", 14], FormatType -> TraditionalForm] <>

         
         ToString[MatrixForm[kmat/(\[CapitalEpsilon] \[CapitalAlpha]/L)], FormatType -> TraditionalForm] <> 
         ToString[MatrixForm[Transpose[{d}]], FormatType -> TraditionalForm] <> " = " <> 
         ToString[Style["c \!\(\*SuperscriptBox[\(L\), \(2\)]\) ", 14], FormatType -> TraditionalForm] <>
         ToString[MatrixForm[(-b - pVector)/(c L^2)], FormatType -> TraditionalForm] <> " + " <> 
         ToString[MatrixForm[pVectorString], FormatType -> TraditionalForm]];
      pVector[[-1]] = p;
      sol = LinearSolve[kmat[[2 ;; -1, 2 ;; -1]], -b[[2 ;; -1]]];
      PrependTo[sol, 0];
      
      
      L = totalLength/numberOfElements;
      \[CapitalEpsilon] = youngModulus;
      c = traction;
      \[CapitalAlpha] = area;
      p = endLoad;
      (*Print["sol=",sol];*)
      data = Table[{(i - 1)*L, sol[[i]]}, {i, 1, numberOfNodes}];
      
      p1 = ListPlot[data, Joined -> True, PlotMarkers -> {Automatic, Medium}, PlotStyle -> {Red, PointSize[Large]}, 
        ImageSize -> {plotWidth, plotHeight}, ImagePadding -> plotImagePadding,
        AspectRatio -> 0.48];
      
      (*Clear[z,y];*)
      exactSolution = 
       First@DSolve[{\[CapitalAlpha] \[CapitalEpsilon] y''[z] == -c z, y[0] == 0, 
          y'[totalLength] == p/(\[CapitalAlpha] \[CapitalEpsilon])}, y[z], z];
      y = y[z] /. exactSolution;
      (*Print["exact solution=",y];*)
      
      p2 = Plot[y, {z, 0, totalLength}, ImageSize -> {plotWidth, plotHeight}, ImagePadding -> plotImagePadding, AspectRatio -> 0.48,
         PlotRange -> All];
      
      p3 = Show[{p2, p1}, Frame -> True, FrameLabel -> {{"U(x)", None}, {"x", "Actual vs. FEM displacement"}}];
      
      Do[d[[i]] = {sol[[i]]}, {i, 1, numberOfNodes}];
      
      (*Print["strain energy before=",N[strainEnergy]];*)
      strain = Table[0, {numberOfElements}];
      
      For[k = 1, k <= numberOfElements, k++,
       strain[[k]] = (\[CapitalBeta].{{sol[[k]]}, {sol[[k + 1]]}})[[1, 1]]
       ];
      
      (*Print["strain =",N[strain]];*)
      
      strainData = Table[{{L (k - 1), strain[[k]]}, {L k, strain[[k]]}}, {k, 1, numberOfElements}];
      
      (*Print["straindata=",N[straindata]];*)
      
      stressData = 
       Table[{{L (k - 1), \[CapitalEpsilon] strain[[k]]}, {L k, \[CapitalEpsilon] strain[[k]]}}, {k, 1, numberOfElements}];
      
      (*Print["stress data=",N[stressData]];*)
      
      p1 = ListPlot[stressData, Joined -> True, AxesOrigin -> {0, 0}, PlotStyle -> {Red}, ImageSize -> {plotWidth, plotHeight}, 
        ImagePadding -> plotImagePadding, AspectRatio -> 0.48];
      
      actualStrain = D[y, z];
      
      actualStress = \[CapitalEpsilon] actualStrain;
      p2 = Plot[actualStress, {z, 0, totalLength}, ImageSize -> {plotWidth, plotHeight}, ImagePadding -> plotImagePadding, 
        AspectRatio -> 0.48];
      
      p4 = Show[{p1, p2}, Frame -> True, PlotRange -> All, FrameLabel -> {{"\[Sigma](x)", None}, {"x", "Actual vs. FEM stress"}}];
      
      Grid[{{output, SpanFromLeft}, {p3 , p4}}, Frame -> None, Alignment -> Center, Spacings -> {1, 0}, 
       ItemSize -> Automatic(*{Scaled[.5],Scaled[.5]}*)]
      ];
   
   )
 ]