Télécharger idmat4.eso

Retour à la liste

Numérotation des lignes :

  1. C IDMAT4 SOURCE BP208322 17/03/01 21:17:40 9325
  2. SUBROUTINE IDMAT4 (NUMP1,NUMP2,NUDIR1,NUDIR2,ANG,ANG2,
  3. . MELEME,MINTE,IPVAL,NPG2,MINTE2)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8(A-H,O-Z)
  6. *
  7. ************************************************************************
  8. *
  9. * I D M A T 4
  10. * -----------
  11. *
  12. * FONCTION:
  13. * ---------
  14. *
  15. * CALCUL DES COS-DIR. DES AXES ORTH. (OU ANISO.)PAR RAPPORT AU
  16. *
  17. * REPERE LOCAL D'ELEMENT , A PARTIR D'UNE DEFINITION GLOBALE
  18. *
  19. *
  20. * MODULES UTILISES:
  21. * -----------------
  22. *
  23. -INC SMCOORD
  24.  
  25. -INC PPARAM
  26. -INC CCOPTIO
  27. -INC CCREEL
  28. -INC SMELEME
  29. -INC SMINTE
  30. *
  31. * PARAMETRES: (E)=ENTREE (S)=SORTIE (+ = CONTENU DANS UN COMMUN)
  32. * -----------
  33. *
  34. * NUMP1 (E) NUMERO DU POINT P1 ASSOCIE A NUDIR1
  35. * NUMP2 (E) NUMERO DU POINT P2 ASSOCIE A NUDIR1
  36. * NUDIR1 (E) NUMERO DE DIRECTIVE UTILISEE DANS LA LISTE:
  37. * "DIRECTION", "RADIAL"
  38. * NUDIR2 (E) NUMERO DE DIRECTIVE UTILISEE DANS LA LISTE:
  39. * "PARALLELE", "PERPENDIC.", "INCLINE",
  40. * POUR LA DEFINITION DES DIRECTIONS D'ORTHOTROPIE.
  41. * ANG ANG2 (E) ANGLEs UTILISEs DANS LA DEFINITION DES DIRECTIONS
  42. * D ORTHOTROPIE (INCLINE)
  43. * MELEME (E) POINTEUR DE "MAILLAGE" A 1 SEUL TYPE D'ELEMENT.
  44. * MINTE (E) SEGMENT CONTENANT LES FONCTIONS DE FORME,ET LEURS
  45. * DERIVEES
  46. * XVAL (S) DIRECTIONS D'ORTHOTROPIE PAR ELEMENT,
  47. * ON FOURNIT LES COSINUS DIRECTEURS DE L'AXE 1
  48. * EN BIDIM ,ET LES COSINUS DIRECTEURS DES AXES 1 ET
  49. * 2 EN TRIDIM,PAR RAPPORT AU REPERE LOCAL DE L'ELEMENT
  50. *
  51. *
  52. * VARIABLES:
  53. *-----------
  54. * VGLOB1,VGLOB2,VGLOB3 = COS.DIRECTEURS DES AXES 1, 2 ET 3 D'ORTH. PAR
  55. * RAPPOT AU REPERE GLOBAL
  56. * TXR = TABLEAU CONTENANT LES COS.DIRECTEURS DES AXES LOCAUX PAR
  57. * RAPPORT AU REPERE GLOBAL
  58. *
  59. *
  60. * REMARQUES:
  61. * ----------
  62. *
  63. * LES DIRECTION P1 ET P2 DEFINISSENT LE PLAN QUI CONTIENT LES AXES 1,2
  64. * NUDIR2 =1 SIGNIFIE QUE L'AXE 1 EST PARALLELE A LA DIRECTION P1,
  65. * NUDIR2 =2 SIGNIFIE QUE L'AXE 1 EST PERPENDICULAIRE A LA DIRECTION P1,
  66. * NUDIR2 =3 SIGNIFIE QUE L'AXE 1 FAIT AN ANGLE DE (ANG) AVEC LA
  67. * DIRECTION P1
  68. *
  69. *
  70. * AUTEUR, DATE DE CREATION:
  71. * -------------------------
  72. *
  73. * P. DOWLATYARI OCT. 1990
  74. *
  75. * LANGAGE:
  76. * --------
  77. * FORTRAN 77 + ESOPE
  78. *
  79. *
  80. ************************************************************************
  81. *
  82. *
  83. DIMENSION VGLOB1(3),VGLOB2(3),VGLOB3(3),VAUX1(3),VAUX2(3)
  84. DIMENSION VPT1(3),VPT2(3)
  85. *
  86. SEGMENT YVAL
  87. REAL*8 VLOC1(IDIM2,NPG2,NBELEM),VLOC2(IDIM2,NPG2,NBELEM)
  88. ENDSEGMENT
  89. *
  90. SEGMENT WORK
  91. REAL*8 XE(3,NBNN),TXR(IDIM,IDIM)
  92. ENDSEGMENT
  93. *
  94. SEGACT MELEME
  95. NBELEM=NUM(/2)
  96. NBNN=NUM(/1)
  97. SEGACT MINTE
  98. IF(NUDIR1.EQ.2) SEGACT MINTE2
  99.  
  100.  
  101. ************************************************************************
  102. * CALCUL DE VGLOB1
  103. ************************************************************************
  104.  
  105. VGLOB1(3)=0.D0
  106. VGLOB2(3)=0.D0
  107.  
  108. ********** CAS 2D *************
  109.  
  110. IF (IDIM.EQ.2)THEN
  111. *
  112. * COORDONNEES DU POINT P1
  113. CALL EXCOO1 (NUMP1,VPT1(1),VPT1(2),VPT1(3),REEL1)
  114. *
  115. * OPTION 'RADIAL' : VPT1 = CENTRE -> ON STOCKE PUIS GOTO 15
  116. IF(NUDIR1.EQ.2) THEN
  117. VAUX1(1)=VPT1(1)
  118. VAUX1(2)=VPT1(2)
  119. GO TO 15
  120. ENDIF
  121. *
  122. * OPTION 'DIRECTION' : VPT1 = DIRECTION 1
  123. * NORMALISATION DU VECTEUR FOURNI
  124. c + CALCUL DE V2 = ORTHOGONAL A V1 DANS LE PLAN 2D
  125. VNORM=SQRT(VPT1(1)*VPT1(1)+VPT1(2)*VPT1(2))
  126. IF (VNORM.EQ.0.) THEN
  127. CALL ERREUR (524)
  128. SEGDES,MINTE,MELEME
  129. RETURN
  130. ENDIF
  131. VGLOB1(1)=VPT1(1)/VNORM
  132. VGLOB1(2)=VPT1(2)/VNORM
  133. VGLOB2(1)=-VGLOB1(2)
  134. VGLOB2(2)=VGLOB1(1)
  135. *
  136. * OPTION PERPENDICULAIRE :
  137. * ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR DE L'AXE 3
  138. IF(NUDIR2.EQ.2)THEN
  139. VGLOB1(1)=VGLOB2(1)
  140. VGLOB1(2)=VGLOB2(2)
  141. VGLOB2(1)=-VGLOB1(2)
  142. VGLOB2(2)=VGLOB1(1)
  143. *
  144. * OPTION INCLINE :
  145. * ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  146. * -> \tilde{[V1,V2,V3]}=[eR eZ eT]*[cosa sina 0 ; -sina cosa 0 ; 0 0 1]
  147. * + eventuellement une rotation de ANG2 autour de \tilde{V1}
  148. * -> [V1,V2,V3]=\tilde{[V1,V2,V3]}*[1 0 0; 0 cosb -sinb; 0 sinb cosb]
  149. ELSEIF(NUDIR2.EQ.3)THEN
  150. COSA=COS(ANG)
  151. SINA=SIN(ANG)
  152. VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  153. VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  154. VGLOB2(1)=-VGLOB1(2)
  155. VGLOB2(2)=VGLOB1(1)
  156. IF(ABS(ANG2).GT.XZPREC) THEN
  157. c VGLOB3=(0 0 1)^T, VGLOB2=(XX1 XX2 0)^T
  158. c VGLOB1 ne change pas, on tourne VGLOB2 (et VGLOB3 implicitement)
  159. VGLOB2(1)=VGLOB2(1)*COS(ANG2)
  160. VGLOB2(2)=VGLOB2(2)*COS(ANG2)
  161. VGLOB2(3)=SIN(ANG2)
  162. ENDIF
  163.  
  164. ENDIF
  165.  
  166.  
  167. ********** CAS 3D *************
  168.  
  169. ELSEIF(IDIM.EQ.3)THEN
  170. *
  171. * COORDONNEES DU POINT P1
  172. CALL EXCOO1 (NUMP1,VPT1(1),VPT1(2),VPT1(3),REEL1)
  173. *
  174. * OPTION 'DIRECTION' : VPT1 = DIRECTION 1
  175. * NORMALISATION
  176. IF(NUDIR1.NE.2) THEN
  177. VNORM=SQRT(VPT1(1)*VPT1(1)+VPT1(2)*VPT1(2)+VPT1(3)*VPT1(3))
  178. IF (VNORM.EQ.0.) THEN
  179. CALL ERREUR (524)
  180. SEGDES,MINTE,MELEME
  181. RETURN
  182. ENDIF
  183. VGLOB1(1)=VPT1(1)/VNORM
  184. VGLOB1(2)=VPT1(2)/VNORM
  185. VGLOB1(3)=VPT1(3)/VNORM
  186. ENDIF
  187. *
  188. * COORDONNEES DU POINT P2
  189. CALL EXCOO1 (NUMP2,VPT2(1),VPT2(2),VPT2(3),REEL1)
  190. *
  191. * OPTION 'RADIAL' : DIRECTION 1 = P1P2
  192. * NORMALISATION PUIS GOTO 15
  193. IF(NUDIR1.EQ.2) THEN
  194. VGLOB1(1)=VPT2(1)-VPT1(1)
  195. VGLOB1(2)=VPT2(2)-VPT1(2)
  196. VGLOB1(3)=VPT2(3)-VPT1(3)
  197. VNORM=SQRT(VGLOB1(1)*VGLOB1(1)+VGLOB1(2)*VGLOB1(2)+
  198. . VGLOB1(3)*VGLOB1(3))
  199. VGLOB1(1)=VGLOB1(1)/VNORM
  200. VGLOB1(2)=VGLOB1(2)/VNORM
  201. VGLOB1(3)=VGLOB1(3)/VNORM
  202. GO TO 15
  203. ENDIF
  204. *
  205. * OPTION 'DIRECTION' :
  206. * DIRECTION 3 = PRODUIT VECTORIEL DES DEUX VECTEURS + NORMALISATION
  207. CALL CROSS2 (VGLOB1(1),VPT2(1),VGLOB3(1),IRR)
  208. IF (IRR.EQ.0)THEN
  209. * DIRECTIONS D'ORTHOTROPIE MAL FOURNIES
  210. CALL ERREUR (524)
  211. SEGDES,MINTE,MELEME
  212. IF(NUDIR1.EQ.2) SEGDES MINTE2
  213. RETURN
  214. ENDIF
  215. * DIRECTION 2 = PRODUIT VECTORIEL DES DEUX VECTEURS + NORMALISATION
  216. CALL CROSS2 (VGLOB3(1),VGLOB1(1),VGLOB2(1),IRR)
  217. *
  218. * OPTION 'PERPENDICULAIRE'
  219. * ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR DE L'AXE 3
  220. IF(NUDIR2.EQ.2)THEN
  221. REEL1=VGLOB1(1)
  222. VGLOB1(1)=VGLOB2(1)
  223. VGLOB2(1)=-REEL1
  224. REEL1=VGLOB1(2)
  225. VGLOB1(2)=VGLOB2(2)
  226. VGLOB2(2)=-REEL1
  227. REEL1=VGLOB1(3)
  228. VGLOB1(3)=VGLOB2(3)
  229. VGLOB2(3)=-REEL1
  230. *
  231. * OPTION 'INCLINE'
  232. * ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  233. ELSEIF(NUDIR2.EQ.3)THEN
  234. COSA=COS(ANG)
  235. SINA=SIN(ANG)
  236. REEL1=VGLOB1(1)
  237. VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  238. VGLOB2(1)=-REEL1*SINA+VGLOB2(1)*COSA
  239. REEL1=VGLOB1(2)
  240. VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  241. VGLOB2(2)=-REEL1*SINA+VGLOB2(2)*COSA
  242. REEL1=VGLOB1(3)
  243. VGLOB1(3)=VGLOB1(3)*COSA+VGLOB2(3)*SINA
  244. VGLOB2(3)=-REEL1*SINA+VGLOB2(3)*COSA
  245. ENDIF
  246.  
  247. ENDIF
  248. ************************************************************************
  249.  
  250. *
  251. 15 CONTINUE
  252. *
  253. IF(IFOUR.EQ.1)THEN
  254. IDIM2=3
  255. ELSE
  256. IDIM2=IDIM
  257. ENDIF
  258. SEGINI YVAL
  259. IPVAL=YVAL
  260. SEGINI WORK
  261. *
  262. ************************************************************************
  263. * BOUCLE SUR LES ELEMENTS
  264. ************************************************************************
  265. *
  266. DO 10 IEL=1,NBELEM
  267. *
  268. DO 20 NC=1,IDIM2
  269. DO 20 NV=1,NPG2
  270. VLOC1(NC,NV,IEL)=0.D0
  271. VLOC2(NC,NV,IEL)=0.D0
  272. 20 CONTINUE
  273. *
  274. * COORDONNEES DES NOEUDS
  275. *
  276. CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XE)
  277. *
  278. * CALCUL DU REPERE LOCAL DE L'ELEMENT
  279. *
  280. NBSH=SHPTOT(/2)
  281. CALL RLOCAL(XE,SHPTOT,NBSH,NBNN,TXR)
  282. if (nbsh.eq.-1) then
  283. call erreur(525)
  284. return
  285. endif
  286. IF(IERR.NE.0)THEN
  287. SEGSUP,WORK,YVAL
  288. SEGDES , MELEME,MINTE
  289. IF(NUDIR1.EQ.2) SEGDES MINTE2
  290. RETURN
  291. ENDIF
  292. *
  293. ********************************************************
  294. * BOUCLE SUR LES POINTS DE GAUSS
  295. ********************************************************
  296. *
  297. DO 80 IGAU=1,NPG2
  298.  
  299. c---- OPTION 'RADIAL' -----
  300. IF (NUDIR1.EQ.2) THEN
  301.  
  302. DO 25 ICO=1,IDIM
  303. VAUX2(ICO)=0.D0
  304. DO 25 ID=1,NBNN
  305. VAUX2(ICO)=VAUX2(ICO)+
  306. . MINTE2.SHPTOT(1,ID,IGAU)*XE(ICO,ID)
  307. 25 CONTINUE
  308. *
  309. * --CAS 2D--
  310. IF(IDIM.EQ.2) THEN
  311.  
  312. VGLOB1(1)=VAUX1(1)-VAUX2(1)
  313. VGLOB1(2)=VAUX1(2)-VAUX2(2)
  314. * NORMALISATION
  315. VNORM=SQRT(VGLOB1(1)*VGLOB1(1)+VGLOB1(2)*VGLOB1(2))
  316. VGLOB1(1)=VGLOB1(1)/VNORM
  317. VGLOB1(2)=VGLOB1(2)/VNORM
  318. VGLOB2(1)=-VGLOB1(2)
  319. VGLOB2(2)=VGLOB1(1)
  320.  
  321. * OPTION 'PERPENDICULAIRE'
  322. * ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR DE L'AXE 3
  323. IF(NUDIR2.EQ.2)THEN
  324. VGLOB1(1)=VGLOB2(1)
  325. VGLOB1(2)=VGLOB2(2)
  326.  
  327. * OPTION 'INCLINE'
  328. * ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  329. ELSEIF(NUDIR2.EQ.3)THEN
  330. COSA=COS(ANG)
  331. SINA=SIN(ANG)
  332. VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  333. VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  334. ENDIF
  335.  
  336. IF(IFOUR.EQ.1) THEN
  337. VGLOB1(3)=0.D0
  338. VGLOB2(1)=-VGLOB1(2)
  339. VGLOB2(2)=VGLOB1(1)
  340. VGLOB2(3)=0.D0
  341. ENDIF
  342.  
  343. ENDIF
  344. * ----------
  345.  
  346. * --CAS 3D--
  347. IF(IDIM.EQ.3) THEN
  348.  
  349. VAUX1(1)=VAUX2(1)-VPT2(1)
  350. VAUX1(2)=VAUX2(2)-VPT2(2)
  351. VAUX1(3)=VAUX2(3)-VPT2(3)
  352. PSC=VAUX1(1)*VGLOB1(1)+VAUX1(2)*VGLOB1(2)+VAUX1(3)*VGLOB1(3)
  353. VGLOB2(1)=VAUX1(1)-PSC*VGLOB1(1)
  354. VGLOB2(2)=VAUX1(2)-PSC*VGLOB1(2)
  355. VGLOB2(3)=VAUX1(3)-PSC*VGLOB1(3)
  356. * NORMALISATION
  357. VNORM=SQRT(VGLOB2(1)*VGLOB2(1)+VGLOB2(2)*VGLOB2(2)+
  358. . VGLOB2(3)*VGLOB2(3))
  359. VGLOB2(1)=VGLOB2(1)/VNORM
  360. VGLOB2(2)=VGLOB2(2)/VNORM
  361. VGLOB2(3)=VGLOB2(3)/VNORM
  362. *
  363. * PRODUIT VECTORIEL DES DEUX VECTEURS ET NORMALISATION
  364. CALL CROSS2 (VGLOB1(1),VGLOB2(1),VGLOB3(1),IRR)
  365. IF (IRR.EQ.0)THEN
  366. CALL ERREUR (524)
  367. SEGDES,MINTE,MELEME
  368. IF(NUDIR1.EQ.2) SEGDES MINTE2
  369. RETURN
  370. ENDIF
  371. *
  372. * OPTION 'PERPENDICULAIRE'
  373. * ON EFFECTUE UNE ROTATION DE 90 DEGRE AUTOUR DE L'AXE 3
  374. IF(NUDIR2.EQ.2)THEN
  375. REEL1=VGLOB1(1)
  376. VGLOB1(1)=VGLOB2(1)
  377. VGLOB2(1)=-REEL1
  378. REEL1=VGLOB1(2)
  379. VGLOB1(2)=VGLOB2(2)
  380. VGLOB2(2)=-REEL1
  381. REEL1=VGLOB1(3)
  382. VGLOB1(3)=VGLOB2(3)
  383. VGLOB2(3)=-REEL1
  384. *
  385. * OPTION 'INCLINE'
  386. * ON EFFECTUE UNE ROTATION DE (ANG) DEGRE AUTOUR DE L'AXE 3
  387. ELSEIF(NUDIR2.EQ.3)THEN
  388. COSA=COS(ANG)
  389. SINA=SIN(ANG)
  390. REEL1=VGLOB1(1)
  391. VGLOB1(1)=VGLOB1(1)*COSA+VGLOB2(1)*SINA
  392. VGLOB2(1)=-REEL1*SINA+VGLOB2(1)*COSA
  393. REEL1=VGLOB1(2)
  394. VGLOB1(2)=VGLOB1(2)*COSA+VGLOB2(2)*SINA
  395. VGLOB2(2)=-REEL1*SINA+VGLOB2(2)*COSA
  396. REEL1=VGLOB1(3)
  397. VGLOB1(3)=VGLOB1(3)*COSA+VGLOB2(3)*SINA
  398. VGLOB2(3)=-REEL1*SINA+VGLOB2(3)*COSA
  399. ENDIF
  400.  
  401. ENDIF
  402. * ----------
  403.  
  404. ENDIF
  405. c---- FIN DE L'OPTION 'RADIAL' -----
  406. *
  407. * CALCUL DES COS.DIRECTEURS DU REPERE LOCAL / REPERE GLOBAL
  408. *
  409. IF(IDIM.EQ.2)THEN
  410. VLOC1(1,IGAU,IEL)=TXR(1,1)*VGLOB1(1)+TXR(2,1)*VGLOB1(2)
  411. VLOC1(2,IGAU,IEL)=TXR(1,2)*VGLOB1(1)+TXR(2,2)*VGLOB1(2)
  412. IF(IFOUR.EQ.1)THEN
  413. c bp: en 2D Fourier, TXR= [ TXR(1,1) TXR(1,2) 0 ] dans (R,Z,T)
  414. c [ TXR(2,1) TXR(2,2) 0 ]
  415. c [ 0 0 1 ]
  416. VLOC1(3,IGAU,IEL)=VGLOB1(3)
  417. VLOC2(1,IGAU,IEL)=TXR(1,1)*VGLOB2(1)+TXR(2,1)*VGLOB2(2)
  418. VLOC2(2,IGAU,IEL)=TXR(1,2)*VGLOB2(1)+TXR(2,2)*VGLOB2(2)
  419. VLOC2(3,IGAU,IEL)=VGLOB2(3)
  420. ENDIF
  421. *
  422. ELSEIF(IDIM.EQ.3)THEN
  423. DO 40 J=1,3
  424. DO 40 I=1,3
  425. VLOC1(J,IGAU,IEL)=VLOC1(J,IGAU,IEL)+VGLOB1(I)*TXR(I,J)
  426. VLOC2(J,IGAU,IEL)=VLOC2(J,IGAU,IEL)+VGLOB2(I)*TXR(I,J)
  427. 40 CONTINUE
  428. ENDIF
  429.  
  430. 80 CONTINUE
  431. ********************************************************
  432. * FIN DE BOUCLE SUR LES POINTS DE GAUSS
  433. ********************************************************
  434.  
  435. 10 CONTINUE
  436. ************************************************************************
  437. * FIN DE BOUCLE SUR LES ELEMENTS
  438. ************************************************************************
  439. *
  440. * DESACTIVATION DES SEGMENTS
  441. *
  442. SEGSUP,WORK
  443. SEGDES , MELEME,MINTE,YVAL
  444. IF(NUDIR1.EQ.2) SEGDES MINTE2
  445. *
  446. RETURN
  447. END
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales