Télécharger dedu.eso

Retour à la liste

Numérotation des lignes :

  1. C DEDU SOURCE BP208322 16/11/18 21:16:11 9177
  2.  
  3. C-----------------------------------------------------------------------
  4. C 1ere possibilite
  5. C OBJ1 = DEDU OBJ2 MAIT_ANC MAIT_NOU (REGU)
  6. C Deduction d'un maillage a partir d'un autre et du deplacement de
  7. C noeuds maitres
  8. C AUTEUR : PV
  9. C-----------------------------------------------------------------------
  10. C 2eme possibilite
  11. C OBJ2 = OBJ1 DEDU 'TRAN' GEO1 GEO2
  12. C 3eme possibilite
  13. C OBJ2 = OBJ1 DEDU FLOT1 POIN1 (POIN2) 'ROTA' GEO1 GEO2
  14. C OBJ1 : type POINT, MAILLAGE, CHPOINT, MCHAML, MMODEL
  15. C Passage du support de OBJ1 a celui de OBJ2 comme celui de GEO1 a GEO2
  16. C Les points du support de OBJ1 doivent appartenir a GEO1
  17. C On realise eventuellement la rotation des composantes (envoi a PROPER)
  18. C AUTEUR : KICH
  19. C-----------------------------------------------------------------------
  20. C 10/2003 : Prise en compte du cas IDIM=1 et des modes associes.
  21. C-----------------------------------------------------------------------
  22. C Remarques :
  23. C -----------
  24. C DEDU est appele par l'operateur DEPL (deplac.eso) avec ICAS=1.
  25. C Ce cas correspond a une transformation affine de l'espace (noeuds
  26. C maitres = un maillage a (IDIM+1) points).
  27. C-----------------------------------------------------------------------
  28.  
  29. SUBROUTINE DEDU(ICAS)
  30.  
  31. IMPLICIT INTEGER(I-N)
  32. IMPLICIT REAL*8 (A-H,O-Z)
  33.  
  34. -INC CCOPTIO
  35. -INC SMCOORD
  36. -INC SMELEME
  37. -INC CCGEOME
  38. -INC SMRIGID
  39. -INC SMCHPOI
  40. -INC SMMODEL
  41. -INC SMCHAML
  42. -INC SMTEXTE
  43. DIMENSION iuniq(10)
  44. CHARACTER*4 mcle(5),strcomp
  45.  
  46. SEGMENT kon(nkon,numnp,2)
  47. SEGMENT xkon(nkon+1,numnp)
  48. SEGMENT ipxkon(nkon+1,numnp)
  49. SEGMENT ikon(numnp)
  50. SEGMENT icpr(xcoor(/1)/(IDIM+1))
  51. SEGMENT icp2(xcoor(/1)/(IDIM+1))
  52. SEGMENT idcp(numnp)
  53. SEGMENT icorr(numnp)
  54. SEGMENT xxtra
  55. REAL*8 xdep(3,4),xarr(3,4),ved1(3),ved2(3),ved3(3),vea1(3),
  56. . vea2(3),vea3(3),pa(3,3),pa2(2,2),pa3(3,3),vb(3)
  57. END SEGMENT
  58.  
  59. DATA mcle / 'REGU','ELIM','TRAN','ROTA','ADAP' /
  60.  
  61. idimp1=IDIM+1
  62.  
  63. C Mot cle optionnel REGU : coordonnees barycentriques identiques
  64. CALL LIRMOT(mcle,5,iregu,0)
  65.  
  66. C 2eme possibilite : option TRAN
  67. C --------------------
  68. IF (iregu.EQ.3) THEN
  69. CALL PROPER(4)
  70. RETURN
  71. ENDIF
  72. C 3eme possibilite : option ROTA
  73. C --------------------
  74. C Option non disponible en dimension 1 (pas de signification)
  75. IF (iregu.EQ.4) THEN
  76. IF (IDIM.EQ.1) THEN
  77. MOTERR(1:4)=mcle(iregu)
  78. INTERR(1)=IDIM
  79. CALL ERREUR(971)
  80. ELSE
  81. CALL PROPER(5)
  82. ENDIF
  83. RETURN
  84. ENDIF
  85. C 5eme possibilite : option ADAP
  86. C On passe la main à la procédure DEDUADAP
  87. C --------------------
  88. IF (iregu.EQ.5) THEN
  89. SEGINI MTEXTE
  90. LTT=8
  91. MTEXT(1:LTT) ='DEDUADAP'
  92. NCART=8
  93. SEGDES MTEXTE
  94. CALL ECROBJ('TEXTE',MTEXTE)
  95. RETURN
  96. ENDIF
  97.  
  98. C Par defaut : option ELIM (ignore les parties concaves)
  99. C Correspond au cas ICAS=1 (operateur DEPL 'DEDU')
  100. IF (iregu.EQ.0) iregu=2
  101.  
  102. C 1ere possibilite : options REGU et ELIM / appel par DEPLAC
  103. C --------------------
  104. C Lecture du maillage initial
  105. CALL LIROBJ('MAILLAGE',ipt1,1,iretou)
  106. IF (ierr.NE.0) RETURN
  107. C Lecture des points maitres initiaux
  108. CALL LIROBJ('MAILLAGE',ipt2,0,iretou)
  109. IF (iretou.NE.0) THEN
  110. ichpo=0
  111. C lecture ddes points maitres finaux
  112. CALL LIROBJ('MAILLAGE',ipt3,1,iretou)
  113. IF (ierr.NE.0) RETURN
  114. ELSE
  115. C Lecture d'un chpoint de deplacement
  116. ichpo=1
  117. CALL LIROBJ('CHPOINT ',mchpoi,1,iretou)
  118. IF (ierr.NE.0) RETURN
  119. C Création de ipt2 et ipt3 pour se ramener au cas précédent
  120. nbpts=xcoor(/1)/idimp1
  121. nbpts0=nbpts
  122. C Comptage des points du chpoint concernes par UX UY UZ UR
  123. SEGACT mchpoi
  124. jatan=jattri(1)
  125. nbpcon=0
  126. DO nso=1,ipchp(/1)
  127. msoupo=ipchp(nso)
  128. SEGACT msoupo
  129. ipt6=igeoc
  130. SEGACT ipt6
  131. ipc=0
  132. DO i=1,nocomp(/2)
  133. strcomp=nocomp(i)
  134. IF ((strcomp.EQ.'UX ').OR.(strcomp.EQ.'UR ').OR.
  135. . ((strcomp.EQ.'UY ').AND.(IDIM.NE.1)).OR.
  136. . ((strcomp.EQ.'UZ ').AND.(IDIM.NE.1))) ipc=ipc+i
  137. ENDDO
  138. IF (ipc.NE.0) nbpcon=nbpcon+ipt6.num(/2)
  139. ENDDO
  140. C Test sur nbpcon=0 a ecrire
  141. nbref=0
  142. nbsous=0
  143. nbelem=nbpcon
  144. nbnn=1
  145. SEGINI ipt2
  146. SEGINI ipt3
  147. ipt2.itypel=1
  148. ipt3.itypel=1
  149. C Remplissage de ipt2 et ipt3 - Creation des points de ipt3
  150. idej=0
  151. DO 9 nso=1,ipchp(/1)
  152. msoupo=ipchp(nso)
  153. ipt6=igeoc
  154. ipox=0
  155. ipoy=0
  156. ipoz=0
  157. DO i=1,nocomp(/2)
  158. strcomp=nocomp(i)
  159. IF ((strcomp.EQ.'UX ').OR.(strcomp.EQ.'UR ')) ipox=i
  160. IF ((strcomp.EQ.'UY ').AND.(IDIM.NE.1)) ipoy=i
  161. IF (strcomp.EQ.'UZ ') THEN
  162. IF (IDIM.EQ.2) THEN
  163. ipoy=i
  164. ELSE IF (IDIM.GE.3) THEN
  165. ipoz=i
  166. ENDIF
  167. ENDIF
  168. ENDDO
  169. IF ((ipox+ipoy+ipoz).EQ.0) GOTO 9
  170. mpoval=ipoval
  171. SEGACT mpoval
  172. nbpan=xcoor(/1)/idimp1
  173. nbpts=nbpts+ipt6.num(/2)
  174. SEGADJ mcoord
  175. DO nel=1,ipt6.num(/2)
  176. iep=ipt6.num(1,nel)
  177. ipt2.num(1,idej+nel)=iep
  178. ipt3.num(1,idej+nel)=nbpan+nel
  179. irefep=(iep-1)*idimp1
  180. iref=(nbpan+nel-1)*idimp1
  181. IF (ipox.NE.0) THEN
  182. xcoor(iref+1)=xcoor(irefep+1)+vpocha(nel,ipox)
  183. ELSE
  184. xcoor(iref+1)=xcoor(irefep+1)
  185. ENDIF
  186. IF (ipoy.NE.0) THEN
  187. xcoor(iref+2)=xcoor(irefep+2)+vpocha(nel,ipoy)
  188. ELSE
  189. IF (IDIM.GE.2) xcoor(iref+2)=xcoor(irefep+2)
  190. ENDIF
  191. IF (ipoz.NE.0) THEN
  192. xcoor(iref+3)=xcoor(irefep+3)+vpocha(nel,ipoz)
  193. ELSE
  194. IF (IDIM.GE.3) xcoor(iref+3)=xcoor(irefep+3)
  195. ENDIF
  196. xcoor(iref+idimp1)=xcoor(irefep+idimp1)
  197. ENDDO
  198. idej=idej+ipt6.num(/2)
  199. SEGDES ipt6,mpoval,msoupo
  200. 9 CONTINUE
  201. SEGDES mchpoi
  202. ENDIF
  203.  
  204. C Verification compatibilite points maitres initiaux et finaux
  205. C Cas ichpo = 0 : cette verification est a faire lorsque la
  206. C transformation est definie par deux maillages.
  207. C Cas ichpo = 1 : les points maitres initiaux et finaux sont crees
  208. C en parallele de la meme maniere a partir du CHPO donne. La veri-
  209. C fication est donc inutile.
  210. IF (ichpo.EQ.1) GOTO 100
  211. SEGACT ipt2,ipt3
  212. IF (ipt2.lisous(/1).NE.ipt3.lisous(/1)) THEN
  213. CALL ERREUR(16)
  214. SEGDES ipt2,ipt3
  215. RETURN
  216. ENDIF
  217. ipt7=ipt2
  218. ipt8=ipt3
  219. DO i=1,MAX(1,ipt2.lisous(/1))
  220. IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(i)
  221. IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(i)
  222. SEGACT ipt7,ipt8
  223. moterr=' '
  224. IF (ipt7.num(/1).NE.ipt8.num(/1)) moterr(1:8)='noeuds'
  225. IF (ipt7.num(/2).NE.ipt8.num(/2)) moterr(1:8)='elements'
  226. SEGDES ipt7,ipt8
  227. IF (moterr.NE.' ') THEN
  228. CALL ERREUR(550)
  229. RETURN
  230. ENDIF
  231. ENDDO
  232.  
  233. 100 SEGACT MCOORD*MOD
  234. C Creation du tableau de correspondance icpr donnant pour chaque
  235. C noeud de ipt2 son noeud image de ipt3
  236. SEGACT ipt2,ipt3
  237. SEGINI icpr
  238. numnp=0
  239. ipt7=ipt2
  240. ipt8=ipt3
  241. DO i=1,MAX(1,ipt2.lisous(/1))
  242. IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(i)
  243. IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(i)
  244. SEGACT ipt7,ipt8
  245. DO k=1,ipt7.num(/2)
  246. DO l=1,ipt7.num(/1)
  247. ip=ipt7.num(l,k)
  248. IF (icpr(ip).EQ.0) THEN
  249. numnp=numnp+1
  250. icpr(ip)=ipt8.num(l,k)
  251. ENDIF
  252. ENDDO
  253. ENDDO
  254. SEGDES ipt7,ipt8
  255. ENDDO
  256. SEGDES ipt2,ipt3
  257.  
  258. C Est-ce une transformation affine ?
  259. IF (numnp.LE.idimp1) THEN
  260. SEGINI xxtra
  261. numnp=0
  262. DO i=1,icpr(/1)
  263. IF (icpr(i).NE.0) THEN
  264. numnp=numnp+1
  265. iref=(i-1)*idimp1
  266. irefep=(icpr(i)-1)*idimp1
  267. DO j=1,idim
  268. xdep(j,numnp)=xcoor(iref+j)
  269. xarr(j,numnp)=xcoor(irefep+j)
  270. ENDDO
  271. ENDIF
  272. ENDDO
  273. SEGSUP icpr
  274. C Faut-il creer de nouveaux points ?
  275. IF (numnp.EQ.1) THEN
  276. C On en cree un par translation en x
  277. numnp=numnp+1
  278. DO j=1,IDIM
  279. xdep(j,numnp)=xdep(j,numnp-1)
  280. xarr(j,numnp)=xarr(j,numnp-1)
  281. ENDDO
  282. xdep(1,numnp)=xdep(1,numnp)+1.
  283. xarr(1,numnp)=xarr(1,numnp)+1.
  284. ENDIF
  285. DO j=1,IDIM
  286. ved1(j)=xdep(j,2)-xdep(j,1)
  287. vea1(j)=xarr(j,2)-xarr(j,1)
  288. ENDDO
  289. IF (numnp.EQ.2) THEN
  290. C En 1D, on ne fait rien.
  291. C En 2D, on cherche un point sur la perpendiculaire a la droite
  292. C definie par les deux premiers points, passant par le premier
  293. C point et a une distance unite (formant un repere orthonorme direct).
  294. C En 3D, on prend un point dans la direction y.
  295. numnp=numnp+1
  296. IF (IDIM.EQ.2) THEN
  297. xdep(1,3)=ved1(2)
  298. xdep(2,3)=-ved1(1)
  299. aa=1./SQRT(xdep(1,3)*xdep(1,3)+xdep(2,3)*xdep(2,3))
  300. xdep(1,3)=xdep(1,1)+xdep(1,3)*aa
  301. xdep(2,3)=xdep(2,1)+xdep(2,3)*aa
  302. xarr(1,3)=vea1(2)
  303. xarr(2,3)=-vea1(2)
  304. aa=1./SQRT(xarr(1,3)*xarr(1,3)+xarr(2,3)*xarr(2,3))
  305. xarr(1,3)=xarr(1,1)+xarr(1,3)*aa
  306. xarr(2,3)=xarr(2,1)+xarr(2,3)*aa
  307. ELSE IF (IDIM.EQ.3) THEN
  308. xdep(1,3)=xdep(1,1)
  309. xdep(2,3)=xdep(2,1)+1.
  310. xdep(3,3)=xdep(3,1)
  311. xarr(1,3)=xarr(1,1)
  312. xarr(2,3)=xarr(2,1)+1.
  313. xarr(3,3)=xarr(3,1)
  314. ENDIF
  315. ENDIF
  316. IF (IDIM.NE.1) THEN
  317. DO j=1,IDIM
  318. ved2(j)=xdep(j,3)-xdep(j,1)
  319. vea2(j)=xarr(j,3)-xarr(j,1)
  320. ENDDO
  321. ENDIF
  322. C En 3D, il faut encore un point hors du plan forme par les 3 points
  323. IF (numnp.LT.idimp1) THEN
  324. ved3(1)=ved1(2)*ved2(3)-ved1(3)*ved2(2)
  325. ved3(2)=ved1(3)*ved2(1)-ved1(1)*ved2(3)
  326. ved3(3)=ved1(1)*ved2(2)-ved1(2)*ved2(1)
  327. aa=1./SQRT(ved3(1)*ved3(1)+ved3(2)*ved3(2)+ved3(3)*ved3(3))
  328. vea3(1)=vea1(2)*vea2(3)-vea1(3)*vea2(2)
  329. vea3(2)=vea1(3)*vea2(1)-vea1(1)*vea2(3)
  330. vea3(3)=vea1(1)*vea2(2)-vea1(2)*vea2(1)
  331. bb=1./SQRT(vea3(1)*vea3(1)+vea3(2)*vea3(2)+vea3(3)*vea3(3))
  332. DO j=1,3
  333. ved3(j)=ved3(j)*aa
  334. vea3(j)=vea3(j)*bb
  335. xdep(j,4)=xdep(j,1)+ved3(j)
  336. xarr(j,4)=xarr(j,1)+vea3(j)
  337. ENDDO
  338. ENDIF
  339. IF (IDIM.EQ.3) THEN
  340. DO j=1,3
  341. ved3(j)=xdep(j,4)-xdep(j,1)
  342. vea3(j)=xarr(j,4)-xarr(j,1)
  343. ENDDO
  344. ENDIF
  345. C Calcul de la matrice de passage pa du repere ved au repere vea
  346. C On traite le cas IDIM=1 a part (sinon reso33 serait a modifier).
  347. IF (IDIM.EQ.1) THEN
  348. IF ((ved1(1).EQ.0.).OR.(vea1(1).EQ.0.)) THEN
  349. CALL ERREUR(21)
  350. SEGSUP xxtra
  351. IF (ichpo.EQ.1) THEN
  352. nbpts=nbpts0
  353. SEGADJ mcoord
  354. SEGSUP ipt2,ipt3
  355. ENDIF
  356. RETURN
  357. ENDIF
  358. pa(1,1)=vea1(1)/ved1(1)
  359. ELSE
  360. DO j=1,IDIM
  361. vb(1)=vea1(j)
  362. vb(2)=vea2(j)
  363. IF (IDIM.EQ.3) vb(3)=vea3(j)
  364. pa3(1,1)=ved1(1)
  365. pa3(1,2)=ved1(2)
  366. pa3(2,1)=ved2(1)
  367. pa3(2,2)=ved2(2)
  368. IF (IDIM.EQ.3) THEN
  369. pa3(1,3)=ved1(3)
  370. pa3(2,3)=ved2(3)
  371. pa3(3,1)=ved3(1)
  372. pa3(3,2)=ved3(2)
  373. pa3(3,3)=ved3(3)
  374. ENDIF
  375. CALL reso33(pa3,vb,3,IDIM,kerr)
  376. C Erreur : en 2D -> 3 points alignes, en 3D -> 4 points coplanaires
  377. IF (kerr.NE.0) THEN
  378. CALL ERREUR(21)
  379. SEGSUP xxtra
  380. IF (ichpo.EQ.1) THEN
  381. nbpts=nbpts0
  382. SEGADJ mcoord
  383. SEGSUP ipt2,ipt3
  384. ENDIF
  385. RETURN
  386. ENDIF
  387. pa(j,1)=vb(1)
  388. pa(j,2)=vb(2)
  389. IF (IDIM.EQ.3) pa(j,3)=vb(3)
  390. ENDDO
  391. ENDIF
  392. C Fabrication de la liste des points a transformer et du maillage
  393. C resultat (segment MELEME)
  394. SEGINI icpr
  395. nbpan=xcoor(/1)/idimp1
  396. SEGACT ipt1
  397. SEGINI,ipt4=ipt1
  398. nbnn=ipt4.num(/1)
  399. nbelem=ipt4.num(/2)
  400. nbsous=ipt4.lisous(/1)
  401. nbref=ipt4.lisref(/1)
  402. IF (ICAS.NE.1) THEN
  403. nbref=0
  404. SEGADJ ipt4
  405. ENDIF
  406. ipt8=ipt1
  407. ipt5=ipt4
  408. numnp=0
  409. DO i=1,MAX(1,ipt1.lisous(/1))
  410. IF (ipt1.lisous(/1).NE.0) THEN
  411. ipt8=ipt1.lisous(i)
  412. SEGACT ipt8
  413. SEGINI,ipt5=ipt8
  414. nbnn=ipt5.num(/1)
  415. nbelem=ipt5.num(/2)
  416. nbsous=ipt5.lisous(/1)
  417. nbref=ipt5.lisref(/1)
  418. IF (ICAS.NE.1) THEN
  419. nbref=0
  420. SEGADJ ipt5
  421. ENDIF
  422. ipt4.lisous(i)=ipt5
  423. ENDIF
  424. DO iel=1,ipt8.num(/2)
  425. DO inn=1,ipt8.num(/1)
  426. ipt=ipt8.num(inn,iel)
  427. IF (icpr(ipt).EQ.0) THEN
  428. IF (ICAS.EQ.1) THEN
  429. icpr(ipt)=ipt
  430. ELSE
  431. numnp=numnp+1
  432. icpr(ipt)=numnp+nbpan
  433. ENDIF
  434. ENDIF
  435. ipt5.num(inn,iel)=icpr(ipt)
  436. ENDDO
  437. ENDDO
  438. SEGDES ipt8,ipt5
  439. ENDDO
  440. SEGDES ipt4,ipt1
  441. IF (ICAS.EQ.0) THEN
  442. nbpts=numnp+nbpan
  443. SEGADJ mcoord
  444. ENDIF
  445. DO ipt=1,nbpan
  446. IF (icpr(ipt).NE.0) THEN
  447. jno=(ipt-1)*idimp1
  448. ino=(icpr(ipt)-1)*idimp1
  449. DO i=1,IDIM
  450. vb(i)=xcoor(jno+i)-xdep(i,1)
  451. ENDDO
  452. DO i=1,IDIM
  453. aa=xarr(i,1)
  454. DO j=1,IDIM
  455. aa=aa+pa(i,j)*vb(j)
  456. ENDDO
  457. xcoor(ino+i)=aa
  458. ENDDO
  459. xcoor(ino+idimp1)=xcoor(jno+idimp1)
  460. ENDIF
  461. ENDDO
  462. SEGSUP icpr,xxtra
  463. CALL ECROBJ('MAILLAGE',ipt4)
  464. IF (ichpo.EQ.1) THEN
  465. SEGSUP ipt2,ipt3
  466. ENDIF
  467. RETURN
  468. ENDIF
  469.  
  470. C Cas d'une transformation non affine
  471. C Interdit dans le cas ICAS=1 (appel par deplac, operateur DEPL 'DEDU')
  472. C On recupere le maillage initial si ichpo=1.
  473. IF (ICAS.EQ.1) THEN
  474. CALL ERREUR(21)
  475. IF (ichpo.EQ.1) THEN
  476. nbpts=nbpts0
  477. SEGADJ mcoord
  478. SEGSUP ipt2,ipt3
  479. ENDIF
  480. RETURN
  481. ENDIF
  482. C CHOIX : Pour l'instant, cas non implemente en DIMEnsion 1.
  483. IF (IDIM.EQ.1) THEN
  484. MOTERR(1:4)=mcle(iregu)
  485. INTERR(1)=IDIM
  486. CALL ERREUR(971)
  487. IF (ichpo.EQ.1) THEN
  488. nbpts=nbpts0
  489. SEGADJ mcoord
  490. SEGSUP ipt2,ipt3
  491. ENDIF
  492. RETURN
  493. ENDIF
  494. C Creation d'un tableau de stockage contenant pour chaque point les
  495. C aretes des elements le touchant sauf celles qui le contiennent
  496. SEGACT ipt1
  497. SEGINI icpr
  498. numnp=0
  499. ipt8=ipt1
  500. DO i=1,MAX(1,ipt1.lisous(/1))
  501. IF (ipt1.lisous(/1).NE.0) ipt8=ipt1.lisous(i)
  502. SEGACT ipt8
  503. DO iel=1,ipt8.num(/2)
  504. DO inn=1,ipt8.num(/1)
  505. ipt=ipt8.num(inn,iel)
  506. IF (icpr(ipt).EQ.0) THEN
  507. numnp=numnp+1
  508. icpr(ipt)=numnp
  509. ENDIF
  510. ENDDO
  511. ENDDO
  512. ENDDO
  513. SEGINI idcp
  514. DO i=1,icpr(/1)
  515. IF (icpr(i).NE.0) idcp(icpr(i))=i
  516. ENDDO
  517. nkon=5
  518. SEGINI kon,ikon
  519. DO i=1,MAX(1,ipt1.lisous(/1))
  520. IF (ipt1.lisous(/1).NE.0) ipt8=ipt1.lisous(i)
  521. DO iel=1,ipt8.num(/2)
  522. DO inn=1,ipt8.num(/1)
  523. ip=ipt8.num(inn,iel)
  524. ipc=icpr(ip)
  525. DO 41 inp=1,ipt8.num(/1)
  526. kpt1=ipt8.num(inp,iel)
  527. IF (kpt1.EQ.ip) GOTO 41
  528. inps=inp+1
  529. IF (inps.GT.ipt8.num(/1)) inps=1
  530. kpt2=ipt8.num(inps,iel)
  531. IF (kpt2.EQ.ip) GOTO 41
  532. DO idj=1,ikon(ipc)
  533. IF ((kon(idj,ipc,1).EQ.kpt1).AND.
  534. . (kon(idj,ipc,2).EQ.kpt2)) GOTO 41
  535. IF ((kon(idj,ipc,1).EQ.kpt2).AND.
  536. . (kon(idj,ipc,2).EQ.kpt1)) GOTO 41
  537. ENDDO
  538. ikon(ipc)=ikon(ipc)+1
  539. IF (ikon(ipc).GT.nkon) THEN
  540. nkon=ikon(ipc)
  541. SEGADJ kon
  542. ENDIF
  543. kon(ikon(ipc),ipc,1)=kpt1
  544. kon(ikon(ipc),ipc,2)=kpt2
  545. 41 CONTINUE
  546. ENDDO
  547. ENDDO
  548. ENDDO
  549. C Fermeture du contour si necessaire (avec 1 segment supplementaire)
  550. DO 80 ipc=1,kon(/2)
  551. iuni=0
  552. DO iko=1,ikon(ipc)
  553. kpt1=kon(iko,ipc,1)
  554. kpt2=kon(iko,ipc,2)
  555. DO ikq=1,ikon(ipc)
  556. IF (ikq.NE.iko) THEN
  557. IF (kpt1.EQ.kon(ikq,ipc,1)) GOTO 86
  558. IF (kpt1.EQ.kon(ikq,ipc,2)) GOTO 86
  559. ENDIF
  560. ENDDO
  561. C kpt1 unique
  562. iuni=iuni+1
  563. iuniq(iuni)=kpt1
  564. 86 DO 88 ikq=1,ikon(ipc)
  565. IF (ikq.EQ.iko) GOTO 88
  566. IF (kpt2.EQ.kon(ikq,ipc,1)) GOTO 90
  567. IF (kpt2.EQ.kon(ikq,ipc,2)) GOTO 90
  568. 88 CONTINUE
  569. C kpt1 unique
  570. iuni=iuni+1
  571. iuniq(iuni)=kpt2
  572. 90 CONTINUE
  573. ENDDO
  574. IF (iuni.EQ.0) GOTO 80
  575. ikon(ipc)=ikon(ipc)+1
  576. IF (ikon(ipc).GT.nkon) THEN
  577. nkon=ikon(ipc)
  578. SEGADJ kon
  579. ENDIF
  580. kon(ikon(ipc),ipc,1)=iuniq(1)
  581. kon(ikon(ipc),ipc,2)=iuniq(2)
  582. 80 CONTINUE
  583. C Calcul des coordonnees barycentriques
  584. SEGINI xkon,ipxkon
  585. DO i=1,kon(/2)
  586. ipcour=0
  587. jp=idcp(i)
  588. ipcini=ipcour
  589. 58 CONTINUE
  590. DO ipc=ipcini+1,ipcour
  591. ipxkon(ipc,i)=0
  592. xkon(ipc,i)=0.
  593. ENDDO
  594. ipcour=ipcini
  595. DO 59 j=1,ikon(i)
  596. C Creation d'un angle
  597. IF (kon(j,i,1).EQ.0) GOTO 59
  598. ip1=kon(j,i,1)
  599. ip2=kon(j,i,2)
  600. DO 60 k=j+1,ikon(i)
  601. IF (kon(k,i,1).EQ.0) GOTO 60
  602. jp1=0
  603. IF (ip1.EQ.kon(k,i,1)) THEN
  604. jp1=ip1
  605. jp2=ip2
  606. jp3=kon(k,i,2)
  607. ELSE IF (ip1.EQ.kon(k,i,2)) THEN
  608. jp1=ip1
  609. jp2=ip2
  610. jp3=kon(k,i,1)
  611. ELSE IF (ip2.EQ.kon(k,i,1)) THEN
  612. jp1=ip2
  613. jp2=ip1
  614. jp3=kon(k,i,2)
  615. ELSE IF (ip2.EQ.kon(k,i,2)) THEN
  616. jp1=ip2
  617. jp2=ip1
  618. jp3=kon(k,i,1)
  619. ELSE
  620. GOTO 60
  621. ENDIF
  622. * jp le point en cours de travail
  623. * jp1 le sommet jp2 et jp3 les deux autres points
  624. xp=xcoor((jp-1)*idimp1+1)
  625. yp=xcoor((jp-1)*idimp1+2)
  626. zp=0.
  627. IF (IDIM.EQ.3) zp=xcoor((jp-1)*idimp1+3)
  628. xp1=xcoor((jp1-1)*idimp1+1)
  629. yp1=xcoor((jp1-1)*idimp1+2)
  630. zp1=0.
  631. IF (IDIM.EQ.3) zp1=xcoor((jp1-1)*idimp1+3)
  632. xp2=xcoor((jp2-1)*idimp1+1)
  633. yp2=xcoor((jp2-1)*idimp1+2)
  634. zp2=0.
  635. IF (IDIM.EQ.3) zp2=xcoor((jp2-1)*idimp1+3)
  636. xp3=xcoor((jp3-1)*idimp1+1)
  637. yp3=xcoor((jp3-1)*idimp1+2)
  638. zp3=0.
  639. IF (IDIM.EQ.3) zp3=xcoor((jp3-1)*idimp1+3)
  640. * le sinus (trafique)
  641. xv1=xp2-xp1
  642. yv1=yp2-yp1
  643. zv1=zp2-zp1
  644. aa=SQRT(xv1*xv1+yv1*yv1+zv1*zv1)
  645. xv1=xv1/aa
  646. yv1=yv1/aa
  647. zv1=zv1/aa
  648. xv2=xp3-xp1
  649. yv2=yp3-yp1
  650. zv2=zp3-zp1
  651. aa=SQRT(xv2*xv2+yv2*yv2+zv2*zv2)
  652. xv2=xv2/aa
  653. yv2=yv2/aa
  654. zv2=zv2/aa
  655. xv=yv1*zv2-zv1*yv2
  656. yv=zv1*xv2-xv1*zv2
  657. zv=xv1*yv2-yv1*xv2
  658. * provisoirement on prend le signe positif
  659. sinp=MAX(1.D-15,SQRT(xv*xv+yv*yv+zv*zv))
  660. * recherche du signe comparaison avec jp jp1 vect jp2 jp3
  661. xw1=xp-xp1
  662. yw1=yp-yp1
  663. zw1=zp-zp1
  664. aa=SQRT(xw1*xw1+yw1*yw1+zw1*zw1)
  665. xw1=xw1/aa
  666. yw1=yw1/aa
  667. zw1=zw1/aa
  668. xw2=xp3-xp2
  669. yw2=yp3-yp2
  670. zw2=zp3-zp2
  671. aa=SQRT(xw2*xw2+yw2*yw2+zw2*zw2)
  672. xw2=xw2/aa
  673. yw2=yw2/aa
  674. zw2=zw2/aa
  675. xw=yw1*zw2-zw1*yw2
  676. yw=zw1*xw2-xw1*zw2
  677. zw=xw1*yw2-yw1*xw2
  678. IF (iregu.NE.2) THEN
  679. IF ((xv*xw+yv*yw+zv*zw).LT.0.) sinp=-sinp
  680. ELSE
  681. IF ((xv*xw+yv*yw+zv*zw).LT.-1.D-15) then
  682. * Essai si le signe est negatif on supprime le coin pour ne travailler
  683. * que dans un convexe
  684. kon(k,i,1)=0
  685. kon(k,i,2)=0
  686. kon(j,i,1)=jp2
  687. kon(j,i,2)=jp3
  688. GOTO 58
  689. ENDIF
  690. ENDIF
  691. ipcour=ipcour+1
  692. ipxkon(ipcour,i)=jp1
  693. * les distances aux cotes
  694. xv=xp-xp1
  695. yv=yp-yp1
  696. zv=zp-zp1
  697. scal=xv*xv1+yv*yv1+zv*zv1
  698. xv=xv-scal*xv1
  699. yv=yv-scal*yv1
  700. zv=zv-scal*zv1
  701. dis1=MAX(1.D-20,SQRT(xv*xv+yv*yv+zv*zv))
  702. xv=xp-xp1
  703. yv=yp-yp1
  704. zv=zp-zp1
  705. scal=xv*xv2+yv*yv2+zv*zv2
  706. xv=xv-scal*xv2
  707. yv=yv-scal*yv2
  708. zv=zv-scal*zv2
  709. dis2=MAX(1.D-20,SQRT(xv*xv+yv*yv+zv*zv))
  710. C la coordonnee barycentrique :
  711. xkon(ipcour,i)=sinp/(dis1*dis2)
  712. IF (iregu.EQ.1) xkon(ipcour,i)=1.
  713. 60 CONTINUE
  714. 59 CONTINUE
  715. * on normalise
  716. skon=0.
  717. DO j=1,ikon(i)+1
  718. skon=skon+xkon(j,i)
  719. ENDDO
  720. DO j=1,ikon(i)+1
  721. xkon(j,i)=xkon(j,i)/skon
  722. ENDDO
  723. ENDDO
  724. SEGSUP kon,ikon
  725. * OK on a les coefficients barycentriques.
  726. * On duplique le maillage puis on cree les nouveaux noeuds
  727. * Mais d'abord voir les points imposes
  728. SEGACT ipt2,ipt3
  729. SEGINI icp2
  730. ipt7=ipt2
  731. ipt8=ipt3
  732. DO isous=1,MAX(1,ipt2.lisous(/1))
  733. IF (ipt2.lisous(/1).NE.0) ipt7=ipt2.lisous(isous)
  734. IF (ipt3.lisous(/1).NE.0) ipt8=ipt3.lisous(isous)
  735. SEGACT ipt7,ipt8
  736. DO i=1,ipt7.num(/1)
  737. DO j=1,ipt7.num(/2)
  738. ip=ipt7.num(i,j)
  739. IF (icp2(ip).EQ.0) icp2(ip)=ipt8.num(i,j)
  740. ENDDO
  741. ENDDO
  742. SEGDES ipt7,ipt8
  743. ENDDO
  744. SEGDES ipt2,ipt3
  745. SEGINI icorr
  746. DO 210 i=1,icpr(/1)
  747. IF (icpr(i).EQ.0) GOTO 210
  748. IF (icp2(i).EQ.0) THEN
  749. nbpts=xcoor(/1)/idimp1+1
  750. SEGADJ mcoord
  751. icorr(icpr(i))=nbpts
  752. DO j=1,idimp1
  753. xcoor((nbpts-1)*idimp1+j)=xcoor((i-1)*idimp1+j)
  754. ENDDO
  755. ELSE
  756. icorr(icpr(i))=icp2(i)
  757. ENDIF
  758. 210 CONTINUE
  759. * On duplique les maillages
  760. SEGINI,ipt4=ipt1
  761. nbref=0
  762. nbsous=ipt4.lisous(/1)
  763. nbnn=ipt4.num(/1)
  764. nbelem=ipt4.num(/2)
  765. SEGADJ ipt4
  766. ipt7=ipt1
  767. ipt8=ipt4
  768. DO isous=1,MAX(1,ipt1.lisous(/1))
  769. IF (ipt1.lisous(/1).NE.0) THEN
  770. ipt7=ipt1.lisous(isous)
  771. SEGINI,ipt8=ipt7
  772. ipt4.lisous(isous)=ipt8
  773. ENDIF
  774. DO i=1,ipt8.num(/1)
  775. DO j=1,ipt8.num(/2)
  776. ipt8.num(i,j)=icorr(icpr(ipt8.num(i,j)))
  777. ENDDO
  778. ENDDO
  779. SEGDES ipt7,ipt8
  780. ENDDO
  781. SEGDES ipt1,ipt4
  782. C Maintenant on peut deplacer les points
  783. nfois=50
  784. CALL LIRENT(nfois,0,iretou)
  785. DO ifois=1,nfois
  786. DO i=1,numnp
  787. IF (icp2(idcp(i)).EQ.0) THEN
  788. DO j=1,IDIM
  789. xaux=0.
  790. DO k=1,nkon+1
  791. IF (ipxkon(k,i).NE.0) xaux=xaux+xkon(k,i)
  792. . *xcoor((icorr(icpr(ipxkon(k,i)))-1)*idimp1+j)
  793. ENDDO
  794. xcoor((icorr(i)-1)*idimp1+j)=xaux
  795. ENDDO
  796. ENDIF
  797. ENDDO
  798. ENDDO
  799.  
  800. IF (ichpo.EQ.0) THEN
  801. CALL ECROBJ('MAILLAGE',ipt4)
  802. ELSE
  803. C Sortie d'un CHPOINT, les points concernes sont reperes dans icpr
  804. SEGSUP ipt2,ipt3
  805. nbelem=icorr(/1)
  806. nbnn=1
  807. nbsous=0
  808. nbref=0
  809. SEGINI ipt2
  810. ipt2.itypel=1
  811. nat=1
  812. nsoupo=1
  813. SEGINI mchpo1
  814. mchpo1.jattri(1)=jatan
  815. IF (IDIM.EQ.1) nc=1
  816. IF (IDIM.EQ.2) nc=2
  817. IF (IDIM.EQ.3) nc=3
  818. SEGINI msoupo
  819. mchpo1.ipchp(1)=msoupo
  820. igeoc=ipt2
  821. n=nbelem
  822. SEGINI mpoval
  823. ipoval=mpoval
  824. IF (IDIM.EQ.1) THEN
  825. IF (NIFOUR.LE.9) THEN
  826. nocomp(1)='UX'
  827. ELSE
  828. nocomp(1)='UR'
  829. ENDIF
  830. ELSE IF (IDIM.EQ.2) THEN
  831. IF (IFOUR.GE.0) THEN
  832. nocomp(1)='UR'
  833. nocomp(2)='UZ'
  834. ELSE
  835. nocomp(1)='UX'
  836. nocomp(2)='UY'
  837. ENDIF
  838. ELSE IF (IDIM.EQ.3) THEN
  839. nocomp(1)='UX'
  840. nocomp(2)='UY'
  841. nocomp(3)='UZ'
  842. ENDIF
  843. SEGDES msoupo
  844. DO i=1,idcp(/1)
  845. ipt2.num(1,i)=idcp(i)
  846. ifi=(icorr(i)-1)*idimp1
  847. ide=(idcp(i)-1)*idimp1
  848. DO j=1,IDIM
  849. vpocha(i,j)=xcoor(ifi+j)-xcoor(ide+j)
  850. ENDDO
  851. ENDDO
  852. SEGDES mpoval,ipt2,mchpo1
  853. CALL ECROBJ('CHPOINT ',mchpo1)
  854. ENDIF
  855. SEGSUP icpr,idcp,xkon,ipxkon,icp2,icorr
  856.  
  857. RETURN
  858. END
  859.  
  860.  
  861.  
  862.  
  863.  
  864.  

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