Télécharger dedu.eso

Retour à la liste

Numérotation des lignes :

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

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