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