OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cbacoor.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| cbacoor ../engine/source/elements/shell/coqueba/cbacoor.F
25!||--- called by ------------------------------------------------------
26!|| cbaforc3 ../engine/source/elements/shell/coqueba/cbaforc3.F
27!||--- calls -----------------------------------------------------
28!|| clskew3 ../engine/source/elements/sh3n/coquedk/cdkcoor3.F
29!|| cortdir3 ../engine/source/elements/shell/coque/cortdir3.F
30!||--- uses -----------------------------------------------------
31!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!|| plyxfem_mod ../engine/share/modules/plyxfem_mod.F
34!||====================================================================
35 SUBROUTINE cbacoor(ELBUF_STR,JFT,JLT,X,V,
36 . VR,IXC,PM,OFFG,LL,
37 1 AREA,VXYZ, RXYZ,VCORE,JAC,HX,HY,VKSI,VETA,
38 2 VQN,VQG,VQ,VJFI,VNRM,VASTN,NPLAT,IPLAT,
39 3 X13_T ,X24_T ,Y13_T,Y24_T,OFF,DI,NLAY,
40 4 IREP ,NPT ,ISMSTR,NEL ,ISROT ,
41 5 SMSTR ,DIR_A ,DIR_B ,FACN ,ZL1 ,
42 6 R11 ,R12 ,R13 ,R21 ,R22 ,R23 ,
43 7 R31 ,R32 ,R33 ,INOD ,RLZ ,
44 8 THK ,IPLYCXFEM ,UX1 ,UX2 ,UX3 ,
45 9 UX4 ,UY1 ,UY2 ,UY3 ,UY4 ,
46 A VL1 ,VL2 ,VL3 ,VL4 ,XL2 ,
47 B XL3 ,XL4 ,YL2 ,YL3 ,YL4 ,XLCOR ,NPINCH)
48C-----------------------------------------------
49C M o d u l e s
50C-----------------------------------------------
51 USE elbufdef_mod
52 USE plyxfem_mod
53 use element_mod , only : nixc
54C-----------------------------------------------
55C I M P L I C I T T Y P E S
56C-----------------------------------------------
57#include "implicit_f.inc"
58#include "comlock.inc"
59c-----------------------------------------------
60c g l o b a l p a r a m e t e r s
61c-----------------------------------------------
62#include "mvsiz_p.inc"
63#include "param_c.inc"
64#include "scr17_c.inc"
65#include "com08_c.inc"
66#include "impl1_c.inc"
67#include "scr14_c.inc"
68C-----------------------------------------------
69C D U M M Y A R G U M E N T S
70C-----------------------------------------------
71 INTEGER JFT,JLT,NNOD,NPLAT,IREP,NPT,NLAY,ISMSTR,ISROT,IPLYCXFEM,NEL,NPINCH
72 INTEGER IXC(NIXC,*),IPLAT(*),INOD(*)
73 PARAMETER (NNOD = 4)
74 my_real x(3,*), v(3,*), vr(3,*),pm(npropm,*),offg(*),off(*),
75 . vcore(mvsiz,3,nnod),vxyz(mvsiz,3,nnod),area(*),vjfi(mvsiz,6,4),
76 . vqn(mvsiz,9,nnod),vqg(mvsiz,9,nnod),vq(mvsiz,3,3),rxyz(mvsiz,2,nnod),
77 . ll(*),vastn(mvsiz,4,nnod),vnrm(mvsiz,3,nnod),
78 . jac(mvsiz,4),hx(mvsiz,4),hy(mvsiz,4),veta(4,4),vksi(4,4),y24_t(*),
79 . x13_t(*),x24_t(*),y13_t(*),off_l, di(mvsiz,6),
80 . dir_a(nel,*),dir_b(nel,*),facn(mvsiz,2),
81 . zl1(mvsiz),r11(mvsiz),r12(mvsiz),r13(mvsiz),
82 . r21(mvsiz),r22(mvsiz),r23(mvsiz),r31(mvsiz),r32(mvsiz),
83 . r33(mvsiz),rlz(mvsiz,4),thk(mvsiz),
84 . ux1(*),ux2(*),ux3(*),ux4(*),uy1(*),uy2(*),uy3(*),uy4(*),
85 . vl1(mvsiz,3),vl2(mvsiz,3),vl3(mvsiz,3),vl4(mvsiz,3),
86 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),yl3(mvsiz),yl4(mvsiz),
87 . xlcor(mvsiz,2,3)
88 double precision
89 . smstr(*)
90C
91 TYPE(elbuf_struct_) :: ELBUF_STR
92C-----------------------------------------------
93c FUNCTION: preliminary compute for QBAT
94c
95c Note:
96c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
97c
98c TYPE NAME FUNCTION
99c I JFT,JLT - element id limit
100c I X, V, VR(3,NUMNOD) - coordinate,velocity, rotational velocity in global system
101c I IXC(NIXC,*NEL) - element connectivity of shell (and other data)
102c I PM(NPROPM,MID) - input Material data
103c I OFFG(NEL),OFF(NEL) - element activation flag, local flag
104c O LL(NEL),FACN(2,NEL)- characteristic length (for dt compute)
105c O AREA(NEL),NTP - Area, num. of integrating points in thickness
106c O VCORE(MVSIZ,3,4) - coordinates of 4 nodes in local system
107c BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
108c for plat element
109c O VXYZ(3,4,NEL) - nodal velocity (4 nodes) in local system
110C V13->VXYZ(,1,),V24->VXYZ(,2,),VHI->VXYZ(,3,)------------
111c for plat element
112c O RXYZ(NEL,2,4) - nodal rotational velocity (4 nodes) in local system 5 dof
113C R13->RXYZ(,1,),R24->RXYZ(,2,),RHI->RXYZ(,3,),RTI->RXYZ(,4,)----
114c for plat element
115c O JAC,HX,HY(4,NEL) - Jacobian,asymmetric part(hourglass) at 4 Gauss points
116c O VKSI,VETA(4,4) - derivation of shape functions------
117c O VQ(NEL,9),VQN,VQG(NEL,9,4) - local system of element, nodal(4) and Gauss points(4)
118c O VJFI,VNRM,VASTN - terms used for assumed strains
119c O NPLAT,IPLAT(NEL) - num. of plat element and indice
120c O Xij_T,Yij_T(NEL) - [B] components used for in-plane shear assumed strain
121c IO XiS,YiS,ZiS(NEL) - reference configurations used for small strain options
122c I ISMSTR,IREP - small strain and orthotrop flags
123c O DIR_A,DIR_B(2,NEL) - actual orthotropic directions
124c O ZL1(NEL) - warped length
125c O Rij(NEL) - local system base
126c I IDRIL - drilling dof flag
127c O RLZ(4,NEL) - nodal Rz (rotational velocity) in local system
128c (used only Idril active )
129c I INOD(NUMNOD) - Node numbers for PLYXFEM
130C-----------------------------------------------
131C L O C A L V A R I A B L E S
132C-----------------------------------------------
133 INTEGER I, II(6), J, K, L, M, I1, I2, I3, I4, EP, SPLAT
134 INTEGER NN(4),JPLAT(MVSIZ)
135 MY_REAL
136 . XX2,XX3,XX4,YY1,YY2,YY3,YY4,ZZ1,ZZ2,ZZ3,ZZ4
137 MY_REAL
138 . C1,C2,L13(MVSIZ),L24(MVSIZ),
139 . RR(3,NNOD)
140 MY_REAL
141 . TOL,PG,J0,J1,J2,DETA,COREL(3,4),S1
142 MY_REAL
143 . X13,X24,Y13,MX13,MX23,MX34,MY13,Z1,Z2,GAMA1,GAMA2,
144 . X21,X34,Y21,Y34,Z21,Z34,X41,X32,Y41,Y32,Z41,Z32,L12,L34,
145 . a_4,sl,sz2,sz,jmx13,jmy13,jmz13,j2myz,lm(mvsiz),y24,
146 . my23,my34,scal,g1x1,g1x3,g1y1,g1y3,g2x1,g2x2,g2y1,g2y2,
147 . ar(3),ad(nnod),btb(6),xx1,yy,zz,xy,xz1,yz,d(6),
148 . alr(3),ald(nnod),dbad(3),abc,xxyz2,zzxy2,yyxz2
149 my_real
150 . lxyz0(3),deta1(mvsiz),rx(mvsiz), ry(mvsiz), rz(mvsiz),
151 . sx(mvsiz),sy(mvsiz),
152 . area_i(mvsiz),ssz(mvsiz)
153C
154 my_real
155 . vg13(3),vg24(3),vghi(3),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
156 . faci,fac1,fac2,
157 . dt05,dt025,exz,eyz,ddrx,ddry,v13x,v24x,vhix,ddrz1,ddrz2,ddrz,
158 . vnx1, vny1, vnz1,vnx2,vny2,vnz2,vnx3,vny3,vnz3,vnx4,vny4,
159 . norm
160 DATA pg/.577350269189626/
161C=======================================================================
162 DO i=1,6
163 ii(i) = nel*(i-1)
164 ENDDO
165C
166 tol=em12
167 IF (isrot > 0 ) tol=em8
168C
169 vksi(1,1)=-fourth*(one+pg)
170 vksi(2,1)=-vksi(1,1)
171 vksi(3,1)= fourth*(one-pg)
172 vksi(4,1)=-vksi(3,1)
173 veta(1,1)=-fourth*(one+pg)
174 veta(2,1)=-fourth*(one-pg)
175 veta(3,1)=-veta(2,1)
176 veta(4,1)=-veta(1,1)
177 vksi(1,2)= vksi(1,1)
178 vksi(2,2)=-vksi(1,2)
179 vksi(3,2)= vksi(3,1)
180 vksi(4,2)=-vksi(3,2)
181 veta(1,2)= veta(2,1)
182 veta(2,2)= veta(1,1)
183 veta(3,2)=-veta(2,2)
184 veta(4,2)=-veta(1,2)
185 vksi(1,3)=-vksi(3,1)
186 vksi(2,3)=-vksi(1,3)
187 vksi(3,3)=-vksi(1,1)
188 vksi(4,3)=-vksi(3,3)
189 veta(1,3)= veta(1,2)
190 veta(2,3)= veta(2,2)
191 veta(3,3)=-veta(2,3)
192 veta(4,3)=-veta(1,3)
193 vksi(1,4)= vksi(1,3)
194 vksi(2,4)=-vksi(1,4)
195 vksi(3,4)= vksi(3,3)
196 vksi(4,4)=-vksi(3,4)
197 veta(1,4)= veta(1,1)
198 veta(2,4)= veta(2,1)
199 veta(3,4)=-veta(2,4)
200 veta(4,4)=-veta(1,4)
201C
202 DO i=jft,jlt
203 rx(i)=x(1,ixc(3,i))+x(1,ixc(4,i))-x(1,ixc(2,i))-x(1,ixc(5,i))
204 sx(i)=x(1,ixc(4,i))+x(1,ixc(5,i))-x(1,ixc(2,i))-x(1,ixc(3,i))
205 ry(i)=x(2,ixc(3,i))+x(2,ixc(4,i))-x(2,ixc(2,i))-x(2,ixc(5,i))
206 sy(i)=x(2,ixc(4,i))+x(2,ixc(5,i))-x(2,ixc(2,i))-x(2,ixc(3,i))
207 rz(i)=x(3,ixc(3,i))+x(3,ixc(4,i))-x(3,ixc(2,i))-x(3,ixc(5,i))
208 ssz(i)=x(3,ixc(4,i))+x(3,ixc(5,i))-x(3,ixc(2,i))-x(3,ixc(3,i))
209 ENDDO
210C----------------------------
211C LOCAL SYSTEM
212C----------------------------
213 k = 0
214 CALL clskew3(jft,jlt,k,
215 . rx, ry, rz,
216 . sx, sy, ssz,
217 . r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
218 DO i=jft,jlt
219 area(i)=fourth*deta1(i)
220 area_i(i)=one/area(i)
221 vq(i,1,1)=r11(i)
222 vq(i,2,1)=r21(i)
223 vq(i,3,1)=r31(i)
224 vq(i,1,2)=r12(i)
225 vq(i,2,2)=r22(i)
226 vq(i,3,2)=r32(i)
227 vq(i,1,3)=r13(i)
228 vq(i,2,3)=r23(i)
229 vq(i,3,3)=r33(i)
230 ENDDO
231c
232 DO i=jft,jlt
233
234 j=ixc(2,i)
235 k=ixc(3,i)
236 l=ixc(4,i)
237 m=ixc(5,i)
238
239 lxyz0(1)=fourth*(x(1,l)+x(1,m)+x(1,j)+x(1,k))
240 lxyz0(2)=fourth*(x(2,l)+x(2,m)+x(2,j)+x(2,k))
241 lxyz0(3)=fourth*(x(3,l)+x(3,m)+x(3,j)+x(3,k))
242
243 xx1=x(1,k)-x(1,j)
244 yy1=x(2,k)-x(2,j)
245 zz1=x(3,k)-x(3,j)
246
247 xl2(i)=r11(i)*xx1+r21(i)*yy1+r31(i)*zz1
248 yl2(i)=r12(i)*xx1+r22(i)*yy1+r32(i)*zz1
249
250 xx2=x(1,j)-lxyz0(1)
251 yy2=x(2,j)-lxyz0(2)
252 zz2=x(3,j)-lxyz0(3)
253 zl1(i)=r13(i)*xx2+r23(i)*yy2+r33(i)*zz2
254
255 xx3=x(1,l)-x(1,j)
256 yy3=x(2,l)-x(2,j)
257 zz3=x(3,l)-x(3,j)
258 xl3(i)=r11(i)*xx3+r21(i)*yy3+r31(i)*zz3
259 yl3(i)=r12(i)*xx3+r22(i)*yy3+r32(i)*zz3
260
261 xx4=x(1,m)-x(1,j)
262 yy4=x(2,m)-x(2,j)
263 zz4=x(3,m)-x(3,j)
264 xl4(i)=r11(i)*xx4+r21(i)*yy4+r31(i)*zz4
265 yl4(i)=r12(i)*xx4+r22(i)*yy4+r32(i)*zz4
266 ENDDO
267C----------------------------
268C SMALL STRAIN
269C----------------------------
270 IF(ismstr==10)THEN
271 DO i=jft,jlt
272 xlcor(i,1,1)=xl2(i)
273 xlcor(i,2,1)=yl2(i)
274 xlcor(i,1,2)=xl3(i)
275 xlcor(i,2,2)=yl3(i)
276 xlcor(i,1,3)=xl4(i)
277 xlcor(i,2,3)=yl4(i)
278 ENDDO
279 ELSEIF(ismstr==11)THEN
280 DO i=jft,jlt
281 IF(abs(offg(i))==one)offg(i)=sign(two,offg(i))
282 ux1(i) = zero
283 uy1(i) = zero
284 ux2(i) = zero
285 uy2(i) = zero
286 ux3(i) = zero
287 uy3(i) = zero
288 ux4(i) = zero
289 uy4(i) = zero
290 IF(abs(offg(i))==two)THEN
291 ux2(i) = xl2(i)-smstr(ii(1)+i)
292 uy2(i) = yl2(i)-smstr(ii(2)+i)
293 ux3(i) = xl3(i)-smstr(ii(3)+i)
294 uy3(i) = yl3(i)-smstr(ii(4)+i)
295 ux4(i) = xl4(i)-smstr(ii(5)+i)
296 uy4(i) = yl4(i)-smstr(ii(6)+i)
297 xl2(i) = smstr(ii(1)+i)
298 yl2(i) = smstr(ii(2)+i)
299 xl3(i) = smstr(ii(3)+i)
300 yl3(i) = smstr(ii(4)+i)
301 xl4(i) = smstr(ii(5)+i)
302 yl4(i) = smstr(ii(6)+i)
303 zl1(i) = zero
304 area(i)=half*((xl2(i)-xl4(i))*yl3(i)-xl3(i)*(yl2(i)-yl4(i)))
305 area_i(i)=one/max(em20,area(i))
306 ELSE
307 smstr(ii(1)+i)=xl2(i)
308 smstr(ii(2)+i)=yl2(i)
309 smstr(ii(3)+i)=xl3(i)
310 smstr(ii(4)+i)=yl3(i)
311 smstr(ii(5)+i)=xl4(i)
312 smstr(ii(6)+i)=yl4(i)
313 ENDIF
314 ENDDO
315 ELSEIF(ismstr==1.OR.ismstr==2)THEN
316 DO i=jft,jlt
317 IF(abs(offg(i))==two)THEN
318 xl2(i)=smstr(ii(1)+i)
319 yl2(i)=smstr(ii(2)+i)
320 xl3(i)=smstr(ii(3)+i)
321 yl3(i)=smstr(ii(4)+i)
322 xl4(i)=smstr(ii(5)+i)
323 yl4(i)=smstr(ii(6)+i)
324 zl1(i)=zero
325 area(i)=half*((xl2(i)-xl4(i))*yl3(i)-xl3(i)*(yl2(i)-yl4(i)))
326 area_i(i)=one/max(em20,area(i))
327 ELSE
328 smstr(ii(1)+i)=xl2(i)
329 smstr(ii(2)+i)=yl2(i)
330 smstr(ii(3)+i)=xl3(i)
331 smstr(ii(4)+i)=yl3(i)
332 smstr(ii(5)+i)=xl4(i)
333 smstr(ii(6)+i)=yl4(i)
334 ENDIF
335 ENDDO
336 ENDIF
337 IF(ismstr==1)THEN
338 DO i=jft,jlt
339 IF(offg(i)==one)offg(i)=two
340 ENDDO
341 ENDIF
342C----------------------------
343C ORTHOTROPY
344C----------------------------
345 CALL cortdir3(elbuf_str,dir_a ,dir_b ,jft ,jlt ,
346 . nlay ,irep ,rx ,ry ,rz ,
347 . sx ,sy ,ssz ,r11 ,r21 ,
348 . r31 ,r12 ,r22 ,r32 ,nel )
349C----------------------------
350 DO ep=jft,jlt
351 lxyz0(1)=fourth*(xl2(ep)+xl3(ep)+xl4(ep))
352 lxyz0(2)=fourth*(yl2(ep)+yl3(ep)+yl4(ep))
353 vcore(ep,1,1)=-lxyz0(1)
354 vcore(ep,1,2)=xl2(ep)-lxyz0(1)
355 vcore(ep,1,3)=xl3(ep)-lxyz0(1)
356 vcore(ep,1,4)=xl4(ep)-lxyz0(1)
357 vcore(ep,2,1)=-lxyz0(2)
358 vcore(ep,2,2)=yl2(ep)-lxyz0(2)
359 vcore(ep,2,3)=yl3(ep)-lxyz0(2)
360 vcore(ep,2,4)=yl4(ep)-lxyz0(2)
361C
362 x13_t(ep) =(vcore(ep,1,1)-vcore(ep,1,3))*half
363 x24_t(ep) =(vcore(ep,1,2)-vcore(ep,1,4))*half
364 y13_t(ep) =(vcore(ep,2,1)-vcore(ep,2,3))*half
365 y24_t(ep) =(vcore(ep,2,2)-vcore(ep,2,4))*half
366 l13(ep)=x13_t(ep)*x13_t(ep)+y13_t(ep)*y13_t(ep)
367 l24(ep)=x24_t(ep)*x24_t(ep)+y24_t(ep)*y24_t(ep)
368 ll(ep)=max(l13(ep),l24(ep))
369 c1 =vcore(ep,1,2)*vcore(ep,2,4)-vcore(ep,2,2)*vcore(ep,1,4)
370 c2 =vcore(ep,1,1)*vcore(ep,2,3)-vcore(ep,2,1)*vcore(ep,1,3)
371 lm(ep) =max(abs(c1),abs(c2))
372 ENDDO
373c
374 DO ep=jft,jlt
375
376 nn(1)=ixc(2,ep)
377 nn(2)=ixc(3,ep)
378 nn(3)=ixc(4,ep)
379 nn(4)=ixc(5,ep)
380
381 vl1(ep,1)=v(1,nn(1))
382 vl1(ep,2)=v(2,nn(1))
383 vl1(ep,3)=v(3,nn(1))
384 vl2(ep,1)=v(1,nn(2))
385 vl2(ep,2)=v(2,nn(2))
386 vl2(ep,3)=v(3,nn(2))
387 vl3(ep,1)=v(1,nn(3))
388 vl3(ep,2)=v(2,nn(3))
389 vl3(ep,3)=v(3,nn(3))
390 vl4(ep,1)=v(1,nn(4))
391 vl4(ep,2)=v(2,nn(4))
392 vl4(ep,3)=v(3,nn(4))
393
394 vg13(1)=vl1(ep,1)-vl3(ep,1)
395 vg24(1)=vl2(ep,1)-vl4(ep,1)
396 vghi(1)=vl1(ep,1)-vl2(ep,1)+vl3(ep,1)-vl4(ep,1)
397C
398 vg13(2)=vl1(ep,2)-vl3(ep,2)
399 vg24(2)=vl2(ep,2)-vl4(ep,2)
400 vghi(2)=vl1(ep,2)-vl2(ep,2)+vl3(ep,2)-vl4(ep,2)
401C
402 vg13(3)=vl1(ep,3)-vl3(ep,3)
403 vg24(3)=vl2(ep,3)-vl4(ep,3)
404 vghi(3)=vl1(ep,3)-vl2(ep,3)+vl3(ep,3)-vl4(ep,3)
405C
406 v13(ep,1)=(vq(ep,1,1)*vg13(1)+vq(ep,2,1)*vg13(2)+vq(ep,3,1)*vg13(3))
407 v24(ep,1)=(vq(ep,1,1)*vg24(1)+vq(ep,2,1)*vg24(2)+vq(ep,3,1)*vg24(3))
408 vhi(ep,1)=(vq(ep,1,1)*vghi(1)+vq(ep,2,1)*vghi(2)+vq(ep,3,1)*vghi(3))
409 v13(ep,2)=(vq(ep,1,2)*vg13(1)+vq(ep,2,2)*vg13(2)+vq(ep,3,2)*vg13(3))
410 v24(ep,2)=(vq(ep,1,2)*vg24(1)+vq(ep,2,2)*vg24(2)+vq(ep,3,2)*vg24(3))
411 vhi(ep,2)=(vq(ep,1,2)*vghi(1)+vq(ep,2,2)*vghi(2)+vq(ep,3,2)*vghi(3))
412 v13(ep,3)=(vq(ep,1,3)*vg13(1)+vq(ep,2,3)*vg13(2)+vq(ep,3,3)*vg13(3))
413 v24(ep,3)=(vq(ep,1,3)*vg24(1)+vq(ep,2,3)*vg24(2)+vq(ep,3,3)*vg24(3))
414 vhi(ep,3)=(vq(ep,1,3)*vghi(1)+vq(ep,2,3)*vghi(2)+vq(ep,3,3)*vghi(3))
415 ENDDO
416C
417 IF (impl_s==0.OR.imp_lr>0) THEN
418 dt05 = half*dt1
419 dt025 =fourth*dt1
420 DO i=jft,jlt
421 exz = y24_t(i)*v13(i,3)-y13_t(i)*v24(i,3)
422 eyz = -x24_t(i)*v13(i,3)+x13_t(i)*v24(i,3)
423 ddry=dt05*exz*area_i(i)
424 ddrx=dt05*eyz*area_i(i)
425 v13x = v13(i,1)
426 v24x = v24(i,1)
427 vhix = vhi(i,1)
428 ddrz1=dt025*(v13(i,2)-v24(i,2))/(x13_t(i)-x24_t(i))
429 IF (abs(x13_t(i)-x24_t(i))<em10) ddrz1 = zero
430 v13(i,1) = v13(i,1)-ddry*v13(i,3)-ddrz1*v13(i,2)
431 v24(i,1) = v24(i,1)-ddry*v24(i,3)-ddrz1*v24(i,2)
432 vhi(i,1) = vhi(i,1)-ddry*vhi(i,3)-ddrz1*vhi(i,2)
433 ddrz2=dt025*(v13x+v24x)/(y13_t(i)+y24_t(i))
434 IF (abs(y13_t(i)+y24_t(i))<em10) ddrz2 = zero
435 v13(i,2) = v13(i,2)-ddrx*v13(i,3)-ddrz2*v13x
436 v24(i,2) = v24(i,2)-ddrx*v24(i,3)-ddrz2*v24x
437 vhi(i,2) = vhi(i,2)-ddrx*vhi(i,3)-ddrz2*vhix
438 ENDDO
439 ENDIF
440 nplat=jft-1
441 splat= 0
442C [PM] in order to make all shells as non-flat shells for pinching
443C [PM] flat shells algorithm not implemented for pinching
444 DO ep=jft,jlt
445 z2=zl1(ep)*zl1(ep)
446 IF ((z2<tol*ll(ep).OR.npt==1).AND.npinch==0) THEN
447 nplat=nplat+1
448 iplat(nplat)=ep
449 ELSE
450 splat=splat+1
451 jplat(splat)=ep
452 ENDIF
453 ENDDO
454 DO ep=nplat+1,jlt
455 iplat(ep)=jplat(ep-nplat)
456 ENDDO
457#include "vectorize.inc"
458 DO i=jft,nplat
459 ep =iplat(i)
460 x13 =x13_t(ep)
461 x24 =x24_t(ep)
462 y13 =y13_t(ep)
463 y24 =y24_t(ep)
464C
465 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
466 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
467 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
468 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
469 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
470 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
471 x13_t(ep) =x13*area_i(ep)
472 x24_t(ep) =x24*area_i(ep)
473 y13_t(ep) =y13*area_i(ep)
474 y24_t(ep) =y24*area_i(ep)
475C--------GAMA(2)
476 gama1=-mx13*y24+my13*x24
477 gama2= mx13*y13-my13*x13
478 z2=zl1(ep)*zl1(ep)
479 nn(1)=ixc(2,ep)
480 nn(2)=ixc(3,ep)
481 nn(3)=ixc(4,ep)
482 nn(4)=ixc(5,ep)
483C-------TRANSMET 3DDL ROTATIONS ----
484 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))
485 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))
486 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))
487 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))
488 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))
489 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))
490 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))
491 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))
492C-------V13(I)->VXYZ(I,1),V24(I)->VXYZ(I,2),VHI(I)->VXYZ(I,3)------------
493 vxyz(ep,1,1)=v13(ep,1)
494 vxyz(ep,2,1)=v13(ep,2)
495 vxyz(ep,3,1)=v13(ep,3)
496C
497 vxyz(ep,1,2)=v24(ep,1)
498 vxyz(ep,2,2)=v24(ep,2)
499 vxyz(ep,3,2)=v24(ep,3)
500C
501 vxyz(ep,1,3)=vhi(ep,1)
502 vxyz(ep,2,3)=vhi(ep,2)
503 vxyz(ep,3,3)=vhi(ep,3)
504C-------R13(I)->RXYZ(I,1),R24(I)->RXYZ(I,2),RHI(I)->RXYZ(I,3),RTI(I)->RXYZ(I,4)------------
505 rxyz(ep,1,1)=rr(1,1)-rr(1,3)
506 rxyz(ep,2,1)=rr(2,1)-rr(2,3)
507 rxyz(ep,1,2)=rr(1,2)-rr(1,4)
508 rxyz(ep,2,2)=rr(2,2)-rr(2,4)
509 rxyz(ep,1,3)=rr(1,1)-rr(1,2)+rr(1,3)-rr(1,4)
510 rxyz(ep,2,3)=rr(2,1)-rr(2,2)+rr(2,3)-rr(2,4)
511 rxyz(ep,1,4)=rr(1,1)+rr(1,2)+rr(1,3)+rr(1,4)
512 rxyz(ep,2,4)=rr(2,1)+rr(2,2)+rr(2,3)+rr(2,4)
513C-------BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
514C-------ARE STORED IN VCORE---
515 vcore(ep,1,1)=y24_t(ep)
516 vcore(ep,2,1)=-y13_t(ep)
517 vcore(ep,3,1)=-x24_t(ep)
518 vcore(ep,1,2)= x13_t(ep)
519 vcore(ep,2,2)=gama1*area_i(ep)
520 vcore(ep,3,2)=gama2*area_i(ep)
521 vcore(ep,1,3)=mx23
522 vcore(ep,2,3)=my23
523 vcore(ep,3,3)=mx34
524 vcore(ep,1,4)=my34
525 vcore(ep,2,4)=mx13
526 vcore(ep,3,4)=my13
527C
528 j1=(mx23*my13-mx13*my23)*pg
529 j2=-(mx13*my34-mx34*my13)*pg
530 j0=area(ep)*fourth
531C---- JAC(jacobiens aux points d'integration) peut etre negative------
532 jac(ep,1)=abs(j0+j2-j1)
533 jac(ep,2)=abs(j0+j2+j1)
534 jac(ep,3)=abs(j0-j2+j1)
535 jac(ep,4)=abs(j0-j2-j1)
536 j1=(my23-my34)*pg
537 j2=-(my23+my34)*pg
538 hx(ep,1)=j1/jac(ep,1)
539 hx(ep,2)=j2/jac(ep,2)
540 hx(ep,3)=-j1/jac(ep,3)
541 hx(ep,4)=-j2/jac(ep,4)
542 j1=(mx34-mx23)*pg
543 j2=(mx34+mx23)*pg
544 hy(ep,1)=j1/jac(ep,1)
545 hy(ep,2)=j2/jac(ep,2)
546 hy(ep,3)=-j1/jac(ep,3)
547 hy(ep,4)=-j2/jac(ep,4)
548C
549 ENDDO
550#include "vectorize.inc"
551 DO i=nplat+1,jlt
552 ep =iplat(i)
553 z1=zl1(ep)
554 z2=z1*z1
555 vcore(ep,3,1)=z1
556 vcore(ep,3,2)=-z1
557 vcore(ep,3,3)=z1
558 vcore(ep,3,4)=-z1
559 corel(1,1)=vcore(ep,1,1)
560 corel(1,2)=vcore(ep,1,2)
561 corel(1,3)=vcore(ep,1,3)
562 corel(1,4)=vcore(ep,1,4)
563 corel(2,1)=vcore(ep,2,1)
564 corel(2,2)=vcore(ep,2,2)
565 corel(2,3)=vcore(ep,2,3)
566 corel(2,4)=vcore(ep,2,4)
567 x13 =x13_t(ep)
568 x24 =x24_t(ep)
569 y13 =y13_t(ep)
570 y24 =y24_t(ep)
571cc
572 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
573 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
574 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
575 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
576 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
577 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
578 x13_t(ep) =x13*area_i(ep)
579 x24_t(ep) =x24*area_i(ep)
580 y13_t(ep) =y13*area_i(ep)
581 y24_t(ep) =y24*area_i(ep)
582C--------GAMA(2)
583 gama1=-mx13*y24+my13*x24
584 gama2= mx13*y13-my13*x13
585 nn(1)=ixc(2,ep)
586 nn(2)=ixc(3,ep)
587 nn(3)=ixc(4,ep)
588 nn(4)=ixc(5,ep)
589C--------------------------
590C BUILD [F], [Q], [NM], [A-S], [AS] AT NODES
591C--------------------------
592C----------------------------------------------------
593C CALCULATION OF [F]
594C----------------------------------------------------
595C---------
596C 2A1
597C---------
598 x21 =mx23-mx13
599 x34 =(corel(1,3)-corel(1,4))*half
600 y21 =my23-my13
601 y34 =(corel(2,3)-corel(2,4))*half
602 z21 = -z1
603 z34 = z1
604 l12 = sqrt(x21*x21+y21*y21+z2)
605 l34 = sqrt(x34*x34+y34*y34+z2)
606C---------
607C 2A2
608C---------
609 x41 =mx34-mx13
610 x32 =(corel(1,3)-corel(1,2))*half
611 y41 =my34-my13
612 y32 =(corel(2,3)-corel(2,2))*half
613 z41 = -z1
614 z32 = z1
615C----------
616C CALCULATION OF [QN] N=1,4
617C----------
618 a_4=area(ep)*fourth
619C---------- N =1----------
620 sl=one/max(l12,em20)
621 vqn(ep,1,1)=x21*sl
622 vqn(ep,2,1)=y21*sl
623 vqn(ep,3,1)=z21*sl
624 sz2=a_4-gama1
625 sz=z2*l24(ep)
626 sl=one/sqrt(sz+sz2*sz2)
627 vqn(ep,7,1)=-z1*y24
628 vqn(ep,8,1)= z1*x24
629 vqn(ep,9,1)= sz2*sl
630C
631 vqn(ep,7,3)=-vqn(ep,7,1)
632 vqn(ep,8,3)=-vqn(ep,8,1)
633 vqn(ep,7,1)= vqn(ep,7,1)*sl
634 vqn(ep,8,1)= vqn(ep,8,1)*sl
635C
636 vqn(ep,4,1)= vqn(ep,8,1)*vqn(ep,3,1)-vqn(ep,9,1)*vqn(ep,2,1)
637 vqn(ep,5,1)= vqn(ep,9,1)*vqn(ep,1,1)-vqn(ep,7,1)*vqn(ep,3,1)
638 vqn(ep,6,1)= vqn(ep,7,1)*vqn(ep,2,1)-vqn(ep,8,1)*vqn(ep,1,1)
639C---------- N =3----------
640 sl=one/max(l34,em20)
641 vqn(ep,1,3)=x34*sl
642 vqn(ep,2,3)=y34*sl
643 vqn(ep,3,3)=z34*sl
644 sz2=a_4+gama1
645 sl=one/sqrt(sz+sz2*sz2)
646 vqn(ep,7,3)= vqn(ep,7,3)*sl
647 vqn(ep,8,3)= vqn(ep,8,3)*sl
648 vqn(ep,9,3)= sz2*sl
649C
650 vqn(ep,4,3)= vqn(ep,8,3)*vqn(ep,3,3)-vqn(ep,9,3)*vqn(ep,2,3)
651 vqn(ep,5,3)= vqn(ep,9,3)*vqn(ep,1,3)-vqn(ep,7,3)*vqn(ep,3,3)
652 vqn(ep,6,3)= vqn(ep,7,3)*vqn(ep,2,3)-vqn(ep,8,3)*vqn(ep,1,3)
653C---------- N =2----------
654 vqn(ep,1,2)=vqn(ep,1,1)
655 vqn(ep,2,2)=vqn(ep,2,1)
656 vqn(ep,3,2)=vqn(ep,3,1)
657 sz2=a_4+gama2
658 sz=z2*l13(ep)
659 sl=one/sqrt(sz+sz2*sz2)
660 vqn(ep,7,2)=-z1*y13
661 vqn(ep,8,2)= z1*x13
662 vqn(ep,9,2)= sz2*sl
663 vqn(ep,7,4)=-vqn(ep,7,2)
664 vqn(ep,8,4)=-vqn(ep,8,2)
665 vqn(ep,7,2)= vqn(ep,7,2)*sl
666 vqn(ep,8,2)= vqn(ep,8,2)*sl
667C
668 vqn(ep,4,2)= vqn(ep,8,2)*vqn(ep,3,2)-vqn(ep,9,2)*vqn(ep,2,2)
669 vqn(ep,5,2)= vqn(ep,9,2)*vqn(ep,1,2)-vqn(ep,7,2)*vqn(ep,3,2)
670 vqn(ep,6,2)= vqn(ep,7,2)*vqn(ep,2,2)-vqn(ep,8,2)*vqn(ep,1,2)
671C---------- N =4----------
672 vqn(ep,1,4)=vqn(ep,1,3)
673 vqn(ep,2,4)=vqn(ep,2,3)
674 vqn(ep,3,4)=vqn(ep,3,3)
675 sz2=a_4-gama2
676 sl=one/sqrt(sz+sz2*sz2)
677 vqn(ep,7,4)= vqn(ep,7,4)*sl
678 vqn(ep,8,4)= vqn(ep,8,4)*sl
679 vqn(ep,9,4)= sz2*sl
680C
681 vqn(ep,4,4)= vqn(ep,8,4)*vqn(ep,3,4)-vqn(ep,9,4)*vqn(ep,2,4)
682 vqn(ep,5,4)= vqn(ep,9,4)*vqn(ep,1,4)-vqn(ep,7,4)*vqn(ep,3,4)
683 vqn(ep,6,4)= vqn(ep,7,4)*vqn(ep,2,4)-vqn(ep,8,4)*vqn(ep,1,4)
684C
685!! IF(ANIM_PLY < 0 ) THEN
686!! I1 = INOD(NN(1))
687!! I2 = INOD(NN(2))
688!! I3 = INOD(NN(3))
689!! I4 = INOD(NN(4))
690C
691!! NORM =SQRT(VQN(EP,7,1)**2 + VQN(EP,8,1)**2 + VQN(EP,9,1)**2)
692!! VN_NOD(1,I1)=VN_NOD(1,I1)+VQN(EP,7,1)/MAX(NORM, EM20)
693!! VN_NOD(2,I1)=VN_NOD(2,I1)+VQN(EP,8,1)/MAX(NORM, EM20)
694!! VN_NOD(3,I1)=VN_NOD(3,I1)+VQN(EP,9,1)/MAX(NORM, EM20)
695C
696!! NORM =SQRT(VQN(EP,7,2)**2 + VQN(EP,8,2)**2 + VQN(EP,9,2)**2)
697!! VN_NOD(1,I2)=VN_NOD(1,I2)+VQN(EP,7,2)/MAX(NORM, EM20)
698!! VN_NOD(2,I2)=VN_NOD(2,I2)+VQN(EP,8,2)/MAX(NORM, EM20)
699!! VN_NOD(3,I2)=VN_NOD(3,I2)+VQN(EP,9,2)/MAX(NORM, EM20)
700C
701!! NORM =SQRT(VQN(EP,7,3)**2 + VQN(EP,8,3)**2 + VQN(EP,9,3)**2)
702!! VN_NOD(1,I3)=VN_NOD(1,I3)+VQN(EP,7,3)/MAX(NORM, EM20)
703!! VN_NOD(2,I3)=VN_NOD(2,I3)+VQN(EP,8,3)/MAX(NORM, EM20)
704!! VN_NOD(3,I3)=VN_NOD(3,I3)+VQN(EP,9,3)/MAX(NORM, EM20)
705c
706!! NORM =SQRT(VQN(EP,7,4)**2 + VQN(EP,8,4)**2 + VQN(EP,9,4)**2)
707!! VN_NOD(1,I4)=VN_NOD(1,I4)+VQN(EP,7,4)/MAX(NORM, EM20)
708!! VN_NOD(2,I4)=VN_NOD(2,I4)+VQN(EP,8,4)/MAX(NORM, EM20)
709!! VN_NOD(3,I4)=VN_NOD(3,I4)+VQN(EP,9,4)/MAX(NORM, EM20)
710C
711!! ENDIF
712C--------------------------------------------------
713C COMPUTE AS N (MIDDLE OF EDGE)
714C--------------------------------------------------
715C J=1 COTE
716 vnrm(ep,1,1)=vqn(ep,7,1)+vqn(ep,7,2)
717 vnrm(ep,2,1)=vqn(ep,8,1)+vqn(ep,8,2)
718 vnrm(ep,3,1)=vqn(ep,9,1)+vqn(ep,9,2)
719 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))
720 c1 = max(em20,c1)
721 vnrm(ep,1,1)=vnrm(ep,1,1)/c1
722 vnrm(ep,2,1)=vnrm(ep,2,1)/c1
723 vnrm(ep,3,1)=vnrm(ep,3,1)/c1
724 vastn(ep,1,1)=zero
725 vastn(ep,2,1)=l12
726 vastn(ep,3,1)=vastn(ep,1,1)
727 vastn(ep,4,1)=vastn(ep,2,1)
728C J=2
729 vnrm(ep,1,2)=vqn(ep,7,4)+vqn(ep,7,3)
730 vnrm(ep,2,2)=vqn(ep,8,4)+vqn(ep,8,3)
731 vnrm(ep,3,2)=vqn(ep,9,4)+vqn(ep,9,3)
732 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))
733 c1 = max(em20,c1)
734 vnrm(ep,1,2)=vnrm(ep,1,2)/c1
735 vnrm(ep,2,2)=vnrm(ep,2,2)/c1
736 vnrm(ep,3,2)=vnrm(ep,3,2)/c1
737 vastn(ep,1,2)=zero
738 vastn(ep,2,2)=l34
739 vastn(ep,3,2)=vastn(ep,1,2)
740 vastn(ep,4,2)=vastn(ep,2,2)
741C J=3
742 vnrm(ep,1,3)=vqn(ep,7,1)+vqn(ep,7,4)
743 vnrm(ep,2,3)=vqn(ep,8,1)+vqn(ep,8,4)
744 vnrm(ep,3,3)=vqn(ep,9,1)+vqn(ep,9,4)
745 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))
746 c1 = max(em20,c1)
747 vnrm(ep,1,3)=vnrm(ep,1,3)/c1
748 vnrm(ep,2,3)=vnrm(ep,2,3)/c1
749 vnrm(ep,3,3)=vnrm(ep,3,3)/c1
750 vastn(ep,1,3)=-(x41*vqn(ep,4,1)+y41*vqn(ep,5,1)+z41*vqn(ep,6,1))
751 vastn(ep,2,3)= x41*vqn(ep,1,1)+y41*vqn(ep,2,1)+z41*vqn(ep,3,1)
752 vastn(ep,3,3)=-(x41*vqn(ep,4,4)+y41*vqn(ep,5,4)+z41*vqn(ep,6,4))
753 vastn(ep,4,3)= x41*vqn(ep,1,4)+y41*vqn(ep,2,4)+z41*vqn(ep,3,4)
754C J=4
755 vnrm(ep,1,4)=vqn(ep,7,2)+vqn(ep,7,3)
756 vnrm(ep,2,4)=vqn(ep,8,2)+vqn(ep,8,3)
757 vnrm(ep,3,4)=vqn(ep,9,2)+vqn(ep,9,3)
758 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))
759 c1 = max(em20,c1)
760 vnrm(ep,1,4)=vnrm(ep,1,4)/c1
761 vnrm(ep,2,4)=vnrm(ep,2,4)/c1
762 vnrm(ep,3,4)=vnrm(ep,3,4)/c1
763 vastn(ep,1,4)=-(x32*vqn(ep,4,2)+y32*vqn(ep,5,2)+z32*vqn(ep,6,2))
764 vastn(ep,2,4)= x32*vqn(ep,1,2)+y32*vqn(ep,2,2)+z32*vqn(ep,3,2)
765 vastn(ep,3,4)=-(x32*vqn(ep,4,3)+y32*vqn(ep,5,3)+z32*vqn(ep,6,3))
766 vastn(ep,4,4)= x32*vqn(ep,1,3)+y32*vqn(ep,2,3)+z32*vqn(ep,3,3)
767C----------
768C CALCULATION OF [QG] NG=1,4
769C----------
770 a_4=a_4/pg
771 jmx13=mx13*pg
772 jmy13=my13*pg
773 jmz13=z1*pg
774 j2myz=jmz13*jmz13
775C---------- NG =1----------
776 sz2=a_4-gama1
777 sz=z2*l24(ep)
778 sl=sqrt(sz+sz2*sz2)
779 jac(ep,1)=sl*pg
780 sl=one/max(sl,em20)
781 vqg(ep,7,1)=-z1*y24
782 vqg(ep,8,1)= z1*x24
783 vqg(ep,9,1)= sz2*sl
784 vqg(ep,7,3)=-vqg(ep,7,1)
785 vqg(ep,8,3)=-vqg(ep,8,1)
786 vqg(ep,7,1)= vqg(ep,7,1)*sl
787 vqg(ep,8,1)= vqg(ep,8,1)*sl
788C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
789 g1x1=mx23-jmx13
790 g1y1=my23-jmy13
791 c1 = sqrt(g1x1*g1x1 + g1y1*g1y1 +j2myz )
792 g2x1=mx34-jmx13
793 g2y1=my34-jmy13
794 c2 = sqrt(g2x1*g2x1 + g2y1*g2y1 +j2myz)
795 vjfi(ep,1,1)=(g2y1*vqg(ep,9,1)+jmz13*vqg(ep,8,1))
796 vjfi(ep,2,1)=(-jmz13*vqg(ep,7,1)-g2x1*vqg(ep,9,1))
797 vjfi(ep,3,1)=(g2x1*vqg(ep,8,1)-g2y1*vqg(ep,7,1))
798 vqg(ep,1,1)= g1x1*c2+vjfi(ep,1,1)*c1
799 vqg(ep,2,1)= g1y1*c2+vjfi(ep,2,1)*c1
800 vqg(ep,3,1)= -jmz13*c2+vjfi(ep,3,1)*c1
801 sl=sl/pg
802 vjfi(ep,1,1)=vjfi(ep,1,1)*sl
803 vjfi(ep,2,1)=vjfi(ep,2,1)*sl
804 vjfi(ep,3,1)=vjfi(ep,3,1)*sl
805C
806 vjfi(ep,4,1)=(-jmz13*vqg(ep,8,1)-g1y1*vqg(ep,9,1))*sl
807 vjfi(ep,5,1)=(g1x1*vqg(ep,9,1)+jmz13*vqg(ep,7,1))*sl
808 vjfi(ep,6,1)=(g1y1*vqg(ep,7,1)-g1x1*vqg(ep,8,1))*sl
809C
810 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))
811 IF ( sl/=zero) sl = one / sl
812 vqg(ep,1,1) = vqg(ep,1,1)*sl
813 vqg(ep,2,1) = vqg(ep,2,1)*sl
814 vqg(ep,3,1) = vqg(ep,3,1)*sl
815C
816 vqg(ep,4,1)= vqg(ep,8,1)*vqg(ep,3,1)-vqg(ep,9,1)*vqg(ep,2,1)
817 vqg(ep,5,1)= vqg(ep,9,1)*vqg(ep,1,1)-vqg(ep,7,1)*vqg(ep,3,1)
818 vqg(ep,6,1)= vqg(ep,7,1)*vqg(ep,2,1)-vqg(ep,8,1)*vqg(ep,1,1)
819C---------- NG =3----------
820 sz2=a_4+gama1
821 sl=sqrt(sz+sz2*sz2)
822 jac(ep,3)=sl*pg
823 sl=one/max(sl,em20)
824 vqg(ep,7,3)= vqg(ep,7,3)*sl
825 vqg(ep,8,3)= vqg(ep,8,3)*sl
826 vqg(ep,9,3)= sz2*sl
827C
828 g1x3=mx23+jmx13
829 g1y3=my23+jmy13
830 j1 = sqrt(g1x3*g1x3 + g1y3*g1y3 +j2myz )
831C--------G2X3=G2X2------
832 g2x2=mx34+jmx13
833 g2y2=my34+jmy13
834 j2 = sqrt(g2x2*g2x2 + g2y2*g2y2 +j2myz)
835 vjfi(ep,1,3)=(g2y2*vqg(ep,9,3)-jmz13*vqg(ep,8,3))
836 vjfi(ep,2,3)=(jmz13*vqg(ep,7,3)-g2x2*vqg(ep,9,3))
837 vjfi(ep,3,3)=(g2x2*vqg(ep,8,3)-g2y2*vqg(ep,7,3))
838 vqg(ep,1,3)= g1x3*j2+vjfi(ep,1,3)*j1
839 vqg(ep,2,3)= g1y3*j2+vjfi(ep,2,3)*j1
840 vqg(ep,3,3)= jmz13*j2+vjfi(ep,3,3)*j1
841 sl=sl/pg
842 vjfi(ep,1,3)=vjfi(ep,1,3)*sl
843 vjfi(ep,2,3)=vjfi(ep,2,3)*sl
844 vjfi(ep,3,3)=vjfi(ep,3,3)*sl
845C
846 vjfi(ep,4,3)=(jmz13*vqg(ep,8,3)-g1y3*vqg(ep,9,3))*sl
847 vjfi(ep,5,3)=(g1x3*vqg(ep,9,3)-jmz13*vqg(ep,7,3))*sl
848 vjfi(ep,6,3)=(g1y3*vqg(ep,7,3)-g1x3*vqg(ep,8,3))*sl
849C
850 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))
851 IF ( sl/=zero) sl = one / sl
852 vqg(ep,1,3) = vqg(ep,1,3)*sl
853 vqg(ep,2,3) = vqg(ep,2,3)*sl
854 vqg(ep,3,3) = vqg(ep,3,3)*sl
855C
856 vqg(ep,4,3)= vqg(ep,8,3)*vqg(ep,3,3)-vqg(ep,9,3)*vqg(ep,2,3)
857 vqg(ep,5,3)= vqg(ep,9,3)*vqg(ep,1,3)-vqg(ep,7,3)*vqg(ep,3,3)
858 vqg(ep,6,3)= vqg(ep,7,3)*vqg(ep,2,3)-vqg(ep,8,3)*vqg(ep,1,3)
859C---------- NG =2----------
860 sz2=a_4+gama2
861 sz=z2*l13(ep)
862 sl=sqrt(sz+sz2*sz2)
863 jac(ep,2)=sl*pg
864 sl=one/max(sl,em20)
865 vqg(ep,7,2)=-z1*y13
866 vqg(ep,8,2)= z1*x13
867 vqg(ep,9,2)= sz2*sl
868 vqg(ep,7,4)=-vqg(ep,7,2)
869 vqg(ep,8,4)=-vqg(ep,8,2)
870 vqg(ep,7,2)= vqg(ep,7,2)*sl
871 vqg(ep,8,2)= vqg(ep,8,2)*sl
872C
873 vjfi(ep,1,2)=(g2y2*vqg(ep,9,2)-jmz13*vqg(ep,8,2))
874 vjfi(ep,2,2)=(jmz13*vqg(ep,7,2)-g2x2*vqg(ep,9,2))
875 vjfi(ep,3,2)=(g2x2*vqg(ep,8,2)-g2y2*vqg(ep,7,2))
876 vqg(ep,1,2)= g1x1*j2+vjfi(ep,1,2)*c1
877 vqg(ep,2,2)= g1y1*j2+vjfi(ep,2,2)*c1
878 vqg(ep,3,2)=-jmz13*j2+vjfi(ep,3,2)*c1
879 sl=sl/pg
880 vjfi(ep,1,2)=vjfi(ep,1,2)*sl
881 vjfi(ep,2,2)=vjfi(ep,2,2)*sl
882 vjfi(ep,3,2)=vjfi(ep,3,2)*sl
883C
884 vjfi(ep,4,2)=(-jmz13*vqg(ep,8,2)-g1y1*vqg(ep,9,2))*sl
885 vjfi(ep,5,2)=(g1x1*vqg(ep,9,2)+jmz13*vqg(ep,7,2))*sl
886 vjfi(ep,6,2)=(g1y1*vqg(ep,7,2)-g1x1*vqg(ep,8,2))*sl
887C
888 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))
889 IF ( sl/=zero) sl = one / sl
890 vqg(ep,1,2) = vqg(ep,1,2)*sl
891 vqg(ep,2,2) = vqg(ep,2,2)*sl
892 vqg(ep,3,2) = vqg(ep,3,2)*sl
893 vqg(ep,4,2)= vqg(ep,8,2)*vqg(ep,3,2)-vqg(ep,9,2)*vqg(ep,2,2)
894 vqg(ep,5,2)= vqg(ep,9,2)*vqg(ep,1,2)-vqg(ep,7,2)*vqg(ep,3,2)
895 vqg(ep,6,2)= vqg(ep,7,2)*vqg(ep,2,2)-vqg(ep,8,2)*vqg(ep,1,2)
896C---------- NG =4----------
897 sz2=a_4-gama2
898 sl=sqrt(sz+sz2*sz2)
899 jac(ep,4)=sl*pg
900 sl=one/max(sl,em20)
901 vqg(ep,7,4)= vqg(ep,7,4)*sl
902 vqg(ep,8,4)= vqg(ep,8,4)*sl
903 vqg(ep,9,4)= sz2*sl
904C
905 vjfi(ep,1,4)=(g2y1*vqg(ep,9,4)+jmz13*vqg(ep,8,4))
906 vjfi(ep,2,4)=(-jmz13*vqg(ep,7,4)-g2x1*vqg(ep,9,4))
907 vjfi(ep,3,4)=(g2x1*vqg(ep,8,4)-g2y1*vqg(ep,7,4))
908 vqg(ep,1,4)= g1x3*c2+vjfi(ep,1,4)*j1
909 vqg(ep,2,4)= g1y3*c2+vjfi(ep,2,4)*j1
910 vqg(ep,3,4)=jmz13*c2+vjfi(ep,3,4)*j1
911 sl=sl/pg
912 vjfi(ep,1,4)=vjfi(ep,1,4)*sl
913 vjfi(ep,2,4)=vjfi(ep,2,4)*sl
914 vjfi(ep,3,4)=vjfi(ep,3,4)*sl
915C
916 vjfi(ep,4,4)=(jmz13*vqg(ep,8,4)-g1y3*vqg(ep,9,4))*sl
917 vjfi(ep,5,4)=(g1x3*vqg(ep,9,4)-jmz13*vqg(ep,7,4))*sl
918 vjfi(ep,6,4)=(g1y3*vqg(ep,7,4)-g1x3*vqg(ep,8,4))*sl
919C
920 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))
921 IF ( sl/=zero) sl = one / sl
922 vqg(ep,1,4) = vqg(ep,1,4)*sl
923 vqg(ep,2,4) = vqg(ep,2,4)*sl
924 vqg(ep,3,4) = vqg(ep,3,4)*sl
925 vqg(ep,4,4)= vqg(ep,8,4)*vqg(ep,3,4)-vqg(ep,9,4)*vqg(ep,2,4)
926 vqg(ep,5,4)= vqg(ep,9,4)*vqg(ep,1,4)-vqg(ep,7,4)*vqg(ep,3,4)
927 vqg(ep,6,4)= vqg(ep,7,4)*vqg(ep,2,4)-vqg(ep,8,4)*vqg(ep,1,4)
928 area(ep)=jac(ep,1)+jac(ep,2)+jac(ep,3)+jac(ep,4)
929 ENDDO
930C ---w/ drilling dof----------
931 IF (isrot>0) THEN
932#include "vectorize.inc"
933 DO i=jft,nplat
934 ep =iplat(i)
935 nn(1)=ixc(2,ep)
936 nn(2)=ixc(3,ep)
937 nn(3)=ixc(4,ep)
938 nn(4)=ixc(5,ep)
939C-------APPLY 3DDL ROTATIONS ----
940 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)))
941 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)))
942 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)))
943 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)))
944 ENDDO
945 END IF !(ISROT>0) THEN
946C
947 DO i=nplat+1,jlt
948 ep =iplat(i)
949 z1=zl1(ep)
950 z2=z1*z1
951 corel(1,1)=vcore(ep,1,1)
952 corel(1,2)=vcore(ep,1,2)
953 corel(1,3)=vcore(ep,1,3)
954 corel(1,4)=vcore(ep,1,4)
955 corel(2,1)=vcore(ep,2,1)
956 corel(2,2)=vcore(ep,2,2)
957 corel(2,3)=vcore(ep,2,3)
958 corel(2,4)=vcore(ep,2,4)
959 nn(1)=ixc(2,ep)
960 nn(2)=ixc(3,ep)
961 nn(3)=ixc(4,ep)
962 nn(4)=ixc(5,ep)
963 x13 =(vcore(ep,1,1)-vcore(ep,1,3))*half
964 x24 =(vcore(ep,1,2)-vcore(ep,1,4))*half
965 y13 =(vcore(ep,2,1)-vcore(ep,2,3))*half
966 y24 =(vcore(ep,2,2)-vcore(ep,2,4))*half
967 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
968 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
969C--------------------------
970C APPLY 6DDL-->5DDL
971C--------------------------
972C-----RXYZ(2,NNOD)=[Q]^T RGXYZ(3,NNOD)-----
973 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))
974 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))
975 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))
976 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))
977 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))
978 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))
979 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))
980 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))
981 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))
982 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))
983 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))
984 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))
985C-----------------------full projection------------------
986 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)
987 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)
988 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)
989C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
990 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)
991 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)
992 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)
993 xz1 = (corel(1,1)-corel(1,2)+corel(1,3)-corel(1,4))*z1
994 yz = (corel(2,1)-corel(2,2)+corel(2,3)-corel(2,4))*z1
995 zz = four*z2
996 d(1)= yy+zz+4
997 d(2)= xx1+zz+4
998 d(3)= xx1+yy+4
999 d(4)= -xy
1000 d(5)= -xz1
1001 d(6)= -yz
1002C
1003 abc = d(1)*d(2)*d(3)
1004 xxyz2 = d(1)*d(6)*d(6)
1005 yyxz2 = d(2)*d(5)*d(5)
1006 zzxy2 = d(3)*d(4)*d(4)
1007 deta = abc+2*d(4)*d(5)*d(6)-xxyz2-yyxz2-zzxy2
1008 IF (deta<em20) THEN
1009 deta=one
1010 ELSE
1011 deta=one/deta
1012 ENDIF
1013 di(ep,1) = (abc-xxyz2)*deta/d(1)
1014 di(ep,2) = (abc-yyxz2)*deta/d(2)
1015 di(ep,3) = (abc-zzxy2)*deta/d(3)
1016 di(ep,4) = (d(5)*d(6)-d(4)*d(3))*deta
1017 di(ep,5) = (d(6)*d(4)-d(5)*d(2))*deta
1018 di(ep,6) = (d(4)*d(5)-d(6)*d(1))*deta
1019C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1020 alr(1) =di(ep,1)*ar(1)+di(ep,4)*ar(2)+di(ep,5)*ar(3)
1021 alr(2) =di(ep,4)*ar(1)+di(ep,2)*ar(2)+di(ep,6)*ar(3)
1022 alr(3) =di(ep,5)*ar(1)+di(ep,6)*ar(2)+di(ep,3)*ar(3)
1023 v13(ep,1)= half*v13(ep,1)+alr(3)*y13
1024 v24(ep,1)= half*v24(ep,1)+alr(3)*y24
1025 vhi(ep,1)= fourth*vhi(ep,1)+(alr(3)*my13-z1*alr(2))
1026 v13(ep,2)= half*v13(ep,2)-alr(3)*x13
1027 v24(ep,2)= half*v24(ep,2)-alr(3)*x24
1028 vhi(ep,2)= fourth*vhi(ep,2)-(alr(3)*mx13-z1*alr(1))
1029 v13(ep,3)= half*v13(ep,3)-(y13*alr(1)-x13*alr(2))
1030 v24(ep,3)= half*v24(ep,3)-(y24*alr(1)-x24*alr(2))
1031 vhi(ep,3)= fourth*vhi(ep,3)+(mx13*alr(2)-my13*alr(1))
1032
1033 vxyz(ep,1,1)= v13(ep,1)+vhi(ep,1)
1034 vxyz(ep,1,2)= v24(ep,1)-vhi(ep,1)
1035 vxyz(ep,1,3)= -v13(ep,1)+vhi(ep,1)
1036 vxyz(ep,1,4)= -v24(ep,1)-vhi(ep,1)
1037
1038 vxyz(ep,2,1)= v13(ep,2)+vhi(ep,2)
1039 vxyz(ep,2,2)= v24(ep,2)-vhi(ep,2)
1040 vxyz(ep,2,3)= -v13(ep,2)+vhi(ep,2)
1041 vxyz(ep,2,4)= -v24(ep,2)-vhi(ep,2)
1042
1043 vxyz(ep,3,1)= v13(ep,3)+vhi(ep,3)
1044 vxyz(ep,3,2)= v24(ep,3)-vhi(ep,3)
1045 vxyz(ep,3,3)= -v13(ep,3)+vhi(ep,3)
1046 vxyz(ep,3,4)= -v24(ep,3)-vhi(ep,3)
1047C
1048 rr(1,1)= rr(1,1)-alr(1)
1049 rr(1,2)= rr(1,2)-alr(1)
1050 rr(1,3)= rr(1,3)-alr(1)
1051 rr(1,4)= rr(1,4)-alr(1)
1052C
1053 rr(2,1)= rr(2,1)-alr(2)
1054 rr(2,2)= rr(2,2)-alr(2)
1055 rr(2,3)= rr(2,3)-alr(2)
1056 rr(2,4)= rr(2,4)-alr(2)
1057C
1058 rr(3,1)= rr(3,1)-alr(3)
1059 rr(3,2)= rr(3,2)-alr(3)
1060 rr(3,3)= rr(3,3)-alr(3)
1061 rr(3,4)= rr(3,4)-alr(3)
1062 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)
1063 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)
1064 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)
1065 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)
1066 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)
1067 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)
1068 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)
1069 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)
1070C--------drilling dof
1071 IF (isrot>0) THEN
1072 rlz(i,1)=(vqn(ep,7,1)*rr(1,1)+vqn(ep,8,1)*rr(2,1)+vqn(ep,9,1)*rr(3,1))
1073 rlz(i,2)=(vqn(ep,7,2)*rr(1,2)+vqn(ep,8,2)*rr(2,2)+vqn(ep,9,2)*rr(3,2))
1074 rlz(i,3)=(vqn(ep,7,3)*rr(1,3)+vqn(ep,8,3)*rr(2,3)+vqn(ep,9,3)*rr(3,3))
1075 rlz(i,4)=(vqn(ep,7,4)*rr(1,4)+vqn(ep,8,4)*rr(2,4)+vqn(ep,9,4)*rr(3,4))
1076 END IF !(ISROT>0) THEN
1077 END DO
1078C
1079C----------------------------
1080C---------Factor for characteristic length
1081C----------------------------
1082 DO i=jft,jlt
1083 rx(i)=xl2(i)+xl3(i)-xl4(i)
1084 ry(i)=yl2(i)+yl3(i)-yl4(i)
1085 sx(i)=-xl2(i)+xl3(i)+xl4(i)
1086 sy(i)=-yl2(i)+yl3(i)+yl4(i)
1087 c1=sqrt(rx(i)*rx(i)+ry(i)*ry(i))
1088 c2=sqrt(sx(i)*sx(i)+sy(i)*sy(i))
1089 s1=fourth*(max(c1,c2)/min(c1,c2)-one)
1090 fac1=min(half,s1)+one
1091 fac2=four*area(i)/(c1*c2)
1092 fac2=3.413*max(zero,fac2-0.7071)
1093 fac2=0.78+0.22*fac2*fac2*fac2
1094 faci=two*fac1*fac2
1095C
1096 facn(i,1)=sqrt(l24(i)/ll(i))
1097 facn(i,2)=sqrt(l13(i)/ll(i))
1098 s1 = sqrt(faci*(four_over_3+lm(i)*area_i(i))*ll(i))
1099 s1 = max(s1,em20)
1100 ll(i) = area(i)/s1
1101 ENDDO
1102
1103 off_l = zero
1104 DO ep=jft,jlt
1105C---------factor for characteristic length---
1106 off(ep) = min(one,abs(offg(ep)))
1107 off_l = min(off_l,offg(ep))
1108 ENDDO
1109C
1110 IF(off_l<zero)THEN
1111 DO ep=jft,jlt
1112 IF(offg(ep)<zero)THEN
1113 vxyz(ep,1,1)= zero
1114 vxyz(ep,2,1)= zero
1115 vxyz(ep,3,1)= zero
1116 vxyz(ep,1,2)= zero
1117 vxyz(ep,2,2)= zero
1118 vxyz(ep,3,2)= zero
1119 vxyz(ep,1,3)= zero
1120 vxyz(ep,2,3)= zero
1121 vxyz(ep,3,3)= zero
1122 vxyz(ep,1,4)= zero
1123 vxyz(ep,2,4)= zero
1124 vxyz(ep,3,4)= zero
1125C
1126 rxyz(ep,1,1)= zero
1127 rxyz(ep,2,1)= zero
1128 rxyz(ep,1,2)= zero
1129 rxyz(ep,2,2)= zero
1130 rxyz(ep,1,3)= zero
1131 rxyz(ep,2,3)= zero
1132 rxyz(ep,1,4)= zero
1133 rxyz(ep,2,4)= zero
1134 ENDIF
1135 ENDDO
1136 ENDIF
1137 IF(anim_ply > 0) THEN
1138#include "lockon.inc"
1139 DO ep=jft,jlt
1140 i1 = inod(ixc(2,ep))
1141 i2 = inod(ixc(3,ep))
1142 i3 = inod(ixc(4,ep))
1143 i4 = inod(ixc(5,ep))
1144cc
1145 vn_nod(1,i1) = vn_nod(1,i1)+vq(ep,1,3)
1146 vn_nod(2,i1) = vn_nod(2,i1)+vq(ep,2,3)
1147 vn_nod(3,i1) = vn_nod(3,i1)+vq(ep,3,3)
1148C
1149 vn_nod(1,i2) =vn_nod(1,i2)+vq(ep,1,3)
1150 vn_nod(2,i2) =vn_nod(2,i2)+vq(ep,2,3)
1151 vn_nod(3,i2) =vn_nod(3,i2)+vq(ep,3,3)
1152C
1153 vn_nod(1,i3) =vn_nod(1,i3)+vq(ep,1,3)
1154 vn_nod(2,i3) =vn_nod(2,i3)+vq(ep,2,3)
1155 vn_nod(3,i3) =vn_nod(3,i3)+vq(ep,3,3)
1156
1157C
1158 vn_nod(1,i4) =vn_nod(1,i4)+vq(ep,1,3)
1159 vn_nod(2,i4) =vn_nod(2,i4)+vq(ep,2,3)
1160 vn_nod(3,i4) =vn_nod(3,i4)+vq(ep,3,3)
1161 ENDDO
1162#include "lockoff.inc"
1163 ENDIF
1164C-------------
1165 RETURN
1166 END
1167!||====================================================================
1168!|| cbacoort ../engine/source/elements/shell/coqueba/cbacoor.F
1169!||--- called by ------------------------------------------------------
1170!|| cbaforc3 ../engine/source/elements/shell/coqueba/cbaforc3.F
1171!||--- calls -----------------------------------------------------
1172!|| clskew3 ../engine/source/elements/sh3n/coquedk/cdkcoor3.F
1173!||--- uses -----------------------------------------------------
1174!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
1175!|| element_mod ../common_source/modules/elements/element_mod.F90
1176!||====================================================================
1177 SUBROUTINE cbacoort(ELBUF_STR,JFT,JLT,X,V,
1178 . VR,DR,IXC,PM,OFFG,AREA,
1179 1 VXYZ, RLZ,VCORE,JAC,HX,
1180 2 HY,VQ,VQG,VJFI,NPLAT,IPLAT,
1181 3 X13_T ,X24_T ,Y13_T,Y24_T,NPT ,
1182 4 SMSTR , ISROT ,XLCOR,ZL ,VQN,NEL)
1183C-----------------------------------------------
1184C M o d u l e s
1185C-----------------------------------------------
1186 USE elbufdef_mod
1187 use element_mod , only : nixc
1188C-----------------------------------------------
1189C I M P L I C I T T Y P E S
1190C-----------------------------------------------
1191#include "implicit_f.inc"
1192#include "comlock.inc"
1193c-----------------------------------------------
1194c g l o b a l p a r a m e t e r s
1195c-----------------------------------------------
1196#include "mvsiz_p.inc"
1197#include "param_c.inc"
1198C-----------------------------------------------
1199C D U M M Y A R G U M E N T S
1200C-----------------------------------------------
1201 INTEGER JFT,JLT,NNOD,NPLAT,NPT,ISROT,NEL
1202 INTEGER IXC(NIXC,*),IPLAT(*)
1203 parameter(nnod = 4)
1204 my_real
1205 . x(3,*), v(3,*), vr(3,*),pm(npropm,*),offg(*),
1206 . vcore(mvsiz,3,nnod),vxyz(mvsiz,3,nnod),area(*),vjfi(mvsiz,6,4),
1207 . vq(mvsiz,3,3),vqg(mvsiz,9,nnod),rlz(mvsiz,nnod),
1208 . jac(mvsiz,4),hx(mvsiz,4),hy(mvsiz,4),y24_t(*),dr(3,*),
1209 . x13_t(*),x24_t(*),y13_t(*),xlcor(mvsiz,2,nnod-1),zl(*),vqn(mvsiz,9,nnod)
1210 double precision
1211 . smstr(*)
1212 TYPE(elbuf_struct_) :: ELBUF_STR
1213C-----RLZ(NNOD,*)
1214C-----------------------------------------------
1215c FUNCTION: preliminary compute for QBAT and Ismstr=10
1216c
1217c Note:
1218c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
1219c
1220c TYPE NAME FUNCTION
1221c I JFT,JLT - element id limit
1222c I X, V, VR(3,NUMNOD) - coordinate,velocity, rotational velocity in global system
1223c I IXC(NIXC,*NEL) - element connectivity of shell (and other data)
1224c I PM(NPROPM,MID) - input Material data
1225c I OFFG(NEL),OFF(NEL) - element activation flag, local flag
1226c O AREA(NEL),NTP - Area, num. of integrating points in thickness
1227c O VCORE(MVSIZ,3,4) - coordinates of 4 nodes in local system
1228c BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
1229c for plat element
1230c O VXYZ(3,4,NEL) - nodal velocity (4 nodes) in local system
1231C V13->VXYZ(,1,),V24->VXYZ(,2,),VHI->VXYZ(,3,)------------
1232c for plat element
1233c O RXYZ(3,4,NEL) - nodal rotational velocity (4 nodes) in local system 5 dof
1234C R13->RXYZ(,1,),R24->RXYZ(,2,),RHI->RXYZ(,3,),RTI->RXYZ(,4,)----
1235c for plat element
1236c O JAC,HX,HY(4,NEL) - Jacobian,asymmetric part(hourglass) at 4 Gauss points
1237c O VQ(9,NEL),VQN,VQG(9,4,) - local system of element, nodal(4) and Gauss points(4)
1238c O Xij_T,Yij_T(NEL) - [B] components used for in-plane shear assumed strain
1239c IO XiS,YiS,ZiS(NEL) - reference configurations used for small strain options
1240c I ISMSTR - small strain
1241c O ZL1(NEL) - warped length
1242c I IDRIL - drilling dof flag
1243c (used only Idril active )
1244C-----------------------------------------------
1245C L O C A L V A R I A B L E S
1246C-----------------------------------------------
1247 INTEGER I,II(9),K,EP,SPLAT
1248 INTEGER NN(4),JPLAT(MVSIZ),IFPROJ
1249 MY_REAL
1250 . XX2,XX3,XX4,YY1,YY2,YY3,YY4,ZZ1,ZZ2,ZZ3,ZZ4
1251 MY_REAL
1252 . C1,C2,L13(MVSIZ),L24(MVSIZ),
1253 . RR(3,NNOD),OFF_L
1254 MY_REAL
1255 . TOL,PG,J0,J1,J2,COREL(3,4),X1,Y1,LL(MVSIZ)
1256 MY_REAL
1257 . X13,X24,Y13,MX13,MX23,MX34,MY13,Z1,Z2,GAMA1,GAMA2,
1258 . X21,X34,Y21,Y34,Z21,Z34,X41,X32,Y41,Y32,Z41,Z32,L12,L34,
1259 . A_4,SL,SZ2,SZ,JMX13,JMY13,JMZ13,J2MYZ,Y24,
1260 . MY23,MY34,G1X1,G1X3,G1Y1,G1Y3,G2X1,G2X2,G2Y1,G2Y2,
1261 . XX1,YY,ZZ,XY,XZ1,YZ,
1262 . XXYZ2,ZZXY2,YYXZ2
1263 my_real
1264 . lxyz0(3),deta1(mvsiz),
1265 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),
1266 . yl3(mvsiz),yl4(mvsiz),xx,area_i(mvsiz),
1267 . x0g2(mvsiz),x0g3(mvsiz),x0g4(mvsiz),y0g2(mvsiz),
1268 . y0g3(mvsiz),y0g4(mvsiz),z0g2(mvsiz),z0g3(mvsiz),z0g4(mvsiz)
1269C
1270 my_real
1271 . vg24(3),vghi(3),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
1272 . vr1_12,vr1_21,
1273 . vhix,
1274 . vnx1, vny1, vnz1,vnx2,vny2,vnz2,vnx3,vny3,vnz3,vnx4,vny4,vnz4,
1275 . norm,zl1(mvsiz),ssz(mvsiz),
1276 . rx(mvsiz),ry(mvsiz),rz(mvsiz),sx(mvsiz),sy(mvsiz),
1277 . r11(mvsiz),r12(mvsiz),r13(mvsiz),r21(mvsiz),r22(mvsiz),
1278 . r23(mvsiz),r31(mvsiz),r32(mvsiz),r33(mvsiz),dirz(mvsiz,2),
1279 . x0l2(mvsiz),x0l3(mvsiz),x0l4(mvsiz),y0l2(mvsiz),
1280 . y0l3(mvsiz),y0l4(mvsiz),vq0(mvsiz,3,3)
1281C-------- [F] will contains only Rz ->2D
1282 my_real
1283 . axyz(mvsiz,3,nnod)
1284 DATA pg/.577350269189626/
1285C=======================================================================
1286 DO i=1,9
1287 ii(i) = nel*(i-1)
1288 ENDDO
1289C
1290 tol=em12
1291 IF (isrot > 0 ) tol=em8
1292 DO i=jft,jlt
1293 IF(abs(offg(i))==one)offg(i)=sign(two,offg(i))
1294 axyz(i,1:3,1:4)= zero
1295C
1296 IF (isrot > 0 ) THEN
1297 nn(1)=ixc(2,i)
1298 nn(2)=ixc(3,i)
1299 nn(3)=ixc(4,i)
1300 nn(4)=ixc(5,i)
1301 axyz(i,1,1) = dr(1,nn(1))
1302 axyz(i,2,1) = dr(2,nn(1))
1303 axyz(i,3,1) = dr(3,nn(1))
1304 axyz(i,1,2) = dr(1,nn(2))
1305 axyz(i,2,2) = dr(2,nn(2))
1306 axyz(i,3,2) = dr(3,nn(2))
1307 axyz(i,1,3) = dr(1,nn(3))
1308 axyz(i,2,3) = dr(2,nn(3))
1309 axyz(i,3,3) = dr(3,nn(3))
1310 axyz(i,1,4) = dr(1,nn(4))
1311 axyz(i,2,4) = dr(2,nn(4))
1312 axyz(i,3,4) = dr(3,nn(4))
1313 END IF !(ISROT > 0 ) THEN
1314
1315 x0g2(i) = smstr(ii(1)+i)
1316 y0g2(i) = smstr(ii(2)+i)
1317 z0g2(i) = smstr(ii(3)+i)
1318 x0g3(i) = smstr(ii(4)+i)
1319 y0g3(i) = smstr(ii(5)+i)
1320 z0g3(i) = smstr(ii(6)+i)
1321 x0g4(i) = smstr(ii(7)+i)
1322 y0g4(i) = smstr(ii(8)+i)
1323 z0g4(i) = smstr(ii(9)+i)
1324 ENDDO
1325C--- in [F](2,2) only Rz is inside, Rx,Ry have been excluded
1326C-bases at T=0 : B0_g : x0g(3,nnod-1) saved in SMSTR, w/ origine at n1
1327C--- B0_l : x0l(2,nnod-1)+Z1,VCORE at only T=0, but can be calculated from x0g(3,nnod-1)
1328C-bases at T>0 : B_G : X(3,numnod), B_l : by rotate B0_l w/ Rz0
1329C---------UL : xl - x0l' (rotated around Rz0)
1330C-- normal in initial conf. to compute Z1
1331 DO i=jft,jlt
1332 rx(i)=x0g2(i)+x0g3(i)-x0g4(i)
1333 sx(i)=x0g3(i)+x0g4(i)-x0g2(i)
1334 ry(i)=y0g2(i)+y0g3(i)-y0g4(i)
1335 sy(i)=y0g3(i)+y0g4(i)-y0g2(i)
1336 rz(i)=z0g2(i)+z0g3(i)-z0g4(i)
1337 ssz(i)=z0g3(i)+z0g4(i)-z0g2(i)
1338 ENDDO
1339C----------------------------
1340C LOCAL SYSTEM VQ0
1341C----------------------------
1342 k = 0
1343 CALL clskew3(jft,jlt,k,
1344 . rx, ry, rz,
1345 . sx, sy, ssz,
1346 . r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
1347 DO i=jft,jlt
1348 area(i)=fourth*deta1(i)
1349 area_i(i)=one/area(i)
1350 vq0(i,1,1)=r11(i)
1351 vq0(i,2,1)=r21(i)
1352 vq0(i,3,1)=r31(i)
1353 vq0(i,1,2)=r12(i)
1354 vq0(i,2,2)=r22(i)
1355 vq0(i,3,2)=r32(i)
1356 vq0(i,1,3)=r13(i)
1357 vq0(i,2,3)=r23(i)
1358 vq0(i,3,3)=r33(i)
1359 ENDDO
1360 DO i=jft,jlt
1361 lxyz0(1)=fourth*(x0g2(i)+x0g3(i)+x0g4(i))
1362 lxyz0(2)=fourth*(y0g2(i)+y0g3(i)+y0g4(i))
1363 lxyz0(3)=fourth*(z0g2(i)+z0g3(i)+z0g4(i))
1364 zl1(i)=-(vq0(i,1,3)*lxyz0(1)+vq0(i,2,3)*lxyz0(2)+vq0(i,3,3)*lxyz0(3))
1365 ENDDO
1366C----------------------------
1367C xl =R1^t x0l; R1^t=(VQ0^t*VQ)^t---extract Rz0 of R1
1368C----------------------------
1369 DO i=jft,jlt
1370 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)
1371 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)
1372 dirz(i,2) = half*(vr1_12-vr1_21)
1373 norm = one-dirz(i,2)*dirz(i,2)
1374 dirz(i,1) = sqrt(max(zero,norm))
1375 ENDDO
1376 DO i=jft,jlt
1377 x0l2(i)=vq0(i,1,1)*x0g2(i)+vq0(i,2,1)*y0g2(i)+vq0(i,3,1)*z0g2(i)
1378 y0l2(i)=vq0(i,1,2)*x0g2(i)+vq0(i,2,2)*y0g2(i)+vq0(i,3,2)*z0g2(i)
1379 x0l3(i)=vq0(i,1,1)*x0g3(i)+vq0(i,2,1)*y0g3(i)+vq0(i,3,1)*z0g3(i)
1380 y0l3(i)=vq0(i,1,2)*x0g3(i)+vq0(i,2,2)*y0g3(i)+vq0(i,3,2)*z0g3(i)
1381 x0l4(i)=vq0(i,1,1)*x0g4(i)+vq0(i,2,1)*y0g4(i)+vq0(i,3,1)*z0g4(i)
1382 y0l4(i)=vq0(i,1,2)*x0g4(i)+vq0(i,2,2)*y0g4(i)+vq0(i,3,2)*z0g4(i)
1383 ENDDO
1384C----------------------------
1385C Rotate x0l of Rz0
1386C----------------------------
1387 DO i=jft,jlt
1388 xl2(i)= x0l2(i)*dirz(i,1)-y0l2(i)*dirz(i,2)
1389 yl2(i)= x0l2(i)*dirz(i,2)+y0l2(i)*dirz(i,1)
1390 xl3(i)= x0l3(i)*dirz(i,1)-y0l3(i)*dirz(i,2)
1391 yl3(i)= x0l3(i)*dirz(i,2)+y0l3(i)*dirz(i,1)
1392 xl4(i)= x0l4(i)*dirz(i,1)-y0l4(i)*dirz(i,2)
1393 yl4(i)= x0l4(i)*dirz(i,2)+y0l4(i)*dirz(i,1)
1394 x0l2(i)=xl2(i)
1395 y0l2(i)=yl2(i)
1396 x0l3(i)=xl3(i)
1397 y0l3(i)=yl3(i)
1398 x0l4(i)=xl4(i)
1399 y0l4(i)=yl4(i)
1400 ENDDO
1401 DO i=jft,jlt
1402 xl2(i)=xlcor(i,1,1)
1403 yl2(i)=xlcor(i,2,1)
1404 xl3(i)=xlcor(i,1,2)
1405 yl3(i)=xlcor(i,2,2)
1406 xl4(i)=xlcor(i,1,3)
1407 yl4(i)=xlcor(i,2,3)
1408 ENDDO
1409C-----------UL : xl - x0l' (rotated of rz0)
1410 DO i=jft,jlt
1411 v13(i,1)=-xl3(i)-(-x0l3(i))
1412 v24(i,1)=xl2(i)-xl4(i)-(x0l2(i)-x0l4(i))
1413 vhi(i,1)=-xl2(i)+xl3(i)-xl4(i)-(-x0l2(i)+x0l3(i)-x0l4(i))
1414 v13(i,2)=-yl3(i)-(-y0l3(i))
1415 v24(i,2)=yl2(i)-yl4(i)-(y0l2(i)-y0l4(i))
1416 vhi(i,2)=-yl2(i)+yl3(i)-yl4(i)-(-y0l2(i)+y0l3(i)-y0l4(i))
1417 v13(i,3)=zero
1418 v24(i,3)=zero
1419 vhi(i,3)=four*(zl(i)-zl1(i))
1420 ENDDO
1421C-------set up B_l by rotating B0_l (Rz0)
1422 DO i=jft,jlt
1423 xl2(i)=x0l2(i)
1424 yl2(i)=y0l2(i)
1425 xl3(i)=x0l3(i)
1426 yl3(i)=y0l3(i)
1427 xl4(i)=x0l4(i)
1428 yl4(i)=y0l4(i)
1429 ENDDO
1430C----------------------------
1431 DO ep=jft,jlt
1432 lxyz0(1)=fourth*(xl2(ep)+xl3(ep)+xl4(ep))
1433 lxyz0(2)=fourth*(yl2(ep)+yl3(ep)+yl4(ep))
1434 vcore(ep,1,1)=-lxyz0(1)
1435 vcore(ep,1,2)=xl2(ep)-lxyz0(1)
1436 vcore(ep,1,3)=xl3(ep)-lxyz0(1)
1437 vcore(ep,1,4)=xl4(ep)-lxyz0(1)
1438 vcore(ep,2,1)=-lxyz0(2)
1439 vcore(ep,2,2)=yl2(ep)-lxyz0(2)
1440 vcore(ep,2,3)=yl3(ep)-lxyz0(2)
1441 vcore(ep,2,4)=yl4(ep)-lxyz0(2)
1442C
1443 x13_t(ep) =(vcore(ep,1,1)-vcore(ep,1,3))*half
1444 x24_t(ep) =(vcore(ep,1,2)-vcore(ep,1,4))*half
1445 y13_t(ep) =(vcore(ep,2,1)-vcore(ep,2,3))*half
1446 y24_t(ep) =(vcore(ep,2,2)-vcore(ep,2,4))*half
1447 l13(ep)=x13_t(ep)*x13_t(ep)+y13_t(ep)*y13_t(ep)
1448 l24(ep)=x24_t(ep)*x24_t(ep)+y24_t(ep)*y24_t(ep)
1449 ll(ep)=max(l13(ep),l24(ep))
1450 ENDDO
1451c
1452 nplat=jft-1
1453 splat= 0
1454 DO ep=jft,jlt
1455 z2=zl1(ep)*zl1(ep)
1456 IF (z2<tol*ll(ep).OR.npt==1) THEN
1457 nplat=nplat+1
1458 iplat(nplat)=ep
1459 ELSE
1460 splat=splat+1
1461 jplat(splat)=ep
1462 ENDIF
1463 ENDDO
1464 DO ep=nplat+1,jlt
1465 iplat(ep)=jplat(ep-nplat)
1466 ENDDO
1467#include "vectorize.inc"
1468 DO i=jft,nplat
1469 ep =iplat(i)
1470 x13 =x13_t(ep)
1471 x24 =x24_t(ep)
1472 y13 =y13_t(ep)
1473 y24 =y24_t(ep)
1474C
1475 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
1476 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
1477 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
1478 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
1479 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
1480 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
1481 x13_t(ep) =x13*area_i(ep)
1482 x24_t(ep) =x24*area_i(ep)
1483 y13_t(ep) =y13*area_i(ep)
1484 y24_t(ep) =y24*area_i(ep)
1485C--------GAMA(2)
1486 gama1=-mx13*y24+my13*x24
1487 gama2= mx13*y13-my13*x13
1488 z2=zl1(ep)*zl1(ep)
1489CC-------V13(I)->VXYZ(I,1),V24(I)->VXYZ(I,2),VHI(I)->VXYZ(I,3)------------
1490 vxyz(ep,1,1)=v13(ep,1)
1491 vxyz(ep,2,1)=v13(ep,2)
1492 vxyz(ep,3,1)=v13(ep,3)
1493C
1494 vxyz(ep,1,2)=v24(ep,1)
1495 vxyz(ep,2,2)=v24(ep,2)
1496 vxyz(ep,3,2)=v24(ep,3)
1497C
1498 vxyz(ep,1,3)=vhi(ep,1)
1499 vxyz(ep,2,3)=vhi(ep,2)
1500 vxyz(ep,3,3)=vhi(ep,3)
1501C-------BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
1502C-------ARE STORED IN VCORE---
1503 vcore(ep,1,1)=y24_t(ep)
1504 vcore(ep,2,1)=-y13_t(ep)
1505 vcore(ep,3,1)=-x24_t(ep)
1506 vcore(ep,1,2)= x13_t(ep)
1507 vcore(ep,2,2)=gama1*area_i(ep)
1508 vcore(ep,3,2)=gama2*area_i(ep)
1509 vcore(ep,1,3)=mx23
1510 vcore(ep,2,3)=my23
1511 vcore(ep,3,3)=mx34
1512 vcore(ep,1,4)=my34
1513 vcore(ep,2,4)=mx13
1514 vcore(ep,3,4)=my13
1515C
1516 j1=(mx23*my13-mx13*my23)*pg
1517 j2=-(mx13*my34-mx34*my13)*pg
1518 j0=area(ep)*fourth
1519C----
1520 jac(ep,1)=abs(j0+j2-j1)
1521 jac(ep,2)=abs(j0+j2+j1)
1522 jac(ep,3)=abs(j0-j2+j1)
1523 jac(ep,4)=abs(j0-j2-j1)
1524 j1=(my23-my34)*pg
1525 j2=-(my23+my34)*pg
1526 hx(ep,1)=j1/jac(ep,1)
1527 hx(ep,2)=j2/jac(ep,2)
1528 hx(ep,3)=-j1/jac(ep,3)
1529 hx(ep,4)=-j2/jac(ep,4)
1530 j1=(mx34-mx23)*pg
1531 j2=(mx34+mx23)*pg
1532 hy(ep,1)=j1/jac(ep,1)
1533 hy(ep,2)=j2/jac(ep,2)
1534 hy(ep,3)=-j1/jac(ep,3)
1535 hy(ep,4)=-j2/jac(ep,4)
1536 ENDDO
1537C ---w/ drilling dof------drilling may not work -> check RLZj w/ Rz0
1538 IF (isrot>0) THEN
1539#include "vectorize.inc"
1540 DO i=jft,nplat
1541 ep =iplat(i)
1542C-------TRANSMET 3DDL ROTATIONS ----
1543 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)
1544 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)
1545 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)
1546 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)
1547 ENDDO
1548 END IF !(ISROT>0) THEN
1549C----------group warped elements , VJFI,VQN,VQG
1550 ifproj = 0
1551#include "vectorize.inc"
1552 DO i=nplat+1,jlt
1553 ep =iplat(i)
1554 z1=zl1(ep)
1555 z2=z1*z1
1556 vcore(ep,3,1)=z1
1557 vcore(ep,3,2)=-z1
1558 vcore(ep,3,3)=z1
1559 vcore(ep,3,4)=-z1
1560 corel(1,1)=vcore(ep,1,1)
1561 corel(1,2)=vcore(ep,1,2)
1562 corel(1,3)=vcore(ep,1,3)
1563 corel(1,4)=vcore(ep,1,4)
1564 corel(2,1)=vcore(ep,2,1)
1565 corel(2,2)=vcore(ep,2,2)
1566 corel(2,3)=vcore(ep,2,3)
1567 corel(2,4)=vcore(ep,2,4)
1568 x13 =x13_t(ep)
1569 x24 =x24_t(ep)
1570 y13 =y13_t(ep)
1571 y24 =y24_t(ep)
1572cc
1573 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
1574 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
1575 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
1576 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
1577 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
1578 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
1579 x13_t(ep) =x13*area_i(ep)
1580 x24_t(ep) =x24*area_i(ep)
1581 y13_t(ep) =y13*area_i(ep)
1582 y24_t(ep) =y24*area_i(ep)
1583C--------GAMA(2)
1584 gama1=-mx13*y24+my13*x24
1585 gama2= mx13*y13-my13*x13
1586C--------------------------
1587C CONSTRUIRE [F], [Q], [NM], [A-S], [AS] AUX NOEUDS
1588C--------------------------
1589 x21 =mx23-mx13
1590 x34 =(corel(1,3)-corel(1,4))*half
1591 y21 =my23-my13
1592 y34 =(corel(2,3)-corel(2,4))*half
1593 z21 = -z1
1594 z34 = z1
1595 l12 = sqrt(x21*x21+y21*y21+z2)
1596 l34 = sqrt(x34*x34+y34*y34+z2)
1597C---------
1598 x41 =mx34-mx13
1599 x32 =(corel(1,3)-corel(1,2))*half
1600 y41 =my34-my13
1601 y32 =(corel(2,3)-corel(2,2))*half
1602 z41 = -z1
1603 z32 = z1
1604C----------
1605 a_4=area(ep)*fourth
1606C
1607C----------
1608C CALCUL DE [QG] NG=1,4
1609C----------
1610 a_4=a_4/pg
1611 jmx13=mx13*pg
1612 jmy13=my13*pg
1613 jmz13=z1*pg
1614 j2myz=jmz13*jmz13
1615C---------- NG =1----------
1616 sz2=a_4-gama1
1617 sz=z2*l24(ep)
1618 sl=sqrt(sz+sz2*sz2)
1619 jac(ep,1)=sl*pg
1620 sl=one/max(sl,em20)
1621 vqg(ep,7,1)=-z1*y24
1622 vqg(ep,8,1)= z1*x24
1623 vqg(ep,9,1)= sz2*sl
1624 vqg(ep,7,3)=-vqg(ep,7,1)
1625 vqg(ep,8,3)=-vqg(ep,8,1)
1626 vqg(ep,7,1)= vqg(ep,7,1)*sl
1627 vqg(ep,8,1)= vqg(ep,8,1)*sl
1628C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1629 g1x1=mx23-jmx13
1630 g1y1=my23-jmy13
1631 c1 = sqrt(g1x1*g1x1 + g1y1*g1y1 +j2myz )
1632 g2x1=mx34-jmx13
1633 g2y1=my34-jmy13
1634 c2 = sqrt(g2x1*g2x1 + g2y1*g2y1 +j2myz)
1635 vjfi(ep,1,1)=(g2y1*vqg(ep,9,1)+jmz13*vqg(ep,8,1))
1636 vjfi(ep,2,1)=(-jmz13*vqg(ep,7,1)-g2x1*vqg(ep,9,1))
1637 vjfi(ep,3,1)=(g2x1*vqg(ep,8,1)-g2y1*vqg(ep,7,1))
1638 vqg(ep,1,1)= g1x1*c2+vjfi(ep,1,1)*c1
1639 vqg(ep,2,1)= g1y1*c2+vjfi(ep,2,1)*c1
1640 vqg(ep,3,1)= -jmz13*c2+vjfi(ep,3,1)*c1
1641 sl=sl/pg
1642 vjfi(ep,1,1)=vjfi(ep,1,1)*sl
1643 vjfi(ep,2,1)=vjfi(ep,2,1)*sl
1644 vjfi(ep,3,1)=vjfi(ep,3,1)*sl
1645C
1646 vjfi(ep,4,1)=(-jmz13*vqg(ep,8,1)-g1y1*vqg(ep,9,1))*sl
1647 vjfi(ep,5,1)=(g1x1*vqg(ep,9,1)+jmz13*vqg(ep,7,1))*sl
1648 vjfi(ep,6,1)=(g1y1*vqg(ep,7,1)-g1x1*vqg(ep,8,1))*sl
1649C
1650 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))
1651 IF ( sl/=zero) sl = one / sl
1652 vqg(ep,1,1) = vqg(ep,1,1)*sl
1653 vqg(ep,2,1) = vqg(ep,2,1)*sl
1654 vqg(ep,3,1) = vqg(ep,3,1)*sl
1655C
1656 vqg(ep,4,1)= vqg(ep,8,1)*vqg(ep,3,1)-vqg(ep,9,1)*vqg(ep,2,1)
1657 vqg(ep,5,1)= vqg(ep,9,1)*vqg(ep,1,1)-vqg(ep,7,1)*vqg(ep,3,1)
1658 vqg(ep,6,1)= vqg(ep,7,1)*vqg(ep,2,1)-vqg(ep,8,1)*vqg(ep,1,1)
1659C---------- NG =3----------
1660 sz2=a_4+gama1
1661 sl=sqrt(sz+sz2*sz2)
1662 jac(ep,3)=sl*pg
1663 sl=one/max(sl,em20)
1664 vqg(ep,7,3)= vqg(ep,7,3)*sl
1665 vqg(ep,8,3)= vqg(ep,8,3)*sl
1666 vqg(ep,9,3)= sz2*sl
1667C
1668 g1x3=mx23+jmx13
1669 g1y3=my23+jmy13
1670 j1 = sqrt(g1x3*g1x3 + g1y3*g1y3 +j2myz )
1671C--------G2X3=G2X2------
1672 g2x2=mx34+jmx13
1673 g2y2=my34+jmy13
1674 j2 = sqrt(g2x2*g2x2 + g2y2*g2y2 +j2myz)
1675 vjfi(ep,1,3)=(g2y2*vqg(ep,9,3)-jmz13*vqg(ep,8,3))
1676 vjfi(ep,2,3)=(jmz13*vqg(ep,7,3)-g2x2*vqg(ep,9,3))
1677 vjfi(ep,3,3)=(g2x2*vqg(ep,8,3)-g2y2*vqg(ep,7,3))
1678 vqg(ep,1,3)= g1x3*j2+vjfi(ep,1,3)*j1
1679 vqg(ep,2,3)= g1y3*j2+vjfi(ep,2,3)*j1
1680 vqg(ep,3,3)= jmz13*j2+vjfi(ep,3,3)*j1
1681 sl=sl/pg
1682 vjfi(ep,1,3)=vjfi(ep,1,3)*sl
1683 vjfi(ep,2,3)=vjfi(ep,2,3)*sl
1684 vjfi(ep,3,3)=vjfi(ep,3,3)*sl
1685C
1686 vjfi(ep,4,3)=(jmz13*vqg(ep,8,3)-g1y3*vqg(ep,9,3))*sl
1687 vjfi(ep,5,3)=(g1x3*vqg(ep,9,3)-jmz13*vqg(ep,7,3))*sl
1688 vjfi(ep,6,3)=(g1y3*vqg(ep,7,3)-g1x3*vqg(ep,8,3))*sl
1689C
1690 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))
1691 IF ( sl/=zero) sl = one / sl
1692 vqg(ep,1,3) = vqg(ep,1,3)*sl
1693 vqg(ep,2,3) = vqg(ep,2,3)*sl
1694 vqg(ep,3,3) = vqg(ep,3,3)*sl
1695C
1696 vqg(ep,4,3)= vqg(ep,8,3)*vqg(ep,3,3)-vqg(ep,9,3)*vqg(ep,2,3)
1697 vqg(ep,5,3)= vqg(ep,9,3)*vqg(ep,1,3)-vqg(ep,7,3)*vqg(ep,3,3)
1698 vqg(ep,6,3)= vqg(ep,7,3)*vqg(ep,2,3)-vqg(ep,8,3)*vqg(ep,1,3)
1699C---------- NG =2----------
1700 sz2=a_4+gama2
1701 sz=z2*l13(ep)
1702 sl=sqrt(sz+sz2*sz2)
1703 jac(ep,2)=sl*pg
1704 sl=one/max(sl,em20)
1705 vqg(ep,7,2)=-z1*y13
1706 vqg(ep,8,2)= z1*x13
1707 vqg(ep,9,2)= sz2*sl
1708 vqg(ep,7,4)=-vqg(ep,7,2)
1709 vqg(ep,8,4)=-vqg(ep,8,2)
1710 vqg(ep,7,2)= vqg(ep,7,2)*sl
1711 vqg(ep,8,2)= vqg(ep,8,2)*sl
1712C
1713 vjfi(ep,1,2)=(g2y2*vqg(ep,9,2)-jmz13*vqg(ep,8,2))
1714 vjfi(ep,2,2)=(jmz13*vqg(ep,7,2)-g2x2*vqg(ep,9,2))
1715 vjfi(ep,3,2)=(g2x2*vqg(ep,8,2)-g2y2*vqg(ep,7,2))
1716 vqg(ep,1,2)= g1x1*j2+vjfi(ep,1,2)*c1
1717 vqg(ep,2,2)= g1y1*j2+vjfi(ep,2,2)*c1
1718 vqg(ep,3,2)=-jmz13*j2+vjfi(ep,3,2)*c1
1719 sl=sl/pg
1720 vjfi(ep,1,2)=vjfi(ep,1,2)*sl
1721 vjfi(ep,2,2)=vjfi(ep,2,2)*sl
1722 vjfi(ep,3,2)=vjfi(ep,3,2)*sl
1723C
1724 vjfi(ep,4,2)=(-jmz13*vqg(ep,8,2)-g1y1*vqg(ep,9,2))*sl
1725 vjfi(ep,5,2)=(g1x1*vqg(ep,9,2)+jmz13*vqg(ep,7,2))*sl
1726 vjfi(ep,6,2)=(g1y1*vqg(ep,7,2)-g1x1*vqg(ep,8,2))*sl
1727C
1728 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))
1729 IF ( sl/=zero) sl = one / sl
1730 vqg(ep,1,2) = vqg(ep,1,2)*sl
1731 vqg(ep,2,2) = vqg(ep,2,2)*sl
1732 vqg(ep,3,2) = vqg(ep,3,2)*sl
1733 vqg(ep,4,2)= vqg(ep,8,2)*vqg(ep,3,2)-vqg(ep,9,2)*vqg(ep,2,2)
1734 vqg(ep,5,2)= vqg(ep,9,2)*vqg(ep,1,2)-vqg(ep,7,2)*vqg(ep,3,2)
1735 vqg(ep,6,2)= vqg(ep,7,2)*vqg(ep,2,2)-vqg(ep,8,2)*vqg(ep,1,2)
1736C---------- NG =4----------
1737 sz2=a_4-gama2
1738 sl=sqrt(sz+sz2*sz2)
1739 jac(ep,4)=sl*pg
1740 sl=one/max(sl,em20)
1741 vqg(ep,7,4)= vqg(ep,7,4)*sl
1742 vqg(ep,8,4)= vqg(ep,8,4)*sl
1743 vqg(ep,9,4)= sz2*sl
1744C
1745 vjfi(ep,1,4)=(g2y1*vqg(ep,9,4)+jmz13*vqg(ep,8,4))
1746 vjfi(ep,2,4)=(-jmz13*vqg(ep,7,4)-g2x1*vqg(ep,9,4))
1747 vjfi(ep,3,4)=(g2x1*vqg(ep,8,4)-g2y1*vqg(ep,7,4))
1748 vqg(ep,1,4)= g1x3*c2+vjfi(ep,1,4)*j1
1749 vqg(ep,2,4)= g1y3*c2+vjfi(ep,2,4)*j1
1750 vqg(ep,3,4)=jmz13*c2+vjfi(ep,3,4)*j1
1751 sl=sl/pg
1752 vjfi(ep,1,4)=vjfi(ep,1,4)*sl
1753 vjfi(ep,2,4)=vjfi(ep,2,4)*sl
1754 vjfi(ep,3,4)=vjfi(ep,3,4)*sl
1755C
1756 vjfi(ep,4,4)=(jmz13*vqg(ep,8,4)-g1y3*vqg(ep,9,4))*sl
1757 vjfi(ep,5,4)=(g1x3*vqg(ep,9,4)-jmz13*vqg(ep,7,4))*sl
1758 vjfi(ep,6,4)=(g1y3*vqg(ep,7,4)-g1x3*vqg(ep,8,4))*sl
1759C
1760 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))
1761 IF ( sl/=zero) sl = one / sl
1762 vqg(ep,1,4) = vqg(ep,1,4)*sl
1763 vqg(ep,2,4) = vqg(ep,2,4)*sl
1764 vqg(ep,3,4) = vqg(ep,3,4)*sl
1765 vqg(ep,4,4)= vqg(ep,8,4)*vqg(ep,3,4)-vqg(ep,9,4)*vqg(ep,2,4)
1766 vqg(ep,5,4)= vqg(ep,9,4)*vqg(ep,1,4)-vqg(ep,7,4)*vqg(ep,3,4)
1767 vqg(ep,6,4)= vqg(ep,7,4)*vqg(ep,2,4)-vqg(ep,8,4)*vqg(ep,1,4)
1768 area(ep)=jac(ep,1)+jac(ep,2)+jac(ep,3)+jac(ep,4)
1769 ENDDO
1770C--------------------------
1771C TRANSMET 6DDL-->5DDL, w/o projection
1772C--------------------------
1773 DO i=nplat+1,jlt
1774 ep =iplat(i)
1775 vxyz(ep,1,1)= half*v13(ep,1)+fourth*vhi(ep,1)
1776 vxyz(ep,1,2)= half*v24(ep,1)-fourth*vhi(ep,1)
1777 vxyz(ep,1,3)= -half*v13(ep,1)+fourth*vhi(ep,1)
1778 vxyz(ep,1,4)= -half*v24(ep,1)-fourth*vhi(ep,1)
1779
1780 vxyz(ep,2,1)= half*v13(ep,2)+fourth*vhi(ep,2)
1781 vxyz(ep,2,2)= half*v24(ep,2)-fourth*vhi(ep,2)
1782 vxyz(ep,2,3)= -half*v13(ep,2)+fourth*vhi(ep,2)
1783 vxyz(ep,2,4)= -half*v24(ep,2)-fourth*vhi(ep,2)
1784C
1785 vxyz(ep,3,1)= half*v13(ep,3)+fourth*vhi(ep,3)
1786 vxyz(ep,3,2)= half*v24(ep,3)-fourth*vhi(ep,3)
1787 vxyz(ep,3,3)= -half*v13(ep,3)+fourth*vhi(ep,3)
1788 vxyz(ep,3,4)= -half*v24(ep,3)-fourth*vhi(ep,3)
1789C--------drilling dof
1790 IF (isrot>0) THEN
1791 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)
1792 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)
1793 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)
1794 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)
1795 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)
1796 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)
1797 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)
1798 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)
1799 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)
1800 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)
1801 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)
1802 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)
1803C----------not right takes the RLZ with cbacoor
1804 rlz(ep,1)=(vqn(ep,7,1)*rr(1,1)+vqn(ep,8,1)*rr(2,1)+vqn(ep,9,1)*rr(3,1))
1805 rlz(ep,2)=(vqn(ep,7,2)*rr(1,2)+vqn(ep,8,2)*rr(2,2)+vqn(ep,9,2)*rr(3,2))
1806 rlz(ep,3)=(vqn(ep,7,3)*rr(1,3)+vqn(ep,8,3)*rr(2,3)+vqn(ep,9,3)*rr(3,3))
1807 rlz(ep,4)=(vqn(ep,7,4)*rr(1,4)+vqn(ep,8,4)*rr(2,4)+vqn(ep,9,4)*rr(3,4))
1808 END IF !(ISROT>0) THEN
1809 END DO
1810C
1811 off_l = zero
1812 DO ep=jft,jlt
1813 off_l = min(off_l,offg(ep))
1814 ENDDO
1815C
1816 IF(off_l<zero)THEN
1817 DO ep=jft,jlt
1818 IF(offg(ep)<zero)THEN
1819 vxyz(ep,1:3,1:4)= zero
1820 rlz(ep,1:4)= zero
1821 ENDIF
1822 ENDDO
1823 ENDIF
1824C-------------
1825 RETURN
1826 END
cbacoor
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)
Definition cbacoor.F:48
cbacoort
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)
Definition cbacoor.F:1183
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