MechanicsExplorers`Beam`

A User-friendly Exploration Program for the Bending of Beams

Stephan Kaufmann
Institute of Mechanics
Swiss Federal Institute of Technology, ETH Zürich
CH-8092 Zürich, Switzerland
Email: kaufmann@ifm.mavt.ethz.ch
URL: http://www.ifm.ethz.ch/~kaufmann/

April 10, 1997

© 1997 by Stephan Kaufmann

Source Code

Package Beginnings

Package Preamble

\(\( (* : Version : \ Mathematica\ 3.0\ \ *) 


     (* : Name : \ MechanicsExplorers`Beam`\ *) 


     (* : Title : \ Beam\ *) 

 (* : Author : \ Stephan\ Kaufmann\ *) 


     (*\( : Date : \ April\ 10\), \ 1997\ *) 


     (* : Package\ \(Version : \ 1.0\)\ *) 


     (* : Mathematica\ \(Version : \ 3.0\)\ *) 


     (* : Copyright : \ ©\ 1997\ by\ Stephan\ Kaufmann\ *) 


     (*\( : Keywords : \ beam\), \ bending, \ mechanics\ *) 


     (* : Requirements : \ Mathematica\ 3.0\ or\ higher\ *) 


     (* : Warnings : 


          The\ code\ in\ this\ file\ is\ in\ no\ way\ meant\ to\ be\ human\ 
            readable.\ Please\ see\ the\ source\ file\ BeamExplorer.nb.
*) 

    
 (* : Source : \ \(None.\)\ *) 


     (* : Limitations : \ \(None.\)\ *) 


     (*\( : Discussion : 


        The\ new\ features\ of\ Mathematica\ Version\ 3.0\ make\ it\ easy\ to
            \ program\ user - 
          friendly\ interfaces\ to\ Mathematica\ packages.\ Palettes\ make\ it
            \ possible\ to\ hide\ or\ automate\ the\ interaction\ with\ the\ 
            Mathematica\ language.\ Such\ programs\ open\ new\ possibilities\ 
            in\ education.

This\ package\ allows\ for\ the\ visualization\ 
            of\ the\ bending\ of\ a\ beam.\ Beams\ are\ ' constructed'\ by\ 
            defining\ supports\ \((simple, \ fixed)\)\), \ hinges, \ 
      and\ loads\ 
        \((discrete\ or\ distributed\ forces, \ moments)\).
		
The\ 
        shearing\ force, \ the\ bending\ moment, \ the\ deflection, \ 
      and\ its\ slope, \ 
      as\ well\ as\ the\ resulting\ forces\ and\ moments\ in\ the\ supports\ 
        are\ calculated\ symbolically.\ After\ the\ calculation, \ 
\ \ 
      the\ results\ can\ be\ visualized\ graphically\ or\ displayed\ as\ 
        formulas.

There\ are\ three\ different\ levels\ of\ simplicity\ for
        \ the\ interaction\ of\ the\ user

      \ \ \ \ \ 
      with\ the\ 
        package.


          \(-\ \(Simple : \ 
              Everything\ is\ evaluated\ automatically\ after\ clicking\ at\ 
                the\ corresponding\ buttons\ of\ a\ palette.\ Numerical\ input
                \ is\ typed\ into\ Input[]\ dialogs\ or\ into\ 
                templates.
\ \ 

                  \(-\ \(Medium : \ 
                      The\ buttons\ paste\ the\ corresponding\ Mathematica\ 
                        expression\ into\ a\ notebook.\ The\ user\ fills\ out
                        \ the\ \ placeholders\ and\ evaluates\ the\ 
                        expressions.
\ \ 

                          \(-\ \(Advanced : \ 
                              Experienced\ users\ of\ Mathematica\ will\ want
                                \ to\ use\ the\ versions\ of\ the\ functions\ 
                                which\ operate\ in\ the\ traditional\ 
                                \(way.\)\)\)\)\)\)\)
*) 



     (* : Code\ \(Description : \ See\ the\ source\ file\ Beam.nb.\)	*) 
      \)\)

Package Begin

BeginPackage["MechanicsExplorers`Beam`", "Calculus`DiracDelta`", 
      "Graphics`Arrow`"];

Turn off Spelling Checker

Off[General::"spell1"]

Usage Statements

Beam::"usage" = 
      "Beam[{leftEnd,rightEnd,bendingStiffnessFunction},{supportsAndLoads}] \
represents a beam with the given left and right end points, the bending \
stiffness (given as a pure function) and the list of supports (Simple[...], \
Fixed[...], and Hinge[...] objects) and loads (Force[...] and Moment[...] \
objects).";

Simple::"usage" = 
      "Simple[position] represents a simple support (no deflection) at the \
given position in a Beam[...] object.";

Fixed::"usage" = 
      "Fixed[position] represents a fixed support (no deflection and \
horizontal slope) at the given position in a Beam[...] object.";

Hinge::"usage" = 
      "Hinge[position] represents a hinge (no bending moment) at the given \
position in a Beam[...] object.";

Force::"usage" = 
      "Force[position,value] represents a discrete force at the given \
position with the given value in a Beam[...] object. \
Force[{left,right},function] defines a distributed force in the interval \
[left,right], given by a pure function. The positive force direction is \
downwards.";

Moment::"usage" = 
      "Moment[position,value] represents a moment at the given position with \
the given value in a Beam[...] object. The positive moment direction points \
into the screen.";

NewWorkBeam::"usage" = 
      "NewWorkBeam[] initializes the work beam with end points {0,1}, \
constant bending stiffness, and no supports or loads. If necessary, a new \
work notebook is opened.";

NewBeam::"usage" = 
      "NewBeam[] returns a default Beam object (with end points {0,1}, \
constant bending stiffness, and no supports or loads).\nNewBeam[beam] defines \
the variable beam as a default Beam object.";

DefineEndPoints::"usage" = 
      "DefineEndPoints[beam,{left,right}] defines the left and right end \
points of the given beam.\nDefineEndPoints[{left,right}] acts on the work \
beam.\nDefineEndPoints[] asks for input and acts on the work beam.";

DefineBendingStiffness::"usage" = 
      "DefineBendingStiffness[beam,function] defines the bending stiffness of \
the given beam. The stiffnes function should either be given as a pure \
function (e.g., (1+#)&), or as a list containing an expression and the \
variable (e.g., {1+x,x}). Discontinuities of the bending stiffness should be \
expressed in term of (linear combinations of) UnitStep functions (e.g., 1+3 \
UnitStep[#-.5]&).\nDefineBendingStiffness[position] acts on the work beam.\n\
DefineBendingStiffness[] asks for input and acts on the work beam.";

AddSimpleSupport::"usage" = 
      "AddSimpleSupport[beam,position] adds a simple support at the given \
position to the given beam.\nAddSimpleSupport[position] acts on the work \
beam.\nAddSimpleSupport[] asks for input and acts on the work beam.";

AddFixedSupport::"usage" = 
      "AddFixedSupport[beam,position] adds a fixed support at the given \
position to the given beam.\nAddFixedSupport[position] acts on the work beam.\
\nAddFixedSupport[] asks for input and acts on the work beam.";

AddHinge::"usage" = 
      "AddHinge[beam,position] adds a hinge at the given position to the \
given beam.\nAddHinge[position] acts on the work beam.\nAddHinge[] asks for \
input and acts on the work beam.";

