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

Go to the source code of this file.

Functions/Subroutines

subroutine czlken3 (jft, jlt, vol, thk0, thk2, hm, hz, a_i, px1, px2, py1, py2, hxx, hyy, hxy, ph1, ph2, z1, nplat, iplat, dhz, k11, k12, k13, k14, k22, k23, k24, k33, k34, k44, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44, mf11, mf12, mf13, mf14, mf22, mf23, mf24, mf33, mf34, mf44, fm12, fm13, fm14, fm23, fm24, fm34, idril)
subroutine czlkenr3 (jft, jlt, vol, thk0, thk2, hm, hz, a_i, px1, px2, py1, py2, hxx, hyy, hxy, ph1, ph2, z1, nplat, iplat, dhz, k11, k12, k13, k14, k22, k23, k24, k33, k34, k44, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44, mf11, mf12, mf13, mf14, mf22, mf23, mf24, mf33, mf34, mf44, fm12, fm13, fm14, fm23, fm24, fm34, phkrx, phkry, phkrxy, pherx, phery, pherxy, phkrz, pherz, phkx, phky, phex, phey)

Function/Subroutine Documentation

◆ czlken3()

subroutine czlken3 ( integer jft,
integer jlt,
vol,
thk0,
thk2,
hm,
hz,
a_i,
px1,
px2,
py1,
py2,
hxx,
hyy,
hxy,
ph1,
ph2,
z1,
integer nplat,
integer, dimension(*) iplat,
dhz,
k11,
k12,
k13,
k14,
k22,
k23,
k24,
k33,
k34,
k44,
m11,
m12,
m13,
m14,
m22,
m23,
m24,
m33,
m34,
m44,
mf11,
mf12,
mf13,
mf14,
mf22,
mf23,
mf24,
mf33,
mf34,
mf44,
fm12,
fm13,
fm14,
fm23,
fm24,
fm34,
integer idril )

Definition at line 28 of file czlken3.F.

