Télécharger orient.eso

Retour à la liste

Numérotation des lignes :

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

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