(* BeginPackage["BeamStatics`", "Calculus`DiracDelta`", "Utilities`FilterOptions`"] *)
BeginPackage["BeamStatics`", "Utilities`FilterOptions`"]
(* Originally written for Mathematica Version 2.2. This version is for
use with Mathematica Version 4.
Changes February 11, 2000 to make BeamStatics work with Mathematica Version 4.
(1) In the BeginPackage statement "Calculus`DiracDelta`" is not
included because Mathematica Version 4 includes the functionality
of the DiracDelta package automatically. Original BeginPackage on line 1
is commented out.
(2) The function SimplifyUnitStep is not called, because it is no longer
available. Lines 433 - 436 commented out, replaced with lines 438 - 439.
The package has no restrictions on use or
copying. Bug reports for the original package can be e-mailed to Levent
Kitis at lk3a@virginia.edu or lk3a@esinet.net *)
Beam::usage = "Beam[x, L, EI, bnd, options] calculates shear, bending
moment, slope, and deflection of a straight beam under static loads; x
measures position from the left end, EI is the list of flexural rigidities,
L is the length, bnd is a list of two strings giving boundary conditions at
the two ends of the beam."
DrawBeam::usage = "DrawBeam[x, L, EI, bnd, options] draws a sketch of the
beam and the loads acting on it. DrawBeam places the figure in a rectangle
with lower left corner at (0, 0) and upper right corner at (Xmax, Ymax),
whose default setting is (100, 130)."
urule::usage = "urule[x, a] replaces UnitStep[x] with 1 and UnitStep[k a] with
UnitStep[k] when k is a number."
uregion::usage = "uregion[x, L, expr] removes all expressions containing
UnitStep in expr and breaks expr into parts that are valid over continuous
subintervals of Interval[{0, L}]. For example, with expr = P UnitStep[x -
L/2], uregion[x, L, expr] gives {{0, P}, {{0, L/2}, {L/2, L}}}."
plots::usage = "plots[f1, f2, x, x0, xf, options] plots f1 and f2 as function of
x from x0 to xf in two frames on one page. Defaults are for f1 = deflection
and f2 = slope."
plotm::usage = "plots[f1, f2, x, x0, xf, options] plots f1 and f2 as function of
x from x0 to xf in two frames on one page. Defaults are for f1 = bending moment
and f2 = shear force."
PointLoad::usage = "PointLoad is an option of Beam and DrawBeam. PointLoad
-> {P1, a1, P2, a2, ...} specifies that a concentrated load of magnitude
P1 acts at x = a1, etc. In DrawBeam, P1 etc. are used to label the
concentrated forces."
LineLoad::usage = "LineLoad is an option of Beam and DrawBeam. LineLoad->
{p1, a, p2, b, q1, c, q2, d,...} specifies linearly distributed
loads of intensity p1 at x = a, p2 at x = b and q1 at x = c, q2 at x = d, etc."
InSpan::usage = "InSpan is an option of Beam and DrawBeam. Spanload -> {a, b,...}
specifies that there are rigid in-span supports at x = a, x = b, etc. InSpan -> a
can be used for a single support at x = a."
PointCouple::usage = "PointCouple is an option of Beam and DrawBeam. PointCouple
-> {M1, a1, M2, a2, ...} specifies that a concentrated moment of magnitude
M1 acts at x = a1, etc. In DrawBeam, M1 etc. are used to label the
concentrated moments."
LoadFunction::usage = "LoadFunction is an option of Beam and
DrawBeam. LoadFunction -> f[x] or LoadFunction -> {f[x], g[x], ...}
specifies distributed loads acting on the beam."
ShearConnection::usage = "ShearConnection is an option of Beam and
DrawBeam. ShearConnection -> {a, b,...} specifies that there are shear
connections, or hinges, at x = a, x = b, etc. ShearConnection -> a can be
used for a single hinge at x = a."
MomentConnection::usage = "MomentConnection is an option of Beam and
DrawBeam. MomentConnection -> {a, b,...} specifies that there are moment
connections at x = a, x = b, etc. MomentConnection -> a can be
used for a single moment connection at x = a."
LoadFunctionRange::usage = "LoadFunctionRange is an option of
DrawBeam. This defines the ranges over which the the list of functions
given by the LoadFunction option are nonzero. For example, LoadFunction ->
{f[x], g[x]} requires LoadFunctionRange -> {{a, b}, {c, d}} to specify that
f(x) is to be drawn between x = a and x = b over the span and g(x) between
x = c and x = d."
LoadFunctionTags::usage = "LoadFunctionTags is an option of DrawBeam. It is
a list of strings defining the labels of the loads specified by
LoadFunction. The labels are centered horizontally at approximately the
midpoints of the function ranges."
LoadFunctionTagLoc::usage = "LoadFunctionTagLoc is an option of
DrawBeam. It defines the vertical distance between the load function curve
and its label for functions specified by the LoadFunction option. The
default is 2.5."
LoadSubdivisions::usage = "LoadSubdivisions is an option of DrawBeam. The
number of subdivisions of the maximum horizontal range of the loads defined
by LineLoad. The other ranges are subdivided such that the length of the
subdivisions approximates that of the load with the maximum range. This
determines the horizontal spacing of the vertical arrows drawn to represent
distributed loads. The default setting is 5."
LoadFunctionSubdivisions::usage = "LoadFunctionSubdivisions is an option of
DrawBeam. It serves the same purpose as LoadSubdivisions for the loads
defined by LoadFunction. The default setting is 5."
LengthScale::usage = "LengthScale is an option of DrawBeam. LengthScale ->
s scales all horizontal lengths x to x s/L units where L is the actual
length of the beam. The default value of s is 80. "
HeightScale::usage = "HeightScale is an option of DrawBeam. If EI is
constant, a rectangle of width LengthScale and height HeightScale is drawn
to represent the beam. If EI is not constant, the portion of the beam with
the maximum EI has height HeightScale. Heights elsewhere are proportional
to the flexural rigidity EI. The default is HeightScale -> 4."
LeftOrigin::usage = "LeftOrigin is an option of DrawBeam. It specifies the
coordinates of the left end of the beam axis. The default is {10, 50}."
AxisLine::usage = "AxisLine is an option of DrawBeam. Axisline -> True draws
a dashed line along the beam axis. The default is Axisline -> False."
GrooveRatio::usage = "GrooveRatio is an option of DrawBeam. This sets the
length of the base of the triangular notch that is drawn at moment and
shear connections equal to 2 HeightScale/GrooveRatio. The default
setting is GrooveRatio -> 4."
BeamLine::usage = "BeamLine is an option of DrawBeam. It sets the
AbsoluteThickness of the lines in the beam outline. The default is 0.75."
ShearRadius::usage = " ShearRadius is an option of DrawBeam. It specifies
the radius of the filled disk in the symbol for a shear connection. The
default is 0.5."
MomentBase::usage = "MomentBase is an option of DrawBeam. It sets the width
of the rectangle drawn to represent a moment connection. The default is
1.25."
ForceHeight::usage = "ForceHeight is an option of DrawBeam. It sets the
height of the vertical arrows representing concentrated forces. The default
is 16."
ForceTagLoc::usage = "ForceTagLoc is an option of DrawBeam. It sets the
vertical distance between a concentrated or distributed load and its
label. The default is 1."
ForceTipHeight::usage = "ForceTipHeight is an option of DrawBeam.
ForceTipHeight -> h sets the base of the filled triangle, or the tip, of
the concentrated force arrows equal to 2 h. The default is 0.5."
ForceTipLength::usage = "ForceTipLength is an option of DrawBeam. It
specifies the height of the filled triangles, or the tips, of the arrows
representing concentrated forces. The default is 1.6."
ForceThickness::usage = "ForceThickness is an option of DrawBeam. The
AbsoluteThickness of the line between the origin and the terminus of
concentrated force arrows. The default is 0.75."
LoadHeight::usage = "LoadHeight is an option of DrawBeam. It sets the
height of the longest arrow among the vertical arrows representing
distributed forces. The default is 12."
LoadTipHeight::usage = " LoadTipHeight is an option of DrawBeam.
LoadTipHeight -> h sets the base of the filled triangle, or the tip, of the
distributed force arrows equal to 2 h The default is 0.5."
LoadTipLength::usage = "LoadTipLength is an option of DrawBeam. It sets the
height of the filled triangles, or the tips, of the vertical arrows
representing distributed forces. The default is 1.6."
LoadThickness::usage = "LoadThickness is an option of DrawBeam. The
AbsoluteThickness of the line between the origin and the terminus of
distributed force arrows. The default is 0.75."
LoadTags::usage = "LoadTags is an option of DrawBeam. This is a list of
strings specifying labels for the loads defined by LineLoad. If LoadTags
is not given as an option, the load intensities specified by LineLoad are
used as labels. For constant loads, only one label is given, which is
centered horizontally at the midpoint of the load range. For nonconstant
linear loads, two labels are given, to be placed at the start and at the
end of the loading curve. If, however, the load function is zero at one
end, no label can be specified for that end."
LengthUnit::usage = "LengthUnit is an option of DrawBeam. It gives the
string to be used as a length unit in labeling distributed loads,
concentrated couples, and locations marked by the scaling line. The default
is LengthUnit -> \"m\". Distributed loads are labeledForceUnit/LengthUnit,
concentrated couples ForceUnit LengthUnit."
ForceUnit::usage = "ForceUnit is an option of DrawBeam. It gives the string
to be used as a force unit in labeling concentrated loads, distributed
loads, and concentrated couples. The default is ForceUnit -> \"kN\". "
ScaleLevel::usage = "ScaleLevel is an option of DrawBeam. If s is
positive, ScaleLevel -> s, places the scaling line at a vertical
distance s below the lower beam line. If s is negative, the scaling
line is at a vertical distance s above the upper beam line. The default
setting is ScaleLevel -> 14."
ScaleTagLoc::usage = "ScaleTagLoc is an option of DrawBeam. If s is
positive, ScaleTagLoc -> s, places the scaling line labels at a vertical
distance s above the double-tipped arrows. If s is negative, the labels are
placed below the arrows. The default setting is ScaleTagLoc -> 1.25"
ScaleTags::usage = "ScaleTags is an option of DrawBeam. DrawBeam places a
length scaling set of double-tipped horizontal arrows above or below the
drawing of the beam. These arrows mark the location of flexural rigidity
changes, in-span supports, concentrated forces and moments, the start and
end points of distributed loads, shear and moment connections, and the two
ends of the beam. In a numerical problem, the values of the lengths are
used to label the arrows. The LengthUnit option can be used to append units
to these numerical labels. The rule ScaleTags -> {\"s1\", \"s2\", ...}
labels the first arrow s1, the second arrow s2 etc., starting at the left
end and proceeding to the right along the beam axis. The number of strings
in the label list must be the same as the number of double-tipped arrows in
the scaling line."
ScaleNear::usage = "ScaleNear is an option of DrawBeam. Vertical lines
passing through the tips of the double-tipped arrows of the scaling line -
see ScaleTags - are drawn to mark length dimensions. The distance between
the lower beam line and the upper end of these vertical lines is defined by
ScaleNear. If ScaleLevel is negative, which indicates that the scaling
line is to be drawn above the beam, then ScaleNear specifies the distance
between the upper beam line and the lower end of the vertical lines. The
default setting is ScaleNear -> 1."
ScaleFar::usage = "ScaleFar is an option of DrawBeam. This gives the
distance between the lower or upper beam line and the far end of the
vertical lines that mark the scaling line (see ScaleNear). The
default setting is ScaleFar -> 15."
ArrowThickness::usage = "ArrowThickness is an option of DrawBeam. It
specifies the AbsoluteThickness of the double-tipped arrows drawn to mark
lengths. The default is 0.25."
TipHeight::usage = "TipHeight is an option of DrawBeam. It
specifies the width of the base of the tip of double-tipped arrows drawn to
mark lengths. The default is 0.30"
TipLength::usage = "TipLength is an option of DrawBeam. It
specifies the length of the tip of double-tipped arrows drawn to
mark lengths. The default is 1.6."
Xmax::usage = "Xmax is an option of DrawBeam, plots, and plotm. It assigns
a maximum width, in the units used in these program, to the display area.
The default is 100."
Ymax::usage = "Ymax is an option of DrawBeam, plots, and plotm. It assigns
a maximum height, in the units used in these program, to the display area.
The default is 130."
XTitle::usage = "XTitle is an option of plots and plotm. It assigns the
x-axis title."
Y1Title::usage = "Y1Title is an option of plots and plotm. It assigns the
y-axis title for the plot on the upper half of the page."
Y2Title::usage = "Y2Title is an option of plots and plotm. It assigns the
y-axis title for the plot on the lower half of the page."
TopTitle1::usage = "TopTitle1 is an option of plots and plotm. It assigns the
label for the top edge of the upper plot."
TopTitle2::usage = "TopTitle2 is an option of plots and plotm. It assigns the
label for the top edge of the lower plot."
Units::usage = "Units is an option of plots and plotm. Units -> {\"m\",
\"mm\", \"mm/m\"} specifies that x is measured in meters, the deflection in
millimeters, and the slope in millimeters per meter. These are the default
values in plots, which are obtained if lengths x, L are in m, forces in kN,
couples in kNm, load intensities in kN/m, and flexural rigidities in MNm^2
in the input data to Beam. In plotm, the default is Units -> {\"m\",
\"kNm\", \"kN\"}}. To supress units use Units -> None."
SideTitle::usage = "SideTitle is an option of plots and plotm. It assigns
a label for the right edge of the plots."
PS::usage = "PS is an option of DrawBeam, plots, and plotm. It gives a
PostScript file name to save a plot. The default, PS -> None, causes
plots to be displayed on the screen only."
Spacing::usage = "Spacing is an option of Spacing -> h leaves a vertical
space of height equal to 100 h % of the maximum height of the display area
between the slope and deflection plots. The default is h = 0.10."
WallDivisions::usage = "WallDivisions is an option of DrawBeam. It
specifies the number of subdivisions of the line indicating a fixed
vertical wall in a fixed support. This determines the number of short lines
drawn at the wall. The default is 8."
TickThickness::usage = "TickThickness is an option of
DrawBeam. TickThickness -> t sets the AbsoluteThickness of the short lines
in a fixed wall symbol equal to WallThickness/t. The setting affects the
walls for fixed, guided, and pinned supports. The default is 3."
WallThickness::usage = "WallThickness is an option of DrawBeam. The
AbsoluteThickness of the wall line in fixed, guided, and pinned
supports. The default is 1.5."
RectThickness::usage = "RectThickness is an option of DrawBeam. It sets the
AbsoluteThickness of the lines forming the rectangles that represent
moment connections. The default is 0.75."
RectTheta::usage = "RectTheta is an option of DrawBeam. This specifies, in
degrees, the angle by which the rectangles that represent moment
connections are to be rotated counterclockwise from the positive
x-axis. The default is 0."
PinHeight::usage = "PinHeight is an option of DrawBeam. It specifies the
height of the triangle drawn to represent a pinned support. The default is
2.5."
PinWidth::usage = "PinWidth is an option of DrawBeam. It sets the half base
length of the triangle drawn to represent a pinned support. The default is
1."
BaseLength::usage = "BaseLength is an option of DrawBeam. It sets the
length of the fixed horizontal wall line drawn at a pinned support. The
default is 3."
PinLine::usage = "PinLine is an option of DrawBeam. It sets the
AbsoluteThickness of the lines for the triangle drawn to represent a pinned
support. The default is 0.75."
BaseDivisions::usage= "BaseDivisions is an option of DrawBeam. It specifies
the number of subdivisions of the line whose length is specified by the
BaseLength option. Short lines drawn at the subdivisions indicate
the fixed horizontal wall of a pinned support. The default is 6."
ArcRadius::usage = "ArcRadius is an option of DrawBeam. It is the radius of
the circular arrow drawn to represent a concentrated moment. The default is
2."
ArcTipHeight::usage = "ArcTipHeight is an option of DrawBeam. ArcTipHeight
-> h sets the base of the filled triangle, or the tip, of the concentrated
moment arrows equal to 2 h. The default is 0.5."
ArcTipLength::usage = "ArcTipLength is an option of DrawBeam. It is the
height of the filled triangles, or the tips, of the circular arrows
representing concentrated moments. The default is 1.6."
ArcThickness::usage = "ArcThickness is an option of DrawBeam. It sets the
AbsoluteThickness of the circular arc drawn to indicate a concentrated
couple. The default is 0.75."
ArcTagLoc::usage = "ArcTagLoc is an option of DrawBeam. ArcTagLoc -> {a, b}
centers the label for the concentrated moments at the coordinates (a, b)
measured from the center of the circular arc. The default is ArcTagLoc ->
{0, -3}."
ArcAngles::usage = "ArcAngles is an option of DrawBeam. ArcAngles -> {a, b}
sets the beginning and ending angles, in degrees, of the circular arc drawn
to represent a concentrated moment. The default setting is ArcAngles ->
{100, 270}"
RollerRadius::usage = "RollerRadius is an option of DrawBeam. It is the
radius of the circles drawn as part of a guided support. The default is
0.5."
GuideWall::usage = "GuideWall is an option of DrawBeam. It sets the number
of subdivisions of the line indicating a fixed vertical wall in a guided
support. This determines the number of short lines drawn at the wall. The
default is 6."
RollerLine::usage = "RollerLine is an option of DrawBeam. It sets the
AbsoluteThickness of the two small circles in a guided support. The default
is 0.75."
Begin["`Private`"]
Fc[{P_, a_}][x_] := -P UnitStep[-a + x]
Fd[{p_, a_, q_, b_}][x_] := -{(-a + x) (-(a p) + 2 b p - a q - p x + q x) UnitStep[-a + x],
(b - x) (b p - 2 a q + b q - p x + q x) UnitStep[-b + x]}/(2 (-a + b))
Fdist[{p_, a_, q_, b_}][x_] := -((b p - a q - p x + q x) (UnitStep[-a + x] - UnitStep[-b + x]))/(-a + b)
Cp[{P_, a_}][x_] := -P DiracDelta[x - a]
Sc[{dtheta_, a_}][x_] := -dtheta DiracDelta'[x - a]
Mc[{dv_, a_}][x_] := -dv DiracDelta''[x - a]
Ef[{e_, a_}][x_] := e UnitStep[x - a]
urule[x_, L_] := {UnitStep[(a_)?(NumberQ[N[#/L]] &)] :> UnitStep[Cancel[a/L]], UnitStep[x] :> 1}
uregion[x_, L_, expr_] :=
Block[{steps, len, stepf},
stepf = expr /. urule[x, L];
steps = Union[Cases[stepf, _UnitStep, Infinity]];
len = Length[steps];
If[len == 0, Return[{{expr}, {{0, L}}}]];
steps = Sort[ Apply[Sequence, steps, 1], OrderedQ[{x - #1 , x - #2}] &];
stepf = Simplify[stepf /. (Thread[(UnitStep /@ steps) -> #]& /@ Prepend[Array[(#1 - #2)&, {len, len}] /.
{_?NonNegative ->1, _?Negative -> 0}, Array[0&, len]])];
{stepf, Append[Partition[Prepend[x - steps, 0], 2, 1], {x - Last[steps], L}]}
]
addsymbol[list_] := Transpose[{Array[Unique[]&, Length[list]], list}]
evalrep[x_, v_, vp_] := ReplaceAll[v, #] & /@ Thread[x -> Last /@ vp]
Options[Beam] = {PointLoad -> {}, LineLoad -> {}, InSpan -> {}, PointCouple -> {},
LoadFunction -> {}, ShearConnection -> {}, MomentConnection -> {}}
Beam[x_Symbol, L_, EI_, bnd:{_String, _String}, options___Rule] :=
Block[{pointload, lineload, inspan, couple, loadfunc, shearcon, momentcon,
flex, loc, V, M, S, U, c1, c2, c3, c4, nhL, nhR, bndL, bndR, inscmc, eqns, soln, res},
{pointload, lineload, inspan, couple, loadfunc, shearcon, momentcon}
= First /@ Options[Beam] /. parameters[Beam, options];
{inspan, shearcon, momentcon} = addsymbol /@ Flatten /@ {{inspan}, {shearcon}, {momentcon}};
{flex, loc} = Transpose[Partition[If[Head[EI] === List, EI, {EI, 0}], 2]];
flex = Prepend[ Apply[ Ef[{1/#2 - 1/#1, #3}][x] &,
MapThread[Append, {Partition[flex, 2,1], Rest[loc]}], 1], 1/First[flex]] /. urule[x, L];
V = Join[
Map[ Fc[#1][x] &, DeleteCases[Join[Partition[pointload, 2], inspan], {_, 0 | 0.0} | {_, L}] ],
Flatten[Map[ Fd[#1][x] &, Partition[lineload, 4] ]] /. UnitStep[x - L] :> 0,
-Integrate[ Flatten[{loadfunc}], x ],
Map[ Sc[#1][x] &, shearcon ],
Map[ Mc[#1][x] &, momentcon ],
Map[ Cp[#1][x] &, DeleteCases[ Partition[couple, 2], {_, 0 | 0.0} | {_, L} ]],
{c1}
];
M = Append[Integrate[V, x], c2];
(*
S = Append[SimplifyUnitStep[ (Integrate[ Flatten[Map[Times[flex, #1]&, M]], x ] /. urule[x, L])
/. UnitStep[x + a_] -> UnitStep[x + a/L] ] /. UnitStep[x + a_] -> UnitStep[x + a L], c3];
*)
S = Append[Simplify[ (Integrate[ Flatten[Map[Times[flex, #1]&, M]], x ] /. urule[x, L])
/. UnitStep[x + a_] -> UnitStep[x + a/L] ] /. UnitStep[x + a_] -> UnitStep[x + a L], c3];
U = Append[Integrate[S, x], c4];
nhL = {Cases[Partition[pointload, 2], {_, 0 | 0.0}], Cases[Partition[couple, 2], {_, 0 | 0.0}]} /.
{{{P_, _}} :> P, {} :> 0};
bndL = Switch[ First[bnd], "Fixed", {U, S}, "Pinned", {U, Append[M, Last[nhL]]},
"Free", {Append[M, Last[nhL]], Append[V, First[nhL]]}, "Guided", {S, Append[V, First[nhL]]}] /. x -> 0;
nhR = {Cases[Partition[pointload, 2], {_, L}], Cases[Partition[couple, 2], {_, L}]} /.
{{{P_, _}} :> P, {} :> 0};
bndR = Switch[ Last[bnd], "Fixed", {U, S}, "Pinned", {U, Append[M, -Last[nhR]]},
"Free", {Append[M, -Last[nhR]], Append[V, -First[nhR]]}, "Guided", {S, Append[V, -First[nhR]]}] /. x -> L;
inscmc = Join[bndR, bndL, evalrep[x, U, inspan], evalrep[x, M, shearcon], evalrep[x, V, momentcon]];
eqns = Thread[ Apply[Plus, inscmc, 1] == 0 ] /. urule[x, L] /. NoDelta;
soln = Flatten[ Solve[eqns, Join[ {c1, c2, c3, c4}, First /@ Join[inspan, shearcon, momentcon]]] ];
res = Apply[Plus, {V, M, S, U} /. NoDelta /. soln /. urule[x, L], 1];
DeleteCases[Append[Simplify[ Collect[res, Cases[res, _UnitStep, Infinity]] ], inspan /. soln], {}]
]
Options[plotf] = {XTitle -> "", YTitle -> "", TopTitle -> "", SideTitle -> ""}
plotf[f_, x_, x0_, xf_, options___Rule] :=
Block[ {htitle, vtitle, toptitle, stitle},
{htitle, vtitle, toptitle, stitle} = First /@ Options[plotf] /. parameters[plotf, options];
Plot[
f, {x, x0, xf},
Evaluate[ FilterOptions[Plot, options] ],
Frame -> True,
GridLines -> None,
PlotRange -> All,
PlotPoints -> 40,
DisplayFunction -> Identity,
FrameLabel -> {htitle, vtitle, toptitle, stitle},
PlotStyle -> {AbsoluteThickness[0.25]},
DefaultFont -> {"Times-Italic", 10},
FrameStyle -> {AbsoluteThickness[0.125]},
Axes -> {True, False},
AxesStyle -> {AbsoluteThickness[0.125]}
]
]
Options[Draw] = {PS -> None}
Draw[graph_, options___Rule] :=
Block[{G, psfile},
psfile = PS /. parameters[Draw, options];
G = Graphics[graph,
FilterOptions[Graphics, options], AspectRatio -> Automatic, PlotRange -> All,
DefaultFont -> {"Times-Italic", 10}];
(* If[psfile === None, Show[G], Display[ StringJoin["!/home/lk3a/mma/psfix > ", psfile], G]] *)
If[psfile === None, Show[G], Display[ StringJoin["!psfix > ", psfile], G]]
]
parameters[name_, options___Rule] := Join[ {FilterOptions[name, options]}, Options[name] ];
NoDelta := DiracDelta -> (0 &)
Options[plots] = {XTitle -> "x", Y1Title -> "DEFLECTION", Y2Title-> "SLOPE",
SideTitle -> " ", TopTitle1 -> " ", TopTitle2 -> " ",
Units -> {"m", "mm", "mm/m"}, Xmax -> 100, Ymax -> 130, Spacing -> 0.1, PS -> None}
plots[f1_, f2_, x_Symbol, x0_, xf_, options___Rule] :=
Block[{xtitle, y1title, y2title, sidetitle, toptitle1,
toptitle2, units, xmax, ymax, spacing, oparen, cparen},
{xtitle, y1title, y2title, sidetitle, toptitle1, toptitle2, units, xmax, ymax, spacing} =
Drop[First /@ Options[plots], -1] /. parameters[plots, options];
If[units === None, oparen = ""; cparen = ""; units = {"", "", ""}, oparen = " ("; cparen = ")"];
Draw[ {Rectangle[{0, (1 + spacing) ymax/2}, {xmax, ymax},
plotf[f1, x, x0, xf, XTitle -> StringJoin[xtitle, oparen, First[units], cparen],
YTitle -> StringJoin[y1title, oparen, Part[units, 2], cparen],
TopTitle -> toptitle1, SideTitle -> sidetitle, options]],
Rectangle[{0, 0}, {xmax, (1 - spacing) ymax/2},
plotf[f2, x, x0, xf, XTitle -> StringJoin[xtitle, oparen, First[units], cparen],
YTitle -> StringJoin[y2title, oparen, Last[units], cparen],
TopTitle -> toptitle2, SideTitle -> sidetitle, options]]}, options]
]
plotm[f1_, f2_, x_Symbol, x0_, xf_, options___Rule] :=
plots[f1, f2, x, x0, xf, options, Y1Title -> "BENDING MOMENT", Y2Title-> "SHEAR FORCE", Units -> {"m", "kNm", "kN"}]
Options[wall] = {WallDivisions -> 8, TickThickness -> 3, WallThickness -> 1.5, WallTheta -> 0, TickDirection -> -1}
wall[x0_, y0_, L_, options___Rule] :=
Block[{nsubdiv, t, tk, tt, theta, pts1, pts2},
{nsubdiv, t, tk, tt, theta} = {WallDivisions, WallThickness, TickDirection, WallThickness/TickThickness, WallTheta}
/. parameters[wall, options];
pts1 = L/nsubdiv Array[{#1, 0}&, nsubdiv + 1, 0];
pts2 = Map[Dot[T[theta], #1] &, Transpose[{pts1, Map[(L/nsubdiv {tk, -1} + #1)&, pts1]}], {2}];
{{AbsoluteThickness[tt], Line /@ Map[(# + {x0, y0})&, Drop[pts2, -tk], {2}]},
{AbsoluteThickness[t], Line[{{x0, y0}, {x0, y0} + T[theta].Last[pts1]}]}}
]
leftwall[lowx_, lowy_, L_, options___Rule] := wall[lowx, lowy + L, L, WallTheta -> -90, TickDirection -> 1, options]
rightwall[lowx_, lowy_, L_, options___Rule] := wall[lowx, lowy, L, WallTheta -> 90, TickDirection -> -1, options]
topwall[leftx_, lefty_, L_, options___Rule] := wall[leftx + L, lefty, L, WallTheta -> 180, TickDirection -> 1, options]
bottomwall[leftx_, lefty_, L_, options___Rule] := wall[leftx, lefty, L, WallTheta -> 0, TickDirection -> -1, options]
points[vectors_List, origin_] := FoldList[ Plus, origin, vectors ]
T[theta_] := With[{x = N[theta Degree]}, {{Cos[x], -Sin[x]},{Sin[x], Cos[x]}}]
T[theta_, vectors_List] := Map[Dot[T[theta], #1]&, vectors]
Options[arrow] = {TipTag -> " ", TipTagLoc -> {0, 0}, TipHeight -> 0.5, TipLength -> 1.6,
ArrowTheta -> 0, ArrowThickness -> 0.75}
arrow[x0_, y0_, L_, options___Rule] :=
Block[{theta, u, plist, length, height, thk, text, textloc},
{text, textloc, height, length, theta, thk} = First /@ Options[arrow] /. parameters[arrow, options];
If[ 0 <= L < length, Return[{}]];
u = {Cos[theta N[Degree]], Sin[theta N[Degree]]};
{height, length} = {height {-Last[u], First[u]}, length u};
plist = points[{u If[L < 0, 0, L] - length, height, -height + length, -height - length, height }, {x0, y0}];
{AbsoluteThickness[thk], Line[Part[plist, {1, 2}]], Polygon[Rest[plist]], Text[text, {x0, y0} + T[theta] . textloc]}
]
doublearrow[x0_, y0_, L_, options___Rule] :=
Block[ {theta, text, textloc},
{theta, text, textloc} = {ArrowTheta, TipTag, TipTagLoc} /. parameters[arrow, options];
{arrow[x0, y0, L, TipTag -> "", options],
arrow[x0, y0, -1, ArrowTheta -> 180 + theta, TipTag -> "", options],
Text[text, {x0, y0} + T[theta] . textloc]}
]
Options[rectangle] = {RectTag -> " ", RectTagLoc -> {0, 0}, RectHeight -> 1, RectThickness -> 0.75, RectTheta -> 0}
rectangle[x0_, y0_, L_, options___Rule] :=
Block[{h, thk, plist, theta, text, textloc},
{h, thk, theta, text, textloc} = {L/RectHeight, RectThickness, RectTheta, RectTag, RectTagLoc}
/. parameters[rectangle, options];
plist = points[ T[theta, {{-L/2, h/2}, {L, 0}, {0, -h}, {-L, 0}, {0, h}}], {x0, y0} ];
{AbsoluteThickness[thk], Line[Rest[plist]], Text[text, {x0, y0} + textloc]}
]
Options[pinsupport] = {PinHeight -> 2.5, PinWidth -> 1, BaseLength -> 3, PinLine -> 0.75, BaseDivisions -> 6}
pinsupport[x0_, y0_, options___Rule] :=
Block[{h, w, L, thk, wd},
{h, w, L, thk, wd} = First /@ Options[pinsupport] /. parameters[pinsupport, options];
{{AbsoluteThickness[thk], Line[{{x0 - w, y0 - h}, {x0, y0}, {x0 + w, y0 - h}}]},
bottomwall[x0 - L/2, y0 - h, L, WallDivisions -> wd, options]}
]
Options[guide] = {RollerRadius -> 0.5, GuideWall -> 6, InterRadial -> 0, RollerLine -> 0.75, LeftRight -> 1}
guide[xl_, yl_, L_, options___Rule] :=
Block[{r, wd, h, thk, s},
{r, wd, h, thk, s} = First /@ Options[guide] /. parameters[guide, options];
{{AbsoluteThickness[thk], Circle[{xl + r s, yl + (L - h)/2}, r], Circle[{xl + r s, yl + (L + h)/2}, r]},
If[ s > 0, leftwall, rightwall][xl, yl, L, WallDivisions -> wd, options]}
]
Options[arc] = {ArcSense -> 1, ArcRadius -> 2, ArcTipHeight -> 0.5, ArcTipLength -> 1.6,
ArcThickness -> 0.75, ArcTag -> "", ArcTagLoc -> {0, 0}, ArcAngles -> {100, 270}}
arc[x0_, y0_, options___Rule] :=
Block[{s, r, th, tl, thk, theta1, theta2, arrowangle, tag, tagloc},
{s, r, th, tl, thk, tag, tagloc, {theta1, theta2}} = First /@ Options[arc] /. parameters[arc, options];
{theta1, theta2} = N[Degree] If[ s == 1, {theta1, theta2}, {theta2, theta1}];
arrowangle = theta2 + If[s == 1, Pi/2, -Pi/2];
{{AbsoluteThickness[thk], Circle[{x0, y0}, r, Sort[{theta1, theta2}]],
Line[{{x0, y0 - r/8}, {x0, y0 + r/8}}], Line[{{x0 - r/8, y0}, {x0 + r/8, y0}}]},
Text[tag, {x0, y0} + tagloc],
arrow[x0 + r Cos[theta2] + tl Cos[arrowangle], y0 + r Sin[theta2] + tl Sin[arrowangle], -1,
TipHeight -> th, TipLength -> tl, ArrowThickness -> thk, ArrowTheta -> arrowangle/N[Degree]]}
]
midpoint[list_] := list /. {a_?NumberQ, b_?NumberQ} -> (a + b)/2
midpts[list_, func_, x_] := Map[(Evaluate[func] & /. x->#), midpoint[list] ]
maxrep[x_List] := Array[Max[x]&, Length[x]]
signfilter[x_?NonNegative ] := {x, 1}
signfilter[x_?Negative] := {-x, 0}
signfilter[x_] := With[{sym = First[Variables[x]]}, signfilter[x /. sym -> 1] {sym, 1} ]
sortfirst := (First[#1] <= First[#2]) &
Options[DrawBeam] = {PointLoad -> {}, LineLoad -> {}, InSpan -> {}, PointCouple -> {},
LoadFunction -> {}, LoadFunctionRange -> {}, LoadFunctionTags -> {},
LoadFunctionTagLoc -> 2.5, LoadFunctionSubdivisions -> 5,
ShearConnection -> {}, MomentConnection -> {},
LengthScale -> 80, HeightScale -> 4, LeftOrigin -> {10, 50}, GrooveRatio -> 4,
BeamLine -> 0.75, ShearRadius -> 0.5, MomentBase -> 1.25,
ForceHeight -> 16, ForceTagLoc -> 1, ForceTipHeight -> 0.5, ForceTipLength -> 1.6,
ForceThickness -> 0.75, LoadHeight -> 12, LoadTipHeight -> 0.5,
LoadTipLength -> 1.6, LoadThickness -> 0.75, LoadSubdivisions -> 5,
LoadTags -> {}, ArcTagLoc -> {0, -3},
LengthUnit -> "m", ForceUnit -> " kN", ScaleLevel -> 14, ArrowThickness -> 0.25,
TipLength -> 1.6, TipHeight -> 0.3,
ScaleTagLoc -> 1.25, ScaleTags -> {}, ScaleNear -> 1, ScaleFar -> 15,
Xmax -> 100, Ymax -> 130, AxisLine -> False}
DrawBeam[x_Symbol, L_, EI_, bnd:{_String, _String}, options___Rule] :=
Block[
{alengths, arcrad, atloc, axisline, beamfig, bndL, bndR, bthk, ccargs,
ccloc, ccouples, cfargs, cfloc,
cforces, cfsigns, cons, couple, cpoints, dathk, datph, datplen,
distargs, ei, eirange, eix, fh, fth, fthk, ftl, ftloc, fu, gr,
graph, hpoints, hs, insp, lfns, lfrange, lftag, lftags, lftloc, lh,
lldiv, llend, lload, lloadfns, llranges, lltag, loadf,
loc, lower, lrend, ls, lsubdiv, ltags, lth, lthk, ltl, lu, maxup, mlen,
momcon, nsubdiv, ptload, ptsfns, ptsll, radius,
rotations, scaleinsp, scfar, scnear, shcon, shradius, silh, slevel,
stloc, ticks, ticktagloc, ticktags, ulend, upper, urend,
vscale, xL, xloc, xmax, xR, xrange, x0, ymax, y0},
{ptload, lload, insp, couple, loadf, lfrange, lftags, lftloc, nsubdiv,
shcon, momcon, ls, hs, {x0, y0}, gr,
bthk, shradius, mlen, fh, ftloc, fth, ftl, fthk, lh, lth, ltl, lthk,
lsubdiv, ltags, atloc,
lu, fu, slevel, dathk, datplen, datph, stloc, ticktags, scnear,
scfar, xmax, ymax, axisline} =
First /@ Options[DrawBeam] /. parameters[DrawBeam, options];
shcon = Flatten[{shcon}];
momcon = Flatten[{momcon}];
radius = RollerRadius /. parameters[guide, options];
arcrad = ArcRadius /. parameters[arc, options];
{ei, loc} = If[Head[EI] === List, Transpose[Partition[ EI, 2]], {{EI}, {0}}];
ei = ei hs /(2 Max[ei]);
loc = loc ls/L;
ei = Plus @@ Map[ Ef[#1][x] &, Transpose[{ei - Drop[ Insert[ei, 0, 1], -1], loc}] ];
upper = Evaluate[(y0 + ei)] & /. x -> #1;
lower = Evaluate[(y0 - ei)] & /. x -> #1;
cons = Join[shcon, momcon] ls/L;
hpoints = (upper /@ cons - lower /@ cons)/gr;
hpoints = Join[cons - hpoints, cons + hpoints];
eix = Sort[Join[Append[loc, ls], Rest[loc], hpoints]] + x0;
eirange = Sort[Join[Append[loc, ls], 0.98 Rest[loc], hpoints]];
cpoints = Map[{# + x0, y0}&, cons];
beamfig = Join[ Sort[ Join[Transpose[{eix, upper /@ eirange}], cpoints], sortfirst ],
Reverse[ Sort[ Join[Transpose[{eix, lower /@ eirange}], cpoints], sortfirst ] ]];
xL = First[eix];
xR = Last[eix];
llend = lower[0];
ulend = upper[0];
lrend = lower[ls];
urend = upper[ls];
bndL = Switch[ First[bnd],
"Free", {},
"Pinned", pinsupport[xL, llend, options],
"Fixed", leftwall[xL, (3 llend - ulend)/2, 2(ulend-llend), options],
"Guided", guide[xL - 2 radius, (3 llend - ulend)/2, 2(ulend-llend), options, InterRadial -> 2(ulend-llend)/3]];
bndR = Switch[ Last[bnd],
"Free", {},
"Pinned", pinsupport[xR, lrend, options],
"Fixed", rightwall[xR, (3 lrend - urend)/2, 2(urend-lrend), options],
"Guided", guide[xR + 2 radius, (3 lrend - urend)/2, 2(urend-lrend), options,
LeftRight -> -1, InterRadial -> 2(urend-lrend)/3]];
scaleinsp = Flatten[{insp}] ls/L;
cfargs = If[ptload === {}, {ptload},
{cforces, cfloc} = Transpose[Partition[ptload, 2]];
{cforces, cfsigns} = Transpose[signfilter /@ cforces];
cfloc = cfloc ls/L;
{x0 + cfloc, fh cfsigns + upper /@ cfloc,
StringJoin[#1, fu] & /@ (ToString /@ InputForm /@ cforces),
(1 - cfsigns) (fh + 2 ftloc) - ftloc, 90 - 180 cfsigns}];
ccargs = If[couple === {}, {couple},
{ccouples, ccloc} = Transpose[Partition[couple, 2]];
{ccouples, ccsigns} = Transpose[signfilter /@ ccouples];
ccloc = ccloc ls/L;
{x0 + ccloc, StringJoin[#1, fu, lu] & /@ (ToString /@ InputForm /@ ccouples), ccsigns}];
lloadfns = Plus @@ Flatten[Map[ Fdist[#1][x] &, Partition[lload, 4] ]];
lfns = Plus @@ Flatten[{-loadf}];
lltag = {};
lftag = {};
{distargs, silh} = If[
lloadfns === 0 && lfns === 0, {{{}},{}},
If[ lfns === 0, xrange = {}; ptsfns = {},
xrange = If[ lfrange === {}, {L/nsubdiv Range[0, nsubdiv]},
llranges = lfrange /. {a_?NumberQ, b_?NumberQ} :> {a, 0.9999 b};
lldiv = llranges /. {a_?NumberQ, b_?NumberQ} :> b - a;
lldiv = Partition[Ceiling /@ (lldiv/ (Max[lldiv]/nsubdiv)), 1];
Apply[(#1 + (Range[0, #3]/#3)(#2 - #1))&, Apply[Join, Transpose[{llranges, lldiv}], 1], 1]
];
ptsfns = Flatten[lfns /. Partition[Thread[x -> Flatten[xrange]], 1]]
];
If[ lloadfns === 0, ptsll = {}; llranges = {}; lldiv = nsubdiv; xloc = {}; lltag = {},
ptsll = First /@ Partition[lload, 2];
lltag = Partition[lload, 4] /. {a_?NumberQ, b_?NumberQ, a_?NumberQ, c_?NumberQ} :> {a, (b + c)/2};
llranges = Partition[Last /@ Partition[lload, 2], 2] /. {a_?NumberQ, b_?NumberQ} :> {a, 0.9999 b};
lldiv = llranges /. {a_?NumberQ, b_?NumberQ} :> b - a;
lldiv = Partition[Ceiling /@ (lldiv/ (Max[lldiv]/lsubdiv)), 1];
xloc = Apply[(#1 + (Range[0, #3]/#3)(#2 - #1))&, Apply[Join, Transpose[{llranges, lldiv}], 1], 1]
];
vscale = lh/Ceiling[ Max[Abs[ptsll], Abs[ptsfns]] ];
alengths = -vscale Join[Map[ (ReplaceAll[ lloadfns, #]) &, Map[Thread[x -> #]&, xloc], {2}],
Map[ (ReplaceAll[ lfns, #]) &, Map[Thread[x -> #]&, xrange], {2}]];
lltag = DeleteCases[Partition[Flatten[lltag], 2], {0, _} | {0.0, _}] /.
{a_?NumberQ, b_?NumberQ} :> {StringJoin[ToString[InputForm[a]], fu, "/", lu],
ls b/L + x0, upper[ls b/L] + a vscale + Sign[a] ftloc};
If[ltags =!= {}, lltag = Transpose[Prepend[Take[Transpose[lltag], -2], ltags]]];
If[lftags =!= {}, lftag = Transpose[{lftags, x0 + ls/L midpoint[lfrange],
- Sign[midpts[lfrange, lfns, x]] lftloc - vscale midpts[lfrange, lfns, x] + upper /@ (ls midpoint[lfrange]/L)}]];
rotations = Map[Last, Map[signfilter, alengths, {2}], {2}] /. {0 -> 90, 1 -> -90};
xloc = x0 + Join[xloc, xrange] ls/L;
maxup = alengths + maxrep /@ Map[upper, xloc, {2}];
{{xloc, maxup, Abs[alengths], rotations}, MapThread[Transpose[{#1, #2}]&, {xloc, maxup}]}
];
ticks = x0 + Join[ Append[loc, ls], scaleinsp, cons,(Last /@ Partition[lload, 2])ls/L, Flatten[lfrange] ls/L ];
ticks = Union[Join[First[cfargs], First[ccargs], ticks]];
ticktagloc = midpoint[Partition[ticks, 2, 1]];
If[ticktags === {}, ticktags = StringJoin[ToString[#1], " ", lu]& /@ (Drop[(RotateLeft[ticks, 1] - ticks), -1] L/ls)];
graph = {
{AbsoluteThickness[bthk], Line[ Append[beamfig, First[beamfig]] ]},
Disk[#, shradius]& /@ Take[cpoints, Length[shcon]],
Apply[ rectangle[#1, #2, mlen, options, RectHeight -> 2/gr] &, Take[cpoints, -Length[momcon]], 1 ],
Apply[ pinsupport[#1, #2, options] &, Transpose[{scaleinsp + x0, lower /@ scaleinsp}], 1 ],
Apply[ arrow[#1, #2, fh, TipHeight -> fth, TipLength -> ftl, ArrowThickness -> fthk,
TipTag -> #3, TipTagLoc -> {#4, 0}, ArrowTheta -> #5]&,
Transpose[cfargs], 1 ],
Apply[ arc[#1, y0, options, ArcTag -> #2, ArcSense -> #3, ArcTagLoc -> atloc] &,
Transpose[ccargs], 1],
Apply[ arrow[#1, #2, #3, ArrowTheta -> #4, TipHeight -> lth, TipLength -> ltl, ArrowThickness -> lthk] &,
Transpose[Flatten /@ distargs], 1],
{AbsoluteThickness[lthk], Line /@ silh},
Apply[ (Text[#1, {##2}] &), lltag, 1 ],
Apply[ (Text[#1, {##2}] &), lftag, 1 ],
Apply[ doublearrow[#1, y0 - slevel - Sign[slevel] hs/2, #2 - #1,
ArrowThickness -> dathk, TipHeight -> datph, TipLength -> datplen]&, Partition[ticks, 2, 1], 1 ],
{AbsoluteThickness[dathk],
Line[{{#1, y0 - Sign[slevel] (hs/2 + scnear)}, {#1, y0 - Sign[slevel] (hs/2 + scfar)}}]& /@ ticks},
Apply[Text[#1, {#2, y0 - Sign[slevel] hs/2 - slevel + stloc}]&, Transpose[{ticktags, ticktagloc}], 1],
If[axisline === False, {},
{AbsoluteThickness[dathk], AbsoluteDashing[{10, 4}], Line[{{x0 - 3 ls/80, y0}, {x0 + 83 ls/80, y0}}]}],
bndL, bndR
};
Draw[Rectangle[{0, 0}, {xmax, ymax},
Graphics[graph, FilterOptions[Graphics, options], AspectRatio -> Automatic, PlotRange -> All]
],
options,
DefaultFont -> {"Times-Italic", 12}]
]
End[]
EndPackage[]