AddForce::"usage" = 
      "AddForce[beam,position,value] adds a discrete force of the given value \
at the given position to the given beam. AddForce[beam,{left,right},function] \
adds a distributed force, starting at position 'left' and ending at position \
'right'. The function must either be a pure function (e.g., #^2&) or a list \
containing an expression and a variable (e.g., {x,x^2}. Positive forces are \
in direction of the y-axis, i.e., vertically downwards.\n\
AddForce[position,value] and AddForce[{left,right},function]act on the work \
beam.\nAddForce[] asks for input and acts on the work beam.";

AddDistributedForce::"usage" = 
      "AddDistributedForce[] asks for input and defines a distributed force \
for the work beam. It is similar to AddForce[].";

AddMoment::"usage" = 
      "AddMoment[beam,position,value] adds a discrete moment of the given \
value at the given position to the given beam. Positive moments are in \
direction of the z-axis, i.e., pointing into the screen.\n\
AddMoment[position,value] acts on the work beam.\nAddMoment[] asks for input \
and acts on the work beam.";

RemoveLastElement::"usage" = 
      "RemoveLastElementElement[] removes the support or load which has been \
defined last. Multiple removals are possible.";

RemoveElement::"usage" = 
      "RemoveElement[beam,element] removes a support or load element from the \
beam. Patterns can be used in the specification of the element(s).";

SolveBeam::"usage" = 
      "SolveBeam[beam] returns a list containing a list of 4 pure functions \
(the shearing force, the bending moment, the slope of the deflection, and the \
deflection), and a Beam object, containing the elements of the unsolved beam \
and the forces and moments at the supports. The positive direction of the \
deflection and the slope is in y-direction, i.e., downwards.\nSolveBeam[] \
solves the work beam and assigns the solution to the work solution.";

PlotShearingForce::"usage" = 
      "PlotShearingForce[result] takes a result of SolveBeam and plots the \
shearing force.\nPlotShearingForce[beam] solves the beam and plots the \
shearing force.\nPlotShearingForce[] calculates the work solution, if \
necessary, and plots the shearing force.";

PlotBendingMoment::"usage" = 
      "PlotBendingMoment[result] takes the result of SolveBeam and plots the \
bending moment.\nPlotBendingMoment[beam] solves the beam and plots the \
bending moment.\nPlotBendingMoment[] calculates the work solution, if \
necessary, and plots the bending moment.";

PlotSlope::"usage" = 
      "PlotShearingForce[result] takes the result of SolveBeam and plots \
the slope of the deflection. Positive values are in the y-direction, i.e., \
downwards.\nPlotSlope[beam] solves the beam and plots the shearing force.\n\
PlotSlope[] calculates the work solution, if necessary, and plots the \
slope.";

PlotDeflection::"usage" = 
      "PlotDeflection[result] takes the result of SolveBeam and plots the \
the deflection. Positive values are in the y-direction, i.e., downwards.\n\
PlotDeflection[beam] solves the beam and plots the deflection.\n\
PlotDeflection[] calculates the work solution, if necessary, and plots the \
deflection.";

PlotBendingStiffness::"usage" = 
      "PlotBendingStiffness[beam] plots the bending stiffness.\n\
PlotBendingStiffness[] plots the bending stiffness of the work beam.\n\
PlotBendingStiffness[result] plots the bending stiffness of a result of \
SolveBeam.";

PrintShearingForce::"usage" = 
      "PrintShearingForce[solution] prints the shearing force of the solution \
as an expression in the variable 'x'. To avoid conflicts, 'x' is a string. \
Therefore, the printed result cannot be used for further calculations.\n\
PrintShearingForce[] prints the shearing force of the work solution.";

ShearingForce::"usage" = 
      "ShearingForce[solution] returns the shearing force of the solution as \
a pure function suitable for further calculations.\nShearingForce[] returns \
the shearing force of the work solution.";

PrintBendingMoment::"usage" = 
      "PrintBendingMoment[solution] prints the bending moment of the solution \
as an expression in the variable 'x'. To avoid conflicts, 'x' is a string. \
Therefore, the printed result cannot be used for further calculations.\n\
PrintBendingMoment[] prints the bending moment of the work solution.";

BendingMoment::"usage" = 
      "BendingMoment[solution] returns the bending moment of the solution as \
a pure function suitable for further calculations.\nBendingMoment[] returns \
the bending moment of the work solution.";

PrintSlope::"usage" = 
      "PrintSlope[solution] prints the slope of the solution as an expression \
in the variable 'x'. The positive slope direction is downwards. To avoid \
conflicts, 'x' is a string. Therefore, the printed result cannot be used for \
further calculations.\nPrintSlope[] prints the slope of the work solution.";

Slope::"usage" = 
      "Slope[solution] returns the slope of the solution as a pure function \
suitable for further calculations. The positive slope direction is downwards.\
\nSlope[] returns the slope of the work solution.";

PrintDeflection::"usage" = 
      "PrintDeflection[solution] prints the deflection function of the \
solution as an expression in the variable 'x'. The positive deflection \
direction is downwards. To avoid conflicts, 'x' is a string. Therefore, the \
printed result cannot be used for further calculations.\nPrintDeflection[] \
prints the deflection of the work solution.";

Deflection::"usage" = 
      "Deflection[solution] returns the bending deflection of the solution \
as a pure function suitable for further calculations. The positive deflection \
direction is downwards.\nDeflection[] returns the deflection of the work \
solution.";

PrintBendingStiffness::"usage" = 
      "PrintBendingStiffness[beam] prints the bending stiffness function of a \
beam (or a solution) as an expression in the variable 'x'. To avoid \
conflicts, 'x' is a string. Therefore, the printed result cannot be used for \
further calculations.\nPrintBendingStiffness[] prints the bending stiffness \
of the work beam.";

BendingStiffness::"usage" = 
      "BendingStiffness[bem] returns the bending stiffness of the beam as a \
pure function suitable for further calculations.\nBendingStiffness[] returns \
the bending stiffness of the work beam.";

DrawBeam::"usage" = 
      "DrawBeam[beam] draws a sketch of the beam.DrawBeam[] draws a sketch of \
the work beam.";

DrawSolution::"usage" = 
      "DrawSolution[solution] draws a sketch of the solution.\nDrawSolution[] \
draws a sketch of the work solution.\nDrawSolution[beam] solved the beam and \
draws the solution. With this version, the solution is lost and has to be \
recalculated for further drawings.";

DrawAll::"usage" = 
      "DrawAlll[beam,solution] draws the original beam and the solved beam \
and plots the bending stiffness, the shearing force, the bending moment, the \
slope, and the deflection.\nDrawAll[] uses the work beam and the work \
solution.\nDrawAll[beam] solved the beam and draws everything (the solution \
is lost).";

SetInputDialogs::"usage" = 
      "SetInputDialogs[value] expects the value True or False. If set True, \
numerical input is given by Input[] dialogs. If set False, numerical input is \
given by inserting values into templates and evaluating the expressions.";

SetDrawEach::"usage" = 
      "SetDrawEach[value] expects the value True or False. If set True, all \
changes of the work beam are sketched automatically.";

SetWorkBeam::"usage" = 
      "SetWorkBeam[beam] sets the work beam to the given beam.";

GetWorkBeam::"usage" = "GetWorkBeam[] returns the work beam.";

SetWorkSolution::"usage" = 
      "SetWorkSolution[solution] sets the work solution to the given \
solution.";

GetWorkSolution::"usage" = 
      "GetWorkSolution[] returns the work beam.";

BeamPasteEvaluate::"usage" = 
      "BeamPasteEvaluate[functions] pastes the function into the work \
notebook and evaluates it.";

