OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
corthdir.F File Reference

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ corthdir()

subroutine corthdir ( type(elbuf_struct_), target elbuf_str,
integer, dimension(npropgi,*) igeo,
geo,
vx,
vy,
vz,
phi1,
phi2,
coor1,
coor2,
coor3,
coor4,
integer, dimension(*) iorthloc,
integer nlay,
integer irep,
integer isubstack,
type (stack_ply) stack,
geo_stack,
integer, dimension(npropgi,*) igeo_stack,
integer ir,
integer is,
integer nel,
integer imat,
integer iprop,
intent(in) x1,
intent(in) x2,
intent(in) x3,
intent(in) x4,
intent(in) y1,
intent(in) y2,
intent(in) y3,
intent(in) y4,
intent(in) z1,
intent(in) z2,
intent(in) z3,
intent(in) z4,
intent(in) e1x,
intent(in) e2x,
intent(in) e3x,
intent(in) e1y,
intent(in) e2y,
intent(in) e3y,
intent(in) e1z,
intent(in) e2z,
intent(in) e3z,
integer npt_all,
integer idrape )

Definition at line 32 of file corthdir.F.

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
my_real
#define my_real
Definition cppsort.cpp:32
stack_mod
Definition stack_mod.F:72