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

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