OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
fsigpini.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine fsigpini (fxbelm, iparg, x, pm, ixp, geo, fxbmod, fxbsig, r, nelp, ibeam_vector, rbeam_vector)
subroutine pevecii (x1, y1, z1, x2, y2, z2, r, rloc, al, nel, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)
subroutine pdefoi (v, exx, exy, exz, al, nel, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)
subroutine pcurvi (v, geo, kxx, kyy, kzz, exy, exz, al, nel, mgm, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)
subroutine pm1inif (pm, for, mom, eint, geo, exx, exy, exz, kxx, kyy, kzz, al, nel, mat, mgm)

Function/Subroutine Documentation

◆ fsigpini()

subroutine fsigpini ( integer, dimension(*) fxbelm,
integer, dimension(nparg,*) iparg,
x,
pm,
integer, dimension(nixp,*) ixp,
geo,
fxbmod,
fxbsig,
r,
integer nelp,
integer, dimension(nelp), intent(in) ibeam_vector,
dimension(3,nelp), intent(in) rbeam_vector )

Definition at line 33 of file fsigpini.F.

36C-----------------------------------------------
37C I m p l i c i t T y p e s
38C-----------------------------------------------
39#include "implicit_f.inc"
40C-----------------------------------------------
41C G l o b a l P a r a m e t e r s
42C-----------------------------------------------
43#include "mvsiz_p.inc"
44C-----------------------------------------------
45C C o m m o n B l o c k s
46C-----------------------------------------------
47#include "param_c.inc"
48C-----------------------------------------------
49C D u m m y A r g u m e n t s
50C-----------------------------------------------
51 INTEGER FXBELM(*), IPARG(NPARG,*), IXP(NIXP,*), NELP
52 INTEGER, INTENT (IN ) :: IBEAM_VECTOR(NELP)
53 my_real
54 . x(3,*), pm(npropm,*), geo(npropg,*), fxbmod(*),
55 . fxbsig(*), r(3,*)
56 my_real , INTENT (IN ) :: rbeam_vector(3,nelp)
57C-----------------------------------------------
58C L o c a l V a r i a b l e s
59C-----------------------------------------------
60 INTEGER IG, OFFSET, LAST, NFT, NFS, I, NG, IEL,
61 . N1, N2
62 INTEGER MAT(MVSIZ), PROP(MVSIZ)
63 my_real
64 . ee1x(mvsiz), ee1y(mvsiz), ee1z(mvsiz),
65 . ee2x(mvsiz), ee2y(mvsiz), ee2z(mvsiz),
66 . ee3x(mvsiz), ee3y(mvsiz), ee3z(mvsiz)
67 my_real
68 . vl(3,2,mvsiz), vrl(3,2,mvsiz)
69 my_real
70 . x1(mvsiz), y1(mvsiz), z1(mvsiz),
71 . x2(mvsiz), y2(mvsiz), z2(mvsiz),
72 . x3(mvsiz), y3(mvsiz), z3(mvsiz)
73 my_real
74 . e2x, e2y, e2z, ee2, rloc(3,mvsiz),
75 . d11, d12, d13, d21, d22, d23,
76 . dr11, dr12, dr13, dr21, dr22, dr23,
77 . al(mvsiz)
78 my_real
79 . for(3,mvsiz), mom(3,mvsiz), eint(2,mvsiz),
80 . exx(mvsiz), exy(mvsiz), exz(mvsiz),
81 . kxx(mvsiz), kyy(mvsiz), kzz(mvsiz)
82C-----------------------------------------------
83C
84 DO ig=1,nelp,mvsiz
85 offset=ig-1
86 last=min(mvsiz,nelp-offset)
87 nft=offset*9
88 nfs=offset*8
89 DO i=1,last
90 ng=fxbelm(nft+9*(i-1)+1)
91 iel=iparg(3,ng)+fxbelm(nft+9*(i-1)+2)
92 mat(i)=ixp(1,iel)
93 prop(i)=ixp(5,iel)
94 x1(i)=x(1,ixp(2,iel))
95 y1(i)=x(2,ixp(2,iel))
96 z1(i)=x(3,ixp(2,iel))
97 x2(i)=x(1,ixp(3,iel))
98 y2(i)=x(2,ixp(3,iel))
99 z2(i)=x(3,ixp(3,iel))
100 x3(i)=x(1,ixp(4,iel))
101 y3(i)=x(2,ixp(4,iel))
102 z3(i)=x(3,ixp(4,iel))
103 IF (ibeam_vector(iel) > 1) THEN
104 e2x=rbeam_vector(1,iel)
105 e2y=rbeam_vector(2,iel)
106 e2z=rbeam_vector(3,iel)
107 ELSE
108 e2x=x3(i)-x1(i)
109 e2y=y3(i)-y1(i)
110 e2z=z3(i)-z1(i)
111 ENDIF
112 ee2=sqrt(e2x**2+e2y**2+e2z**2)
113 rloc(1,i)=e2x/ee2
114 rloc(2,i)=e2y/ee2
115 rloc(3,i)=e2z/ee2
116 n1=fxbelm(nft+9*(i-1)+3)
117 n2=fxbelm(nft+9*(i-1)+4)
118 d11=fxbmod(6*(n1-1)+1)
119 d12=fxbmod(6*(n1-1)+2)
120 d13=fxbmod(6*(n1-1)+3)
121 d21=fxbmod(6*(n2-1)+1)
122 d22=fxbmod(6*(n2-1)+2)
123 d23=fxbmod(6*(n2-1)+3)
124 vl(1,1,i)=r(1,1)*d11+r(1,2)*d12+r(1,3)*d13
125 vl(2,1,i)=r(2,1)*d11+r(2,2)*d12+r(2,3)*d13
126 vl(3,1,i)=r(3,1)*d11+r(3,2)*d12+r(3,3)*d13
127 vl(1,2,i)=r(1,1)*d21+r(1,2)*d22+r(1,3)*d23
128 vl(2,2,i)=r(2,1)*d21+r(2,2)*d22+r(2,3)*d23
129 vl(3,2,i)=r(3,1)*d21+r(3,2)*d22+r(3,3)*d23
130 dr11=fxbmod(6*(n1-1)+4)
131 dr12=fxbmod(6*(n1-1)+5)
132 dr13=fxbmod(6*(n1-1)+6)
133 dr21=fxbmod(6*(n2-1)+4)
134 dr22=fxbmod(6*(n2-1)+5)
135 dr23=fxbmod(6*(n2-1)+6)
136 vrl(1,1,i)=r(1,1)*dr11+r(1,2)*dr12+r(1,3)*dr13
137 vrl(2,1,i)=r(2,1)*dr11+r(2,2)*dr12+r(2,3)*dr13
138 vrl(3,1,i)=r(3,1)*dr11+r(3,2)*dr12+r(3,3)*dr13
139 vrl(1,2,i)=r(1,1)*dr21+r(1,2)*dr22+r(1,3)*dr23
140 vrl(2,2,i)=r(2,1)*dr21+r(2,2)*dr22+r(2,3)*dr23
141 vrl(3,2,i)=r(3,1)*dr21+r(3,2)*dr22+r(3,3)*dr23
142 for(1,i)=zero
143 for(2,i)=zero
144 for(3,i)=zero
145 mom(1,i)=zero
146 mom(2,i)=zero
147 mom(3,i)=zero
148 ENDDO
149C
150 CALL pevecii(x1, y1, z1, x2, y2,
151 . z2, vrl, rloc, al, last,
152 . ee1x, ee1y, ee1z,
153 . ee2x, ee2y, ee2z,
154 . ee3x, ee3y, ee3z)
155C
156 CALL pdefoi(vl , exx , exy, exz, al, last,
157 . ee1x, ee1y, ee1z,
158 . ee2x, ee2y, ee2z,
159 . ee3x, ee3y, ee3z)
160 CALL pcurvi(vrl, geo , kxx , kyy , kzz ,
161 . exy , exz , al , last, prop,
162 . ee1x, ee1y, ee1z,
163 . ee2x, ee2y, ee2z,
164 . ee3x, ee3y, ee3z)
165C
166 CALL pm1inif(pm, for, mom , eint, geo,
167 . exx, exy, exz , kxx , kyy,
168 . kzz, al , last, mat , prop)
169C
170 DO i=1,last
171 fxbsig(nfs+8*(i-1)+1)=for(1,i)
172 fxbsig(nfs+8*(i-1)+2)=for(2,i)
173 fxbsig(nfs+8*(i-1)+3)=for(3,i)
174 fxbsig(nfs+8*(i-1)+4)=mom(1,i)
175 fxbsig(nfs+8*(i-1)+5)=mom(2,i)
176 fxbsig(nfs+8*(i-1)+6)=mom(3,i)
177 fxbsig(nfs+8*(i-1)+7)=eint(1,i)
178 fxbsig(nfs+8*(i-1)+8)=eint(2,i)
179 ENDDO
180 ENDDO
181C
182 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
pevecii
subroutine pevecii(x1, y1, z1, x2, y2, z2, r, rloc, al, nel, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)
Definition fsigpini.F:194
pm1inif
subroutine pm1inif(pm, for, mom, eint, geo, exx, exy, exz, kxx, kyy, kzz, al, nel, mat, mgm)
Definition fsigpini.F:490
pdefoi
subroutine pdefoi(v, exx, exy, exz, al, nel, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)
Definition fsigpini.F:320
pcurvi
subroutine pcurvi(v, geo, kxx, kyy, kzz, exy, exz, al, nel, mgm, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z)
Definition fsigpini.F:383
min
#define min(a, b)
Definition macros.h:20
for
for(i8=*sizetab-1;i8 >=0;i8--)
Definition mumps_common.c:91

