OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
corthdir.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!|| corthdir ../starter/source/elements/shell/coque/corthdir.f
25!||--- called by ------------------------------------------------------
26!|| c3init3 ../starter/source/elements/sh3n/coque3n/c3init3.F
27!|| cinit3 ../starter/source/elements/shell/coque/cinit3.F
28!|| cmaini3 ../starter/source/elements/sh3n/coquedk/cmaini3.F
29!||--- uses -----------------------------------------------------
30!|| stack_mod ../starter/share/modules1/stack_mod.F
31!||====================================================================
32 SUBROUTINE corthdir(ELBUF_STR,
33 . IGEO ,GEO ,VX ,VY ,VZ ,
34 . PHI1 ,PHI2 ,COOR1 ,COOR2 ,COOR3 ,
35 . COOR4 ,IORTHLOC ,NLAY ,IREP ,ISUBSTACK,
36 . STACK ,GEO_STACK ,IGEO_STACK,IR ,IS ,
37 . NEL ,IMAT ,IPROP ,
38 . X1 ,X2 ,X3 ,X4 ,Y1 ,Y2 ,
39 . Y3 ,Y4 ,Z1 ,Z2 ,Z3 ,Z4 ,
40 . E1X, E2X, E3X, E1Y, E2Y, E3Y ,E1Z, E2Z, E3Z ,
41 . NPT_ALL, IDRAPE)
42C-----------------------------------------------
43C M o d u l e s
44C-----------------------------------------------
45 USE elbufdef_mod
46 USE stack_mod
47C-----------------------------------------------
48C INITIALISE LES COORDONEES LOCALES DES AXES D'ORTHOTROPIE
49C INITIALISE LES EPAISSEURS ET MATERIAUX DES COUCHES
50C-----------------------------------------------
51C I m p l i c i t T y p e s
52C-----------------------------------------------
53#include "implicit_f.inc"
54C-----------------------------------------------
55C G l o b a l P a r a m e t e r s
56C-----------------------------------------------
57#include "mvsiz_p.inc"
58C-----------------------------------------------
59C C o m m o n B l o c k s
60C-----------------------------------------------
61#include "vect01_c.inc"
62#include "param_c.inc"
63C-----------------------------------------------
64C D u m m y A r g u m e n t s
65C-----------------------------------------------
66 INTEGER NLAY,IREP,ISUBSTACK,IS,IR,NEL,IMAT,IPROP,NPT_ALL,IDRAPE
67 INTEGER IGEO(NPROPGI,*),IORTHLOC(*),IGEO_STACK(NPROPGI,*)
68 my_real, DIMENSION(MVSIZ), INTENT(IN) :: E1X,E2X,E3X,E1Y,E2Y,E3Y,E1Z,E2Z,E3Z,
69 . X1,X2,X3,X4,Y1,Y2,Y3,Y4,Z1,Z2,Z3,Z4
70 my_real
71 . GEO(NPROPG,*),VX(*),VY(*),VZ(*),PHI1(NPT_ALL,*),PHI2(NPT_ALL,*),
72 . COOR1(NPT_ALL,MVSIZ),COOR2(NPT_ALL,MVSIZ),
73 . COOR3(NPT_ALL,MVSIZ),COOR4(NPT_ALL,MVSIZ),
74 . GEO_STACK(NPROPG,*)
75 TYPE(elbuf_struct_), TARGET :: ELBUF_STR
76 TYPE (STACK_PLY):: STACK
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER I,II,J,I1,I2,I3,IL,IGTYP,IPTHK,IPMAT,ILAW,L_DMG,NPTT,
81 . IPT ,IPT_ALL,IT,IPID_PLY,IRP
82 my_real R,S,D1,D2,D11,D12,D21,D22,U1X,U1Y,U2X,U2Y,
83 . DET,W1X,W2X,W1Y,W2Y,CSP,SNP
84 my_real, DIMENSION(MVSIZ) :: e11,e12,e13,e21,e22,e23,vr,vs
85 my_real, DIMENSION(:), POINTER :: dir1,dir2,dir_dmg
86C
87 TYPE(buf_lay_) ,POINTER :: BUFLY
88 TYPE(L_BUFEL_) ,POINTER :: LBUF
89 TYPE(l_bufel_dir_) ,POINTER :: LBUF_DIR
90C======================================================================|
91 igtyp = igeo(11,iprop)
92C
93 IF (igtyp == 1) THEN
94C non orthotropic property
95 IF( elbuf_str%BUFLY(1)%LY_DIRA /= 0) THEN
96 dir1 => elbuf_str%BUFLY(1)%DIRA
97 DO i=lft,llt
98 dir1(i) = one
99 dir1(i+nel) = zero
100 ENDDO
101 ENDIF
102 ELSE
103 ipmat = 100
104 ipthk = 300
105C
106C--- isoparametric (material) axes
107 IF (ity == 3) THEN
108C--- shell 4N
109 DO i=lft,llt
110 e11(i)= x2(i)+x3(i)-x1(i)-x4(i)
111 e12(i)= y2(i)+y3(i)-y1(i)-y4(i)
112 e13(i)= z2(i)+z3(i)-z1(i)-z4(i)
113 e21(i)= x3(i)+x4(i)-x1(i)-x2(i)
114 e22(i)= y3(i)+y4(i)-y1(i)-y2(i)
115 e23(i)= z3(i)+z4(i)-z1(i)-z2(i)
116 ENDDO
117 ELSE
118C--- shell 3N
119 DO i=lft,llt
120 e11(i)= x2(i)-x1(i)
121 e12(i)= y2(i)-y1(i)
122 e13(i)= z2(i)-z1(i)
123 e21(i)= x3(i)-x1(i)
124 e22(i)= y3(i)-y1(i)
125 e23(i)= z3(i)-z1(i)
126 ENDDO
127 ENDIF
128C--- vecteur de projection dans repere elementaire
129 irp = igeo(14,iprop)
130 IF(irp == 26 ) THEN
131 DO i=lft,llt
132 vr(i)= one
133 vs(i)= zero
134 ENDDO
135 ELSE
136 DO i=lft,llt
137 vr(i)=vx(i)*e1x(i)+vy(i)*e1y(i)+vz(i)*e1z(i)
138 vs(i)=vx(i)*e2x(i)+vy(i)*e2y(i)+vz(i)*e2z(i)
139 ENDDO
140 ENDIF
141C-------
142 IF (igtyp == 9) THEN
143 dir1 => elbuf_str%BUFLY(1)%DIRA
144C-------
145 DO i=lft,llt
146 csp=cos(phi1(1,i))
147 snp=sin(phi1(1,i))
148 dir1(i) = vr(i)*csp-vs(i)*snp
149 dir1(i+nel) = vs(i)*csp+vr(i)*snp
150 IF (iorthloc(i) /= 0) THEN
151 dir1(i) = coor1(1,i)
152 dir1(i+nel) = coor2(1,i)
153 ENDIF
154 ENDDO
155C-------
156 ELSEIF (igtyp == 10) THEN
157C-------
158 DO il=1,nlay
159 dir1 => elbuf_str%BUFLY(il)%DIRA
160 DO i=lft,llt
161 csp=cos(phi1(il,i))
162 snp=sin(phi1(il,i))
163 dir1(i) = vr(i)*csp-vs(i)*snp
164 dir1(i+nel) = vs(i)*csp+vr(i)*snp
165 IF (iorthloc(i) /= 0) THEN
166 dir1(i) = coor1(il,i)
167 dir1(i+nel) = coor2(il,i)
168 ENDIF
169 ENDDO
170 ENDDO
171 DO il=1,nlay
172 i2=ipthk+il
173 i3=ipmat+il
174 DO i=lft,llt
175 geo(i2,iprop) = one/nlay
176 igeo(i3,iprop) = imat
177 ENDDO
178 ENDDO
179C-------
180 ELSEIF (igtyp == 11 .OR. igtyp == 17) THEN
181C-------
182 DO il=1,nlay
183 dir1 => elbuf_str%BUFLY(il)%DIRA
184 DO i=lft,llt
185 csp=cos(phi1(il,i))
186 snp=sin(phi1(il,i))
187 dir1(i) = vr(i)*csp-vs(i)*snp
188 dir1(i+nel) = vs(i)*csp+vr(i)*snp
189 IF (iorthloc(i) /= 0) THEN
190 dir1(i) = coor1(il,i)
191 dir1(i+nel) = coor2(il,i)
192 ENDIF
193C
194 ENDDO
195 ENDDO
196 l_dmg = elbuf_str%BUFLY(1)%L_DMG
197 IF(l_dmg > 0) THEN
198 DO il=1,nlay
199 ilaw = elbuf_str%BUFLY(il)%ILAW
200 IF (ilaw == 27) THEN
201 lbuf => elbuf_str%BUFLY(il)%LBUF(ir,is,1)
202 dir_dmg => lbuf%DMG(1:l_dmg*llt)
203 dir_dmg(lft:llt) = one
204 IF(l_dmg > 1) dir_dmg(llt+1:l_dmg*llt) = zero
205 ENDIF
206 ENDDO
207 ENDIF
208 IF (irep == 1) THEN
209 DO i=lft,llt
210 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
211 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
212 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
213 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
214 det = u1x*u2y-u1y*u2x
215 w1x = u2y/det
216 w2y = u1x/det
217 w1y = -u1y/det
218 w2x = -u2x/det
219 DO il=1,nlay
220 dir1 => elbuf_str%BUFLY(il)%DIRA
221 IF ( iorthloc(i) /= 0 ) THEN
222 dir1(i) = coor1(il,i)
223 dir1(i+nel) = coor2(il,i)
224 ENDIF
225 d1 = dir1(i)
226 d2 = dir1(i+nel)
227 dir1(i) = w1x*d1 + w2x*d2
228 dir1(i+nel) = w1y*d1 + w2y*d2
229 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
230 dir1(i) = dir1(i)/s
231 dir1(i+nel) = dir1(i+nel)/s
232 ENDDO
233 ENDDO
234 ENDIF
235C-------
236 ELSEIF (igtyp == 51 .OR. igtyp == 52) THEN
237 IF(idrape > 0) THEN
238 ipt_all = 0
239 DO il=1,nlay
240 nptt = elbuf_str%BUFLY(il)%NPTT
241 ipid_ply=stack%IGEO(2 + il,isubstack)
242 IF(ipid_ply > 0)THEN
243 DO it = 1, nptt
244 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
245 dir1 => lbuf_dir%DIRA
246 ipt = ipt_all + it
247 DO i=lft,llt
248 csp=cos(phi1(ipt,i))
249 snp=sin(phi1(ipt,i))
250 dir1(i) = vr(i)*csp-vs(i)*snp
251 dir1(i+nel) = vs(i)*csp+vr(i)*snp
252 IF (iorthloc(i) /= 0) THEN
253 dir1(i) = coor1(ipt,i)
254 dir1(i+nel) = coor2(ipt,i)
255 ENDIF
256 ENDDO
257 ENDDO ! NPPT
258 ENDIF !
259 ipt_all = ipt_all + nptt
260 ENDDO ! NLAY
261 IF(irep == 1) THEN
262 DO i=lft,llt
263 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
264 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
265 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
266 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
267 det = u1x*u2y-u1y*u2x
268 w1x = u2y/det
269 w2y = u1x/det
270 w1y = -u1y/det
271 w2x = -u2x/det
272 ipt_all = 0
273 DO il=1,nlay
274 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
275 dir1 => lbuf_dir%DIRA
276 nptt = elbuf_str%BUFLY(il)%NPTT
277 DO it = 1,nptt
278 ipt = ipt_all + it
279 IF (iorthloc(i) /= 0) THEN
280 dir1(i) = coor1(ipt,i) ! il --> It
281 dir1(i+nel) = coor2(ipt,i)
282 ENDIF
283 d1 = dir1(i)
284 d2 = dir1(i+nel)
285 dir1(i) = w1x*d1 + w2x*d2
286 dir1(i+nel) = w1y*d1 + w2y*d2
287 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
288 dir1(i) = dir1(i)/s
289 dir1(i+nel) = dir1(i+nel)/s
290 ENDDO ! nptt
291 ENDDO ! nlay
292 ENDDO
293 ELSEIF (irep == 2) THEN
294C--- Axe I d'anisotropie
295 DO il=1,nlay
296 nptt = elbuf_str%BUFLY(il)%NPTT
297 DO it = 1, nptt
298 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
299 dir1 => lbuf_dir%DIRA
300 dir2 => lbuf_dir%DIRB
301 ipt = ipt_all + it
302C--- Axe I d'anisotropie
303 DO i=lft,llt
304 csp=cos(phi1(ipt,i))
305 snp=sin(phi1(ipt,i))
306 dir1(i) = vr(i)*csp-vs(i)*snp
307 dir1(i+nel) = vs(i)*csp+vr(i)*snp
308 ENDDO
309C--- Axe II d'anisotropie
310 DO i=lft,llt
311 csp=cos(phi2(ipt,i))
312 snp=sin(phi2(ipt,i))
313 r = dir1(i)
314 s = dir1(i+nel)
315 dir2(i) = r*csp-s*snp
316 dir2(i+nel) = s*csp+r*snp
317 ENDDO
318 ENDDO ! NPTT
319 ipt_all = ipt_all + nptt
320 ENDDO ! DO IL=1,NLAY
321C---
322 ipt_all = 0
323 DO il=1,nlay
324 nptt = elbuf_str%BUFLY(il)%NPTT
325 DO it = 1, nptt
326 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
327 dir1 => lbuf_dir%DIRA
328 dir2 => lbuf_dir%DIRB
329 ipt = ipt_all + it
330 DO i=lft,llt
331 IF (iorthloc(i) /= 0) THEN
332 dir1(i) = coor1(ipt,i) !! IL ou ipt
333 dir1(i+nel) = coor2(ipt,i)
334 dir2(i) = coor3(ipt,i)
335 dir2(i+nel) = coor4(ipt,i)
336 ENDIF
337 ENDDO
338 ENDDO ! NPTT
339 ipt_all = ipt_all + nptt
340 ENDDO ! DO IL=1,NLAY
341C---
342 ipt_all = 0
343 DO i=lft,llt
344 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
345 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
346 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
347 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
348 det = u1x*u2y-u1y*u2x
349 w1x = u2y/det
350 w2y = u1x/det
351 w1y = -u1y/det
352 w2x = -u2x/det
353 DO il=1,nlay
354 nptt = elbuf_str%BUFLY(il)%NPTT
355 DO it = 1, nptt
356 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
357 dir1 => lbuf_dir%DIRA
358 dir2 => lbuf_dir%DIRB
359 ipt = ipt_all + it
360C--- dir I
361 d11 = dir1(i)
362 d21 = dir1(i+nel)
363 dir1(i) = w1x*d11 + w2x*d21
364 dir1(i+nel) = w1y*d11 + w2y*d21
365 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
366 dir1(i) = dir1(i)/s
367 dir1(i+nel) = dir1(i+nel)/s
368C--- dir II
369 d12 = dir2(i)
370 d22= dir2(i+nel)
371 dir2(i) = w1x*d12 + w2x*d22
372 dir2(i+nel) = w1y*d12 + w2y*d22
373 s = sqrt(dir2(i)**2 + dir2(i+nel)**2)
374 dir2(i) = dir2(i)/s
375 dir2(i+nel) = dir2(i+nel)/s
376 ENDDO
377 ENDDO
378 ENDDO
379 ELSEIF (irep == 3) THEN
380C mixing law58 with other user laws with IREP = 0 within PID51
381 ipt_all = 0
382 DO il=1,nlay
383 ilaw = elbuf_str%BUFLY(il)%ILAW
384 nptt = elbuf_str%BUFLY(il)%NPTT
385 IF (ilaw == 58) THEN
386 DO it = 1, nptt
387 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
388 dir1 => lbuf_dir%DIRA
389 dir2 => lbuf_dir%DIRB
390 ipt = ipt_all + it
391C--- Axe I d'anisotropie
392 DO i=lft,llt
393 csp=cos(phi1(ipt,i))
394 snp=sin(phi1(ipt,i))
395 dir1(i) = vr(i)*csp-vs(i)*snp
396 dir1(i+nel) = vs(i)*csp+vr(i)*snp
397 ENDDO
398C--- Axe II d'anisotropie
399 DO i=lft,llt
400 csp=cos(phi2(ipt,i))
401 snp=sin(phi2(ipt,i))
402 r = dir1(i)
403 s = dir1(i+nel)
404 dir2(i) = r*csp-s*snp
405 dir2(i+nel) = s*csp+r*snp
406 ENDDO
407C---
408 DO i=lft,llt
409 IF (iorthloc(i) /= 0) THEN
410 dir1(i) = coor1(ipt,i) ! to check
411 dir1(i+nel) = coor2(ipt,i)
412 dir2(i) = coor3(ipt,i)
413 dir2(i+nel) = coor4(ipt,i)
414 ENDIF
415 ENDDO
416C---
417 DO i=lft,llt
418 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
419 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
420 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
421 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
422 det = u1x*u2y-u1y*u2x
423 w1x = u2y/det
424 w2y = u1x/det
425 w1y = -u1y/det
426 w2x = -u2x/det
427C--- dir I
428 d11 = dir1(i)
429 d21 = dir1(i+nel)
430 dir1(i) = w1x*d11 + w2x*d21
431 dir1(i+nel) = w1y*d11 + w2y*d21
432 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
433 dir1(i) = dir1(i)/s
434 dir1(i+nel) = dir1(i+nel)/s
435C--- dir II
436 d12 = dir2(i)
437 d22= dir2(i+nel)
438 dir2(i) = w1x*d12 + w2x*d22
439 dir2(i+nel) = w1y*d12 + w2y*d22
440 s = sqrt(dir2(i)**2 + dir2(i+nel)**2)
441 dir2(i) = dir2(i)/s
442 dir2(i+nel) = dir2(i+nel)/s
443 ENDDO
444 ENDDO ! NPTT
445 ipt_all = ipt_all + nptt
446 ELSE ! IREP = 0 within PID51
447 DO it = 1, nptt
448 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
449 dir1 => lbuf_dir%DIRA
450 dir2 => lbuf_dir%DIRB
451 ipt = ipt_all + it
452 DO i=lft,llt
453 csp=cos(phi1(ipt,i))
454 snp=sin(phi1(ipt,i))
455 dir1(i) = vr(i)*csp-vs(i)*snp
456 dir1(i+nel) = vs(i)*csp+vr(i)*snp
457 IF (iorthloc(i) /= 0) THEN
458 dir1(i) = coor1(ipt,i)
459 dir1(i+nel) = coor2(ipt,i)
460 ENDIF
461 ENDDO
462 ENDDO
463 ipt_all = ipt_all + nptt
464 ENDIF ! IF (ILAW == 58) THEN
465 ENDDO ! DO N=1,NLAY
466 ELSEIF (irep == 4) THEN
467C mixing law58 with other user laws with IREP = 1 within PID51
468 ipt_all = 0
469 DO il=1,nlay
470 ilaw = elbuf_str%BUFLY(il)%ILAW
471 IF (ilaw == 58) THEN
472 DO it = 1, nptt
473 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
474 dir1 => lbuf_dir%DIRA
475 dir2 => lbuf_dir%DIRB
476 ipt = ipt_all + it
477C--- Axe I d'anisotropie
478 DO i=lft,llt
479 csp=cos(phi1(ipt,i))
480 snp=sin(phi1(ipt,i))
481 dir1(i) = vr(i)*csp-vs(i)*snp
482 dir1(i+nel) = vs(i)*csp+vr(i)*snp
483 ENDDO
484C--- Axe II d'anisotropie
485 DO i=lft,llt
486 csp=cos(phi2(ipt,i))
487 snp=sin(phi2(ipt,i))
488 r = dir1(i)
489 s = dir1(i+nel)
490 dir2(i) = r*csp-s*snp
491 dir2(i+nel) = s*csp+r*snp
492 ENDDO
493C---
494 DO i=lft,llt
495 IF (iorthloc(i) /= 0) THEN
496 dir1(i) = coor1(ipt,i)
497 dir1(i+nel) = coor2(ipt,i)
498 dir2(i) = coor3(ipt,i)
499 dir2(i+nel) = coor4(ipt,i)
500 ENDIF
501 ENDDO
502C---
503 DO i=lft,llt
504 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
505 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
506 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
507 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
508 det = u1x*u2y-u1y*u2x
509 w1x = u2y/det
510 w2y = u1x/det
511 w1y = -u1y/det
512 w2x = -u2x/det
513C--- dir I
514 d11 = dir1(i)
515 d21 = dir1(i+nel)
516 dir1(i) = w1x*d11 + w2x*d21
517 dir1(i+nel) = w1y*d11 + w2y*d21
518 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
519 dir1(i) = dir1(i)/s
520 dir1(i+nel) = dir1(i+nel)/s
521C--- dir II
522 d12 = dir2(i)
523 d22= dir2(i+nel)
524 dir2(i) = w1x*d12 + w2x*d22
525 dir2(i+nel) = w1y*d12 + w2y*d22
526 s = sqrt(dir2(i)**2 + dir2(i+nel)**2)
527 dir2(i) = dir2(i)/s
528 dir2(i+nel) = dir2(i+nel)/s
529 ENDDO
530 ENDDO ! NPTT
531 ipt_all = ipt_all + nptt
532 ELSE ! Ilaw
533 DO it = 1, nptt
534 lbuf_dir => elbuf_str%BUFLY(il)%LBUF_DIR(it)
535 dir1 => lbuf_dir%DIRA
536 dir2 => lbuf_dir%DIRB
537 ipt = ipt_all + it
538 DO i=lft,llt
539 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
540 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
541 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
542 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
543 det = u1x*u2y-u1y*u2x
544 w1x = u2y/det
545 w2y = u1x/det
546 w1y = -u1y/det
547 w2x = -u2x/det
548 IF (iorthloc(i) /= 0) THEN
549 dir1(i) = coor1(ipt,i)
550 dir1(i+nel) = coor2(ipt,i)
551 ENDIF
552 d1 = dir1(i)
553 d2 = dir1(i+nel)
554 dir1(i) = w1x*d1 + w2x*d2
555 dir1(i+nel) = w1y*d1 + w2y*d2
556 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
557 dir1(i) = dir1(i)/s
558 dir1(i+nel) = dir1(i+nel)/s
559 ENDDO
560 ENDDO
561 ipt_all = ipt_all + nptt
562 ENDIF ! IF (ILAW == 58) THEN
563 ENDDO ! DO N=1,NLAY
564 ENDIF ! IF (IREP == 1)
565 ELSE ! idrape ==0
566 DO il=1,nlay
567 dir1 => elbuf_str%BUFLY(il)%DIRA
568 ipid_ply = stack%IGEO(2+il,isubstack)
569 DO i=lft,llt
570 csp=cos(phi1(il,i))
571 snp=sin(phi1(il,i))
572 IF(ipid_ply > 0)THEN
573 dir1(i) = vr(i)*csp-vs(i)*snp
574 dir1(i+nel) = vs(i)*csp+vr(i)*snp
575 IF (iorthloc(i) /= 0) THEN
576 dir1(i) = coor1(il,i)
577 dir1(i+nel) = coor2(il,i)
578 ENDIF
579 ENDIF
580 ENDDO
581 ENDDO
582C---
583 IF (irep == 1) THEN
584 DO i=lft,llt
585 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
586 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
587 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
588 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
589 det = u1x*u2y-u1y*u2x
590 w1x = u2y/det
591 w2y = u1x/det
592 w1y = -u1y/det
593 w2x = -u2x/det
594 DO il=1,nlay
595 dir1 => elbuf_str%BUFLY(il)%DIRA
596 IF (iorthloc(i) /= 0) THEN
597 dir1(i) = coor1(il,i)
598 dir1(i+nel) = coor2(il,i)
599 ENDIF
600 d1 = dir1(i)
601 d2 = dir1(i+nel)
602 dir1(i) = w1x*d1 + w2x*d2
603 dir1(i+nel) = w1y*d1 + w2y*d2
604 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
605 dir1(i) = dir1(i)/s
606 dir1(i+nel) = dir1(i+nel)/s
607 ENDDO
608 ENDDO
609 ELSEIF (irep == 2) THEN
610C--- Axe I d'anisotropie
611 DO il=1,nlay
612 dir1 => elbuf_str%BUFLY(il)%DIRA
613 dir2 => elbuf_str%BUFLY(il)%DIRB
614C--- Axe I d'anisotropie
615 DO i=lft,llt
616 csp=cos(phi1(il,i))
617 snp=sin(phi1(il,i))
618 dir1(i) = vr(i)*csp-vs(i)*snp
619 dir1(i+nel) = vs(i)*csp+vr(i)*snp
620 ENDDO
621C--- Axe II d'anisotropie
622 DO i=lft,llt
623 csp=cos(phi2(il,i))
624 snp=sin(phi2(il,i))
625 r = dir1(i)
626 s = dir1(i+nel)
627 dir2(i) = r*csp-s*snp
628 dir2(i+nel) = s*csp+r*snp
629 ENDDO
630 ENDDO ! DO IL=1,NLAY
631C---
632 DO il=1,nlay
633 dir1 => elbuf_str%BUFLY(il)%DIRA
634 dir2 => elbuf_str%BUFLY(il)%DIRB
635 DO i=lft,llt
636 IF (iorthloc(i) /= 0) THEN
637 dir1(i) = coor1(il,i)
638 dir1(i+nel) = coor2(il,i)
639 dir2(i) = coor3(il,i)
640 dir2(i+nel) = coor4(il,i)
641 ENDIF
642 ENDDO
643 ENDDO ! DO IL=1,NLAY
644C---
645 DO i=lft,llt
646 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
647 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
648 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
649 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
650 det = u1x*u2y-u1y*u2x
651 w1x = u2y/det
652 w2y = u1x/det
653 w1y = -u1y/det
654 w2x = -u2x/det
655 DO il=1,nlay
656 dir1 => elbuf_str%BUFLY(il)%DIRA
657 dir2 => elbuf_str%BUFLY(il)%DIRB
658C--- dir I
659 d11 = dir1(i)
660 d21 = dir1(i+nel)
661 dir1(i) = w1x*d11 + w2x*d21
662 dir1(i+nel) = w1y*d11 + w2y*d21
663 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
664 dir1(i) = dir1(i)/s
665 dir1(i+nel) = dir1(i+nel)/s
666C--- dir II
667 d12 = dir2(i)
668 d22= dir2(i+nel)
669 dir2(i) = w1x*d12 + w2x*d22
670 dir2(i+nel) = w1y*d12 + w2y*d22
671 s = sqrt(dir2(i)**2 + dir2(i+nel)**2)
672 dir2(i) = dir2(i)/s
673 dir2(i+nel) = dir2(i+nel)/s
674 ENDDO
675 ENDDO
676 ELSEIF (irep == 3) THEN
677C mixing law58 with other user laws with IREP = 0 within PID51
678 DO il=1,nlay
679 ilaw = elbuf_str%BUFLY(il)%ILAW
680 IF (ilaw == 58) THEN
681 dir1 => elbuf_str%BUFLY(il)%DIRA
682 dir2 => elbuf_str%BUFLY(il)%DIRB
683C--- Axe I d'anisotropie
684 DO i=lft,llt
685 csp=cos(phi1(il,i))
686 snp=sin(phi1(il,i))
687 dir1(i) = vr(i)*csp-vs(i)*snp
688 dir1(i+nel) = vs(i)*csp+vr(i)*snp
689 ENDDO
690C--- Axe II d'anisotropie
691 DO i=lft,llt
692 csp=cos(phi2(il,i))
693 snp=sin(phi2(il,i))
694 r = dir1(i)
695 s = dir1(i+nel)
696 dir2(i) = r*csp-s*snp
697 dir2(i+nel) = s*csp+r*snp
698 ENDDO
699C---
700 DO i=lft,llt
701 IF (iorthloc(i) /= 0) THEN
702 dir1(i) = coor1(il,i)
703 dir1(i+nel) = coor2(il,i)
704 dir2(i) = coor3(il,i)
705 dir2(i+nel) = coor4(il,i)
706 ENDIF
707 ENDDO
708C---
709 DO i=lft,llt
710 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
711 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
712 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
713 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
714 det = u1x*u2y-u1y*u2x
715 w1x = u2y/det
716 w2y = u1x/det
717 w1y = -u1y/det
718 w2x = -u2x/det
719C--- dir I
720 d11 = dir1(i)
721 d21 = dir1(i+nel)
722 dir1(i) = w1x*d11 + w2x*d21
723 dir1(i+nel) = w1y*d11 + w2y*d21
724 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
725 dir1(i) = dir1(i)/s
726 dir1(i+nel) = dir1(i+nel)/s
727C--- dir II
728 d12 = dir2(i)
729 d22= dir2(i+nel)
730 dir2(i) = w1x*d12 + w2x*d22
731 dir2(i+nel) = w1y*d12 + w2y*d22
732 s = sqrt(dir2(i)**2 + dir2(i+nel)**2)
733 dir2(i) = dir2(i)/s
734 dir2(i+nel) = dir2(i+nel)/s
735 ENDDO
736 ELSE ! IREP = 0 within PID51
737 dir1 => elbuf_str%BUFLY(il)%DIRA
738 DO i=lft,llt
739 csp=cos(phi1(il,i))
740 snp=sin(phi1(il,i))
741 dir1(i) = vr(i)*csp-vs(i)*snp
742 dir1(i+nel) = vs(i)*csp+vr(i)*snp
743 IF (iorthloc(i) /= 0) THEN
744 dir1(i) = coor1(il,i)
745 dir1(i+nel) = coor2(il,i)
746 ENDIF
747 ENDDO
748 ENDIF ! IF (ILAW == 58) THEN
749 ENDDO ! DO N=1,NLAY
750 ELSEIF (irep == 4) THEN
751C mixing law58 with other user laws with IREP = 1 within PID51
752 DO il=1,nlay
753 ilaw = elbuf_str%BUFLY(il)%ILAW
754 IF (ilaw == 58) THEN
755 dir1 => elbuf_str%BUFLY(il)%DIRA
756 dir2 => elbuf_str%BUFLY(il)%DIRB
757C--- Axe I d'anisotropie
758 DO i=lft,llt
759 csp=cos(phi1(il,i))
760 snp=sin(phi1(il,i))
761 dir1(i) = vr(i)*csp-vs(i)*snp
762 dir1(i+nel) = vs(i)*csp+vr(i)*snp
763 ENDDO
764C--- Axe II d'anisotropie
765 DO i=lft,llt
766 csp=cos(phi2(il,i))
767 snp=sin(phi2(il,i))
768 r = dir1(i)
769 s = dir1(i+nel)
770 dir2(i) = r*csp-s*snp
771 dir2(i+nel) = s*csp+r*snp
772 ENDDO
773C---
774 DO i=lft,llt
775 IF (iorthloc(i) /= 0) THEN
776 dir1(i) = coor1(il,i)
777 dir1(i+nel) = coor2(il,i)
778 dir2(i) = coor3(il,i)
779 dir2(i+nel) = coor4(il,i)
780 ENDIF
781 ENDDO
782C---
783 DO i=lft,llt
784 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
785 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
786 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
787 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
788 det = u1x*u2y-u1y*u2x
789 w1x = u2y/det
790 w2y = u1x/det
791 w1y = -u1y/det
792 w2x = -u2x/det
793C--- dir I
794 d11 = dir1(i)
795 d21 = dir1(i+nel)
796 dir1(i) = w1x*d11 + w2x*d21
797 dir1(i+nel) = w1y*d11 + w2y*d21
798 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
799 dir1(i) = dir1(i)/s
800 dir1(i+nel) = dir1(i+nel)/s
801C--- dir II
802 d12 = dir2(i)
803 d22= dir2(i+nel)
804 dir2(i) = w1x*d12 + w2x*d22
805 dir2(i+nel) = w1y*d12 + w2y*d22
806 s = sqrt(dir2(i)**2 + dir2(i+nel)**2)
807 dir2(i) = dir2(i)/s
808 dir2(i+nel) = dir2(i+nel)/s
809 ENDDO
810 ELSE ! IREP = 1 within PID51
811 dir1 => elbuf_str%BUFLY(il)%DIRA
812 DO i=lft,llt
813 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
814 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
815 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
816 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
817 det = u1x*u2y-u1y*u2x
818 w1x = u2y/det
819 w2y = u1x/det
820 w1y = -u1y/det
821 w2x = -u2x/det
822 IF (iorthloc(i) /= 0) THEN
823 dir1(i) = coor1(il,i)
824 dir1(i+nel) = coor2(il,i)
825 ENDIF
826 d1 = dir1(i)
827 d2 = dir1(i+nel)
828 dir1(i) = w1x*d1 + w2x*d2
829 dir1(i+nel) = w1y*d1 + w2y*d2
830 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
831 dir1(i) = dir1(i)/s
832 dir1(i+nel) = dir1(i+nel)/s
833 ENDDO
834 ENDIF ! IF (ILAW == 58) THEN
835 ENDDO ! DO N=1,NLAY
836 ENDIF ! IF (IREP == 1)
837 ENDIF ! IDRAPE
838C-------
839 ELSEIF (igtyp == 16) THEN
840C-------
841C--- Axe I d'anisotropie
842 DO il=1,nlay
843 dir1 => elbuf_str%BUFLY(il)%DIRA
844 DO i=lft,llt
845 csp=cos(phi1(il,i))
846 snp=sin(phi1(il,i))
847 dir1(i) = vr(i)*csp-vs(i)*snp
848 dir1(i+nel) = vs(i)*csp+vr(i)*snp
849 ENDDO
850 ENDDO
851C--- Axe II d'anisotropie
852 DO il=1,nlay
853 dir1 => elbuf_str%BUFLY(il)%DIRA
854 dir2 => elbuf_str%BUFLY(il)%DIRB
855 DO i=lft,llt
856 csp=cos(phi2(il,i))
857 snp=sin(phi2(il,i))
858 r = dir1(i)
859 s = dir1(i+nel)
860 dir2(i) = r*csp-s*snp
861 dir2(i+nel) = s*csp+r*snp
862 ENDDO
863 ENDDO
864 DO il=1,nlay
865 dir1 => elbuf_str%BUFLY(il)%DIRA
866 dir2 => elbuf_str%BUFLY(il)%DIRB
867 DO i=lft,llt
868 IF (iorthloc(i) /= 0) THEN
869 dir1(i) = coor1(il,i)
870 dir1(i+nel) = coor2(il,i)
871 dir2(i) = coor3(il,i)
872 dir2(i+nel) = coor4(il,i)
873 ENDIF
874 ENDDO
875 ENDDO
876C------
877 DO i=lft,llt
878 u1x = e11(i)*e1x(i)+e12(i)*e1y(i)+e13(i)*e1z(i)
879 u1y = e11(i)*e2x(i)+e12(i)*e2y(i)+e13(i)*e2z(i)
880 u2x = e21(i)*e1x(i)+e22(i)*e1y(i)+e23(i)*e1z(i)
881 u2y = e21(i)*e2x(i)+e22(i)*e2y(i)+e23(i)*e2z(i)
882 det = u1x*u2y-u1y*u2x
883 w1x = u2y/det
884 w2y = u1x/det
885 w1y = -u1y/det
886 w2x = -u2x/det
887 DO il=1,nlay
888 dir1 => elbuf_str%BUFLY(il)%DIRA
889 dir2 => elbuf_str%BUFLY(il)%DIRB
890C--- dir I
891 d11 = dir1(i)
892 d21 = dir1(i+nel)
893 dir1(i) = w1x*d11 + w2x*d21
894 dir1(i+nel) = w1y*d11 + w2y*d21
895 s = sqrt(dir1(i)**2 + dir1(i+nel)**2)
896 dir1(i) = dir1(i)/s
897 dir1(i+nel) = dir1(i+nel)/s
898C--- dir II
899 d12 = dir2(i)
900 d22= dir2(i+nel)
901 dir2(i) = w1x*d12 + w2x*d22
902 dir2(i+nel) = w1y*d12 + w2y*d22
903 s = sqrt(dir2(i)**2 + dir2(i+nel)**2)
904 dir2(i) = dir2(i)/s
905 dir2(i+nel) = dir2(i+nel)/s
906 ENDDO
907 ENDDO
908C----
909 ENDIF
910C----
911 ENDIF ! IF (IGTYP == 1)
912C-----------
913 RETURN
914 END
corthdir
subroutine corthdir(elbuf_str, igeo, geo, vx, vy, vz, phi1, phi2, coor1, coor2, coor3, coor4, iorthloc, nlay, irep, isubstack, stack, geo_stack, igeo_stack, ir, is, nel, imat, iprop, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z, npt_all, idrape)
Definition corthdir.F:42
stack_mod
Definition stack_mod.F:72
starter
program starter
Definition starter.F:39