OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cbacoor.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "scr17_c.inc"
#include "com08_c.inc"
#include "impl1_c.inc"
#include "scr14_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine cbacoor (elbuf_str, jft, jlt, x, v, vr, ixc, pm, offg, ll, area, vxyz, rxyz, vcore, jac, hx, hy, vksi, veta, vqn, vqg, vq, vjfi, vnrm, vastn, nplat, iplat, x13_t, x24_t, y13_t, y24_t, off, di, nlay, irep, npt, ismstr, nel, isrot, smstr, dir_a, dir_b, facn, zl1, r11, r12, r13, r21, r22, r23, r31, r32, r33, inod, rlz, thk, iplycxfem, ux1, ux2, ux3, ux4, uy1, uy2, uy3, uy4, vl1, vl2, vl3, vl4, xl2, xl3, xl4, yl2, yl3, yl4, xlcor, npinch)
subroutine cbacoort (elbuf_str, jft, jlt, x, v, vr, dr, ixc, pm, offg, area, vxyz, rlz, vcore, jac, hx, hy, vq, vqg, vjfi, nplat, iplat, x13_t, x24_t, y13_t, y24_t, npt, smstr, isrot, xlcor, zl, vqn, nel)

Function/Subroutine Documentation

◆ cbacoor()

subroutine cbacoor ( type(elbuf_struct_) elbuf_str,
integer jft,
integer jlt,
x,
v,
vr,
integer, dimension(nixc,*) ixc,
pm,
offg,
ll,
area,
vxyz,
rxyz,
vcore,
jac,
hx,
hy,
vksi,
veta,
vqn,
vqg,
vq,
vjfi,
vnrm,
vastn,
integer nplat,
integer, dimension(*) iplat,
x13_t,
x24_t,
y13_t,
y24_t,
off,
di,
integer nlay,
integer irep,
integer npt,
integer ismstr,
integer nel,
integer isrot,
double precision, dimension(*) smstr,
dir_a,
dir_b,
facn,
zl1,
r11,
r12,
r13,
r21,
r22,
r23,
r31,
r32,
r33,
integer, dimension(*) inod,
rlz,
thk,
integer iplycxfem,
ux1,
ux2,
ux3,
ux4,
uy1,
uy2,
uy3,
uy4,
vl1,
vl2,
vl3,
vl4,
xl2,
xl3,
xl4,
yl2,
yl3,
yl4,
xlcor,
integer npinch )

Definition at line 34 of file cbacoor.F.