37C--------------------------------------------------------------------------------------------------
38C-----------------------------------------------
39C I M P L I C I T T Y P E S
40C-----------------------------------------------
41#include "implicit_f.inc"
42#include "mvsiz_p.inc"
43C-----------------------------------------------
44C D U M M Y A R G U M E N T S
45C-----------------------------------------------
46C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
47 INTEGER JFT,JLT,NPLAT,IPLAT(*),IDRIL
48 my_real
49 . vol(*),thk0(*),thk2(*),a_i(*),z1(*),hm(mvsiz,4),hz(*),
50 . px1(*) ,px2(*) ,py1(*) ,py2(*) ,ph1(*) ,ph2(*) ,
51 . hxx(*),hyy(*),hxy(*),
52 . k11(3,3,*),k12(3,3,*),k13(3,3,*),k14(3,3,*),
53 . k22(3,3,*),k23(3,3,*),k24(3,3,*),k33(3,3,*),
54 . m11(3,3,*),m12(3,3,*),m13(3,3,*),m14(3,3,*),
55 . m22(3,3,*),m23(3,3,*),m24(3,3,*),m33(3,3,*),
56 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),mf14(3,3,*),
57 . mf22(3,3,*),mf23(3,3,*),mf24(3,3,*),mf33(3,3,*),
58 . fm12(3,3,*),fm13(3,3,*),fm14(3,3,*),
59 . fm23(3,3,*),fm24(3,3,*),fm34(3,3,*),
60 . k34(3,3,*),k44(3,3,*),m34(3,3,*),m44(3,3,*),
61 . mf34(3,3,*),mf44(3,3,*),dhz(*)
62C---------------|[KIJ][MFIJ]|----
63C-----KE(6x6)= | |
64C---------------|[FMIJ]{MIJ]|----
65C-----------------------------------------------
66C L O C A L V A R I A B L E S
67C-----------------------------------------------
68 INTEGER EP,I,J,NF,M
69 my_real
70 . dm(2,2,mvsiz),c1,c2,gm(mvsiz),
71 . g11(mvsiz),g12(mvsiz),g13(mvsiz),
72 . g14(mvsiz),g22(mvsiz),g23(mvsiz),
73 . g24(mvsiz),g33(mvsiz),g34(mvsiz),g44(mvsiz),
74 . chm(2,2,mvsiz),chf(2,2,mvsiz),facf(mvsiz),fac,zr(mvsiz),
75 . gama(4,mvsiz),cxx,cyy,cxy,c11,c22,c12,cx1,cx2,cy1,cy2
76C-----------gama(I)=hI/4-PH(I), PH:anti-sym comme bxI------
77 nf=nplat+1
78C-----------Attention Matrice sym Kii ne calcul que la moitie---------72
79#include "vectorize.inc"
80 DO m=jft,jlt
81 ep=iplat(m)
82 c2=vol(ep)
83 c1=thk2(ep)*c2
84 dm(1,1,m)=hm(ep,1)*c2
85 dm(2,2,m)=hm(ep,2)*c2
86 dm(1,2,m)=hm(ep,3)*c2
87 dm(2,1,m)=dm(1,2,m)
88 gm(m) =hm(ep,4)*c2
89 facf(m)=one_over_12*thk2(ep)
90C DHZ(M)= HZ(EP)*C1
91 zr(m) = z1(ep)
92 ENDDO
93C------------------shear est constante----
94 DO ep=jft,jlt
95 chm(1,1,ep)=dm(1,1,ep)*hxx(ep)
96C+GM(EP)*HYY(EP)
97 chm(2,2,ep)=dm(2,2,ep)*hyy(ep)
98C+GM(EP)*HXX(EP)
99 chm(1,2,ep)=dm(1,2,ep)*hxy(ep)
100C+GM(EP)*HXY(EP)
101 chm(2,1,ep)=chm(1,2,ep)
102C------------------pour etre coherant avec QEPH explicite---
103 chf(1,1,ep)=facf(ep)*chm(2,2,ep)
104C+GM(EP)*HXX(EP))
105 chf(2,2,ep)=facf(ep)*chm(1,1,ep)
106C+GM(EP)*HYY(EP))
107 chf(1,2,ep)=-facf(ep)*chm(1,2,ep)
108C+GM(EP)*HXY(EP))
109c FAC = FACF(EP)*GM(EP)
110c CHF(1,1,EP)=CHF(1,1,EP)+FAC*HXX(EP)
111c CHF(2,2,EP)=CHF(2,2,EP)+FAC*HYY(EP)
112c CHF(1,2,EP)=CHF(1,2,EP)-FAC*HXY(EP)
113C
114 chf(2,1,ep)=chf(1,2,ep)
115 ENDDO
116C
117 IF (idril>0) THEN
118 DO ep=jft,jlt
119 chm(1,1,ep)=chm(1,1,ep)+gm(ep)*hyy(ep)
120 chm(2,2,ep)=chm(2,2,ep)+gm(ep)*hxx(ep)
121 chm(1,2,ep)=chm(1,2,ep)+gm(ep)*hxy(ep)
122 chm(2,1,ep)=chm(1,2,ep)
123C----------------------------
124 fac = facf(ep)*gm(ep)
125 chf(1,1,ep)=chf(1,1,ep)+fac*hxx(ep)
126 chf(2,2,ep)=chf(2,2,ep)+fac*hyy(ep)
127 chf(1,2,ep)=chf(1,2,ep)-fac*hxy(ep)
128 chf(2,1,ep)=chf(1,2,ep)
129 ENDDO
130 END IF !(IDRIL>0) THEN
131C--------------gamaI--------
132 DO ep=jft,jlt
133 gama(1,ep)=fourth-ph1(ep)
134 gama(2,ep)=-fourth-ph2(ep)
135 gama(3,ep)=fourth+ph1(ep)
136 gama(4,ep)=-fourth+ph2(ep)
137 g11(ep) =gama(1,ep)*gama(1,ep)
138 g12(ep) =gama(1,ep)*gama(2,ep)
139 g13(ep) =gama(1,ep)*gama(3,ep)
140 g14(ep) =gama(1,ep)*gama(4,ep)
141 g22(ep) =gama(2,ep)*gama(2,ep)
142 g23(ep) =gama(2,ep)*gama(3,ep)
143 g24(ep) =gama(2,ep)*gama(4,ep)
144 g33(ep) =gama(3,ep)*gama(3,ep)
145 g34(ep) =gama(3,ep)*gama(4,ep)
146 g44(ep) =gama(4,ep)*gama(4,ep)
147 ENDDO
148C-------
149 DO i=1,2
150 DO j=i,2
151 DO ep=jft,jlt
152 k11(i,j,ep)=k11(i,j,ep)+chm(i,j,ep)*g11(ep)
153 k22(i,j,ep)=k22(i,j,ep)+chm(i,j,ep)*g22(ep)
154 k33(i,j,ep)=k33(i,j,ep)+chm(i,j,ep)*g33(ep)
155 k44(i,j,ep)=k44(i,j,ep)+chm(i,j,ep)*g44(ep)
156 m11(i,j,ep)=m11(i,j,ep)+chf(i,j,ep)*g11(ep)
157 m22(i,j,ep)=m22(i,j,ep)+chf(i,j,ep)*g22(ep)
158 m33(i,j,ep)=m33(i,j,ep)+chf(i,j,ep)*g33(ep)
159 m44(i,j,ep)=m44(i,j,ep)+chf(i,j,ep)*g44(ep)
160 ENDDO
161 ENDDO
162 ENDDO
163 DO i=1,2
164 DO j=1,2
165 DO ep=jft,jlt
166 k12(i,j,ep)=k12(i,j,ep)+chm(i,j,ep)*g12(ep)
167 k13(i,j,ep)=k13(i,j,ep)+chm(i,j,ep)*g13(ep)
168 k14(i,j,ep)=k14(i,j,ep)+chm(i,j,ep)*g14(ep)
169 k23(i,j,ep)=k23(i,j,ep)+chm(i,j,ep)*g23(ep)
170 k24(i,j,ep)=k24(i,j,ep)+chm(i,j,ep)*g24(ep)
171 k34(i,j,ep)=k34(i,j,ep)+chm(i,j,ep)*g34(ep)
172 m12(i,j,ep)=m12(i,j,ep)+chf(i,j,ep)*g12(ep)
173 m13(i,j,ep)=m13(i,j,ep)+chf(i,j,ep)*g13(ep)
174 m14(i,j,ep)=m14(i,j,ep)+chf(i,j,ep)*g14(ep)
175 m23(i,j,ep)=m23(i,j,ep)+chf(i,j,ep)*g23(ep)
176 m24(i,j,ep)=m24(i,j,ep)+chf(i,j,ep)*g24(ep)
177 m34(i,j,ep)=m34(i,j,ep)+chf(i,j,ep)*g34(ep)
178 ENDDO
179 ENDDO
180 ENDDO
181C---+---------+---------+warped elements-------------
182 DO ep=nf,jlt
183 cxx =zr(ep)*chm(1,1,ep)
184 cyy =zr(ep)*chm(2,2,ep)
185 cxy =zr(ep)*chm(1,2,ep)
186 cx1 = cxx*px1(ep)+cxy*py1(ep)
187 cx2 = cxx*px2(ep)+cxy*py2(ep)
188 cy1 = cxy*px1(ep)+cyy*py1(ep)
189 cy2 = cxy*px2(ep)+cyy*py2(ep)
190 k11(1,3,ep) = cx1*gama(1,ep)
191 k11(2,3,ep) = cy1*gama(1,ep)
192 k11(3,1,ep) = k11(1,3,ep)
193 k11(3,2,ep) = k11(2,3,ep)
194 k22(1,3,ep) = cx2*gama(2,ep)
195 k22(2,3,ep) = cy2*gama(2,ep)
196 k22(3,1,ep) = k22(1,3,ep)
197 k22(3,2,ep) = k22(2,3,ep)
198 k33(1,3,ep) = -cx1*gama(3,ep)
199 k33(2,3,ep) = -cy1*gama(3,ep)
200 k33(3,1,ep) = k33(1,3,ep)
201 k33(3,2,ep) = k33(2,3,ep)
202 k44(1,3,ep) = -cx2*gama(4,ep)
203 k44(2,3,ep) = -cy2*gama(4,ep)
204 k44(3,1,ep) = k44(1,3,ep)
205 k44(3,2,ep) = k44(2,3,ep)
206 k12(1,3,ep) = cx2*gama(1,ep)
207 k12(2,3,ep) = cy2*gama(1,ep)
208 k12(3,1,ep) = cx1*gama(2,ep)
209 k12(3,2,ep) = cy1*gama(2,ep)
210 k13(1,3,ep) = -k11(1,3,ep)
211 k13(2,3,ep) = -k11(2,3,ep)
212 k13(3,1,ep) = -k33(1,3,ep)
213 k13(3,2,ep) = -k33(2,3,ep)
214 k14(1,3,ep) = -k12(1,3,ep)
215 k14(2,3,ep) = -k12(2,3,ep)
216 k14(3,1,ep) = cx1*gama(4,ep)
217 k14(3,2,ep) = cy1*gama(4,ep)
218 k23(1,3,ep) = -k12(3,1,ep)
219 k23(2,3,ep) = -k12(3,2,ep)
220 k23(3,1,ep) = cx2*gama(3,ep)
221 k23(3,2,ep) = cy2*gama(3,ep)
222 k24(1,3,ep) = -k22(1,3,ep)
223 k24(2,3,ep) = -k22(2,3,ep)
224 k24(3,1,ep) = -k44(1,3,ep)
225 k24(3,2,ep) = -k44(2,3,ep)
226 k34(1,3,ep) = -k23(3,1,ep)
227 k34(2,3,ep) = -k23(3,2,ep)
228 k34(3,1,ep) = -k14(3,1,ep)
229 k34(3,2,ep) = -k14(3,2,ep)
230 ENDDO
231 DO ep=nf,jlt
232 fac = zr(ep)*zr(ep)
233 cxx =fac*chm(1,1,ep)
234 cyy =fac*chm(2,2,ep)
235 cxy =fac*chm(1,2,ep)
236 c11 = cxx*px1(ep)*px1(ep)+cyy*py1(ep)*py1(ep)+
237 1 cxy*(px1(ep)*py1(ep)+py1(ep)*px1(ep))
238 c12 = cxx*px1(ep)*px2(ep)+cyy*py1(ep)*py2(ep)+
239 1 cxy*(px1(ep)*py2(ep)+py1(ep)*px2(ep))
240 c22 = cxx*px2(ep)*px2(ep)+cyy*py2(ep)*py2(ep)+
241 1 cxy*(px2(ep)*py2(ep)+py2(ep)*px2(ep))
242 k11(3,3,ep) = k11(3,3,ep)+c11
243 k12(3,3,ep) = k12(3,3,ep)+c12
244 k13(3,3,ep) = k13(3,3,ep)-c11
245 k14(3,3,ep) = k14(3,3,ep)-c12
246 k22(3,3,ep) = k22(3,3,ep)+c22
247 k23(3,3,ep) = k23(3,3,ep)-c12
248 k24(3,3,ep) = k24(3,3,ep)-c22
249 k33(3,3,ep) = k33(3,3,ep)+c11
250 k34(3,3,ep) = k34(3,3,ep)+c12
251 k44(3,3,ep) = k44(3,3,ep)+c22
252 ENDDO
253C------ Mzz ----------
254 IF (idril == 0) THEN
255 DO ep=jft,jlt
256 fac =dhz(ep)*(hxx(ep)+hyy(ep))
257 m11(3,3,ep)=m11(3,3,ep)+fac*g11(ep)
258 m12(3,3,ep)=m12(3,3,ep)+fac*g12(ep)
259 m13(3,3,ep)=m13(3,3,ep)+fac*g13(ep)
260 m14(3,3,ep)=m14(3,3,ep)+fac*g14(ep)
261 m22(3,3,ep)=m22(3,3,ep)+fac*g22(ep)
262 m23(3,3,ep)=m23(3,3,ep)+fac*g23(ep)
263 m24(3,3,ep)=m24(3,3,ep)+fac*g24(ep)
264 m33(3,3,ep)=m33(3,3,ep)+fac*g33(ep)
265 m34(3,3,ep)=m34(3,3,ep)+fac*g34(ep)
266 m44(3,3,ep)=m44(3,3,ep)+fac*g44(ep)
267 ENDDO
268 END IF !(IDRIL == 0) THEN
269C
270 RETURN
my_real
#define my_real
Definition cppsort.cpp:32

