Télécharger dedu.eso

Retour à la liste

Numérotation des lignes :

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

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