OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
czcorc.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!|| czcorc1 ../engine/source/elements/shell/coquez/czcorc.F
25!||--- called by ------------------------------------------------------
26!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
27!|| czforc3_crk ../engine/source/elements/xfem/czforc3_crk.F
28!||--- calls -----------------------------------------------------
29!|| clskew3 ../engine/source/elements/sh3n/coquedk/cdkcoor3.F
30!|| cortdir3 ../engine/source/elements/shell/coque/cortdir3.F
31!|| czcorp5 ../engine/source/elements/shell/coquez/czcorp5.F
32!|| czcorp5x ../engine/source/elements/shell/coquez/czcorc.F
33!||--- uses -----------------------------------------------------
34!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
35!||====================================================================
36 SUBROUTINE czcorc1(NUMNOD ,NUMELC ,ELBUF_STR,
37 1 JFT ,JLT ,X ,V ,VR ,
38 2 IXC ,PM ,PLAT ,AREA ,
39 3 AREA_I ,V13 ,V24 ,VHI ,RLXYZ ,
40 4 VQN ,VQ ,LL ,L13 ,L24 ,
41 5 X13 ,X24 ,Y13 ,Y24 ,MX13 ,
42 6 MX23 ,MX34 ,MY13 ,MY23 ,MY34 ,
43 7 Z1 ,COREL ,DI ,DB ,SMSTR ,
44 9 IREP ,NPT ,NLAY ,ISMSTR ,
45 A DIR_A ,DIR_B ,OFFG ,RLXYZV ,CORELV ,
46 B FACN ,PY1 ,PX2 ,PY2 ,R11 ,
47 C R12 ,R13 ,R21 ,R22 ,R23 ,
48 D R31 ,R32 ,R33 ,RLZ ,IDRIL ,
49 E IXFEM ,VX1 ,VX2 ,VX3 ,VX4 ,
50 F VY1 ,VY2 ,VY3 ,VY4 ,VZ1 ,
51 G VZ2 ,VZ3 ,VZ4 ,VRX1 ,VRX2 ,
52 H VRX3 ,VRX4 ,VRY1 ,VRY2 ,VRY3 ,
53 I VRY4 ,VRZ1 ,VRZ2 ,VRZ3 ,VRZ4 ,
54 J X1G ,X2G ,X3G ,X4G ,Y1G ,
55 K Y2G ,Y3G ,Y4G ,Z1G ,Z2G ,
56 L Z3G ,Z4G ,THK ,DIZ ,UX1 ,
57 M UX2 ,UX3 ,UX4 ,UY1 ,UY2 ,
58 N UY3 ,UY4 ,XL2 ,XL3 ,XL4 ,
59 O YL2 ,YL3 ,YL4 ,VL1 ,VL2 ,
60 P VL3 ,VL4 ,NEL ,Z2 )
61C-----------------------------------------------
62C M o d u l e s
63C-----------------------------------------------
64 USE elbufdef_mod
65C-----------------------------------------------
66#include "implicit_f.inc"
67c-----------------------------------------------
68c g l o b a l p a r a m e t e r s
69c-----------------------------------------------
70#include "mvsiz_p.inc"
71#include "param_c.inc"
72#include "parit_c.inc"
73#include "scr05_c.inc"
74#include "scr17_c.inc"
75#include "com08_c.inc"
76#include "impl1_c.inc"
77C-----------------------------------------------
78C D U M M Y A R G U M E N T S
79C-----------------------------------------------
80 LOGICAL PLAT(*)
81 INTEGER, INTENT(IN) :: NUMNOD,NUMELC
82 INTEGER IXC(NIXC,*),JFT,JLT,IREP,NPT,NLAY,ISMSTR,IDRIL,IXFEM,NEL
83 my_real
84 . PM(NPROPM,*), X(3,*), V(3,*), VR(3,*),RLXYZ(MVSIZ,2,4),
85 . V13(MVSIZ,3),V24(MVSIZ,3),VHI(MVSIZ,3),MX23(NEL),MY13(NEL),MY23(NEL),MY34(NEL),
86 . LL(NEL),L13(NEL),L24(NEL),X13(NEL),X24(NEL),Y13(NEL),Y24(NEL),MX13(NEL),
87 . VQ(MVSIZ,3,3),AREA(NEL),Z1(NEL),MX34(NEL),VQN(MVSIZ,3,4),AREA_I(NEL),
88 . corel(mvsiz,2,4),di(mvsiz,6),db(mvsiz,3,4),dir_a(nel,*),dir_b(nel,*),offg(nel),
89 . rlxyzv(mvsiz,2,4),corelv(mvsiz,2,4),facn(mvsiz,2),px2(*),
90 . py1(nel), py2(nel),r11(mvsiz),r12(mvsiz),r13(mvsiz),
91 . r21(mvsiz),r22(mvsiz),r23(mvsiz),r31(mvsiz),r32(mvsiz),
92 . r33(mvsiz),rlz(mvsiz,4),diz(mvsiz,3),thk(*),
93 . vx1(mvsiz),vx2(mvsiz),vx3(mvsiz),vx4(mvsiz),
94 . vy1(mvsiz),vy2(mvsiz),vy3(mvsiz),vy4(mvsiz),
95 . vz1(mvsiz),vz2(mvsiz),vz3(mvsiz),vz4(mvsiz),
96 . vrx1(mvsiz),vrx2(mvsiz),vrx3(mvsiz),vrx4(mvsiz),
97 . vry1(mvsiz),vry2(mvsiz),vry3(mvsiz),vry4(mvsiz),
98 . vrz1(mvsiz),vrz2(mvsiz),vrz3(mvsiz),vrz4(mvsiz),
99 . x1g(mvsiz),x2g(mvsiz),x3g(mvsiz),x4g(mvsiz),
100 . y1g(mvsiz),y2g(mvsiz),y3g(mvsiz),y4g(mvsiz),
101 . z1g(mvsiz),z2g(mvsiz),z3g(mvsiz),z4g(mvsiz),
102 . ux1(mvsiz),ux2(mvsiz),ux3(mvsiz),ux4(mvsiz),
103 . uy1(mvsiz),uy2(mvsiz),uy3(mvsiz),uy4(mvsiz),
104 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),
105 . yl2(mvsiz),yl3(mvsiz),yl4(mvsiz),
106 . z2(mvsiz),vl1(mvsiz,3),vl2(mvsiz,3),vl3(mvsiz,3),vl4(mvsiz,3)
107 TYPE (ELBUF_STRUCT_) :: ELBUF_STR
108 DOUBLE PRECISION
109 . SMSTR(*)
110C-----------------------------------------------
111C L O C A L V A R I A B L E S
112C-----------------------------------------------
113 INTEGER IXCTMP2,IXCTMP3,IXCTMP4,IXCTMP5
114 INTEGER NNOD,I,J,K,L,II(6)
115 PARAMETER (NNOD = 4)
116 my_real
117 . lxyz0(3),deta,deta1(mvsiz),rx(mvsiz), ry(mvsiz), rz(mvsiz),
118 . sx(mvsiz),sy(mvsiz),ssz(mvsiz),
119 . vcore(3,nnod),rrxyz(3,nnod),vg13(3),vg24(3),vghi(3),
120 . tol,x13_2(mvsiz),y13_2(mvsiz),x24_2(mvsiz),y24_2(mvsiz),
121 . a_4,sz,sz1,sz2,sl,c1,c2,s1,
122 . vlxyz(3,nnod),ar(3),ad(nnod),btb(6),xx,yy,zz,xy,xz,yz,
123 . abc,xxyz2,yyxz2,zzxy2,d(6),diz1(6),
124 . alr(3),ald(nnod),dbad(3),btb_c,
125 . hs(mvsiz),faci,fac1,fac2,lm(mvsiz),facdt,
126 . dt05,dt025,exz,eyz,ddrx,ddry,v13x,v24x,vhix,ddrz1,ddrz2
127 my_real :: off_loc
128C----------------------------------------------
129
130 IF(iresp == 1)THEN
131 tol=em7
132 ELSE
133 tol=em8
134 END IF
135C
136 DO i=1,6
137 ii(i) = nel*(i-1)
138 ENDDO
139C
140 IF(ixfem==0)THEN
141 DO i=jft,jlt
142 ixctmp2=ixc(2,i)
143 ixctmp3=ixc(3,i)
144 ixctmp4=ixc(4,i)
145 ixctmp5=ixc(5,i)
146
147 rx(i)=x(1,ixctmp3)+x(1,ixctmp4)-x(1,ixctmp2)-x(1,ixctmp5)
148 sx(i)=x(1,ixctmp4)+x(1,ixctmp5)-x(1,ixctmp2)-x(1,ixctmp3)
149 ry(i)=x(2,ixctmp3)+x(2,ixctmp4)-x(2,ixctmp2)-x(2,ixctmp5)
150 sy(i)=x(2,ixctmp4)+x(2,ixctmp5)-x(2,ixctmp2)-x(2,ixctmp3)
151 rz(i)=x(3,ixctmp3)+x(3,ixctmp4)-x(3,ixctmp2)-x(3,ixctmp5)
152 ssz(i)=x(3,ixctmp4)+x(3,ixctmp5)-x(3,ixctmp2)-x(3,ixctmp3)
153 ENDDO
154 ELSE IF(ixfem==1)THEN
155 DO i=jft,jlt
156 rx(i)=x2g(i)+x3g(i)-x1g(i)-x4g(i)
157 sx(i)=x3g(i)+x4g(i)-x1g(i)-x2g(i)
158 ry(i)=y2g(i)+y3g(i)-y1g(i)-y4g(i)
159 sy(i)=y3g(i)+y4g(i)-y1g(i)-y2g(i)
160 rz(i)=z2g(i)+z3g(i)-z1g(i)-z4g(i)
161 ssz(i)=z3g(i)+z4g(i)-z1g(i)-z2g(i)
162 ENDDO
163 END IF
164C----------------------------
165C LOCAL SYSTEM
166C----------------------------
167 k = 0
168 CALL clskew3(jft,jlt,k,
169 . rx, ry, rz,
170 . sx, sy, ssz,
171 . r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
172 DO i=jft,jlt
173 area(i)=fourth*deta1(i)
174 off_loc = zero
175 IF(abs(offg(i))/=zero) off_loc = one
176 area_i(i)=off_loc/area(i)
177 area_i(i) = max(area_i(i),em20)
178 vq(i,1,1)=r11(i)
179 vq(i,2,1)=r21(i)
180 vq(i,3,1)=r31(i)
181 vq(i,1,2)=r12(i)
182 vq(i,2,2)=r22(i)
183 vq(i,3,2)=r32(i)
184 vq(i,1,3)=r13(i)
185 vq(i,2,3)=r23(i)
186 vq(i,3,3)=r33(i)
187 ENDDO
188C--------------------------
189C TRANSFERT GLOBAL-->LOCAL
190C--------------------------
191 IF(ixfem==0)THEN
192 DO i=jft,jlt
193 ixctmp2=ixc(2,i)
194 ixctmp3=ixc(3,i)
195 ixctmp4=ixc(4,i)
196 ixctmp5=ixc(5,i)
197
198 lxyz0(1)=fourth*(x(1,ixctmp4)+x(1,ixctmp5) + x(1,ixctmp2)+x(1,ixctmp3))
199 lxyz0(2)=fourth*(x(2,ixctmp4)+x(2,ixctmp5) + x(2,ixctmp2)+x(2,ixctmp3))
200 lxyz0(3)=fourth*(x(3,ixctmp4)+x(3,ixctmp5) + x(3,ixctmp2)+x(3,ixctmp3))
201
202 j=ixctmp2
203 k=ixctmp3
204 xx=x(1,k)-x(1,j)
205 yy=x(2,k)-x(2,j)
206 zz=x(3,k)-x(3,j)
207 xl2(i)=r11(i)*xx+r21(i)*yy+r31(i)*zz
208 yl2(i)=r12(i)*xx+r22(i)*yy+r32(i)*zz
209 xx=x(1,j)-lxyz0(1)
210 yy=x(2,j)-lxyz0(2)
211 zz=x(3,j)-lxyz0(3)
212 z1(i)=r13(i)*xx+r23(i)*yy+r33(i)*zz
213 k=ixctmp4
214 xx=x(1,k)-x(1,j)
215 yy=x(2,k)-x(2,j)
216 zz=x(3,k)-x(3,j)
217 xl3(i)=r11(i)*xx+r21(i)*yy+r31(i)*zz
218 yl3(i)=r12(i)*xx+r22(i)*yy+r32(i)*zz
219 k=ixctmp5
220 xx=x(1,k)-x(1,j)
221 yy=x(2,k)-x(2,j)
222 zz=x(3,k)-x(3,j)
223 xl4(i)=r11(i)*xx+r21(i)*yy+r31(i)*zz
224 yl4(i)=r12(i)*xx+r22(i)*yy+r32(i)*zz
225 z2(i) = z1(i)*z1(i)
226 ENDDO
227 ELSE IF(ixfem==1)THEN
228 DO i=jft,jlt
229 lxyz0(1)=fourth*(x3g(i)+x4g(i)+x1g(i)+x2g(i))
230 lxyz0(2)=fourth*(y3g(i)+y4g(i)+y1g(i)+y2g(i))
231 lxyz0(3)=fourth*(z3g(i)+z4g(i)+z1g(i)+z2g(i))
232 xx=x2g(i)-x1g(i)
233 yy=y2g(i)-y1g(i)
234 zz=z2g(i)-z1g(i)
235 xl2(i)=r11(i)*xx+r21(i)*yy+r31(i)*zz
236 yl2(i)=r12(i)*xx+r22(i)*yy+r32(i)*zz
237 xx=x1g(i)-lxyz0(1)
238 yy=y1g(i)-lxyz0(2)
239 zz=z1g(i)-lxyz0(3)
240 z1(i)=r13(i)*xx+r23(i)*yy+r33(i)*zz
241 xx=x3g(i)-x1g(i)
242 yy=y3g(i)-y1g(i)
243 zz=z3g(i)-z1g(i)
244 xl3(i)=r11(i)*xx+r21(i)*yy+r31(i)*zz
245 yl3(i)=r12(i)*xx+r22(i)*yy+r32(i)*zz
246 xx=x4g(i)-x1g(i)
247 yy=y4g(i)-y1g(i)
248 zz=z4g(i)-z1g(i)
249 xl4(i)=r11(i)*xx+r21(i)*yy+r31(i)*zz
250 yl4(i)=r12(i)*xx+r22(i)*yy+r32(i)*zz
251 z2(i) = z1(i)*z1(i)
252 ENDDO
253 END IF
254C----------------------------
255C SMALL STRAIN
256C----------------------------
257 IF(ismstr == 11)THEN
258#include "vectorize.inc"
259 DO i=jft,jlt
260 IF(abs(offg(i)) == one)offg(i)=sign(two,offg(i))
261 ux1(i) = zero
262 uy1(i) = zero
263 ux2(i) = zero
264 uy2(i) = zero
265 ux3(i) = zero
266 uy3(i) = zero
267 ux4(i) = zero
268 uy4(i) = zero
269 IF(abs(offg(i)) == two)THEN
270 ux2(i) = xl2(i)-smstr(ii(1)+i)
271 uy2(i) = yl2(i)-smstr(ii(2)+i)
272 ux3(i) = xl3(i)-smstr(ii(3)+i)
273 uy3(i) = yl3(i)-smstr(ii(4)+i)
274 ux4(i) = xl4(i)-smstr(ii(5)+i)
275 uy4(i) = yl4(i)-smstr(ii(6)+i)
276 xl2(i) = smstr(ii(1)+i)
277 yl2(i) = smstr(ii(2)+i)
278 xl3(i) = smstr(ii(3)+i)
279 yl3(i) = smstr(ii(4)+i)
280 xl4(i) = smstr(ii(5)+i)
281 yl4(i) = smstr(ii(6)+i)
282 z1(i) = zero
283 area(i) = half*((xl2(i)-xl4(i))*yl3(i)-xl3(i)*(yl2(i)-yl4(i)))
284 area_i(i) = one/max(em20,area(i))
285 ELSE
286 smstr(ii(1)+i) = xl2(i)
287 smstr(ii(2)+i) = yl2(i)
288 smstr(ii(3)+i) = xl3(i)
289 smstr(ii(4)+i) = yl3(i)
290 smstr(ii(5)+i) = xl4(i)
291 smstr(ii(6)+i) = yl4(i)
292 ENDIF
293 ENDDO
294 ELSEIF(ismstr == 1.OR.ismstr == 2)THEN
295#include "vectorize.inc"
296 DO i=jft,jlt
297 IF(abs(offg(i)) == two)THEN
298 xl2(i) = smstr(ii(1)+i)
299 yl2(i) = smstr(ii(2)+i)
300 xl3(i) = smstr(ii(3)+i)
301 yl3(i) = smstr(ii(4)+i)
302 xl4(i) = smstr(ii(5)+i)
303 yl4(i) = smstr(ii(6)+i)
304 z1(i) = zero
305 area(i) = half*((xl2(i)-xl4(i))*yl3(i)-xl3(i)*(yl2(i)-yl4(i)))
306 area_i(i) = one/max(em20,area(i))
307 ELSE
308 smstr(ii(1)+i)=xl2(i)
309 smstr(ii(2)+i)=yl2(i)
310 smstr(ii(3)+i)=xl3(i)
311 smstr(ii(4)+i)=yl3(i)
312 smstr(ii(5)+i)=xl4(i)
313 smstr(ii(6)+i)=yl4(i)
314 ENDIF
315 ENDDO
316 ENDIF
317 IF(ismstr == 1)THEN
318 DO i=jft,jlt
319 IF(offg(i) == one)offg(i)=two
320 ENDDO
321 ENDIF
322C----------------------------
323C ORTHOTROPY
324C----------------------------
325 CALL cortdir3(elbuf_str,dir_a ,dir_b ,jft ,jlt ,
326 . nlay ,irep ,rx ,ry ,rz ,
327 . sx ,sy ,ssz ,r11 ,r21 ,
328 . r31 ,r12 ,r22 ,r32 ,nel )
329C----------------------------
330 IF (ivector == 1)THEN
331 ELSE
332C
333 DO i=jft,jlt
334 lxyz0(1)=fourth*(xl2(i)+xl3(i)+xl4(i))
335 lxyz0(2)=fourth*(yl2(i)+yl3(i)+yl4(i))
336 corel(i,1,1)=-lxyz0(1)
337 corel(i,1,2)=xl2(i)-lxyz0(1)
338 corel(i,1,3)=xl3(i)-lxyz0(1)
339 corel(i,1,4)=xl4(i)-lxyz0(1)
340 corel(i,2,1)=-lxyz0(2)
341 corel(i,2,2)=yl2(i)-lxyz0(2)
342 corel(i,2,3)=yl3(i)-lxyz0(2)
343 corel(i,2,4)=yl4(i)-lxyz0(2)
344 x13(i)=(corel(i,1,1)-corel(i,1,3))*half
345 x24(i)=(corel(i,1,2)-corel(i,1,4))*half
346 y13(i)=(corel(i,2,1)-corel(i,2,3))*half
347 y24(i)=(corel(i,2,2)-corel(i,2,4))*half
348C
349 mx13(i)=(corel(i,1,1)+corel(i,1,3))*half
350 mx23(i)=(corel(i,1,2)+corel(i,1,3))*half
351 mx34(i)=(corel(i,1,3)+corel(i,1,4))*half
352 my13(i)=(corel(i,2,1)+corel(i,2,3))*half
353 my23(i)=(corel(i,2,2)+corel(i,2,3))*half
354 my34(i)=(corel(i,2,3)+corel(i,2,4))*half
355C
356 py1(i) = -x24(i)
357 px2(i) = -y13(i)
358 py2(i) = x13(i)
359C-
360 x13_2(i) =x13(i)*x13(i)
361 y13_2(i) =y13(i)*y13(i)
362 x24_2(i) =x24(i)*x24(i)
363 y24_2(i) =y24(i)*y24(i)
364 l13(i)=x13_2(i)+y13_2(i)
365C-------------------------------
366 l24(i)=x24_2(i)+y24_2(i)
367C ---------Factor for characteristic length---
368 c1 =corel(i,1,2)*corel(i,2,4)-corel(i,2,2)*corel(i,1,4)
369 c2 =corel(i,1,1)*corel(i,2,3)-corel(i,2,1)*corel(i,1,3)
370 hs(i) =max(abs(c1),abs(c2))*area_i(i)
371C
372 ENDDO
373 ENDIF
374C----------------------------
375C---------characteristic length ratio ---
376C----------------------------
377 facdt = five_over_4
378C-------same than QBAT
379 IF (idt1sh>0) facdt =four_over_3
380 DO i=jft,jlt
381 rx(i)=xl2(i)+xl3(i)-xl4(i)
382 ry(i)=yl2(i)+yl3(i)-yl4(i)
383 sx(i)=-xl2(i)+xl3(i)+xl4(i)
384 sy(i)=-yl2(i)+yl3(i)+yl4(i)
385 c1=sqrt(rx(i)*rx(i)+ry(i)*ry(i))
386 c2=sqrt(sx(i)*sx(i)+sy(i)*sy(i))
387 s1=fourth*(max(c1,c2)/min(c1,c2)-one)
388 fac1=min(half,s1)+one
389 fac2=four*area(i)/(c1*c2)
390 fac2=3.413*max(zero,fac2-0.7071)
391 fac2=0.78+0.22*fac2*fac2*fac2
392 faci=two*fac1*fac2
393C
394 ll(i)=max(l13(i),l24(i))
395 lm(i)=half*(l13(i)+l24(i))
396 facn(i,1)=sqrt(l24(i)/ll(i))
397 facn(i,2)=sqrt(l13(i)/ll(i))
398 s1 = sqrt(faci*(facdt+hs(i))*ll(i))
399 s1 = max(s1,em10)
400 ll(i) = area(i)/s1
401 ENDDO
402C-------cas thk >>L----
403 IF (impl_s>0) THEN
404 DO i=jft,jlt
405 s1=em01*thk(i)*thk(i)
406 lm(i)=max(lm(i),s1)
407 ENDDO
408 END IF
409 IF (ivector == 1)THEN
410 ELSE
411 IF(ixfem==0)THEN
412 DO i=jft,jlt
413 k=ixc(2,i)
414 rlxyz(i,1,1) =vq(i,1,1)*vr(1,k)+vq(i,2,1)*vr(2,k)+vq(i,3,1)*vr(3,k)
415 rlxyz(i,2,1) =vq(i,1,2)*vr(1,k)+vq(i,2,2)*vr(2,k)+vq(i,3,2)*vr(3,k)
416 k=ixc(3,i)
417 rlxyz(i,1,2) =vq(i,1,1)*vr(1,k)+vq(i,2,1)*vr(2,k)+vq(i,3,1)*vr(3,k)
418 rlxyz(i,2,2) =vq(i,1,2)*vr(1,k)+vq(i,2,2)*vr(2,k)+vq(i,3,2)*vr(3,k)
419 k=ixc(4,i)
420 rlxyz(i,1,3) =vq(i,1,1)*vr(1,k)+vq(i,2,1)*vr(2,k)+vq(i,3,1)*vr(3,k)
421 rlxyz(i,2,3) =vq(i,1,2)*vr(1,k)+vq(i,2,2)*vr(2,k)+vq(i,3,2)*vr(3,k)
422 k=ixc(5,i)
423 rlxyz(i,1,4) =vq(i,1,1)*vr(1,k)+vq(i,2,1)*vr(2,k)+vq(i,3,1)*vr(3,k)
424 rlxyz(i,2,4) =vq(i,1,2)*vr(1,k)+vq(i,2,2)*vr(2,k)+vq(i,3,2)*vr(3,k)
425 ENDDO
426 ELSE IF(ixfem==1)THEN
427 DO i=jft,jlt
428
429 rlxyz(i,1,1) =vq(i,1,1)*vrx1(i)+vq(i,2,1)*vry1(i)+vq(i,3,1)*vrz1(i)
430 rlxyz(i,2,1) =vq(i,1,2)*vrx1(i)+vq(i,2,2)*vry1(i)+vq(i,3,2)*vrz1(i)
431
432 rlxyz(i,1,2) =vq(i,1,1)*vrx2(i)+vq(i,2,1)*vry2(i)+vq(i,3,1)*vrz2(i)
433 rlxyz(i,2,2) =vq(i,1,2)*vrx2(i)+vq(i,2,2)*vry2(i)+vq(i,3,2)*vrz2(i)
434
435 rlxyz(i,1,3) =vq(i,1,1)*vrx3(i)+vq(i,2,1)*vry3(i)+vq(i,3,1)*vrz3(i)
436 rlxyz(i,2,3) =vq(i,1,2)*vrx3(i)+vq(i,2,2)*vry3(i)+vq(i,3,2)*vrz3(i)
437
438 rlxyz(i,1,4) =vq(i,1,1)*vrx4(i)+vq(i,2,1)*vry4(i)+vq(i,3,1)*vrz4(i)
439 rlxyz(i,2,4) =vq(i,1,2)*vrx4(i)+vq(i,2,2)*vry4(i)+vq(i,3,2)*vrz4(i)
440 ENDDO
441 ENDIF
442 ENDIF
443C
444 IF(ixfem==0)THEN
445 DO i=jft,jlt
446
447 ixctmp2=ixc(2,i)
448 ixctmp3=ixc(3,i)
449 ixctmp4=ixc(4,i)
450 ixctmp5=ixc(5,i)
451
452 vl1(i,1)=v(1,ixctmp2)
453 vl1(i,2)=v(2,ixctmp2)
454 vl1(i,3)=v(3,ixctmp2)
455 vl2(i,1)=v(1,ixctmp3)
456 vl2(i,2)=v(2,ixctmp3)
457 vl2(i,3)=v(3,ixctmp3)
458 vl3(i,1)=v(1,ixctmp4)
459 vl3(i,2)=v(2,ixctmp4)
460 vl3(i,3)=v(3,ixctmp4)
461 vl4(i,1)=v(1,ixctmp5)
462 vl4(i,2)=v(2,ixctmp5)
463 vl4(i,3)=v(3,ixctmp5)
464 ENDDO
465 DO i=jft,jlt
466 vg13(1)=vl1(i,1)-vl3(i,1)
467 vg24(1)=vl2(i,1)-vl4(i,1)
468 vghi(1)=vl1(i,1)-vl2(i,1)+vl3(i,1)-vl4(i,1)
469C
470 vg13(2)=vl1(i,2)-vl3(i,2)
471 vg24(2)=vl2(i,2)-vl4(i,2)
472 vghi(2)=vl1(i,2)-vl2(i,2)+vl3(i,2)-vl4(i,2)
473C
474 vg13(3)=vl1(i,3)-vl3(i,3)
475 vg24(3)=vl2(i,3)-vl4(i,3)
476 vghi(3)=vl1(i,3)-vl2(i,3)+vl3(i,3)-vl4(i,3)
477C
478 v13(i,1)=(vq(i,1,1)*vg13(1)+vq(i,2,1)*vg13(2)+vq(i,3,1)*vg13(3))
479 v24(i,1)=(vq(i,1,1)*vg24(1)+vq(i,2,1)*vg24(2)+vq(i,3,1)*vg24(3))
480 vhi(i,1)=(vq(i,1,1)*vghi(1)+vq(i,2,1)*vghi(2)+vq(i,3,1)*vghi(3))
481 v13(i,2)=(vq(i,1,2)*vg13(1)+vq(i,2,2)*vg13(2)+vq(i,3,2)*vg13(3))
482 v24(i,2)=(vq(i,1,2)*vg24(1)+vq(i,2,2)*vg24(2)+vq(i,3,2)*vg24(3))
483 vhi(i,2)=(vq(i,1,2)*vghi(1)+vq(i,2,2)*vghi(2)+vq(i,3,2)*vghi(3))
484 v13(i,3)=(vq(i,1,3)*vg13(1)+vq(i,2,3)*vg13(2)+vq(i,3,3)*vg13(3))
485 v24(i,3)=(vq(i,1,3)*vg24(1)+vq(i,2,3)*vg24(2)+vq(i,3,3)*vg24(3))
486 vhi(i,3)=(vq(i,1,3)*vghi(1)+vq(i,2,3)*vghi(2)+vq(i,3,3)*vghi(3))
487 ENDDO
488 ELSE IF(ixfem==1)THEN
489 DO i=jft,jlt
490 vg13(1)=vx1(i)-vx3(i)
491 vg24(1)=vx2(i)-vx4(i)
492 vghi(1)=vx1(i)-vx2(i)+vx3(i)-vx4(i)
493C
494 vg13(2)=vy1(i)-vy3(i)
495 vg24(2)=vy2(i)-vy4(i)
496 vghi(2)=vy1(i)-vy2(i)+vy3(i)-vy4(i)
497C
498 vg13(3)=vz1(i)-vz3(i)
499 vg24(3)=vz2(i)-vz4(i)
500 vghi(3)=vz1(i)-vz2(i)+vz3(i)-vz4(i)
501C
502 v13(i,1)=(vq(i,1,1)*vg13(1)+vq(i,2,1)*vg13(2)+vq(i,3,1)*vg13(3))
503 v24(i,1)=(vq(i,1,1)*vg24(1)+vq(i,2,1)*vg24(2)+vq(i,3,1)*vg24(3))
504 vhi(i,1)=(vq(i,1,1)*vghi(1)+vq(i,2,1)*vghi(2)+vq(i,3,1)*vghi(3))
505 v13(i,2)=(vq(i,1,2)*vg13(1)+vq(i,2,2)*vg13(2)+vq(i,3,2)*vg13(3))
506 v24(i,2)=(vq(i,1,2)*vg24(1)+vq(i,2,2)*vg24(2)+vq(i,3,2)*vg24(3))
507 vhi(i,2)=(vq(i,1,2)*vghi(1)+vq(i,2,2)*vghi(2)+vq(i,3,2)*vghi(3))
508 v13(i,3)=(vq(i,1,3)*vg13(1)+vq(i,2,3)*vg13(2)+vq(i,3,3)*vg13(3))
509 v24(i,3)=(vq(i,1,3)*vg24(1)+vq(i,2,3)*vg24(2)+vq(i,3,3)*vg24(3))
510 vhi(i,3)=(vq(i,1,3)*vghi(1)+vq(i,2,3)*vghi(2)+vq(i,3,3)*vghi(3))
511 ENDDO
512 END IF
513C--------------------------
514 IF (idril>0) THEN
515 IF(ixfem==0)THEN
516 DO i=jft,jlt
517 k=ixc(2,i)
518 rlz(i,1) =(vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k))*area_i(i)
519 k=ixc(3,i)
520 rlz(i,2) =(vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k))*area_i(i)
521 k=ixc(4,i)
522 rlz(i,3) =(vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k))*area_i(i)
523 k=ixc(5,i)
524 rlz(i,4) =(vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k))*area_i(i)
525 ENDDO
526 ELSE IF(ixfem==1)THEN
527 DO i=jft,jlt
528 rlz(i,1) =(vq(i,1,3)*vrx1(i)+vq(i,2,3)*vry1(i)+vq(i,3,3)*vrz1(i))*area_i(i)
529 rlz(i,2) =(vq(i,1,3)*vrx2(i)+vq(i,2,3)*vry2(i)+vq(i,3,3)*vrz2(i))*area_i(i)
530 rlz(i,3) =(vq(i,1,3)*vrx3(i)+vq(i,2,3)*vry3(i)+vq(i,3,3)*vrz3(i))*area_i(i)
531 rlz(i,4) =(vq(i,1,3)*vrx4(i)+vq(i,2,3)*vry4(i)+vq(i,3,3)*vrz4(i))*area_i(i)
532 ENDDO
533 ENDIF
534 END IF
535C--------------------------
536C-------Correction 2nd order rigid rotation due to V(t+dt/2),X(t+dt)----
537C--------------------------
538 IF (impl_s == 0.OR.imp_lr>0) THEN
539 dt05 = half*dt1
540 dt025 =fourth*dt1
541 DO i=jft,jlt
542 exz = y24(i)*v13(i,3)-y13(i)*v24(i,3)
543 eyz = -x24(i)*v13(i,3)+x13(i)*v24(i,3)
544 ddry=dt05*exz*area_i(i)
545 ddrx=dt05*eyz*area_i(i)
546 v13x = v13(i,1)
547 v24x = v24(i,1)
548 vhix = vhi(i,1)
549 IF (abs(x13(i)-x24(i)) < em10) THEN
550 ddrz1 = zero
551 ELSE
552 ddrz1 = dt025*(v13(i,2)-v24(i,2))/(x13(i)-x24(i))
553 ENDIF
554 v13(i,1) = v13(i,1)-ddry*v13(i,3)-ddrz1*v13(i,2)
555 v24(i,1) = v24(i,1)-ddry*v24(i,3)-ddrz1*v24(i,2)
556 vhi(i,1) = vhi(i,1)-ddry*vhi(i,3)-ddrz1*vhi(i,2)
557 IF (abs(y13(i)+y24(i)) < em10) THEN
558 ddrz2 = zero
559 ELSE
560 ddrz2 = dt025*(v13x+v24x)/(y13(i)+y24(i))
561 ENDIF
562 v13(i,2) = v13(i,2)-ddrx*v13(i,3)-ddrz2*v13x
563 v24(i,2) = v24(i,2)-ddrx*v24(i,3)-ddrz2*v24x
564 vhi(i,2) = vhi(i,2)-ddrx*vhi(i,3)-ddrz2*vhix
565 ENDDO
566 ENDIF
567C--------------------------
568C------- 3 ROTATIONS to 2 ROTATIONS----
569C--------------------------
570 IF (ivector == 1)THEN
571 ELSE
572 IF(ixfem==0)THEN
573 CALL czcorp5(numnod ,jlt ,numelc ,vr ,npt ,tol ,
574 2 ixc ,plat ,area ,area_i ,v13 ,
575 3 v24 ,vhi ,rlxyz ,vqn ,vq ,
576 4 x13 ,x24 ,y13 ,y24 ,mx13 ,
577 6 my13 ,
578 7 z1 ,di ,db ,corel ,rlz ,
579 8 lm ,
580 9 l13 ,l24 ,idril ,diz )
581
582 ELSE IF(ixfem==1)THEN
583 CALL czcorp5x(jft ,jlt ,vr ,npt ,tol ,
584 2 ixc ,plat ,area ,area_i ,v13 ,
585 3 v24 ,vhi ,rlxyz ,vqn ,vq ,
586 4 x13 ,x24 ,y13 ,y24 ,mx13 ,
587 6 mx23 ,mx34 ,my13 ,my23 ,my34 ,
588 7 z1 ,di ,db ,corel ,rlz ,
589 8 lm ,x13_2 ,y13_2 ,x24_2 ,y24_2 ,
590 9 l13 ,l24 ,vrx1 ,vrx2 ,vrx3 ,
591 a vrx4 ,vry1 ,vry2 ,vry3 ,vry4 ,
592 b vrz1 ,vrz2 ,vrz3 ,vrz4 )
593 END IF
594 END IF !IF (IVECTOR == 1)THEN
595C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
596 DO i=jft,jlt
597 v13(i,1)=v13(i,1)*area_i(i)
598 v24(i,1)=v24(i,1)*area_i(i)
599 vhi(i,1)=vhi(i,1)*fourth
600C
601 v13(i,2)=v13(i,2)*area_i(i)
602 v24(i,2)=v24(i,2)*area_i(i)
603 vhi(i,2)=vhi(i,2)*fourth
604C
605 v13(i,3)=v13(i,3)*area_i(i)
606 v24(i,3)=v24(i,3)*area_i(i)
607 vhi(i,3)=vhi(i,3)*fourth
608 ENDDO
609C
610 RETURN
611 END
612!||====================================================================
613!|| a3invdp ../engine/source/elements/shell/coquez/czcorc.f
614!||--- called by ------------------------------------------------------
615!|| czcorp5v ../engine/source/elements/shell/coquez/czcorc.F
616!|| czcorp5x ../engine/source/elements/shell/coquez/czcorc.F
617!||====================================================================
618 SUBROUTINE a3invdp(D,DI)
619C-----------------------------------------------
620C I m p l i c i t T y p e s
621C-----------------------------------------------
622#include "implicit_f.inc"
623C-----------------------------------------------
624C D u m m y A r g u m e n t s
625C-----------------------------------------------
626 my_real d(6), di(6)
627C-----------------------------------------------
628C L o c a l V a r i a b l e s
629C-----------------------------------------------
630 double precision
631 . abc,xxyz2,yyxz2,zzxy2,deta
632C-----------------------------------------------
633 abc = d(1)*d(2)*d(3)
634 xxyz2 = d(1)*d(6)*d(6)
635 yyxz2 = d(2)*d(5)*d(5)
636 zzxy2 = d(3)*d(4)*d(4)
637 deta = abs(abc+two*d(4)*d(5)*d(6)-xxyz2-yyxz2-zzxy2)
638 deta = one/max(deta,em20)
639 di(1) = (abc-xxyz2)*deta/max(d(1),em20)
640 di(2) = (abc-yyxz2)*deta/max(d(2),em20)
641 di(3) = (abc-zzxy2)*deta/max(d(3),em20)
642 di(4) = (d(5)*d(6)-d(4)*d(3))*deta
643 di(5) = (d(6)*d(4)-d(5)*d(2))*deta
644 di(6) = (d(4)*d(5)-d(6)*d(1))*deta
645C
646 RETURN
647 END
648!||====================================================================
649!|| czcorp4v ../engine/source/elements/shell/coquez/czcorc.F
650!||====================================================================
651 SUBROUTINE czcorp4v(JFT ,JLT ,VR ,NPT ,TOL ,
652 2 IXC ,PLAT ,AREA ,AREA_I ,V13 ,
653 3 V24 ,VHI ,RLXYZV ,VQN ,VQ ,
654 4 X13 ,X24 ,Y13 ,Y24 ,MX13 ,
655 6 MX23 ,MX34 ,MY13 ,MY23 ,MY34 ,
656 7 Z1 ,DI ,DB ,CORELV ,RLZ ,
657 8 LL ,X13_2 ,Y13_2 ,X24_2 ,Y24_2 ,
658 9 L13 ,L24 )
659C-----------------------------------------------
660C I m p l i c i t T y p e s
661C-----------------------------------------------
662#include "implicit_f.inc"
663C-----------------------------------------------
664C G l o b a l P a r a m e t e r s
665C-----------------------------------------------
666#include "mvsiz_p.inc"
667C-----------------------------------------------
668C D u m m y A r g u m e n t s
669C-----------------------------------------------
670 LOGICAL PLAT(*)
671 INTEGER IXC(NIXC,*),JFT,JLT,NPT
672 my_real
673 . vr(3,*),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
674 . mx23(*),my13(*),my23(*),my34(*),
675 . x13(*),x24(*),y13(*),y24(*),mx13(*),
676 . vq(mvsiz,3,3),area(*),z1(*),mx34(*),vqn(mvsiz,3,4),area_i(*),
677 . di(mvsiz,6),db(mvsiz,3,4),ll(*),l13(*),l24(*),
678 . tol,x13_2(*),y13_2(*),x24_2(*),y24_2(*),
679 . rlxyzv(mvsiz,2,4),corelv(mvsiz,2,4),rlz(mvsiz,4)
680C-----------------------------------------------
681C L O C A L V A R I A B L E S
682C-----------------------------------------------
683 INTEGER NNOD,I,J,K,L
684 PARAMETER (NNOD = 4)
685 my_real
686 . deta,
687 . vcore(3,nnod),rrxyz(3,nnod),vg13(3),vg24(3),vghi(3),
688 . z2(mvsiz),a_4,sz,sz1,sz2,sl,c1,c2,s1,
689 . ar(3),ad(nnod),btb(6),xx,yy,zz,xy,xz,yz,
690 . abc,xxyz2,yyxz2,zzxy2,d(6),
691 . alr(3),ald(nnod),dbad(3),btb_c
692C-----------------------------------------------
693 DO i=jft,jlt
694 z2(i)=z1(i)*z1(i)
695 IF (z2(i)<ll(i)*tol.OR.npt == 1) THEN
696 z1(i)=zero
697 plat(i)=.true.
698 ELSE
699 plat(i)=.false.
700C--------------------------------------------------
701C WARPING SPECIAL TREATMENT
702C full projection eliminer drilling rotations and rigid rotations
703C--------------------------------------------------------------------------
704 a_4=area(i)*fourth
705C---------- node N ----------
706 sz1=mx13(i)*y24(i)-my13(i)*x24(i)
707 sz2=a_4+sz1
708 sz=z2(i)*l24(i)
709 sl=one/sqrt(sz+sz2*sz2)
710 vqn(i,1,1)=-z1(i)*y24(i)
711 vqn(i,2,1)= z1(i)*x24(i)
712 vqn(i,3,1)= sz2*sl
713 vqn(i,1,3)=-vqn(i,1,1)
714 vqn(i,2,3)=-vqn(i,2,1)
715 vqn(i,1,1)= vqn(i,1,1)*sl
716 vqn(i,2,1)= vqn(i,2,1)*sl
717C
718 sz2=a_4-sz1
719 sl=one/sqrt(sz+sz2*sz2)
720 vqn(i,1,3)= vqn(i,1,3)*sl
721 vqn(i,2,3)= vqn(i,2,3)*sl
722 vqn(i,3,3)= sz2*sl
723C
724 sz1=mx13(i)*y13(i)-my13(i)*x13(i)
725 sz2=a_4+sz1
726 sz=z2(i)*l13(i)
727 sl=one/sqrt(sz+sz2*sz2)
728 vqn(i,1,2)=-z1(i)*y13(i)
729 vqn(i,2,2)= z1(i)*x13(i)
730 vqn(i,3,2)= sz2*sl
731 vqn(i,1,4)=-vqn(i,1,2)
732 vqn(i,2,4)=-vqn(i,2,2)
733 vqn(i,1,2)= vqn(i,1,2)*sl
734 vqn(i,2,2)= vqn(i,2,2)*sl
735C
736 sz2=a_4-sz1
737 sl=one/sqrt(sz+sz2*sz2)
738 vqn(i,1,4)= vqn(i,1,4)*sl
739 vqn(i,2,4)= vqn(i,2,4)*sl
740 vqn(i,3,4)= sz2*sl
741C
742 k=ixc(2,i)
743 rrxyz(1,1) =rlxyzv(i,1,1)
744 rrxyz(2,1) =rlxyzv(i,2,1)
745 rrxyz(3,1) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k)
746
747 k=ixc(3,i)
748 rrxyz(1,2) =rlxyzv(i,1,2)
749 rrxyz(2,2) =rlxyzv(i,2,2)
750 rrxyz(3,2) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k)
751
752 k=ixc(4,i)
753 rrxyz(1,3) =rlxyzv(i,1,3)
754 rrxyz(2,3) =rlxyzv(i,2,3)
755 rrxyz(3,3) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k)
756
757 k=ixc(5,i)
758 rrxyz(1,4) =rlxyzv(i,1,4)
759 rrxyz(2,4) =rlxyzv(i,2,4)
760 rrxyz(3,4) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k)
761C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
762C-----------------------full projection------------------
763 ar(1)=-z1(i)*vhi(i,2)+y13(i)*v13(i,3)+y24(i)*v24(i,3)
764 1 +my13(i)*vhi(i,3)
765 2 +rrxyz(1,1)+rrxyz(1,2)+rrxyz(1,3)+rrxyz(1,4)
766 ar(2)= z1(i)*vhi(i,1)-x13(i)*v13(i,3)-x24(i)*v24(i,3)
767 1 -mx13(i)*vhi(i,3)
768 2 +rrxyz(2,1)+rrxyz(2,2)+rrxyz(2,3)+rrxyz(2,4)
769 ar(3)= x13(i)*v13(i,2)+x24(i)*v24(i,2)+mx13(i)*vhi(i,2)
770 1 -y13(i)*v13(i,1)-y24(i)*v24(i,1)-my13(i)*vhi(i,1)
771 2 +rrxyz(3,1)+rrxyz(3,2)+rrxyz(3,3)+rrxyz(3,4)
772 ad(1)= vqn(i,1,1)*rrxyz(1,1)+vqn(i,2,1)*rrxyz(2,1)+vqn(i,3,1)*rrxyz(3,1)
773 ad(2)= vqn(i,1,2)*rrxyz(1,2)+vqn(i,2,2)*rrxyz(2,2)+vqn(i,3,2)*rrxyz(3,2)
774 ad(3)= vqn(i,1,3)*rrxyz(1,3)+vqn(i,2,3)*rrxyz(2,3)+vqn(i,3,3)*rrxyz(3,3)
775 ad(4)= vqn(i,1,4)*rrxyz(1,4)+vqn(i,2,4)*rrxyz(2,4)+vqn(i,3,4)*rrxyz(3,4)
776C
777C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
778 xx = corelv(i,1,1)*corelv(i,1,1)+corelv(i,1,2)*corelv(i,1,2)+corelv(i,1,3)*corelv(i,1,3)+corelv(i,1,4)*corelv(i,1,4)
779 yy = corelv(i,2,1)*corelv(i,2,1)+corelv(i,2,2)*corelv(i,2,2)+corelv(i,2,3)*corelv(i,2,3)+corelv(i,2,4)*corelv(i,2,4)
780 xy = corelv(i,1,1)*corelv(i,2,1)+corelv(i,1,2)*corelv(i,2,2)+corelv(i,1,3)*corelv(i,2,3)+corelv(i,1,4)*corelv(i,2,4)
781 zz = four*z2(i)
782 c1 = z1(i)/a_4
783 btb_c = two *c1*c1
784 btb(1)= btb_c*(y24_2(i)+y13_2(i))
785 btb(2)= btb_c*(x24_2(i)+x13_2(i))
786 btb(3)= btb_c*(x24(i)*y24(i)+x13(i)*y13(i))
787C
788 d(1)= yy+zz+4-btb(1)
789 d(2)= xx+zz+4-btb(2)
790 d(3)= -xy+btb(3)
791 deta = d(1)*d(2)-d(3)*d(3)
792 IF (deta<=em20) THEN
793 z1(i)=zero
794 plat(i)=.true.
795 ELSE
796 deta = one/deta
797 di(i,1) = d(2)*deta
798 di(i,2) = d(1)*deta
799 di(i,4) =-d(3)*deta
800 di(i,3) = one/max(xx+yy,em20)
801 di(i,5) = zero
802 di(i,6) = zero
803C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
804 DO j=1,nnod
805 db(i,1,j)= di(i,1)*vqn(i,1,j)+di(i,4)*vqn(i,2,j)
806 db(i,2,j)= di(i,4)*vqn(i,1,j)+di(i,2)*vqn(i,2,j)
807 db(i,3,j)= di(i,3)*vqn(i,3,j)
808 ENDDO
809C
810 dbad(1)= db(i,1,1)*ad(1)+db(i,1,2)*ad(2)+db(i,1,3)*ad(3)+db(i,1,4)*ad(4)
811 dbad(2)= db(i,2,1)*ad(1)+db(i,2,2)*ad(2)+db(i,2,3)*ad(3)+db(i,2,4)*ad(4)
812 dbad(3)= db(i,3,1)*ad(1)+db(i,3,2)*ad(2)+db(i,3,3)*ad(3)+db(i,3,4)*ad(4)
813C
814 alr(1) =di(i,1)*ar(1)+di(i,4)*ar(2)-dbad(1)
815 alr(2) =di(i,4)*ar(1)+di(i,2)*ar(2)-dbad(2)
816 alr(3) = di(i,3)*ar(3)-dbad(3)
817C
818 ald(1) = ad(1)+vqn(i,1,1)*dbad(1)+vqn(i,2,1)*dbad(2)
819 1 +vqn(i,3,1)*dbad(3)
820 2 -db(i,1,1)*ar(1)-db(i,2,1)*ar(2)-db(i,3,1)*ar(3)
821 ald(2) = ad(2)+vqn(i,1,2)*dbad(1)+vqn(i,2,2)*dbad(2)
822 1 +vqn(i,3,2)*dbad(3)
823 2 -db(i,1,2)*ar(1)-db(i,2,2)*ar(2)-db(i,3,2)*ar(3)
824 ald(3) = ad(3)+vqn(i,1,3)*dbad(1)+vqn(i,2,3)*dbad(2)
825 1 +vqn(i,3,3)*dbad(3)
826 2 -db(i,1,3)*ar(1)-db(i,2,3)*ar(2)-db(i,3,3)*ar(3)
827 ald(4) = ad(4)+vqn(i,1,4)*dbad(1)+vqn(i,2,4)*dbad(2)
828 1 +vqn(i,3,4)*dbad(3)
829 2 -db(i,1,4)*ar(1)-db(i,2,4)*ar(2)-db(i,3,4)*ar(3)
830C
831 c1 = 2.0*alr(3)
832 v13(i,1)= v13(i,1)+c1*y13(i)
833 v24(i,1)= v24(i,1)+c1*y24(i)
834 vhi(i,1)= vhi(i,1)+four*(alr(3)*my13(i)-z1(i)*alr(2))
835 v13(i,2)= v13(i,2)-c1*x13(i)
836 v24(i,2)= v24(i,2)-c1*x24(i)
837 vhi(i,2)= vhi(i,2)-four*(alr(3)*mx13(i)-z1(i)*alr(1))
838 v13(i,3)= v13(i,3)-two*(y13(i)*alr(1)-x13(i)*alr(2))
839 v24(i,3)= v24(i,3)-two*(y24(i)*alr(1)-x24(i)*alr(2))
840 vhi(i,3)= vhi(i,3)+four*(mx13(i)*alr(2)-my13(i)*alr(1))
841C
842 rlxyzv(i,1,1)= rrxyz(1,1)-alr(1)-vqn(i,1,1)*ald(1)
843 rlxyzv(i,1,2)= rrxyz(1,2)-alr(1)-vqn(i,1,2)*ald(2)
844 rlxyzv(i,1,3)= rrxyz(1,3)-alr(1)-vqn(i,1,3)*ald(3)
845 rlxyzv(i,1,4)= rrxyz(1,4)-alr(1)-vqn(i,1,4)*ald(4)
846C
847 rlxyzv(i,2,1)= rrxyz(2,1)-alr(2)-vqn(i,2,1)*ald(1)
848 rlxyzv(i,2,2)= rrxyz(2,2)-alr(2)-vqn(i,2,2)*ald(2)
849 rlxyzv(i,2,3)= rrxyz(2,3)-alr(2)-vqn(i,2,3)*ald(3)
850 rlxyzv(i,2,4)= rrxyz(2,4)-alr(2)-vqn(i,2,4)*ald(4)
851 ENDIF
852C
853 END IF !IF(Z2(I)<LL(I)*TOL.OR.NPT == 1)
854C
855 ENDDO
856C
857 RETURN
858 END
859!||====================================================================
860!|| czcorp4 ../engine/source/elements/shell/coquez/czcorc.F
861!||====================================================================
862 SUBROUTINE czcorp4(JFT ,JLT ,VR ,NPT ,TOL ,
863 2 IXC ,PLAT ,AREA ,AREA_I ,V13 ,
864 3 V24 ,VHI ,RLXYZ ,VQN ,VQ ,
865 4 X13 ,X24 ,Y13 ,Y24 ,MX13 ,
866 6 MX23 ,MX34 ,MY13 ,MY23 ,MY34 ,
867 7 Z1 ,DI ,DB ,COREL ,RLZ ,
868 8 LL ,X13_2 ,Y13_2 ,X24_2 ,Y24_2 ,
869 9 L13 ,L24 )
870C-----------------------------------------------
871C I m p l i c i t T y p e s
872C-----------------------------------------------
873#include "implicit_f.inc"
874C-----------------------------------------------
875C G l o b a l P a r a m e t e r s
876C-----------------------------------------------
877#include "mvsiz_p.inc"
878C-----------------------------------------------
879C D u m m y A r g u m e n t s
880C-----------------------------------------------
881 LOGICAL PLAT(*)
882 INTEGER IXC(NIXC,*),JFT,JLT,NPT
883 my_real
884 . vr(3,*),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
885 . mx23(*),my13(*),my23(*),my34(*),
886 . x13(*),x24(*),y13(*),y24(*),mx13(*),
887 . vq(mvsiz,3,3),area(*),z1(*),mx34(*),vqn(mvsiz,3,4),area_i(*),
888 . di(mvsiz,6),db(mvsiz,3,4),ll(*),l13(*),l24(*),
889 . tol,x13_2(*),y13_2(*),x24_2(*),y24_2(*),
890 . rlxyz(mvsiz,2,4),corel(mvsiz,2,4),rlz(mvsiz,4)
891C-----------------------------------------------
892C L O C A L V A R I A B L E S
893C-----------------------------------------------
894 INTEGER NNOD,I,J,K,L
895 parameter(nnod = 4)
896 my_real
897 . deta,
898 . vcore(3,nnod),rrxyz(3,nnod),vg13(3),vg24(3),vghi(3),
899 . z2(mvsiz),a_4,sz,sz1,sz2,sl,c1,c2,s1,
900 . ar(3),ad(nnod),btb(6),xx,yy,zz,xy,xz,yz,
901 . abc,xxyz2,yyxz2,zzxy2,d(6),
902 . alr(3),ald(nnod),dbad(3),btb_c
903C----------------------------------
904 DO i=jft,jlt
905 z2(i)=z1(i)*z1(i)
906 IF (z2(i)<ll(i)*tol.OR.npt == 1) THEN
907 z1(i)=zero
908 plat(i)=.true.
909 ELSE
910 plat(i)=.false.
911C--------------------------------------------------
912C WARPING SPECIAL TREATMENT
913C full projection eliminer drilling rotations and rigid rotations
914C--------------------------------------------------------------------------
915 a_4=area(i)*fourth
916C---------- node N ----------
917 sz1=mx13(i)*y24(i)-my13(i)*x24(i)
918 sz2=a_4+sz1
919 sz=z2(i)*l24(i)
920 sl=one/sqrt(sz+sz2*sz2)
921 vqn(i,1,1)=-z1(i)*y24(i)
922 vqn(i,2,1)= z1(i)*x24(i)
923 vqn(i,3,1)= sz2*sl
924 vqn(i,1,3)=-vqn(i,1,1)
925 vqn(i,2,3)=-vqn(i,2,1)
926 vqn(i,1,1)= vqn(i,1,1)*sl
927 vqn(i,2,1)= vqn(i,2,1)*sl
928C
929 sz2=a_4-sz1
930 sl=one/sqrt(sz+sz2*sz2)
931 vqn(i,1,3)= vqn(i,1,3)*sl
932 vqn(i,2,3)= vqn(i,2,3)*sl
933 vqn(i,3,3)= sz2*sl
934C
935 sz1=mx13(i)*y13(i)-my13(i)*x13(i)
936 sz2=a_4+sz1
937 sz=z2(i)*l13(i)
938 sl=one/sqrt(sz+sz2*sz2)
939 vqn(i,1,2)=-z1(i)*y13(i)
940 vqn(i,2,2)= z1(i)*x13(i)
941 vqn(i,3,2)= sz2*sl
942 vqn(i,1,4)=-vqn(i,1,2)
943 vqn(i,2,4)=-vqn(i,2,2)
944 vqn(i,1,2)= vqn(i,1,2)*sl
945 vqn(i,2,2)= vqn(i,2,2)*sl
946C
947 sz2=a_4-sz1
948 sl=one/sqrt(sz+sz2*sz2)
949 vqn(i,1,4)= vqn(i,1,4)*sl
950 vqn(i,2,4)= vqn(i,2,4)*sl
951 vqn(i,3,4)= sz2*sl
952C
953 k=ixc(2,i)
954 rrxyz(1,1) =rlxyz(i,1,1)
955 rrxyz(2,1) =rlxyz(i,2,1)
956 rrxyz(3,1) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k)
957
958 k=ixc(3,i)
959 rrxyz(1,2) =rlxyz(i,1,2)
960 rrxyz(2,2) =rlxyz(i,2,2)
961 rrxyz(3,2) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k)
962
963 k=ixc(4,i)
964 rrxyz(1,3) =rlxyz(i,1,3)
965 rrxyz(2,3) =rlxyz(i,2,3)
966 rrxyz(3,3) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k)
967
968 k=ixc(5,i)
969 rrxyz(1,4) =rlxyz(i,1,4)
970 rrxyz(2,4) =rlxyz(i,2,4)
971 rrxyz(3,4) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)+vq(i,3,3)*vr(3,k)
972
973C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
974C-----------------------full projection------------------
975 ar(1)=-z1(i)*vhi(i,2)+y13(i)*v13(i,3)+y24(i)*v24(i,3)
976 1 +my13(i)*vhi(i,3)
977 2 +rrxyz(1,1)+rrxyz(1,2)+rrxyz(1,3)+rrxyz(1,4)
978 ar(2)= z1(i)*vhi(i,1)-x13(i)*v13(i,3)-x24(i)*v24(i,3)
979 1 -mx13(i)*vhi(i,3)
980 2 +rrxyz(2,1)+rrxyz(2,2)+rrxyz(2,3)+rrxyz(2,4)
981 ar(3)= x13(i)*v13(i,2)+x24(i)*v24(i,2)+mx13(i)*vhi(i,2)
982 1 -y13(i)*v13(i,1)-y24(i)*v24(i,1)-my13(i)*vhi(i,1)
983 2 +rrxyz(3,1)+rrxyz(3,2)+rrxyz(3,3)+rrxyz(3,4)
984 ad(1)= vqn(i,1,1)*rrxyz(1,1)+vqn(i,2,1)*rrxyz(2,1)+
985 1 vqn(i,3,1)*rrxyz(3,1)
986 ad(2)= vqn(i,1,2)*rrxyz(1,2)+vqn(i,2,2)*rrxyz(2,2)+
987 1 vqn(i,3,2)*rrxyz(3,2)
988 ad(3)= vqn(i,1,3)*rrxyz(1,3)+vqn(i,2,3)*rrxyz(2,3)+
989 1 vqn(i,3,3)*rrxyz(3,3)
990 ad(4)= vqn(i,1,4)*rrxyz(1,4)+vqn(i,2,4)*rrxyz(2,4)+
991 1 vqn(i,3,4)*rrxyz(3,4)
992C
993C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
994 xx = corel(i,1,1)*corel(i,1,1)+corel(i,1,2)*corel(i,1,2)
995 1 +corel(i,1,3)*corel(i,1,3)+corel(i,1,4)*corel(i,1,4)
996 yy = corel(i,2,1)*corel(i,2,1)+corel(i,2,2)*corel(i,2,2)
997 1 +corel(i,2,3)*corel(i,2,3)+corel(i,2,4)*corel(i,2,4)
998 xy = corel(i,1,1)*corel(i,2,1)+corel(i,1,2)*corel(i,2,2)
999 1 +corel(i,1,3)*corel(i,2,3)+corel(i,1,4)*corel(i,2,4)
1000 zz = four*z2(i)
1001 c1 = z1(i)/a_4
1002 btb_c = two *c1*c1
1003 btb(1)= btb_c*(y24_2(i)+y13_2(i))
1004 btb(2)= btb_c*(x24_2(i)+x13_2(i))
1005 btb(3)= btb_c*(x24(i)*y24(i)+x13(i)*y13(i))
1006C
1007 d(1)= yy+zz+4-btb(1)
1008 d(2)= xx+zz+4-btb(2)
1009 d(3)= -xy+btb(3)
1010 deta = d(1)*d(2)-d(3)*d(3)
1011 IF (deta<=em20) THEN
1012 z1(i)=zero
1013 plat(i)=.true.
1014 ELSE
1015 deta = one/deta
1016 di(i,1) = d(2)*deta
1017 di(i,2) = d(1)*deta
1018 di(i,4) =-d(3)*deta
1019 di(i,3) = one/max(xx+yy,em20)
1020 di(i,5) = zero
1021 di(i,6) = zero
1022C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1023 DO j=1,nnod
1024 db(i,1,j)= di(i,1)*vqn(i,1,j)+di(i,4)*vqn(i,2,j)
1025 db(i,2,j)= di(i,4)*vqn(i,1,j)+di(i,2)*vqn(i,2,j)
1026 db(i,3,j)= di(i,3)*vqn(i,3,j)
1027 ENDDO
1028C
1029 dbad(1)= db(i,1,1)*ad(1)+db(i,1,2)*ad(2)
1030 1 +db(i,1,3)*ad(3)+db(i,1,4)*ad(4)
1031 dbad(2)= db(i,2,1)*ad(1)+db(i,2,2)*ad(2)
1032 1 +db(i,2,3)*ad(3)+db(i,2,4)*ad(4)
1033 dbad(3)= db(i,3,1)*ad(1)+db(i,3,2)*ad(2)
1034 1 +db(i,3,3)*ad(3)+db(i,3,4)*ad(4)
1035C
1036 alr(1) =di(i,1)*ar(1)+di(i,4)*ar(2)-dbad(1)
1037 alr(2) =di(i,4)*ar(1)+di(i,2)*ar(2)-dbad(2)
1038 alr(3) = di(i,3)*ar(3)-dbad(3)
1039C
1040 ald(1) = ad(1)+vqn(i,1,1)*dbad(1)+vqn(i,2,1)*dbad(2)
1041 1 +vqn(i,3,1)*dbad(3)
1042 2 -db(i,1,1)*ar(1)-db(i,2,1)*ar(2)-db(i,3,1)*ar(3)
1043 ald(2) = ad(2)+vqn(i,1,2)*dbad(1)+vqn(i,2,2)*dbad(2)
1044 1 +vqn(i,3,2)*dbad(3)
1045 2 -db(i,1,2)*ar(1)-db(i,2,2)*ar(2)-db(i,3,2)*ar(3)
1046 ald(3) = ad(3)+vqn(i,1,3)*dbad(1)+vqn(i,2,3)*dbad(2)
1047 1 +vqn(i,3,3)*dbad(3)
1048 2 -db(i,1,3)*ar(1)-db(i,2,3)*ar(2)-db(i,3,3)*ar(3)
1049 ald(4) = ad(4)+vqn(i,1,4)*dbad(1)+vqn(i,2,4)*dbad(2)
1050 1 +vqn(i,3,4)*dbad(3)
1051 2 -db(i,1,4)*ar(1)-db(i,2,4)*ar(2)-db(i,3,4)*ar(3)
1052C
1053C
1054 c1 = two*alr(3)
1055 v13(i,1)= v13(i,1)+c1*y13(i)
1056 v24(i,1)= v24(i,1)+c1*y24(i)
1057 vhi(i,1)= vhi(i,1)+four*(alr(3)*my13(i)-z1(i)*alr(2))
1058 v13(i,2)= v13(i,2)-c1*x13(i)
1059 v24(i,2)= v24(i,2)-c1*x24(i)
1060 vhi(i,2)= vhi(i,2)-four*(alr(3)*mx13(i)-z1(i)*alr(1))
1061 v13(i,3)= v13(i,3)-two*(y13(i)*alr(1)-x13(i)*alr(2))
1062 v24(i,3)= v24(i,3)-two*(y24(i)*alr(1)-x24(i)*alr(2))
1063 vhi(i,3)= vhi(i,3)+four*(mx13(i)*alr(2)-my13(i)*alr(1))
1064 rlxyz(i,1,1)= rrxyz(1,1)-alr(1)-vqn(i,1,1)*ald(1)
1065 rlxyz(i,1,2)= rrxyz(1,2)-alr(1)-vqn(i,1,2)*ald(2)
1066 rlxyz(i,1,3)= rrxyz(1,3)-alr(1)-vqn(i,1,3)*ald(3)
1067 rlxyz(i,1,4)= rrxyz(1,4)-alr(1)-vqn(i,1,4)*ald(4)
1068C
1069 rlxyz(i,2,1)= rrxyz(2,1)-alr(2)-vqn(i,2,1)*ald(1)
1070 rlxyz(i,2,2)= rrxyz(2,2)-alr(2)-vqn(i,2,2)*ald(2)
1071 rlxyz(i,2,3)= rrxyz(2,3)-alr(2)-vqn(i,2,3)*ald(3)
1072 rlxyz(i,2,4)= rrxyz(2,4)-alr(2)-vqn(i,2,4)*ald(4)
1073 ENDIF
1074C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1075C
1076 ENDIF
1077C
1078 ENDDO
1079 RETURN
1080 END
1081!||====================================================================
1082!|| czcorp4x ../engine/source/elements/shell/coquez/czcorc.F
1083!||====================================================================
1084 SUBROUTINE czcorp4x(JFT ,JLT ,VR ,NPT ,TOL ,
1085 2 IXC ,PLAT ,AREA ,AREA_I ,V13 ,
1086 3 V24 ,VHI ,RLXYZ ,VQN ,VQ ,
1087 4 X13 ,X24 ,Y13 ,Y24 ,MX13 ,
1088 6 MX23 ,MX34 ,MY13 ,MY23 ,MY34 ,
1089 7 Z1 ,DI ,DB ,COREL ,RLZ ,
1090 8 LL ,X13_2 ,Y13_2 ,X24_2 ,Y24_2 ,
1091 9 L13 ,L24 ,VRX1 ,VRX2 ,VRX3 ,
1092 A VRX4 ,VRY1 ,VRY2 ,VRY3 ,VRY4 ,
1093 B VRZ1 ,VRZ2 ,VRZ3 ,VRZ4 )
1094C-----------------------------------------------
1095C I m p l i c i t T y p e s
1096C-----------------------------------------------
1097#include "implicit_f.inc"
1098C-----------------------------------------------
1099C G l o b a l P a r a m e t e r s
1100C-----------------------------------------------
1101#include "mvsiz_p.inc"
1102C-----------------------------------------------
1103C D u m m y A r g u m e n t s
1104C-----------------------------------------------
1105 LOGICAL PLAT(*)
1106 INTEGER IXC(NIXC,*),JFT,JLT,NPT
1107 my_real
1108 . vr(3,*),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
1109 . mx23(*),my13(*),my23(*),my34(*),
1110 . x13(*),x24(*),y13(*),y24(*),mx13(*),
1111 . vq(mvsiz,3,3),area(*),z1(*),mx34(*),vqn(mvsiz,3,4),area_i(*),
1112 . di(mvsiz,6),db(mvsiz,3,4),ll(*),l13(*),l24(*),
1113 . tol,x13_2(*),y13_2(*),x24_2(*),y24_2(*),
1114 . rlxyz(mvsiz,2,4),corel(mvsiz,2,4),rlz(mvsiz,4),
1115 . vrx1(*),vrx2(*),vrx3(*),vrx4(*),
1116 . vry1(*),vry2(*),vry3(*),vry4(*),
1117 . vrz1(*),vrz2(*),vrz3(*),vrz4(*)
1118C-----------------------------------------------
1119C L O C A L V A R I A B L E S
1120C-----------------------------------------------
1121 INTEGER NNOD,I,J,K,L
1122 PARAMETER (NNOD = 4)
1123 my_real
1124 . deta,
1125 . vcore(3,nnod),rrxyz(3,nnod),vg13(3),vg24(3),vghi(3),
1126 . z2(mvsiz),a_4,sz,sz1,sz2,sl,c1,c2,s1,
1127 . ar(3),ad(nnod),btb(6),xx,yy,zz,xy,xz,yz,
1128 . abc,xxyz2,yyxz2,zzxy2,d(6),
1129 . alr(3),ald(nnod),dbad(3),btb_c
1130C----------------------------------
1131 DO i=jft,jlt
1132 z2(i)=z1(i)*z1(i)
1133 IF (z2(i)<ll(i)*tol.OR.npt == 1) THEN
1134 z1(i)=zero
1135 plat(i)=.true.
1136 ELSE
1137 plat(i)=.false.
1138C--------------------------------------------------
1139C WAPRING SPECIAL TREATMENT
1140C full projection eliminer drilling rotations and rigid rotations
1141C--------------------------------------------------------------------------
1142 a_4=area(i)*fourth
1143C
1144C---------- node N ----------
1145 sz1=mx13(i)*y24(i)-my13(i)*x24(i)
1146 sz2=a_4+sz1
1147 sz=z2(i)*l24(i)
1148 sl=one/sqrt(sz+sz2*sz2)
1149 vqn(i,1,1)=-z1(i)*y24(i)
1150 vqn(i,2,1)= z1(i)*x24(i)
1151 vqn(i,3,1)= sz2*sl
1152 vqn(i,1,3)=-vqn(i,1,1)
1153 vqn(i,2,3)=-vqn(i,2,1)
1154 vqn(i,1,1)= vqn(i,1,1)*sl
1155 vqn(i,2,1)= vqn(i,2,1)*sl
1156C
1157 sz2=a_4-sz1
1158 sl=one/sqrt(sz+sz2*sz2)
1159 vqn(i,1,3)= vqn(i,1,3)*sl
1160 vqn(i,2,3)= vqn(i,2,3)*sl
1161 vqn(i,3,3)= sz2*sl
1162C
1163 sz1=mx13(i)*y13(i)-my13(i)*x13(i)
1164 sz2=a_4+sz1
1165 sz=z2(i)*l13(i)
1166 sl=one/sqrt(sz+sz2*sz2)
1167 vqn(i,1,2)=-z1(i)*y13(i)
1168 vqn(i,2,2)= z1(i)*x13(i)
1169 vqn(i,3,2)= sz2*sl
1170 vqn(i,1,4)=-vqn(i,1,2)
1171 vqn(i,2,4)=-vqn(i,2,2)
1172 vqn(i,1,2)= vqn(i,1,2)*sl
1173 vqn(i,2,2)= vqn(i,2,2)*sl
1174C
1175 sz2=a_4-sz1
1176 sl=one/sqrt(sz+sz2*sz2)
1177 vqn(i,1,4)= vqn(i,1,4)*sl
1178 vqn(i,2,4)= vqn(i,2,4)*sl
1179 vqn(i,3,4)= sz2*sl
1180C
1181 k=ixc(2,i)
1182 rrxyz(1,1) =rlxyz(i,1,1)
1183 rrxyz(2,1) =rlxyz(i,2,1)
1184 rrxyz(3,1) =vq(i,1,3)*vrx1(i)+vq(i,2,3)*vry1(i)
1185 1 +vq(i,3,3)*vrz1(i)
1186 k=ixc(3,i)
1187 rrxyz(1,2) =rlxyz(i,1,2)
1188 rrxyz(2,2) =rlxyz(i,2,2)
1189 rrxyz(3,2) =vq(i,1,3)*vrx2(i)+vq(i,2,3)*vry2(i)
1190 1 +vq(i,3,3)*vrz2(i)
1191 k=ixc(4,i)
1192 rrxyz(1,3) =rlxyz(i,1,3)
1193 rrxyz(2,3) =rlxyz(i,2,3)
1194 rrxyz(3,3) =vq(i,1,3)*vrx3(i)+vq(i,2,3)*vry3(i)
1195 1 +vq(i,3,3)*vrz3(i)
1196 k=ixc(5,i)
1197 rrxyz(1,4) =rlxyz(i,1,4)
1198 rrxyz(2,4) =rlxyz(i,2,4)
1199 rrxyz(3,4) =vq(i,1,3)*vrx4(i)+vq(i,2,3)*vry4(i)
1200 1 +vq(i,3,3)*vrz4(i)
1201C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1202C-----------------------full projection------------------
1203 ar(1)=-z1(i)*vhi(i,2)+y13(i)*v13(i,3)+y24(i)*v24(i,3)
1204 1 +my13(i)*vhi(i,3)
1205 2 +rrxyz(1,1)+rrxyz(1,2)+rrxyz(1,3)+rrxyz(1,4)
1206 ar(2)= z1(i)*vhi(i,1)-x13(i)*v13(i,3)-x24(i)*v24(i,3)
1207 1 -mx13(i)*vhi(i,3)
1208 2 +rrxyz(2,1)+rrxyz(2,2)+rrxyz(2,3)+rrxyz(2,4)
1209 ar(3)= x13(i)*v13(i,2)+x24(i)*v24(i,2)+mx13(i)*vhi(i,2)
1210 1 -y13(i)*v13(i,1)-y24(i)*v24(i,1)-my13(i)*vhi(i,1)
1211 2 +rrxyz(3,1)+rrxyz(3,2)+rrxyz(3,3)+rrxyz(3,4)
1212 ad(1)= vqn(i,1,1)*rrxyz(1,1)+vqn(i,2,1)*rrxyz(2,1)+
1213 1 vqn(i,3,1)*rrxyz(3,1)
1214 ad(2)= vqn(i,1,2)*rrxyz(1,2)+vqn(i,2,2)*rrxyz(2,2)+
1215 1 vqn(i,3,2)*rrxyz(3,2)
1216 ad(3)= vqn(i,1,3)*rrxyz(1,3)+vqn(i,2,3)*rrxyz(2,3)+
1217 1 vqn(i,3,3)*rrxyz(3,3)
1218 ad(4)= vqn(i,1,4)*rrxyz(1,4)+vqn(i,2,4)*rrxyz(2,4)+
1219 1 vqn(i,3,4)*rrxyz(3,4)
1220C
1221C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1222 xx = corel(i,1,1)*corel(i,1,1)+corel(i,1,2)*corel(i,1,2)
1223 1 +corel(i,1,3)*corel(i,1,3)+corel(i,1,4)*corel(i,1,4)
1224 yy = corel(i,2,1)*corel(i,2,1)+corel(i,2,2)*corel(i,2,2)
1225 1 +corel(i,2,3)*corel(i,2,3)+corel(i,2,4)*corel(i,2,4)
1226 xy = corel(i,1,1)*corel(i,2,1)+corel(i,1,2)*corel(i,2,2)
1227 1 +corel(i,1,3)*corel(i,2,3)+corel(i,1,4)*corel(i,2,4)
1228 zz = four*z2(i)
1229 c1 = z1(i)/a_4
1230 btb_c = two *c1*c1
1231 btb(1)= btb_c*(y24_2(i)+y13_2(i))
1232 btb(2)= btb_c*(x24_2(i)+x13_2(i))
1233 btb(3)= btb_c*(x24(i)*y24(i)+x13(i)*y13(i))
1234C
1235 d(1)= yy+zz+4-btb(1)
1236 d(2)= xx+zz+4-btb(2)
1237 d(3)= -xy+btb(3)
1238 deta = d(1)*d(2)-d(3)*d(3)
1239 IF (deta<=em20) THEN
1240 z1(i)=zero
1241 plat(i)=.true.
1242 ELSE
1243 deta = one/deta
1244 di(i,1) = d(2)*deta
1245 di(i,2) = d(1)*deta
1246 di(i,4) =-d(3)*deta
1247 di(i,3) = one/max(xx+yy,em20)
1248 di(i,5) = zero
1249 di(i,6) = zero
1250C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1251 DO j=1,nnod
1252 db(i,1,j)= di(i,1)*vqn(i,1,j)+di(i,4)*vqn(i,2,j)
1253 db(i,2,j)= di(i,4)*vqn(i,1,j)+di(i,2)*vqn(i,2,j)
1254 db(i,3,j)= di(i,3)*vqn(i,3,j)
1255 ENDDO
1256C
1257 dbad(1)= db(i,1,1)*ad(1)+db(i,1,2)*ad(2)
1258 1 +db(i,1,3)*ad(3)+db(i,1,4)*ad(4)
1259 dbad(2)= db(i,2,1)*ad(1)+db(i,2,2)*ad(2)
1260 1 +db(i,2,3)*ad(3)+db(i,2,4)*ad(4)
1261 dbad(3)= db(i,3,1)*ad(1)+db(i,3,2)*ad(2)
1262 1 +db(i,3,3)*ad(3)+db(i,3,4)*ad(4)
1263C
1264 alr(1) =di(i,1)*ar(1)+di(i,4)*ar(2)-dbad(1)
1265 alr(2) =di(i,4)*ar(1)+di(i,2)*ar(2)-dbad(2)
1266 alr(3) = di(i,3)*ar(3)-dbad(3)
1267C
1268 ald(1) = ad(1)+vqn(i,1,1)*dbad(1)+vqn(i,2,1)*dbad(2)
1269 1 +vqn(i,3,1)*dbad(3)
1270 2 -db(i,1,1)*ar(1)-db(i,2,1)*ar(2)-db(i,3,1)*ar(3)
1271 ald(2) = ad(2)+vqn(i,1,2)*dbad(1)+vqn(i,2,2)*dbad(2)
1272 1 +vqn(i,3,2)*dbad(3)
1273 2 -db(i,1,2)*ar(1)-db(i,2,2)*ar(2)-db(i,3,2)*ar(3)
1274 ald(3) = ad(3)+vqn(i,1,3)*dbad(1)+vqn(i,2,3)*dbad(2)
1275 1 +vqn(i,3,3)*dbad(3)
1276 2 -db(i,1,3)*ar(1)-db(i,2,3)*ar(2)-db(i,3,3)*ar(3)
1277 ald(4) = ad(4)+vqn(i,1,4)*dbad(1)+vqn(i,2,4)*dbad(2)
1278 1 +vqn(i,3,4)*dbad(3)
1279 2 -db(i,1,4)*ar(1)-db(i,2,4)*ar(2)-db(i,3,4)*ar(3)
1280C
1281C
1282 c1 = two*alr(3)
1283 v13(i,1)= v13(i,1)+c1*y13(i)
1284 v24(i,1)= v24(i,1)+c1*y24(i)
1285 vhi(i,1)= vhi(i,1)+four*(alr(3)*my13(i)-z1(i)*alr(2))
1286 v13(i,2)= v13(i,2)-c1*x13(i)
1287 v24(i,2)= v24(i,2)-c1*x24(i)
1288 vhi(i,2)= vhi(i,2)-four*(alr(3)*mx13(i)-z1(i)*alr(1))
1289 v13(i,3)= v13(i,3)-two*(y13(i)*alr(1)-x13(i)*alr(2))
1290 v24(i,3)= v24(i,3)-two*(y24(i)*alr(1)-x24(i)*alr(2))
1291 vhi(i,3)= vhi(i,3)+four*(mx13(i)*alr(2)-my13(i)*alr(1))
1292 rlxyz(i,1,1)= rrxyz(1,1)-alr(1)-vqn(i,1,1)*ald(1)
1293 rlxyz(i,1,2)= rrxyz(1,2)-alr(1)-vqn(i,1,2)*ald(2)
1294 rlxyz(i,1,3)= rrxyz(1,3)-alr(1)-vqn(i,1,3)*ald(3)
1295 rlxyz(i,1,4)= rrxyz(1,4)-alr(1)-vqn(i,1,4)*ald(4)
1296C
1297 rlxyz(i,2,1)= rrxyz(2,1)-alr(2)-vqn(i,2,1)*ald(1)
1298 rlxyz(i,2,2)= rrxyz(2,2)-alr(2)-vqn(i,2,2)*ald(2)
1299 rlxyz(i,2,3)= rrxyz(2,3)-alr(2)-vqn(i,2,3)*ald(3)
1300 rlxyz(i,2,4)= rrxyz(2,4)-alr(2)-vqn(i,2,4)*ald(4)
1301 ENDIF
1302C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1303C
1304 ENDIF
1305C
1306 ENDDO
1307C
1308 RETURN
1309 END
1310!||====================================================================
1311!|| czcorp5v ../engine/source/elements/shell/coquez/czcorc.F
1312!||--- calls -----------------------------------------------------
1313!|| a3invdp ../engine/source/elements/shell/coquez/czcorc.F
1314!||====================================================================
1315 SUBROUTINE czcorp5v(JFT ,JLT ,VR ,NPT ,TOL ,
1316 2 IXC ,PLAT ,AREA ,AREA_I ,V13 ,
1317 3 V24 ,VHI ,RLXYZV ,VQN ,VQ ,
1318 4 X13 ,X24 ,Y13 ,Y24 ,MX13 ,
1319 6 MX23 ,MX34 ,MY13 ,MY23 ,MY34 ,
1320 7 Z1 ,DI ,DB ,CORELV ,RLZ ,
1321 8 LL ,X13_2 ,Y13_2 ,X24_2 ,Y24_2 ,
1322 9 L13 ,L24 ,IDRIL ,DIZ )
1323C-----------------------------------------------
1324C I m p l i c i t T y p e s
1325C-----------------------------------------------
1326#include "implicit_f.inc"
1327C-----------------------------------------------
1328C G l o b a l P a r a m e t e r s
1329C-----------------------------------------------
1330#include "mvsiz_p.inc"
1331#include "impl1_c.inc"
1332#include "scr05_c.inc"
1333C-----------------------------------------------
1334C D u m m y A r g u m e n t s
1335C-----------------------------------------------
1336 LOGICAL PLAT(*)
1337 INTEGER IXC(NIXC,*),JFT,JLT,IDRIL,NPT
1338 my_real
1339 . vr(3,*),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
1340 . mx23(*),my13(*),my23(*),my34(*),
1341 . x13(*),x24(*),y13(*),y24(*),mx13(*),
1342 . vq(mvsiz,3,3),area(*),z1(*),mx34(*),vqn(mvsiz,3,4),area_i(*),
1343 . di(mvsiz,6),db(mvsiz,3,4),ll(*),l13(*),l24(*),
1344 . tol,x13_2(*),y13_2(*),x24_2(*),y24_2(*),
1345 . rlxyzv(mvsiz,2,4),corelv(mvsiz,2,4),rlz(mvsiz,4),diz(mvsiz,3)
1346C-----------------------------------------------
1347C L O C A L V A R I A B L E S
1348C-----------------------------------------------
1349 INTEGER NNOD,I,J,K,L
1350 PARAMETER (NNOD = 4)
1351 my_real
1352 . deta,
1353 . vcore(3,nnod),rrxyz(3,nnod),vg13(3),vg24(3),vghi(3),
1354 . z2(mvsiz),a_4,sz,sz1,sz2,sl,c1,c2,s1,
1355 . ar(3),ad(nnod),btb(6),xx,yy,zz,xy,xz,yz,
1356 . abc,xxyz2,yyxz2,zzxy2,d(6),diz1(6),diz2(6),
1357 . alr(3),ald(nnod),dbad(3),btb_c,alrz
1358C-----------------------------------------------
1359 DO i=jft,jlt
1360 z2(i)=z1(i)*z1(i)
1361 IF (z2(i)<ll(i)*tol.OR.npt == 1) THEN
1362 z1(i)=zero
1363 plat(i)=.true.
1364 ELSE
1365 plat(i)=.false.
1366C--------------------------------------------------
1367C WARPING SPECIAL TREATMENT
1368C full projection eliminer drilling rotations and rigid rotations
1369C--------------------------------------------------------------------------
1370 a_4=area(i)*fourth
1371C
1372C---------- node N ----------
1373 sz1=mx13(i)*y24(i)-my13(i)*x24(i)
1374 sz2=a_4+sz1
1375 sz=z2(i)*l24(i)
1376 sl=one/sqrt(sz+sz2*sz2)
1377 vqn(i,1,1)=-z1(i)*y24(i)
1378 vqn(i,2,1)= z1(i)*x24(i)
1379 vqn(i,3,1)= sz2*sl
1380 vqn(i,1,3)=-vqn(i,1,1)
1381 vqn(i,2,3)=-vqn(i,2,1)
1382 vqn(i,1,1)= vqn(i,1,1)*sl
1383 vqn(i,2,1)= vqn(i,2,1)*sl
1384C
1385 sz2=a_4-sz1
1386 sl=one/sqrt(sz+sz2*sz2)
1387 vqn(i,1,3)= vqn(i,1,3)*sl
1388 vqn(i,2,3)= vqn(i,2,3)*sl
1389 vqn(i,3,3)= sz2*sl
1390C
1391 sz1=mx13(i)*y13(i)-my13(i)*x13(i)
1392 sz2=a_4+sz1
1393 sz=z2(i)*l13(i)
1394 sl=one/sqrt(sz+sz2*sz2)
1395 vqn(i,1,2)=-z1(i)*y13(i)
1396 vqn(i,2,2)= z1(i)*x13(i)
1397 vqn(i,3,2)= sz2*sl
1398 vqn(i,1,4)=-vqn(i,1,2)
1399 vqn(i,2,4)=-vqn(i,2,2)
1400 vqn(i,1,2)= vqn(i,1,2)*sl
1401 vqn(i,2,2)= vqn(i,2,2)*sl
1402C
1403 sz2=a_4-sz1
1404 sl=one/sqrt(sz+sz2*sz2)
1405 vqn(i,1,4)= vqn(i,1,4)*sl
1406 vqn(i,2,4)= vqn(i,2,4)*sl
1407 vqn(i,3,4)= sz2*sl
1408C
1409 k=ixc(2,i)
1410 rrxyz(1,1) =rlxyzv(i,1,1)
1411 rrxyz(2,1) =rlxyzv(i,2,1)
1412 rrxyz(3,1) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)
1413 1 +vq(i,3,3)*vr(3,k)
1414 k=ixc(3,i)
1415 rrxyz(1,2) =rlxyzv(i,1,2)
1416 rrxyz(2,2) =rlxyzv(i,2,2)
1417 rrxyz(3,2) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)
1418 1 +vq(i,3,3)*vr(3,k)
1419 k=ixc(4,i)
1420 rrxyz(1,3) =rlxyzv(i,1,3)
1421 rrxyz(2,3) =rlxyzv(i,2,3)
1422 rrxyz(3,3) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)
1423 1 +vq(i,3,3)*vr(3,k)
1424 k=ixc(5,i)
1425 rrxyz(1,4) =rlxyzv(i,1,4)
1426 rrxyz(2,4) =rlxyzv(i,2,4)
1427 rrxyz(3,4) =vq(i,1,3)*vr(1,k)+vq(i,2,3)*vr(2,k)
1428 1 +vq(i,3,3)*vr(3,k)
1429 IF (impl_s>0.AND.ikproj<0) THEN
1430C-------------------------------------
1431C DRILLING PROJECTION ONLY
1432C-------------------------------------
1433 IF (idril>0) THEN
1434 rlz(i,1)=area_i(i)*(vqn(i,1,1)*rrxyz(1,1)+
1435 1 vqn(i,2,1)*rrxyz(2,1)+vqn(i,3,1)*rrxyz(3,1))
1436 rlz(i,2)=area_i(i)*(vqn(i,1,2)*rrxyz(1,2)+
1437 1 vqn(i,2,2)*rrxyz(2,2)+vqn(i,3,2)*rrxyz(3,2))
1438 rlz(i,3)=area_i(i)*(vqn(i,1,3)*rrxyz(1,3)+
1439 1 vqn(i,2,3)*rrxyz(2,3)+vqn(i,3,3)*rrxyz(3,3))
1440 rlz(i,4)=area_i(i)*(vqn(i,1,4)*rrxyz(1,4)+
1441 1 vqn(i,2,4)*rrxyz(2,4)+vqn(i,3,4)*rrxyz(3,4))
1442 END IF
1443C
1444 rlxyzv(i,1,1)=(one-vqn(i,1,1)*vqn(i,1,1))*rrxyz(1,1)
1445 1 -vqn(i,1,1)*vqn(i,2,1) *rrxyz(2,1)
1446 2 -vqn(i,1,1)*vqn(i,3,1) *rrxyz(3,1)
1447 rlxyzv(i,2,1)=(one-vqn(i,2,1)*vqn(i,2,1))*rrxyz(2,1)
1448 1 -vqn(i,1,1)*vqn(i,2,1) *rrxyz(1,1)
1449 2 -vqn(i,2,1)*vqn(i,3,1) *rrxyz(3,1)
1450C
1451 rlxyzv(i,1,2)=(one-vqn(i,1,2)*vqn(i,1,2))*rrxyz(1,2)
1452 1 -vqn(i,1,2)*vqn(i,2,2) *rrxyz(2,2)
1453 2 -vqn(i,1,2)*vqn(i,3,2) *rrxyz(3,2)
1454 rlxyzv(i,2,2)=(one-vqn(i,2,2)*vqn(i,2,2))*rrxyz(2,2)
1455 1 -vqn(i,1,2)*vqn(i,2,2) *rrxyz(1,2)
1456 2 -vqn(i,2,2)*vqn(i,3,2) *rrxyz(3,2)
1457C
1458 rlxyzv(i,1,3)=(one-vqn(i,1,3)*vqn(i,1,3))*rrxyz(1,3)
1459 1 -vqn(i,1,3)*vqn(i,2,3) *rrxyz(2,3)
1460 2 -vqn(i,1,3)*vqn(i,3,3) *rrxyz(3,3)
1461 rlxyzv(i,2,3)=(one-vqn(i,2,3)*vqn(i,2,3))*rrxyz(2,3)
1462 1 -vqn(i,1,3)*vqn(i,2,3) *rrxyz(1,3)
1463 2 -vqn(i,2,3)*vqn(i,3,3) *rrxyz(3,3)
1464C
1465 rlxyzv(i,1,4)=(one-vqn(i,1,4)*vqn(i,1,4))*rrxyz(1,4)
1466 1 -vqn(i,1,4)*vqn(i,2,4) *rrxyz(2,4)
1467 2 -vqn(i,1,4)*vqn(i,3,4) *rrxyz(3,4)
1468 rlxyzv(i,2,4)=(one-vqn(i,2,4)*vqn(i,2,4))*rrxyz(2,4)
1469 1 -vqn(i,1,4)*vqn(i,2,4) *rrxyz(1,4)
1470 2 -vqn(i,2,4)*vqn(i,3,4) *rrxyz(3,4)
1471 ELSE
1472C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1473C-----------------------full projection------------------
1474 ar(1)=-z1(i)*vhi(i,2)+y13(i)*v13(i,3)+y24(i)*v24(i,3)
1475 1 +my13(i)*vhi(i,3)
1476 2 +rrxyz(1,1)+rrxyz(1,2)+rrxyz(1,3)+rrxyz(1,4)
1477 ar(2)= z1(i)*vhi(i,1)-x13(i)*v13(i,3)-x24(i)*v24(i,3)
1478 1 -mx13(i)*vhi(i,3)
1479 2 +rrxyz(2,1)+rrxyz(2,2)+rrxyz(2,3)+rrxyz(2,4)
1480 ar(3)= x13(i)*v13(i,2)+x24(i)*v24(i,2)+mx13(i)*vhi(i,2)
1481 1 -y13(i)*v13(i,1)-y24(i)*v24(i,1)-my13(i)*vhi(i,1)
1482 2 +rrxyz(3,1)+rrxyz(3,2)+rrxyz(3,3)+rrxyz(3,4)
1483 ad(1)= vqn(i,1,1)*rrxyz(1,1)+vqn(i,2,1)*rrxyz(2,1)+
1484 1 vqn(i,3,1)*rrxyz(3,1)
1485 ad(2)= vqn(i,1,2)*rrxyz(1,2)+vqn(i,2,2)*rrxyz(2,2)+
1486 1 vqn(i,3,2)*rrxyz(3,2)
1487 ad(3)= vqn(i,1,3)*rrxyz(1,3)+vqn(i,2,3)*rrxyz(2,3)+
1488 1 vqn(i,3,3)*rrxyz(3,3)
1489 ad(4)= vqn(i,1,4)*rrxyz(1,4)+vqn(i,2,4)*rrxyz(2,4)+
1490 1 vqn(i,3,4)*rrxyz(3,4)
1491C
1492C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1493 xx = corelv(i,1,1)*corelv(i,1,1)+corelv(i,1,2)*corelv(i,1,2)
1494 1 +corelv(i,1,3)*corelv(i,1,3)+corelv(i,1,4)*corelv(i,1,4)
1495 yy = corelv(i,2,1)*corelv(i,2,1)+corelv(i,2,2)*corelv(i,2,2)
1496 1 +corelv(i,2,3)*corelv(i,2,3)+corelv(i,2,4)*corelv(i,2,4)
1497 xy = corelv(i,1,1)*corelv(i,2,1)+corelv(i,1,2)*corelv(i,2,2)
1498 1 +corelv(i,1,3)*corelv(i,2,3)+corelv(i,1,4)*corelv(i,2,4)
1499 xz =(corelv(i,1,1)-corelv(i,1,2)+corelv(i,1,3)-corelv(i,1,4))
1500 . *z1(i)
1501 yz =(corelv(i,2,1)-corelv(i,2,2)+corelv(i,2,3)-corelv(i,2,4))
1502 . *z1(i)
1503 zz = four*z2(i)
1504C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1505 btb(1)= vqn(i,1,1)*vqn(i,1,1)+vqn(i,1,2)*vqn(i,1,2)
1506 1 +vqn(i,1,3)*vqn(i,1,3)+vqn(i,1,4)*vqn(i,1,4)
1507 btb(2)= vqn(i,2,1)*vqn(i,2,1)+vqn(i,2,2)*vqn(i,2,2)
1508 1 +vqn(i,2,3)*vqn(i,2,3)+vqn(i,2,4)*vqn(i,2,4)
1509 btb(3)= vqn(i,3,1)*vqn(i,3,1)+vqn(i,3,2)*vqn(i,3,2)
1510 1 +vqn(i,3,3)*vqn(i,3,3)+vqn(i,3,4)*vqn(i,3,4)
1511 btb(4)= vqn(i,1,1)*vqn(i,2,1)+vqn(i,1,2)*vqn(i,2,2)
1512 1 +vqn(i,1,3)*vqn(i,2,3)+vqn(i,1,4)*vqn(i,2,4)
1513 btb(5)= vqn(i,1,1)*vqn(i,3,1)+vqn(i,1,2)*vqn(i,3,2)
1514 1 +vqn(i,1,3)*vqn(i,3,3)+vqn(i,1,4)*vqn(i,3,4)
1515 btb(6)= vqn(i,2,1)*vqn(i,3,1)+vqn(i,2,2)*vqn(i,3,2)
1516 1 +vqn(i,2,3)*vqn(i,3,3)+vqn(i,2,4)*vqn(i,3,4)
1517 d(1)= yy+zz+four-btb(1)
1518 d(2)= xx+zz+four-btb(2)
1519 d(3)= xx+yy+four-btb(3)
1520 d(4)= -xy-btb(4)
1521 d(5)= -xz-btb(5)
1522 d(6)= -yz-btb(6)
1523 IF(iresp == 1)THEN
1524 CALL a3invdp(d,diz2)
1525 di(i,1:6) = diz2(1:6)
1526 ELSE
1527 abc = d(1)*d(2)*d(3)
1528 xxyz2 = d(1)*d(6)*d(6)
1529 yyxz2 = d(2)*d(5)*d(5)
1530 zzxy2 = d(3)*d(4)*d(4)
1531 deta = abs(abc+two*d(4)*d(5)*d(6)-xxyz2-yyxz2-zzxy2)
1532 deta = one/max(deta,em20)
1533 di(i,1) = (abc-xxyz2)*deta/max(d(1),em20)
1534 di(i,2) = (abc-yyxz2)*deta/max(d(2),em20)
1535 di(i,3) = (abc-zzxy2)*deta/max(d(3),em20)
1536 di(i,4) = (d(5)*d(6)-d(4)*d(3))*deta
1537 di(i,5) = (d(6)*d(4)-d(5)*d(2))*deta
1538 di(i,6) = (d(4)*d(5)-d(6)*d(1))*deta
1539 END IF !(IRESP == 1)THEN
1540 DO j=1,nnod
1541 db(i,1,j)= di(i,1)*vqn(i,1,j)+di(i,4)*vqn(i,2,j)
1542 1 +di(i,5)*vqn(i,3,j)
1543 db(i,2,j)= di(i,4)*vqn(i,1,j)+di(i,2)*vqn(i,2,j)
1544 1 +di(i,6)*vqn(i,3,j)
1545 db(i,3,j)= di(i,5)*vqn(i,1,j)+di(i,6)*vqn(i,2,j)
1546 1 +di(i,3)*vqn(i,3,j)
1547 ENDDO
1548C
1549 dbad(1)= db(i,1,1)*ad(1)+db(i,1,2)*ad(2)
1550 1 +db(i,1,3)*ad(3)+db(i,1,4)*ad(4)
1551 dbad(2)= db(i,2,1)*ad(1)+db(i,2,2)*ad(2)
1552 1 +db(i,2,3)*ad(3)+db(i,2,4)*ad(4)
1553 dbad(3)= db(i,3,1)*ad(1)+db(i,3,2)*ad(2)
1554 1 +db(i,3,3)*ad(3)+db(i,3,4)*ad(4)
1555C
1556 alr(1) =di(i,1)*ar(1)+di(i,4)*ar(2)+di(i,5)*ar(3)-dbad(1)
1557 alr(2) =di(i,4)*ar(1)+di(i,2)*ar(2)+di(i,6)*ar(3)-dbad(2)
1558 alr(3) =di(i,5)*ar(1)+di(i,6)*ar(2)+di(i,3)*ar(3)-dbad(3)
1559C
1560 ald(1) = ad(1)+vqn(i,1,1)*dbad(1)+vqn(i,2,1)*dbad(2)
1561 1 +vqn(i,3,1)*dbad(3)
1562 2 -db(i,1,1)*ar(1)-db(i,2,1)*ar(2)-db(i,3,1)*ar(3)
1563 ald(2) = ad(2)+vqn(i,1,2)*dbad(1)+vqn(i,2,2)*dbad(2)
1564 1 +vqn(i,3,2)*dbad(3)
1565 2 -db(i,1,2)*ar(1)-db(i,2,2)*ar(2)-db(i,3,2)*ar(3)
1566 ald(3) = ad(3)+vqn(i,1,3)*dbad(1)+vqn(i,2,3)*dbad(2)
1567 1 +vqn(i,3,3)*dbad(3)
1568 2 -db(i,1,3)*ar(1)-db(i,2,3)*ar(2)-db(i,3,3)*ar(3)
1569 ald(4) = ad(4)+vqn(i,1,4)*dbad(1)+vqn(i,2,4)*dbad(2)
1570 1 +vqn(i,3,4)*dbad(3)
1571 2 -db(i,1,4)*ar(1)-db(i,2,4)*ar(2)-db(i,3,4)*ar(3)
1572C
1573 c1 = two*alr(3)
1574 v13(i,1)= v13(i,1)+c1*y13(i)
1575 v24(i,1)= v24(i,1)+c1*y24(i)
1576 vhi(i,1)= vhi(i,1)+four*(alr(3)*my13(i)-z1(i)*alr(2))
1577 v13(i,2)= v13(i,2)-c1*x13(i)
1578 v24(i,2)= v24(i,2)-c1*x24(i)
1579 vhi(i,2)= vhi(i,2)-four*(alr(3)*mx13(i)-z1(i)*alr(1))
1580 v13(i,3)= v13(i,3)-two*(y13(i)*alr(1)-x13(i)*alr(2))
1581 v24(i,3)= v24(i,3)-two*(y24(i)*alr(1)-x24(i)*alr(2))
1582 vhi(i,3)= vhi(i,3)+four*(mx13(i)*alr(2)-my13(i)*alr(1))
1583C
1584 rlxyzv(i,1,1)= rrxyz(1,1)-alr(1)-vqn(i,1,1)*ald(1)
1585 rlxyzv(i,1,2)= rrxyz(1,2)-alr(1)-vqn(i,1,2)*ald(2)
1586 rlxyzv(i,1,3)= rrxyz(1,3)-alr(1)-vqn(i,1,3)*ald(3)
1587 rlxyzv(i,1,4)= rrxyz(1,4)-alr(1)-vqn(i,1,4)*ald(4)
1588C
1589 rlxyzv(i,2,1)= rrxyz(2,1)-alr(2)-vqn(i,2,1)*ald(1)
1590 rlxyzv(i,2,2)= rrxyz(2,2)-alr(2)-vqn(i,2,2)*ald(2)
1591 rlxyzv(i,2,3)= rrxyz(2,3)-alr(2)-vqn(i,2,3)*ald(3)
1592 rlxyzv(i,2,4)= rrxyz(2,4)-alr(2)-vqn(i,2,4)*ald(4)
1593 IF (idril>0) THEN
1594 d(1)= yy+zz+four
1595 d(2)= xx+zz+four
1596 d(3)= xx+yy+four
1597 d(4)= -xy
1598 d(5)= -xz
1599 d(6)= -yz
1600 IF(iresp == 1)THEN
1601 CALL a3invdp(d,diz1)
1602 diz(i,3) = diz1(3)
1603 diz(i,1) = diz1(5)
1604 diz(i,2) = diz1(6)
1605 ELSE
1606 abc = d(1)*d(2)*d(3)
1607 xxyz2 = d(1)*d(6)*d(6)
1608 yyxz2 = d(2)*d(5)*d(5)
1609 zzxy2 = d(3)*d(4)*d(4)
1610 deta = abs(abc+two*d(4)*d(5)*d(6)-xxyz2-yyxz2-zzxy2)
1611 deta = one/max(deta,em20)
1612 diz(i,3) = (abc-zzxy2)*deta/max(d(3),em20)
1613 diz(i,1) = (d(6)*d(4)-d(5)*d(2))*deta
1614 diz(i,2) = (d(4)*d(5)-d(6)*d(1))*deta
1615 END IF !(IRESP == 1)THEN
1616C
1617 alrz=area_i(i)*(diz(i,1)*ar(1)+diz(i,2)*ar(2)+diz(i,3)*ar(3))
1618 rlz(i,1)=rlz(i,1)-alrz
1619 rlz(i,2)=rlz(i,2)-alrz
1620 rlz(i,3)=rlz(i,3)-alrz
1621 rlz(i,4)=rlz(i,4)-alrz
1622 END IF !IF (IDRIL>0) THEN
1623C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1624 END IF !((IMPL_S>0.AND.IKPROJ<0).OR.IDRIL>0) THEN
1625C
1626 ENDIF
1627C
1628 ENDDO
1629C
1630 RETURN
1631 END
1632!||====================================================================
1633!|| czcorp5x ../engine/source/elements/shell/coquez/czcorc.F
1634!||--- called by ------------------------------------------------------
1635!|| czcorc1 ../engine/source/elements/shell/coquez/czcorc.F
1636!||--- calls -----------------------------------------------------
1637!|| a3invdp ../engine/source/elements/shell/coquez/czcorc.F
1638!||====================================================================
1639 SUBROUTINE czcorp5x(JFT ,JLT ,VR ,NPT ,TOL ,
1640 2 IXC ,PLAT ,AREA ,AREA_I ,V13 ,
1641 3 V24 ,VHI ,RLXYZ ,VQN ,VQ ,
1642 4 X13 ,X24 ,Y13 ,Y24 ,MX13 ,
1643 6 MX23 ,MX34 ,MY13 ,MY23 ,MY34 ,
1644 7 Z1 ,DI ,DB ,COREL ,RLZ ,
1645 8 LL ,X13_2 ,Y13_2 ,X24_2 ,Y24_2 ,
1646 9 L13 ,L24 ,VRX1 ,VRX2 ,VRX3 ,
1647 A VRX4 ,VRY1 ,VRY2 ,VRY3 ,VRY4 ,
1648 B VRZ1 ,VRZ2 ,VRZ3 ,VRZ4 )
1649C-----------------------------------------------
1650C I m p l i c i t T y p e s
1651C-----------------------------------------------
1652#include "implicit_f.inc"
1653C-----------------------------------------------
1654C G l o b a l P a r a m e t e r s
1655C-----------------------------------------------
1656#include "mvsiz_p.inc"
1657#include "impl1_c.inc"
1658#include "scr05_c.inc"
1659C-----------------------------------------------
1660C D u m m y A r g u m e n t s
1661C-----------------------------------------------
1662 LOGICAL PLAT(*)
1663 INTEGER IXC(NIXC,*),JFT,JLT,IDRIL,NPT
1664 my_real
1665 . vr(3,*),v13(mvsiz,3),v24(mvsiz,3),vhi(mvsiz,3),
1666 . mx23(*),my13(*),my23(*),my34(*),
1667 . x13(*),x24(*),y13(*),y24(*),mx13(*),
1668 . vq(mvsiz,3,3),area(*),z1(*),mx34(*),vqn(mvsiz,3,4),area_i(*),
1669 . di(mvsiz,6),db(mvsiz,3,4),ll(*),l13(*),l24(*),
1670 . tol,x13_2(*),y13_2(*),x24_2(*),y24_2(*),
1671 . rlxyz(mvsiz,2,4),corel(mvsiz,2,4),rlz(mvsiz,4),
1672 . vrx1(*),vrx2(*),vrx3(*),vrx4(*),
1673 . vry1(*),vry2(*),vry3(*),vry4(*),
1674 . vrz1(*),vrz2(*),vrz3(*),vrz4(*)
1675C-----------------------------------------------
1676C L O C A L V A R I A B L E S
1677C-----------------------------------------------
1678 INTEGER NNOD,I,J,K,L
1679 PARAMETER (NNOD = 4)
1680 my_real
1681 . deta,
1682 . vcore(3,nnod),rrxyz(3,nnod),vg13(3),vg24(3),vghi(3),
1683 . z2(mvsiz),a_4,sz,sz1,sz2,sl,c1,c2,s1,
1684 . ar(3),ad(nnod),btb(6),xx,yy,zz,xy,xz,yz,
1685 . abc,xxyz2,yyxz2,zzxy2,d(6),diz2(6),
1686 . alr(3),ald(nnod),dbad(3),btb_c
1687C=======================================================================
1688 DO i=jft,jlt
1689 z2(i)=z1(i)*z1(i)
1690 IF (z2(i)<ll(i)*tol.OR.npt == 1) THEN
1691 z1(i)=zero
1692 plat(i)=.true.
1693 ELSE
1694 plat(i)=.false.
1695C--------------------------------------------------
1696C WARPING SPECIAL TREATMENT
1697C full projection eliminer drilling rotations and rigid rotations
1698C--------------------------------------------------------------------------
1699 a_4=area(i)*fourth
1700C
1701C---------- node N ----------
1702 sz1=mx13(i)*y24(i)-my13(i)*x24(i)
1703 sz2=a_4+sz1
1704 sz=z2(i)*l24(i)
1705 sl=one/sqrt(sz+sz2*sz2)
1706 vqn(i,1,1)=-z1(i)*y24(i)
1707 vqn(i,2,1)= z1(i)*x24(i)
1708 vqn(i,3,1)= sz2*sl
1709 vqn(i,1,3)=-vqn(i,1,1)
1710 vqn(i,2,3)=-vqn(i,2,1)
1711 vqn(i,1,1)= vqn(i,1,1)*sl
1712 vqn(i,2,1)= vqn(i,2,1)*sl
1713C
1714 sz2=a_4-sz1
1715 sl=one/sqrt(sz+sz2*sz2)
1716 vqn(i,1,3)= vqn(i,1,3)*sl
1717 vqn(i,2,3)= vqn(i,2,3)*sl
1718 vqn(i,3,3)= sz2*sl
1719C
1720 sz1=mx13(i)*y13(i)-my13(i)*x13(i)
1721 sz2=a_4+sz1
1722 sz=z2(i)*l13(i)
1723 sl=one/sqrt(sz+sz2*sz2)
1724 vqn(i,1,2)=-z1(i)*y13(i)
1725 vqn(i,2,2)= z1(i)*x13(i)
1726 vqn(i,3,2)= sz2*sl
1727 vqn(i,1,4)=-vqn(i,1,2)
1728 vqn(i,2,4)=-vqn(i,2,2)
1729 vqn(i,1,2)= vqn(i,1,2)*sl
1730 vqn(i,2,2)= vqn(i,2,2)*sl
1731C
1732 sz2=a_4-sz1
1733 sl=one/sqrt(sz+sz2*sz2)
1734 vqn(i,1,4)= vqn(i,1,4)*sl
1735 vqn(i,2,4)= vqn(i,2,4)*sl
1736 vqn(i,3,4)= sz2*sl
1737C
1738 k=ixc(2,i)
1739 rrxyz(1,1) =rlxyz(i,1,1)
1740 rrxyz(2,1) =rlxyz(i,2,1)
1741 rrxyz(3,1) =vq(i,1,3)*vrx1(i)+vq(i,2,3)*vry1(i)
1742 1 +vq(i,3,3)*vrz1(i)
1743 k=ixc(3,i)
1744 rrxyz(1,2) =rlxyz(i,1,2)
1745 rrxyz(2,2) =rlxyz(i,2,2)
1746 rrxyz(3,2) =vq(i,1,3)*vrx2(i)+vq(i,2,3)*vry2(i)
1747 1 +vq(i,3,3)*vrz2(i)
1748 k=ixc(4,i)
1749 rrxyz(1,3) =rlxyz(i,1,3)
1750 rrxyz(2,3) =rlxyz(i,2,3)
1751 rrxyz(3,3) =vq(i,1,3)*vrx3(i)+vq(i,2,3)*vry3(i)
1752 1 +vq(i,3,3)*vrz3(i)
1753 k=ixc(5,i)
1754 rrxyz(1,4) =rlxyz(i,1,4)
1755 rrxyz(2,4) =rlxyz(i,2,4)
1756 rrxyz(3,4) =vq(i,1,3)*vrx4(i)+vq(i,2,3)*vry4(i)
1757 1 +vq(i,3,3)*vrz4(i)
1758 IF (impl_s>0.AND.ikproj<0) THEN
1759C-------------------------------------
1760C DRILLING PROJECTION ONLY
1761C-------------------------------------
1762 rlxyz(i,1,1)=(1.-vqn(i,1,1)*vqn(i,1,1))*rrxyz(1,1)
1763 1 -vqn(i,1,1)*vqn(i,2,1) *rrxyz(2,1)
1764 2 -vqn(i,1,1)*vqn(i,3,1) *rrxyz(3,1)
1765 rlxyz(i,2,1)=(1.-vqn(i,2,1)*vqn(i,2,1))*rrxyz(2,1)
1766 1 -vqn(i,1,1)*vqn(i,2,1) *rrxyz(1,1)
1767 2 -vqn(i,2,1)*vqn(i,3,1) *rrxyz(3,1)
1768C
1769 rlxyz(i,1,2)=(1.-vqn(i,1,2)*vqn(i,1,2))*rrxyz(1,2)
1770 1 -vqn(i,1,2)*vqn(i,2,2) *rrxyz(2,2)
1771 2 -vqn(i,1,2)*vqn(i,3,2) *rrxyz(3,2)
1772 rlxyz(i,2,2)=(1.-vqn(i,2,2)*vqn(i,2,2))*rrxyz(2,2)
1773 1 -vqn(i,1,2)*vqn(i,2,2) *rrxyz(1,2)
1774 2 -vqn(i,2,2)*vqn(i,3,2) *rrxyz(3,2)
1775C
1776 rlxyz(i,1,3)=(1.-vqn(i,1,3)*vqn(i,1,3))*rrxyz(1,3)
1777 1 -vqn(i,1,3)*vqn(i,2,3) *rrxyz(2,3)
1778 2 -vqn(i,1,3)*vqn(i,3,3) *rrxyz(3,3)
1779 rlxyz(i,2,3)=(1.-vqn(i,2,3)*vqn(i,2,3))*rrxyz(2,3)
1780 1 -vqn(i,1,3)*vqn(i,2,3) *rrxyz(1,3)
1781 2 -vqn(i,2,3)*vqn(i,3,3) *rrxyz(3,3)
1782C
1783 rlxyz(i,1,4)=(1.-vqn(i,1,4)*vqn(i,1,4))*rrxyz(1,4)
1784 1 -vqn(i,1,4)*vqn(i,2,4) *rrxyz(2,4)
1785 2 -vqn(i,1,4)*vqn(i,3,4) *rrxyz(3,4)
1786 rlxyz(i,2,4)=(1.-vqn(i,2,4)*vqn(i,2,4))*rrxyz(2,4)
1787 1 -vqn(i,1,4)*vqn(i,2,4) *rrxyz(1,4)
1788 2 -vqn(i,2,4)*vqn(i,3,4) *rrxyz(3,4)
1789 ELSE
1790C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1791C-----------------------full projection------------------
1792 ar(1)=-z1(i)*vhi(i,2)+y13(i)*v13(i,3)+y24(i)*v24(i,3)
1793 1 +my13(i)*vhi(i,3)
1794 2 +rrxyz(1,1)+rrxyz(1,2)+rrxyz(1,3)+rrxyz(1,4)
1795 ar(2)= z1(i)*vhi(i,1)-x13(i)*v13(i,3)-x24(i)*v24(i,3)
1796 1 -mx13(i)*vhi(i,3)
1797 2 +rrxyz(2,1)+rrxyz(2,2)+rrxyz(2,3)+rrxyz(2,4)
1798 ar(3)= x13(i)*v13(i,2)+x24(i)*v24(i,2)+mx13(i)*vhi(i,2)
1799 1 -y13(i)*v13(i,1)-y24(i)*v24(i,1)-my13(i)*vhi(i,1)
1800 2 +rrxyz(3,1)+rrxyz(3,2)+rrxyz(3,3)+rrxyz(3,4)
1801 ad(1)= vqn(i,1,1)*rrxyz(1,1)+vqn(i,2,1)*rrxyz(2,1)+
1802 1 vqn(i,3,1)*rrxyz(3,1)
1803 ad(2)= vqn(i,1,2)*rrxyz(1,2)+vqn(i,2,2)*rrxyz(2,2)+
1804 1 vqn(i,3,2)*rrxyz(3,2)
1805 ad(3)= vqn(i,1,3)*rrxyz(1,3)+vqn(i,2,3)*rrxyz(2,3)+
1806 1 vqn(i,3,3)*rrxyz(3,3)
1807 ad(4)= vqn(i,1,4)*rrxyz(1,4)+vqn(i,2,4)*rrxyz(2,4)+
1808 1 vqn(i,3,4)*rrxyz(3,4)
1809C
1810C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1811 xx = corel(i,1,1)*corel(i,1,1)+corel(i,1,2)*corel(i,1,2)
1812 1 +corel(i,1,3)*corel(i,1,3)+corel(i,1,4)*corel(i,1,4)
1813 yy = corel(i,2,1)*corel(i,2,1)+corel(i,2,2)*corel(i,2,2)
1814 1 +corel(i,2,3)*corel(i,2,3)+corel(i,2,4)*corel(i,2,4)
1815 xy = corel(i,1,1)*corel(i,2,1)+corel(i,1,2)*corel(i,2,2)
1816 1 +corel(i,1,3)*corel(i,2,3)+corel(i,1,4)*corel(i,2,4)
1817 xz =(corel(i,1,1)-corel(i,1,2)+corel(i,1,3)-corel(i,1,4))
1818 . *z1(i)
1819 yz =(corel(i,2,1)-corel(i,2,2)+corel(i,2,3)-corel(i,2,4))
1820 . *z1(i)
1821 zz = four*z2(i)
1822C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1823 btb(1)= vqn(i,1,1)*vqn(i,1,1)+vqn(i,1,2)*vqn(i,1,2)
1824 1 +vqn(i,1,3)*vqn(i,1,3)+vqn(i,1,4)*vqn(i,1,4)
1825 btb(2)= vqn(i,2,1)*vqn(i,2,1)+vqn(i,2,2)*vqn(i,2,2)
1826 1 +vqn(i,2,3)*vqn(i,2,3)+vqn(i,2,4)*vqn(i,2,4)
1827 btb(3)= vqn(i,3,1)*vqn(i,3,1)+vqn(i,3,2)*vqn(i,3,2)
1828 1 +vqn(i,3,3)*vqn(i,3,3)+vqn(i,3,4)*vqn(i,3,4)
1829 btb(4)= vqn(i,1,1)*vqn(i,2,1)+vqn(i,1,2)*vqn(i,2,2)
1830 1 +vqn(i,1,3)*vqn(i,2,3)+vqn(i,1,4)*vqn(i,2,4)
1831 btb(5)= vqn(i,1,1)*vqn(i,3,1)+vqn(i,1,2)*vqn(i,3,2)
1832 1 +vqn(i,1,3)*vqn(i,3,3)+vqn(i,1,4)*vqn(i,3,4)
1833 btb(6)= vqn(i,2,1)*vqn(i,3,1)+vqn(i,2,2)*vqn(i,3,2)
1834 1 +vqn(i,2,3)*vqn(i,3,3)+vqn(i,2,4)*vqn(i,3,4)
1835 d(1)= yy+zz+four-btb(1)
1836 d(2)= xx+zz+four-btb(2)
1837 d(3)= xx+yy+four-btb(3)
1838 d(4)= -xy-btb(4)
1839 d(5)= -xz-btb(5)
1840 d(6)= -yz-btb(6)
1841 IF(iresp == 1)THEN
1842 CALL a3invdp(d,diz2)
1843 di(i,1:6) = diz2(1:6)
1844 ELSE
1845 abc = d(1)*d(2)*d(3)
1846 xxyz2 = d(1)*d(6)*d(6)
1847 yyxz2 = d(2)*d(5)*d(5)
1848 zzxy2 = d(3)*d(4)*d(4)
1849 deta = abs(abc+two*d(4)*d(5)*d(6)-xxyz2-yyxz2-zzxy2)
1850 deta = one/max(deta,em20)
1851 di(i,1) = (abc-xxyz2)*deta/max(d(1),em20)
1852 di(i,2) = (abc-yyxz2)*deta/max(d(2),em20)
1853 di(i,3) = (abc-zzxy2)*deta/max(d(3),em20)
1854 di(i,4) = (d(5)*d(6)-d(4)*d(3))*deta
1855 di(i,5) = (d(6)*d(4)-d(5)*d(2))*deta
1856 di(i,6) = (d(4)*d(5)-d(6)*d(1))*deta
1857 END IF !(IRESP == 1)THEN
1858 DO j=1,nnod
1859 db(i,1,j)= di(i,1)*vqn(i,1,j)+di(i,4)*vqn(i,2,j)
1860 1 +di(i,5)*vqn(i,3,j)
1861 db(i,2,j)= di(i,4)*vqn(i,1,j)+di(i,2)*vqn(i,2,j)
1862 1 +di(i,6)*vqn(i,3,j)
1863 db(i,3,j)= di(i,5)*vqn(i,1,j)+di(i,6)*vqn(i,2,j)
1864 1 +di(i,3)*vqn(i,3,j)
1865 ENDDO
1866C
1867 dbad(1)= db(i,1,1)*ad(1)+db(i,1,2)*ad(2)
1868 1 +db(i,1,3)*ad(3)+db(i,1,4)*ad(4)
1869 dbad(2)= db(i,2,1)*ad(1)+db(i,2,2)*ad(2)
1870 1 +db(i,2,3)*ad(3)+db(i,2,4)*ad(4)
1871 dbad(3)= db(i,3,1)*ad(1)+db(i,3,2)*ad(2)
1872 1 +db(i,3,3)*ad(3)+db(i,3,4)*ad(4)
1873C
1874 alr(1) =di(i,1)*ar(1)+di(i,4)*ar(2)+di(i,5)*ar(3)-dbad(1)
1875 alr(2) =di(i,4)*ar(1)+di(i,2)*ar(2)+di(i,6)*ar(3)-dbad(2)
1876 alr(3) =di(i,5)*ar(1)+di(i,6)*ar(2)+di(i,3)*ar(3)-dbad(3)
1877C
1878 ald(1) = ad(1)+vqn(i,1,1)*dbad(1)+vqn(i,2,1)*dbad(2)
1879 1 +vqn(i,3,1)*dbad(3)
1880 2 -db(i,1,1)*ar(1)-db(i,2,1)*ar(2)-db(i,3,1)*ar(3)
1881 ald(2) = ad(2)+vqn(i,1,2)*dbad(1)+vqn(i,2,2)*dbad(2)
1882 1 +vqn(i,3,2)*dbad(3)
1883 2 -db(i,1,2)*ar(1)-db(i,2,2)*ar(2)-db(i,3,2)*ar(3)
1884 ald(3) = ad(3)+vqn(i,1,3)*dbad(1)+vqn(i,2,3)*dbad(2)
1885 1 +vqn(i,3,3)*dbad(3)
1886 2 -db(i,1,3)*ar(1)-db(i,2,3)*ar(2)-db(i,3,3)*ar(3)
1887 ald(4) = ad(4)+vqn(i,1,4)*dbad(1)+vqn(i,2,4)*dbad(2)
1888 1 +vqn(i,3,4)*dbad(3)
1889 2 -db(i,1,4)*ar(1)-db(i,2,4)*ar(2)-db(i,3,4)*ar(3)
1890C
1891 c1 = two*alr(3)
1892 v13(i,1)= v13(i,1)+c1*y13(i)
1893 v24(i,1)= v24(i,1)+c1*y24(i)
1894 vhi(i,1)= vhi(i,1)+four*(alr(3)*my13(i)-z1(i)*alr(2))
1895 v13(i,2)= v13(i,2)-c1*x13(i)
1896 v24(i,2)= v24(i,2)-c1*x24(i)
1897 vhi(i,2)= vhi(i,2)-four*(alr(3)*mx13(i)-z1(i)*alr(1))
1898 v13(i,3)= v13(i,3)-two*(y13(i)*alr(1)-x13(i)*alr(2))
1899 v24(i,3)= v24(i,3)-two*(y24(i)*alr(1)-x24(i)*alr(2))
1900 vhi(i,3)= vhi(i,3)+four*(mx13(i)*alr(2)-my13(i)*alr(1))
1901 rlxyz(i,1,1)= rrxyz(1,1)-alr(1)-vqn(i,1,1)*ald(1)
1902 rlxyz(i,1,2)= rrxyz(1,2)-alr(1)-vqn(i,1,2)*ald(2)
1903 rlxyz(i,1,3)= rrxyz(1,3)-alr(1)-vqn(i,1,3)*ald(3)
1904 rlxyz(i,1,4)= rrxyz(1,4)-alr(1)-vqn(i,1,4)*ald(4)
1905C
1906 rlxyz(i,2,1)= rrxyz(2,1)-alr(2)-vqn(i,2,1)*ald(1)
1907 rlxyz(i,2,2)= rrxyz(2,2)-alr(2)-vqn(i,2,2)*ald(2)
1908 rlxyz(i,2,3)= rrxyz(2,3)-alr(2)-vqn(i,2,3)*ald(3)
1909 rlxyz(i,2,4)= rrxyz(2,4)-alr(2)-vqn(i,2,4)*ald(4)
1910C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1911 END IF !(IMPL_S>0.AND.IKPROJ<0) THEN
1912C
1913 ENDIF
1914 ENDDO
1915C
1916 RETURN
1917 END
1918!||====================================================================
1919!|| czcorct ../engine/source/elements/shell/coquez/czcorc.F
1920!||--- called by ------------------------------------------------------
1921!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
1922!||--- calls -----------------------------------------------------
1923!|| clskew3 ../engine/source/elements/sh3n/coquedk/cdkcoor3.F
1924!||--- uses -----------------------------------------------------
1925!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.f90
1926!||====================================================================
1927 SUBROUTINE czcorct(ELBUF_STR,
1928 1 JFT ,JLT ,X ,V ,VR ,
1929 2 IXC ,PM ,OFFG ,AREA ,AREA_I ,
1930 3 V13 ,V24 ,DR ,RLXYZ ,VQ ,
1931 4 X13_T ,X24_T ,Y13_T ,Y24_T ,MX13 ,
1932 5 MX23 ,MX34 ,MY13 ,MY23 ,MY34 ,
1933 6 Z1 ,SMSTR ,THK ,NPT ,ISMSTR ,
1934 7 IDRIL ,XLCOR ,ZL ,VQN ,NEL )
1935C-----------------------------------------------
1936C M o d u l e s
1937C-----------------------------------------------
1938 USE elbufdef_mod
1939C-----------------------------------------------
1940#include "implicit_f.inc"
1941c-----------------------------------------------
1942c g l o b a l p a r a m e t e r s
1943c-----------------------------------------------
1944#include "mvsiz_p.inc"
1945#include "param_c.inc"
1946#include "scr05_c.inc"
1947#include "impl1_c.inc"
1948C-----------------------------------------------
1949C D U M M Y A R G U M E N T S
1950C-----------------------------------------------
1951 INTEGER IXC(NIXC,*),JFT,JLT,IREP,NPT,ISMSTR,IDRIL,NEL
1952 my_real
1953 . PM(NPROPM,*), X(3,*), V(3,*), VR(3,*),RLXYZ(MVSIZ,4),
1954 . v13(mvsiz,2),v24(mvsiz,2),offg(*),dr(3,*),
1955 . mx23(*),my13(*),my23(*),my34(*),
1956 . x13_t(*),x24_t(*),y13_t(*),y24_t(*),mx13(*),
1957 . vq(mvsiz,3,3),area(*),z1(*),mx34(*),area_i(*),
1958 . thk(*),xlcor(mvsiz,2,4),zl(*),vqn(mvsiz,9,4)
1959 double precision
1960 . smstr(*)
1961 TYPE (ELBUF_STRUCT_) :: ELBUF_STR
1962C-----------------------------------------------
1963C L O C A L V A R I A B L E S
1964C-----------------------------------------------
1965 INTEGER NNOD,I,J,K,L,II(9),IFPROJ
1966 INTEGER IXCTMP2,IXCTMP3,IXCTMP4,IXCTMP5
1967 parameter(nnod = 4)
1968 my_real
1969 . lxyz0(3),deta,deta1(mvsiz),alrz,corel(mvsiz,2,4),
1970 . di(mvsiz,6),db(mvsiz,3,4),l13(mvsiz),l24(mvsiz),
1971 . rrxyz(3,nnod),vg13(3),vg24(3),vghi(3),
1972 . tol,z2(mvsiz),a_4,sz,sz1,sz2,sl,c1,c2,s1,
1973 . vlxyz(3,nnod),ar(3),ad(nnod),btb(6),xx,yy,zz,xy,xz,yz,
1974 . abc,xxyz2,yyxz2,zzxy2,d(6),diz1(6),diz(3),
1975 . alr(3),ald(nnod),dbad(3),btb_c,lm(mvsiz),
1976 . facdt,ll(mvsiz),axyz(mvsiz,3,nnod),
1977 . vhi(mvsiz,3),v13x,v24x,vhix,ddrz1,ddrz2
1978 my_real
1979 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),
1980 . yl3(mvsiz),yl4(mvsiz),zl1(mvsiz),off_l,
1981 . x13(mvsiz),x24(mvsiz),y13(mvsiz),y24(mvsiz),
1982 . x0g2(mvsiz),x0g3(mvsiz),x0g4(mvsiz),y0g2(mvsiz),
1983 . y0g3(mvsiz),y0g4(mvsiz),z0g2(mvsiz),z0g3(mvsiz),z0g4(mvsiz),
1984 . rx(mvsiz),ry(mvsiz),rz(mvsiz),sx(mvsiz),sy(mvsiz),
1985 . r11(mvsiz),r12(mvsiz),r13(mvsiz),r21(mvsiz),r22(mvsiz),
1986 . r23(mvsiz),r31(mvsiz),r32(mvsiz),r33(mvsiz),dirz(mvsiz,2),
1987 . ssz(mvsiz),vq0(mvsiz,3,3),det_1,
1988 . x0l2(mvsiz),x0l3(mvsiz),x0l4(mvsiz),y0l2(mvsiz),
1989 . y0l3(mvsiz),y0l4(mvsiz),vr1_12,vr1_21
1990C----------------------------------------------
1991 DO i=1,9
1992 ii(i) = nel*(i-1)
1993 ENDDO
1994C
1995 IF(iresp == 1)THEN
1996 tol=em7
1997 ELSE
1998 tol=em8
1999 END IF
2000C
2001 DO i=jft,jlt
2002 IF(abs(offg(i))==one)offg(i)=sign(two,offg(i))
2003 axyz(i,1:3,1:4)= zero
2004C
2005 IF (idril > 0 ) THEN
2006 ixctmp2=ixc(2,i)
2007 ixctmp3=ixc(3,i)
2008 ixctmp4=ixc(4,i)
2009 ixctmp5=ixc(5,i)
2010 axyz(i,1,1) = dr(1,ixctmp2)
2011 axyz(i,2,1) = dr(2,ixctmp2)
2012 axyz(i,3,1) = dr(3,ixctmp2)
2013 axyz(i,1,2) = dr(1,ixctmp3)
2014 axyz(i,2,2) = dr(2,ixctmp3)
2015 axyz(i,3,2) = dr(3,ixctmp3)
2016 axyz(i,1,3) = dr(1,ixctmp4)
2017 axyz(i,2,3) = dr(2,ixctmp4)
2018 axyz(i,3,3) = dr(3,ixctmp4)
2019 axyz(i,1,4) = dr(1,ixctmp5)
2020 axyz(i,2,4) = dr(2,ixctmp5)
2021 axyz(i,3,4) = dr(3,ixctmp5)
2022 END IF !(ISROT > 0 ) THEN
2023
2024 x0g2(i) = smstr(ii(1)+i)
2025 y0g2(i) = smstr(ii(2)+i)
2026 z0g2(i) = smstr(ii(3)+i)
2027 x0g3(i) = smstr(ii(4)+i)
2028 y0g3(i) = smstr(ii(5)+i)
2029 z0g3(i) = smstr(ii(6)+i)
2030 x0g4(i) = smstr(ii(7)+i)
2031 y0g4(i) = smstr(ii(8)+i)
2032 z0g4(i) = smstr(ii(9)+i)
2033 ENDDO
2034C-- normal in initial conf. to compute Z1
2035 DO i=jft,jlt
2036 rx(i)=x0g2(i)+x0g3(i)-x0g4(i)
2037 sx(i)=x0g3(i)+x0g4(i)-x0g2(i)
2038 ry(i)=y0g2(i)+y0g3(i)-y0g4(i)
2039 sy(i)=y0g3(i)+y0g4(i)-y0g2(i)
2040 rz(i)=z0g2(i)+z0g3(i)-z0g4(i)
2041 ssz(i)=z0g3(i)+z0g4(i)-z0g2(i)
2042 ENDDO
2043C----------------------------
2044C LOCAL SYSTEM VQ0
2045C----------------------------
2046 k = 0
2047 CALL clskew3(jft,jlt,k,
2048 . rx, ry, rz,
2049 . sx, sy, ssz,
2050 . r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
2051 DO i=jft,jlt
2052 area(i)=fourth*deta1(i)
2053 area_i(i)=one/area(i)
2054 vq0(i,1,1)=r11(i)
2055 vq0(i,2,1)=r21(i)
2056 vq0(i,3,1)=r31(i)
2057 vq0(i,1,2)=r12(i)
2058 vq0(i,2,2)=r22(i)
2059 vq0(i,3,2)=r32(i)
2060 vq0(i,1,3)=r13(i)
2061 vq0(i,2,3)=r23(i)
2062 vq0(i,3,3)=r33(i)
2063 ENDDO
2064 DO i=jft,jlt
2065 lxyz0(1)=fourth*(x0g2(i)+x0g3(i)+x0g4(i))
2066 lxyz0(2)=fourth*(y0g2(i)+y0g3(i)+y0g4(i))
2067 lxyz0(3)=fourth*(z0g2(i)+z0g3(i)+z0g4(i))
2068 z1(i)=-(vq0(i,1,3)*lxyz0(1)+vq0(i,2,3)*lxyz0(2)+vq0(i,3,3)*lxyz0(3))
2069 ENDDO
2070C----------------------------
2071C xl =R1^t x0l; R1^t=(VQ0^t*VQ)^t---extract Rz0 of R1
2072C----------------------------
2073 DO i=jft,jlt
2074 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)
2075 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)
2076 dirz(i,2) = half*(vr1_12-vr1_21)
2077 det_1 = one-dirz(i,2)*dirz(i,2)
2078 dirz(i,1) = sqrt(max(zero,det_1))
2079 ENDDO
2080 DO i=jft,jlt
2081 x0l2(i)=vq0(i,1,1)*x0g2(i)+vq0(i,2,1)*y0g2(i)+vq0(i,3,1)*z0g2(i)
2082 y0l2(i)=vq0(i,1,2)*x0g2(i)+vq0(i,2,2)*y0g2(i)+vq0(i,3,2)*z0g2(i)
2083 x0l3(i)=vq0(i,1,1)*x0g3(i)+vq0(i,2,1)*y0g3(i)+vq0(i,3,1)*z0g3(i)
2084 y0l3(i)=vq0(i,1,2)*x0g3(i)+vq0(i,2,2)*y0g3(i)+vq0(i,3,2)*z0g3(i)
2085 x0l4(i)=vq0(i,1,1)*x0g4(i)+vq0(i,2,1)*y0g4(i)+vq0(i,3,1)*z0g4(i)
2086 y0l4(i)=vq0(i,1,2)*x0g4(i)+vq0(i,2,2)*y0g4(i)+vq0(i,3,2)*z0g4(i)
2087 ENDDO
2088C----------------------------
2089C Rotate x0l of Rz0
2090C----------------------------
2091 DO i=jft,jlt
2092 xl2(i)= x0l2(i)*dirz(i,1)-y0l2(i)*dirz(i,2)
2093 yl2(i)= x0l2(i)*dirz(i,2)+y0l2(i)*dirz(i,1)
2094 xl3(i)= x0l3(i)*dirz(i,1)-y0l3(i)*dirz(i,2)
2095 yl3(i)= x0l3(i)*dirz(i,2)+y0l3(i)*dirz(i,1)
2096 xl4(i)= x0l4(i)*dirz(i,1)-y0l4(i)*dirz(i,2)
2097 yl4(i)= x0l4(i)*dirz(i,2)+y0l4(i)*dirz(i,1)
2098 x0l2(i)=xl2(i)
2099 y0l2(i)=yl2(i)
2100 x0l3(i)=xl3(i)
2101 y0l3(i)=yl3(i)
2102 x0l4(i)=xl4(i)
2103 y0l4(i)=yl4(i)
2104 ENDDO
2105 DO i=jft,jlt
2106 xl2(i)=xlcor(i,1,2)-xlcor(i,1,1)
2107 yl2(i)=xlcor(i,2,2)-xlcor(i,2,1)
2108 xl3(i)=xlcor(i,1,3)-xlcor(i,1,1)
2109 yl3(i)=xlcor(i,2,3)-xlcor(i,2,1)
2110 xl4(i)=xlcor(i,1,4)-xlcor(i,1,1)
2111 yl4(i)=xlcor(i,2,4)-xlcor(i,2,1)
2112 ENDDO
2113C-----------UL : xl - x0l' (rotated of rz0)
2114 DO i=jft,jlt
2115 v13(i,1)=-xl3(i)-(-x0l3(i))
2116 v24(i,1)=xl2(i)-xl4(i)-(x0l2(i)-x0l4(i))
2117 vhi(i,1)=-xl2(i)+xl3(i)-xl4(i)-(-x0l2(i)+x0l3(i)-x0l4(i))
2118 v13(i,2)=-yl3(i)-(-y0l3(i))
2119 v24(i,2)=yl2(i)-yl4(i)-(y0l2(i)-y0l4(i))
2120 vhi(i,2)=-yl2(i)+yl3(i)-yl4(i)-(-y0l2(i)+y0l3(i)-y0l4(i))
2121 vhi(i,3)=four*(zl(i)-z1(i))
2122 ENDDO
2123C-------set up B_l by rotating B0_l (Rz0)
2124 DO i=jft,jlt
2125 xl2(i)=x0l2(i)
2126 yl2(i)=y0l2(i)
2127 xl3(i)=x0l3(i)
2128 yl3(i)=y0l3(i)
2129 xl4(i)=x0l4(i)
2130 yl4(i)=y0l4(i)
2131 ENDDO
2132C----------------------------
2133 DO i=jft,jlt
2134 lxyz0(1)=fourth*(xl2(i)+xl3(i)+xl4(i))
2135 lxyz0(2)=fourth*(yl2(i)+yl3(i)+yl4(i))
2136 corel(i,1,1)=-lxyz0(1)
2137 corel(i,1,2)=xl2(i)-lxyz0(1)
2138 corel(i,1,3)=xl3(i)-lxyz0(1)
2139 corel(i,1,4)=xl4(i)-lxyz0(1)
2140 corel(i,2,1)=-lxyz0(2)
2141 corel(i,2,2)=yl2(i)-lxyz0(2)
2142 corel(i,2,3)=yl3(i)-lxyz0(2)
2143 corel(i,2,4)=yl4(i)-lxyz0(2)
2144 x13(i) =(corel(i,1,1)-corel(i,1,3))*half
2145 x24(i) =(corel(i,1,2)-corel(i,1,4))*half
2146 y13(i) =(corel(i,2,1)-corel(i,2,3))*half
2147 y24(i) =(corel(i,2,2)-corel(i,2,4))*half
2148 x13_t(i) =x13(i)*area_i(i)
2149 x24_t(i) =x24(i)*area_i(i)
2150 y13_t(i) =y13(i)*area_i(i)
2151 y24_t(i) =y24(i)*area_i(i)
2152C
2153 mx13(i)=(corel(i,1,1)+corel(i,1,3))*half
2154 mx23(i)=(corel(i,1,2)+corel(i,1,3))*half
2155 mx34(i)=(corel(i,1,3)+corel(i,1,4))*half
2156 my13(i)=(corel(i,2,1)+corel(i,2,3))*half
2157 my23(i)=(corel(i,2,2)+corel(i,2,3))*half
2158 my34(i)=(corel(i,2,3)+corel(i,2,4))*half
2159C-
2160 l13(i) = x13(i)*x13(i)+y13(i)*y13(i)
2161C-------------------------------
2162 l24(i) = x24(i)*x24(i)+y24(i)*y24(i)
2163 lm(i)=half*(l13(i)+l24(i))
2164C
2165 ENDDO
2166C----------------------------
2167C-------cas thk >>L----
2168 IF (impl_s>0) THEN
2169 DO i=jft,jlt
2170 s1=em01*thk(i)*thk(i)
2171 lm(i)=max(lm(i),s1)
2172 ENDDO
2173 END IF
2174C--------------------------
2175 IF (idril>0) THEN
2176 DO i=jft,jlt
2177 rlxyz(i,1) =vq(i,1,3)*axyz(i,1,1)+vq(i,2,3)*axyz(i,2,1)
2178 1 +vq(i,3,3)*axyz(i,3,1)
2179 rlxyz(i,2) =vq(i,1,3)*axyz(i,1,2)+vq(i,2,3)*axyz(i,2,2)
2180 1 +vq(i,3,3)*axyz(i,3,2)
2181 rlxyz(i,3) =vq(i,1,3)*axyz(i,1,3)+vq(i,2,3)*axyz(i,2,3)
2182 1 +vq(i,3,3)*axyz(i,3,3)
2183 rlxyz(i,4) =vq(i,1,3)*axyz(i,1,4)+vq(i,2,3)*axyz(i,2,4)
2184 1 +vq(i,3,3)*axyz(i,3,4)
2185 ENDDO
2186 ELSE
2187 rlxyz(jft:jlt,1:4) = zero
2188 END IF
2189C--------------------------
2190 DO i=jft,jlt
2191 z2(i)=z1(i)*z1(i)
2192 IF (z2(i)<lm(i)*tol.OR.npt == 1) THEN
2193 z1(i)=zero
2194 ELSE
2195C--------------------------------------------------
2196C WARPING SPECIAL TREATMENT
2197C full projection eliminer drilling rotations and rigid rotations
2198C--------------------------------------------------------------------------
2199 a_4=area(i)*fourth
2200C
2201C-------------------------------------
2202C DRILLING PROJECTION ONLY
2203C-------------------------------------
2204 IF (idril>0.0 ) THEN
2205 rrxyz(1,1) =vq(i,1,1)*axyz(i,1,1)+vq(i,2,1)*axyz(i,2,1)
2206 + +vq(i,3,1)*axyz(i,3,1)
2207 rrxyz(1,2) =vq(i,1,1)*axyz(i,1,2)+vq(i,2,1)*axyz(i,2,2)
2208 + +vq(i,3,1)*axyz(i,3,2)
2209 rrxyz(1,3) =vq(i,1,1)*axyz(i,1,3)+vq(i,2,1)*axyz(i,2,3)
2210 + +vq(i,3,1)*axyz(i,3,3)
2211 rrxyz(1,4) =vq(i,1,1)*axyz(i,1,4)+vq(i,2,1)*axyz(i,2,4)
2212 + +vq(i,3,1)*axyz(i,3,4)
2213 rrxyz(2,1) =vq(i,1,2)*axyz(i,1,1)+vq(i,2,2)*axyz(i,2,1)
2214 + +vq(i,3,2)*axyz(i,3,1)
2215 rrxyz(2,2) =vq(i,1,2)*axyz(i,1,2)+vq(i,2,2)*axyz(i,2,2)
2216 + +vq(i,3,2)*axyz(i,3,2)
2217 rrxyz(2,3) =vq(i,1,2)*axyz(i,1,3)+vq(i,2,2)*axyz(i,2,3)
2218 + +vq(i,3,2)*axyz(i,3,3)
2219 rrxyz(2,4) =vq(i,1,2)*axyz(i,1,4)+vq(i,2,2)*axyz(i,2,4)
2220 + +vq(i,3,2)*axyz(i,3,4)
2221 rlxyz(i,1)=vqn(i,1,1)*rrxyz(1,1)+
2222 1 vqn(i,2,1)*rrxyz(2,1)+vqn(i,3,1)*rrxyz(3,1)
2223 rlxyz(i,2)=vqn(i,1,2)*rrxyz(1,2)+
2224 1 vqn(i,2,2)*rrxyz(2,2)+vqn(i,3,2)*rrxyz(3,2)
2225 rlxyz(i,3)=vqn(i,1,3)*rrxyz(1,3)+
2226 1 vqn(i,2,3)*rrxyz(2,3)+vqn(i,3,3)*rrxyz(3,3)
2227 rlxyz(i,4)=vqn(i,1,4)*rrxyz(1,4)+
2228 1 vqn(i,2,4)*rrxyz(2,4)+vqn(i,3,4)*rrxyz(3,4)
2229 END IF
2230C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2231C
2232 END IF
2233 ENDDO
2234 off_l = zero
2235 DO i=jft,jlt
2236 off_l = min(off_l,offg(i))
2237 ENDDO
2238C
2239 IF(off_l<zero)THEN
2240 DO i=jft,jlt
2241 IF(offg(i)<zero)THEN
2242 v13(i,1:2)= zero
2243 v24(i,1:2)= zero
2244 rlxyz(i,1:4)= zero
2245 ENDIF
2246 ENDDO
2247 ENDIF
2248C
2249 RETURN
2250 END
2251!||====================================================================
2252!|| czcorcht ../engine/source/elements/shell/coquez/czcorc.F
2253!||--- called by ------------------------------------------------------
2254!|| czforc3 ../engine/source/elements/shell/coquez/czforc3.F
2255!||--- calls -----------------------------------------------------
2256!|| clskew3 ../engine/source/elements/sh3n/coquedk/cdkcoor3.F
2257!||--- uses -----------------------------------------------------
2258!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
2259!||====================================================================
2260 SUBROUTINE czcorcht(ELBUF_STR,
2261 1 JFT ,JLT ,X ,V ,VR ,
2262 2 IXC ,PM ,OFFG ,
2263 3 VQ ,HOURG ,THK ,NPT ,ISMSTR ,
2264 7 XLCOR ,ZL2 ,IINT ,NEL )
2265C-----------------------------------------------
2266C M o d u l e s
2267C-----------------------------------------------
2268 USE elbufdef_mod
2269C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
2270#include "implicit_f.inc"
2271c-----------------------------------------------
2272c g l o b a l p a r a m e t e r s
2273c-----------------------------------------------
2274#include "mvsiz_p.inc"
2275#include "param_c.inc"
2276#include "scr05_c.inc"
2277#include "com08_c.inc"
2278C-----------------------------------------------
2279C D U M M Y A R G U M E N T S
2280C-----------------------------------------------
2281 INTEGER IXC(NIXC,*),JFT,JLT,NPT,ISMSTR,NEL,IINT
2282 my_real
2283 . pm(npropm,*), x(3,*), v(3,*), vr(3,*),
2284 . offg(*),vq(mvsiz,3,3),thk(*),xlcor(mvsiz,2,4),zl2(*)
2285 my_real
2286 . hourg(nel,12)
2287 TYPE (ELBUF_STRUCT_) :: ELBUF_STR
2288C-----------------------------------------------
2289C L O C A L V A R I A B L E S
2290C-----------------------------------------------
2291 INTEGER NNOD,I,J,K,L,II(9),IFPROJ
2292 my_real
2293 . lxyz0(3),deta,deta1(mvsiz),alrz,corel(mvsiz,2,4),
2294 . l13(mvsiz),l24(mvsiz),
2295 . vg13(3),vg24(3),vghi(3),
2296 . tol,z2(mvsiz),a_2,sz,sz1,sz2,sl,c1,c2,s1,
2297 . xx,yy,zz,xy,xz,yz,lm(mvsiz),
2298 . my13(mvsiz),mx13(mvsiz),
2299 . v13(mvsiz,2),v24(mvsiz,2),area(mvsiz),area_i(mvsiz),
2300 . vhi(mvsiz,3),v13x,v24x,vhix,bxv2,byv1,tolh
2301 my_real
2302 . xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),
2303 . yl3(mvsiz),yl4(mvsiz),z1(mvsiz),off_l,
2304 . x13(mvsiz),x24(mvsiz),y13(mvsiz),y24(mvsiz),
2305 . x0g2(mvsiz),x0g3(mvsiz),x0g4(mvsiz),y0g2(mvsiz),
2306 . y0g3(mvsiz),y0g4(mvsiz),z0g2(mvsiz),z0g3(mvsiz),z0g4(mvsiz),
2307 . rx(mvsiz),ry(mvsiz),rz(mvsiz),sx(mvsiz),sy(mvsiz),
2308 . r11(mvsiz),r12(mvsiz),r13(mvsiz),r21(mvsiz),r22(mvsiz),
2309 . r23(mvsiz),r31(mvsiz),r32(mvsiz),r33(mvsiz),dirz(mvsiz,2),
2310 . ssz(mvsiz),vq0(mvsiz,3,3),det_1,
2311 . x0l2(mvsiz),x0l3(mvsiz),x0l4(mvsiz),y0l2(mvsiz),
2312 . y0l3(mvsiz),y0l4(mvsiz),vr1_12,vr1_21,hg1,hg2,vdef1,vdef2
2313C----------------------------------------------
2314 IF(iresp == 1)THEN
2315 tol=em7
2316 tolh= two*em6
2317 ELSE
2318 tol=em8
2319 tolh= em12
2320 END IF
2321C
2322 IF (iint==0.AND.tt==zero) THEN
2323 DO i=jft,jlt
2324 hourg(i,1) = x(1,ixc(3,i))-x(1,ixc(2,i))
2325 hourg(i,2) = x(2,ixc(3,i))-x(2,ixc(2,i))
2326 hourg(i,3) = x(3,ixc(3,i))-x(3,ixc(2,i))
2327 hourg(i,4) = x(1,ixc(4,i))-x(1,ixc(2,i))
2328 hourg(i,5) = x(2,ixc(4,i))-x(2,ixc(2,i))
2329 hourg(i,6) = x(3,ixc(4,i))-x(3,ixc(2,i))
2330 hourg(i,7) = x(1,ixc(5,i))-x(1,ixc(2,i))
2331 hourg(i,8) = x(2,ixc(5,i))-x(2,ixc(2,i))
2332 hourg(i,9) = x(3,ixc(5,i))-x(3,ixc(2,i))
2333 ENDDO
2334 END IF
2335C
2336 DO i=jft,jlt
2337 x0g2(i) = hourg(i,1)
2338 y0g2(i) = hourg(i,2)
2339 z0g2(i) = hourg(i,3)
2340 x0g3(i) = hourg(i,4)
2341 y0g3(i) = hourg(i,5)
2342 z0g3(i) = hourg(i,6)
2343 x0g4(i) = hourg(i,7)
2344 y0g4(i) = hourg(i,8)
2345 z0g4(i) = hourg(i,9)
2346 ENDDO
2347C-- normal in initial conf. to compute Z1
2348 DO i=jft,jlt
2349 rx(i)=x0g2(i)+x0g3(i)-x0g4(i)
2350 sx(i)=x0g3(i)+x0g4(i)-x0g2(i)
2351 ry(i)=y0g2(i)+y0g3(i)-y0g4(i)
2352 sy(i)=y0g3(i)+y0g4(i)-y0g2(i)
2353 rz(i)=z0g2(i)+z0g3(i)-z0g4(i)
2354 ssz(i)=z0g3(i)+z0g4(i)-z0g2(i)
2355 ENDDO
2356C----------------------------
2357C LOCAL SYSTEM VQ0
2358C----------------------------
2359 k = 0
2360 CALL clskew3(jft,jlt,k,
2361 . rx, ry, rz,
2362 . sx, sy, ssz,
2363 . r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
2364 DO i=jft,jlt
2365 area(i)=fourth*deta1(i)
2366 area_i(i)=one/area(i)
2367 vq0(i,1,1)=r11(i)
2368 vq0(i,2,1)=r21(i)
2369 vq0(i,3,1)=r31(i)
2370 vq0(i,1,2)=r12(i)
2371 vq0(i,2,2)=r22(i)
2372 vq0(i,3,2)=r32(i)
2373 vq0(i,1,3)=r13(i)
2374 vq0(i,2,3)=r23(i)
2375 vq0(i,3,3)=r33(i)
2376 ENDDO
2377 DO i=jft,jlt
2378 lxyz0(1)=fourth*(x0g2(i)+x0g3(i)+x0g4(i))
2379 lxyz0(2)=fourth*(y0g2(i)+y0g3(i)+y0g4(i))
2380 lxyz0(3)=fourth*(z0g2(i)+z0g3(i)+z0g4(i))
2381 z1(i)=-(vq0(i,1,3)*lxyz0(1)+vq0(i,2,3)*lxyz0(2)+vq0(i,3,3)*lxyz0(3))
2382 ENDDO
2383C----------------------------
2384C xl =R1^t x0l; R1^t=(VQ0^t*VQ)^t---extract Rz0 of R1
2385C----------------------------
2386 DO i=jft,jlt
2387 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)
2388 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)
2389 dirz(i,2) = half*(vr1_12-vr1_21)
2390 det_1 = one-dirz(i,2)*dirz(i,2)
2391 dirz(i,1) = sqrt(max(zero,det_1))
2392 ENDDO
2393 DO i=jft,jlt
2394 x0l2(i)=vq0(i,1,1)*x0g2(i)+vq0(i,2,1)*y0g2(i)+vq0(i,3,1)*z0g2(i)
2395 y0l2(i)=vq0(i,1,2)*x0g2(i)+vq0(i,2,2)*y0g2(i)+vq0(i,3,2)*z0g2(i)
2396 x0l3(i)=vq0(i,1,1)*x0g3(i)+vq0(i,2,1)*y0g3(i)+vq0(i,3,1)*z0g3(i)
2397 y0l3(i)=vq0(i,1,2)*x0g3(i)+vq0(i,2,2)*y0g3(i)+vq0(i,3,2)*z0g3(i)
2398 x0l4(i)=vq0(i,1,1)*x0g4(i)+vq0(i,2,1)*y0g4(i)+vq0(i,3,1)*z0g4(i)
2399 y0l4(i)=vq0(i,1,2)*x0g4(i)+vq0(i,2,2)*y0g4(i)+vq0(i,3,2)*z0g4(i)
2400 ENDDO
2401C----------------------------
2402C Rotate x0l of Rz0
2403C----------------------------
2404 DO i=jft,jlt
2405 xl2(i)= x0l2(i)*dirz(i,1)-y0l2(i)*dirz(i,2)
2406 yl2(i)= x0l2(i)*dirz(i,2)+y0l2(i)*dirz(i,1)
2407 xl3(i)= x0l3(i)*dirz(i,1)-y0l3(i)*dirz(i,2)
2408 yl3(i)= x0l3(i)*dirz(i,2)+y0l3(i)*dirz(i,1)
2409 xl4(i)= x0l4(i)*dirz(i,1)-y0l4(i)*dirz(i,2)
2410 yl4(i)= x0l4(i)*dirz(i,2)+y0l4(i)*dirz(i,1)
2411 x0l2(i)=xl2(i)
2412 y0l2(i)=yl2(i)
2413 x0l3(i)=xl3(i)
2414 y0l3(i)=yl3(i)
2415 x0l4(i)=xl4(i)
2416 y0l4(i)=yl4(i)
2417 ENDDO
2418 DO i=jft,jlt
2419 xl2(i)=xlcor(i,1,2)-xlcor(i,1,1)
2420 yl2(i)=xlcor(i,2,2)-xlcor(i,2,1)
2421 xl3(i)=xlcor(i,1,3)-xlcor(i,1,1)
2422 yl3(i)=xlcor(i,2,3)-xlcor(i,2,1)
2423 xl4(i)=xlcor(i,1,4)-xlcor(i,1,1)
2424 yl4(i)=xlcor(i,2,4)-xlcor(i,2,1)
2425 ENDDO
2426C-----------UL : xl - x0l' (rotated of rz0)
2427 DO i=jft,jlt
2428 v13(i,1)=-xl3(i)-(-x0l3(i))
2429 v24(i,1)=xl2(i)-xl4(i)-(x0l2(i)-x0l4(i))
2430 vhi(i,1)=-xl2(i)+xl3(i)-xl4(i)-(-x0l2(i)+x0l3(i)-x0l4(i))
2431 v13(i,2)=-yl3(i)-(-y0l3(i))
2432 v24(i,2)=yl2(i)-yl4(i)-(y0l2(i)-y0l4(i))
2433 vhi(i,2)=-yl2(i)+yl3(i)-yl4(i)-(-y0l2(i)+y0l3(i)-y0l4(i))
2434C----zl -zl zl -zl w 1/4
2435 vhi(i,3)=sqrt(zl2(i))-z1(i)
2436 ENDDO
2437C-------set up B_l by rotating B0_l (Rz0)
2438 DO i=jft,jlt
2439 xl2(i)=x0l2(i)
2440 yl2(i)=y0l2(i)
2441 xl3(i)=x0l3(i)
2442 yl3(i)=y0l3(i)
2443 xl4(i)=x0l4(i)
2444 yl4(i)=y0l4(i)
2445 ENDDO
2446C----------------------------
2447 DO i=jft,jlt
2448 lxyz0(1)=fourth*(xl2(i)+xl3(i)+xl4(i))
2449 lxyz0(2)=fourth*(yl2(i)+yl3(i)+yl4(i))
2450 corel(i,1,1)=-lxyz0(1)
2451 corel(i,1,2)=xl2(i)-lxyz0(1)
2452 corel(i,1,3)=xl3(i)-lxyz0(1)
2453 corel(i,1,4)=xl4(i)-lxyz0(1)
2454 corel(i,2,1)=-lxyz0(2)
2455 corel(i,2,2)=yl2(i)-lxyz0(2)
2456 corel(i,2,3)=yl3(i)-lxyz0(2)
2457 corel(i,2,4)=yl4(i)-lxyz0(2)
2458 x13(i) =(corel(i,1,1)-corel(i,1,3))*half
2459 x24(i) =(corel(i,1,2)-corel(i,1,4))*half
2460 y13(i) =(corel(i,2,1)-corel(i,2,3))*half
2461 y24(i) =(corel(i,2,2)-corel(i,2,4))*half
2462C
2463 mx13(i)=(corel(i,1,1)+corel(i,1,3))*half
2464 my13(i)=(corel(i,2,1)+corel(i,2,3))*half
2465C----1/4*hx
2466c HX(I) = FOURTH*(COREL(I,1,1)-COREL(I,1,2)+COREL(I,1,3)-COREL(I,1,4))
2467c HY(I) = FOURTH*(COREL(I,2,1)-COREL(I,2,2)+COREL(I,2,3)-COREL(I,2,4))
2468C-
2469 l13(i) = x13(i)*x13(i)+y13(i)*y13(i)
2470C-------------------------------
2471 l24(i) = x24(i)*x24(i)+y24(i)*y24(i)
2472 lm(i)=half*(l13(i)+l24(i))
2473C
2474 ENDDO
2475 DO i=jft,jlt
2476 x13(i) =x13(i)*area_i(i)
2477 x24(i) =x24(i)*area_i(i)
2478 y13(i) =y13(i)*area_i(i)
2479 y24(i) =y24(i)*area_i(i)
2480 ENDDO
2481C--------------------------
2482 DO i=jft,jlt
2483 z2(i)=z1(i)*z1(i)
2484 IF (z2(i)<lm(i)*tol.OR.npt == 1) THEN
2485 z1(i)=zero
2486 END IF
2487 ENDDO
2488 DO i=jft,jlt
2489 vdef1= y24(i)*v13(i,1)-y13(i)*v24(i,1)
2490 vdef2=-x24(i)*v13(i,2)+x13(i)*v24(i,2)
2491 bxv2= y24(i)*v13(i,2)-y13(i)*v24(i,2)
2492 byv1=-x24(i)*v13(i,1)+x13(i)*v24(i,1)
2493C----------------------------------
2494C HOURGLASS STRAIN RATE
2495C----------------------------------
2496 hg1=vhi(i,1)-mx13(i)*vdef1-my13(i)*byv1
2497 hg2=vhi(i,2)-mx13(i)*bxv2 -my13(i)*vdef2
2498 hourg(i,11) = hg1
2499 hourg(i,12) = hg2
2500 IF (abs(hg1)<tolh) hourg(i,11) = zero
2501 IF (abs(hg2)<tolh) hourg(i,12) = zero
2502 ENDDO
2503C
2504 RETURN
2505 END
2506
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
czcorc1
subroutine czcorc1(numnod, numelc, elbuf_str, jft, jlt, x, v, vr, ixc, pm, plat, area, area_i, v13, v24, vhi, rlxyz, vqn, vq, ll, l13, l24, x13, x24, y13, y24, mx13, mx23, mx34, my13, my23, my34, z1, corel, di, db, smstr, irep, npt, nlay, ismstr, dir_a, dir_b, offg, rlxyzv, corelv, facn, py1, px2, py2, r11, r12, r13, r21, r22, r23, r31, r32, r33, rlz, idril, ixfem, vx1, vx2, vx3, vx4, vy1, vy2, vy3, vy4, vz1, vz2, vz3, vz4, vrx1, vrx2, vrx3, vrx4, vry1, vry2, vry3, vry4, vrz1, vrz2, vrz3, vrz4, x1g, x2g, x3g, x4g, y1g, y2g, y3g, y4g, z1g, z2g, z3g, z4g, thk, diz, ux1, ux2, ux3, ux4, uy1, uy2, uy3, uy4, xl2, xl3, xl4, yl2, yl3, yl4, vl1, vl2, vl3, vl4, nel, z2)
Definition czcorc.F:61
czcorp4v
subroutine czcorp4v(jft, jlt, vr, npt, tol, ixc, plat, area, area_i, v13, v24, vhi, rlxyzv, vqn, vq, x13, x24, y13, y24, mx13, mx23, mx34, my13, my23, my34, z1, di, db, corelv, rlz, ll, x13_2, y13_2, x24_2, y24_2, l13, l24)
Definition czcorc.F:659
czcorp5x
subroutine czcorp5x(jft, jlt, vr, npt, tol, ixc, plat, area, area_i, v13, v24, vhi, rlxyz, vqn, vq, x13, x24, y13, y24, mx13, mx23, mx34, my13, my23, my34, z1, di, db, corel, rlz, ll, x13_2, y13_2, x24_2, y24_2, l13, l24, vrx1, vrx2, vrx3, vrx4, vry1, vry2, vry3, vry4, vrz1, vrz2, vrz3, vrz4)
Definition czcorc.F:1649
czcorp5v
subroutine czcorp5v(jft, jlt, vr, npt, tol, ixc, plat, area, area_i, v13, v24, vhi, rlxyzv, vqn, vq, x13, x24, y13, y24, mx13, mx23, mx34, my13, my23, my34, z1, di, db, corelv, rlz, ll, x13_2, y13_2, x24_2, y24_2, l13, l24, idril, diz)
Definition czcorc.F:1323
czcorcht
subroutine czcorcht(elbuf_str, jft, jlt, x, v, vr, ixc, pm, offg, vq, hourg, thk, npt, ismstr, xlcor, zl2, iint, nel)
Definition czcorc.F:2265
czcorp4
subroutine czcorp4(jft, jlt, vr, npt, tol, ixc, plat, area, area_i, v13, v24, vhi, rlxyz, vqn, vq, x13, x24, y13, y24, mx13, mx23, mx34, my13, my23, my34, z1, di, db, corel, rlz, ll, x13_2, y13_2, x24_2, y24_2, l13, l24)
Definition czcorc.F:870
czcorp4x
subroutine czcorp4x(jft, jlt, vr, npt, tol, ixc, plat, area, area_i, v13, v24, vhi, rlxyz, vqn, vq, x13, x24, y13, y24, mx13, mx23, mx34, my13, my23, my34, z1, di, db, corel, rlz, ll, x13_2, y13_2, x24_2, y24_2, l13, l24, vrx1, vrx2, vrx3, vrx4, vry1, vry2, vry3, vry4, vrz1, vrz2, vrz3, vrz4)
Definition czcorc.F:1094
czcorct
subroutine czcorct(elbuf_str, jft, jlt, x, v, vr, ixc, pm, offg, area, area_i, v13, v24, dr, rlxyz, vq, x13_t, x24_t, y13_t, y24_t, mx13, mx23, mx34, my13, my23, my34, z1, smstr, thk, npt, ismstr, idril, xlcor, zl, vqn, nel)
Definition czcorc.F:1935
a3invdp
subroutine a3invdp(d, di)
Definition czcorc.F:619
czcorp5
subroutine czcorp5(numnod, nel, numelc, vr, npt, tol, ixc, plat, area, area_i, v13, v24, vhi, rlxyz, vqn, vq, x13, x24, y13, y24, mx13, my13, z1, di, db, corel, rlz, ll, l13, l24, idril, diz)
Definition czcorp5.F:38
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