◆ czlkenr3()

subroutine czlkenr3 ( integer jft,
integer jlt,
vol,
thk0,
thk2,
hm,
hz,
a_i,
px1,
px2,
py1,
py2,
hxx,
hyy,
hxy,
ph1,
ph2,
z1,
integer nplat,
integer, dimension(*) iplat,
dhz,
k11,
k12,
k13,
k14,
k22,
k23,
k24,
k33,
k34,
k44,
m11,
m12,
m13,
m14,
m22,
m23,
m24,
m33,
m34,
m44,
mf11,
mf12,
mf13,
mf14,
mf22,
mf23,
mf24,
mf33,
mf34,
mf44,
fm12,
fm13,
fm14,
fm23,
fm24,
fm34,
phkrx,
phkry,
phkrxy,
pherx,
phery,
pherxy,
phkrz,
pherz,
phkx,
phky,
phex,
phey )

Definition at line 277 of file czlken3.F.

287C--------------------------------------------------------------------------------------------------
288C-----------------------------------------------
289C I M P L I C I T T Y P E S
290C-----------------------------------------------
291#include "implicit_f.inc"
292#include "mvsiz_p.inc"
293C-----------------------------------------------
294C D U M M Y A R G U M E N T S
295C-----------------------------------------------
296C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
297 INTEGER JFT,JLT,NPLAT,IPLAT(*)
298 my_real
299 . vol(*),thk0(*),thk2(*),a_i(*),z1(*),hm(mvsiz,4),hz(*),
300 . px1(*) ,px2(*) ,py1(*) ,py2(*) ,ph1(*) ,ph2(*) ,
301 . hxx(*),hyy(*),hxy(*),
302 . k11(3,3,*),k12(3,3,*),k13(3,3,*),k14(3,3,*),
303 . k22(3,3,*),k23(3,3,*),k24(3,3,*),k33(3,3,*),
304 . m11(3,3,*),m12(3,3,*),m13(3,3,*),m14(3,3,*),
305 . m22(3,3,*),m23(3,3,*),m24(3,3,*),m33(3,3,*),
306 . mf11(3,3,*),mf12(3,3,*),mf13(3,3,*),mf14(3,3,*),
307 . mf22(3,3,*),mf23(3,3,*),mf24(3,3,*),mf33(3,3,*),
308 . fm12(3,3,*),fm13(3,3,*),fm14(3,3,*),
309 . fm23(3,3,*),fm24(3,3,*),fm34(3,3,*),
310 . k34(3,3,*),k44(3,3,*),m34(3,3,*),m44(3,3,*),
311 . mf34(3,3,*),mf44(3,3,*),dhz(*)
312 my_real
313 . phkrx(4,*),phkry(4,*),
314 . phkrxy(4,*),pherx(4,*),phery(4,*),pherxy(4,*),
315 . phkrz(4,*),pherz(4,*),phkx(*) ,phky(*) ,phex(*) ,phey(*)
316C---------------|[KIJ][MFIJ]|----
317C-----KE(6x6)= | |
318C---------------|[FMIJ]{MIJ]|----
319C-----------------------------------------------
320C L O C A L V A R I A B L E S
321C-----------------------------------------------
322 INTEGER EP,I,J,NF,M
323 my_real
324 . dm(3,3,mvsiz),c1,c2,gm(mvsiz),
325 . g11(mvsiz),g12(mvsiz),g13(mvsiz),
326 . g14(mvsiz),g22(mvsiz),g23(mvsiz),
327 . g24(mvsiz),g33(mvsiz),g34(mvsiz),g44(mvsiz),
328 . chm(2,2,mvsiz),chf(2,2,mvsiz),facf(mvsiz),fac,
329 . gama(4,mvsiz),cxx,cyy,cxy,c11,c22,c12,rz_3,
330 . cbkrx(4,mvsiz),cbkry(4,mvsiz),cbkrz(4,mvsiz),
331 . cberx(4,mvsiz),cbery(4,mvsiz),cberz(4,mvsiz),
332 . grz3(4,mvsiz),zr(mvsiz),c3,cx1,cx2,cy1,cy2,zr2
333C-----------gama(I)=hI/4-PH(I), PH:anti-sym comme bxI------
334C-----------Attention Matrice sym Kii ne calcul que la moitie---------72
335#include "vectorize.inc"
336 DO m=jft,jlt
337 ep=iplat(m)
338 c2=vol(ep)
339 c1= third*c2
340 dm(1,1,m)=hm(ep,1)*c1
341 dm(2,2,m)=hm(ep,2)*c1
342 dm(1,2,m)=hm(ep,3)*c1
343c GM(M) =HM(EP,4)*C2
344 dm(1,3,m)=zero
345 dm(2,3,m)=zero
346 dm(3,3,m)=hm(ep,4)*c1
347 zr(m) = z1(ep)
348c zr(m)=zero
349 ENDDO
350 DO m=jft,jlt
351 dm(2,1,m)=dm(1,2,m)
352 dm(3,1,m)=dm(1,3,m)
353 dm(3,2,m)=dm(2,3,m)
354 ENDDO
355C------------------shear est non constante----
356C-------0.5*[-By Bx 0 0 0 BRZ]^tKG[-By Bx 0 0 0 BRZ]*0.5
357C-------DHZ=0.25*Kg*V---------
358 DO ep=jft,jlt
359 chm(1,1,ep)=dhz(ep)*hyy(ep)
360 chm(2,2,ep)=dhz(ep)*hxx(ep)
361 chm(1,2,ep)=-dhz(ep)*hxy(ep)
362 chm(2,1,ep)=chm(1,2,ep)
363 ENDDO
364C--------------gamaI--------
365 DO ep=jft,jlt
366 gama(1,ep)=fourth-ph1(ep)
367 gama(2,ep)=-fourth-ph2(ep)
368 gama(3,ep)=fourth+ph1(ep)
369 gama(4,ep)=-fourth+ph2(ep)
370 g11(ep) =gama(1,ep)*gama(1,ep)
371 g12(ep) =gama(1,ep)*gama(2,ep)
372 g13(ep) =gama(1,ep)*gama(3,ep)
373 g14(ep) =gama(1,ep)*gama(4,ep)
374 g22(ep) =gama(2,ep)*gama(2,ep)
375 g23(ep) =gama(2,ep)*gama(3,ep)
376 g24(ep) =gama(2,ep)*gama(4,ep)
377 g33(ep) =gama(3,ep)*gama(3,ep)
378 g34(ep) =gama(3,ep)*gama(4,ep)
379 g44(ep) =gama(4,ep)*gama(4,ep)
380 ENDDO
381C-------
382 DO i=1,2
383 DO j=i,2
384 DO ep=jft,jlt
385 k11(i,j,ep)=k11(i,j,ep)+chm(i,j,ep)*g11(ep)
386 k22(i,j,ep)=k22(i,j,ep)+chm(i,j,ep)*g22(ep)
387 k33(i,j,ep)=k33(i,j,ep)+chm(i,j,ep)*g33(ep)
388 k44(i,j,ep)=k44(i,j,ep)+chm(i,j,ep)*g44(ep)
389 ENDDO
390 ENDDO
391 ENDDO
392 DO i=1,2
393 DO j=1,2
394 DO ep=jft,jlt
395 k12(i,j,ep)=k12(i,j,ep)+chm(i,j,ep)*g12(ep)
396 k13(i,j,ep)=k13(i,j,ep)+chm(i,j,ep)*g13(ep)
397 k14(i,j,ep)=k14(i,j,ep)+chm(i,j,ep)*g14(ep)
398 k23(i,j,ep)=k23(i,j,ep)+chm(i,j,ep)*g23(ep)
399 k24(i,j,ep)=k24(i,j,ep)+chm(i,j,ep)*g24(ep)
400 k34(i,j,ep)=k34(i,j,ep)+chm(i,j,ep)*g34(ep)
401 ENDDO
402 ENDDO
403 ENDDO
404C
405 DO j=1,4
406 DO m=jft,jlt
407 rz_3= dhz(m)*third
408 grz3(j,m)=gama(j,m)*rz_3
409 ENDDO
410 ENDDO
411C
412 DO m=jft,jlt
413 i=1
414 j=1
415 c1 = -(phky(m)*phkrz(j,m)+phey(m)*pherz(j,m))
416 c2 = (phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m))
417 mf11(1,3,m)=mf11(1,3,m)+grz3(i,m)*c1
418 mf11(2,3,m)=mf11(2,3,m)+grz3(i,m)*c2
419 i=2
420 fm12(3,1,m)=fm12(3,1,m)+ grz3(i,m)*c1
421 fm12(3,2,m)=fm12(3,2,m)+ grz3(i,m)*c2
422 i=3
423 fm13(3,1,m)=fm13(3,1,m)+grz3(i,m)*c1
424 fm13(3,2,m)=fm13(3,2,m)+grz3(i,m)*c2
425 i=4
426 fm14(3,1,m)=fm14(3,1,m)+ grz3(i,m)*c1
427 fm14(3,2,m)=fm14(3,2,m)+ grz3(i,m)*c2
428 j=2
429 c1 = -(phky(m)*phkrz(j,m)+phey(m)*pherz(j,m))
430 c2 = (phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m))
431 i=1
432 mf12(1,3,m)=mf12(1,3,m)+grz3(i,m)*c1
433 mf12(2,3,m)=mf12(2,3,m)+grz3(i,m)*c2
434 i=2
435 mf22(1,3,m)=mf22(1,3,m)+grz3(i,m)*c1
436 mf22(2,3,m)=mf22(2,3,m)+grz3(i,m)*c2
437 i=3
438 fm23(3,1,m)=fm23(3,1,m)+grz3(i,m)*c1
439 fm23(3,2,m)=fm23(3,2,m)+grz3(i,m)*c2
440 i=4
441 fm24(3,1,m)=fm24(3,1,m)+ grz3(i,m)*c1
442 fm24(3,2,m)=fm24(3,2,m)+ grz3(i,m)*c2
443 j=3
444 c1 = -(phky(m)*phkrz(j,m)+phey(m)*pherz(j,m))
445 c2 = (phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m))
446 i=1
447 mf13(1,3,m)=mf13(1,3,m)+ grz3(i,m)*c1
448 mf13(2,3,m)=mf13(2,3,m)+ grz3(i,m)*c2
449 i=2
450 mf23(1,3,m)=mf23(1,3,m)+grz3(i,m)*c1
451 mf23(2,3,m)=mf23(2,3,m)+grz3(i,m)*c2
452 i=3
453 mf33(1,3,m)=mf33(1,3,m)+ grz3(i,m)*c1
454 mf33(2,3,m)=mf33(2,3,m)+ grz3(i,m)*c2
455 i=4
456 fm34(3,1,m)=fm34(3,1,m)+ grz3(i,m)*c1
457 fm34(3,2,m)=fm34(3,2,m)+ grz3(i,m)*c2
458 j=4
459 c1 = -(phky(m)*phkrz(j,m)+phey(m)*pherz(j,m))
460 c2 = (phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m))
461 i=1
462 mf14(1,3,m)=mf14(1,3,m)+ grz3(i,m)*c1
463 mf14(2,3,m)=mf14(2,3,m)+ grz3(i,m)*c2
464 i=2
465 mf24(1,3,m)=mf24(1,3,m)+ grz3(i,m)*c1
466 mf24(2,3,m)=mf24(2,3,m)+ grz3(i,m)*c2
467 i=3
468 mf34(1,3,m)=mf34(1,3,m)+ grz3(i,m)*c1
469 mf34(2,3,m)=mf34(2,3,m)+ grz3(i,m)*c2
470 i=4
471 mf44(1,3,m)=mf44(1,3,m)+ grz3(i,m)*c1
472 mf44(2,3,m)=mf44(2,3,m)+ grz3(i,m)*c2
473 ENDDO
474C------ Mzz ----------
475 DO ep=jft,jlt
476 rz_3= dhz(ep)*third
477 i=1
478 j=1
479 m11(3,3,ep)=m11(3,3,ep)+
480 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
481 j=2
482 m12(3,3,ep)=m12(3,3,ep)+
483 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
484 j=3
485 m13(3,3,ep)=m13(3,3,ep)+
486 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
487 j=4
488 m14(3,3,ep)=m14(3,3,ep)+
489 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
490 i=2
491 j=2
492 m22(3,3,ep)=m22(3,3,ep)+
493 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
494 j=3
495 m23(3,3,ep)=m23(3,3,ep)+
496 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
497 j=4
498 m24(3,3,ep)=m24(3,3,ep)+
499 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
500 i=3
501 j=3
502 m33(3,3,ep)=m33(3,3,ep)+
503 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
504 j=4
505 m34(3,3,ep)=m34(3,3,ep)+
506 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
507 i=4
508 m44(3,3,ep)=m44(3,3,ep)+
509 . rz_3*(phkrz(i,ep)*phkrz(j,ep)+pherz(i,ep)*pherz(j,ep))
510 ENDDO
511C-------[MFIJ]=[Bm]^t[C][BRm]; [MIJ]=[BRm]^t[C][BRm];----
512C---------------|0 0 BRX |----
513C-----BR(3x6)= |0 0 BRY |
514C---------------|0 0 BRXY+BRYX |----
515C---------------Ksi first---
516 DO j=1,4
517 DO m=jft,jlt
518 cbkrx(j,m) =dm(1,1,m)*phkrx(j,m)+dm(1,2,m)*phkry(j,m)+
519 . dm(1,3,m)*phkrxy(j,m)
520 cbkry(j,m) =dm(2,1,m)*phkrx(j,m)+dm(2,2,m)*phkry(j,m)+
521 . dm(2,3,m)*phkrxy(j,m)
522 cbkrz(j,m) =dm(3,1,m)*phkrx(j,m)+dm(3,2,m)*phkry(j,m)+
523 . dm(3,3,m)*phkrxy(j,m)
524 cberx(j,m) =dm(1,1,m)*pherx(j,m)+dm(1,2,m)*phery(j,m)+
525 . dm(1,3,m)*pherxy(j,m)
526 cbery(j,m) =dm(2,1,m)*pherx(j,m)+dm(2,2,m)*phery(j,m)+
527 . dm(2,3,m)*pherxy(j,m)
528 cberz(j,m) =dm(3,1,m)*pherx(j,m)+dm(3,2,m)*phery(j,m)+
529 . dm(3,3,m)*pherxy(j,m)
530 ENDDO
531 ENDDO
532C
533 DO m=jft,jlt
534 i=1
535 j=1
536 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
537 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
538 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
539 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
540 mf11(1,3,m)=mf11(1,3,m)+gama(i,m)*c1
541 mf11(2,3,m)=mf11(2,3,m)+gama(i,m)*c2
542 i=2
543 fm12(3,1,m)=fm12(3,1,m)+ gama(i,m)*c1
544 fm12(3,2,m)=fm12(3,2,m)+ gama(i,m)*c2
545 i=3
546 fm13(3,1,m)=fm13(3,1,m)+gama(i,m)*c1
547 fm13(3,2,m)=fm13(3,2,m)+gama(i,m)*c2
548 i=4
549 fm14(3,1,m)=fm14(3,1,m)+ gama(i,m)*c1
550 fm14(3,2,m)=fm14(3,2,m)+ gama(i,m)*c2
551 j=2
552 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
553 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
554 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
555 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
556 i=1
557 mf12(1,3,m)=mf12(1,3,m)+gama(i,m)*c1
558 mf12(2,3,m)=mf12(2,3,m)+gama(i,m)*c2
559 i=2
560 mf22(1,3,m)=mf22(1,3,m)+gama(i,m)*c1
561 mf22(2,3,m)=mf22(2,3,m)+gama(i,m)*c2
562 i=3
563 fm23(3,1,m)=fm23(3,1,m)+ gama(i,m)*c1
564 fm23(3,2,m)=fm23(3,2,m)+ gama(i,m)*c2
565 i=4
566 fm24(3,1,m)=fm24(3,1,m)+ gama(i,m)*c1
567 fm24(3,2,m)=fm24(3,2,m)+ gama(i,m)*c2
568 j=3
569 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
570 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
571 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
572 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
573 i=1
574 mf13(1,3,m)=mf13(1,3,m)+ gama(i,m)*c1
575 mf13(2,3,m)=mf13(2,3,m)+ gama(i,m)*c2
576 i=2
577 mf23(1,3,m)=mf23(1,3,m)+gama(i,m)*c1
578 mf23(2,3,m)=mf23(2,3,m)+gama(i,m)*c2
579 i=3
580 mf33(1,3,m)=mf33(1,3,m)+ gama(i,m)*c1
581 mf33(2,3,m)=mf33(2,3,m)+ gama(i,m)*c2
582 i=4
583 fm34(3,1,m)=fm34(3,1,m)+ gama(i,m)*c1
584 fm34(3,2,m)=fm34(3,2,m)+ gama(i,m)*c2
585 j=4
586 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
587 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
588 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
589 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
590 i=1
591 mf14(1,3,m)=mf14(1,3,m)+ gama(i,m)*c1
592 mf14(2,3,m)=mf14(2,3,m)+ gama(i,m)*c2
593 i=2
594 mf24(1,3,m)=mf24(1,3,m)+ gama(i,m)*c1
595 mf24(2,3,m)=mf24(2,3,m)+ gama(i,m)*c2
596 i=3
597 mf34(1,3,m)=mf34(1,3,m)+ gama(i,m)*c1
598 mf34(2,3,m)=mf34(2,3,m)+ gama(i,m)*c2
599 i=4
600 mf44(1,3,m)=mf44(1,3,m)+ gama(i,m)*c1
601 mf44(2,3,m)=mf44(2,3,m)+ gama(i,m)*c2
602 ENDDO
603C------ Mzz ----------
604 DO m=jft,jlt
605 i=1
606 j=1
607 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
608 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
609 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
610 m11(3,3,m)=m11(3,3,m)+c1
611 j=2
612 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
613 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
614 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
615 m12(3,3,m)=m12(3,3,m)+c1
616 j=3
617 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
618 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
619 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
620 m13(3,3,m)=m13(3,3,m)+c1
621 j=4
622 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
623 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
624 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
625 m14(3,3,m)=m14(3,3,m)+c1
626 i=2
627 j=2
628 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
629 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
630 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
631 m22(3,3,m)=m22(3,3,m)+c1
632 j=3
633 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
634 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
635 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
636 m23(3,3,m)=m23(3,3,m)+c1
637 j=4
638 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
639 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
640 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
641 m24(3,3,m)=m24(3,3,m)+c1
642 i=3
643 j=3
644 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
645 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
646 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
647 m33(3,3,m)=m33(3,3,m)+c1
648 j=4
649 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
650 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
651 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
652 m34(3,3,m)=m34(3,3,m)+c1
653 i=4
654 c1 = phkrx(i,m)*cbkrx(j,m)+pherx(i,m)*cberx(j,m)
655 . +phkry(i,m)*cbkry(j,m)+phery(i,m)*cbery(j,m)
656 . +phkrxy(i,m)*cbkrz(j,m)+pherxy(i,m)*cberz(j,m)
657 m44(3,3,m)=m44(3,3,m)+c1
658 ENDDO
659C---+---------+---------+warped elements-------------
660 nf=nplat+1
661 DO m=nf,jlt
662 j=1
663 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
664 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
665 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
666 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
667C I=1,3
668 c3 = zr(m)*(px1(m)*c1+py1(m)*c2)
669 mf11(3,3,m)=mf11(3,3,m)+ c3
670 fm13(3,3,m)=fm13(3,3,m)- c3
671C I=2,4
672 c3 = zr(m)*(px2(m)*c1+py2(m)*c2)
673 fm12(3,3,m)=fm12(3,3,m)+ c3
674 fm14(3,3,m)=fm14(3,3,m)- c3
675 j=2
676 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
677 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
678 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
679 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
680C I=1,3
681 c3 = zr(m)*(px1(m)*c1+py1(m)*c2)
682 mf12(3,3,m)=mf12(3,3,m)+ c3
683 fm23(3,3,m)=fm23(3,3,m)- c3
684C I=2,4
685 c3 = zr(m)*(px2(m)*c1+py2(m)*c2)
686 mf22(3,3,m)=mf22(3,3,m)+ c3
687 fm24(3,3,m)=fm24(3,3,m)- c3
688 j=3
689 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
690 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
691 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
692 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
693C I=1,3
694 c3 = zr(m)*(px1(m)*c1+py1(m)*c2)
695 mf13(3,3,m)=mf13(3,3,m)+ c3
696 mf33(3,3,m)=mf33(3,3,m)- c3
697C I=2,4
698 c3 = zr(m)*(px2(m)*c1+py2(m)*c2)
699 mf23(3,3,m)=mf23(3,3,m)+ c3
700 fm34(3,3,m)=fm34(3,3,m)- c3
701 j=4
702 c1 = (phkx(m)*cbkrx(j,m)+phex(m)*cberx(j,m)+
703 . phky(m)*cbkrz(j,m)+phey(m)*cberz(j,m))
704 c2 = (phky(m)*cbkry(j,m)+phey(m)*cbery(j,m)+
705 . phkx(m)*cbkrz(j,m)+phex(m)*cberz(j,m))
706C I=1,3
707 c3 = zr(m)*(px1(m)*c1+py1(m)*c2)
708 mf14(3,3,m)=mf14(3,3,m)+ c3
709 mf34(3,3,m)=mf34(3,3,m)- c3
710C I=2,4
711 c3 = zr(m)*(px2(m)*c1+py2(m)*c2)
712 mf24(3,3,m)=mf24(3,3,m)+ c3
713 mf44(3,3,m)=mf44(3,3,m)- c3
714 ENDDO
715C------------------add coulpling terms w/ warped----
716C-------0.5*[-By Bx Zr* 0 0 BRZ]^tKG[-By Bx Zr* 0 0 BRZ]*0.5
717C-------DHZ=0.25*KG*V---------
718C-----KIJ(3,1)=KJI(1,3)=DHZ*Zr*(HYY*BxI*GamaJ-HXY*ByI*GamaJ)----------
719C-----KIJ(3,2)=KJI(2,3)=DHZ*Zr*(HXX*ByI*GamaJ-HXY*BxI*GamaJ)----------
720 DO m=nf,jlt
721C-------------------------
722 j=1
723C-----------------------
724 c11 = zr(m)*gama(j,m)
725C I=1,3
726 c1 = (chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))*c11
727 c2 = (chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))*c11
728 k11(1,3,m)=k11(1,3,m)+ c1
729 k13(1,3,m)=k13(1,3,m)- c1
730 k11(2,3,m)=k11(2,3,m)+ c2
731 k13(2,3,m)=k13(2,3,m)- c2
732C I=2,4
733 c1 = (chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))*c11
734 c2 = (chm(1,2,m)*px2(m)+chm(2,2,m)*py2(m))*c11
735 k12(1,3,m)=k12(1,3,m)+ c1
736 k14(1,3,m)=k14(1,3,m)- c1
737 k12(2,3,m)=k12(2,3,m)+ c2
738 k14(2,3,m)=k14(2,3,m)- c2
739C-------------------------
740 j=2
741C---------------------------
742 c11 = zr(m)*gama(j,m)
743C I=1,3
744 c1 = (chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))*c11
745 c2 = (chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))*c11
746 k12(3,1,m)=k12(3,1,m)+ c1
747 k23(1,3,m)=k23(1,3,m)- c1
748 k12(3,2,m)=k12(3,2,m)+ c2
749 k23(2,3,m)=k23(2,3,m)- c2
750C I=2,4
751 c1 = (chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))*c11
752 c2 = (chm(1,2,m)*px2(m)+chm(2,2,m)*py2(m))*c11
753 k22(1,3,m)=k22(1,3,m)+ c1
754 k24(1,3,m)=k24(1,3,m)- c1
755 k22(2,3,m)=k22(2,3,m)+ c2
756 k24(2,3,m)=k24(2,3,m)- c2
757C-------------------------
758 j=3
759C-------------------------
760 c11 = zr(m)*gama(j,m)
761C I=1,3
762 c1 = (chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))*c11
763 c2 = (chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))*c11
764 k13(3,1,m)=k13(3,1,m)+ c1
765 k33(1,3,m)=k33(1,3,m)- c1
766 k13(3,2,m)=k13(3,2,m)+ c2
767 k33(2,3,m)=k33(2,3,m)- c2
768C I=2,4
769 c1 = (chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))*c11
770 c2 = (chm(1,2,m)*px2(m)+chm(2,2,m)*py2(m))*c11
771 k23(3,1,m)=k23(3,1,m)+ c1
772 k34(1,3,m)=k34(1,3,m)- c1
773 k23(3,2,m)=k23(3,2,m)+ c2
774 k34(2,3,m)=k34(2,3,m)- c2
775C-------------------------
776 j=4
777C-------------------------
778 c11 = zr(m)*gama(j,m)
779C I=1,3
780 c1 = (chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))*c11
781 c2 = (chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))*c11
782 k14(3,1,m)=k14(3,1,m)+ c1
783 k34(3,1,m)=k34(3,1,m)- c1
784 k14(3,2,m)=k14(3,2,m)+ c2
785 k34(3,2,m)=k34(3,2,m)- c2
786C I=2,4
787 c1 = (chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))*c11
788 c2 = (chm(1,2,m)*px2(m)+chm(2,2,m)*py2(m))*c11
789 k24(3,1,m)=k24(3,1,m)+ c1
790 k44(1,3,m)=k44(1,3,m)- c1
791 k24(3,2,m)=k24(3,2,m)+ c2
792 k44(2,3,m)=k44(2,3,m)- c2
793 ENDDO
794C-----KIJ(3,3)=DHZ*Zr*Zr*(HYY*BxI*BxJ+HXX*ByI*ByJ-HXY*BxI*ByJ-HXY*ByI*BxJ)----------
795 DO m=nf,jlt
796 zr2= zr(m)*zr(m)
797C-------------------------
798 j=1
799C-----------------------
800 cx1=zr2*(chm(1,1,m)*px1(m)+chm(1,2,m)*py1(m))
801 cy1=zr2*(chm(1,2,m)*px1(m)+chm(2,2,m)*py1(m))
802C I=1,3
803 c3=px1(m)*cx1+py1(m)*cy1
804 k11(3,3,m)=k11(3,3,m)+ c3
805c J=3
806 k13(3,3,m)=k13(3,3,m)- c3
807 k33(3,3,m)=k33(3,3,m)+ c3
808C I=2,4
809 c3=px2(m)*cx1+py2(m)*cy1
810c J=3
811 k23(3,3,m)=k23(3,3,m)- c3
812C-------------------------
813 j=2
814C---------------------------
815 cx1=zr2*(chm(1,1,m)*px2(m)+chm(1,2,m)*py2(m))
816 cy1=zr2*(chm(1,2,m)*px2(m)+chm(2,2,m)*py2(m))
817C I=1,3
818 c3=px1(m)*cx1+py1(m)*cy1
819 k12(3,3,m)=k12(3,3,m)+ c3
820c J=4
821 k14(3,3,m)=k14(3,3,m)- c3
822 k34(3,3,m)=k34(3,3,m)+ c3
823C I=2,4
824 c3=px2(m)*cx1+py2(m)*cy1
825 k22(3,3,m)=k22(3,3,m)+ c3
826c J=4
827 k24(3,3,m)=k24(3,3,m)- c3
828 k44(3,3,m)=k44(3,3,m)+ c3
829 ENDDO
830C-----MFIJ(3,3)=FMJI(3,3)=DHZ*Zr(-PHy*BxI*BRZJ+PHx*ByI*BRZJ)----------
831 DO m=nf,jlt
832 rz_3= zr(m)*dhz(m)*third
833C-------------------------
834 j=1
835C-----------------------
836 cx2 = phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m)
837 cy2 = phky(m)*phkrz(j,m)+phey(m)*pherz(j,m)
838C I=1,3
839 c3 = rz_3*(-px1(m)*cy2+py1(m)*cx2)
840 mf11(3,3,m)=mf11(3,3,m)+ c3
841 fm13(3,3,m)=fm13(3,3,m)- c3
842C I=2,4
843 c3 = rz_3*(-px2(m)*cy2+py2(m)*cx2)
844 fm12(3,3,m)=fm12(3,3,m)+ c3
845 fm14(3,3,m)=fm14(3,3,m)- c3
846C-------------------------
847 j=2
848C---------------------------
849 cx2 = phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m)
850 cy2 = phky(m)*phkrz(j,m)+phey(m)*pherz(j,m)
851C I=1,3
852 c3 = rz_3*(-px1(m)*cy2+py1(m)*cx2)
853 mf12(3,3,m)=mf12(3,3,m)+ c3
854 fm23(3,3,m)=fm23(3,3,m)- c3
855C I=2,4
856 c3 = rz_3*(-px2(m)*cy2+py2(m)*cx2)
857 mf22(3,3,m)=mf22(3,3,m)+ c3
858 fm24(3,3,m)=fm24(3,3,m)- c3
859C-------------------------
860 j=3
861C-------------------------
862 cx2 = phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m)
863 cy2 = phky(m)*phkrz(j,m)+phey(m)*pherz(j,m)
864C I=1,3
865 c3 = rz_3*(-px1(m)*cy2+py1(m)*cx2)
866 mf13(3,3,m)=mf13(3,3,m)+ c3
867 mf33(3,3,m)=mf33(3,3,m)- c3
868C I=2,4
869 c3 = rz_3*(-px2(m)*cy2+py2(m)*cx2)
870 mf23(3,3,m)=mf23(3,3,m)+ c3
871 fm34(3,3,m)=fm34(3,3,m)- c3
872C-------------------------
873 j=4
874C-------------------------
875 cx2 = phkx(m)*phkrz(j,m)+phex(m)*pherz(j,m)
876 cy2 = phky(m)*phkrz(j,m)+phey(m)*pherz(j,m)
877C I=1,3
878 c3 = rz_3*(-px1(m)*cy2+py1(m)*cx2)
879 mf14(3,3,m)=mf14(3,3,m)+ c3
880 mf34(3,3,m)=mf34(3,3,m)- c3
881C I=2,4
882 c3 = rz_3*(-px2(m)*cy2+py2(m)*cx2)
883 mf24(3,3,m)=mf24(3,3,m)+ c3
884 mf44(3,3,m)=mf44(3,3,m)- c3
885 ENDDO
886C
887 RETURN