◆ pcurvi()

subroutine pcurvi ( v,
geo,
kxx,
kyy,
kzz,
exy,
exz,
al,
integer nel,
integer, dimension(*) mgm,
e1x,
e1y,
e1z,
e2x,
e2y,
e2z,
e3x,
e3y,
e3z )

Definition at line 378 of file fsigpini.F.

383C-----------------------------------------------
384C I m p l i c i t T y p e s
385C-----------------------------------------------
386#include "implicit_f.inc"
387C-----------------------------------------------
388C G l o b a l P a r a m e t e r s
389C-----------------------------------------------
390#include "mvsiz_p.inc"
391C-----------------------------------------------
392C C o m m o n B l o c k s
393C-----------------------------------------------
394#include "param_c.inc"
395C-----------------------------------------------
396C D u m m y A r g u m e n t s
397C-----------------------------------------------
398 INTEGER :: NEL, MGM(*)
399 my_real
400 . v(3,2,*), geo(npropg,*), kxx(*), kyy(*), kzz(*),
401 . exy(*), exz(*), al(*),
402 . e1x(*), e1y(*), e1z(*),
403 . e2x(*), e2y(*), e2z(*),
404 . e3x(*), e3y(*), e3z(*)
405C-----------------------------------------------
406C L o c a l V a r i a b l e s
407C-----------------------------------------------
408 INTEGER I, IG, IRX, IR1Y, IR1Z, IR2Y, IR2Z, IRY, IRZ
409 my_real
410 . rx1g(mvsiz), ry1g(mvsiz), rz1g(mvsiz),
411 . rx2g(mvsiz), ry2g(mvsiz), rz2g(mvsiz),
412 . rx1(mvsiz), ry1(mvsiz), rz1(mvsiz),
413 . rx2(mvsiz), ry2(mvsiz), rz2(mvsiz),
414 . rxav(mvsiz), ryav(mvsiz), rzav(mvsiz)
415C
416 DO i=1,nel
417 rx1g(i)=v(1,1,i)
418 ry1g(i)=v(2,1,i)
419 rz1g(i)=v(3,1,i)
420 rx2g(i)=v(1,2,i)
421 ry2g(i)=v(2,2,i)
422 rz2g(i)=v(3,2,i)
423 ENDDO
424C
425 DO i=1,nel
426 rx1(i)=e1x(i)*rx1g(i)+e1y(i)*ry1g(i)+e1z(i)*rz1g(i)
427 ry1(i)=e2x(i)*rx1g(i)+e2y(i)*ry1g(i)+e2z(i)*rz1g(i)
428 rz1(i)=e3x(i)*rx1g(i)+e3y(i)*ry1g(i)+e3z(i)*rz1g(i)
429 rx2(i)=e1x(i)*rx2g(i)+e1y(i)*ry2g(i)+e1z(i)*rz2g(i)
430 ry2(i)=e2x(i)*rx2g(i)+e2y(i)*ry2g(i)+e2z(i)*rz2g(i)
431 rz2(i)=e3x(i)*rx2g(i)+e3y(i)*ry2g(i)+e3z(i)*rz2g(i)
432 ENDDO
433C---------------------------------------------------
434C Free rotations
435C---------------------------------------------------
436 DO i=1,nel
437 ig=mgm(i)
438 irx =nint(geo(7 ,ig))
439 ir1y=nint(geo(8 ,ig))
440 ir1z=nint(geo(9 ,ig))
441 ir2y=nint(geo(10,ig))
442 ir2z=nint(geo(11,ig))
443 iry =min(1,ir1y+ir2y)
444 irz =min(1,ir1z+ir2z)
445 rx1(i)=rx1(i)*irx
446 ry1(i)=ry1(i)*iry
447 rz1(i)=rz1(i)*irz
448 rx2(i)=rx2(i)*irx
449 ry2(i)=ry2(i)*iry
450 rz2(i)=rz2(i)*irz
451 exz(i)=exz(i)*iry
452 exy(i)=exy(i)*irz
453 ry1(i)=ir1y*ry1(i)
454 + -(one -ir1y)*(three_half*exz(i)+half*ry2(i))
455 ry2(i)=ir2y*ry2(i)
456 + -(one -ir2y)*(three_half*exz(i)+half*ry1(i))
457 rz1(i)=ir1z*rz1(i)
458 + +(one-ir1z)*(three_half*exy(i)-half*rz2(i))
459 rz2(i)=ir2z*rz2(i)
460 + +(one -ir2z)*(three_half*exy(i)-half*rz1(i))
461 ENDDO
462C
463 DO i=1,nel
464 kxx(i)=(rx2(i)-rx1(i))/al(i)
465 kyy(i)=(ry2(i)-ry1(i))/al(i)
466 kzz(i)=(rz2(i)-rz1(i))/al(i)
467 ENDDO
468C
469 DO i=1,nel
470 rxav(i)=rx1(i)+rx2(i)
471 rzav(i)=rz1(i)+rz2(i)
472 ryav(i)=ry1(i)+ry2(i)
473 ENDDO
474C
475 DO i=1,nel
476 exz(i)=exz(i) + half*ryav(i)
477 exy(i)=exy(i) - half*rzav(i)
478 ENDDO
479C
480 RETURN

