1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
|
SUBROUTINE BARD
C
C THIS SUBROUTINE PROCESSES BAR ELEMENT DATA TO PRODUCE STIFFNESS
C AND MASS MATRICES. IF THE HEAT TRANSFER OPTION IS ON, CONDUCTIVITY
C AND CAPACITY MATRICES ARE PRODUCED.
C
C DOUBLE PRECISION VERSION
C THIS ROUTINE WILL PRODUCE MASS MATRICES BY EITHER THE CONSISTENT
C OR CONVENTIONAL MASS METHODS.
C THE ECPT/EST ENTRIES FOR THE BAR (ELEMENT TYPE 34) ARE
C
C
C ECPT( 1) - IELID ELEMENT ID. NUMBER
C ECPT( 2) - ISILNO(2) * SCALAR INDEX NOS. OF THE GRID POINTS
C ECPT( 3) - ... *
C ECPT( 4) - SMALLV(3) $ REFERENCE VECTOR
C ECPT( 5) - ... $
C ECPT( 6) - ... $
C ECPT( 7) - ICSSV COOR. SYS. ID FOR SMALLV VECTOR
C ECPT( 8) - IPINFL(2) * PIN FLAGS
C ECPT( 9) - ... *
C ECPT(10) - ZA(3) $ OFFSET VECTOR FOR POINT A
C ECPT(11) - ... $
C ECPT(12) - ... $
C ECPT(13) - ZB(3) * OFFSET VECTOR FOR POINT B
C ECPT(14) - ... *
C ECPT(15) - ... *
C ECPT(16) - IMATID MATERIAL ID.
C ECPT(17) - A CROSS-SECTIONAL AREA
C ECPT(18) - I1 $ AREA MOMENTS OF INERTIA
C ECPT(19) - I2 $
C ECPT(20) - FJ TORSIONAL CONSTANT
C ECPT(21) - NSM NON-STRUCTURAL MASS
C ECPT(22) - FE FORCE ELEM DESCRIPTIONS (FORCE METHOD)
C ECPT(23) - C1 * STRESS RECOVERY COEFFICIENTS
C ECPT(24) - C2 *
C ECPT(25) - D1 *
C ECPT(26) - D2 *
C ECPT(27) - F1 *
C ECPT(28) - F2 *
C ECPT(29) - G1 *
C ECPT(30) - G2 *
C ECPT(31) - K1 $ AREA FACTORS FOR SHEAR
C ECPT(32) - K2 $
C ECPT(33) - I12 AREA MOMENT OF INERTIA
C ECPT(34) - MCSIDA COOR. SYS. ID. FOR GRID POINT A
C ECPT(35) - GPA(3) * BASIC COORDINATES FOR GRID POINT A
C ECPT(36) - ... *
C ECPT(37) - ... *
C ECPT(38) - MCSIDB COOR. SYS. ID. FOR GRID POINT B
C ECPT(39) - GPB(3) $ BASIC COORDINATES FOR GRID POINT B
C ECPT(40) - ... $
C ECPT(41) - ... $
C ECPT(42) - ELTEMP AVG. ELEMENT TEMPERATURE
C
LOGICAL BASIC,OFFSET,NOGO,AOFSET,BOFSET
INTEGER DICT(7),IS12OR(4),IS21OR(4),GSUBE,ESTID,IECPT(38)
REAL K1,K2,I1,I2,I12,NSM
DIMENSION ECPT(42),IPIN(10),IKK(4)
DOUBLE PRECISION CONST,BL22,BLSQ3,FM,KE,KK,KEP,M,MEP,ME,
1 LR1,LR2,LB,L2B3,L2B6,VECI(3),VECJ(3),VECK(3),
2 SMALVN,TA,TB,VEC,DELA,DELB,FL,SMALLV(3),
3 FLL,BL,BLSQ,BLCUBE,EI1,EI2,R1,R2,SK1,SK2,
4 SK3,SK4,AEL,GJL,BETA,BL13,BLSQ4,A2B,LIMIT,
5 EPSI,EPSI2
CHARACTER UFM*23,UWM*25
COMMON /XMSSG / UFM,UWM
COMMON /SYSTEM/ KSYSTM(100)
COMMON /EMGEST/ IELID,ISILNO(2),SMALV(3),ICSSV,IPINFL(2),ZA(3),
1 ZB(3),IMATID,A,I1,I2,FJ,NSM,FE,C1,C2,D1,D2,
2 F1,F2,G1,G2,K1,K2,I12,MCSIDA,GPA(3),MCSIDB,
3 GPB(3),TEMPEL
COMMON /EMGPRM/ IXTRA,JCORE,NCORE,DUM(12),ISTIF,IMASS,IDAMP,
1 IPREC,NOGO,HEAT,ICMBAR,LCSTM,LMAT,LHMAT
COMMON /EMGDIC/ IDUMM, LDICT,NGRIDS,ELID,ESTID
COMMON /EMGTRX/ KE(144),KEP(144),M(12,12),ME(144),MEP(144),
1 KK(144),SMALVN(6),TA(18),TB(9),VEC(10),
2 DELA(6),DELB(6)
COMMON /MATIN / MATIDC,MATFLG,ELTEMP,STRESS,SINTH,COSTH
COMMON /MATOUT/ E,G,NU,RHO,ALPHA,TSUBO,GSUBE,SIGT,SIGC,SIGS
COMMON /HMTOUT/ FK
EQUIVALENCE (KSYSTM(2),IOUTPT),(KSYSTM(56),IHEAT),
1 (ECPT(1),IECPT(1),IELID),(KSYSTM(87),KSY87),
2 (VEC(1),VECI(1)),(VEC(4),VECJ(1)),
3 (VEC(7),VECK(1))
DATA IKK / 1,7,73,79 /, EPSI,EPSI2 / 1.0D-18,1.0D-7 /
DATA IS12OR/ 1,37,109,73 /, IS21OR / 73,109,37,1 /
C
C
DICT(1) = ESTID
C
C SET UP POINTERS TO COOR. SYS. IDS., OFFSET VECTORS, AND PIN FLAGS.
C ICSIDA AND ICSIDB ARE COOR. SYS. IDS.
C
JCSIDA = 34
JCSIDB = 38
JOFSTA = 10
JOFSTB = 13
JPINA = 8
JPINB = 9
ICSIDA = IECPT(34)
ICSIDB = IECPT(38)
LIMIT = DBLE(IABS(KSY87))*.01D0
C
C NORMALIZE THE REFERENCE VECTOR WHICH LIES IN THE FIRST PRINCIPAL
C AXIS PLANE (FMMS - 36 P. 4)
C
FL = 0.D0
DO 40 I = 1,3
SMALLV(I) = SMALV(I)
40 FL = FL + SMALLV(I)**2
FL = DSQRT(FL)
IF (DABS(FL) .LT. EPSI) GO TO 7770
DO 50 I = 1,3
50 SMALVN(I) = SMALLV(I)/FL
C
C DETERMINE IF POINT A AND B ARE IN BASIC COORDINATES OR NOT.
C COMPUTE THE TRANSFORMATION MATRICES TA AND TB IF NECESSARY
C
IF (ICSIDA .NE. 0) CALL TRANSD(ECPT(JCSIDA),TA)
IF (ICSIDB .NE. 0) CALL TRANSD(ECPT(JCSIDB),TB)
C
C DETERMINE IF WE HAVE NON-ZERO OFFSET VECTORS.
C
AOFSET = .TRUE.
J = JOFSTA - 1
DO 70 I = 1,3
J = J + 1
IF (ECPT(J) .NE. 0.0) GO TO 80
70 CONTINUE
AOFSET = .FALSE.
80 BOFSET = .TRUE.
J = JOFSTB - 1
DO 90 I = 1,3
J = J + 1
IF (ECPT(J) .NE. 0.0) GO TO 100
90 CONTINUE
BOFSET = .FALSE.
C
C FORM THE CENTER AXIS OF THE BEAM WITHOUT OFFSETS.
C
100 DO 105 I = 1,3
JTA = I + JCSIDA
JTB = I + JCSIDB
105 VECI(I) = ECPT(JTA) - ECPT(JTB)
C
C SAVE IN A2B THE LENGTH OF BAR, WITHOUT OFFSET, FROM GRID PT. A
C TO B
C
A2B = DSQRT(VECI(1)**2 + VECI(2)**2 + VECI(3)**2)
C
C TRANSFORM THE OFFSET VECTORS IF NECESSARY
C
IF (.NOT.AOFSET .AND. .NOT.BOFSET) GO TO 150
C
C TRANSFORM THE OFFSET VECTOR FOR POINT A IF NECESSARY.
C
IDELA = 1
J = JOFSTA - 1
DO 110 I = 1,3
J = J + 1
110 DELA(I) = ECPT(J)
IF (ICSIDA .EQ. 0) GO TO 120
IDELA = 4
CALL GMMATD (TA, 3,3,0, DELA(1),3,1,0, DELA(4))
C
C TRANSFORM THE OFFSET VECTOR FOR POINT B IF NECESSARY
C
120 IDELB = 1
J = JOFSTB - 1
DO 130 I = 1,3
J = J + 1
130 DELB(I) = ECPT(J)
IF (ICSIDB .EQ. 0) GOTO 140
IDELB = 4
CALL GMMATD (TB, 3,3,0, DELB(1),3,1,0, DELB(4))
C
C SINCE THERE WAS AT LEAST ONE NON-ZERO OFFSET VECTOR RECOMPUTE VECI
C
140 DO 145 I = 1,3
JTA = I - 1 + IDELA
JTB = I - 1 + IDELB
145 VECI(I) = VECI(I)+DELA(JTA) - DELB(JTB)
C
C COMPUTE THE LENGTH OF THE BIG V (VECI) VECTOR AND NORMALIZE
C
150 FL = 0.D0
DO 155 I = 1,3
VECI(I) = -VECI(I)
155 FL = FL + VECI(I)**2
FL = DSQRT(FL)
IF (DABS(FL) .LT. EPSI) GO TO 7770
DO 160 I = 1,3
160 VECI(I) = VECI(I)/FL
C
C NOW THAT LENGTH HAS BEEN COMPUTED, CHECK POSSIBLE OFFSET ERROR
C ISSUE WARNING MESSAGE IF OFFSET EXCEEDS A2B BY 'LIMIT' PERCENT.
C (DEFAULT IS 15 PERCENT, KSYSTM(87) WORD)
C
IF (DABS(FL-A2B)/A2B .LE. LIMIT) GO TO 170
WRITE (IOUTPT,165) UWM,IELID
165 FORMAT (A25,' - UNUSUALLY LARGE OFFSET IS DETECTED FOR CBAR ',
1 'ELEMENT ID =',I8)
IF (KSY87 .LE. 0) GO TO 170
WRITE (IOUTPT,167) KSY87
167 FORMAT (/5X,'(OFFSET BAR LENGTH EXCEEDS NON-OFFSET LENGTH BY',
1 I4,' PERCENT, SET BY SYSTEM 87TH WORD)')
KSY87 = -KSY87
C
C BRANCH IF THIS IS A -HEAT- FORMULATION.
C
170 IF (IHEAT .EQ.1) GO TO 500
C
C COMPUTE THE SMALV0 VECTOR
C
ISV = 1
IF (ICSSV .EQ. 0) GO TO 180
ISV = 4
CALL GMMATD (TA,3,3,0, SMALVN(1),3,1,0, SMALVN(4))
C
C COMPUTE THE K VECTOR, VECK = VECI X SMALV0, AND NORMALIZE
C
180 VECK(1) = VECI(2)*SMALVN(ISV+2) - VECI(3)*SMALVN(ISV+1)
VECK(2) = VECI(3)*SMALVN(ISV ) - VECI(1)*SMALVN(ISV+2)
VECK(3) = VECI(1)*SMALVN(ISV+1) - VECI(2)*SMALVN(ISV )
FLL = DSQRT(VECK(1)**2 + VECK(2)**2 + VECK(3)**2)
IF (DABS(FLL) .LT. EPSI2) GO TO 7770
DO 190 I = 1,3
190 VECK(I) = VECK(I)/FLL
C
C COMPUTE THE J VECTOR, VECJ = VECK X VECI, AND NORMALIZE
C
VECJ(1) = VECK(2)*VECI(3) - VECK(3)*VECI(2)
VECJ(2) = VECK(3)*VECI(1) - VECK(1)*VECI(3)
VECJ(3) = VECK(1)*VECI(2) - VECK(2)*VECI(1)
FLL = DSQRT (VECJ(1)**2 + VECJ(2)**2 + VECJ(3)**2)
IF (DABS(FLL) .LT. EPSI2) GO TO 7770
VECJ(1) = VECJ(1)/FLL
VECJ(2) = VECJ(2)/FLL
VECJ(3) = VECJ(3)/FLL
C
C SEARCH THE MATERIAL PROPERTIES TABLE FOR E,G AND THE DAMPING
C CONSTANT.
C
MATIDC = IMATID
MATFLG = 1
ELTEMP = TEMPEL
CALL MAT (IECPT(1))
C
IF (ISTIF .EQ. 0) GOTO 600
C
C IF ELASTICITY AND SHEAR MODULES BOTH ZERO, SKIP STIFFNESS
C CALCULATION
C
IF (E.EQ.0. .AND. G.EQ.0.) GO TO 600
C
C SET UP INTERMEDIATE VARIABLES FOR ELEMENT STIFFNESS MATRIX
C CALCULATION
C
ASSIGN 305 TO K OR M
205 BL = FL
BLSQ = FL**2
BLCUBE= BLSQ*BL
C
C COMPUTE SOME TERMS TO BE USED IN STIFFNESS MATRIX KE
C
EI1 = DBLE(E)*DBLE(I1)
EI2 = DBLE(E)*DBLE(I2)
IF (K1.EQ.0.0 .OR. I12.NE.0.0) GO TO 210
GAK1 = DBLE(G)*DBLE(A)*DBLE(K1)
R1 = (12.D0*EI1*GAK1)/(GAK1*BLCUBE + 12.D0*BL*EI1)
GO TO 220
210 R1 = 12.D0*EI1/BLCUBE
220 IF (K2.EQ.0.0 .OR. I12.NE.0.0) GO TO 230
GAK2 = DBLE(G)*DBLE(A)*DBLE(K2)
R2 = (12.D0*EI2*GAK2)/(GAK2*BLCUBE + 12.D0*BL*EI2)
GO TO 240
230 R2 = 12.D0*EI2/BLCUBE
C
240 SK1 = .25D0*R1*BLSQ + EI1/BL
SK2 = .25D0*R2*BLSQ + EI2/BL
SK3 = .25D0*R1*BLSQ - EI1/BL
SK4 = .25D0*R2*BLSQ - EI2/BL
C
AEL = DBLE(A)*DBLE(E)/BL
LR1 = BL*R1/2.D0
LR2 = BL*R2/2.D0
GJL = DBLE(G)*DBLE(FJ)/BL
C
C
C CONSTRUCT THE GENERAL 12X12 MATRIX FOR THE BAR ELEMENT
C
C ** **
C * K K *
C * AA AB*
C K = * T *
C * K K *
C * AB BB*
C ** **
C
C
C
C FIRST SET THE COMPONENT CODE AND THE DOF
C
ICODE = 63
NDOF = 12
NSQ = NDOF**2
C
C CONSTRUCT THE 12 X 12 MATRIX KE
C
C ** **
C * 1 73 *
C * 14 62 86 134 *
C * 27 51 99 123 *
C * 40 112 *
C * 29 53 101 125 *
C * 18 66 90 138 *
C * 7 79 *
C * 20 68 92 140 *
C * 33 57 105 129 *
C * 46 118 *
C * 35 59 107 131 *
C * 24 72 96 144 *
C ** **
C
DO 300 I = 1,144
300 KE(I) = 0.D0
KE( 1) = AEL
KE( 7) = -AEL
KE( 14) = R1
KE( 18) = LR1
KE( 20) = -R1
KE( 24) = LR1
KE( 27) = R2
KE( 29) = -LR2
KE( 33) = -R2
KE( 35) = -LR2
KE( 40) = GJL
KE( 46) = -GJL
KE( 51) = -LR2
KE( 53) = SK2
KE( 57) = LR2
KE( 59) = SK4
KE( 62) = LR1
KE( 66) = SK1
KE( 68) = -LR1
KE( 72) = SK3
KE( 73) = -AEL
KE( 79) = AEL
KE( 86) = -R1
KE( 90) = -LR1
KE( 92) = R1
KE( 96) = -LR1
KE( 99) = -R2
KE(101) = LR2
KE(105) = R2
KE(107) = LR2
KE(112) = -GJL
KE(118) = GJL
KE(123) = -LR2
KE(125) = SK4
KE(129) = LR2
KE(131) = SK2
KE(134) = LR1
KE(138) = SK3
KE(140) = -LR1
KE(144) = SK1
IF (I12 .EQ. 0.) GOTO 303
BETA = 12.D0*DBLE(E)*DBLE(I12)/BLCUBE
LB = BL*BETA/2.0D0
L2B3 = BLSQ*BETA/3.0D0
L2B6 = BLSQ*BETA/6.0D0
KE( 15) = BETA
KE( 17) = -LB
KE( 21) = -BETA
KE( 23) = -LB
KE( 26) = BETA
KE( 30) = LB
KE( 32) = -BETA
KE( 36) = LB
KE( 50) = -LB
KE( 54) = -L2B3
KE( 56) = LB
KE( 60) = -L2B6
KE( 63) = LB
KE( 65) = -L2B3
KE( 69) = -LB
KE( 71) = -L2B6
KE( 87) = -BETA
KE( 89) = LB
KE( 93) = BETA
KE( 95) = LB
KE( 98) = -BETA
KE(102) = -LB
KE(104) = BETA
KE(108) = -LB
KE(122) = -LB
KE(126) = -L2B6
KE(128) = LB
KE(132) = -L2B3
KE(135) = LB
KE(137) = -L2B6
KE(141) = -LB
KE(143) = -L2B3
303 GO TO K OR M, (305,640,465,750)
C
C DETERMINE IF THERE ARE NON-ZERO PIN FLAGS.
C
305 KA = IECPT(JPINA)
KB = IECPT(JPINB)
IF (KA.EQ.0 .AND. KB.EQ.0) GOTO 345
C
C SET UP THE IPIN ARRAY
C
DO 310 I = 1,5
IPIN(I ) = MOD(KA,10)
IPIN(I+5) = MOD(KB,10) + 6
IF (IPIN(I+5) .EQ. 6) IPIN(I+5) = 0
KA = KA/10
310 KB = KB/10
C
C ALTER KE MATRIX DUE TO PIN FLAGS.
C
DO 340 I = 1,10
IF (IPIN(I) .EQ.0) GO TO 340
II = 13*IPIN(I) - 12
IF (KE(II) .NE. 0.D0) GO TO 320
IL = IPIN(I)
II = II - IL
DO 315 J = 1,12
II = II + 1
KE(II) = 0.D0
KE(IL) = 0.D0
IL = IL + 12
315 CONTINUE
GO TO 340
320 DO 330 J = 1,12
JI = 12*(J-1) + IPIN(I)
IJ = 12*(IPIN(I)-1) + J
DO 325 LL = 1,12
JLL = 12*(J-1) + LL
ILL = 12*(IPIN(I)-1) + LL
KEP(JLL) = KE(JLL) - (KE(ILL)/KE(II))*KE(JI)
325 CONTINUE
KEP(IJ) = 0.D0
KEP(JI) = 0.D0
330 CONTINUE
DO 335 K = 1,144
335 KE(K) = KEP(K)
340 CONTINUE
C
C DIVIDE KE INTO FOUR SUBMATRICES AND STORE IN OPEN CORE
C
C E E E
C K = KK(1 TO 36) K = KK(37 TO 72) K = KK(73 TO 108)
C AA AB BA
C
C E
C K = KK(109 TO 144)
C BB
C
C
345 J = 0
DO 355 I = 1,72,12
LOW = I
LIM = I + 5
DO 350 K = LOW,LIM
J = J + 1
KK(J ) = KE(K )
KK(J+36) = KE(K+ 6)
KK(J+72) = KE(K+72)
350 KK(J+108)= KE(K+78)
355 CONTINUE
C
ASSIGN 465 TO K OR M
C
C ZERO OUT THE ARRAY WHERE THE 3X3 MATRIX H AND THE W AND W 6X6
C MATRICES WILL RESIDE. A B
C T
C A MATRIX NOW STORED IN KE
C
358 DO 357 I = 1,9
357 KE(I) = VEC(I)
C
DO 360 I = 28,144
360 KE(I) = 0.D0
C
C
C SET UP POINTERS
C
BASIC = ICSIDA.EQ.0
JCSID = JCSIDA
OFFSET = AOFSET
JOFSET = JOFSTA
DO 395 I = 1,2
IWBEG = I*36
C
C SET UP THE -G- MATRIX. IG POINTS TO THE BEGINNING OF THE G MATRIX
C G = AT X TI
C
IG = 1
IF (BASIC) GO TO 380
CALL TRANSD (ECPT(JCSID),KE(10))
CALL GMMATD (KE(1), 3,3,0, KE(10), 3,3,0, KE(19))
IG = 19
C
C IF THERE IS A NON-ZERO OFFSET FOR THE POINT, SET UP THE D 3X3
C MATRIX.
C
380 IF (.NOT.OFFSET) GO TO 385
KE(10) = 0.D0
KE(11) = ECPT(JOFSET+2)
KE(12) = -ECPT(JOFSET+1)
KE(13) = -KE(11)
KE(14) = 0.D0
KE(15) = ECPT(JOFSET)
KE(16) = -KE(12)
KE(17) = -KE(15)
KE(18) = 0.D0
C
C FORM THE 3X3 PRODUCT H = G X D, I.E., KE(28) = KE(IG) X KE(10)
C
CALL GMMATD (KE(IG), 3,3,0, KE(10), 3,3,0, KE(28))
C
C
C FORM THE W SUBMATRICES IN KE(37) AND KE(73)
C
C
385 KE(IWBEG+ 1) = KE(IG )
KE(IWBEG+ 2) = KE(IG+1)
KE(IWBEG+ 3) = KE(IG+2)
KE(IWBEG+ 7) = KE(IG+3)
KE(IWBEG+ 8) = KE(IG+4)
KE(IWBEG+ 9) = KE(IG+5)
KE(IWBEG+13) = KE(IG+6)
KE(IWBEG+14) = KE(IG+7)
KE(IWBEG+15) = KE(IG+8)
KE(IWBEG+22) = KE(IG )
KE(IWBEG+23) = KE(IG+1)
KE(IWBEG+24) = KE(IG+2)
KE(IWBEG+28) = KE(IG+3)
KE(IWBEG+29) = KE(IG+4)
KE(IWBEG+30) = KE(IG+5)
KE(IWBEG+34) = KE(IG+6)
KE(IWBEG+35) = KE(IG+7)
KE(IWBEG+36) = KE(IG+8)
IF (.NOT.OFFSET) GO TO 390
KE(IWBEG+ 4) = KE(28)
KE(IWBEG+ 5) = KE(29)
KE(IWBEG+ 6) = KE(30)
KE(IWBEG+10) = KE(31)
KE(IWBEG+11) = KE(32)
KE(IWBEG+12) = KE(33)
KE(IWBEG+16) = KE(34)
KE(IWBEG+17) = KE(35)
KE(IWBEG+18) = KE(36)
390 BASIC = ICSIDB.EQ.0
JCSID = JCSIDB
OFFSET = BOFSET
JOFSET = JOFSTB
395 CONTINUE
C
C CONVERT THE K PARTITIONS TO GLOBAL COORDINATES AND STORE IN KEP
C
IAFT = 37
DO 400 I = 1,4
IKX = (I-1)*36 + 1
IK = IKX
IF (I .GE. 3) IKX = (7-I-1)*36 + 1
IFORE = ((I-1)/2)*36 + 37
CALL GMMATD (KE(IFORE), 6,6,1, KK(IKX), 6,6,0, KE(109))
CALL GMMATD (KE(109), 6,6,0, KE(IAFT), 6,6,0, KEP(IK))
IAFT = 73
IF (I .EQ. 3) IAFT = 37
400 CONTINUE
C
C REFORM THE K MATRIX (12X12) FROM THE FOUR SUBMATRICES (6X6) AND
C ORDER THE SUBMATRICES BY INCREASING SIL VALUE
C
DO 460 II = 1,4
IX1 = IKK(II)
IX2 = IX1 + 60
IS = IS12OR(II)
IF (ISILNO(1) .GT. ISILNO(2)) IS = IS21OR(II)
DO 450 I = IX1,IX2,12
IP5 = I + 5
DO 440 J = I,IP5
KE(J) = KEP(IS)
440 IS = IS + 1
450 CONTINUE
460 CONTINUE
C
GO TO K OR M, (305,640,465,750)
C
C OUTPUT THE STIFFNESS MATRIX
C
465 DICT(2) = 1
DICT(3) = NDOF
DICT(4) = ICODE
DICT(5) = GSUBE
CALL EMGOUT (KE(1),KE(1),NSQ,1,DICT,1,IPREC)
GO TO 600
C
C THE MASS MATRIX IS GENERATED HERE. IF THE PARAMETER ICMBAR IS
C .LT. 0, CALL THE CONVENTIONAL MASS MATRIX GENERATION ROUTINE FOR
C THE BAR. OTHERWISE CALL THE ROUTINE TO GENERATE CONSISTENT MASS
C MATRICES FOR THE BAR.
C
600 CONST = (FL*(DBLE(RHO)*DBLE(A) + DBLE(NSM)))/420.D0
IF (IMASS.EQ.0 .OR. CONST.EQ.0.D0) RETURN
IF (ICMBAR .LT. 0) GO TO 800
C
C CALCULATE THE CONSISTENT/CONVENTIONAL MASS MATRIX
C
C CALL THE MAT ROUTINE TO FETCH SINGLE PRECISION MATERIAL PROPERTIES
C
MATIDC = IMATID
MATFLG = 1
ELTEMP = TEMPEL
CALL MAT (IECPT(1))
C
C COMPUTE TERMS OF THE ELEMENT MASS MATRIX
C
BL22 = 22.D0*FL
BL13 = 13.D0*FL
BLSQ4 = 4.D0*FL**2
BLSQ3 = 3.D0*FL**2
C
C CONSTRUCT THE ELEMENT MASS MATRIX.
C
DO 610 I = 1,12
DO 610 J = 1,12
610 M( I, J) = 0.D0
M( 1, 1) = 175.D0
M( 1, 7) = 35.D0
M( 2, 2) = 156.D0
M( 2, 6) = BL22
M( 2, 8) = 54.D0
M( 2,12) =-BL13
M( 3, 3) = 156.D0
M( 3, 5) =-BL22
M( 3, 9) = 54.D0
M( 3,11) = BL13
M( 5, 5) = BLSQ4
M( 5, 9) =-BL13
M( 5,11) =-BLSQ3
M( 6, 6) = BLSQ4
M( 6, 8) = BL13
M( 6,12) =-BLSQ3
M( 7, 7) = 175.D0
M( 8, 8) = 156.D0
M( 8,12) =-BL22
M( 9, 9) = 156.D0
M( 9,11) = BL22
M(11,11) = BLSQ4
M(12,12) = BLSQ4
C
C STORE THE UPPER TRIANGULAR PART OF THE MATRIX IN THE LOWER PART.
C
DO 625 I = 1,10
LOW = I + 1
DO 620 J = LOW,12
M(J,I) = M(I,J)
620 CONTINUE
625 CONTINUE
C
C MULTIPLY BY CONSTANT AND STORE ROW-WISE IN THE ARRAY ME
C
K = 0
DO 630 I = 1,12
DO 630 J = 1,12
K = K + 1
630 ME(K) = CONST*M(I,J)
C
C IF THERE ARE NO PIN FLAGS THERE IS NO NEED TO CALCULATE THE
C ELEMENT STIFFNESS MATRIX
C
KA = IECPT(JPINA)
KB = IECPT(JPINB)
IF (KA.EQ.0 .AND. KB.EQ.0) GO TO 705
C
C COMPUTE THE STIFFNESS MATRIX KE
C
C
ASSIGN 640 TO K OR M
GO TO 205
C
C RETURN HERE AFTER COMPUTING THE STIFFNESS MATRIX
C
C
C SET UP TNHE IPIN ARRAY
C
640 DO 645 I = 1,5
IPIN(I ) = MOD(KA,10)
IPIN(I+5) = MOD(KB,10)+6
IF (IPIN(I+5) .EQ. 6) IPIN(I+5) = 0
KA = KA/10
645 KB = KB/10
C
C ALTER THE ELEMENT MASS MATRIX DUE TO PIN FLAGS. NOTE THAT THE
C FOLLOWING CODE IS CONGRUENT AS IT WERE TO THE CODE IN SUBROUTINE
C DBEAM IN THE DSMG1 MODULE.
C
DO 700 J = 1,10
IF (IPIN(J) .EQ. 0) GO TO 700
JJ = 12*(IPIN(J)-1) + IPIN(J)
IF (KE(JJ) .EQ. 0.) GO TO 680
DO 660 I = 1,12
JI = 12*(IPIN(J)-1) + I
IJ = 12*(I-1) + IPIN(J)
DO 650 L1 = 1,12
IL = 12*(I -1) + L1
LJ = 12*(L1-1) + IPIN(J)
MEP(IL) = ME(IL) - KE(LJ)*ME(JI)/KE(JJ) - KE(JI)*ME(LJ)/KE(JJ)
1 + KE(LJ)*KE(JI)*ME(JJ)/KE(JJ)**2
650 CONTINUE
660 CONTINUE
DO 670 K = 1,144
670 ME(K) = MEP(K)
C
C ZERO OUT THE IPIN(J) TH ROW AND COLUMN OF ME
C
680 J1 = JJ - IPIN(J)
J2 = IPIN(J)
DO 690 K = 1,12
J1 = J1 + 1
ME(J1) = 0.D0
ME(J2) = 0.D0
690 J2 = J2 + 12
700 CONTINUE
C
C E E E
C STORE M AT KK(1 TO 36), M AT KK (37 TO 72), M AT KK(73 TO 108)
C AA AB BA
C
C E
C M AT KK(109 TO 144)
C BB
C
705 J = 0
DO 720 I = 1,72,12
LOW = I
LIM = LOW + 5
DO 710 K = LOW,LIM
J = J + 1
KK(J) = ME(K)
KK(J+ 36) = ME(K+ 6)
KK(J+ 72) = ME(K+72)
710 KK(J+108) = ME(K+78)
720 CONTINUE
C
C CALCULATE THE TRANSFORMATION VECTORS
C
ASSIGN 750 TO K OR M
GO TO 358
C
C OUTPUT THE CONSISTENT MASS MATRIX
C
750 DICT(2) = 1
DICT(3) = NDOF
DICT(4) = ICODE
DICT(5) = 0
CALL EMGOUT (KE(1),KE(1),144,1,DICT,2,IPREC)
RETURN
C
C CALCULATE THE CONVENTIONAL MASS MATRIX HERE
C
C
C GET RHO FROM MPT BY CALLING MAT
C
800 MATIDC = IMATID
MATFLG = 4
ELTEMP = TEMPEL
CALL MAT (ECPT(1))
DO 810 I = 1,72
810 MEP(I) = 0.D0
FM = .5D0*FL*(DBLE(RHO)*DBLE(A) + DBLE(NSM))
C
C DETERMINE IF THE GRID POINT IS ASSOCIATED WITH A NON-ZERO OFFSET.
C
JOFSET = 9
DO 850 II = 1,2
IX = (II-1)*36
J = JOFSET
DO 815 I = 1,3
J = J + 1
IF (ECPT(J) .NE. 0.) GO TO 820
815 CONTINUE
GO TO 840
C
C FORM UPPER RIGHT CORNER OF THE MATRIX
C
820 MEP(IX+ 1) = 1.D0
MEP(IX+ 8) = 1.D0
MEP(IX+15) = 1.D0
MEP(IX+ 5) = ECPT(JOFSET+3)
MEP(IX+ 6) =-ECPT(JOFSET+2)
MEP(IX+12) = ECPT(JOFSET+1)
MEP(IX+10) =-MEP(IX+ 5)
MEP(IX+16) =-MEP(IX+ 6)
MEP(IX+17) =-MEP(IX+12)
MEP(IX+20) =-MEP(IX+ 5)
MEP(IX+21) =-MEP(IX+ 6)
MEP(IX+25) =-MEP(IX+10)
MEP(IX+27) =-MEP(IX+12)
MEP(IX+31) =-MEP(IX+16)
MEP(IX+32) =-MEP(IX+17)
MEP(IX+22) = ECPT(JOFSET+3)**2 + ECPT(JOFSET+2)**2
MEP(IX+29) = ECPT(JOFSET+3)**2 + ECPT(JOFSET+1)**2
MEP(IX+36) = ECPT(JOFSET+2)**2 + ECPT(JOFSET+1)**2
MEP(IX+23) =-ECPT(JOFSET+1)*ECPT(JOFSET+2)
MEP(IX+24) =-ECPT(JOFSET+1)*ECPT(JOFSET+3)
MEP(IX+30) =-ECPT(JOFSET+2)*ECPT(JOFSET+3)
MEP(IX+28) = MEP(IX+23)
MEP(IX+34) = MEP(IX+24)
MEP(IX+35) = MEP(IX+30)
C
C MULTIPLY M BY THE CONSTANT FL
C
DO 830 I = 1,36
IS = IX + I
830 MEP(IS) = MEP(IS)*FM
GO TO 850
C
C HERE WE HAVE A ZERO OFFSET VECTOR
C
840 MEP(IX+ 1) = FM
MEP(IX+ 8) = FM
MEP(IX+15) = FM
850 JOFSET = 12
C
C INSERT THE M AND M SUBMATRICES INTO M ACCORDING TO INCREASING
C SIL A B
C
DO 860 I = 1,144
860 ME(I) = 0.D0
C
IF (ISILNO(1) - ISILNO(2)) 870,870,880
870 IX1 = 1
IX2 = 37
GO TO 890
880 IX1 = 37
IX2 = 1
890 CONTINUE
DO 900 JJ = 1,36
MM = MOD(JJ,6)
IF (MM .EQ. 0) MM = 6
I = ((JJ-1)/6)*12 + MM
J = I + 78
ME(I) = MEP(IX1)
ME(J) = MEP(IX2)
IX1 = IX1 + 1
900 IX2 = IX2 + 1
C
C OUTPUT THE CONVENTIONAL MASS MATRIX
C
DICT(2) = 1
DICT(3) = NDOF
DICT(4) = ICODE
DICT(5) = 0
C
CALL EMGOUT (ME,ME,144,1,DICT,2,IPREC)
C
RETURN
C
C HEAT FORMULATION CONTINUES HERE. GET MATERIAL PROPERTY -K- FROM
C HMAT
C
500 MATFLG = 1
MATIDC = IECPT(16)
ELTEMP = ECPT(42)
DICT(2) = 1
DICT(3) = 2
DICT(4) = 1
DICT(5) = 0
IF (ISTIF .EQ. 0) GO TO 540
CALL HMAT (IELID)
C
KK(1) = DBLE(FK)*DBLE(ECPT(17))/FL
IF (KK(1).EQ. 0.D0) GO TO 520
KK(2) =-KK(1)
KK(3) = KK(2)
KK(4) = KK(1)
CALL EMGOUT (KK(1),KK(1),4,1,DICT,1,IPREC)
C
520 MATFLG = 4
C
C ERROR IN NEXT CARD FOR HEAT FORMULATION. REMOVED BY G.CHAN/SPERRY,
C 1984. ALSO, CHANGE GO TO 520 TO 540, 11-TH CARD ABOVE, AND
C CALL EMGOUT BELOW AND WRITE TO THE 3-RD FILE INSTEAD OF THE 2-ND.
C
CALL HMAT (IELID)
KK(1) = (DBLE(FK)*DBLE(ECPT(17)))*FL/2.D0
IF (KK(1) .EQ. 0.D0) RETURN
KK(2) = KK(1)
DICT(2) = 2
CALL EMGOUT (KK(1),KK(1),2,1,DICT,3,IPREC)
540 RETURN
C
C ERROR RETURNS
C
7770 CONTINUE
WRITE (IOUTPT,7775) UFM,IELID
7775 FORMAT (A23,' 3176, BAR ELEMENT ID',I9,
1 ' HAS ILLEGAL GEOMETRY OR CONNECTIONS.')
NOGO = .TRUE.
RETURN
END
|
