code 33303;
procedure FEM HERM SYM(X, Y, N, P, Q, R, F, ORDER, E);
value N, ORDER; integer N, ORDER;
array X, Y, E;
real procedure P, Q, R, F;
begin integer L, N2, V, W;
array A[1:8*(N - 1)], EM[2:3];
real A11, A12, A13, A14, A22, A23, A24, A33, A34, A44,
YA, YB, ZA, ZB,
B1, B2, B3, B4, D1, D2, E1, R1, R2, XL1, XL;
procedure ELEMENTMATVECEVALUATION;
if ORDER=4then
begin real X2, H, H2, H3, P1, P2,
Q1, Q2, R1, R2, F1, F2,
B11, B12, B13, B14, B22, B23, B24, B33, B34, B44,
S11, S12, S13, S14, S22, S23, S24, S33, S34, S44,
M11, M12, M13, M14, M22, M23, M24, M33, M34, M44;
own real P3, Q3, R3, F3;
H:= XL - XL1; H2:= H*H; H3:= H*H2;
X2:= (XL1 + XL)/2;
if L=1then
begin P3:= P(XL1); Q3:= Q(XL1); R3:= R(XL1); F3:= F(XL1)
end;
comment ELEMENT BENDING MATRIX;
P1:= P3; P2:= P(X2); P3:= P(XL);
B11:= 6*(P1 + P3); B12:= 4*P1 + 2*P3;
B13:= - B11; B14:= B11 - B12;
B22:= (4*P1 + P2 + P3)/1.5; B23:= - B12; B24:= B12 - B22;
B33:= B11; B34:= - B14; B44:= B14 - B24;
comment ELEMENT STIFFNESS MATRIX;
Q1:= Q3; Q2:= Q(X2); Q3:= Q(XL);
S11:= 1.5*Q2; S12:= Q2/4; S13:= - S11; S14:= S12;
S24:= Q2/24; S22:= Q1/6 + S24; S23:= - S12;
S33:= S11; S34:= - S12; S44:= S24 + Q3/6;
comment ELEMENT MASS MATRIX;
R1:= R3; R2:= R(X2); R3:= R(XL);
M11:= (R1 + R2)/6; M12:= R2/24; M13:= R2/6; M14:= - M12;
M22:= R2/96; M23:= - M14; M24:= - M22;
M33:= (R2 + R3)/6; M34:= M14; M44:= M22;
comment ELEMENT LOAD VECTOR;
F1:= F3; F2:= F(X2); F3:= F(XL);
B1:= H*(F1 + 2*F2)/6; B3:= H*(F3 + 2*F2)/6;
B2:= H2*F2/12; B4:= - B2;
A11:= B11/H3 + S11/H + M11*H; A12:= B12/H2 + S12 + M12*H2;
A13:= B13/H3 + S13/H + M13*H; A14:= B14/H2 + S14 + M14*H2;
A22:= B22/H + S22*H + M22*H3; A23:= B23/H2 + S23 + M23*H2;
A24:= B24/H + S24*H + M24*H3; A34:= B34/H2 + S34 + M34*H2;
A33:= B33/H3 + S33/H + M33*H; A44:= B44/H + S44*H + M44*H3
end else if ORDER=6then
begin own real P4, Q4, R4, F4;
real H, H2, H3, X2, X3,
P1, P2, P3, Q1, Q2, Q3,
R1, R2, R3, F1, F2, F3,
B11, B12, B13, B14, B15, B22, B23, B24, B25,
B33, B34, B35, B44, B45, B55,
S11, S12, S13, S14, S15, S22, S23, S24, S25,
S33, S34, S35, S44, S45, S55,
M11, M12, M13, M14, M15, M22, M23, M24, M25,
M33, M34, M35, M44, M45, M55,
A15, A25, A35, A45, A55, C1, C2, C3, C4, B5;
if L=1then
begin P4:= P(XL1); Q4:= Q(XL1); R4:= R(XL1); F4:= F(XL1)
end;
H:= XL - XL1; H2:= H*H; H3:= H*H2;
X2:= 0.27639320225*H + XL1; X3:= XL1 + XL - X2;
comment ELEMENT BENDING MATRIX;
P1:= P4; P2:= P(X2); P3:= P(X3); P4:= P(XL);
B11:= + 4.0333333333333"+1*P1 + 1.1124913866738"-1*P2
+ 1.4422084194664"+1*P3 + 8.3333333333333"+0*P4;
B12:= + 1.4666666666667"+1*P1 - 3.3191425091659"-1*P2
+ 2.7985809175818"+0*P3 + 1.6666666666667"+0*P4;
B13:= + 1.8333333333333"+1*(P1+P4)
+ 1.2666666666667"+0*(P2+P3);
B15:= - (B11 + B13); B14:= - (B12 + B13 + B15/2);
B22:= + 5.3333333333333"+0*P1 + 9.9027346441674"-1*P2
+ 5.4305986891624"-1*P3 + 3.3333333333333"-1*P4;
B23:= + 6.6666666666667"+0*P1 - 3.7791278464167"+0*P2
+ 2.4579451308295"-1*P3 + 3.6666666666667"+0*P4;
B25:= - (B12 + B23); B24:= - (B22 + B23 + B25/2);
B33:= + 8.3333333333333"+0*P1 + 1.4422084194666"+1*P2
+ 1.1124913866726"-1*P3 + 4.0333333333333"+1*P4;
B35:= - (B13 + B33); B34:= - (B23 + B33 + B35/2);
B45:= - (B14 + B34); B44:= - (B24 + B34 + B45/2);
B55:= - (B15 + B35);
comment ELEMENT STIFFNESS MATRIX;
Q1:= Q4; Q2:= Q(X2); Q3:= Q(X3); Q4:= Q(XL);
S11:= + 2.8844168389330"+0*Q2 + 2.2249827733448"-2*Q3;
S12:= + 2.5671051872498"-1*Q2 + 3.2894812749994"-3*Q3;
S13:= + 2.5333333333333"-1*(Q2+Q3);
S14:= - 3.7453559925005"-2*Q2 - 2.2546440074988"-2*Q3;
S15:= - (S13 + S11);
S22:= + 8.3333333333333"-2*Q1 + 2.2847006554164"-2*Q2
+ 4.8632677916445"-4*Q3;
S23:= + 2.2546440075002"-2*Q2 + 3.7453559924873"-2*Q3;
S24:= - 3.3333333333333"-3*(Q2+Q3);
S25:= - (S12 + S23);
S33:= + 2.2249827733471"-2*Q2 + 2.8844168389330"+0*Q3;
S34:= - 3.2894812750127"-3*Q2 - 2.5671051872496"-1*Q3;
S35:= - (S13 + S33);
S44:= + 4.8632677916788"-4*Q2
+ 2.2847006554161"-2*Q3 + 8.3333333333338"-2*Q4;
S45:= - (S14 + S34);
S55:= - (S15 + S35);
comment ELEMENT MASS MATRIX;
R1:= R4; R2:= R(X2); R3:= R(X3); R4:= R(XL);
M11:= + 8.3333333333333"-2*R1 + 1.0129076086083"-1*R2
+ 7.3759058058380"-3*R3;
M12:= + 1.3296181273333"-2*R2 + 1.3704853933353"-3*R3;
M13:= - 2.7333333333333"-2*(R2+R3);
M14:= + 5.0786893258335"-3*R2 + 3.5879773408333"-3*R3;
M15:= + 1.3147987115999"-1*R2 - 3.5479871159991"-2*R3;
M22:= + 1.7453559925000"-3*R2 + 2.5464400750059"-4*R3;
M23:= - 3.5879773408336"-3*R2 - 5.0786893258385"-3*R3;
M24:= + 6.6666666666667"-4*(R2+R3);
M25:= + 1.7259029213333"-2*R2 - 6.5923625466719"-3*R3;
M33:= + 7.3759058058380"-3*R2
+ 1.0129076086083"-1*R3 + 8.3333333333333"-2*R4;
M34:= - 1.3704853933333"-3*R2 - 1.3296181273333"-2*R3;
M35:= - 3.5479871159992"-2*R2 + 1.3147987115999"-1*R3;
M44:= + 2.5464400750008"-4*R2 + 1.7453559924997"-3*R3;
M45:= + 6.5923625466656"-3*R2 - 1.7259029213330"-2*R3;
M55:= + .17066666666667"+0*(R2+R3);
comment ELEMENT LOAD VECTOR;
F1:= F4; F2:= F(X2); F3:= F(X3); F4:= F(XL);
B1:= + 8.3333333333333"-2*F1 + 2.0543729868749"-1*F2
- 5.5437298687489"-2*F3;
B2:= + 2.6967233145832"-2*F2 - 1.0300566479175"-2*F3;
B3:= - 5.5437298687489"-2*F2
+ 2.0543729868749"-1*F3 + 8.3333333333333"-2*F4;
B4:= + 1.0300566479165"-2*F2 - 2.6967233145830"-2*F3;
B5:= + 2.6666666666667"-1*(F2+F3);
A11:= H2*(H2*M11 + S11) + B11; A12:= H2*(H2*M12 + S12) + B12;
A13:= H2*(H2*M13 + S13) + B13; A14:= H2*(H2*M14 + S14) + B14;
A15:= H2*(H2*M15 + S15) + B15; A22:= H2*(H2*M22 + S22) + B22;
A23:= H2*(H2*M23 + S23) + B23; A24:= H2*(H2*M24 + S24) + B24;
A25:= H2*(H2*M25 + S25) + B25; A33:= H2*(H2*M33 + S33) + B33;
A34:= H2*(H2*M34 + S34) + B34; A35:= H2*(H2*M35 + S35) + B35;
A44:= H2*(H2*M44 + S44) + B44; A45:= H2*(H2*M45 + S45) + B45;
A55:= H2*(H2*M55 + S55) + B55;
comment STATIC CONDENSATION;
C1:= A15/A55; C2:= A25/A55; C3:= A35/A55; C4:= A45/A55;
B1:= (B1 - C1*B5)*H; B2:= (B2 - C2*B5)*H2;
B3:= (B3 - C3*B5)*H; B4:= (B4 - C4*B5)*H2;
A11:= (A11 - C1*A15)/H3; A12:= (A12 - C1*A25)/H2;
A13:= (A13 - C1*A35)/H3; A14:= (A14 - C1*A45)/H2;
A22:= (A22 - C2*A25)/H; A23:= (A23 - C2*A35)/H2;
A24:= (A24 - C2*A45)/H; A33:= (A33 - C3*A35)/H3;
A34:= (A34 - C3*A45)/H2; A44:= (A44 - C4*A45)/H;
end else
begin own real P5, Q5, R5, F5;
real X2, X3, X4, H, H2, H3,
P1, P2, P3, P4, Q1, Q2, Q3, Q4,
R1, R2, R3, R4, F1, F2, F3, F4,
B11, B12, B13, B14, B15, B16, B22, B23, B24, B25, B26,
B33, B34, B35, B36, B44, B45, B46, B55, B56, B66,
S11, S12, S13, S14, S15, S16, S22, S23, S24, S25, S26,
S33, S34, S35, S36, S44, S45, S46, S55, S56, S66,
M11, M12, M13, M14, M15, M16, M22, M23, M24, M25, M26,
M33, M34, M35, M36, M44, M45, M46, M55, M56, M66,
C15, C16, C25, C26, C35, C36, C45, C46, B5, B6,
A15, A16, A25, A26, A35, A36, A45, A46, A55, A56, A66, DET;
if L=1then
begin P5:= P(XL1); Q5:= Q(XL1); R5:= R(XL1); F5:= F(XL1)
end;
H:= XL - XL1; H2:= H*H; H3:= H*H2;
X2:= XL1 + H*.172673164646; X3:= XL1 + H/2; X4:= XL1 + XL - X2;
comment ELEMENT BENDING MATRIX;
P1:= P5; P2:= P(X2); P3:= P(X3); P4:= P(X4); P5:= P(XL);
B11:= + 105.8*P1 + 9.8*P5 + 7.3593121303513"-2*P2
+ 2.2755555555556"+1*P3 + 7.0565656088553"+0*P4;
B12:= + 27.6*P1 + 1.4*P5 - 3.41554824811"-1*P2
+ 2.8444444444444"+0*P3 + 1.0113960946522"+0*P4;
B13:= - 32.2*(P1 + P5) - 7.2063492063505"-1*(P2 + P4)
+ 2.2755555555556"+1*P3;
B14:= + 4.6*P1 + 8.4*P5 + 1.0328641222944"-1*P2
- 2.8444444444444"+0*P3 - 3.3445562534992"+0*P4;
B15:= - (B11 + B13); B16:= - (B12 + B13 + B14 + B15/2);
B22:= + 7.2*P1 + 0.2*P5 + 1.5851984028581"+0*P2
+ 3.5555555555556"-1*P3 + 1.4496032730059"-1*P4;
B23:= - 8.4*P1 - 4.6*P5 + 3.3445562534992"+0*P2
+ 2.8444444444444"+0*P3 - 1.0328641222944"-1*P4;
B24:= + 1.2*(P1 + P5) - 4.7936507936508"-1*(P2 + P4)
- 3.5555555555556"-1*P3;
B25:= - (B12 + B23); B26:= - (B22 + B23 + B24 + B25/2);
B33:= + 7.0565656088553"+0*P2 + 2.2755555555556"+1*P3
+ 7.3593121303513"-2*P4 + 105.8*P5 + 9.8*P1;
B34:= - 1.4*P1 - 27.6*P5 - 1.0113960946522"+0*P2
- 2.8444444444444"+0*P3 + 3.4155482481100"-1*P4;
B35:= - (B13 + B33); B36:= - (B23 + B33 + B34 + B35/2);
B44:= +7.2*P5 + P1/5 + 1.4496032730059"-1*P2
+ 3.5555555555556"-1*P3 + 1.5851984028581"+0*P4;
B45:= - (B14 + B34); B46:= - (B24 + B34 + B44 + B45/2);
B55:= - (B15 + B35); B56:= - (B16 + B36);
B66:= - (B26 + B36 + B46 + B56/2);
comment ELEMENT STIFFNESS MATRIX;
Q1:= Q5; Q2:= Q(X2); Q3:= Q(X3); Q4:= Q(X4); Q5:= Q(XL);
S11:= + 3.0242424037951"+0*Q2 + 3.1539909130065"-2*Q4;
S12:= + 1.2575525581744"-1*Q2 + 4.1767169716742"-3*Q4;
S13:= - 3.0884353741496"-1*(Q2+Q4);
S14:= + 4.0899041243062"-2*Q2 + 1.2842455355577"-2*Q4;
S15:= - (S13 + S11);
S16:= + 5.9254861177068"-1*Q2 + 6.0512612719116"-2*Q4;
S22:= + 5.2292052865422"-3*Q2 + 5.5310763862796"-4*Q4 + Q1/20;
S23:= - 1.2842455355577"-2*Q2 - 4.0899041243062"-2*Q4;
S24:= + 1.7006802721088"-3*(Q2+Q4);
S25:= - (S12 + S23);
S26:= + 2.4639593097426"-2*Q2 + 8.0134681270641"-3*Q4;
S33:= + 3.1539909130065"-2*Q2 + 3.0242424037951"+0*Q4;
S34:= - 4.1767169716742"-3*Q2 - 1.2575525581744"-1*Q4;
S35:= - (S13 + S33);
S36:= - 6.0512612719116"-2*Q2 - 5.9254861177068"-1*Q4;
S44:= + 5.5310763862796"-4*Q2 + 5.2292052865422"-3*Q4 + Q5/20;
S45:= - (S14 + S34);
S46:= + 8.0134681270641"-3*Q2 + 2.4639593097426"-2*Q4;
S55:= - (S15 + S35); S56:= -(S16 + S36);
S66:= + 1.1609977324263"-1*(Q2+Q4) + 3.5555555555556"-1*Q3;
comment ELEMENT MASS MATRIX;
R1:= R5; R2:= R(X2); R3:= R(X3); R4:= R(X4); R5:= R(XL);
M11:= + 9.7107020727310"-2*R2 + 1.5810259199180"-3*R4 + R1/20;
M12:= + 8.2354889460254"-3*R2 + 2.1932154960071"-4*R4;
M13:= + 1.2390670553936"-2*(R2+R4);
M14:= - 1.7188466249968"-3*R2 - 1.0508326752939"-3*R4;
M15:= + 5.3089789712119"-2*R2 + 6.7741558661060"-3*R4;
M16:= - 1.7377712856076"-2*R2 + 2.2173630018466"-3*R4;
M22:= + 6.9843846173145"-4*R2 + 3.0424512029349"-5*R4;
M23:= + 1.0508326752947"-3*R2 + 1.7188466249936"-3*R4;
M24:= - 1.4577259475206"-4*(R2+R4);
M25:= + 4.5024589679127"-3*R2 + 9.3971790283374"-4*R4;
M26:= - 1.4737756452780"-3*R2 + 3.0759488725998"-4*R4;
M33:= + 1.5810259199209"-3*R2 + 9.7107020727290"-2*R4 + R5/20;
M34:= - 2.1932154960131"-4*R2 - 8.2354889460254"-3*R4;
M35:= + 6.7741558661123"-3*R2 + 5.3089789712112"-2*R4;
M36:= - 2.2173630018492"-3*R2 + 1.7377712856071"-2*R4;
M44:= + 3.0424512029457"-5*R2 + 6.9843846173158"-4*R4;
M45:= - 9.3971790283542"-4*R2 - 4.5024589679131"-3*R4;
M46:= + 3.0759488726060"-4*R2 - 1.4737756452778"-3*R4;
M55:= + 2.9024943310657"-2*(R2+R4) + 3.5555555555556"-1*R3;
M56:= + 9.5006428402050"-3*(R4-R2);
M66:= + 3.1098153547125"-3*(R2+R4);
comment ELEMENT LOAD VECTOR;
F1:= F5; F2:= F(X2); F3:= F(X3); F4:= F(X4); F5:= F(XL);
B1:= + 1.6258748099336"-1*F2 + 2.0745852339969"-2*F4 + F1/20;
B2:= + 1.3788780589233"-2*F2 + 2.8778860774335"-3*F4;
B3:= + 2.0745852339969"-2*F2 + 1.6258748099336"-1*F4 + F5/20;
B4:= - 2.8778860774335"-3*F2 - 1.3788780589233"-2*F4;
B5:= + (F2 + F4)/11.25 + 3.5555555555556"-1*F3;
B6:= + 2.9095718698132"-2*(F4-F2);
A11:= H2*(H2*M11 + S11) + B11; A12:= H2*(H2*M12 + S12) + B12;
A13:= H2*(H2*M13 + S13) + B13; A14:= H2*(H2*M14 + S14) + B14;
A15:= H2*(H2*M15 + S15) + B15; A16:= H2*(H2*M16 + S16) + B16;
A22:= H2*(H2*M22 + S22) + B22; A23:= H2*(H2*M23 + S23) + B23;
A24:= H2*(H2*M24 + S24) + B24; A25:= H2*(H2*M25 + S25) + B25;
A26:= H2*(H2*M26 + S26) + B26; A33:= H2*(H2*M33 + S33) + B33;
A34:= H2*(H2*M34 + S34) + B34; A35:= H2*(H2*M35 + S35) + B35;
A36:= H2*(H2*M36 + S36) + B36; A44:= H2*(H2*M44 + S44) + B44;
A45:= H2*(H2*M45 + S45) + B45; A46:= H2*(H2*M46 + S46) + B46;
A55:= H2*(H2*M55 + S55) + B55; A56:= H2*(H2*M56 + S56) + B56;
A66:= H2*(H2*M66 + S66) + B66;
comment STATIC CONDENSATION;
DET:= - A55*A66 + A56*A56;
C15:= (A15*A66 - A16*A56)/DET; C16:= (A16*A55 - A15*A56)/DET;
C25:= (A25*A66 - A26*A56)/DET; C26:= (A26*A55 - A25*A56)/DET;
C35:= (A35*A66 - A36*A56)/DET; C36:= (A36*A55 - A35*A56)/DET;
C45:= (A45*A66 - A46*A56)/DET; C46:= (A46*A55 - A45*A56)/DET;
A11:= (A11 + C15*A15 + C16*A16)/H3;
A12:= (A12 + C15*A25 + C16*A26)/H2;
A13:= (A13 + C15*A35 + C16*A36)/H3;
A14:= (A14 + C15*A45 + C16*A46)/H2;
A22:= (A22 + C25*A25 + C26*A26)/H;
A23:= (A23 + C25*A35 + C26*A36)/H2;
A24:= (A24 + C25*A45 + C26*A46)/H;
A33:= (A33 + C35*A35 + C36*A36)/H3;
A34:= (A34 + C35*A45 + C36*A46)/H2;
A44:= (A44 + C45*A45 + C46*A46)/H;
B1:= (B1 + C15*B5 + C16*B6)*H; B2:= (B2 + C25*B5 + C26*B6)*H2;
B3:= (B3 + C35*B5 + C36*B6)*H; B4:= (B4 + C45*B5 + C46*B6)*H2;
end EL.MATVECEVAL.;
L:= 1; W:= V:= 0; N2:= N + N - 2; XL1:= X[0]; XL:= X[1];
YA:= E[1]; ZA:= E[2]; YB:= E[3]; ZB:= E[4];
ELEMENTMATVECEVALUATION; EM[2]:= "-12;
R1:= B3 - A13*YA - A23*ZA; D1:= A33; D2:= A44;
R2:= B4 - A14*YA - A24*ZA; E1:= A34;
for L:= L + 1while L<N"DO"
begin XL1:= XL; XL:= X[L];
ELEMENTMATVECEVALUATION;
A[W + 1]:= D1 + A11; A[W + 4]:= E1 + A12;
A[W + 7]:= A13; A[W + 10]:= A14;
A[W + 5]:= D2 + A22; A[W + 8]:= A23;
A[W + 11]:= A24; A[W + 14]:= 0;
Y[V + 1]:= R1 + B1; Y[V + 2]:= R2 + B2;
R1:= B3; R2:= B4; V:= V + 2; W:= W + 8;
D1:= A33; D2:= A44; E1:= A34
end;
L:= N; XL1:= XL; XL:= X[L]; ELEMENTMATVECEVALUATION;
Y[N2 - 1]:= R1 + B1 - A13*YB - A14*ZB;
Y[N2]:= R2 + B2 - A23*YB - A24*ZB;
A[W + 1]:= D1 + A11; A[W + 4]:= E1 + A12; A[W + 5]:= D2 + A22;
CHLDECSOLBND(A, N2, 3, EM, Y)
end FEMHERM;
eop