◆ pdefoi()

subroutine pdefoi ( v,
exx,
exy,
exz,
al,
integer nel,
e1x,
e1y,
e1z,
e2x,
e2y,
e2z,
e3x,
e3y,
e3z )

Definition at line 316 of file fsigpini.F.

320C-----------------------------------------------
321C I m p l i c i t T y p e s
322C-----------------------------------------------
323#include "implicit_f.inc"
324C-----------------------------------------------
325C G l o b a l P a r a m e t e r s
326C-----------------------------------------------
327#include "mvsiz_p.inc"
328C-----------------------------------------------
329C D u m m y A r g u m e n t s
330C-----------------------------------------------
331 INTEGER :: NEL
332 my_real
333 . v(3,2,*), exx(*), exy(*), exz(*), al(*),
334 . e1x(*), e1y(*), e1z(*),
335 . e2x(*), e2y(*), e2z(*),
336 . e3x(*), e3y(*), e3z(*)
337C-----------------------------------------------
338C L o c a l V a r i a b l e s
339C-----------------------------------------------
340 INTEGER I
341 my_real
342 . vx1g(mvsiz), vy1g(mvsiz), vz1g(mvsiz),
343 . vx2g(mvsiz), vy2g(mvsiz), vz2g(mvsiz),
344 . vx1(mvsiz), vy1(mvsiz), vz1(mvsiz),
345 . vx2(mvsiz), vy2(mvsiz), vz2(mvsiz)
346C
347 DO i=1,nel
348 vx1g(i)=v(1,1,i)
349 vy1g(i)=v(2,1,i)
350 vz1g(i)=v(3,1,i)
351 vx2g(i)=v(1,2,i)
352 vy2g(i)=v(2,2,i)
353 vz2g(i)=v(3,2,i)
354 ENDDO
355C
356 DO i=1,nel
357 vx1(i)=e1x(i)*vx1g(i)+e1y(i)*vy1g(i)+e1z(i)*vz1g(i)
358 vy1(i)=e2x(i)*vx1g(i)+e2y(i)*vy1g(i)+e2z(i)*vz1g(i)
359 vz1(i)=e3x(i)*vx1g(i)+e3y(i)*vy1g(i)+e3z(i)*vz1g(i)
360 vx2(i)=e1x(i)*vx2g(i)+e1y(i)*vy2g(i)+e1z(i)*vz2g(i)
361 vy2(i)=e2x(i)*vx2g(i)+e2y(i)*vy2g(i)+e2z(i)*vz2g(i)
362 vz2(i)=e3x(i)*vx2g(i)+e3y(i)*vy2g(i)+e3z(i)*vz2g(i)
363 ENDDO
364C
365 DO i=1,nel
366 exx(i)=(vx2(i)-vx1(i))/al(i)
367 exy(i)=(vy2(i)-vy1(i))/al(i)
368 exz(i)=(vz2(i)-vz1(i))/al(i)
369 ENDDO
370C
371 RETURN

