Télécharger orient.eso

Retour à la liste

Numérotation des lignes :

  1. C ORIENT SOURCE PV 20/04/14 21:15:04 10584
  2. C CET OPERATEUR ORIENTE LES ELEMENTS ORIENTABLES EN FONCTION DE XP
  3. C
  4. SUBROUTINE ORIENT(IPT1,IPT2,XP,ICLE)
  5. *
  6. C IPT1 (E) MAILLAGE A ORIENTER (segment ACTIF)
  7. C IPT2 (S) MAILLAGE ORIENTE (segment ACTIF)
  8. C IPT2 peut être égal à IPT1 si on ne savait pas orienter ou si le
  9. C maillage
  10. C est deja orienté...
  11. *
  12. * ICLE=1 ORIEN DIRE
  13. * ICLE=2 ORIEN POINT
  14. *
  15. C SG : 2016/05/17 ajout orientation elements massifs
  16. *
  17. *
  18. IMPLICIT INTEGER(I-N)
  19. IMPLICIT REAL*8 (A-H,O-z)
  20.  
  21. -INC CCREEL
  22. -INC SMELEME
  23. -INC SMLENTI
  24. -INC CCGEOME
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC SMCOORD
  29.  
  30. DIMENSION XP(3)
  31. DIMENSION XQ(3)
  32. LOGICAL LQUAF
  33. PARAMETER (NLMASS=22)
  34. * NLNOMA=2*2 + 8*3 + 12*4
  35. PARAMETER (NLNOMA=76)
  36. * Tableau donnant la liste des types d'elements susceptibles d'etre
  37. * traites dans le cas massif
  38. INTEGER LTYMAS(NLMASS)
  39. * Pour chaque element susceptible d'etre traite, ce tableau donne
  40. * l'adresse dans le tableau LNOMAS de la description des noeuds
  41. * servant à calculer l'orientation
  42. INTEGER LADMAS(NLMASS)
  43. * Ce tableau donne les idim+1 noeuds (triedre) de chaque element massif
  44. * servant pour calculer l'orientation
  45. * On aurait pu le simplifier sachant que les QUAD et QUAF sont
  46. * similaires...
  47. INTEGER LNOMAS(NLNOMA)
  48. * SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8 QUA9 10
  49. DATA LTYMAS/ 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
  50. C CUB8 CU20 PRI6 PR15 TET4 TE10 PYR5 PY13 18
  51. $ 14, 15, 16, 17, 23, 24, 25, 26,
  52. C CU27 PR21 TE15 PY19 22
  53. $ 33, 34, 35, 36/
  54. * SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4 QUA5 QUA8 QUA9 10
  55. DATA LADMAS/ 1, 3, 5, 8, 11, 14, 17, 20, 23, 26,
  56. C CUB8 CU20 PRI6 PR15 TET4 TE10 PYR5 PY13 18
  57. $ 29, 33, 37, 41, 45, 49, 53, 57,
  58. C CU27 PR21 TE15 PY19 22
  59. $ 61, 65, 69, 73/
  60. * SEG2 SEG3 TRI3 TRI4 TRI6 TRI7 QUA4
  61. DATA LNOMAS/ 1,2, 1,3, 1,2,3, 1,2,3, 1,3,5, 1,3,5, 1,2,4,
  62. C QUA5 QUA8 QUA9 CUB8 CU20 PRI6
  63. $ 1,2,4, 1,3,7, 1,3,7, 1,2,4,5, 1,3,7,13, 1,2,3,4,
  64. C PR15 TET4 TE10 PYR5 PY13
  65. $ 1,3,5,10, 1,2,3,4, 1,3,5,10, 1,2,4,5, 1,3,7,13,
  66. C CU27 PR21 TE15 PY19
  67. $ 1,3,7,13, 1,3,5,10, 1,3,5,10, 1,3,7,13/
  68. segact mcoord
  69. *
  70. * Cas de l'orientation des éléments massifs
  71. * On souhaite orienter positivement les éléments (ie même
  72. * orientation que l'élément de référence)
  73. *
  74. * Est-on dans le cas des elements massifs ?
  75. ilmass=0
  76. ITY=IPT1.ITYPEL
  77. DO i=1,nlmass
  78. if (ity.eq.ltymas(i)) then
  79. ilmass=i
  80. goto 666
  81. endif
  82. enddo
  83. 666 continue
  84. *dbg write(ioimp,*) 'ilmass=',ilmass
  85. *dbg write(ioimp,*) 'idim=',idim
  86. if (ilmass.eq.0) goto 1000
  87. if (LDLR(ity).ne.idim) goto 1000
  88. *dbg write(ioimp,*) 'Cas massif'
  89. * Si oui, il faut que ICLE=0 sauf pour le cub8 (shb8) qui est traité
  90. * plus loin
  91. if (icle.ne.0) then
  92. if (ity.eq.14) goto 1000
  93. if (icle.eq.1) moterr(1:8)='DIRE '
  94. if (icle.eq.2) moterr(1:8)='POIN '
  95. *dbg write(ioimp,*) 'icle=',icle
  96. * 803 2
  97. * Option %m1:8 incompatible avec les donnees
  98. call erreur(803)
  99. return
  100. endif
  101. * Ici, on teste l'orientation des éléments et on la stocke dans une
  102. * liste d'entiers valant +1 ou -1.
  103. * On stocke le nombre d'orientation négative dans ninv
  104. *
  105. * oubli des references (suivant la logique de invers)
  106. *eff nbref=0
  107. *eff nbsous=0
  108. nbnn=ipt1.num(/1)
  109. nbelem=ipt1.num(/2)
  110. *eff segini meleme
  111. *eff itypel=ipt1.itypel
  112. iadr=ladmas(ilmass)-1
  113. *dbg write(ioimp,*) 'iadr=',iadr
  114. jg=nbelem
  115. segini mlenti
  116. ninv=0
  117. do il=1,nbelem
  118. * test orientation
  119. if (idim.eq.1) then
  120. ia=ipt1.num(lnomas(iadr+1),il)
  121. ib=ipt1.num(lnomas(iadr+2),il)
  122. *dbg write(ioimp,*) 'ia=',ia,' ib=',ib
  123. ira=(idim+1)*(ia-1)
  124. irb=(idim+1)*(ib-1)
  125. x1=xcoor(irb+1)-xcoor(ira+1)
  126. pmix=x1
  127. elseif (idim.eq.2) then
  128. ia=ipt1.num(lnomas(iadr+1),il)
  129. ib=ipt1.num(lnomas(iadr+2),il)
  130. ic=ipt1.num(lnomas(iadr+3),il)
  131. ira=(idim+1)*(ia-1)
  132. irb=(idim+1)*(ib-1)
  133. irc=(idim+1)*(ic-1)
  134. x1=xcoor(irb+1)-xcoor(ira+1)
  135. y1=xcoor(irb+2)-xcoor(ira+2)
  136. x2=xcoor(irc+1)-xcoor(ira+1)
  137. y2=xcoor(irc+2)-xcoor(ira+2)
  138. pmix=x1*y2-y1*x2
  139. elseif (idim.eq.3) then
  140. *dbg write(ioimp,*) 'lnomas1=',lnomas(iadr+1)
  141. ia=ipt1.num(lnomas(iadr+1),il)
  142. ib=ipt1.num(lnomas(iadr+2),il)
  143. ic=ipt1.num(lnomas(iadr+3),il)
  144. id=ipt1.num(lnomas(iadr+4),il)
  145. *dbg write(ioimp,*) 'ia=',ia,' ib=',ib,' ic=',ic,' id=',id
  146. ira=(idim+1)*(ia-1)
  147. irb=(idim+1)*(ib-1)
  148. irc=(idim+1)*(ic-1)
  149. ird=(idim+1)*(id-1)
  150. x1=xcoor(irb+1)-xcoor(ira+1)
  151. y1=xcoor(irb+2)-xcoor(ira+2)
  152. z1=xcoor(irb+3)-xcoor(ira+3)
  153. x2=xcoor(irc+1)-xcoor(ira+1)
  154. y2=xcoor(irc+2)-xcoor(ira+2)
  155. z2=xcoor(irc+3)-xcoor(ira+3)
  156. x3=xcoor(ird+1)-xcoor(ira+1)
  157. y3=xcoor(ird+2)-xcoor(ira+2)
  158. z3=xcoor(ird+3)-xcoor(ira+3)
  159. pmix=x1*(y2*z3-y3*z2)+y1*(z2*x3-z3*x2)+z1*(x2*y3-x3*y2)
  160. endif
  161. *dbg write(ioimp,*) 'pmix=',pmix
  162. if (pmix.ge.0.d0) then
  163. lect(il)=+1
  164. else
  165. ninv=ninv+1
  166. lect(il)=-1
  167. endif
  168. enddo
  169. if (ninv.gt.0) then
  170. call inver4(ipt1,ipt2,3,mlenti)
  171. if (ierr.ne.0) return
  172. else
  173. ipt2=ipt1
  174. endif
  175. segsup mlenti
  176. * fin du cas massif
  177. return
  178. *
  179. * Cas de l'orientation des éléments non massifs (sauf cub8 shb8)
  180. *
  181. *
  182. 1000 continue
  183. IF (IDIM.EQ.3) THEN
  184. IF (ICLE.NE.1.AND.ICLE.NE.2) CALL ERREUR(1074)
  185. ELSE
  186. * SG Ceci est un residu de l'ancienne programmation en 2D,
  187. * ou on oriente par rapport a un vecteur (0. 0. 1.). En fait, comme
  188. * on n'orientait que les TRI et les QUA, ceux-ci ont deja ete
  189. * traites dans le cas massif. Du coup, on ne passera normalement ici
  190. * que pour les elements que l'on ne sait pas orienter (SEG en 2D par ex)
  191. * Voir IPT2=IPT1 plus loin
  192. ICLE=1
  193. ENDIF
  194.  
  195. IF (ICLE.EQ.1) THEN
  196. * NORMALISATION DU VECTEUR
  197. DP=SQRT(XP(1)**2+XP(2)**2+XP(3)**2)
  198. IF (DP.LE.XPETIT) THEN
  199. CALL ERREUR(277)
  200. DP=1.D0
  201. ENDIF
  202. XQ(1)=XP(1)/DP
  203. XQ(2)=XP(2)/DP
  204. XQ(3)=XP(3)/DP
  205. ENDIF
  206. NBREF=IPT1.LISREF(/1)
  207. NBSOUS=0
  208. NBELEM=IPT1.NUM(/2)
  209. NBNN=IPT1.NUM(/1)
  210. ITYP=KSURF(IPT1.ITYPEL)
  211. IF (ITYP.NE.0.AND.ITYP.EQ.IPT1.ITYPEL) GOTO 1
  212. *
  213. * ce n'est pas des éléments de surface
  214. * est-on en présence de cub8 (shb8)
  215. *
  216. if(ipt1.itypel.eq.14) then
  217. nbref=2
  218. segini ipt2
  219. nbnn=4
  220. nbref=0
  221. segini ipt3,ipt4
  222. ipt2.lisref(1)=ipt3
  223. ipt2.lisref(2)=ipt4
  224. ipt3.itypel=8
  225. ipt4.itypel=8
  226. ipt2.itypel=14
  227. idim1=idim+1
  228. do i=1,ipt1.num(/2)
  229. ia=ipt1.num(1,i) - 1
  230. ib=ipt1.num(5,i) - 1
  231. ic=ipt1.num(2,i) - 1
  232. id=ipt1.num(3,i) - 1
  233. xab=xcoor(ib*idim1+1)-xcoor(ia*idim1+1)
  234. yab=xcoor(ib*idim1+2)-xcoor(ia*idim1+2)
  235. zab=xcoor(ib*idim1+3)-xcoor(ia*idim1+3)
  236. xac=xcoor(ic*idim1+1)-xcoor(ia*idim1+1)
  237. yac=xcoor(ic*idim1+2)-xcoor(ia*idim1+2)
  238. zac=xcoor(ic*idim1+3)-xcoor(ia*idim1+3)
  239. xad=xcoor(id*idim1+1)-xcoor(ia*idim1+1)
  240. yad=xcoor(id*idim1+2)-xcoor(ia*idim1+2)
  241. zad=xcoor(id*idim1+3)-xcoor(ia*idim1+3)
  242. xpvec= yac*zad - zac*yad
  243. ypvec= zac*xad - xac*zad
  244. zpvec= xac*yad - yac*xad
  245. pmix=xpvec*xab+ypvec*yab+zpvec*zab
  246. isens=0
  247. if(pmix.ge.0.d0) isens=1
  248. if(icle.eq.1) then
  249. xsc=xab*xq(1)+yab*xq(2)+zab*xq(3)
  250. else
  251. xsc= xab*(xq(1)-xcoor(ib*idim1+1))+
  252. $ yab*(xq(2)-xcoor(ib*idim1+2))+
  253. $ zab*(xq(3)-xcoor(ib*idim1+3))
  254. endif
  255. if(xsc.ge.0.d0) then
  256. do j=1,8
  257. ipt2.num(j,i)=ipt1.num(j,i)
  258. enddo
  259. if(isens.eq.2) then
  260. iu=ipt2.num(2,i)
  261. ipt2.num(2,i)=ipt2.num(4,i)
  262. ipt2.num(4,i)=iu
  263. iu=ipt2.num(6,i)
  264. ipt2.num(6,i)=ipt2.num(8,i)
  265. ipt2.num(8,i)=iu
  266. endif
  267. else
  268. ipt2.num(1,i)=ipt1.num(5,i)
  269. ipt2.num(2,i)=ipt1.num(8,i)
  270. ipt2.num(3,i)=ipt1.num(7,i)
  271. ipt2.num(4,i)=ipt1.num(6,i)
  272. ipt2.num(5,i)=ipt1.num(1,i)
  273. ipt2.num(6,i)=ipt1.num(4,i)
  274. ipt2.num(7,i)=ipt1.num(3,i)
  275. ipt2.num(8,i)=ipt1.num(2,i)
  276. if(isens.eq.2) then
  277. iu=ipt2.num(2,i)
  278. ipt2.num(2,i)=ipt2.num(4,i)
  279. ipt2.num(4,i)=iu
  280. iu=ipt2.num(6,i)
  281. ipt2.num(6,i)=ipt2.num(8,i)
  282. ipt2.num(8,i)=iu
  283. endif
  284. endif
  285. ipt2.icolor(i)=ipt1.icolor(i)
  286. ipt3.num(1,i)=ipt2.num(1,i)
  287. ipt3.num(2,i)=ipt2.num(2,i)
  288. ipt3.num(3,i)=ipt2.num(3,i)
  289. ipt3.num(4,i)=ipt2.num(4,i)
  290. ipt3.icolor(i)=ipt2.icolor(i)
  291. ipt4.num(1,i)=ipt2.num(5,i)
  292. ipt4.num(2,i)=ipt2.num(6,i)
  293. ipt4.num(3,i)=ipt2.num(7,i)
  294. ipt4.num(4,i)=ipt2.num(8,i)
  295. ipt4.icolor(i)=ipt2.icolor(i)
  296. enddo
  297. segdes ipt3,ipt4
  298. return
  299. endif
  300. *
  301. * Sinon, on ne savait pas orienter les éléments de ipt1
  302. *
  303. IPT2=IPT1
  304. RETURN
  305. *
  306. * on a des vrais éléments de surfaces
  307. *
  308. 1 CONTINUE
  309. * Cas QUAF TRI7 ou QUA9
  310. LQUAF=(ITYP.EQ.11.OR.ITYP.EQ.7)
  311. SEGINI IPT2
  312. IPT2.ITYPEL=ITYP
  313. IF (NBREF.EQ.0) GOTO 3
  314. DO 2 I=1,NBREF
  315. IPT2.LISREF(I)=IPT1.LISREF(I)
  316. 2 CONTINUE
  317. 3 CONTINUE
  318. SEGACT MCOORD
  319. IB=NSPOS(ITYP)-1
  320. IS1=IBSOM(IB+1)
  321. IS2=IBSOM(IB+2)
  322. IS3=IBSOM(IB+3)
  323.  
  324.  
  325. DO 4 J=1,NBELEM
  326. IP1=IPT1.NUM(IS1,J)
  327. IP2=IPT1.NUM(IS2,J)
  328. IP3=IPT1.NUM(IS3,J)
  329. IREF1=(IDIM+1)*(IP1-1)
  330. IREF2=(IDIM+1)*(IP2-1)
  331. IREF3=(IDIM+1)*(IP3-1)
  332. XV1=XCOOR(IREF2+1)-XCOOR(IREF1+1)
  333. YV1=XCOOR(IREF2+2)-XCOOR(IREF1+2)
  334. ZV1=XCOOR(IREF2+3)-XCOOR(IREF1+3)
  335. IF (IDIM.NE.3) ZV1=0.D0
  336. XV2=XCOOR(IREF2+1)-XCOOR(IREF3+1)
  337. YV2=XCOOR(IREF2+2)-XCOOR(IREF3+2)
  338. ZV2=XCOOR(IREF2+3)-XCOOR(IREF3+3)
  339. IF (IDIM.NE.3) ZV2=0.D0
  340. XV3=ZV1*YV2-ZV2*YV1
  341. YV3=XV1*ZV2-XV2*ZV1
  342. ZV3=YV1*XV2-YV2*XV1
  343. XVN=SQRT(XV3**2+YV3**2+ZV3**2)
  344. IF (ICLE.EQ.2) THEN
  345. XQ(1)=0.D0
  346. XQ(2)=0.D0
  347. XQ(3)=0.D0
  348. DO 10 L=1,NBNN
  349. IREF=(IPT1.NUM(L,J)-1)*(IDIM+1)
  350. DO 11 K=1,3
  351. XQ(K)=XQ(K)+XCOOR(IREF+K)
  352. 11 CONTINUE
  353. 10 CONTINUE
  354. XQ(1)=XQ(1)/NBNN
  355. XQ(2)=XQ(2)/NBNN
  356. XQ(3)=XQ(3)/NBNN
  357. XQ(1)=XP(1)-XQ(1)
  358. XQ(2)=XP(2)-XQ(2)
  359. XQ(3)=XP(3)-XQ(3)
  360. DP=SQRT(XQ(1)**2+XQ(2)**2+XQ(3)**2)
  361. IF (DP.LE.XPETIT) THEN
  362. CALL ERREUR(277)
  363. DP=1.D0
  364. ENDIF
  365. XQ(1)=XQ(1)/DP
  366. XQ(2)=XQ(2)/DP
  367. XQ(3)=XQ(3)/DP
  368. ENDIF
  369. TEST=XV3*XQ(1)+YV3*XQ(2)+ZV3*XQ(3)
  370. IF (ABS(TEST).LE.0.0175D0*XVN) CALL ERREUR(278)
  371. IF (TEST.GE.0.D0) THEN
  372. DO 6 I=1,NBNN
  373. IPT2.NUM(I,J)=IPT1.NUM(I,J)
  374. 6 CONTINUE
  375. ELSE
  376. IF (.NOT.LQUAF) THEN
  377. DO 7 I=1,NBNN
  378. IPT2.NUM(MOD(NBNN+1-I,NBNN)+1,J)=IPT1.NUM(I,J)
  379. 7 CONTINUE
  380. ELSE
  381. NBNN1=NBNN-1
  382. DO 8 I=1,NBNN1
  383. IPT2.NUM(MOD(NBNN1+1-I,NBNN1)+1,J)=IPT1.NUM(I,J)
  384. 8 CONTINUE
  385. IPT2.NUM(NBNN,J)=IPT1.NUM(NBNN,J)
  386. ENDIF
  387. ENDIF
  388. IPT2.ICOLOR(J)=IPT1.ICOLOR(J)
  389. 4 CONTINUE
  390. RETURN
  391. END
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  

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