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