◆ pevecii()

subroutine pevecii ( x1,
y1,
z1,
x2,
y2,
z2,
r,
rloc,
al,
integer nel,
e1x,
e1y,
e1z,
e2x,
e2y,
e2z,
e3x,
e3y,
e3z )

Definition at line 189 of file fsigpini.F.

194C-----------------------------------------------
195C I m p l i c i t T y p e s
196C-----------------------------------------------
197#include "implicit_f.inc"
198C-----------------------------------------------
199C G l o b a l P a r a m e t e r s
200C-----------------------------------------------
201#include "mvsiz_p.inc"
202C-----------------------------------------------
203C D u m m y A r g u m e n t s
204C-----------------------------------------------
205 INTEGER :: NEL
206 my_real
207 . x1(*), y1(*), z1(*), x2(*), y2(*), z2(*),
208 . r(3,2,*), rloc(3,*), al(*),
209 . e1x(*), e1y(*), e1z(*),
210 . e2x(*), e2y(*), e2z(*),
211 . e3x(*), e3y(*), e3z(*)
212C-----------------------------------------------
213C L o c a l V a r i a b l e s
214C-----------------------------------------------
215 INTEGER I
216 my_real
217 . rx1g(mvsiz), ry1g(mvsiz), rz1g(mvsiz),
218 . rx2g(mvsiz), ry2g(mvsiz), rz2g(mvsiz),
219 . rx1(mvsiz),
220 . rx2(mvsiz),
221 . theta, sum2(mvsiz), sum3(mvsiz), sum(mvsiz),
222 . cost(mvsiz), sint(mvsiz)
223
224C
225 DO i=1,nel
226 rx1g(i)=r(1,1,i)
227 ry1g(i)=r(2,1,i)
228 rz1g(i)=r(3,1,i)
229 rx2g(i)=r(1,2,i)
230 ry2g(i)=r(2,2,i)
231 rz2g(i)=r(3,2,i)
232 ENDDO
233C
234 DO i=1,nel
235 e2x(i)=rloc(1,i)
236 e2y(i)=rloc(2,i)
237 e2z(i)=rloc(3,i)
238 ENDDO
239C
240 DO i=1,nel
241 e1x(i)=x2(i)-x1(i)
242 e1y(i)=y2(i)-y1(i)
243 e1z(i)=z2(i)-z1(i)
244 ENDDO
245C
246 DO i=1,nel
247 al(i)=sqrt(e1x(i)**2+e1y(i)**2+e1z(i)**2)
248 ENDDO
249C
250 DO i=1,nel
251 e1x(i)=e1x(i)/al(i)
252 e1y(i)=e1y(i)/al(i)
253 e1z(i)=e1z(i)/al(i)
254 ENDDO
255C
256 DO i=1,nel
257 e3x(i)=e1y(i)*e2z(i)-e1z(i)*e2y(i)
258 e3y(i)=e1z(i)*e2x(i)-e1x(i)*e2z(i)
259 e3z(i)=e1x(i)*e2y(i)-e1y(i)*e2x(i)
260 ENDDO
261C
262 DO i=1,nel
263 e2x(i)=e3y(i)*e1z(i)-e3z(i)*e1y(i)
264 e2y(i)=e3z(i)*e1x(i)-e3x(i)*e1z(i)
265 e2z(i)=e3x(i)*e1y(i)-e3y(i)*e1x(i)
266 ENDDO
267C--------------------------------------------
268C Average torsion in global coordinates
269C--------------------------------------------
270 DO i=1,nel
271 rx1(i)=e1x(i)*rx1g(i)+e1y(i)*ry1g(i)+e1z(i)*rz1g(i)
272 rx2(i)=e1x(i)*rx2g(i)+e1y(i)*ry2g(i)+e1z(i)*rz2g(i)
273 theta=half*(rx1(i)+rx2(i))
274 sum2(i)=sqrt(e2x(i)**2+e2y(i)**2+e2z(i)**2)
275 sum3(i)=sqrt(e3x(i)**2+e3y(i)**2+e3z(i)**2)
276 cost(i)=cos(theta)/sum2(i)
277 sint(i)=sin(theta)/sum3(i)
278 ENDDO
279C
280 DO i=1,nel
281 e2x(i)=e2x(i)*cost(i)+e3x(i)*sint(i)
282 e2y(i)=e2y(i)*cost(i)+e3y(i)*sint(i)
283 e2z(i)=e2z(i)*cost(i)+e3z(i)*sint(i)
284 ENDDO
285C
286 DO i=1,nel
287 sum(i)=sqrt(e2x(i)**2+e2y(i)**2+e2z(i)**2)
288 ENDDO
289C
290 DO i=1,nel
291 e2x(i)=e2x(i)/sum(i)
292 e2y(i)=e2y(i)/sum(i)
293 e2z(i)=e2z(i)/sum(i)
294 ENDDO
295C
296 DO i=1,nel
297 e3x(i)=e1y(i)*e2z(i)-e1z(i)*e2y(i)
298 e3y(i)=e1z(i)*e2x(i)-e1x(i)*e2z(i)
299 e3z(i)=e1x(i)*e2y(i)-e1y(i)*e2x(i)
300 ENDDO
301C
302 DO i=1,nel
303 sum(i)=sqrt(e3x(i)**2+e3y(i)**2+e3z(i)**2)
304 e3x(i)=e3x(i)/sum(i)
305 e3y(i)=e3y(i)/sum(i)
306 e3z(i)=e3z(i)/sum(i)
307 ENDDO
308C
309 RETURN