Messages

Beam::"undef" = 
      "Your beam seems to be kinematically undefinded. Consider adding more \
supports.";

Beam::"nosol" = "Your beam could not be solved.";

Beam::"badends" = 
      "Warning: Your beam contains elements outside of its end points. The \
end points have been automatically adjusted to contain all elements.";

Beam::"nobeam" = "This function only works for Beam objects.";

Beam::"multelem" = 
      "The beam contained multiple copies of elements. Only one of each is \
kept.";

Beam::"illdef" = "Your beam is ill-defined.";

Beam::"symbends" = 
      "Warning: Your beam contains symbolic end points. All supports and \
loads are assumed to be between the end points. You will not be able to plot \
the results without defining numerical end points.";

Beam::"numends" = 
      "The draw and plot functions can only work for numerical end points.";

Beam::"fixedhinge" = 
      "A fixed support and a hinge at the same position are equivalent to a \
simple support and a hinge. They have been replaced accordingly.";

Beam::"simplefixed" = 
      "It does not make much sense to have a fixed support and a fixed \
support at the same position. The supports have been replaced by a single \
fixed support.";

Beam::"unitstep" = 
      "Warning: Equations containing UnitStep functions have to be solved. \
This might take quite long. If no result is returned within reasonable time, \
consider using numerical values for the positions. Results containing \
UnitStep functions in the denominator can not be trusted.";

Beam::"badstiffness" = 
      "The bending stiffness is ill-defined. Either use a constant function \
(e.g., '1&'), or a function containing UnitSteps(e.g., '1+3 \
UnitStep[#-.5]&'), or a list containing an expression and the variable (e.g., \
'{1,x}' or '{1+3 UnitStep[x-.5],x}').";

Beam::"baddistforce" = 
      "The distributed force is ill-defined. Either give a pure function \
(e.g., '#^2&'), or a list containing an expression and the variable (e.g., \
'{x^2,x}'.";

Beam::"nolast" = 
      "The last element cannot be removed. It is only possible to remove the \
last definition of a support or a load.";

Definitions

Begin Private Context

Begin["`Private`"];

Variable for Work Objects

$WorkBeam is the work beam.

$WorkBeam = Beam[{0, 1, 1&}, {Fixed[0], Force[1, 1]}];

$WorkSolution is the work solution.

