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 34 of file fsigpini.F.

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

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

◆ pdefoi()

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

Definition at line 318 of file fsigpini.F.

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

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

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