◆ pm1inif()

subroutine pm1inif ( pm,
for,
mom,
eint,
geo,
exx,
exy,
exz,
kxx,
kyy,
kzz,
al,
integer nel,
integer, dimension(*) mat,
integer, dimension(*) mgm )

Definition at line 487 of file fsigpini.F.

490C-----------------------------------------------
491C I m p l i c i t T y p e s
492C-----------------------------------------------
493#include "implicit_f.inc"
494C-----------------------------------------------
495C G l o b a l P a r a m e t e r s
496C-----------------------------------------------
497#include "mvsiz_p.inc"
498C-----------------------------------------------
499C C o m m o n B l o c k s
500C-----------------------------------------------
501#include "param_c.inc"
502C-----------------------------------------------
503C D u m m y A r g u m e n t s
504C-----------------------------------------------
505 INTEGER :: NEL, MAT(*), MGM(*)
506 my_real
507 . pm(npropm,*), for(3,*), mom(3,*), eint(2,*),
508 . geo(npropg,*), exx(*), exy(*), exz(*), kxx(*),
509 . kyy(*), kzz(*), al(*)
510C-----------------------------------------------
511C L o c a l V a r i a b l e s
512C-----------------------------------------------
513 INTEGER I
514 my_real
515 . rho(mvsiz), g(mvsiz), ym(mvsiz), a1(mvsiz), b1(mvsiz),
516 . b2(mvsiz), b3(mvsiz), shf(mvsiz), sh(mvsiz),
517 . yma2(mvsiz), sh10(mvsiz), sh20(mvsiz), sh0(mvsiz),
518 . sh1(mvsiz), sh2(mvsiz), degmb(mvsiz), degfx(mvsiz)
519C
520 DO i=1,nel
521 rho(i) =pm( 1,mat(i))
522 g(i) =pm(22,mat(i))
523 ym(i) =pm(20,mat(i))
524 a1(i) =geo(1,mgm(i))
525 b1(i) =geo(2,mgm(i))
526 b2(i) =geo(18,mgm(i))
527 b3(i) =geo(4,mgm(i))
528 shf(i) =geo(37,mgm(i))
529 ENDDO
530C
531C Transverse shear computed with K1=12EI/L**2 K2=5/6GA
532 DO i=1,nel
533 sh(i)=five_over_6*g(i)*a1(i)
534 yma2(i)=twelve*ym(i)/al(i)**2
535 sh10(i)=yma2(i)*b1(i)
536 sh20(i)=yma2(i)*b2(i)
537 sh0(i)=(one-shf(i))*sh(i)
538 sh1(i)=sh0(i)*sh10(i)/(sh(i)+sh10(i)) + shf(i)*sh10(i)
539 sh2(i)=sh0(i)*sh20(i)/(sh(i)+sh20(i)) + shf(i)*sh20(i)
540C
541 for(1,i)=for(1,i)+ exx(i)*a1(i)*ym(i)
542 for(2,i)=for(2,i)+ exy(i)*sh2(i)
543 for(3,i)=for(3,i)+ exz(i)*sh1(i)
544 mom(1,i)=mom(1,i)+ kxx(i)*g(i)*b3(i)
545 mom(2,i)=mom(2,i)+ kyy(i)*ym(i)*b1(i)
546 mom(3,i)=mom(3,i)+ kzz(i)*ym(i)*b2(i)
547 ENDDO
548C
549 DO i=1,nel
550 degmb(i) = for(1,i)*exx(i)+for(2,i)*exy(i)+for(3,i)*exz(i)
551 degfx(i) = mom(1,i)*kxx(i)+mom(2,i)*kyy(i)+mom(3,i)*kzz(i)
552 ENDDO
553C
554 DO i=1,nel
555 eint(1,i) = degmb(i)*al(i)*half
556 eint(2,i) = degfx(i)*al(i)*half
557 ENDDO
558C
559 RETURN