47C-----------------------------------------------
48C M o d u l e s
49C-----------------------------------------------
50 USE elbufdef_mod
51 USE plyxfem_mod
52C-----------------------------------------------
53C I M P L I C I T T Y P E S
54C-----------------------------------------------
55#include "implicit_f.inc"
56#include "comlock.inc"
57c-----------------------------------------------
58c g l o b a l p a r a m e t e r s
59c-----------------------------------------------
60#include "mvsiz_p.inc"
61#include "param_c.inc"
62#include "scr17_c.inc"
63#include "com08_c.inc"
64#include "impl1_c.inc"
65#include "scr14_c.inc"
66C-----------------------------------------------
67C D U M M Y A R G U M E N T S
68C-----------------------------------------------
69 INTEGER JFT,JLT,NNOD,NPLAT,IREP,NPT,NLAY,ISMSTR,ISROT,IPLYCXFEM,NEL,NPINCH
70 INTEGER IXC(NIXC,*),IPLAT(*),INOD(*)
71 parameter(nnod = 4)
72 my_real x(3,*), v(3,*), vr(3,*),pm(npropm,*),offg(*),off(*),
73 . vcore(mvsiz,3,nnod),vxyz(mvsiz,3,nnod),area(*),vjfi(mvsiz,6,4),
74 . vqn(mvsiz,9,nnod),vqg(mvsiz,9,nnod),vq(mvsiz,3,3),rxyz(mvsiz,2,nnod),
75 . ll(*),vastn(mvsiz,4,nnod),vnrm(mvsiz,3,nnod),
76 . jac(mvsiz,4),hx(mvsiz,4),hy(mvsiz,4),veta(4,4),vksi(4,4),y24_t(*),
77 . x13_t(*),x24_t(*),y13_t(*),off_l, di(mvsiz,6),
78 . dir_a(nel,*),dir_b(nel,*),facn(mvsiz,2),
79 . zl1(mvsiz),r11(mvsiz),r12(mvsiz),r13(mvsiz),
80 . r21(mvsiz),r22(mvsiz),r23(mvsiz),r31(mvsiz),r32(mvsiz),
81 . r33(mvsiz),rlz(mvsiz,4),thk(mvsiz),
82 . ux1(*),ux2(*),ux3(*),ux4(*),uy1(*),uy2(*),uy3(*),uy4(*),
83 . vl1(mvsiz,3),vl2(mvsiz,3),vl3(mvsiz,3),vl4(mvsiz,3),
84 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),yl3(mvsiz),yl4(mvsiz),
85 . xlcor(mvsiz,2,3)
86 double precision
87 . smstr(*)
88C
89 TYPE(ELBUF_STRUCT_) :: ELBUF_STR
90C-----------------------------------------------
91c FUNCTION: premilitary compute for QBAT
92c
93c Note:
94c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
95c
96c TYPE NAME FUNCTION
97c I JFT,JLT - element id limit
98c I X, V, VR(3,NUMNOD) - coordinate,velocity, rotational velocity in global system
99c I IXC(NIXC,*NEL) - element connectivity of shell (and other data)
100c I PM(NPROPM,MID) - input Material data
101c I OFFG(NEL),OFF(NEL) - element activation flag, local flag
102c O LL(NEL),FACN(2,NEL)- characteristic length (for dt compute)
103c O AREA(NEL),NTP - Area, num. of integrating points in thickness
104c O VCORE(MVSIZ,3,4) - coordinates of 4 nodes in local system
105c BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
106c for plat element
107c O VXYZ(3,4,NEL) - nodal velocity (4 nodes) in local system
108C V13->VXYZ(,1,),V24->VXYZ(,2,),VHI->VXYZ(,3,)------------
109c for plat element
110c O RXYZ(NEL,2,4) - nodal rotational velocity (4 nodes) in local system 5 dof
111C R13->RXYZ(,1,),R24->RXYZ(,2,),RHI->RXYZ(,3,),RTI->RXYZ(,4,)----
112c for plat element
113c O JAC,HX,HY(4,NEL) - Jacobien,asymmetric part(hourglass) at 4 Gauss points
114c O VKSI,VETA(4,4) - derivation of shape functions------
115c O VQ(NEL,9),VQN,VQG(NEL,9,4) - local system of element, nodal(4) and Gauss points(4)
116c O VJFI,VNRM,VASTN - terms used for assumed strains
117c O NPLAT,IPLAT(NEL) - num. of plat element and indice
118c O Xij_T,Yij_T(NEL) - [B] components used for in-plane shear assumed strain
119c IO XiS,YiS,ZiS(NEL) - reference configurations used for small strain options
120c I ISMSTR,IREP - small strain and orthotrop flags
121c O DIR_A,DIR_B(2,NEL) - actual orthotropic directions
122c O ZL1(NEL) - warped length
123c O Rij(NEL) - local system base
124c I IDRIL - drilling dof flag
125c O RLZ(4,NEL) - nodal Rz (rotational velocity) in local system
126c (used only Idril active )
127c I INOD(NUMNOD) - Node numbers for PLYXFEM
128C-----------------------------------------------
129C L O C A L V A R I A B L E S
130C-----------------------------------------------
131 INTEGER I,II(6),J,K,L,M,I1,I2,I3,I4,EP,SPLAT,JJ
132 INTEGER NN(4),JPLAT(MVSIZ)
133 my_real
134 . xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4
135 my_real
136 . c1,c2,l13(mvsiz),l24(mvsiz),
137 . aa,vv(3,nnod),rr(3,nnod)
138 my_real
139 . tol,pg,j0,j1,j2,deta,corel(3,4),x1,y1,s1
140 my_real
141 . x13,x24,y13,mx13,mx23,mx34,my13,y13_2,z1,z2,gama1,gama2,
142 . x21,x34,y21,y34,z21,z34,x41,x32,y41,y32,z_2,z41,z32,l12,l34,
143 . a_4,sl,sz2,sz,jmx13,jmy13,jmz13,j2myz,lm(mvsiz),y24,y24_2,
144 . my23,my34,scal,g1x1,g1x3,g1y1,g1y3,g2x1,g2x2,g2y1,g2y2,
145 . ar(3),ad(nnod),btb(6),xx1,yy,zz,xy,xz1,yz,d(6),
146 . alr(3),ald(nnod),dbad(3),abc,xxyz2,zzxy2,yyxz2
147 my_real
148 . lxyz0(3),deta1(mvsiz),rx(mvsiz), ry(mvsiz), rz(mvsiz),
149 . sx(mvsiz),sy(mvsiz),
150 . xx,area_i(mvsiz),ssz(mvsiz)
151C
152 my_real
153 . vg13(3),vg24(3),vghi(3),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
154 . faci,fac1,fac2,
155 . dt05,dt025,exz,eyz,ddrx,ddry,v13x,v24x,vhix,ddrz1,ddrz2,ddrz,
156 . vnx1, vny1, vnz1,vnx2,vny2,vnz2,vnx3,vny3,vnz3,vnx4,vny4,vnz4,
157 . norm
158 DATA pg/.577350269189626/
159C=======================================================================
160 DO i=1,6
161 ii(i) = nel*(i-1)
162 ENDDO
163C
164 tol=em12
165 IF (isrot > 0 ) tol=em8
166C
167 vksi(1,1)=-fourth*(one+pg)
168 vksi(2,1)=-vksi(1,1)
169 vksi(3,1)= fourth*(one-pg)
170 vksi(4,1)=-vksi(3,1)
171 veta(1,1)=-fourth*(one+pg)
172 veta(2,1)=-fourth*(one-pg)
173 veta(3,1)=-veta(2,1)
174 veta(4,1)=-veta(1,1)
175 vksi(1,2)= vksi(1,1)
176 vksi(2,2)=-vksi(1,2)
177 vksi(3,2)= vksi(3,1)
178 vksi(4,2)=-vksi(3,2)
179 veta(1,2)= veta(2,1)
180 veta(2,2)= veta(1,1)
181 veta(3,2)=-veta(2,2)
182 veta(4,2)=-veta(1,2)
183 vksi(1,3)=-vksi(3,1)
184 vksi(2,3)=-vksi(1,3)
185 vksi(3,3)=-vksi(1,1)
186 vksi(4,3)=-vksi(3,3)
187 veta(1,3)= veta(1,2)
188 veta(2,3)= veta(2,2)
189 veta(3,3)=-veta(2,3)
190 veta(4,3)=-veta(1,3)
191 vksi(1,4)= vksi(1,3)
192 vksi(2,4)=-vksi(1,4)
193 vksi(3,4)= vksi(3,3)
194 vksi(4,4)=-vksi(3,4)
195 veta(1,4)= veta(1,1)
196 veta(2,4)= veta(2,1)
197 veta(3,4)=-veta(2,4)
198 veta(4,4)=-veta(1,4)
199C
200 DO i=jft,jlt
201 rx(i)=x(1,ixc(3,i))+x(1,ixc(4,i))-x(1,ixc(2,i))-x(1,ixc(5,i))
202 sx(i)=x(1,ixc(4,i))+x(1,ixc(5,i))-x(1,ixc(2,i))-x(1,ixc(3,i))
203 ry(i)=x(2,ixc(3,i))+x(2,ixc(4,i))-x(2,ixc(2,i))-x(2,ixc(5,i))
204 sy(i)=x(2,ixc(4,i))+x(2,ixc(5,i))-x(2,ixc(2,i))-x(2,ixc(3,i))
205 rz(i)=x(3,ixc(3,i))+x(3,ixc(4,i))-x(3,ixc(2,i))-x(3,ixc(5,i))
206 ssz(i)=x(3,ixc(4,i))+x(3,ixc(5,i))-x(3,ixc(2,i))-x(3,ixc(3,i))
207 ENDDO
208C----------------------------
209C LOCAL SYSTEM
210C----------------------------
211 k = 0
212 CALL clskew3(jft,jlt,k,
213 . rx, ry, rz,
214 . sx, sy, ssz,
215 . r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
216 DO i=jft,jlt
217 area(i)=fourth*deta1(i)
218 area_i(i)=one/area(i)
219 vq(i,1,1)=r11(i)
220 vq(i,2,1)=r21(i)
221 vq(i,3,1)=r31(i)
222 vq(i,1,2)=r12(i)
223 vq(i,2,2)=r22(i)
224 vq(i,3,2)=r32(i)
225 vq(i,1,3)=r13(i)
226 vq(i,2,3)=r23(i)
227 vq(i,3,3)=r33(i)
228 ENDDO
229c
230 DO i=jft,jlt
231
232 j=ixc(2,i)
233 k=ixc(3,i)
234 l=ixc(4,i)
235 m=ixc(5,i)
236
237 lxyz0(1)=fourth*(x(1,l)+x(1,m)+x(1,j)+x(1,k))
238 lxyz0(2)=fourth*(x(2,l)+x(2,m)+x(2,j)+x(2,k))
239 lxyz0(3)=fourth*(x(3,l)+x(3,m)+x(3,j)+x(3,k))
240
241 xx1=x(1,k)-x(1,j)
242 yy1=x(2,k)-x(2,j)
243 zz1=x(3,k)-x(3,j)
244
245 xl2(i)=r11(i)*xx1+r21(i)*yy1+r31(i)*zz1
246 yl2(i)=r12(i)*xx1+r22(i)*yy1+r32(i)*zz1
247
248 xx2=x(1,j)-lxyz0(1)
249 yy2=x(2,j)-lxyz0(2)
250 zz2=x(3,j)-lxyz0(3)
251 zl1(i)=r13(i)*xx2+r23(i)*yy2+r33(i)*zz2
252
253 xx3=x(1,l)-x(1,j)
254 yy3=x(2,l)-x(2,j)
255 zz3=x(3,l)-x(3,j)
256 xl3(i)=r11(i)*xx3+r21(i)*yy3+r31(i)*zz3
257 yl3(i)=r12(i)*xx3+r22(i)*yy3+r32(i)*zz3
258
259 xx4=x(1,m)-x(1,j)
260 yy4=x(2,m)-x(2,j)
261 zz4=x(3,m)-x(3,j)
262 xl4(i)=r11(i)*xx4+r21(i)*yy4+r31(i)*zz4
263 yl4(i)=r12(i)*xx4+r22(i)*yy4+r32(i)*zz4
264 ENDDO
265C----------------------------
266C SMALL STRAIN
267C----------------------------
268 IF(ismstr==10)THEN
269 DO i=jft,jlt
270 xlcor(i,1,1)=xl2(i)
271 xlcor(i,2,1)=yl2(i)
272 xlcor(i,1,2)=xl3(i)
273 xlcor(i,2,2)=yl3(i)
274 xlcor(i,1,3)=xl4(i)
275 xlcor(i,2,3)=yl4(i)
276 ENDDO
277 ELSEIF(ismstr==11)THEN
278 DO i=jft,jlt
279 IF(abs(offg(i))==one)offg(i)=sign(two,offg(i))
280 ux1(i) = zero
281 uy1(i) = zero
282 ux2(i) = zero
283 uy2(i) = zero
284 ux3(i) = zero
285 uy3(i) = zero
286 ux4(i) = zero
287 uy4(i) = zero
288 IF(abs(offg(i))==two)THEN
289 ux2(i) = xl2(i)-smstr(ii(1)+i)
290 uy2(i) = yl2(i)-smstr(ii(2)+i)
291 ux3(i) = xl3(i)-smstr(ii(3)+i)
292 uy3(i) = yl3(i)-smstr(ii(4)+i)
293 ux4(i) = xl4(i)-smstr(ii(5)+i)
294 uy4(i) = yl4(i)-smstr(ii(6)+i)
295 xl2(i) = smstr(ii(1)+i)
296 yl2(i) = smstr(ii(2)+i)
297 xl3(i) = smstr(ii(3)+i)
298 yl3(i) = smstr(ii(4)+i)
299 xl4(i) = smstr(ii(5)+i)
300 yl4(i) = smstr(ii(6)+i)
301 zl1(i) = zero
302 area(i)=half*((xl2(i)-xl4(i))*yl3(i)-xl3(i)*(yl2(i)-yl4(i)))
303 area_i(i)=one/max(em20,area(i))
304 ELSE
305 smstr(ii(1)+i)=xl2(i)
306 smstr(ii(2)+i)=yl2(i)
307 smstr(ii(3)+i)=xl3(i)
308 smstr(ii(4)+i)=yl3(i)
309 smstr(ii(5)+i)=xl4(i)
310 smstr(ii(6)+i)=yl4(i)
311 ENDIF
312 ENDDO
313 ELSEIF(ismstr==1.OR.ismstr==2)THEN
314 DO i=jft,jlt
315 IF(abs(offg(i))==two)THEN
316 xl2(i)=smstr(ii(1)+i)
317 yl2(i)=smstr(ii(2)+i)
318 xl3(i)=smstr(ii(3)+i)
319 yl3(i)=smstr(ii(4)+i)
320 xl4(i)=smstr(ii(5)+i)
321 yl4(i)=smstr(ii(6)+i)
322 zl1(i)=zero
323 area(i)=half*((xl2(i)-xl4(i))*yl3(i)-xl3(i)*(yl2(i)-yl4(i)))
324 area_i(i)=one/max(em20,area(i))
325 ELSE
326 smstr(ii(1)+i)=xl2(i)
327 smstr(ii(2)+i)=yl2(i)
328 smstr(ii(3)+i)=xl3(i)
329 smstr(ii(4)+i)=yl3(i)
330 smstr(ii(5)+i)=xl4(i)
331 smstr(ii(6)+i)=yl4(i)
332 ENDIF
333 ENDDO
334 ENDIF
335 IF(ismstr==1)THEN
336 DO i=jft,jlt
337 IF(offg(i)==one)offg(i)=two
338 ENDDO
339 ENDIF
340C----------------------------
341C ORTHOTROPY
342C----------------------------
343 CALL cortdir3(elbuf_str,dir_a ,dir_b ,jft ,jlt ,
344 . nlay ,irep ,rx ,ry ,rz ,
345 . sx ,sy ,ssz ,r11 ,r21 ,
346 . r31 ,r12 ,r22 ,r32 ,nel )
347C----------------------------
348 DO ep=jft,jlt
349 lxyz0(1)=fourth*(xl2(ep)+xl3(ep)+xl4(ep))
350 lxyz0(2)=fourth*(yl2(ep)+yl3(ep)+yl4(ep))
351 vcore(ep,1,1)=-lxyz0(1)
352 vcore(ep,1,2)=xl2(ep)-lxyz0(1)
353 vcore(ep,1,3)=xl3(ep)-lxyz0(1)
354 vcore(ep,1,4)=xl4(ep)-lxyz0(1)
355 vcore(ep,2,1)=-lxyz0(2)
356 vcore(ep,2,2)=yl2(ep)-lxyz0(2)
357 vcore(ep,2,3)=yl3(ep)-lxyz0(2)
358 vcore(ep,2,4)=yl4(ep)-lxyz0(2)
359C
360 x13_t(ep) =(vcore(ep,1,1)-vcore(ep,1,3))*half
361 x24_t(ep) =(vcore(ep,1,2)-vcore(ep,1,4))*half
362 y13_t(ep) =(vcore(ep,2,1)-vcore(ep,2,3))*half
363 y24_t(ep) =(vcore(ep,2,2)-vcore(ep,2,4))*half
364 l13(ep)=x13_t(ep)*x13_t(ep)+y13_t(ep)*y13_t(ep)
365 l24(ep)=x24_t(ep)*x24_t(ep)+y24_t(ep)*y24_t(ep)
366 ll(ep)=max(l13(ep),l24(ep))
367 c1 =vcore(ep,1,2)*vcore(ep,2,4)-vcore(ep,2,2)*vcore(ep,1,4)
368 c2 =vcore(ep,1,1)*vcore(ep,2,3)-vcore(ep,2,1)*vcore(ep,1,3)
369 lm(ep) =max(abs(c1),abs(c2))
370 ENDDO
371c
372 DO ep=jft,jlt
373
374 nn(1)=ixc(2,ep)
375 nn(2)=ixc(3,ep)
376 nn(3)=ixc(4,ep)
377 nn(4)=ixc(5,ep)
378
379 vl1(ep,1)=v(1,nn(1))
380 vl1(ep,2)=v(2,nn(1))
381 vl1(ep,3)=v(3,nn(1))
382 vl2(ep,1)=v(1,nn(2))
383 vl2(ep,2)=v(2,nn(2))
384 vl2(ep,3)=v(3,nn(2))
385 vl3(ep,1)=v(1,nn(3))
386 vl3(ep,2)=v(2,nn(3))
387 vl3(ep,3)=v(3,nn(3))
388 vl4(ep,1)=v(1,nn(4))
389 vl4(ep,2)=v(2,nn(4))
390 vl4(ep,3)=v(3,nn(4))
391
392 vg13(1)=vl1(ep,1)-vl3(ep,1)
393 vg24(1)=vl2(ep,1)-vl4(ep,1)
394 vghi(1)=vl1(ep,1)-vl2(ep,1)+vl3(ep,1)-vl4(ep,1)
395C
396 vg13(2)=vl1(ep,2)-vl3(ep,2)
397 vg24(2)=vl2(ep,2)-vl4(ep,2)
398 vghi(2)=vl1(ep,2)-vl2(ep,2)+vl3(ep,2)-vl4(ep,2)
399C
400 vg13(3)=vl1(ep,3)-vl3(ep,3)
401 vg24(3)=vl2(ep,3)-vl4(ep,3)
402 vghi(3)=vl1(ep,3)-vl2(ep,3)+vl3(ep,3)-vl4(ep,3)
403C
404 v13(ep,1)=(vq(ep,1,1)*vg13(1)+vq(ep,2,1)*vg13(2)+vq(ep,3,1)*vg13(3))
405 v24(ep,1)=(vq(ep,1,1)*vg24(1)+vq(ep,2,1)*vg24(2)+vq(ep,3,1)*vg24(3))
406 vhi(ep,1)=(vq(ep,1,1)*vghi(1)+vq(ep,2,1)*vghi(2)+vq(ep,3,1)*vghi(3))
407 v13(ep,2)=(vq(ep,1,2)*vg13(1)+vq(ep,2,2)*vg13(2)+vq(ep,3,2)*vg13(3))
408 v24(ep,2)=(vq(ep,1,2)*vg24(1)+vq(ep,2,2)*vg24(2)+vq(ep,3,2)*vg24(3))
409 vhi(ep,2)=(vq(ep,1,2)*vghi(1)+vq(ep,2,2)*vghi(2)+vq(ep,3,2)*vghi(3))
410 v13(ep,3)=(vq(ep,1,3)*vg13(1)+vq(ep,2,3)*vg13(2)+vq(ep,3,3)*vg13(3))
411 v24(ep,3)=(vq(ep,1,3)*vg24(1)+vq(ep,2,3)*vg24(2)+vq(ep,3,3)*vg24(3))
412 vhi(ep,3)=(vq(ep,1,3)*vghi(1)+vq(ep,2,3)*vghi(2)+vq(ep,3,3)*vghi(3))
413 ENDDO
414C
415 IF (impl_s==0.OR.imp_lr>0) THEN
416 dt05 = half*dt1
417 dt025 =fourth*dt1
418 DO i=jft,jlt
419 exz = y24_t(i)*v13(i,3)-y13_t(i)*v24(i,3)
420 eyz = -x24_t(i)*v13(i,3)+x13_t(i)*v24(i,3)
421 ddry=dt05*exz*area_i(i)
422 ddrx=dt05*eyz*area_i(i)
423 v13x = v13(i,1)
424 v24x = v24(i,1)
425 vhix = vhi(i,1)
426 ddrz1=dt025*(v13(i,2)-v24(i,2))/(x13_t(i)-x24_t(i))
427 IF (abs(x13_t(i)-x24_t(i))<em10) ddrz1 = zero
428 v13(i,1) = v13(i,1)-ddry*v13(i,3)-ddrz1*v13(i,2)
429 v24(i,1) = v24(i,1)-ddry*v24(i,3)-ddrz1*v24(i,2)
430 vhi(i,1) = vhi(i,1)-ddry*vhi(i,3)-ddrz1*vhi(i,2)
431 ddrz2=dt025*(v13x+v24x)/(y13_t(i)+y24_t(i))
432 IF (abs(y13_t(i)+y24_t(i))<em10) ddrz2 = zero
433 v13(i,2) = v13(i,2)-ddrx*v13(i,3)-ddrz2*v13x
434 v24(i,2) = v24(i,2)-ddrx*v24(i,3)-ddrz2*v24x
435 vhi(i,2) = vhi(i,2)-ddrx*vhi(i,3)-ddrz2*vhix
436 ENDDO
437 ENDIF
438 nplat=jft-1
439 splat= 0
440C [PM] in order to make all shells as non-flat shells for pinching
441C [PM] flat shells algorithm not implemented for pinching
442 DO ep=jft,jlt
443 z2=zl1(ep)*zl1(ep)
444 IF ((z2<tol*ll(ep).OR.npt==1).AND.npinch==0) THEN
445 nplat=nplat+1
446 iplat(nplat)=ep
447 ELSE
448 splat=splat+1
449 jplat(splat)=ep
450 ENDIF
451 ENDDO
452 DO ep=nplat+1,jlt
453 iplat(ep)=jplat(ep-nplat)
454 ENDDO
455#include "vectorize.inc"
456 DO i=jft,nplat
457 ep =iplat(i)
458 x13 =x13_t(ep)
459 x24 =x24_t(ep)
460 y13 =y13_t(ep)
461 y24 =y24_t(ep)
462C
463 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
464 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
465 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
466 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
467 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
468 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
469 x13_t(ep) =x13*area_i(ep)
470 x24_t(ep) =x24*area_i(ep)
471 y13_t(ep) =y13*area_i(ep)
472 y24_t(ep) =y24*area_i(ep)
473C--------GAMA(2)
474 gama1=-mx13*y24+my13*x24
475 gama2= mx13*y13-my13*x13
476 z2=zl1(ep)*zl1(ep)
477 nn(1)=ixc(2,ep)
478 nn(2)=ixc(3,ep)
479 nn(3)=ixc(4,ep)
480 nn(4)=ixc(5,ep)
481C-------TRANSMET 3DDL ROTATIONS ----
482 rr(1,1) =vq(ep,1,1)*vr(1,nn(1))+vq(ep,2,1)*vr(2,nn(1))+vq(ep,3,1)*vr(3,nn(1))
483 rr(1,2) =vq(ep,1,1)*vr(1,nn(2))+vq(ep,2,1)*vr(2,nn(2))+vq(ep,3,1)*vr(3,nn(2))
484 rr(1,3) =vq(ep,1,1)*vr(1,nn(3))+vq(ep,2,1)*vr(2,nn(3))+vq(ep,3,1)*vr(3,nn(3))
485 rr(1,4) =vq(ep,1,1)*vr(1,nn(4))+vq(ep,2,1)*vr(2,nn(4))+vq(ep,3,1)*vr(3,nn(4))
486 rr(2,1) =vq(ep,1,2)*vr(1,nn(1))+vq(ep,2,2)*vr(2,nn(1))+vq(ep,3,2)*vr(3,nn(1))
487 rr(2,2) =vq(ep,1,2)*vr(1,nn(2))+vq(ep,2,2)*vr(2,nn(2))+vq(ep,3,2)*vr(3,nn(2))
488 rr(2,3) =vq(ep,1,2)*vr(1,nn(3))+vq(ep,2,2)*vr(2,nn(3))+vq(ep,3,2)*vr(3,nn(3))
489 rr(2,4) =vq(ep,1,2)*vr(1,nn(4))+vq(ep,2,2)*vr(2,nn(4))+vq(ep,3,2)*vr(3,nn(4))
490C-------V13(I)->VXYZ(I,1),V24(I)->VXYZ(I,2),VHI(I)->VXYZ(I,3)------------
491 vxyz(ep,1,1)=v13(ep,1)
492 vxyz(ep,2,1)=v13(ep,2)
493 vxyz(ep,3,1)=v13(ep,3)
494C
495 vxyz(ep,1,2)=v24(ep,1)
496 vxyz(ep,2,2)=v24(ep,2)
497 vxyz(ep,3,2)=v24(ep,3)
498C
499 vxyz(ep,1,3)=vhi(ep,1)
500 vxyz(ep,2,3)=vhi(ep,2)
501 vxyz(ep,3,3)=vhi(ep,3)
502C-------R13(I)->RXYZ(I,1),R24(I)->RXYZ(I,2),RHI(I)->RXYZ(I,3),RTI(I)->RXYZ(I,4)------------
503 rxyz(ep,1,1)=rr(1,1)-rr(1,3)
504 rxyz(ep,2,1)=rr(2,1)-rr(2,3)
505 rxyz(ep,1,2)=rr(1,2)-rr(1,4)
506 rxyz(ep,2,2)=rr(2,2)-rr(2,4)
507 rxyz(ep,1,3)=rr(1,1)-rr(1,2)+rr(1,3)-rr(1,4)
508 rxyz(ep,2,3)=rr(2,1)-rr(2,2)+rr(2,3)-rr(2,4)
509 rxyz(ep,1,4)=rr(1,1)+rr(1,2)+rr(1,3)+rr(1,4)
510 rxyz(ep,2,4)=rr(2,1)+rr(2,2)+rr(2,3)+rr(2,4)
511C-------BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
512C-------SONT STOCKE DANS VCORE---
513 vcore(ep,1,1)=y24_t(ep)
514 vcore(ep,2,1)=-y13_t(ep)
515 vcore(ep,3,1)=-x24_t(ep)
516 vcore(ep,1,2)= x13_t(ep)
517 vcore(ep,2,2)=gama1*area_i(ep)
518 vcore(ep,3,2)=gama2*area_i(ep)
519 vcore(ep,1,3)=mx23
520 vcore(ep,2,3)=my23
521 vcore(ep,3,3)=mx34
522 vcore(ep,1,4)=my34
523 vcore(ep,2,4)=mx13
524 vcore(ep,3,4)=my13
525C
526 j1=(mx23*my13-mx13*my23)*pg
527 j2=-(mx13*my34-mx34*my13)*pg
528 j0=area(ep)*fourth
529C---- JAC(jacobiens aux points d'integration)peut etre negative------
530 jac(ep,1)=abs(j0+j2-j1)
531 jac(ep,2)=abs(j0+j2+j1)
532 jac(ep,3)=abs(j0-j2+j1)
533 jac(ep,4)=abs(j0-j2-j1)
534 j1=(my23-my34)*pg
535 j2=-(my23+my34)*pg
536 hx(ep,1)=j1/jac(ep,1)
537 hx(ep,2)=j2/jac(ep,2)
538 hx(ep,3)=-j1/jac(ep,3)
539 hx(ep,4)=-j2/jac(ep,4)
540 j1=(mx34-mx23)*pg
541 j2=(mx34+mx23)*pg
542 hy(ep,1)=j1/jac(ep,1)
543 hy(ep,2)=j2/jac(ep,2)
544 hy(ep,3)=-j1/jac(ep,3)
545 hy(ep,4)=-j2/jac(ep,4)
546C
547 ENDDO
548#include "vectorize.inc"
549 DO i=nplat+1,jlt
550 ep =iplat(i)
551 z1=zl1(ep)
552 z2=z1*z1
553 vcore(ep,3,1)=z1
554 vcore(ep,3,2)=-z1
555 vcore(ep,3,3)=z1
556 vcore(ep,3,4)=-z1
557 corel(1,1)=vcore(ep,1,1)
558 corel(1,2)=vcore(ep,1,2)
559 corel(1,3)=vcore(ep,1,3)
560 corel(1,4)=vcore(ep,1,4)
561 corel(2,1)=vcore(ep,2,1)
562 corel(2,2)=vcore(ep,2,2)
563 corel(2,3)=vcore(ep,2,3)
564 corel(2,4)=vcore(ep,2,4)
565 x13 =x13_t(ep)
566 x24 =x24_t(ep)
567 y13 =y13_t(ep)
568 y24 =y24_t(ep)
569cc
570 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
571 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
572 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
573 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
574 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
575 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
576 x13_t(ep) =x13*area_i(ep)
577 x24_t(ep) =x24*area_i(ep)
578 y13_t(ep) =y13*area_i(ep)
579 y24_t(ep) =y24*area_i(ep)
580C--------GAMA(2)
581 gama1=-mx13*y24+my13*x24
582 gama2= mx13*y13-my13*x13
583 nn(1)=ixc(2,ep)
584 nn(2)=ixc(3,ep)
585 nn(3)=ixc(4,ep)
586 nn(4)=ixc(5,ep)
587C--------------------------
588C BUILD [F], [Q], [NM], [A-S], [AS] AT NODES
589C--------------------------
590C----------------------------------------------------
591C CALCULATION OF [F]
592C----------------------------------------------------
593C---------
594C 2A1
595C---------
596 x21 =mx23-mx13
597 x34 =(corel(1,3)-corel(1,4))*half
598 y21 =my23-my13
599 y34 =(corel(2,3)-corel(2,4))*half
600 z21 = -z1
601 z34 = z1
602 l12 = sqrt(x21*x21+y21*y21+z2)
603 l34 = sqrt(x34*x34+y34*y34+z2)
604C---------
605C 2A2
606C---------
607 x41 =mx34-mx13
608 x32 =(corel(1,3)-corel(1,2))*half
609 y41 =my34-my13
610 y32 =(corel(2,3)-corel(2,2))*half
611 z41 = -z1
612 z32 = z1
613C----------
614C CALCULATION OF [QN] N=1,4
615C----------
616 a_4=area(ep)*fourth
617C---------- N =1----------
618 sl=one/max(l12,em20)
619 vqn(ep,1,1)=x21*sl
620 vqn(ep,2,1)=y21*sl
621 vqn(ep,3,1)=z21*sl
622 sz2=a_4-gama1
623 sz=z2*l24(ep)
624 sl=one/sqrt(sz+sz2*sz2)
625 vqn(ep,7,1)=-z1*y24
626 vqn(ep,8,1)= z1*x24
627 vqn(ep,9,1)= sz2*sl
628C
629 vqn(ep,7,3)=-vqn(ep,7,1)
630 vqn(ep,8,3)=-vqn(ep,8,1)
631 vqn(ep,7,1)= vqn(ep,7,1)*sl
632 vqn(ep,8,1)= vqn(ep,8,1)*sl
633C
634 vqn(ep,4,1)= vqn(ep,8,1)*vqn(ep,3,1)-vqn(ep,9,1)*vqn(ep,2,1)
635 vqn(ep,5,1)= vqn(ep,9,1)*vqn(ep,1,1)-vqn(ep,7,1)*vqn(ep,3,1)
636 vqn(ep,6,1)= vqn(ep,7,1)*vqn(ep,2,1)-vqn(ep,8,1)*vqn(ep,1,1)
637C---------- N =3----------
638 sl=one/max(l34,em20)
639 vqn(ep,1,3)=x34*sl
640 vqn(ep,2,3)=y34*sl
641 vqn(ep,3,3)=z34*sl
642 sz2=a_4+gama1
643 sl=one/sqrt(sz+sz2*sz2)
644 vqn(ep,7,3)= vqn(ep,7,3)*sl
645 vqn(ep,8,3)= vqn(ep,8,3)*sl
646 vqn(ep,9,3)= sz2*sl
647C
648 vqn(ep,4,3)= vqn(ep,8,3)*vqn(ep,3,3)-vqn(ep,9,3)*vqn(ep,2,3)
649 vqn(ep,5,3)= vqn(ep,9,3)*vqn(ep,1,3)-vqn(ep,7,3)*vqn(ep,3,3)
650 vqn(ep,6,3)= vqn(ep,7,3)*vqn(ep,2,3)-vqn(ep,8,3)*vqn(ep,1,3)
651C---------- N =2----------
652 vqn(ep,1,2)=vqn(ep,1,1)
653 vqn(ep,2,2)=vqn(ep,2,1)
654 vqn(ep,3,2)=vqn(ep,3,1)
655 sz2=a_4+gama2
656 sz=z2*l13(ep)
657 sl=one/sqrt(sz+sz2*sz2)
658 vqn(ep,7,2)=-z1*y13
659 vqn(ep,8,2)= z1*x13
660 vqn(ep,9,2)= sz2*sl
661 vqn(ep,7,4)=-vqn(ep,7,2)
662 vqn(ep,8,4)=-vqn(ep,8,2)
663 vqn(ep,7,2)= vqn(ep,7,2)*sl
664 vqn(ep,8,2)= vqn(ep,8,2)*sl
665C
666 vqn(ep,4,2)= vqn(ep,8,2)*vqn(ep,3,2)-vqn(ep,9,2)*vqn(ep,2,2)
667 vqn(ep,5,2)= vqn(ep,9,2)*vqn(ep,1,2)-vqn(ep,7,2)*vqn(ep,3,2)
668 vqn(ep,6,2)= vqn(ep,7,2)*vqn(ep,2,2)-vqn(ep,8,2)*vqn(ep,1,2)
669C---------- N =4----------
670 vqn(ep,1,4)=vqn(ep,1,3)
671 vqn(ep,2,4)=vqn(ep,2,3)
672 vqn(ep,3,4)=vqn(ep,3,3)
673 sz2=a_4-gama2
674 sl=one/sqrt(sz+sz2*sz2)
675 vqn(ep,7,4)= vqn(ep,7,4)*sl
676 vqn(ep,8,4)= vqn(ep,8,4)*sl
677 vqn(ep,9,4)= sz2*sl
678C
679 vqn(ep,4,4)= vqn(ep,8,4)*vqn(ep,3,4)-vqn(ep,9,4)*vqn(ep,2,4)
680 vqn(ep,5,4)= vqn(ep,9,4)*vqn(ep,1,4)-vqn(ep,7,4)*vqn(ep,3,4)
681 vqn(ep,6,4)= vqn(ep,7,4)*vqn(ep,2,4)-vqn(ep,8,4)*vqn(ep,1,4)
682C
683!! IF(ANIM_PLY < 0 ) THEN
684!! I1 = INOD(NN(1))
685!! I2 = INOD(NN(2))
686!! I3 = INOD(NN(3))
687!! I4 = INOD(NN(4))
688C
689!! NORM =SQRT(VQN(EP,7,1)**2 + VQN(EP,8,1)**2 + VQN(EP,9,1)**2)
690!! VN_NOD(1,I1)=VN_NOD(1,I1)+VQN(EP,7,1)/MAX(NORM, EM20)
691!! VN_NOD(2,I1)=VN_NOD(2,I1)+VQN(EP,8,1)/MAX(NORM, EM20)
692!! VN_NOD(3,I1)=VN_NOD(3,I1)+VQN(EP,9,1)/MAX(NORM, EM20)
693C
694!! NORM =SQRT(VQN(EP,7,2)**2 + VQN(EP,8,2)**2 + VQN(EP,9,2)**2)
695!! VN_NOD(1,I2)=VN_NOD(1,I2)+VQN(EP,7,2)/MAX(NORM, EM20)
696!! VN_NOD(2,I2)=VN_NOD(2,I2)+VQN(EP,8,2)/MAX(NORM, EM20)
697!! VN_NOD(3,I2)=VN_NOD(3,I2)+VQN(EP,9,2)/MAX(NORM, EM20)
698C
699!! NORM =SQRT(VQN(EP,7,3)**2 + VQN(EP,8,3)**2 + VQN(EP,9,3)**2)
700!! VN_NOD(1,I3)=VN_NOD(1,I3)+VQN(EP,7,3)/MAX(NORM, EM20)
701!! VN_NOD(2,I3)=VN_NOD(2,I3)+VQN(EP,8,3)/MAX(NORM, EM20)
702!! VN_NOD(3,I3)=VN_NOD(3,I3)+VQN(EP,9,3)/MAX(NORM, EM20)
703c
704!! NORM =SQRT(VQN(EP,7,4)**2 + VQN(EP,8,4)**2 + VQN(EP,9,4)**2)
705!! VN_NOD(1,I4)=VN_NOD(1,I4)+VQN(EP,7,4)/MAX(NORM, EM20)
706!! VN_NOD(2,I4)=VN_NOD(2,I4)+VQN(EP,8,4)/MAX(NORM, EM20)
707!! VN_NOD(3,I4)=VN_NOD(3,I4)+VQN(EP,9,4)/MAX(NORM, EM20)
708C
709!! ENDIF
710C--------------------------------------------------
711C COMPUTE AS N (MIDDLE OF EDGE)
712C--------------------------------------------------
713C J=1 COTE
714 vnrm(ep,1,1)=vqn(ep,7,1)+vqn(ep,7,2)
715 vnrm(ep,2,1)=vqn(ep,8,1)+vqn(ep,8,2)
716 vnrm(ep,3,1)=vqn(ep,9,1)+vqn(ep,9,2)
717 c1=sqrt(vnrm(ep,1,1)*vnrm(ep,1,1)+vnrm(ep,2,1)*vnrm(ep,2,1)+vnrm(ep,3,1)*vnrm(ep,3,1))
718 c1 = max(em20,c1)
719 vnrm(ep,1,1)=vnrm(ep,1,1)/c1
720 vnrm(ep,2,1)=vnrm(ep,2,1)/c1
721 vnrm(ep,3,1)=vnrm(ep,3,1)/c1
722 vastn(ep,1,1)=zero
723 vastn(ep,2,1)=l12
724 vastn(ep,3,1)=vastn(ep,1,1)
725 vastn(ep,4,1)=vastn(ep,2,1)
726C J=2
727 vnrm(ep,1,2)=vqn(ep,7,4)+vqn(ep,7,3)
728 vnrm(ep,2,2)=vqn(ep,8,4)+vqn(ep,8,3)
729 vnrm(ep,3,2)=vqn(ep,9,4)+vqn(ep,9,3)
730 c1=sqrt(vnrm(ep,1,2)*vnrm(ep,1,2)+vnrm(ep,2,2)*vnrm(ep,2,2)+vnrm(ep,3,2)*vnrm(ep,3,2))
731 c1 = max(em20,c1)
732 vnrm(ep,1,2)=vnrm(ep,1,2)/c1
733 vnrm(ep,2,2)=vnrm(ep,2,2)/c1
734 vnrm(ep,3,2)=vnrm(ep,3,2)/c1
735 vastn(ep,1,2)=zero
736 vastn(ep,2,2)=l34
737 vastn(ep,3,2)=vastn(ep,1,2)
738 vastn(ep,4,2)=vastn(ep,2,2)
739C J=3
740 vnrm(ep,1,3)=vqn(ep,7,1)+vqn(ep,7,4)
741 vnrm(ep,2,3)=vqn(ep,8,1)+vqn(ep,8,4)
742 vnrm(ep,3,3)=vqn(ep,9,1)+vqn(ep,9,4)
743 c1=sqrt(vnrm(ep,1,3)*vnrm(ep,1,3)+vnrm(ep,2,3)*vnrm(ep,2,3)+vnrm(ep,3,3)*vnrm(ep,3,3))
744 c1 = max(em20,c1)
745 vnrm(ep,1,3)=vnrm(ep,1,3)/c1
746 vnrm(ep,2,3)=vnrm(ep,2,3)/c1
747 vnrm(ep,3,3)=vnrm(ep,3,3)/c1
748 vastn(ep,1,3)=-(x41*vqn(ep,4,1)+y41*vqn(ep,5,1)+z41*vqn(ep,6,1))
749 vastn(ep,2,3)= x41*vqn(ep,1,1)+y41*vqn(ep,2,1)+z41*vqn(ep,3,1)
750 vastn(ep,3,3)=-(x41*vqn(ep,4,4)+y41*vqn(ep,5,4)+z41*vqn(ep,6,4))
751 vastn(ep,4,3)= x41*vqn(ep,1,4)+y41*vqn(ep,2,4)+z41*vqn(ep,3,4)
752C J=4
753 vnrm(ep,1,4)=vqn(ep,7,2)+vqn(ep,7,3)
754 vnrm(ep,2,4)=vqn(ep,8,2)+vqn(ep,8,3)
755 vnrm(ep,3,4)=vqn(ep,9,2)+vqn(ep,9,3)
756 c1=sqrt(vnrm(ep,1,4)*vnrm(ep,1,4)+vnrm(ep,2,4)*vnrm(ep,2,4)+vnrm(ep,3,4)*vnrm(ep,3,4))
757 c1 = max(em20,c1)
758 vnrm(ep,1,4)=vnrm(ep,1,4)/c1
759 vnrm(ep,2,4)=vnrm(ep,2,4)/c1
760 vnrm(ep,3,4)=vnrm(ep,3,4)/c1
761 vastn(ep,1,4)=-(x32*vqn(ep,4,2)+y32*vqn(ep,5,2)+z32*vqn(ep,6,2))
762 vastn(ep,2,4)= x32*vqn(ep,1,2)+y32*vqn(ep,2,2)+z32*vqn(ep,3,2)
763 vastn(ep,3,4)=-(x32*vqn(ep,4,3)+y32*vqn(ep,5,3)+z32*vqn(ep,6,3))
764 vastn(ep,4,4)= x32*vqn(ep,1,3)+y32*vqn(ep,2,3)+z32*vqn(ep,3,3)
765C----------
766C CALCULATION OF [QG] NG=1,4
767C----------
768 a_4=a_4/pg
769 jmx13=mx13*pg
770 jmy13=my13*pg
771 jmz13=z1*pg
772 j2myz=jmz13*jmz13
773C---------- NG =1----------
774 sz2=a_4-gama1
775 sz=z2*l24(ep)
776 sl=sqrt(sz+sz2*sz2)
777 jac(ep,1)=sl*pg
778 sl=one/max(sl,em20)
779 vqg(ep,7,1)=-z1*y24
780 vqg(ep,8,1)= z1*x24
781 vqg(ep,9,1)= sz2*sl
782 vqg(ep,7,3)=-vqg(ep,7,1)
783 vqg(ep,8,3)=-vqg(ep,8,1)
784 vqg(ep,7,1)= vqg(ep,7,1)*sl
785 vqg(ep,8,1)= vqg(ep,8,1)*sl
786C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
787 g1x1=mx23-jmx13
788 g1y1=my23-jmy13
789 c1 = sqrt(g1x1*g1x1 + g1y1*g1y1 +j2myz )
790 g2x1=mx34-jmx13
791 g2y1=my34-jmy13
792 c2 = sqrt(g2x1*g2x1 + g2y1*g2y1 +j2myz)
793 vjfi(ep,1,1)=(g2y1*vqg(ep,9,1)+jmz13*vqg(ep,8,1))
794 vjfi(ep,2,1)=(-jmz13*vqg(ep,7,1)-g2x1*vqg(ep,9,1))
795 vjfi(ep,3,1)=(g2x1*vqg(ep,8,1)-g2y1*vqg(ep,7,1))
796 vqg(ep,1,1)= g1x1*c2+vjfi(ep,1,1)*c1
797 vqg(ep,2,1)= g1y1*c2+vjfi(ep,2,1)*c1
798 vqg(ep,3,1)= -jmz13*c2+vjfi(ep,3,1)*c1
799 sl=sl/pg
800 vjfi(ep,1,1)=vjfi(ep,1,1)*sl
801 vjfi(ep,2,1)=vjfi(ep,2,1)*sl
802 vjfi(ep,3,1)=vjfi(ep,3,1)*sl
803C
804 vjfi(ep,4,1)=(-jmz13*vqg(ep,8,1)-g1y1*vqg(ep,9,1))*sl
805 vjfi(ep,5,1)=(g1x1*vqg(ep,9,1)+jmz13*vqg(ep,7,1))*sl
806 vjfi(ep,6,1)=(g1y1*vqg(ep,7,1)-g1x1*vqg(ep,8,1))*sl
807C
808 sl = sqrt(vqg(ep,1,1)*vqg(ep,1,1) + vqg(ep,2,1)*vqg(ep,2,1) + vqg(ep,3,1)*vqg(ep,3,1))
809 IF ( sl/=zero) sl = one / sl
810 vqg(ep,1,1) = vqg(ep,1,1)*sl
811 vqg(ep,2,1) = vqg(ep,2,1)*sl
812 vqg(ep,3,1) = vqg(ep,3,1)*sl
813C
814 vqg(ep,4,1)= vqg(ep,8,1)*vqg(ep,3,1)-vqg(ep,9,1)*vqg(ep,2,1)
815 vqg(ep,5,1)= vqg(ep,9,1)*vqg(ep,1,1)-vqg(ep,7,1)*vqg(ep,3,1)
816 vqg(ep,6,1)= vqg(ep,7,1)*vqg(ep,2,1)-vqg(ep,8,1)*vqg(ep,1,1)
817C---------- NG =3----------
818 sz2=a_4+gama1
819 sl=sqrt(sz+sz2*sz2)
820 jac(ep,3)=sl*pg
821 sl=one/max(sl,em20)
822 vqg(ep,7,3)= vqg(ep,7,3)*sl
823 vqg(ep,8,3)= vqg(ep,8,3)*sl
824 vqg(ep,9,3)= sz2*sl
825C
826 g1x3=mx23+jmx13
827 g1y3=my23+jmy13
828 j1 = sqrt(g1x3*g1x3 + g1y3*g1y3 +j2myz )
829C--------G2X3=G2X2------
830 g2x2=mx34+jmx13
831 g2y2=my34+jmy13
832 j2 = sqrt(g2x2*g2x2 + g2y2*g2y2 +j2myz)
833 vjfi(ep,1,3)=(g2y2*vqg(ep,9,3)-jmz13*vqg(ep,8,3))
834 vjfi(ep,2,3)=(jmz13*vqg(ep,7,3)-g2x2*vqg(ep,9,3))
835 vjfi(ep,3,3)=(g2x2*vqg(ep,8,3)-g2y2*vqg(ep,7,3))
836 vqg(ep,1,3)= g1x3*j2+vjfi(ep,1,3)*j1
837 vqg(ep,2,3)= g1y3*j2+vjfi(ep,2,3)*j1
838 vqg(ep,3,3)= jmz13*j2+vjfi(ep,3,3)*j1
839 sl=sl/pg
840 vjfi(ep,1,3)=vjfi(ep,1,3)*sl
841 vjfi(ep,2,3)=vjfi(ep,2,3)*sl
842 vjfi(ep,3,3)=vjfi(ep,3,3)*sl
843C
844 vjfi(ep,4,3)=(jmz13*vqg(ep,8,3)-g1y3*vqg(ep,9,3))*sl
845 vjfi(ep,5,3)=(g1x3*vqg(ep,9,3)-jmz13*vqg(ep,7,3))*sl
846 vjfi(ep,6,3)=(g1y3*vqg(ep,7,3)-g1x3*vqg(ep,8,3))*sl
847C
848 sl = sqrt(vqg(ep,1,3)*vqg(ep,1,3) + vqg(ep,2,3)*vqg(ep,2,3) + vqg(ep,3,3)*vqg(ep,3,3))
849 IF ( sl/=zero) sl = one / sl
850 vqg(ep,1,3) = vqg(ep,1,3)*sl
851 vqg(ep,2,3) = vqg(ep,2,3)*sl
852 vqg(ep,3,3) = vqg(ep,3,3)*sl
853C
854 vqg(ep,4,3)= vqg(ep,8,3)*vqg(ep,3,3)-vqg(ep,9,3)*vqg(ep,2,3)
855 vqg(ep,5,3)= vqg(ep,9,3)*vqg(ep,1,3)-vqg(ep,7,3)*vqg(ep,3,3)
856 vqg(ep,6,3)= vqg(ep,7,3)*vqg(ep,2,3)-vqg(ep,8,3)*vqg(ep,1,3)
857C---------- NG =2----------
858 sz2=a_4+gama2
859 sz=z2*l13(ep)
860 sl=sqrt(sz+sz2*sz2)
861 jac(ep,2)=sl*pg
862 sl=one/max(sl,em20)
863 vqg(ep,7,2)=-z1*y13
864 vqg(ep,8,2)= z1*x13
865 vqg(ep,9,2)= sz2*sl
866 vqg(ep,7,4)=-vqg(ep,7,2)
867 vqg(ep,8,4)=-vqg(ep,8,2)
868 vqg(ep,7,2)= vqg(ep,7,2)*sl
869 vqg(ep,8,2)= vqg(ep,8,2)*sl
870C
871 vjfi(ep,1,2)=(g2y2*vqg(ep,9,2)-jmz13*vqg(ep,8,2))
872 vjfi(ep,2,2)=(jmz13*vqg(ep,7,2)-g2x2*vqg(ep,9,2))
873 vjfi(ep,3,2)=(g2x2*vqg(ep,8,2)-g2y2*vqg(ep,7,2))
874 vqg(ep,1,2)= g1x1*j2+vjfi(ep,1,2)*c1
875 vqg(ep,2,2)= g1y1*j2+vjfi(ep,2,2)*c1
876 vqg(ep,3,2)=-jmz13*j2+vjfi(ep,3,2)*c1
877 sl=sl/pg
878 vjfi(ep,1,2)=vjfi(ep,1,2)*sl
879 vjfi(ep,2,2)=vjfi(ep,2,2)*sl
880 vjfi(ep,3,2)=vjfi(ep,3,2)*sl
881C
882 vjfi(ep,4,2)=(-jmz13*vqg(ep,8,2)-g1y1*vqg(ep,9,2))*sl
883 vjfi(ep,5,2)=(g1x1*vqg(ep,9,2)+jmz13*vqg(ep,7,2))*sl
884 vjfi(ep,6,2)=(g1y1*vqg(ep,7,2)-g1x1*vqg(ep,8,2))*sl
885C
886 sl = sqrt(vqg(ep,1,2)*vqg(ep,1,2) + vqg(ep,2,2)*vqg(ep,2,2) + vqg(ep,3,2)*vqg(ep,3,2))
887 IF ( sl/=zero) sl = one / sl
888 vqg(ep,1,2) = vqg(ep,1,2)*sl
889 vqg(ep,2,2) = vqg(ep,2,2)*sl
890 vqg(ep,3,2) = vqg(ep,3,2)*sl
891 vqg(ep,4,2)= vqg(ep,8,2)*vqg(ep,3,2)-vqg(ep,9,2)*vqg(ep,2,2)
892 vqg(ep,5,2)= vqg(ep,9,2)*vqg(ep,1,2)-vqg(ep,7,2)*vqg(ep,3,2)
893 vqg(ep,6,2)= vqg(ep,7,2)*vqg(ep,2,2)-vqg(ep,8,2)*vqg(ep,1,2)
894C---------- NG =4----------
895 sz2=a_4-gama2
896 sl=sqrt(sz+sz2*sz2)
897 jac(ep,4)=sl*pg
898 sl=one/max(sl,em20)
899 vqg(ep,7,4)= vqg(ep,7,4)*sl
900 vqg(ep,8,4)= vqg(ep,8,4)*sl
901 vqg(ep,9,4)= sz2*sl
902C
903 vjfi(ep,1,4)=(g2y1*vqg(ep,9,4)+jmz13*vqg(ep,8,4))
904 vjfi(ep,2,4)=(-jmz13*vqg(ep,7,4)-g2x1*vqg(ep,9,4))
905 vjfi(ep,3,4)=(g2x1*vqg(ep,8,4)-g2y1*vqg(ep,7,4))
906 vqg(ep,1,4)= g1x3*c2+vjfi(ep,1,4)*j1
907 vqg(ep,2,4)= g1y3*c2+vjfi(ep,2,4)*j1
908 vqg(ep,3,4)=jmz13*c2+vjfi(ep,3,4)*j1
909 sl=sl/pg
910 vjfi(ep,1,4)=vjfi(ep,1,4)*sl
911 vjfi(ep,2,4)=vjfi(ep,2,4)*sl
912 vjfi(ep,3,4)=vjfi(ep,3,4)*sl
913C
914 vjfi(ep,4,4)=(jmz13*vqg(ep,8,4)-g1y3*vqg(ep,9,4))*sl
915 vjfi(ep,5,4)=(g1x3*vqg(ep,9,4)-jmz13*vqg(ep,7,4))*sl
916 vjfi(ep,6,4)=(g1y3*vqg(ep,7,4)-g1x3*vqg(ep,8,4))*sl
917C
918 sl = sqrt(vqg(ep,1,4)*vqg(ep,1,4) + vqg(ep,2,4)*vqg(ep,2,4) + vqg(ep,3,4)*vqg(ep,3,4))
919 IF ( sl/=zero) sl = one / sl
920 vqg(ep,1,4) = vqg(ep,1,4)*sl
921 vqg(ep,2,4) = vqg(ep,2,4)*sl
922 vqg(ep,3,4) = vqg(ep,3,4)*sl
923 vqg(ep,4,4)= vqg(ep,8,4)*vqg(ep,3,4)-vqg(ep,9,4)*vqg(ep,2,4)
924 vqg(ep,5,4)= vqg(ep,9,4)*vqg(ep,1,4)-vqg(ep,7,4)*vqg(ep,3,4)
925 vqg(ep,6,4)= vqg(ep,7,4)*vqg(ep,2,4)-vqg(ep,8,4)*vqg(ep,1,4)
926 area(ep)=jac(ep,1)+jac(ep,2)+jac(ep,3)+jac(ep,4)
927 ENDDO
928C ---w/ drilling dof----------
929 IF (isrot>0) THEN
930#include "vectorize.inc"
931 DO i=jft,nplat
932 ep =iplat(i)
933 nn(1)=ixc(2,ep)
934 nn(2)=ixc(3,ep)
935 nn(3)=ixc(4,ep)
936 nn(4)=ixc(5,ep)
937C-------APPLY 3DDL ROTATIONS ----
938 rlz(i,1) =(vq(ep,1,3)*vr(1,nn(1))+vq(ep,2,3)*vr(2,nn(1)) +vq(ep,3,3)*vr(3,nn(1)))
939 rlz(i,2) =(vq(ep,1,3)*vr(1,nn(2))+vq(ep,2,3)*vr(2,nn(2)) +vq(ep,3,3)*vr(3,nn(2)))
940 rlz(i,3) =(vq(ep,1,3)*vr(1,nn(3))+vq(ep,2,3)*vr(2,nn(3)) +vq(ep,3,3)*vr(3,nn(3)))
941 rlz(i,4) =(vq(ep,1,3)*vr(1,nn(4))+vq(ep,2,3)*vr(2,nn(4)) +vq(ep,3,3)*vr(3,nn(4)))
942 ENDDO
943 END IF !(ISROT>0) THEN
944C
945 DO i=nplat+1,jlt
946 ep =iplat(i)
947 z1=zl1(ep)
948 z2=z1*z1
949 corel(1,1)=vcore(ep,1,1)
950 corel(1,2)=vcore(ep,1,2)
951 corel(1,3)=vcore(ep,1,3)
952 corel(1,4)=vcore(ep,1,4)
953 corel(2,1)=vcore(ep,2,1)
954 corel(2,2)=vcore(ep,2,2)
955 corel(2,3)=vcore(ep,2,3)
956 corel(2,4)=vcore(ep,2,4)
957 nn(1)=ixc(2,ep)
958 nn(2)=ixc(3,ep)
959 nn(3)=ixc(4,ep)
960 nn(4)=ixc(5,ep)
961 x13 =(vcore(ep,1,1)-vcore(ep,1,3))*half
962 x24 =(vcore(ep,1,2)-vcore(ep,1,4))*half
963 y13 =(vcore(ep,2,1)-vcore(ep,2,3))*half
964 y24 =(vcore(ep,2,2)-vcore(ep,2,4))*half
965 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
966 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
967C--------------------------
968C APPLY 6DDL-->5DDL
969C--------------------------
970C-----RXYZ(2,NNOD)=[Q]^T RGXYZ(3,NNOD)-----
971 rr(1,1) =vq(ep,1,1)*vr(1,nn(1))+vq(ep,2,1)*vr(2,nn(1))+vq(ep,3,1)*vr(3,nn(1))
972 rr(1,2) =vq(ep,1,1)*vr(1,nn(2))+vq(ep,2,1)*vr(2,nn(2))+vq(ep,3,1)*vr(3,nn(2))
973 rr(1,3) =vq(ep,1,1)*vr(1,nn(3))+vq(ep,2,1)*vr(2,nn(3))+vq(ep,3,1)*vr(3,nn(3))
974 rr(1,4) =vq(ep,1,1)*vr(1,nn(4))+vq(ep,2,1)*vr(2,nn(4))+vq(ep,3,1)*vr(3,nn(4))
975 rr(2,1) =vq(ep,1,2)*vr(1,nn(1))+vq(ep,2,2)*vr(2,nn(1))+vq(ep,3,2)*vr(3,nn(1))
976 rr(2,2) =vq(ep,1,2)*vr(1,nn(2))+vq(ep,2,2)*vr(2,nn(2))+vq(ep,3,2)*vr(3,nn(2))
977 rr(2,3) =vq(ep,1,2)*vr(1,nn(3))+vq(ep,2,2)*vr(2,nn(3))+vq(ep,3,2)*vr(3,nn(3))
978 rr(2,4) =vq(ep,1,2)*vr(1,nn(4))+vq(ep,2,2)*vr(2,nn(4))+vq(ep,3,2)*vr(3,nn(4))
979 rr(3,1) =vq(ep,1,3)*vr(1,nn(1))+vq(ep,2,3)*vr(2,nn(1))+vq(ep,3,3)*vr(3,nn(1))
980 rr(3,2) =vq(ep,1,3)*vr(1,nn(2))+vq(ep,2,3)*vr(2,nn(2))+vq(ep,3,3)*vr(3,nn(2))
981 rr(3,3) =vq(ep,1,3)*vr(1,nn(3))+vq(ep,2,3)*vr(2,nn(3))+vq(ep,3,3)*vr(3,nn(3))
982 rr(3,4) =vq(ep,1,3)*vr(1,nn(4))+vq(ep,2,3)*vr(2,nn(4))+vq(ep,3,3)*vr(3,nn(4))
983C-----------------------full projection------------------
984 ar(1)=-z1*vhi(ep,2)+y13*v13(ep,3)+y24*v24(ep,3)+my13*vhi(ep,3)+rr(1,1)+rr(1,2)+rr(1,3)+rr(1,4)
985 ar(2)= z1*vhi(ep,1)-x13*v13(ep,3)-x24*v24(ep,3)-mx13*vhi(ep,3)+rr(2,1)+rr(2,2)+rr(2,3)+rr(2,4)
986 ar(3)= x13*v13(ep,2)+x24*v24(ep,2)+mx13*vhi(ep,2)-y13*v13(ep,1)-y24*v24(ep,1)-my13*vhi(ep,1)+rr(3,1)+rr(3,2)+rr(3,3)+rr(3,4)
987C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
988 xx1 = corel(1,1)*corel(1,1)+corel(1,2)*corel(1,2)+corel(1,3)*corel(1,3)+corel(1,4)*corel(1,4)
989 yy = corel(2,1)*corel(2,1)+corel(2,2)*corel(2,2)+corel(2,3)*corel(2,3)+corel(2,4)*corel(2,4)
990 xy = corel(1,1)*corel(2,1)+corel(1,2)*corel(2,2)+corel(1,3)*corel(2,3)+corel(1,4)*corel(2,4)
991 xz1 = (corel(1,1)-corel(1,2)+corel(1,3)-corel(1,4))*z1
992 yz = (corel(2,1)-corel(2,2)+corel(2,3)-corel(2,4))*z1
993 zz = four*z2
994 d(1)= yy+zz+4
995 d(2)= xx1+zz+4
996 d(3)= xx1+yy+4
997 d(4)= -xy
998 d(5)= -xz1
999 d(6)= -yz
1000C
1001 abc = d(1)*d(2)*d(3)
1002 xxyz2 = d(1)*d(6)*d(6)
1003 yyxz2 = d(2)*d(5)*d(5)
1004 zzxy2 = d(3)*d(4)*d(4)
1005 deta = abc+2*d(4)*d(5)*d(6)-xxyz2-yyxz2-zzxy2
1006 IF (deta<em20) THEN
1007 deta=one
1008 ELSE
1009 deta=one/deta
1010 ENDIF
1011 di(ep,1) = (abc-xxyz2)*deta/d(1)
1012 di(ep,2) = (abc-yyxz2)*deta/d(2)
1013 di(ep,3) = (abc-zzxy2)*deta/d(3)
1014 di(ep,4) = (d(5)*d(6)-d(4)*d(3))*deta
1015 di(ep,5) = (d(6)*d(4)-d(5)*d(2))*deta
1016 di(ep,6) = (d(4)*d(5)-d(6)*d(1))*deta
1017C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1018 alr(1) =di(ep,1)*ar(1)+di(ep,4)*ar(2)+di(ep,5)*ar(3)
1019 alr(2) =di(ep,4)*ar(1)+di(ep,2)*ar(2)+di(ep,6)*ar(3)
1020 alr(3) =di(ep,5)*ar(1)+di(ep,6)*ar(2)+di(ep,3)*ar(3)
1021 v13(ep,1)= half*v13(ep,1)+alr(3)*y13
1022 v24(ep,1)= half*v24(ep,1)+alr(3)*y24
1023 vhi(ep,1)= fourth*vhi(ep,1)+(alr(3)*my13-z1*alr(2))
1024 v13(ep,2)= half*v13(ep,2)-alr(3)*x13
1025 v24(ep,2)= half*v24(ep,2)-alr(3)*x24
1026 vhi(ep,2)= fourth*vhi(ep,2)-(alr(3)*mx13-z1*alr(1))
1027 v13(ep,3)= half*v13(ep,3)-(y13*alr(1)-x13*alr(2))
1028 v24(ep,3)= half*v24(ep,3)-(y24*alr(1)-x24*alr(2))
1029 vhi(ep,3)= fourth*vhi(ep,3)+(mx13*alr(2)-my13*alr(1))
1030
1031 vxyz(ep,1,1)= v13(ep,1)+vhi(ep,1)
1032 vxyz(ep,1,2)= v24(ep,1)-vhi(ep,1)
1033 vxyz(ep,1,3)= -v13(ep,1)+vhi(ep,1)
1034 vxyz(ep,1,4)= -v24(ep,1)-vhi(ep,1)
1035
1036 vxyz(ep,2,1)= v13(ep,2)+vhi(ep,2)
1037 vxyz(ep,2,2)= v24(ep,2)-vhi(ep,2)
1038 vxyz(ep,2,3)= -v13(ep,2)+vhi(ep,2)
1039 vxyz(ep,2,4)= -v24(ep,2)-vhi(ep,2)
1040
1041 vxyz(ep,3,1)= v13(ep,3)+vhi(ep,3)
1042 vxyz(ep,3,2)= v24(ep,3)-vhi(ep,3)
1043 vxyz(ep,3,3)= -v13(ep,3)+vhi(ep,3)
1044 vxyz(ep,3,4)= -v24(ep,3)-vhi(ep,3)
1045C
1046 rr(1,1)= rr(1,1)-alr(1)
1047 rr(1,2)= rr(1,2)-alr(1)
1048 rr(1,3)= rr(1,3)-alr(1)
1049 rr(1,4)= rr(1,4)-alr(1)
1050C
1051 rr(2,1)= rr(2,1)-alr(2)
1052 rr(2,2)= rr(2,2)-alr(2)
1053 rr(2,3)= rr(2,3)-alr(2)
1054 rr(2,4)= rr(2,4)-alr(2)
1055C
1056 rr(3,1)= rr(3,1)-alr(3)
1057 rr(3,2)= rr(3,2)-alr(3)
1058 rr(3,3)= rr(3,3)-alr(3)
1059 rr(3,4)= rr(3,4)-alr(3)
1060 rxyz(ep,1,1)=vqn(ep,1,1)*rr(1,1)+vqn(ep,2,1)*rr(2,1)+vqn(ep,3,1)*rr(3,1)
1061 rxyz(ep,2,1)=vqn(ep,4,1)*rr(1,1)+vqn(ep,5,1)*rr(2,1)+vqn(ep,6,1)*rr(3,1)
1062 rxyz(ep,1,2)=vqn(ep,1,2)*rr(1,2)+vqn(ep,2,2)*rr(2,2)+vqn(ep,3,2)*rr(3,2)
1063 rxyz(ep,2,2)=vqn(ep,4,2)*rr(1,2)+vqn(ep,5,2)*rr(2,2)+vqn(ep,6,2)*rr(3,2)
1064 rxyz(ep,1,3)=vqn(ep,1,3)*rr(1,3)+vqn(ep,2,3)*rr(2,3)+vqn(ep,3,3)*rr(3,3)
1065 rxyz(ep,2,3)=vqn(ep,4,3)*rr(1,3)+vqn(ep,5,3)*rr(2,3)+vqn(ep,6,3)*rr(3,3)
1066 rxyz(ep,1,4)=vqn(ep,1,4)*rr(1,4)+vqn(ep,2,4)*rr(2,4)+vqn(ep,3,4)*rr(3,4)
1067 rxyz(ep,2,4)=vqn(ep,4,4)*rr(1,4)+vqn(ep,5,4)*rr(2,4)+vqn(ep,6,4)*rr(3,4)
1068C--------drilling dof
1069 IF (isrot>0) THEN
1070 rlz(i,1)=(vqn(ep,7,1)*rr(1,1)+vqn(ep,8,1)*rr(2,1)+vqn(ep,9,1)*rr(3,1))
1071 rlz(i,2)=(vqn(ep,7,2)*rr(1,2)+vqn(ep,8,2)*rr(2,2)+vqn(ep,9,2)*rr(3,2))
1072 rlz(i,3)=(vqn(ep,7,3)*rr(1,3)+vqn(ep,8,3)*rr(2,3)+vqn(ep,9,3)*rr(3,3))
1073 rlz(i,4)=(vqn(ep,7,4)*rr(1,4)+vqn(ep,8,4)*rr(2,4)+vqn(ep,9,4)*rr(3,4))
1074 END IF !(ISROT>0) THEN
1075 END DO
1076C
1077C----------------------------
1078C---------Factor for characteristic length
1079C----------------------------
1080 DO i=jft,jlt
1081 rx(i)=xl2(i)+xl3(i)-xl4(i)
1082 ry(i)=yl2(i)+yl3(i)-yl4(i)
1083 sx(i)=-xl2(i)+xl3(i)+xl4(i)
1084 sy(i)=-yl2(i)+yl3(i)+yl4(i)
1085 c1=sqrt(rx(i)*rx(i)+ry(i)*ry(i))
1086 c2=sqrt(sx(i)*sx(i)+sy(i)*sy(i))
1087 s1=fourth*(max(c1,c2)/min(c1,c2)-one)
1088 fac1=min(half,s1)+one
1089 fac2=four*area(i)/(c1*c2)
1090 fac2=3.413*max(zero,fac2-0.7071)
1091 fac2=0.78+0.22*fac2*fac2*fac2
1092 faci=two*fac1*fac2
1093C
1094 facn(i,1)=sqrt(l24(i)/ll(i))
1095 facn(i,2)=sqrt(l13(i)/ll(i))
1096 s1 = sqrt(faci*(four_over_3+lm(i)*area_i(i))*ll(i))
1097 s1 = max(s1,em20)
1098 ll(i) = area(i)/s1
1099 ENDDO
1100
1101 off_l = zero
1102 DO ep=jft,jlt
1103C---------factor for characteristic length---
1104 off(ep) = min(one,abs(offg(ep)))
1105 off_l = min(off_l,offg(ep))
1106 ENDDO
1107C
1108 IF(off_l<zero)THEN
1109 DO ep=jft,jlt
1110 IF(offg(ep)<zero)THEN
1111 vxyz(ep,1,1)= zero
1112 vxyz(ep,2,1)= zero
1113 vxyz(ep,3,1)= zero
1114 vxyz(ep,1,2)= zero
1115 vxyz(ep,2,2)= zero
1116 vxyz(ep,3,2)= zero
1117 vxyz(ep,1,3)= zero
1118 vxyz(ep,2,3)= zero
1119 vxyz(ep,3,3)= zero
1120 vxyz(ep,1,4)= zero
1121 vxyz(ep,2,4)= zero
1122 vxyz(ep,3,4)= zero
1123C
1124 rxyz(ep,1,1)= zero
1125 rxyz(ep,2,1)= zero
1126 rxyz(ep,1,2)= zero
1127 rxyz(ep,2,2)= zero
1128 rxyz(ep,1,3)= zero
1129 rxyz(ep,2,3)= zero
1130 rxyz(ep,1,4)= zero
1131 rxyz(ep,2,4)= zero
1132 ENDIF
1133 ENDDO
1134 ENDIF
1135 IF(anim_ply > 0) THEN
1136#include "lockon.inc"
1137 DO ep=jft,jlt
1138 i1 = inod(ixc(2,ep))
1139 i2 = inod(ixc(3,ep))
1140 i3 = inod(ixc(4,ep))
1141 i4 = inod(ixc(5,ep))
1142cc
1143 vn_nod(1,i1) = vn_nod(1,i1)+vq(ep,1,3)
1144 vn_nod(2,i1) = vn_nod(2,i1)+vq(ep,2,3)
1145 vn_nod(3,i1) = vn_nod(3,i1)+vq(ep,3,3)
1146C
1147 vn_nod(1,i2) =vn_nod(1,i2)+vq(ep,1,3)
1148 vn_nod(2,i2) =vn_nod(2,i2)+vq(ep,2,3)
1149 vn_nod(3,i2) =vn_nod(3,i2)+vq(ep,3,3)
1150C
1151 vn_nod(1,i3) =vn_nod(1,i3)+vq(ep,1,3)
1152 vn_nod(2,i3) =vn_nod(2,i3)+vq(ep,2,3)
1153 vn_nod(3,i3) =vn_nod(3,i3)+vq(ep,3,3)
1154
1155C
1156 vn_nod(1,i4) =vn_nod(1,i4)+vq(ep,1,3)
1157 vn_nod(2,i4) =vn_nod(2,i4)+vq(ep,2,3)
1158 vn_nod(3,i4) =vn_nod(3,i4)+vq(ep,3,3)
1159 ENDDO
1160#include "lockoff.inc"
1161 ENDIF
1162C-------------
1163 RETURN
clskew3
subroutine clskew3(jft, jlt, irep, rx, ry, rz, sx, sy, sz, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, det)
Definition clskew.F:34
cortdir3
subroutine cortdir3(elbuf_str, dir_a, dir_b, jft, jlt, nlay, irep, rx, ry, rz, sx, sy, sz, e1x, e1y, e1z, e2x, e2y, e2z, nel)
Definition cortdir3.F:45
my_real
#define my_real
Definition cppsort.cpp:32
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
plyxfem_mod
Definition plyxfem_mod.F:31

◆ cbacoort()

subroutine cbacoort ( type(elbuf_struct_) elbuf_str,
integer jft,
integer jlt,
x,
v,
vr,
dr,
integer, dimension(nixc,*) ixc,
pm,
offg,
area,
vxyz,
rlz,
vcore,
jac,
hx,
hy,
vq,
vqg,
vjfi,
integer nplat,
integer, dimension(*) iplat,
x13_t,
x24_t,
y13_t,
y24_t,
integer npt,
double precision, dimension(*) smstr,
integer isrot,
xlcor,
zl,
vqn,
integer nel )

Definition at line 1174 of file cbacoor.F.

1180C-----------------------------------------------
1181C M o d u l e s
1182C-----------------------------------------------
1183 USE elbufdef_mod
1184C-----------------------------------------------
1185C I M P L I C I T T Y P E S
1186C-----------------------------------------------
1187#include "implicit_f.inc"
1188#include "comlock.inc"
1189c-----------------------------------------------
1190c g l o b a l p a r a m e t e r s
1191c-----------------------------------------------
1192#include "mvsiz_p.inc"
1193#include "param_c.inc"
1194C-----------------------------------------------
1195C D U M M Y A R G U M E N T S
1196C-----------------------------------------------
1197 INTEGER JFT,JLT,NNOD,NPLAT,NPT,ISROT,NEL
1198 INTEGER IXC(NIXC,*),IPLAT(*)
1199 parameter(nnod = 4)
1200 my_real
1201 . x(3,*), v(3,*), vr(3,*),pm(npropm,*),offg(*),
1202 . vcore(mvsiz,3,nnod),vxyz(mvsiz,3,nnod),area(*),vjfi(mvsiz,6,4),
1203 . vq(mvsiz,3,3),vqg(mvsiz,9,nnod),rlz(mvsiz,nnod),
1204 . jac(mvsiz,4),hx(mvsiz,4),hy(mvsiz,4),y24_t(*),dr(3,*),
1205 . x13_t(*),x24_t(*),y13_t(*),xlcor(mvsiz,2,nnod-1),zl(*),vqn(mvsiz,9,nnod)
1206 double precision
1207 . smstr(*)
1208 TYPE(ELBUF_STRUCT_) :: ELBUF_STR
1209C-----RLZ(NNOD,*)
1210C-----------------------------------------------
1211c FUNCTION: premilitary compute for QBAT and Ismstr=10
1212c
1213c Note:
1214c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
1215c
1216c TYPE NAME FUNCTION
1217c I JFT,JLT - element id limit
1218c I X, V, VR(3,NUMNOD) - coordinate,velocity, rotational velocity in global system
1219c I IXC(NIXC,*NEL) - element connectivity of shell (and other data)
1220c I PM(NPROPM,MID) - input Material data
1221c I OFFG(NEL),OFF(NEL) - element activation flag, local flag
1222c O AREA(NEL),NTP - Area, num. of integrating points in thickness
1223c O VCORE(MVSIZ,3,4) - coordinates of 4 nodes in local system
1224c BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
1225c for plat element
1226c O VXYZ(3,4,NEL) - nodal velocity (4 nodes) in local system
1227C V13->VXYZ(,1,),V24->VXYZ(,2,),VHI->VXYZ(,3,)------------
1228c for plat element
1229c O RXYZ(3,4,NEL) - nodal rotational velocity (4 nodes) in local system 5 dof
1230C R13->RXYZ(,1,),R24->RXYZ(,2,),RHI->RXYZ(,3,),RTI->RXYZ(,4,)----
1231c for plat element
1232c O JAC,HX,HY(4,NEL) - Jacobien,asymmetric part(hourglass) at 4 Gauss points
1233c O VQ(9,NEL),VQN,VQG(9,4,) - local system of element, nodal(4) and Gauss points(4)
1234c O Xij_T,Yij_T(NEL) - [B] components used for in-plane shear assumed strain
1235c IO XiS,YiS,ZiS(NEL) - reference configurations used for small strain options
1236c I ISMSTR - small strain
1237c O ZL1(NEL) - warped length
1238c I IDRIL - drilling dof flag
1239c (used only Idril active )
1240C-----------------------------------------------
1241C L O C A L V A R I A B L E S
1242C-----------------------------------------------
1243 INTEGER I,II(9),J,K,L,M,I1,I2,I3,I4,EP,SPLAT,JJ
1244 INTEGER NN(4),JPLAT(MVSIZ),IFPROJ
1245 my_real
1246 . xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4
1247 my_real
1248 . c1,c2,l13(mvsiz),l24(mvsiz),
1249 . aa,vv(3,nnod),rr(3,nnod),off_l
1250 my_real
1251 . tol,pg,j0,j1,j2,deta,corel(3,4),x1,y1,s1,ll(mvsiz)
1252 my_real
1253 . x13,x24,y13,mx13,mx23,mx34,my13,y13_2,z1,z2,gama1,gama2,
1254 . x21,x34,y21,y34,z21,z34,x41,x32,y41,y32,z_2,z41,z32,l12,l34,
1255 . a_4,sl,sz2,sz,jmx13,jmy13,jmz13,j2myz,lm(mvsiz),y24,y24_2,
1256 . my23,my34,scal,g1x1,g1x3,g1y1,g1y3,g2x1,g2x2,g2y1,g2y2,
1257 . ar(3),ad(nnod),btb(6),xx1,yy,zz,xy,xz1,yz,d(6),
1258 . alr(3),ald(nnod),dbad(3),abc,xxyz2,zzxy2,yyxz2
1259 my_real
1260 . lxyz0(3),deta1(mvsiz),
1261 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),
1262 . yl3(mvsiz),yl4(mvsiz),xx,area_i(mvsiz),
1263 . x0g2(mvsiz),x0g3(mvsiz),x0g4(mvsiz),y0g2(mvsiz),
1264 . y0g3(mvsiz),y0g4(mvsiz),z0g2(mvsiz),z0g3(mvsiz),z0g4(mvsiz)
1265C
1266 my_real
1267 . vg13(3),vg24(3),vghi(3),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
1268 . faci,fac1,fac2,vr1_12,vr1_21,
1269 . dt05,dt025,exz,eyz,ddrx,ddry,v13x,v24x,vhix,ddrz1,ddrz2,ddrz,
1270 . vnx1, vny1, vnz1,vnx2,vny2,vnz2,vnx3,vny3,vnz3,vnx4,vny4,vnz4,
1271 . norm,di(mvsiz,6),zl1(mvsiz),ssz(mvsiz),
1272 . rx(mvsiz),ry(mvsiz),rz(mvsiz),sx(mvsiz),sy(mvsiz),
1273 . r11(mvsiz),r12(mvsiz),r13(mvsiz),r21(mvsiz),r22(mvsiz),
1274 . r23(mvsiz),r31(mvsiz),r32(mvsiz),r33(mvsiz),dirz(mvsiz,2),
1275 . x0l2(mvsiz),x0l3(mvsiz),x0l4(mvsiz),y0l2(mvsiz),
1276 . y0l3(mvsiz),y0l4(mvsiz),nx,ny,nz,det,vq0(mvsiz,3,3),det_1
1277C-------- [F] will contains only Rz ->2D
1278 my_real
1279 . axyz(mvsiz,3,nnod)
1280 DATA pg/.577350269189626/
1281C=======================================================================
1282 DO i=1,9
1283 ii(i) = nel*(i-1)
1284 ENDDO
1285C
1286 tol=em12
1287 IF (isrot > 0 ) tol=em8
1288 DO i=jft,jlt
1289 IF(abs(offg(i))==one)offg(i)=sign(two,offg(i))
1290 axyz(i,1:3,1:4)= zero
1291C
1292 IF (isrot > 0 ) THEN
1293 nn(1)=ixc(2,i)
1294 nn(2)=ixc(3,i)
1295 nn(3)=ixc(4,i)
1296 nn(4)=ixc(5,i)
1297 axyz(i,1,1) = dr(1,nn(1))
1298 axyz(i,2,1) = dr(2,nn(1))
1299 axyz(i,3,1) = dr(3,nn(1))
1300 axyz(i,1,2) = dr(1,nn(2))
1301 axyz(i,2,2) = dr(2,nn(2))
1302 axyz(i,3,2) = dr(3,nn(2))
1303 axyz(i,1,3) = dr(1,nn(3))
1304 axyz(i,2,3) = dr(2,nn(3))
1305 axyz(i,3,3) = dr(3,nn(3))
1306 axyz(i,1,4) = dr(1,nn(4))
1307 axyz(i,2,4) = dr(2,nn(4))
1308 axyz(i,3,4) = dr(3,nn(4))
1309 END IF !(ISROT > 0 ) THEN
1310
1311 x0g2(i) = smstr(ii(1)+i)
1312 y0g2(i) = smstr(ii(2)+i)
1313 z0g2(i) = smstr(ii(3)+i)
1314 x0g3(i) = smstr(ii(4)+i)
1315 y0g3(i) = smstr(ii(5)+i)
1316 z0g3(i) = smstr(ii(6)+i)
1317 x0g4(i) = smstr(ii(7)+i)
1318 y0g4(i) = smstr(ii(8)+i)
1319 z0g4(i) = smstr(ii(9)+i)
1320 ENDDO
1321C--- in [F](2,2) only Rz is inside, Rx,Ry have been excluded
1322C-bases at T=0 : B0_g : x0g(3,nnod-1) saved in SMSTR, w/ origine at n1
1323C--- B0_l : x0l(2,nnod-1)+Z1,VCORE at only T=0, but can be calculated from x0g(3,nnod-1)
1324C-bases at T>0 : B_G : X(3,numnod), B_l : by rotate B0_l w/ Rz0
1325C---------UL : xl - x0l' (rotated around Rz0)
1326C-- normal in initial conf. to compute Z1
1327 DO i=jft,jlt
1328 rx(i)=x0g2(i)+x0g3(i)-x0g4(i)
1329 sx(i)=x0g3(i)+x0g4(i)-x0g2(i)
1330 ry(i)=y0g2(i)+y0g3(i)-y0g4(i)
1331 sy(i)=y0g3(i)+y0g4(i)-y0g2(i)
1332 rz(i)=z0g2(i)+z0g3(i)-z0g4(i)
1333 ssz(i)=z0g3(i)+z0g4(i)-z0g2(i)
1334 ENDDO
1335C----------------------------
1336C LOCAL SYSTEM VQ0
1337C----------------------------
1338 k = 0
1339 CALL clskew3(jft,jlt,k,
1340 . rx, ry, rz,
1341 . sx, sy, ssz,
1342 . r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
1343 DO i=jft,jlt
1344 area(i)=fourth*deta1(i)
1345 area_i(i)=one/area(i)
1346 vq0(i,1,1)=r11(i)
1347 vq0(i,2,1)=r21(i)
1348 vq0(i,3,1)=r31(i)
1349 vq0(i,1,2)=r12(i)
1350 vq0(i,2,2)=r22(i)
1351 vq0(i,3,2)=r32(i)
1352 vq0(i,1,3)=r13(i)
1353 vq0(i,2,3)=r23(i)
1354 vq0(i,3,3)=r33(i)
1355 ENDDO
1356 DO i=jft,jlt
1357 lxyz0(1)=fourth*(x0g2(i)+x0g3(i)+x0g4(i))
1358 lxyz0(2)=fourth*(y0g2(i)+y0g3(i)+y0g4(i))
1359 lxyz0(3)=fourth*(z0g2(i)+z0g3(i)+z0g4(i))
1360 zl1(i)=-(vq0(i,1,3)*lxyz0(1)+vq0(i,2,3)*lxyz0(2)+vq0(i,3,3)*lxyz0(3))
1361 ENDDO
1362C----------------------------
1363C xl =R1^t x0l; R1^t=(VQ0^t*VQ)^t---extract Rz0 of R1
1364C----------------------------
1365 DO i=jft,jlt
1366 vr1_12=vq0(i,1,1)*vq(i,1,2)+vq0(i,2,1)*vq(i,2,2)+vq0(i,3,1)*vq(i,3,2)
1367 vr1_21=vq0(i,1,2)*vq(i,1,1)+vq0(i,2,2)*vq(i,2,1)+vq0(i,3,1)*vq(i,3,2)
1368 dirz(i,2) = half*(vr1_12-vr1_21)
1369 norm = one-dirz(i,2)*dirz(i,2)
1370 dirz(i,1) = sqrt(max(zero,norm))
1371 ENDDO
1372 DO i=jft,jlt
1373 x0l2(i)=vq0(i,1,1)*x0g2(i)+vq0(i,2,1)*y0g2(i)+vq0(i,3,1)*z0g2(i)
1374 y0l2(i)=vq0(i,1,2)*x0g2(i)+vq0(i,2,2)*y0g2(i)+vq0(i,3,2)*z0g2(i)
1375 x0l3(i)=vq0(i,1,1)*x0g3(i)+vq0(i,2,1)*y0g3(i)+vq0(i,3,1)*z0g3(i)
1376 y0l3(i)=vq0(i,1,2)*x0g3(i)+vq0(i,2,2)*y0g3(i)+vq0(i,3,2)*z0g3(i)
1377 x0l4(i)=vq0(i,1,1)*x0g4(i)+vq0(i,2,1)*y0g4(i)+vq0(i,3,1)*z0g4(i)
1378 y0l4(i)=vq0(i,1,2)*x0g4(i)+vq0(i,2,2)*y0g4(i)+vq0(i,3,2)*z0g4(i)
1379 ENDDO
1380C----------------------------
1381C Rotate x0l of Rz0
1382C----------------------------
1383 DO i=jft,jlt
1384 xl2(i)= x0l2(i)*dirz(i,1)-y0l2(i)*dirz(i,2)
1385 yl2(i)= x0l2(i)*dirz(i,2)+y0l2(i)*dirz(i,1)
1386 xl3(i)= x0l3(i)*dirz(i,1)-y0l3(i)*dirz(i,2)
1387 yl3(i)= x0l3(i)*dirz(i,2)+y0l3(i)*dirz(i,1)
1388 xl4(i)= x0l4(i)*dirz(i,1)-y0l4(i)*dirz(i,2)
1389 yl4(i)= x0l4(i)*dirz(i,2)+y0l4(i)*dirz(i,1)
1390 x0l2(i)=xl2(i)
1391 y0l2(i)=yl2(i)
1392 x0l3(i)=xl3(i)
1393 y0l3(i)=yl3(i)
1394 x0l4(i)=xl4(i)
1395 y0l4(i)=yl4(i)
1396 ENDDO
1397 DO i=jft,jlt
1398 xl2(i)=xlcor(i,1,1)
1399 yl2(i)=xlcor(i,2,1)
1400 xl3(i)=xlcor(i,1,2)
1401 yl3(i)=xlcor(i,2,2)
1402 xl4(i)=xlcor(i,1,3)
1403 yl4(i)=xlcor(i,2,3)
1404 ENDDO
1405C-----------UL : xl - x0l' (rotated of rz0)
1406 DO i=jft,jlt
1407 v13(i,1)=-xl3(i)-(-x0l3(i))
1408 v24(i,1)=xl2(i)-xl4(i)-(x0l2(i)-x0l4(i))
1409 vhi(i,1)=-xl2(i)+xl3(i)-xl4(i)-(-x0l2(i)+x0l3(i)-x0l4(i))
1410 v13(i,2)=-yl3(i)-(-y0l3(i))
1411 v24(i,2)=yl2(i)-yl4(i)-(y0l2(i)-y0l4(i))
1412 vhi(i,2)=-yl2(i)+yl3(i)-yl4(i)-(-y0l2(i)+y0l3(i)-y0l4(i))
1413 v13(i,3)=zero
1414 v24(i,3)=zero
1415 vhi(i,3)=four*(zl(i)-zl1(i))
1416 ENDDO
1417C-------set up B_l by rotating B0_l (Rz0)
1418 DO i=jft,jlt
1419 xl2(i)=x0l2(i)
1420 yl2(i)=y0l2(i)
1421 xl3(i)=x0l3(i)
1422 yl3(i)=y0l3(i)
1423 xl4(i)=x0l4(i)
1424 yl4(i)=y0l4(i)
1425 ENDDO
1426C----------------------------
1427 DO ep=jft,jlt
1428 lxyz0(1)=fourth*(xl2(ep)+xl3(ep)+xl4(ep))
1429 lxyz0(2)=fourth*(yl2(ep)+yl3(ep)+yl4(ep))
1430 vcore(ep,1,1)=-lxyz0(1)
1431 vcore(ep,1,2)=xl2(ep)-lxyz0(1)
1432 vcore(ep,1,3)=xl3(ep)-lxyz0(1)
1433 vcore(ep,1,4)=xl4(ep)-lxyz0(1)
1434 vcore(ep,2,1)=-lxyz0(2)
1435 vcore(ep,2,2)=yl2(ep)-lxyz0(2)
1436 vcore(ep,2,3)=yl3(ep)-lxyz0(2)
1437 vcore(ep,2,4)=yl4(ep)-lxyz0(2)
1438C
1439 x13_t(ep) =(vcore(ep,1,1)-vcore(ep,1,3))*half
1440 x24_t(ep) =(vcore(ep,1,2)-vcore(ep,1,4))*half
1441 y13_t(ep) =(vcore(ep,2,1)-vcore(ep,2,3))*half
1442 y24_t(ep) =(vcore(ep,2,2)-vcore(ep,2,4))*half
1443 l13(ep)=x13_t(ep)*x13_t(ep)+y13_t(ep)*y13_t(ep)
1444 l24(ep)=x24_t(ep)*x24_t(ep)+y24_t(ep)*y24_t(ep)
1445 ll(ep)=max(l13(ep),l24(ep))
1446 ENDDO
1447c
1448 nplat=jft-1
1449 splat= 0
1450 DO ep=jft,jlt
1451 z2=zl1(ep)*zl1(ep)
1452 IF (z2<tol*ll(ep).OR.npt==1) THEN
1453 nplat=nplat+1
1454 iplat(nplat)=ep
1455 ELSE
1456 splat=splat+1
1457 jplat(splat)=ep
1458 ENDIF
1459 ENDDO
1460 DO ep=nplat+1,jlt
1461 iplat(ep)=jplat(ep-nplat)
1462 ENDDO
1463#include "vectorize.inc"
1464 DO i=jft,nplat
1465 ep =iplat(i)
1466 x13 =x13_t(ep)
1467 x24 =x24_t(ep)
1468 y13 =y13_t(ep)
1469 y24 =y24_t(ep)
1470C
1471 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
1472 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
1473 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
1474 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
1475 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
1476 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
1477 x13_t(ep) =x13*area_i(ep)
1478 x24_t(ep) =x24*area_i(ep)
1479 y13_t(ep) =y13*area_i(ep)
1480 y24_t(ep) =y24*area_i(ep)
1481C--------GAMA(2)
1482 gama1=-mx13*y24+my13*x24
1483 gama2= mx13*y13-my13*x13
1484 z2=zl1(ep)*zl1(ep)
1485CC-------V13(I)->VXYZ(I,1),V24(I)->VXYZ(I,2),VHI(I)->VXYZ(I,3)------------
1486 vxyz(ep,1,1)=v13(ep,1)
1487 vxyz(ep,2,1)=v13(ep,2)
1488 vxyz(ep,3,1)=v13(ep,3)
1489C
1490 vxyz(ep,1,2)=v24(ep,1)
1491 vxyz(ep,2,2)=v24(ep,2)
1492 vxyz(ep,3,2)=v24(ep,3)
1493C
1494 vxyz(ep,1,3)=vhi(ep,1)
1495 vxyz(ep,2,3)=vhi(ep,2)
1496 vxyz(ep,3,3)=vhi(ep,3)
1497C-------BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
1498C-------SONT STOCKE DANS VCORE---
1499 vcore(ep,1,1)=y24_t(ep)
1500 vcore(ep,2,1)=-y13_t(ep)
1501 vcore(ep,3,1)=-x24_t(ep)
1502 vcore(ep,1,2)= x13_t(ep)
1503 vcore(ep,2,2)=gama1*area_i(ep)
1504 vcore(ep,3,2)=gama2*area_i(ep)
1505 vcore(ep,1,3)=mx23
1506 vcore(ep,2,3)=my23
1507 vcore(ep,3,3)=mx34
1508 vcore(ep,1,4)=my34
1509 vcore(ep,2,4)=mx13
1510 vcore(ep,3,4)=my13
1511C
1512 j1=(mx23*my13-mx13*my23)*pg
1513 j2=-(mx13*my34-mx34*my13)*pg
1514 j0=area(ep)*fourth
1515C----
1516 jac(ep,1)=abs(j0+j2-j1)
1517 jac(ep,2)=abs(j0+j2+j1)
1518 jac(ep,3)=abs(j0-j2+j1)
1519 jac(ep,4)=abs(j0-j2-j1)
1520 j1=(my23-my34)*pg
1521 j2=-(my23+my34)*pg
1522 hx(ep,1)=j1/jac(ep,1)
1523 hx(ep,2)=j2/jac(ep,2)
1524 hx(ep,3)=-j1/jac(ep,3)
1525 hx(ep,4)=-j2/jac(ep,4)
1526 j1=(mx34-mx23)*pg
1527 j2=(mx34+mx23)*pg
1528 hy(ep,1)=j1/jac(ep,1)
1529 hy(ep,2)=j2/jac(ep,2)
1530 hy(ep,3)=-j1/jac(ep,3)
1531 hy(ep,4)=-j2/jac(ep,4)
1532 ENDDO
1533C ---w/ drilling dof------drilling may not work -> check RLZj w/ Rz0
1534 IF (isrot>0) THEN
1535#include "vectorize.inc"
1536 DO i=jft,nplat
1537 ep =iplat(i)
1538C-------TRANSMET 3DDL ROTATIONS ----
1539 rlz(ep,1) =vq(ep,1,3)*axyz(ep,1,1)+vq(ep,2,3)*axyz(ep,2,1)+vq(ep,3,3)*axyz(ep,3,1)
1540 rlz(ep,2) =vq(ep,1,3)*axyz(ep,1,2)+vq(ep,2,3)*axyz(ep,2,2)+vq(ep,3,3)*axyz(ep,3,2)
1541 rlz(ep,3) =vq(ep,1,3)*axyz(ep,1,3)+vq(ep,2,3)*axyz(ep,2,3)+vq(ep,3,3)*axyz(ep,3,3)
1542 rlz(ep,4) =vq(ep,1,3)*axyz(ep,1,4)+vq(ep,2,3)*axyz(ep,2,4)+vq(ep,3,3)*axyz(ep,3,4)
1543 ENDDO
1544 END IF !(ISROT>0) THEN
1545C----------group warped elements , VJFI,VQN,VQG
1546 ifproj = 0
1547#include "vectorize.inc"
1548 DO i=nplat+1,jlt
1549 ep =iplat(i)
1550 z1=zl1(ep)
1551 z2=z1*z1
1552 vcore(ep,3,1)=z1
1553 vcore(ep,3,2)=-z1
1554 vcore(ep,3,3)=z1
1555 vcore(ep,3,4)=-z1
1556 corel(1,1)=vcore(ep,1,1)
1557 corel(1,2)=vcore(ep,1,2)
1558 corel(1,3)=vcore(ep,1,3)
1559 corel(1,4)=vcore(ep,1,4)
1560 corel(2,1)=vcore(ep,2,1)
1561 corel(2,2)=vcore(ep,2,2)
1562 corel(2,3)=vcore(ep,2,3)
1563 corel(2,4)=vcore(ep,2,4)
1564 x13 =x13_t(ep)
1565 x24 =x24_t(ep)
1566 y13 =y13_t(ep)
1567 y24 =y24_t(ep)
1568cc
1569 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
1570 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
1571 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
1572 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
1573 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
1574 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
1575 x13_t(ep) =x13*area_i(ep)
1576 x24_t(ep) =x24*area_i(ep)
1577 y13_t(ep) =y13*area_i(ep)
1578 y24_t(ep) =y24*area_i(ep)
1579C--------GAMA(2)
1580 gama1=-mx13*y24+my13*x24
1581 gama2= mx13*y13-my13*x13
1582C--------------------------
1583C CONSTRUIRE [F], [Q], [NM], [A-S], [AS] AUX NOEUDS
1584C--------------------------
1585 x21 =mx23-mx13
1586 x34 =(corel(1,3)-corel(1,4))*half
1587 y21 =my23-my13
1588 y34 =(corel(2,3)-corel(2,4))*half
1589 z21 = -z1
1590 z34 = z1
1591 l12 = sqrt(x21*x21+y21*y21+z2)
1592 l34 = sqrt(x34*x34+y34*y34+z2)
1593C---------
1594 x41 =mx34-mx13
1595 x32 =(corel(1,3)-corel(1,2))*half
1596 y41 =my34-my13
1597 y32 =(corel(2,3)-corel(2,2))*half
1598 z41 = -z1
1599 z32 = z1
1600C----------
1601 a_4=area(ep)*fourth
1602C
1603C----------
1604C CALCUL DE [QG] NG=1,4
1605C----------
1606 a_4=a_4/pg
1607 jmx13=mx13*pg
1608 jmy13=my13*pg
1609 jmz13=z1*pg
1610 j2myz=jmz13*jmz13
1611C---------- NG =1----------
1612 sz2=a_4-gama1
1613 sz=z2*l24(ep)
1614 sl=sqrt(sz+sz2*sz2)
1615 jac(ep,1)=sl*pg
1616 sl=one/max(sl,em20)
1617 vqg(ep,7,1)=-z1*y24
1618 vqg(ep,8,1)= z1*x24
1619 vqg(ep,9,1)= sz2*sl
1620 vqg(ep,7,3)=-vqg(ep,7,1)
1621 vqg(ep,8,3)=-vqg(ep,8,1)
1622 vqg(ep,7,1)= vqg(ep,7,1)*sl
1623 vqg(ep,8,1)= vqg(ep,8,1)*sl
1624C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1625 g1x1=mx23-jmx13
1626 g1y1=my23-jmy13
1627 c1 = sqrt(g1x1*g1x1 + g1y1*g1y1 +j2myz )
1628 g2x1=mx34-jmx13
1629 g2y1=my34-jmy13
1630 c2 = sqrt(g2x1*g2x1 + g2y1*g2y1 +j2myz)
1631 vjfi(ep,1,1)=(g2y1*vqg(ep,9,1)+jmz13*vqg(ep,8,1))
1632 vjfi(ep,2,1)=(-jmz13*vqg(ep,7,1)-g2x1*vqg(ep,9,1))
1633 vjfi(ep,3,1)=(g2x1*vqg(ep,8,1)-g2y1*vqg(ep,7,1))
1634 vqg(ep,1,1)= g1x1*c2+vjfi(ep,1,1)*c1
1635 vqg(ep,2,1)= g1y1*c2+vjfi(ep,2,1)*c1
1636 vqg(ep,3,1)= -jmz13*c2+vjfi(ep,3,1)*c1
1637 sl=sl/pg
1638 vjfi(ep,1,1)=vjfi(ep,1,1)*sl
1639 vjfi(ep,2,1)=vjfi(ep,2,1)*sl
1640 vjfi(ep,3,1)=vjfi(ep,3,1)*sl
1641C
1642 vjfi(ep,4,1)=(-jmz13*vqg(ep,8,1)-g1y1*vqg(ep,9,1))*sl
1643 vjfi(ep,5,1)=(g1x1*vqg(ep,9,1)+jmz13*vqg(ep,7,1))*sl
1644 vjfi(ep,6,1)=(g1y1*vqg(ep,7,1)-g1x1*vqg(ep,8,1))*sl
1645C
1646 sl = sqrt(vqg(ep,1,1)*vqg(ep,1,1) + vqg(ep,2,1)*vqg(ep,2,1) + vqg(ep,3,1)*vqg(ep,3,1))
1647 IF ( sl/=zero) sl = one / sl
1648 vqg(ep,1,1) = vqg(ep,1,1)*sl
1649 vqg(ep,2,1) = vqg(ep,2,1)*sl
1650 vqg(ep,3,1) = vqg(ep,3,1)*sl
1651C
1652 vqg(ep,4,1)= vqg(ep,8,1)*vqg(ep,3,1)-vqg(ep,9,1)*vqg(ep,2,1)
1653 vqg(ep,5,1)= vqg(ep,9,1)*vqg(ep,1,1)-vqg(ep,7,1)*vqg(ep,3,1)
1654 vqg(ep,6,1)= vqg(ep,7,1)*vqg(ep,2,1)-vqg(ep,8,1)*vqg(ep,1,1)
1655C---------- NG =3----------
1656 sz2=a_4+gama1
1657 sl=sqrt(sz+sz2*sz2)
1658 jac(ep,3)=sl*pg
1659 sl=one/max(sl,em20)
1660 vqg(ep,7,3)= vqg(ep,7,3)*sl
1661 vqg(ep,8,3)= vqg(ep,8,3)*sl
1662 vqg(ep,9,3)= sz2*sl
1663C
1664 g1x3=mx23+jmx13
1665 g1y3=my23+jmy13
1666 j1 = sqrt(g1x3*g1x3 + g1y3*g1y3 +j2myz )
1667C--------G2X3=G2X2------
1668 g2x2=mx34+jmx13
1669 g2y2=my34+jmy13
1670 j2 = sqrt(g2x2*g2x2 + g2y2*g2y2 +j2myz)
1671 vjfi(ep,1,3)=(g2y2*vqg(ep,9,3)-jmz13*vqg(ep,8,3))
1672 vjfi(ep,2,3)=(jmz13*vqg(ep,7,3)-g2x2*vqg(ep,9,3))
1673 vjfi(ep,3,3)=(g2x2*vqg(ep,8,3)-g2y2*vqg(ep,7,3))
1674 vqg(ep,1,3)= g1x3*j2+vjfi(ep,1,3)*j1
1675 vqg(ep,2,3)= g1y3*j2+vjfi(ep,2,3)*j1
1676 vqg(ep,3,3)= jmz13*j2+vjfi(ep,3,3)*j1
1677 sl=sl/pg
1678 vjfi(ep,1,3)=vjfi(ep,1,3)*sl
1679 vjfi(ep,2,3)=vjfi(ep,2,3)*sl
1680 vjfi(ep,3,3)=vjfi(ep,3,3)*sl
1681C
1682 vjfi(ep,4,3)=(jmz13*vqg(ep,8,3)-g1y3*vqg(ep,9,3))*sl
1683 vjfi(ep,5,3)=(g1x3*vqg(ep,9,3)-jmz13*vqg(ep,7,3))*sl
1684 vjfi(ep,6,3)=(g1y3*vqg(ep,7,3)-g1x3*vqg(ep,8,3))*sl
1685C
1686 sl = sqrt(vqg(ep,1,3)*vqg(ep,1,3) + vqg(ep,2,3)*vqg(ep,2,3) + vqg(ep,3,3)*vqg(ep,3,3))
1687 IF ( sl/=zero) sl = one / sl
1688 vqg(ep,1,3) = vqg(ep,1,3)*sl
1689 vqg(ep,2,3) = vqg(ep,2,3)*sl
1690 vqg(ep,3,3) = vqg(ep,3,3)*sl
1691C
1692 vqg(ep,4,3)= vqg(ep,8,3)*vqg(ep,3,3)-vqg(ep,9,3)*vqg(ep,2,3)
1693 vqg(ep,5,3)= vqg(ep,9,3)*vqg(ep,1,3)-vqg(ep,7,3)*vqg(ep,3,3)
1694 vqg(ep,6,3)= vqg(ep,7,3)*vqg(ep,2,3)-vqg(ep,8,3)*vqg(ep,1,3)
1695C---------- NG =2----------
1696 sz2=a_4+gama2
1697 sz=z2*l13(ep)
1698 sl=sqrt(sz+sz2*sz2)
1699 jac(ep,2)=sl*pg
1700 sl=one/max(sl,em20)
1701 vqg(ep,7,2)=-z1*y13
1702 vqg(ep,8,2)= z1*x13
1703 vqg(ep,9,2)= sz2*sl
1704 vqg(ep,7,4)=-vqg(ep,7,2)
1705 vqg(ep,8,4)=-vqg(ep,8,2)
1706 vqg(ep,7,2)= vqg(ep,7,2)*sl
1707 vqg(ep,8,2)= vqg(ep,8,2)*sl
1708C
1709 vjfi(ep,1,2)=(g2y2*vqg(ep,9,2)-jmz13*vqg(ep,8,2))
1710 vjfi(ep,2,2)=(jmz13*vqg(ep,7,2)-g2x2*vqg(ep,9,2))
1711 vjfi(ep,3,2)=(g2x2*vqg(ep,8,2)-g2y2*vqg(ep,7,2))
1712 vqg(ep,1,2)= g1x1*j2+vjfi(ep,1,2)*c1
1713 vqg(ep,2,2)= g1y1*j2+vjfi(ep,2,2)*c1
1714 vqg(ep,3,2)=-jmz13*j2+vjfi(ep,3,2)*c1
1715 sl=sl/pg
1716 vjfi(ep,1,2)=vjfi(ep,1,2)*sl
1717 vjfi(ep,2,2)=vjfi(ep,2,2)*sl
1718 vjfi(ep,3,2)=vjfi(ep,3,2)*sl
1719C
1720 vjfi(ep,4,2)=(-jmz13*vqg(ep,8,2)-g1y1*vqg(ep,9,2))*sl
1721 vjfi(ep,5,2)=(g1x1*vqg(ep,9,2)+jmz13*vqg(ep,7,2))*sl
1722 vjfi(ep,6,2)=(g1y1*vqg(ep,7,2)-g1x1*vqg(ep,8,2))*sl
1723C
1724 sl = sqrt(vqg(ep,1,2)*vqg(ep,1,2) + vqg(ep,2,2)*vqg(ep,2,2) + vqg(ep,3,2)*vqg(ep,3,2))
1725 IF ( sl/=zero) sl = one / sl
1726 vqg(ep,1,2) = vqg(ep,1,2)*sl
1727 vqg(ep,2,2) = vqg(ep,2,2)*sl
1728 vqg(ep,3,2) = vqg(ep,3,2)*sl
1729 vqg(ep,4,2)= vqg(ep,8,2)*vqg(ep,3,2)-vqg(ep,9,2)*vqg(ep,2,2)
1730 vqg(ep,5,2)= vqg(ep,9,2)*vqg(ep,1,2)-vqg(ep,7,2)*vqg(ep,3,2)
1731 vqg(ep,6,2)= vqg(ep,7,2)*vqg(ep,2,2)-vqg(ep,8,2)*vqg(ep,1,2)
1732C---------- NG =4----------
1733 sz2=a_4-gama2
1734 sl=sqrt(sz+sz2*sz2)
1735 jac(ep,4)=sl*pg
1736 sl=one/max(sl,em20)
1737 vqg(ep,7,4)= vqg(ep,7,4)*sl
1738 vqg(ep,8,4)= vqg(ep,8,4)*sl
1739 vqg(ep,9,4)= sz2*sl
1740C
1741 vjfi(ep,1,4)=(g2y1*vqg(ep,9,4)+jmz13*vqg(ep,8,4))
1742 vjfi(ep,2,4)=(-jmz13*vqg(ep,7,4)-g2x1*vqg(ep,9,4))
1743 vjfi(ep,3,4)=(g2x1*vqg(ep,8,4)-g2y1*vqg(ep,7,4))
1744 vqg(ep,1,4)= g1x3*c2+vjfi(ep,1,4)*j1
1745 vqg(ep,2,4)= g1y3*c2+vjfi(ep,2,4)*j1
1746 vqg(ep,3,4)=jmz13*c2+vjfi(ep,3,4)*j1
1747 sl=sl/pg
1748 vjfi(ep,1,4)=vjfi(ep,1,4)*sl
1749 vjfi(ep,2,4)=vjfi(ep,2,4)*sl
1750 vjfi(ep,3,4)=vjfi(ep,3,4)*sl
1751C
1752 vjfi(ep,4,4)=(jmz13*vqg(ep,8,4)-g1y3*vqg(ep,9,4))*sl
1753 vjfi(ep,5,4)=(g1x3*vqg(ep,9,4)-jmz13*vqg(ep,7,4))*sl
1754 vjfi(ep,6,4)=(g1y3*vqg(ep,7,4)-g1x3*vqg(ep,8,4))*sl
1755C
1756 sl = sqrt(vqg(ep,1,4)*vqg(ep,1,4) + vqg(ep,2,4)*vqg(ep,2,4) + vqg(ep,3,4)*vqg(ep,3,4))
1757 IF ( sl/=zero) sl = one / sl
1758 vqg(ep,1,4) = vqg(ep,1,4)*sl
1759 vqg(ep,2,4) = vqg(ep,2,4)*sl
1760 vqg(ep,3,4) = vqg(ep,3,4)*sl
1761 vqg(ep,4,4)= vqg(ep,8,4)*vqg(ep,3,4)-vqg(ep,9,4)*vqg(ep,2,4)
1762 vqg(ep,5,4)= vqg(ep,9,4)*vqg(ep,1,4)-vqg(ep,7,4)*vqg(ep,3,4)
1763 vqg(ep,6,4)= vqg(ep,7,4)*vqg(ep,2,4)-vqg(ep,8,4)*vqg(ep,1,4)
1764 area(ep)=jac(ep,1)+jac(ep,2)+jac(ep,3)+jac(ep,4)
1765 ENDDO
1766C--------------------------
1767C TRANSMET 6DDL-->5DDL, w/o projection
1768C--------------------------
1769 DO i=nplat+1,jlt
1770 ep =iplat(i)
1771 vxyz(ep,1,1)= half*v13(ep,1)+fourth*vhi(ep,1)
1772 vxyz(ep,1,2)= half*v24(ep,1)-fourth*vhi(ep,1)
1773 vxyz(ep,1,3)= -half*v13(ep,1)+fourth*vhi(ep,1)
1774 vxyz(ep,1,4)= -half*v24(ep,1)-fourth*vhi(ep,1)
1775
1776 vxyz(ep,2,1)= half*v13(ep,2)+fourth*vhi(ep,2)
1777 vxyz(ep,2,2)= half*v24(ep,2)-fourth*vhi(ep,2)
1778 vxyz(ep,2,3)= -half*v13(ep,2)+fourth*vhi(ep,2)
1779 vxyz(ep,2,4)= -half*v24(ep,2)-fourth*vhi(ep,2)
1780C
1781 vxyz(ep,3,1)= half*v13(ep,3)+fourth*vhi(ep,3)
1782 vxyz(ep,3,2)= half*v24(ep,3)-fourth*vhi(ep,3)
1783 vxyz(ep,3,3)= -half*v13(ep,3)+fourth*vhi(ep,3)
1784 vxyz(ep,3,4)= -half*v24(ep,3)-fourth*vhi(ep,3)
1785C--------drilling dof
1786 IF (isrot>0) THEN
1787 rr(1,1) =vq(ep,1,1)*axyz(ep,1,1)+vq(ep,2,1)*axyz(ep,2,1)+vq(ep,3,1)*axyz(ep,3,1)
1788 rr(1,2) =vq(ep,1,1)*axyz(ep,1,2)+vq(ep,2,1)*axyz(ep,2,2)+vq(ep,3,1)*axyz(ep,3,2)
1789 rr(1,3) =vq(ep,1,1)*axyz(ep,1,3)+vq(ep,2,1)*axyz(ep,2,3)+vq(ep,3,1)*axyz(ep,3,3)
1790 rr(1,4) =vq(ep,1,1)*axyz(ep,1,4)+vq(ep,2,1)*axyz(ep,2,4)+vq(ep,3,1)*axyz(ep,3,4)
1791 rr(2,1) =vq(ep,1,2)*axyz(ep,1,1)+vq(ep,2,2)*axyz(ep,2,1)+vq(ep,3,2)*axyz(ep,3,1)
1792 rr(2,2) =vq(ep,1,2)*axyz(ep,1,2)+vq(ep,2,2)*axyz(ep,2,2)+vq(ep,3,2)*axyz(ep,3,2)
1793 rr(2,3) =vq(ep,1,2)*axyz(ep,1,3)+vq(ep,2,2)*axyz(ep,2,3)+vq(ep,3,2)*axyz(ep,3,3)
1794 rr(2,4) =vq(ep,1,2)*axyz(ep,1,4)+vq(ep,2,2)*axyz(ep,2,4)+vq(ep,3,2)*axyz(ep,3,4)
1795 rr(3,1) =vq(ep,1,3)*axyz(ep,1,1)+vq(ep,2,3)*axyz(ep,2,1)+vq(ep,3,3)*axyz(ep,3,1)
1796 rr(3,2) =vq(ep,1,3)*axyz(ep,1,2)+vq(ep,2,3)*axyz(ep,2,2)+vq(ep,3,3)*axyz(ep,3,2)
1797 rr(3,3) =vq(ep,1,3)*axyz(ep,1,3)+vq(ep,2,3)*axyz(ep,2,3)+vq(ep,3,3)*axyz(ep,3,3)
1798 rr(3,4) =vq(ep,1,3)*axyz(ep,1,4)+vq(ep,2,3)*axyz(ep,2,4)+vq(ep,3,3)*axyz(ep,3,4)
1799C----------not right takes the RLZ with cbacoor
1800 rlz(ep,1)=(vqn(ep,7,1)*rr(1,1)+vqn(ep,8,1)*rr(2,1)+vqn(ep,9,1)*rr(3,1))
1801 rlz(ep,2)=(vqn(ep,7,2)*rr(1,2)+vqn(ep,8,2)*rr(2,2)+vqn(ep,9,2)*rr(3,2))
1802 rlz(ep,3)=(vqn(ep,7,3)*rr(1,3)+vqn(ep,8,3)*rr(2,3)+vqn(ep,9,3)*rr(3,3))
1803 rlz(ep,4)=(vqn(ep,7,4)*rr(1,4)+vqn(ep,8,4)*rr(2,4)+vqn(ep,9,4)*rr(3,4))
1804 END IF !(ISROT>0) THEN
1805 END DO
1806C
1807 off_l = zero
1808 DO ep=jft,jlt
1809 off_l = min(off_l,offg(ep))
1810 ENDDO
1811C
1812 IF(off_l<zero)THEN
1813 DO ep=jft,jlt
1814 IF(offg(ep)<zero)THEN
1815 vxyz(ep,1:3,1:4)= zero
1816 rlz(ep,1:4)= zero
1817 ENDIF
1818 ENDDO
1819 ENDIF
1820C-------------
1821 RETURN