$WorkSolution
       = {{1&, 1 - #1&, -(1/2) (-2 + #1) #1&, 
          -(1/6) (-3 + #1) #(1^2)&}, 
        Beam[{0, 1, 1&}, {Fixed[0], Force[0, -1], Moment[0, -1], 
            Force[1, 1]}]};

$LastElements is the support or load defined last.

$LastElements = {};

$WorkNotebook is the work notebook.

$WorkNotebook = Null;

$CalculationNotebook is the notebook for the information about a runnung calculation.

$CalculationNotebook = Null;

$Calculated is a global variable. If set True, the solution $WorkSolution of the beam $WorkBeam has already been calculated.

$Calculated = True;

$Calculated is a global variable. If True, each definition automatically calls a DrawBeam.

$DrawEach = True;

$DrawEach is a global variable. If True, all input is read with Input[] dialogs

$InputDialogs = True;

General Auxiliary Functions

MakePure[expression,x] constructs a pure function for the given expression in the variable x.

MakePure[expr_, x_] := Evaluate[expr /. x -> #1]&

MakeExpr[purefunction] applies the pure function to the string 'x'. This is useful only for printing results.

MakeExpr[fct_] := fct["x"]

Neg[fct] returns the negative pure function.

Neg[fct_Function] := Module[{x}, MakePure[-fct[x], x]]

InvertNumber[num] calculates -num for numbers.

SetAttributes[InvertNumber, Listable]

InvertNumber[x_] := If[NumberQ[x], -x, x]

InsertUnknowns[beam] inserts unknown forces and moments at the supports.

InsertUnknowns[beam_Beam] := 
      Beam[beam[[1]], 
        Flatten[beam[[2]] /. {
              Fixed[x_] -> {Fixed[x], Force[x], Moment[x]}, 
              Simple[x_] -> {Simple[x], Force[x]}}]]

InvertBeta[expression,x] handles expressions containing UnitStep functions linearly. It calculates 1/expression in terms of UnitStep functions.

\(InvertBeta[expr_, x_] := 
      Module[{e = Collect[expr, UnitStep[_]], free, us, ic, in}, 
        free = e /. UnitStep[_] -> 0; 
		
        us = Sort[
            Cases[e, UnitStep[_], 
              ∞], #1[[1, 1
                  ]] > #2[[1, 1
                  ]]&]; 
		
        ic = 1\/FoldList[Plus, free, \((Coefficient[e, #1]&)\)/@us]; 
        in = RotateRight[ReplacePart[ic, 0, \(-1\)]]; us = Prepend[us, 1]; 
        Simplify[Plus@@\((\((ic - in)\)\ us)\)]]\)

CheckBeam[beam] performs some consistency tests on the endpoints and the positions of supports and loads. If possible, inconsistencies are removed automatically.

\(CheckBeam[beam_Beam] := 
      Module[{ep = Flatten[beam[[1]]], 
          sl = Flatten[beam[[2]]], 
          mima}, Which[MatchQ[ep, {_?AtomQ, _?AtomQ, _Function}], Null, 
          MatchQ[ep, {_?AtomQ, _?AtomQ}], AppendTo[ep, 1&], True, 
          Message[Beam::"illdef"]; Throw[Null, "error"]]; 
        If[Length[sl] =!= Length[sl = Union[sl]], Message[Beam::"multelem"]]; 
        If[Length[sl] =!= Length[Union[First/@sl]], 
          If[sl =!= 
              \((sl = Union[
                    sl //. {a___, Fixed[x_], b___, Hinge[x_], c___} 
                        -> {a, b, c, Hinge[x], Simple[x]}])\), 
            Message[Beam::"fixedhinge"]]; 
          If[sl =!= 
              \((sl = Union[
                    sl //. {a___, Fixed[x_], b___, Simple[x_], c___} 
                        -> {a, b, c, Fixed[x]}])\), 
            Message[Beam::"simplefixed"]]]; 
		
        sl = DeleteCases[sl, 
            Hinge[ep[[1]]] | 
              Hinge[ep[[2]]]]; 
        If[\(\(! 
            \((And@@\(NumberQ/@Drop[ep, \(-1\)]\))\)\)\), 
          Message[Beam::"symbends"]]; 
        If[\(\(! \((Or@@\(NumberQ/@ep\))\)\)\), 
          Return[Beam[ep, sl]]]; 
        mima = Through[{Min, Max}[
              Select[Flatten[{First/@sl, ep}], NumberQ]]]; 
        If[\(\(! 
            \((ep[[1]] ≤ 
                mima[[1]])\)\)\), 
          Message[Beam::"badends"]; 
          ep[[1]] = 
            mima[[1]]]; 
		
        If[\(\(! 
            \((ep[[2]] ≥ 
                mima[[2]])\)\)\), 
          Message[Beam::"badends"]; 
          ep[[2]] = 
            mima[[2]]]; 
		
        sl = DeleteCases[sl, 
            Hinge[ep[[1]]] | 
              Hinge[ep[[2]]]]; 
        Beam[ep, sl]]\)

GetInput[dialogStrings,notebookString,template] pastes the notebookString and the template into the work notebook ($InputDialogs==False), or uses the dialogStrings to get the input via Input[], pastes the template into the work notebook, replaces the placeholdes by the input values (in the sequence given), and evaluates the expression.

\(GetInput[dialogStrings_List, notebookString_String, template_String] := 
      
	Module[{values}, 
		CheckNotebooks[]; 
		
        If[$InputDialogs, 
          values = \(ToString[#, InputForm]&\)/@\(Input/@dialogStrings\); 
          NotebookWrite[$WorkNotebook, template, Placeholder]; 
          \((NotebookWrite[$WorkNotebook, #1, Placeholder]&)\)/@values; 
          SelectionMove[$WorkNotebook, All, Cell]; 
          SelectionEvaluateCreateCell[$WorkNotebook], 
          NotebookWrite[$WorkNotebook, Cell[notebookString, "Text"]]; 
          NotebookWrite[$WorkNotebook, template, Placeholder]]]\)

AddElement[beam,element] adds the element to the beam.

SetAttributes[AddElement, HoldFirst]

AddElement[beam_, element_] := 
      Module[{res}, 
        If[Head[beam] =!= Beam, Message[Beam::"nobeam"]; Return[]]; 
        res = If[MatchQ[beam, Beam[_List, _List]], 
            Beam[beam[[1]], 
              Union[Append[beam[[2]], 
                  element]]], Beam[{}, element]]; 
        If[ValueQ[beam] && FreeQ[Hold[beam], Out], beam = res, res]]

AddToWorkBeam[fct,params] applies fct[params] to the work beam, sets $Calculated to False, and draws the work beam, if necessary.

\(AddToWorkBeam[fct_, params___] := 
      \((fct[$WorkBeam, params]; 
		$Calculated = False; 
		
        AppendTo[$LastElements, 
			
          Switch[fct, 
AddSimpleSupport, Simple[params], 
AddFixedSupport, 
            Fixed[params], 
AddHinge, Hinge[params], 
AddForce, 
            If[FreeQ[{params}, List], Force[params], 
					
              Force[{params}[[1]], 
                If[MatchQ[{params}[[2
                      ]], {_, _?AtomQ}], 
                  MakePure@@{params}[[2
                      ]], {params}[[2
                    ]]]]], 
AddMoment, Moment[params]]]; 
        
			If[$DrawEach, DrawBeam[]])\)\)

If[MatchQ[fct, {_, _?AtomQ}], f = MakePure@@fct, f = fct]

SolveFirst[beam,fct] opens the calculation window, solves the beam, closes the calcuation window, and applies fct to the result of the calculation.

SolveFirst[beam_Beam, fct_, opts___] := 
      Module[{s}, InitCalculationWindow[]; 
        s = Catch[DoSolveBeam[beam], "error"]; CloseCalculationWindow[]; 
        If[s =!= Null, fct[s, opts]]]

SolveAndPlot[plotfunction,opts] solves the work beam, if necessary, and applies the the given plotfunction.

SolveAndPlot[fct_, opts___] := 
      (If[! $Calculated, DoSolveBeam[]]; 
        If[$WorkSolution =!= Null, fct[$WorkSolution, opts]])

SolveAndEvaluate[resultType] solves the work beam, if necessary, and evaluates resultType[].

SolveAndEvaluate[resultType_] := 
      (If[! $Calculated, DoSolveBeam[]]; 
        If[$WorkSolution =!= Null, resultType[$WorkSolution]])

UnitStepToLists[fct,leftend,rightend] takes a pure function and two end points. It returns a list containing lists of values of the function between step intervals and the intervals themselves.

\(UnitStepToLists[fct_Function, x_, l_, r_] := 
      Module[{us, e}, e = fct[x]; 
		
        us = Sort[
            Union[Cases[e, UnitStep[_?\((\(! FreeQ[#, x]\)&)\)], 
                ∞]], #1[[1, 1
                  ]] > #2[[1, 1
                  ]]&]; 
		
        Transpose[{
            \(\(FoldList[Plus, 
                    e /. UnitStep[_?\((\(! FreeQ[#, x]\)&)\)] -> 0, 
                    \((Coefficient[e, #1]&)\)/@us] // Expand\) // Chop\) // 
              Simplify, 
            Partition[
              Join[{l}, 
                us /. UnitStep[a__?\((\(! FreeQ[#, x]\)&)\)] :> 
                    \((\(-a\) /. x -> 0)\), {r}], 2, 1]}]]\)

PrintResult[function,name,left,right] prints the result function in a nice form. The left and right end points must be given.

PrintResult[fct_, name_String, l_, r_] := 
      Module[{x}, 
        If[FreeQ[fct, UnitStep[_]], 
          CellPrint[
            Cell[BoxData[
                GridBox[{{name, "=", ToBoxes[fct["x"]], ",", l, "<", "x", 
                      "<", r}}]], "StandardForm"]], 
          CellPrint[
            Cell[BoxData[
                RowBox[{name, "=", "\[LeftBracketingBar]", 
                    GridBox[
                      ({ToBoxes[#1[[1
                                    ]] /. x -> "x"], 
                              ",", ToBoxes[#1[[2, 1
                                  ]]], "<", "x", "<", 
                              ToBoxes[#1[[2, 2
                                  ]]]}&)/@
                        UnitStepToLists[fct, x, l, r]]}]], "StandardForm"]]]]

CheckNotebooks[] opens the work notebook $WorkNotebook if no InputNotebook is available.

CheckNotebooks[] := 
      (If[! 
            (MemberQ[Notebooks[], InputNotebook[]]), $WorkNotebook = 
            NotebookCreate[WindowTitle -> "Work Beam"], $WorkNotebook = 
            InputNotebook[]]; )

InitCalculationWindow[] sets up a window indicating that the solution is being calculated.

InitCalculationWindow[] := 
      ($CalculationNotebook = 
          NotebookCreate[
            WindowMargins -> {{0, Automatic}, {Automatic, 0}}, 
            WindowSize -> {300, 100}, 
            WindowTitle -> "Calculating..."]; 
        NotebookWrite[$CalculationNotebook, 
          Cell["Solving the beam. Please hold on...", "Text"]]; )

CloseCalculationWindow[] closes the calculation window.

CloseCalculationWindow[] := NotebookClose[$CalculationNotebook];

Access to Local Variables

SetInputDialogs[value_?(MemberQ[{True, False}, #1]&)] := 
      $InputDialogs = value;

SetDrawEach[value_?(MemberQ[{True, False}, #1]&)] := 
      $DrawEach = value;

SetWorkBeam[beam_Beam] := 
      ($Calculated = False; $WorkBeam = beam)

GetWorkBeam[] := $WorkBeam

SetWorkSolution[x : {{_, _, _, _}, _Beam}] := $WorkSolution = x

GetWorkSolution[] := $WorkSolution

Constructing Beams

SetAttributes[{NewBeam, DefineEndPoints, DefineBendingStiffness, 
        AddSimpleSupport, AddFixedSupport, AddHinge, AddForce, AddMoment, 
        RemoveElement}, {HoldFirst}]

BeamPasteEvaluate[fct_String] := 
      (CheckNotebooks[]; NotebookWrite[$WorkNotebook, fct, All]; 
        SelectionEvaluateCreateCell[$WorkNotebook])

NewBeam[] := Beam[{0, 1, 1&}, {}]

NewBeam[beam_] := 
      If[(ValueQ[beam] || AtomQ[beam]) && FreeQ[Hold[beam], Out], 
        beam = Beam[{0, 1, 1&}, {}], Beam[{0, 1, 1&}, {}]]

NewWorkBeam[] := 
      ($WorkBeam = Beam[{0, 1, 1&}, {}]; $Calculated = False; $LastElements
           = {}; )

DefineEndPoints[beam_, {left_, right_}] := 
      Module[{res}, 
        If[Head[beam] =!= Beam, Message[Beam::"nobeam"]; Return[]]; 
        res = If[MatchQ[beam, Beam[{_, _, _}, _List]], 
            ReplacePart[
              beam, {left, right, 
                beam[[1, 3]]}, 1], 
            Beam[{left, right, 1&}, 
              beam[[2]]]]; 
        If[ValueQ[beam] && 
            ! 
              (MemberQ[Attributes[beam], Protected]), beam = res, 
          res]]

DefineEndPoints[{left_, right_}] := 
      (DefineEndPoints[$WorkBeam, {left, right}]; $Calculated = False; )

DefineEndPoints[] := 
      GetInput[{"Enter the position of the left end of the beam", 
          "Enter the position of the right end of the beam"}, 
        "Replace the placeholders in the following input cell by the \
positions of the left and the right end points of the beam. Then evaluate the \
cell.", "DefineEndPoints[{\[Placeholder],\[Placeholder]}]"]

DefineBendingStiffness[beam_, sf_] := 
      Module[{fct, res}, 
        If[Head[beam] =!= Beam, Message[Beam::"nobeam"]; Return[]]; 
        If[MatchQ[sf, {_, _?AtomQ}], fct = MakePure@@sf, fct = sf]; 
        If[Head[fct] =!= Function, Message[Beam::"badstiffness"]; Return[b]]; 
        res = If[MatchQ[beam, Beam[{_, _, _}, _List]], 
            ReplacePart[beam, fct, {1, 3}], 
            Beam[{0, 1, fct}, 
              beam[[2]]]]; 
        If[ValueQ[beam] && 
            ! 
              (MemberQ[Attributes[beam], Protected]), beam = res, 
          res]]

DefineBendingStiffness[sf_] := 
      (DefineBendingStiffness[$WorkBeam, sf]; $Calculated = False; )

DefineBendingStiffness[] := 
      GetInput[{
          "Enter a pure function defining the bending stiffness (e.g., \
(1+.2UnitStep[#-.5])&) or a list containing an expression defining the \
bending stiffness and the variable (e.g., {1+.2UnitStep[x-.5], x})"}, 
        "Replace the placeholder in the following input cell either by a pure \
function defining the bending stiffness (e.g., (1+.2UnitStep[#-.5])&), or by \
a list containing an expressing defining the bending stiffness and the \
variable (e.g., {1+.2UnitStep[x-.5], x}). Then evaluate the cell.", 
        "DefineBendingStiffness[\[Placeholder]]"]

AddSimpleSupport[beam_, pos_] := AddElement[beam, Simple[pos]]

AddSimpleSupport[pos_] := AddToWorkBeam[AddSimpleSupport, pos]

AddSimpleSupport[] := 
      GetInput[{"Enter the position of the simple support"}, 
        "Replace the placeholder in the following input cell by the position \
of the simple support. Then evaluate the cell.", 
        "AddSimpleSupport[\[Placeholder]]"]

AddFixedSupport[beam_, pos_] := AddElement[beam, Fixed[pos]]

AddFixedSupport[pos_] := AddToWorkBeam[AddFixedSupport, pos]

AddFixedSupport[] := 
      GetInput[{"Enter the position of the fixed support"}, 
        "Replace the placeholder in the following input cell by the position \
of the fixed support. Then evaluate the cell.", 
        "AddFixedSupport[\[Placeholder]]"]

AddHinge[beam_, pos_] := AddElement[beam, Hinge[pos]]

AddHinge[pos_] := AddToWorkBeam[AddHinge, pos]

AddHinge[] := 
      GetInput[{"Enter the position of the hinge"}, 
        "Replace the placeholder in the following input cell by the position \
of the hinge. Then evaluate the cell.", "AddHinge[\[Placeholder]]"]

\(AddForce[beam_, {l_, r_}, fct_] := 
      Module[{f}, If[MatchQ[fct, {_, _?AtomQ}], f = MakePure@@fct, f = fct]; 
        
		If[Head[f] =!= Function, Message[Beam::"baddistforce"]; 
          Return[beam]]; AddElement[beam, Force[{l, r}, f]]]\)

AddForce[{l_, r_}, fct_] := AddToWorkBeam[AddForce, {l, r}, fct]

AddDistributedForce[] := 
      GetInput[{"Enter the the left end of the distributed force", 
          "Enter the right end of the distributed force", 
          "Enter the function defining the force. Either use a pure function \
(e.g., #^2& or Function[x,x^2] for a quadratic force distribution), or use a \
list with an expression and the variable (e.g., {x^2,x}); forces pointing \
downwards are positive"}, 
        "Replace the first two placehoders by the left and the right end of \
the distributed force. The third placeholder should either contain a pure \
function defining the force (e.g., #^2& or Function[x,x^2] for a quadratic \
distribution), or a list containing an expression and the variable (e.g., \
{x^2,x}). Forces pointing downwards are positive. Then evaluate the cell.", 
        "AddForce[{\[Placeholder],\[Placeholder]},\[Placeholder]]"]

AddForce[beam_, pos_, value_] := AddElement[beam, Force[pos, value]]

AddForce[pos_, value_] := AddToWorkBeam[AddForce, pos, value]

AddForce[] := 
      GetInput[{"Enter the position of the force", 
          "Enter the value of the force (the positive force direction is \
downwards)"}, 
        "Replace the placeholders in the following input cell by the \
position and the value of the force. Forces pointing downwards are positive. \
Then evaluate the cell.", "AddForce[\[Placeholder],\[Placeholder]]"]

AddMoment[beam_, pos_, value_] := 
      AddElement[beam, Moment[pos, value]]

AddMoment[pos_, value_] := AddToWorkBeam[AddMoment, pos, value]

AddMoment[] := 
      GetInput[{"Please enter the position of the moment", 
          "Enter the value of the moment (moment vectors pointing into the \
screen are positive)"}, 
        "Replace the placeholders in the following input cell by the position \
and the value of the moment. Moment vectors pointing into the screen are \
positive. Then evaluate the cell.", 
        "AddMoment[\[Placeholder],\[Placeholder]]"]

\(RemoveLastElement[] := 
	
      If[$LastElements === {}, Message[Beam::"nolast"], 

        		$WorkBeam = 
          DeleteCases[$WorkBeam, $LastElements
              [[\(-1\)]], 
            ∞]; $LastElements = Drop[$LastElements, \(-1\)]; 

        		$Calculated = False; 
		If[$DrawEach, DrawBeam[]]]\)

RemoveElement[beam_, element_] := 
      Module[{res}, 
        If[Head[beam] =!= Beam, Message[Beam::"nobeam"]; Return[]]; 
        res = Beam@@
            DeleteCases[beam /. Beam -> List, element, ∞]; 
        If[ValueQ[beam] && 
            ! 
              (MemberQ[Attributes[beam], Protected]), beam = res, 
          res]]

Solving Beams

DoSolveBeam is the actual function which solves a beam. It is defined with a separate name to simplify catching errors.

\(DoSolveBeam[beam_Beam] := 
      Module[{b, S, M, d, s, x, l, r, ef, if, betaInverse, c, unknownFM, 
          newUnknownFM, equations, sol, res, cleanRule}, 
        b = InsertUnknowns[CheckBeam[beam]]; 
		
        l = b[[1, 1]]; 
        r = b[[1, 2]]; 
		
        ef = Cases[b[[2]], 
            Force[l, ___] | Force[r, ___] | Moment[l, ___] | Moment[r, ___]]; 
        if = Complement[b[[2]], ef]; 
        cleanRule
           = {UnitStep[a_?Positive] -> 1, 
            UnitStep[b_?Negative] -> 0, 
            UnitStep[a_?Positive\ b_?Negative] -> 0, 
            DiracDelta[_?\((FreeQ[#1, x]&)\)] -> 0}; 
        betaInverse = 
          If[FreeQ[b[[1, 3]], 
              UnitStep[_]], 
            1\/\(b[[1, 3]]\)[x], 
            InvertBeta[\(b[[1, 3]]\)[x], 
              x]]; S = 
          \(-Plus@@Cases[ef, Force[l, v_: Force[l]] -> v]\) - 
            ∫Plus@@
                  \((Cases[if, Force[{_, _}, _]] /. {
                        Force[{p1_, r}, f_] -> UnitStep[x - p1]\ f[x], 
                        Force[{p1_, p2_}, f_] -> 
                          UnitStep[x - p1]\ f[x] - UnitStep[x - p2]\ f[x]})
                    \)\[DifferentialD]x - 
            Plus@@Join[
                Cases[if, Force[a_] -> Force[a]\ UnitStep[x - a]], 
                Cases[if, 
                  Force[a_?\((\(\(! \((ListQ[#1])\)\)\)&)\), 
                      v_] -> v\ UnitStep[x - a]]]; 
        equations
           = {\((S /. x -> r)\) - 
                Plus@@Cases[ef, Force[r, v_: Force[r]] -> v] == 0}; 
        M = \(-Plus@@Cases[ef, Moment[l, v_: Moment[l]] -> v]\) - 
            ∫S\[DifferentialD]x - 
            Plus@@Join[
                Cases[if, Moment[a_] -> Moment[a]\ UnitStep[x - a]], 
                Cases[if, Moment[a_, v_] -> v\ UnitStep[x - a]]]; 
        AppendTo[equations, 
          \((M /. x -> r)\) - 
              Plus@@Cases[ef, Moment[r, v_: Moment[r]] -> v] == 0]; 
        s = ∫\(M\ betaInverse\)\[DifferentialD]x + c[1] + 
            Plus@@Cases[if, Hinge[a_] -> c[a, a]\ UnitStep[x - a]]; 
        d = ∫s\[DifferentialD]x + c[2]; 
        equations = 
          Join[equations, 
              Flatten[Cases[if, Simple[_] | Fixed[_] | Hinge[_]] /. {
                    Simple[a_] -> \((d /. x -> a)\) == 0, 
                    Fixed[a_] -> 
                      Thread[\((
                            NestList[∂\_x #1&, d, 1] /. x -> a)
                            \) == 0], 
                    Hinge[a_] -> \((M /. x -> a)\) == 0}]] /. 
            cleanRule; 
        If[\(\(! \((FreeQ[equations, UnitStep])\)\)\), 
          Message[Beam::"unitstep"]]; 
        sol = Solve[equations, 
            Union[Flatten[{Cases[d, c[__], ∞], 
                  Cases[b[[2]], 
                    Force[_] | Moment[_], ∞]}]]]; 
        If[sol == {}, Message[Beam::"undef"]; Throw[Null, "error"]]; 
        unknownFM = 
          Cases[b[[2]], 
            Force[_] | Moment[_]]; 
        newUnknownFM = 
          \(unknownFM
               /. {Force[a_] -> Force[a, Force[a]], 
                Moment[a_] -> Moment[a, Moment[a]]}\) /. 
            Chop[sol[[1]]]; 
        res = Simplify[
            Chop[Expand[{S, M, s, d} /. 
                  sol[[1]]]]]; 
        If[\(\(! \((FreeQ[res, c[__]])\)\)\), 
          Message[Beam::"nosol"]; Throw[Null, "error"]]; 

        		{\((MakePure[#1, x]&)\)/@res, 
          b /. Thread[unknownFM -> newUnknownFM]}]\)

DoSolveBeam[] := 
      (InitCalculationWindow[]; $WorkSolution = 
          Catch[DoSolveBeam[$WorkBeam], "error"]; 
        If[$WorkSolution =!= Null, $Calculated = True]; 
        CloseCalculationWindow[]; )

SolveBeam[beam_Beam] := Catch[DoSolveBeam[beam], "error"]

SolveBeam[] := 
      (DoSolveBeam[]; If[$WorkSolution =!= Null, $WorkSolution])

Plots of Results

PlotSolution contains common code of the plot functions.

\(PlotSolution[solution_, l_, r_, xLabel_String, title_String, opts___]
         := Module[{options, fs}, 
        If[\(\(! \((NumberQ[l] && NumberQ[r])\)\)\), 
          Message[Beam::"numends"]; Return[Null]]; 
        fs = \((\(DefaultFont /. Flatten[{opts}]\) /. 
                  DefaultFont -> $DefaultFont)\)[[2
              ]] + 2; 
        options = 
          Flatten[{opts, 
					
              PlotLabel -> 
                StyleForm[xLabel, FontFamily -> "Times", 
                  FontSize -> fs, FontWeight -> Bold], 
              FrameLabel 
                -> {"", "", "", 
                  StyleForm[title, FontFamily -> "Times", 
                    FontSlant -> Italic, FontSize -> fs]}, 
              Frame -> True, 
					RotateLabel -> False, 

              					PlotRange -> All}]; 
        Plot[solution[x], {x, l, r}, Evaluate[options]]]\)

PlotShearingForce[{{S_, _, _, _}, Beam[{l_, r_, _}, _]}, opts___] := 
      PlotSolution[S, l, r, "Shearing Force", "\!\(S\_y\)", opts]

PlotShearingForce[beam_Beam, opts___] := 
      SolveFirst[beam, PlotShearingForce, opts]

PlotShearingForce[opts___] := SolveAndPlot[PlotShearingForce, opts]

PlotBendingMoment[{{_, M_, _, _}, Beam[{l_, r_, _}, _]}, opts___] := 
      PlotSolution[M, l, r, "Bending Moment", "\!\(M\_z\)", opts]

PlotBendingMoment[beam_Beam, opts___] := 
      SolveFirst[beam, PlotBendingMoment, opts]

PlotBendingMoment[opts___] := SolveAndPlot[PlotBendingMoment, opts]

PlotSlope[{{_, _, s_, _}, Beam[{l_, r_, _}, _]}, opts___] := 
      PlotSolution[Neg[$WorkSolution], l, r, "Slope", "-d'", opts]

PlotSlope[{{_, _, s_, _}, Beam[{l_, r_, _}, _]}, opts___] := 
      Module[{p, t, df}, 
        df = DisplayFunction /. Flatten[{opts}] /. 
            DisplayFunction -> $DisplayFunction; 
        p = PlotSolution[Neg[s], l, r, "Slope", "\!\(d'\)", 
            Flatten[{opts, DisplayFunction -> Identity}]]; 
        t = FrameTicks /. FullOptions[p]; 
        Show[p, FrameTicks -> 
            ReplacePart[t, 
              Transpose[
                MapAt[InvertNumber, 
                  Transpose[t[[2]]], 
                  2]], 2], DisplayFunction -> df]]

PlotSlope[beam_Beam, opts___] := SolveFirst[beam, PlotSlope, opts]

PlotSlope[opts___] := SolveAndPlot[PlotSlope, opts]

\(PlotDeflection[{{_, _, _, d_}, Beam[{l_, r_, _}, _]}, opts___] := 
      Module[{p, t, df}, 
        df = \(DisplayFunction /. Flatten[{opts}]\) /. 
            DisplayFunction -> $DisplayFunction; 
        p = PlotSolution[Neg[d], l, r, "Deflection", "d", 
            Flatten[{opts, DisplayFunction -> Identity}]]; 
        t = FrameTicks /. FullOptions[p]; 
	
        t = FrameTicks /. FullOptions[p]; 
        Show[p, FrameTicks -> 
            ReplacePart[t, 
              Transpose[
                MapAt[InvertNumber, 
                  Transpose[t[[2]]], 
                  2]], 2], DisplayFunction -> df]]\)

PlotDeflection[beam_Beam, opts___] := 
      SolveFirst[beam, PlotDeflection, opts]

PlotDeflection[opts___] := 
      (If[! $Calculated, DoSolveBeam[]]; 
        PlotDeflection[$WorkSolution, opts])

PlotBendingStiffness[Beam[{l_, r_, sf_}, _], opts___] := 
      PlotSolution[sf, l, r, "Bending Stiffness", "\!\(E I\_z\)", opts]

PlotBendingStiffness[{{_, _, _, _}, beam_}, opts___] := 
      PlotBendingStiffness[beam, opts]

PlotBendingStiffness[opts___] := 
      PlotBendingStiffness[$WorkBeam, opts]

Symbolic Results

PrintShearingForce[{{S_, _, _, _}, Beam[{l_, r_, _}, _]}] := 
      PrintResult[S, "\!\(S\_y\)", l, r]

PrintShearingForce[beam_Beam] := 
      SolveFirst[beam, PrintShearingForce]

PrintShearingForce[] := SolveAndEvaluate[PrintShearingForce]

ShearingForce[{{S_, _, _, _}, Beam[{_, _, _}, _]}] := S

ShearingForce[beam_Beam] := SolveFirst[beam, ShearingForce]

ShearingForce[] := SolveAndEvaluate[ShearingForce]

PrintBendingMoment[{{_, M_, _, _}, Beam[{l_, r_, _}, _]}] := 
      PrintResult[M, "\!\(M\_z\)", l, r]

PrintBendingMoment[beam_Beam] := 
      SolveFirst[beam, PrintBendingMoment]

PrintBendingMoment[] := SolveAndEvaluate[PrintBendingMoment]

BendingMoment[{{_, M_, _, _}, Beam[{_, _, _}, _]}] := M

BendingMoment[bean_Beam] := SolveFirst[beam, BendingMoment]

BendingMoment[] := SolveAndEvaluate[BendingMoment]

PrintSlope[{{_, _, s_, _}, Beam[{l_, r_, _}, _]}] := 
      PrintResult[s, "d'", l, r]

PrintSlope[beam_Beam] := SolveFirst[beam, PrintSlope]

PrintSlope[] := SolveAndEvaluate[PrintSlope]

Slope[{{_, _, s_, _}, Beam[{_, _, _}, _]}] := s

Slope[beam_Beam] := SolveFirst[beam, Slope]

Slope[] := SolveAndEvaluate[Slope]

PrintDeflection[{{_, _, _, d_}, Beam[{l_, r_, _}, _]}] := 
      PrintResult[d, "d", l, r]

PrintDeflection[beam_Beam] := SolveFirst[beam, PrintDeflection]

PrintDeflection[] := SolveAndEvaluate[PrintDeflection]

Deflection[{{_, _, _, d_}, Beam[{_, _, _}, _]}] := d

Deflection[beam_Beam] := SolveFirst[beam, Deflection]

Deflection[] := SolveAndEvaluate[Deflection]

PrintBendingStiffness[Beam[{l_, r_, bs_}, _]] := 
      PrintResult[bs, "E\!\(I\_z\)", l, r]

PrintBendingStiffness[{{_, _, _, _}, Beam[{l_, r_, bs_}, _]}] := 
      PrintResult[bs, "E\!\(I\_z\)", l, r]

PrintBendingStiffness[] := PrintBendingStiffness[$WorkBeam]

BendingStiffness[Beam[{_, _, bs_}, _]] := bs

BendingStiffness[{{_, _, _, _}, Beam[{_, _, bs_}, _]}] := bs

BendingStiffness[] := BendingStiffness[$WorkBeam]

Auxiliary Functions for Drawings

MakeFixed[pos,as] returns graphics objects for a fixed support at the given position with a scaling as.

MakeFixed[pos_?NumberQ, as_] := {GrayLevel[.5], 
        Rectangle[{pos - as, -2 as}, {pos + as, 2 as}]}

MakeFixed[pos_, as_] := {Text[Fixed[pos], Scaled[{.5, .1}]]}

MakeSimple[pos, as] returns graphics objects for a simple support at the given position with a scaling as.

MakeSimple[pos_?NumberQ, as_] := {
        Line[{{pos - 2 as, -2 as}, {pos + 2 as, 2 as}}], 
        Line[{{pos - 2 as, 2 as}, {pos + 2 as, -2 as}}], {
          GrayLevel[.5], 
          Rectangle[{pos - 2 as, -4 as}, {pos + 2 as, -2 as}], 
          Rectangle[{pos - 2 as, 4 as}, {pos + 2 as, 2 as}]}}

MakeSimple[pos_, as_] := {Text[Simple[pos], Scaled[{.5, .2}]]}

MakeHinge[pos,as] returns graphics objects for a hinge at the given position with a scaling as and an optional displacement dpos.

MakeHinge[pos_?NumberQ, as_, dpos_: 0]
         := {{GrayLevel[1], Disk[{pos, dpos}, as]}, Circle[{pos, dpos}, as]}

MakeHinge[pos_, as_, dpos_: 0] := {Text[Hinge[pos], Scaled[{.5, .3}]]}

MakeForce[pos,value,as,fs,dpos] returns graphics objects for a force at pos with value, absolute scaling as, force scaling fs an optional displacement dpos (for solutions).

MakeForce[pos_?NumberQ, value_, as_, fs_, dpos_: 0] := 
      Module[{v}, 
        v = If[NumberQ[value], value, 1]; {
          Arrow[{pos, dpos + fs v}, {pos, dpos}, 
            HeadScaling -> Relative], 
          Text[value, {pos, dpos + fs v + 2 Sign[v] as}]}]

MakeForce[pos_, value_, as_, fs_, dpos_: 0] := {
        Text[Force[pos, value], Scaled[{.5, .9}]]}

MakeDistributedForce[{left,right},ffct,as,fs,dpos] returns graphics objects for a disributed force in the interval [left,right] with force function ffct, absolute scaling as, force scaling fs an optional displacement dpos (for solutions).

MakeDistributedForce[{l_?NumberQ, r_?NumberQ}, 
        ffct_?(NumberQ[#[1]]&), as_, fs_, dfct_:(0&)] := 
      Module[{x, p = ((l + r)/2), fp}, 
        fp = ffct[p]; {Line[{{l, dfct[l] + fs ffct[l]}, {l, dfct[l]}}], 
          Line[{{r, dfct[r] + fs ffct[r]}, {r, dfct[r]}}], 
          Plot[fs ffct[x] + dfct[x], {x, l, r}, 
              DisplayFunction -> Identity][[1, 1
            ]], 
          Text[TraditionalForm[ffct[x] /. x -> "x"], {p, 
              dfct[p] + fs fp + 3 as Sign[fp]}]}]

MakeDistributedForce[{l_, r_}, ffct_, as_, fs_, dfct_:(0&)] := {
        Text[Force[{l, r}, ffct], Scaled[{.5, .8}]]}

MakeMoment[pos,value,as,fs,dpos] returns graphics objects for a moment at pos with value, absolute scaling as, force scaling fs an optional displacement dpos (for solutions).

MakeMoment[pos_?NumberQ, value_, as_, dpos_: 0] := {
        Circle[{pos, dpos}, as], 
        If[! (NumberQ[value]) || value > 0, {
            Line[{{pos - 0.7 as, dpos - 0.7 as}, {pos + 0.7 as, 
                  dpos + 0.7 as}}], 
            Line[{{pos - 0.7 as, dpos + 0.7 as}, {pos + 0.7 as, 
                  dpos - 0.7 as}}]}, Disk[{pos, dpos}, .3 as]], 
        Text[value, {pos + 1.8 as, dpos + 1.8 as}]}

MakeMoment[pos_, value_, as_, dpos_: 0] := {
        Text[Moment[pos, value], Scaled[{.5, .7}]]}

BeamGraphics[beam] returns all graphics objects for a beam.

\(BeamGraphics[Beam[{l_, r_, _ : \((1&)\)}, obj_], opts___] := 
      Module[{as, fs, options}, 
		
        If[\(! \((NumberQ[l] && NumberQ[r])\)\), Message[Beam::"numends"]; 
          Return[Null]]; 
        options = 
          Flatten[{opts, Axes -> None, AspectRatio -> Automatic, 
              PlotRange -> All}]; as = .02\ Abs[r - l]; 
        fs = \((10\ as)\)/
            Max[Flatten[{
                  Cases[obj, Force[__]] /. {
                      Force[x_?NumberQ, v_?NumberQ] -> Abs[v], 
                      Force[x_?NumberQ, _] -> 1, 
                      Force[{left_?NumberQ, right_?NumberQ}, 
                          fct_?\((NumberQ[#[1]]&)\)] -> 
                        Abs[{fct[left], fct[right]}], 
                      Force[___] -> 0}, .0001}]]; 
        Graphics[{{AbsoluteThickness[3], Line[{{l, 0}, {r, 0}}]}, 
            obj /. {Fixed[x_] :> MakeFixed[x, as], 
                Simple[x_] :> MakeSimple[x, as], 
                Hinge[x_] :> MakeHinge[x, as], 
                Force[ends_List, fct_] :> 
                  MakeDistributedForce[ends, fct, as, fs], 
                Force[x_, v_] :> MakeForce[x, v, as, fs], 
                Moment[x_, v_] :> MakeMoment[x, v, as]}}, 
          options]]\)

Drawings of Beams

\(DrawBeam[beam_Beam, opts___] := 
	
      Module[{res = Catch[BeamGraphics[CheckBeam[beam], opts], "error"]}, 
        If[res =!= Null, Show[res], Null]]\)

DrawBeam[{{_, _, _, _}, beam_Beam}, opts___] := 
      DrawBeam[N[beam], opts]

DrawBeam[] := 
      Module[{res = 
            Catch[BeamGraphics[$WorkBeam = CheckBeam[$WorkBeam]], "error"]}, 
        If[res =!= Null, Show[res], Null]]

DrawSolution[{{_, _, _, d_}, Beam[{l_, r_, _ : (1&)}, obj_]}, 
        opts___] := 
      Module[{ar, as, fs, p, options, max, scaledD}, 
        If[! (NumberQ[l] && NumberQ[r]), 
          Message[Beam::"numends"]; Return[Null]]; 
        options = 
          Flatten[{opts, Axes -> None, AspectRatio -> Automatic, 
              PlotRange -> All, 
              DisplayFunction -> $DisplayFunction}]; 
        as = .02 Abs[r - l]; 
        fs = (10 as)/
            Max[Flatten[{
                  Cases[obj, Force[__]] /. {
                      Force[x_?NumberQ, v_?NumberQ] :> Abs[v], 
                      Force[x_?NumberQ, _] -> 1, 
                      Force[{left_?NumberQ, right_?NumberQ}, 
                          fct_?(NumberQ[#[1]]&)] :> 
                        Abs[{fct[left], fct[right]}], 
                      Force[___] -> 0}, .0001}]]; 
        ar = (.3 (r - l))/
            Plus@@Through[{Max, -Min[#1]&}[
                  Append[Table[d[x], {x, l, r, ((r - l)/20)}], .000001]]]; 
        scaledD = Function[x, -.3 ar Evaluate[d[x]]]; 
        p = Plot[Evaluate[scaledD[x]], {x, l, r}, 
            Evaluate[
              Flatten[{PlotStyle -> AbsoluteThickness[3], 
                  DisplayFunction -> Identity, options}]]]; 
        Show[{p, 
            Graphics[{
                Cases[obj, Force[__] | Moment[__] | Hinge[__], ∞]
                     /. {Hinge[x_] :> 
                      MakeHinge[x, as, scaledD[x]], 
                    Force[ends_List, fct_] :> 
                      MakeDistributedForce[ends, fct, as, fs, scaledD], 
                    Force[x_, v_] :> 
                      MakeForce[x, N[v, Max[2, N[Log[10, Abs[v]]] + 1]], as, 
                        fs, scaledD[x]], 
                    Moment[x_, v_] :> 
                      MakeMoment[x, N[v, Max[2, N[Log[10, Abs[v]]] + 1]], as, 
                        scaledD[x]]}}]}, options]]

DrawSolution[beam_Beam] := SolveFirst[beam, DrawSolution]

DrawSolution[] := SolveAndEvaluate[DrawSolution]

Draw and Plot all Results

\(DrawAll[beam_Beam, sol : {{_, _, _, _}, _Beam}] := 
      Module[{df
             = {$DefaultFont[[1
                ]], $DefaultFont[[2
                  ]] - 1}}, 
        If[\(\(! 
            \((NumberQ[beam[[1, 1]]] && 
                NumberQ[beam[[1, 2]]])
              \)\)\), Message[Beam::"numends"]; Return[Null]]; 
        Show[GraphicsArray[{

              						{
                DrawBeam[beam, DisplayFunction -> Identity, 
                  DefaultFont -> df], 
							
                PlotShearingForce[sol, DisplayFunction -> Identity, 
                  AspectRatio -> .5, DefaultFont -> df]}, {
                DrawSolution[sol, DisplayFunction -> Identity, 
                  DefaultFont -> df], 
							
                PlotBendingMoment[sol, DisplayFunction -> Identity, 
                  AspectRatio -> .5, DefaultFont -> df]}, {
                PlotBendingStiffness[beam, DisplayFunction -> Identity, 
                  AspectRatio -> .5, DefaultFont -> df], 
                PlotSlope[sol, DisplayFunction -> Identity, 
                  AspectRatio -> .5, DefaultFont -> df]}}], 
          DisplayFunction -> $DisplayFunction, AspectRatio -> 1]]
        \)

DrawAll[beam_Beam] := SolveFirst[beam, DrawAll[beam, #1]&]

DrawAll[] := 
      (If[! $Calculated, DoSolveBeam[]]; 
        If[$WorkSolution =!= Null, DrawAll[$WorkBeam, $WorkSolution]])

Package Endings

End Private Context

End[];

Turn on Spelling Checker

On[General::"spell1"]

End Package

EndPackage[];