OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cbacoork.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!|| cbacoork ../engine/source/elements/shell/coqueba/cbacoork.F
25!||--- called by ------------------------------------------------------
26!|| cbake3 ../engine/source/elements/shell/coqueba/cbake3.F
27!||--- calls -----------------------------------------------------
28!|| clskew3 ../engine/source/elements/sh3n/coquedk/cdkcoor3.F
29!|| cortdir3 ../engine/source/elements/shell/coque/cortdir3.F
30!||--- uses -----------------------------------------------------
31!|| elbufdef_mod ../common_source/modules/mat_elem/elbufdef_mod.F90
32!|| element_mod ../common_source/modules/elements/element_mod.F90
33!||====================================================================
34 SUBROUTINE cbacoork(JFT,JLT,X,IXC,PM,OFFG,
35 1 GEO,AREA,VCORE,JAC,HX,HY,
36 2 VQN,VQG,VQ,VJFI,VNRM,VASTN,NPLAT,IPLAT,
37 3 X13_T ,X24_T ,Y13_T,Y24_T,
38 4 ELBUF_STR,NLAY, SMSTR,
39 5 IREP,NPT,ISMSTR,DIR_A ,DIR_B,
40 6 PID,MAT,NGL,OFF,ISROT ,NEL)
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE elbufdef_mod
45 use element_mod , only : nixc
46C-----------------------------------------------
47C I M P L I C I T T Y P E S
48C-----------------------------------------------
49#include "implicit_f.inc"
50c-----------------------------------------------
51c g l o b a l p a r a m e t e r s
52c-----------------------------------------------
53#include "mvsiz_p.inc"
54#include "param_c.inc"
55#include "impl1_c.inc"
56#include "comlock.inc"
57#include "units_c.inc"
58#include "scr17_c.inc"
59C-----------------------------------------------
60C D U M M Y A R G U M E N T S
61C-----------------------------------------------
62 INTEGER IXC(NIXC,*),JFT,JLT,NNOD,NLAY,NPLAT,IPLAT(*),ISROT,NEL
63 MY_REAL
64 . X(3,*), PM(NPROPM,*),GEO(NPROPG,*),OFFG(*)
65 PARAMETER (NNOD = 4)
66 my_real
67 . vcore(mvsiz,3,nnod),area(*),vjfi(mvsiz,6,4),
68 . vqn(mvsiz,9,nnod),vqg(mvsiz,9,nnod),vq(mvsiz,3,3),
69 . vastn(mvsiz,4,nnod),vnrm(mvsiz,3,nnod),
70 . jac(mvsiz,4),hx(mvsiz,4),hy(mvsiz,4),y24_t(*),
71 . dir_a(nel,*),dir_b(nel,*),off(*),
72 . x13_t(*),x24_t(*),y13_t(*)
73 double precision
74 . smstr(*)
75 INTEGER IREP,NPT,ISMSTR,PID(*),MAT(*),NGL(*)
76 TYPE(ELBUF_STRUCT_) :: ELBUF_STR
77C-----------------------------------------------
78c FUNCTION: preliminary compute for QBAT [K] build
79c
80c Note:
81c ARGUMENTS: (I: input, O: output, IO: input * output, W: workspace)
82c
83c TYPE NAME FUNCTION
84c I JFT,JLT - element id limit
85c I X - coordinate x,y,z in global system (basic)
86c I IXC(NIXC,*NEL) - element connectivity of shell (and other data)
87c I PM(NPROPM,MID),GEO - input Material and geometric property data
88c I OFFG(NEL),OFF(NEL) - element activation flag, local flag
89c O AREA(NEL) ,NPT - Area (element), num. of integrating points in thickness
90c O VCORE(3,4,NEL) - coordinates of 4 nodes in local system
91c BX0(2),BY0(2),GAMA(2),MX23,MY23,MX34,MY34,MX13,MY13
92c for plat element
93c O JAC,HX,HY(4,NEL) - Jacobian,asymmetric part(hourglass) at 4 Gauss points
94c O VQ(9,NEL),VQN,VQG(9,4,NEL) - local system of element, nodal(4) and Gauss points(4)
95c O VJFI,VNRM,VASTN - terms used for assumed strains
96c O NPLAT,IPLAT(NEL) - num. of plat element and indice
97c O Xij_T,Yij_T(NEL) - [B] components used for in-plane shear assumed strain
98c IO XiS,YiS,ZiS(NEL) - reference configurations used for small strain options
99c I ISMSTR,IREP - small strain and orthotrop flags
100c O DIR_A,DIR_B(2,NEL) - orthotropic directions
101c I ISROT - drilling dof flag
102c O PID,MAT,NGL(NEL) - geometric property id, material id and users' element id
103C-----------------------------------------------
104C L O C A L V A R I A B L E S
105C-----------------------------------------------
106 INTEGER J,I,II(6),K,EP,SPLAT,JPLAT(MVSIZ),L,M,MAT_1
107 MY_REAL
108 . J0,J1,J2,PG
109 MY_REAL
110 . X13,X24,Y13,MX13,MX23,MX34,MY13,Z1,Z2,GAMA1,GAMA2,
111 . X21,X34,Y21,Y34,Z21,Z34,X41,X32,Y41,Y32,Z41,Z32,L12,L34,
112 . a_4,sl,sz2,sz,jmx13,jmy13,jmz13,j2myz,y24,
113 . my23,my34,g1x1,g1x3,g1y1,g1y3,g2x1,g2x2,g2y1,g2y2
114 my_real
115 . lxyz0(3),deta1(mvsiz),rx(mvsiz), ry(mvsiz), rz(mvsiz),
116 . sx(mvsiz),sy(mvsiz),r11(mvsiz),r12(mvsiz),r13(mvsiz),
117 . r21(mvsiz),r22(mvsiz),r23(mvsiz),r31(mvsiz),r32(mvsiz),
118 . r33(mvsiz),xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),
119 . yl3(mvsiz),yl4(mvsiz),zl1(mvsiz),area_i(mvsiz),ssz(mvsiz),
120 . l13(mvsiz),l24(mvsiz),c1,c2,ll(mvsiz)
121 my_real
122 . xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4,tol
123 DATA pg/.577350269189626/
124C-----------------------------------------------
125 DO i=1,6
126 ii(i) = nel*(i-1)
127 ENDDO
128C
129 tol=em12
130 IF (isrot > 0 ) tol=em8
131 mat_1 = ixc(1,jft)
132 DO i=jft,jlt
133 mat(i) = mat_1
134 pid(i) = ixc(6,i)
135 ngl(i) = ixc(7,i)
136 ENDDO
137C----------------------------
138 DO i=jft,jlt
139 rx(i)=x(1,ixc(3,i))+x(1,ixc(4,i))-x(1,ixc(2,i))-x(1,ixc(5,i))
140 sx(i)=x(1,ixc(4,i))+x(1,ixc(5,i))-x(1,ixc(2,i))-x(1,ixc(3,i))
141 ry(i)=x(2,ixc(3,i))+x(2,ixc(4,i))-x(2,ixc(2,i))-x(2,ixc(5,i))
142 sy(i)=x(2,ixc(4,i))+x(2,ixc(5,i))-x(2,ixc(2,i))-x(2,ixc(3,i))
143 rz(i)=x(3,ixc(3,i))+x(3,ixc(4,i))-x(3,ixc(2,i))-x(3,ixc(5,i))
144 ssz(i)=x(3,ixc(4,i))+x(3,ixc(5,i))-x(3,ixc(2,i))-x(3,ixc(3,i))
145 ENDDO
146C----------------------------
147C LOCAL SYSTEM
148C----------------------------
149 k = 0
150 CALL clskew3(jft,jlt,k,
151 . rx, ry, rz,
152 . sx, sy, ssz,
153 . r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
154 DO i=jft,jlt
155 area(i)=fourth*deta1(i)
156 area_i(i)=one/area(i)
157 vq(i,1,1)=r11(i)
158 vq(i,2,1)=r21(i)
159 vq(i,3,1)=r31(i)
160 vq(i,1,2)=r12(i)
161 vq(i,2,2)=r22(i)
162 vq(i,3,2)=r32(i)
163 vq(i,1,3)=r13(i)
164 vq(i,2,3)=r23(i)
165 vq(i,3,3)=r33(i)
166 ENDDO
167C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
168 DO i=jft,jlt
169 j=ixc(2,i)
170 k=ixc(3,i)
171 l=ixc(4,i)
172 m=ixc(5,i)
173 lxyz0(1)=fourth*(x(1,l)+x(1,m)+x(1,j)+x(1,k))
174 lxyz0(2)=fourth*(x(2,l)+x(2,m)+x(2,j)+x(2,k))
175 lxyz0(3)=fourth*(x(3,l)+x(3,m)+x(3,j)+x(3,k))
176C
177 xx1=x(1,k)-x(1,j)
178 yy1=x(2,k)-x(2,j)
179 zz1=x(3,k)-x(3,j)
180C
181 xl2(i)=r11(i)*xx1+r21(i)*yy1+r31(i)*zz1
182 yl2(i)=r12(i)*xx1+r22(i)*yy1+r32(i)*zz1
183C
184 xx2=x(1,j)-lxyz0(1)
185 yy2=x(2,j)-lxyz0(2)
186 zz2=x(3,j)-lxyz0(3)
187 zl1(i)=r13(i)*xx2+r23(i)*yy2+r33(i)*zz2
188C
189 xx3=x(1,l)-x(1,j)
190 yy3=x(2,l)-x(2,j)
191 zz3=x(3,l)-x(3,j)
192 xl3(i)=r11(i)*xx3+r21(i)*yy3+r31(i)*zz3
193 yl3(i)=r12(i)*xx3+r22(i)*yy3+r32(i)*zz3
194C
195 xx4=x(1,m)-x(1,j)
196 yy4=x(2,m)-x(2,j)
197 zz4=x(3,m)-x(3,j)
198 xl4(i)=r11(i)*xx4+r21(i)*yy4+r31(i)*zz4
199 yl4(i)=r12(i)*xx4+r22(i)*yy4+r32(i)*zz4
200 ENDDO
201C----------------------------
202C SMALL STRAIN
203C----------------------------
204 IF(ismstr==1.OR.ismstr==2)THEN
205 DO i=jft,jlt
206 IF(abs(offg(i))==two)THEN
207 xl2(i)=smstr(ii(1)+i)
208 yl2(i)=smstr(ii(2)+i)
209 xl3(i)=smstr(ii(3)+i)
210 yl3(i)=smstr(ii(4)+i)
211 xl4(i)=smstr(ii(5)+i)
212 yl4(i)=smstr(ii(6)+i)
213 zl1(i)=zero
214 area(i)=half*
215 . ((xl2(i)-xl4(i))*yl3(i)-xl3(i)*(yl2(i)-yl4(i)))
216 area_i(i)=one/max(em20,area(i))
217 ELSE
218 smstr(ii(1)+i)=xl2(i)
219 smstr(ii(2)+i)=yl2(i)
220 smstr(ii(3)+i)=xl3(i)
221 smstr(ii(4)+i)=yl3(i)
222 smstr(ii(5)+i)=xl4(i)
223 smstr(ii(6)+i)=yl4(i)
224 ENDIF
225 ENDDO
226 ENDIF
227 IF(ismstr==1)THEN
228 DO i=jft,jlt
229 IF(offg(i)==one)offg(i)=two
230 ENDDO
231 ENDIF
232C----------------------------
233C Orthotropy (later)
234C----------------------------
235 IF (irep > 0) THEN
236 CALL cortdir3(elbuf_str,dir_a,dir_b ,jft ,jlt ,
237 . nlay ,irep ,rx ,ry ,rz ,
238 . sx ,sy ,ssz ,r11 ,r21 ,
239 . r31 ,r12 ,r22 ,r32 ,nel )
240 ENDIF
241c
242 DO ep=jft,jlt
243 lxyz0(1)=fourth*(xl2(ep)+xl3(ep)+xl4(ep))
244 lxyz0(2)=fourth*(yl2(ep)+yl3(ep)+yl4(ep))
245 vcore(ep,1,1)=-lxyz0(1)
246 vcore(ep,1,2)=xl2(ep)-lxyz0(1)
247 vcore(ep,1,3)=xl3(ep)-lxyz0(1)
248 vcore(ep,1,4)=xl4(ep)-lxyz0(1)
249 vcore(ep,2,1)=-lxyz0(2)
250 vcore(ep,2,2)=yl2(ep)-lxyz0(2)
251 vcore(ep,2,3)=yl3(ep)-lxyz0(2)
252 vcore(ep,2,4)=yl4(ep)-lxyz0(2)
253 x13_t(ep) =(vcore(ep,1,1)-vcore(ep,1,3))*half
254 x24_t(ep) =(vcore(ep,1,2)-vcore(ep,1,4))*half
255 y13_t(ep) =(vcore(ep,2,1)-vcore(ep,2,3))*half
256 y24_t(ep) =(vcore(ep,2,2)-vcore(ep,2,4))*half
257 l13(ep)=x13_t(ep)*x13_t(ep)+y13_t(ep)*y13_t(ep)
258 l24(ep)=x24_t(ep)*x24_t(ep)+y24_t(ep)*y24_t(ep)
259 ll(ep)=max(l13(ep),l24(ep))
260 ENDDO
261 IF (imp_chk > 0) THEN
262 DO ep=jft,jlt
263 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
264 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
265 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
266 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
267 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
268 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
269 j1=(mx23*my13-mx13*my23)*pg
270 j2=-(mx13*my34-mx34*my13)*pg
271 j0=area(ep)*fourth
272 jac(ep,1)=j0+j2-j1
273 jac(ep,2)=j0+j2+j1
274 jac(ep,3)=j0-j2+j1
275 jac(ep,4)=j0-j2-j1
276 IF(offg(ep)/=zero)THEN
277 IF(jac(ep,1)<=zero.OR.jac(ep,2)<=zero.OR.
278 . jac(ep,3)<=zero.OR.jac(ep,4)<=zero)THEN
279#include "lockon.inc"
280 WRITE(iout ,2001) ngl(ep)
281#include "lockoff.inc"
282 idel7nok = 1
283 imp_ir = imp_ir + 1
284 ENDIF
285 ENDIF
286 ENDDO
287 ENDIF
288C ------- Vectorize- a day-to-day JFT, NPLAT, JLT for all--
289 nplat=jft-1
290 splat= 0
291 DO ep=jft,jlt
292 z2=zl1(ep)*zl1(ep)
293 IF (z2<tol*ll(ep)) THEN
294 nplat=nplat+1
295 iplat(nplat)=ep
296 ELSE
297 splat=splat+1
298 jplat(splat)=ep
299 ENDIF
300 ENDDO
301 DO ep=nplat+1,jlt
302 iplat(ep)=jplat(ep-nplat)
303 ENDDO
304#include "vectorize.inc"
305 DO i=jft,nplat
306 ep =iplat(i)
307 x13 =x13_t(ep)
308 x24 =x24_t(ep)
309 y13 =y13_t(ep)
310 y24 =y24_t(ep)
311 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
312 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
313 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
314 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
315 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
316 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
317 x13_t(ep) =x13*area_i(ep)
318 x24_t(ep) =x24*area_i(ep)
319 y13_t(ep) =y13*area_i(ep)
320 y24_t(ep) =y24*area_i(ep)
321C--------GAMA(2)
322 gama1=-mx13*y24+my13*x24
323 gama2= mx13*y13-my13*x13
324 vcore(ep,1,1)=y24_t(ep)
325 vcore(ep,2,1)=-y13_t(ep)
326 vcore(ep,3,1)=-x24_t(ep)
327 vcore(ep,1,2)= x13_t(ep)
328 vcore(ep,2,2)=gama1*area_i(ep)
329 vcore(ep,3,2)=gama2*area_i(ep)
330 vcore(ep,1,3)=mx23
331 vcore(ep,2,3)=my23
332 vcore(ep,3,3)=mx34
333 vcore(ep,1,4)=my34
334 vcore(ep,2,4)=mx13
335 vcore(ep,3,4)=my13
336 j1=(mx23*my13-mx13*my23)*pg
337 j2=-(mx13*my34-mx34*my13)*pg
338 j0=area(ep)*fourth
339 jac(ep,1)=abs(j0+j2-j1)
340 jac(ep,2)=abs(j0+j2+j1)
341 jac(ep,3)=abs(j0-j2+j1)
342 jac(ep,4)=abs(j0-j2-j1)
343 j1=(my23-my34)*pg
344 j2=-(my23+my34)*pg
345 hx(ep,1)=j1/jac(ep,1)
346 hx(ep,2)=j2/jac(ep,2)
347 hx(ep,3)=-j1/jac(ep,3)
348 hx(ep,4)=-j2/jac(ep,4)
349 j1=(mx34-mx23)*pg
350 j2=(mx34+mx23)*pg
351 hy(ep,1)=j1/jac(ep,1)
352 hy(ep,2)=j2/jac(ep,2)
353 hy(ep,3)=-j1/jac(ep,3)
354 hy(ep,4)=-j2/jac(ep,4)
355 ENDDO
356#include "vectorize.inc"
357 DO i=nplat+1,jlt
358 ep =iplat(i)
359 z1=zl1(ep)
360 z2=z1*z1
361 vcore(ep,3,1)=z1
362 vcore(ep,3,2)=-z1
363 vcore(ep,3,3)=z1
364 vcore(ep,3,4)=-z1
365 x13 =x13_t(ep)
366 x24 =x24_t(ep)
367 y13 =y13_t(ep)
368 y24 =y24_t(ep)
369 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
370 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
371 mx23=(vcore(ep,1,2)+vcore(ep,1,3))*half
372 my23=(vcore(ep,2,2)+vcore(ep,2,3))*half
373 mx34=(vcore(ep,1,3)+vcore(ep,1,4))*half
374 my34=(vcore(ep,2,3)+vcore(ep,2,4))*half
375 x13_t(ep) =x13*area_i(ep)
376 x24_t(ep) =x24*area_i(ep)
377 y13_t(ep) =y13*area_i(ep)
378 y24_t(ep) =y24*area_i(ep)
379C--------GAMA(2)
380 gama1=-mx13*y24+my13*x24
381 gama2= mx13*y13-my13*x13
382C--------------------------
383C CONSTRUIRE [F], [Q], [NM], [A-S], [AS] AUX NOEUDS
384C--------------------------
385C----------------------------------------------------
386C CALCULATION OF [F]
387C----------------------------------------------------
388C---------
389C 2A1
390C---------
391 x21 =mx23-mx13
392 x34 =(vcore(ep,1,3)-vcore(ep,1,4))*half
393 y21 =my23-my13
394 y34 =(vcore(ep,2,3)-vcore(ep,2,4))*half
395 z21 = -z1
396 z34 = z1
397 l12 = sqrt(x21*x21+y21*y21+z2)
398 l34 = sqrt(x34*x34+y34*y34+z2)
399C---------
400C 2A2
401C---------
402 x41 =mx34-mx13
403 x32 =(vcore(ep,1,3)-vcore(ep,1,2))*half
404 y41 =my34-my13
405 y32 =(vcore(ep,2,3)-vcore(ep,2,2))*half
406 z41 = -z1
407 z32 = z1
408C----------
409C CALCUL DE [QN] N=1,4
410C----------
411 a_4=area(ep)*fourth
412C
413C---------- N =1----------
414 sl=one/max(l12,em20)
415 vqn(ep,1,1)=x21*sl
416 vqn(ep,2,1)=y21*sl
417 vqn(ep,3,1)=z21*sl
418 sz2=a_4-gama1
419 sz=z2*l24(ep)
420 sl=one/sqrt(sz+sz2*sz2)
421 vqn(ep,7,1)=-z1*y24
422 vqn(ep,8,1)= z1*x24
423 vqn(ep,9,1)= sz2*sl
424C
425 vqn(ep,7,3)=-vqn(ep,7,1)
426 vqn(ep,8,3)=-vqn(ep,8,1)
427 vqn(ep,7,1)= vqn(ep,7,1)*sl
428 vqn(ep,8,1)= vqn(ep,8,1)*sl
429C
430 vqn(ep,4,1)= vqn(ep,8,1)*vqn(ep,3,1)-vqn(ep,9,1)*vqn(ep,2,1)
431 vqn(ep,5,1)= vqn(ep,9,1)*vqn(ep,1,1)-vqn(ep,7,1)*vqn(ep,3,1)
432 vqn(ep,6,1)= vqn(ep,7,1)*vqn(ep,2,1)-vqn(ep,8,1)*vqn(ep,1,1)
433C---------- N =3----------
434 sl=one/max(l34,em20)
435 vqn(ep,1,3)=x34*sl
436 vqn(ep,2,3)=y34*sl
437 vqn(ep,3,3)=z34*sl
438 sz2=a_4+gama1
439 sl=one/sqrt(sz+sz2*sz2)
440 vqn(ep,7,3)= vqn(ep,7,3)*sl
441 vqn(ep,8,3)= vqn(ep,8,3)*sl
442 vqn(ep,9,3)= sz2*sl
443C
444 vqn(ep,4,3)= vqn(ep,8,3)*vqn(ep,3,3)-vqn(ep,9,3)*vqn(ep,2,3)
445 vqn(ep,5,3)= vqn(ep,9,3)*vqn(ep,1,3)-vqn(ep,7,3)*vqn(ep,3,3)
446 vqn(ep,6,3)= vqn(ep,7,3)*vqn(ep,2,3)-vqn(ep,8,3)*vqn(ep,1,3)
447C---------- N =2----------
448 vqn(ep,1,2)=vqn(ep,1,1)
449 vqn(ep,2,2)=vqn(ep,2,1)
450 vqn(ep,3,2)=vqn(ep,3,1)
451 sz2=a_4+gama2
452 sz=z2*l13(ep)
453 sl=one/sqrt(sz+sz2*sz2)
454 vqn(ep,7,2)=-z1*y13
455 vqn(ep,8,2)= z1*x13
456 vqn(ep,9,2)= sz2*sl
457 vqn(ep,7,4)=-vqn(ep,7,2)
458 vqn(ep,8,4)=-vqn(ep,8,2)
459 vqn(ep,7,2)= vqn(ep,7,2)*sl
460 vqn(ep,8,2)= vqn(ep,8,2)*sl
461C
462 vqn(ep,4,2)= vqn(ep,8,2)*vqn(ep,3,2)-vqn(ep,9,2)*vqn(ep,2,2)
463 vqn(ep,5,2)= vqn(ep,9,2)*vqn(ep,1,2)-vqn(ep,7,2)*vqn(ep,3,2)
464 vqn(ep,6,2)= vqn(ep,7,2)*vqn(ep,2,2)-vqn(ep,8,2)*vqn(ep,1,2)
465C---------- N =4----------
466 vqn(ep,1,4)=vqn(ep,1,3)
467 vqn(ep,2,4)=vqn(ep,2,3)
468 vqn(ep,3,4)=vqn(ep,3,3)
469 sz2=a_4-gama2
470 sl=one/sqrt(sz+sz2*sz2)
471 vqn(ep,7,4)= vqn(ep,7,4)*sl
472 vqn(ep,8,4)= vqn(ep,8,4)*sl
473 vqn(ep,9,4)= sz2*sl
474C
475 vqn(ep,4,4)= vqn(ep,8,4)*vqn(ep,3,4)-vqn(ep,9,4)*vqn(ep,2,4)
476 vqn(ep,5,4)= vqn(ep,9,4)*vqn(ep,1,4)-vqn(ep,7,4)*vqn(ep,3,4)
477 vqn(ep,6,4)= vqn(ep,7,4)*vqn(ep,2,4)-vqn(ep,8,4)*vqn(ep,1,4)
478C--------------------------------------------------
479C CALCULATION OF AS N AT THE MIDDLE OF SIDES
480C--------------------------------------------------
481C J=1 COTE
482 vnrm(ep,1,1)=vqn(ep,7,1)+vqn(ep,7,2)
483 vnrm(ep,2,1)=vqn(ep,8,1)+vqn(ep,8,2)
484 vnrm(ep,3,1)=vqn(ep,9,1)+vqn(ep,9,2)
485 c1=sqrt(vnrm(ep,1,1)*vnrm(ep,1,1)+
486 1 vnrm(ep,2,1)*vnrm(ep,2,1)+vnrm(ep,3,1)*vnrm(ep,3,1))
487 c1 = max(em20,c1)
488 vnrm(ep,1,1)=vnrm(ep,1,1)/c1
489 vnrm(ep,2,1)=vnrm(ep,2,1)/c1
490 vnrm(ep,3,1)=vnrm(ep,3,1)/c1
491 vastn(ep,1,1)=zero
492 vastn(ep,2,1)=l12
493 vastn(ep,3,1)=vastn(ep,1,1)
494 vastn(ep,4,1)=vastn(ep,2,1)
495C J=2
496 vnrm(ep,1,2)=vqn(ep,7,4)+vqn(ep,7,3)
497 vnrm(ep,2,2)=vqn(ep,8,4)+vqn(ep,8,3)
498 vnrm(ep,3,2)=vqn(ep,9,4)+vqn(ep,9,3)
499 c1=sqrt(vnrm(ep,1,2)*vnrm(ep,1,2)+
500 1 vnrm(ep,2,2)*vnrm(ep,2,2)+vnrm(ep,3,2)*vnrm(ep,3,2))
501 c1 = max(em20,c1)
502 vnrm(ep,1,2)=vnrm(ep,1,2)/c1
503 vnrm(ep,2,2)=vnrm(ep,2,2)/c1
504 vnrm(ep,3,2)=vnrm(ep,3,2)/c1
505 vastn(ep,1,2)=zero
506 vastn(ep,2,2)=l34
507 vastn(ep,3,2)=vastn(ep,1,2)
508 vastn(ep,4,2)=vastn(ep,2,2)
509C J=3
510 vnrm(ep,1,3)=vqn(ep,7,1)+vqn(ep,7,4)
511 vnrm(ep,2,3)=vqn(ep,8,1)+vqn(ep,8,4)
512 vnrm(ep,3,3)=vqn(ep,9,1)+vqn(ep,9,4)
513 c1=sqrt(vnrm(ep,1,3)*vnrm(ep,1,3)+
514 1 vnrm(ep,2,3)*vnrm(ep,2,3)+vnrm(ep,3,3)*vnrm(ep,3,3))
515 c1 = max(em20,c1)
516 vnrm(ep,1,3)=vnrm(ep,1,3)/c1
517 vnrm(ep,2,3)=vnrm(ep,2,3)/c1
518 vnrm(ep,3,3)=vnrm(ep,3,3)/c1
519 vastn(ep,1,3)=-(x41*vqn(ep,4,1)+y41*vqn(ep,5,1)+z41*vqn(ep,6,1))
520 vastn(ep,2,3)= x41*vqn(ep,1,1)+y41*vqn(ep,2,1)+z41*vqn(ep,3,1)
521 vastn(ep,3,3)=-(x41*vqn(ep,4,4)+y41*vqn(ep,5,4)+z41*vqn(ep,6,4))
522 vastn(ep,4,3)= x41*vqn(ep,1,4)+y41*vqn(ep,2,4)+z41*vqn(ep,3,4)
523C J=4
524 vnrm(ep,1,4)=vqn(ep,7,2)+vqn(ep,7,3)
525 vnrm(ep,2,4)=vqn(ep,8,2)+vqn(ep,8,3)
526 vnrm(ep,3,4)=vqn(ep,9,2)+vqn(ep,9,3)
527 c1=sqrt(vnrm(ep,1,4)*vnrm(ep,1,4)+
528 1 vnrm(ep,2,4)*vnrm(ep,2,4)+vnrm(ep,3,4)*vnrm(ep,3,4))
529 c1 = max(em20,c1)
530 vnrm(ep,1,4)=vnrm(ep,1,4)/c1
531 vnrm(ep,2,4)=vnrm(ep,2,4)/c1
532 vnrm(ep,3,4)=vnrm(ep,3,4)/c1
533 vastn(ep,1,4)=-(x32*vqn(ep,4,2)+y32*vqn(ep,5,2)+z32*vqn(ep,6,2))
534 vastn(ep,2,4)= x32*vqn(ep,1,2)+y32*vqn(ep,2,2)+z32*vqn(ep,3,2)
535 vastn(ep,3,4)=-(x32*vqn(ep,4,3)+y32*vqn(ep,5,3)+z32*vqn(ep,6,3))
536 vastn(ep,4,4)= x32*vqn(ep,1,3)+y32*vqn(ep,2,3)+z32*vqn(ep,3,3)
537C----------
538C CALCUL DE [QG] NG=1,4
539C----------
540 a_4=a_4/pg
541 jmx13=mx13*pg
542 jmy13=my13*pg
543 jmz13=z1*pg
544 j2myz=jmz13*jmz13
545C---------- NG =1----------
546 sz2=a_4-gama1
547 sz=z2*l24(ep)
548 sl=sqrt(sz+sz2*sz2)
549 jac(ep,1)=sl*pg
550 sl=one/max(sl,em20)
551 vqg(ep,7,1)=-z1*y24
552 vqg(ep,8,1)= z1*x24
553 vqg(ep,9,1)= sz2*sl
554 vqg(ep,7,3)=-vqg(ep,7,1)
555 vqg(ep,8,3)=-vqg(ep,8,1)
556 vqg(ep,7,1)= vqg(ep,7,1)*sl
557 vqg(ep,8,1)= vqg(ep,8,1)*sl
558C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
559 g1x1=mx23-jmx13
560 g1y1=my23-jmy13
561 c1 = sqrt(g1x1*g1x1 + g1y1*g1y1 +j2myz )
562 g2x1=mx34-jmx13
563 g2y1=my34-jmy13
564 c2 = sqrt(g2x1*g2x1 + g2y1*g2y1 +j2myz)
565 vjfi(ep,1,1)=(g2y1*vqg(ep,9,1)+jmz13*vqg(ep,8,1))
566 vjfi(ep,2,1)=(-jmz13*vqg(ep,7,1)-g2x1*vqg(ep,9,1))
567 vjfi(ep,3,1)=(g2x1*vqg(ep,8,1)-g2y1*vqg(ep,7,1))
568 vqg(ep,1,1)= g1x1*c2+vjfi(ep,1,1)*c1
569 vqg(ep,2,1)= g1y1*c2+vjfi(ep,2,1)*c1
570 vqg(ep,3,1)= -jmz13*c2+vjfi(ep,3,1)*c1
571 sl=sl/pg
572 vjfi(ep,1,1)=vjfi(ep,1,1)*sl
573 vjfi(ep,2,1)=vjfi(ep,2,1)*sl
574 vjfi(ep,3,1)=vjfi(ep,3,1)*sl
575C
576 vjfi(ep,4,1)=(-jmz13*vqg(ep,8,1)-g1y1*vqg(ep,9,1))*sl
577 vjfi(ep,5,1)=(g1x1*vqg(ep,9,1)+jmz13*vqg(ep,7,1))*sl
578 vjfi(ep,6,1)=(g1y1*vqg(ep,7,1)-g1x1*vqg(ep,8,1))*sl
579C
580 sl = sqrt(vqg(ep,1,1)*vqg(ep,1,1) + vqg(ep,2,1)*vqg(ep,2,1)
581 1 + vqg(ep,3,1)*vqg(ep,3,1))
582 IF ( sl/=zero) sl = one / sl
583 vqg(ep,1,1) = vqg(ep,1,1)*sl
584 vqg(ep,2,1) = vqg(ep,2,1)*sl
585 vqg(ep,3,1) = vqg(ep,3,1)*sl
586C
587 vqg(ep,4,1)= vqg(ep,8,1)*vqg(ep,3,1)-vqg(ep,9,1)*vqg(ep,2,1)
588 vqg(ep,5,1)= vqg(ep,9,1)*vqg(ep,1,1)-vqg(ep,7,1)*vqg(ep,3,1)
589 vqg(ep,6,1)= vqg(ep,7,1)*vqg(ep,2,1)-vqg(ep,8,1)*vqg(ep,1,1)
590C---------- NG =3----------
591 sz2=a_4+gama1
592 sl=sqrt(sz+sz2*sz2)
593 jac(ep,3)=sl*pg
594 sl=one/max(sl,em20)
595 vqg(ep,7,3)= vqg(ep,7,3)*sl
596 vqg(ep,8,3)= vqg(ep,8,3)*sl
597 vqg(ep,9,3)= sz2*sl
598C
599 g1x3=mx23+jmx13
600 g1y3=my23+jmy13
601 j1 = sqrt(g1x3*g1x3 + g1y3*g1y3 +j2myz )
602C--------G2X3=G2X2------
603 g2x2=mx34+jmx13
604 g2y2=my34+jmy13
605 j2 = sqrt(g2x2*g2x2 + g2y2*g2y2 +j2myz)
606 vjfi(ep,1,3)=(g2y2*vqg(ep,9,3)-jmz13*vqg(ep,8,3))
607 vjfi(ep,2,3)=(jmz13*vqg(ep,7,3)-g2x2*vqg(ep,9,3))
608 vjfi(ep,3,3)=(g2x2*vqg(ep,8,3)-g2y2*vqg(ep,7,3))
609 vqg(ep,1,3)= g1x3*j2+vjfi(ep,1,3)*j1
610 vqg(ep,2,3)= g1y3*j2+vjfi(ep,2,3)*j1
611 vqg(ep,3,3)= jmz13*j2+vjfi(ep,3,3)*j1
612 sl=sl/pg
613 vjfi(ep,1,3)=vjfi(ep,1,3)*sl
614 vjfi(ep,2,3)=vjfi(ep,2,3)*sl
615 vjfi(ep,3,3)=vjfi(ep,3,3)*sl
616C
617 vjfi(ep,4,3)=(jmz13*vqg(ep,8,3)-g1y3*vqg(ep,9,3))*sl
618 vjfi(ep,5,3)=(g1x3*vqg(ep,9,3)-jmz13*vqg(ep,7,3))*sl
619 vjfi(ep,6,3)=(g1y3*vqg(ep,7,3)-g1x3*vqg(ep,8,3))*sl
620C
621 sl = sqrt(vqg(ep,1,3)*vqg(ep,1,3) + vqg(ep,2,3)*vqg(ep,2,3)
622 1 + vqg(ep,3,3)*vqg(ep,3,3))
623 IF ( sl/=zero) sl = one / sl
624 vqg(ep,1,3) = vqg(ep,1,3)*sl
625 vqg(ep,2,3) = vqg(ep,2,3)*sl
626 vqg(ep,3,3) = vqg(ep,3,3)*sl
627C
628 vqg(ep,4,3)= vqg(ep,8,3)*vqg(ep,3,3)-vqg(ep,9,3)*vqg(ep,2,3)
629 vqg(ep,5,3)= vqg(ep,9,3)*vqg(ep,1,3)-vqg(ep,7,3)*vqg(ep,3,3)
630 vqg(ep,6,3)= vqg(ep,7,3)*vqg(ep,2,3)-vqg(ep,8,3)*vqg(ep,1,3)
631C---------- NG =2----------
632 sz2=a_4+gama2
633 sz=z2*l13(ep)
634 sl=sqrt(sz+sz2*sz2)
635 jac(ep,2)=sl*pg
636 sl=one/max(sl,em20)
637 vqg(ep,7,2)=-z1*y13
638 vqg(ep,8,2)= z1*x13
639 vqg(ep,9,2)= sz2*sl
640 vqg(ep,7,4)=-vqg(ep,7,2)
641 vqg(ep,8,4)=-vqg(ep,8,2)
642 vqg(ep,7,2)= vqg(ep,7,2)*sl
643 vqg(ep,8,2)= vqg(ep,8,2)*sl
644C
645 vjfi(ep,1,2)=(g2y2*vqg(ep,9,2)-jmz13*vqg(ep,8,2))
646 vjfi(ep,2,2)=(jmz13*vqg(ep,7,2)-g2x2*vqg(ep,9,2))
647 vjfi(ep,3,2)=(g2x2*vqg(ep,8,2)-g2y2*vqg(ep,7,2))
648 vqg(ep,1,2)= g1x1*j2+vjfi(ep,1,2)*c1
649 vqg(ep,2,2)= g1y1*j2+vjfi(ep,2,2)*c1
650 vqg(ep,3,2)=-jmz13*j2+vjfi(ep,3,2)*c1
651 sl=sl/pg
652 vjfi(ep,1,2)=vjfi(ep,1,2)*sl
653 vjfi(ep,2,2)=vjfi(ep,2,2)*sl
654 vjfi(ep,3,2)=vjfi(ep,3,2)*sl
655C
656 vjfi(ep,4,2)=(-jmz13*vqg(ep,8,2)-g1y1*vqg(ep,9,2))*sl
657 vjfi(ep,5,2)=(g1x1*vqg(ep,9,2)+jmz13*vqg(ep,7,2))*sl
658 vjfi(ep,6,2)=(g1y1*vqg(ep,7,2)-g1x1*vqg(ep,8,2))*sl
659C
660 sl = sqrt(vqg(ep,1,2)*vqg(ep,1,2) + vqg(ep,2,2)*vqg(ep,2,2)
661 1 + vqg(ep,3,2)*vqg(ep,3,2))
662 IF ( sl/=0.) sl = 1. / sl
663 vqg(ep,1,2) = vqg(ep,1,2)*sl
664 vqg(ep,2,2) = vqg(ep,2,2)*sl
665 vqg(ep,3,2) = vqg(ep,3,2)*sl
666 vqg(ep,4,2)= vqg(ep,8,2)*vqg(ep,3,2)-vqg(ep,9,2)*vqg(ep,2,2)
667 vqg(ep,5,2)= vqg(ep,9,2)*vqg(ep,1,2)-vqg(ep,7,2)*vqg(ep,3,2)
668 vqg(ep,6,2)= vqg(ep,7,2)*vqg(ep,2,2)-vqg(ep,8,2)*vqg(ep,1,2)
669C---------- NG =4----------
670 sz2=a_4-gama2
671 sl=sqrt(sz+sz2*sz2)
672 jac(ep,4)=sl*pg
673 sl=one/max(sl,em20)
674 vqg(ep,7,4)= vqg(ep,7,4)*sl
675 vqg(ep,8,4)= vqg(ep,8,4)*sl
676 vqg(ep,9,4)= sz2*sl
677C
678 vjfi(ep,1,4)=(g2y1*vqg(ep,9,4)+jmz13*vqg(ep,8,4))
679 vjfi(ep,2,4)=(-jmz13*vqg(ep,7,4)-g2x1*vqg(ep,9,4))
680 vjfi(ep,3,4)=(g2x1*vqg(ep,8,4)-g2y1*vqg(ep,7,4))
681 vqg(ep,1,4)= g1x3*c2+vjfi(ep,1,4)*j1
682 vqg(ep,2,4)= g1y3*c2+vjfi(ep,2,4)*j1
683 vqg(ep,3,4)=jmz13*c2+vjfi(ep,3,4)*j1
684 sl=sl/pg
685 vjfi(ep,1,4)=vjfi(ep,1,4)*sl
686 vjfi(ep,2,4)=vjfi(ep,2,4)*sl
687 vjfi(ep,3,4)=vjfi(ep,3,4)*sl
688C
689 vjfi(ep,4,4)=(jmz13*vqg(ep,8,4)-g1y3*vqg(ep,9,4))*sl
690 vjfi(ep,5,4)=(g1x3*vqg(ep,9,4)-jmz13*vqg(ep,7,4))*sl
691 vjfi(ep,6,4)=(g1y3*vqg(ep,7,4)-g1x3*vqg(ep,8,4))*sl
692C
693 sl = sqrt(vqg(ep,1,4)*vqg(ep,1,4) + vqg(ep,2,4)*vqg(ep,2,4)
694 1 + vqg(ep,3,4)*vqg(ep,3,4))
695 IF ( sl/=zero) sl = one / sl
696 vqg(ep,1,4) = vqg(ep,1,4)*sl
697 vqg(ep,2,4) = vqg(ep,2,4)*sl
698 vqg(ep,3,4) = vqg(ep,3,4)*sl
699 vqg(ep,4,4)= vqg(ep,8,4)*vqg(ep,3,4)-vqg(ep,9,4)*vqg(ep,2,4)
700 vqg(ep,5,4)= vqg(ep,9,4)*vqg(ep,1,4)-vqg(ep,7,4)*vqg(ep,3,4)
701 vqg(ep,6,4)= vqg(ep,7,4)*vqg(ep,2,4)-vqg(ep,8,4)*vqg(ep,1,4)
702 area(ep)=jac(ep,1)+jac(ep,2)+jac(ep,3)+jac(ep,4)
703 ENDDO
704C
705 DO i=jft,jlt
706 off(i)=offg(i)
707 ENDDO
708C
709 RETURN
710 2001 FORMAT(/' ZERO OR NEGATIVE SHELL SUB-AREA : ELEMENT NB:',
711 . i8/)
712 END
cbacoork
subroutine cbacoork(jft, jlt, x, ixc, pm, offg, geo, area, vcore, jac, hx, hy, vqn, vqg, vq, vjfi, vnrm, vastn, nplat, iplat, x13_t, x24_t, y13_t, y24_t, elbuf_str, nlay, smstr, irep, npt, ismstr, dir_a, dir_b, pid, mat, ngl, off, isrot, nel)
Definition cbacoork.F:41
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
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
max
#define max(a, b)
Definition macros.h:21