Télécharger ligne.eso

Retour à la liste

Numérotation des lignes :

ligne
  1. C LIGNE SOURCE CB215821 23/10/11 21:15:03 11749
  2.  
  3. C-----------------------------------------------------------------------
  4. C FABRIQUE DES LIGNES CONSTITUES D ELEM A 2 OU 3 NOEUDS
  5. C
  6. C ILIGGI indique quel est le type de la ligne support :
  7. C 1 = SEGMENT DE DROITE (operateur DROI)
  8. C 2 =
  9. C 3 = ARC DE CERCLE (operateur CERC)
  10. C 4 =
  11. C 5 = <-- indisponible
  12. C 6 = <-- indisponible
  13. C 7 = ARC DE PARABOLE (operateur PARA)
  14. C 8 = <-- indisponible : ELLIPSE (operateur ELLI)
  15. C 9 = PORTION DE CUBIQUE (operateur CUBP)
  16. C 10 = PORTION DE CUBIQUE (operateur CUBT)
  17. C 11 = ARC DE CERCLE (operateur CER3)
  18. C 12 =
  19. C
  20. C 09/2003 :
  21. C Modifications suite a la mise en place du cas IDIM=1. Ajout de
  22. C tests de coherence. Indisponibilite des operateurs PARA,CERC,CER3,
  23. C CUBP,CUBT dans le cas IDIM=1. Ajout de commentaires.
  24. C
  25. C BP (2017-03-27) : 3 facons de creer un CERCle
  26. C -CENT (option par defaut de CERC) = via pt init + centre + pt fin
  27. C -PASS (ancien CER3) = via pt ini + pt entre ini et fin + pt fin
  28. C -ROTA (nouveau inspiré de ROTA) = via pt ini + 2 points definissant
  29. C l'axe de rotation + angle
  30. C
  31. C-----------------------------------------------------------------------
  32.  
  33. SUBROUTINE LIGNE(ILIGGI,IFONC,DEN1OX,DEN2OX,INBR)
  34.  
  35. IMPLICIT INTEGER(I-N)
  36. IMPLICIT REAL*8 (A-H,O-Z)
  37.  
  38. -INC CCREEL
  39.  
  40. -INC PPARAM
  41. -INC CCOPTIO
  42. -INC CCGEOME
  43. -INC SMELEME
  44. -INC SMCOORD
  45.  
  46. PARAMETER (NITER=5)
  47. logical ltelq
  48. SEGMENT TABPAR(NBELEM)
  49. SEGMENT / STRAV / (DIST2(NBNO+1),XPOIN(NBNO+1),YPOIN(NBNO+1),
  50. . ZPOIN(NBNO+1),DIST(NBNO+1),U(NSTRAV,NBNO+1))
  51.  
  52. C :::::: AJOUT PIFFARD : 16/01/86 POUR CUBIQUES ET CER3 ::::::::
  53. DIMENSION NDP(12)
  54. DIMENSION TU(4),CU(4,4),DLONGT(NITER+1),DMOY(NITER)
  55. DIMENSION Q(4)
  56. CHARACTER*(4) ITTEMP*8
  57. CHARACTER*4 MCLE(2),MCLE2(3)
  58. CHARACTER*4 MUNIF(1),MELIM(2)
  59. REAL*8 NVECT
  60. c nombre de points a lire
  61. c 2 pour DROI, 3 pour CERC PARA et CER3, 4 pour CUBP et CUBT
  62. DATA NDP / 2,2,3,3,0,0,3,0,4,4,3,3 /
  63. DATA MCLE / 'DINI','DFIN' /
  64. DATA MUNIF / 'UNIF' /
  65. DATA MCLE2 / 'CENT','PASS','ROTA' /
  66. DATA MELIM / 'NOEL','ELIM' /
  67. C :::::: AJOUT PIFFARD : 16/01/86 POUR CUBIQUES ::::::::
  68. DATA CU / 1.D0, 0.D0, -5.5D0, -1.0D0,
  69. . 0.D0, 0.D0, 9.0D0, 4.5D0,
  70. . 0.D0, 0.D0, -4.5D0, -9.0D0,
  71. . 0.D0, 1.D0, 1.0D0, 5.5D0 /
  72.  
  73. C-----------------------------------------------------------------------
  74. C INITIALISATIONS ET VERIFICATIONS
  75. C-----------------------------------------------------------------------
  76. XUN =1.D0
  77. XDEUX=2.D0
  78. XDEMI=0.5D0
  79. XQUAR=0.25D0
  80. X180 =180.D0
  81.  
  82. C Initialisations necessaires sinon SIGFPE en DEBUG...
  83. XINI =XZERO
  84. XINT1=XZERO
  85. XINT2=XZERO
  86. XFIN =XZERO
  87. YINI =XZERO
  88. YINT1=XZERO
  89. YINT2=XZERO
  90. YFIN =XZERO
  91. ZINI =XZERO
  92. ZINT1=XZERO
  93. ZINT2=XZERO
  94. ZFIN =XZERO
  95. F1 =XZERO
  96. F2 =XZERO
  97. F3 =XZERO
  98. F4 =XZERO
  99. segact mcoord*mod
  100.  
  101. C Recup des donnees en entree
  102. ILIG=ILIGGI
  103. DEN1=DEN1OX
  104. DEN2=DEN2OX
  105.  
  106. C En DIMEnsion 1, seul l'operateur DROIte peut etre utilise
  107. IF ((IDIM.EQ.1).AND.(ILIG.NE.1)) THEN
  108. INTERR(1)=IDIM
  109. CALL ERREUR(709)
  110. RETURN
  111. ENDIF
  112.  
  113. C Quelques verifications
  114. IF (ILCOUR.EQ.0) THEN
  115. CALL ERREUR(16)
  116. RETURN
  117. ENDIF
  118. C INOMB = nombre de "points" a lire pour definir la ligne
  119. INOMB=NDP(ILIG)
  120. IF (INOMB.EQ.0) THEN
  121. CALL ERREUR(19)
  122. RETURN
  123. ENDIF
  124. ITYPELT=KDEGRE(ILCOUR)
  125. C Pas de verif. sur ITYPELT (2 ou 3 uniquement)
  126. IF (ITYPELT.EQ.0) THEN
  127. CALL ERREUR(16)
  128. RETURN
  129. ENDIF
  130. C Arc de CUBIQUE (CUBP et CUBT) uniquement en SEG2
  131. IF (((ILIG.EQ.9).OR.(ILIG.EQ.10)).AND.(ITYPELT.EQ.3)) THEN
  132. CALL ERREUR(16)
  133. RETURN
  134. ENDIF
  135. IF ((ILIG.EQ.12).AND.(ITYPELT.NE.3)) THEN
  136. CALL ERREUR(16)
  137. RETURN
  138. ENDIF
  139. NBNN=NBNNE(ITYPELT)
  140. IF (NBNN.EQ.0) THEN
  141. CALL ERREUR(16)
  142. RETURN
  143. ENDIF
  144. IF (IERR.NE.0) RETURN
  145.  
  146. C Initialisations dans le cas de CUBP et CUBT
  147. IF (ILIG.EQ.9.OR.ILIG.EQ.10) THEN
  148. DO i=1,NITER
  149. DLONGT(i)=XZERO
  150. DMOY(i) =XZERO
  151. ENDDO
  152. DLONGT(NITER+1)=XZERO
  153. ENDIF
  154.  
  155.  
  156. C-----------------------------------------------------------------------
  157. C LECTURE DES DONNEES FACULTATIVES (DECOUPAGE, DENSITES...)
  158. C-----------------------------------------------------------------------
  159.  
  160. C Usage en tant que fonction -> on saute
  161. IF (IFONC.EQ.0) GOTO 400
  162.  
  163. C Y a-t-il un facteur de decoupage ?
  164. IF (ILIG.NE.10) THEN
  165. INBR=0
  166. CALL MESLIR(-172)
  167. CALL LIRENT(INBR,0,IRETOU)
  168. ELSE
  169. C Si CUBIque a vecteurs TGTS, il FAUT un facteur de decoupage.
  170. CALL MESLIR(-172)
  171. CALL LIRENT(INBR,1,IRETOU)
  172. IF (IERR.NE.0) RETURN
  173. ENDIF
  174.  
  175. C Y a t-il des densites imposees ?
  176. IMPOI=0
  177. IMPOF=0
  178. 413 CALL MESLIR(-171)
  179. CALL LIRMOT(MCLE,2,IRETOU,0)
  180. IF (IRETOU.EQ.1) THEN
  181. CALL MESLIR(-170)
  182. CALL LIRREE(DEN1,1,IRETOU)
  183. IF (IERR.NE.0) RETURN
  184. IMPOI=1
  185. GOTO 413
  186. ELSEIF (IRETOU.EQ.2) THEN
  187. CALL MESLIR(-169)
  188. CALL LIRREE(DEN2,1,IRETOU)
  189. IF (IERR.NE.0) RETURN
  190. IMPOF=1
  191. GOTO 413
  192. ENDIF
  193.  
  194. 400 CONTINUE
  195.  
  196.  
  197. C-----------------------------------------------------------------------
  198. C LECTURE DES DONNEES OBLIGATOIRES (MAILLAGES, POINTS ...)
  199. C-----------------------------------------------------------------------
  200.  
  201. C Y-a-t-il le mot UNIF pour les cubiques ?
  202. IUNIF=0
  203. IF (ILIG.EQ.9) THEN
  204. CALL LIRMOT(MUNIF,1,IUNIF,0)
  205. ENDIF
  206.  
  207. C Y-a-t-il un mot cle ('CERC', 'PASS' ou 'ROTA') pour les CERCles ?
  208. ICLE2=0
  209. IELIM=0
  210. IF (ILIG.EQ.3) THEN
  211. CALL LIRMOT(MCLE2,3,ICLE2,0)
  212. c cas CERC ROTA : il faut obligatoirement un angle
  213. c IF(ICLE2.NE.0) WRITE(IOIMP,*) 'MOT CLE LU =',MCLE2(ICLE2)
  214. IF(ICLE2.EQ.3) THEN
  215. CALL MESLIR(-235)
  216. CALL LIRREE(ANG,1,IRETOU)
  217. IF (IERR.NE.0) RETURN
  218. c si angle=360 degres, on veut probablement eliminer le dernier point
  219. IF (A_EGALE_B(ANG,360.D0)) THEN
  220. CALL LIRMOT(MELIM,2,IELIM,0)
  221. ENDIF
  222. ANG=ANG*XPI/X180
  223. IF(IDIM.LE.2) INOMB=INOMB-1
  224. ENDIF
  225. ENDIF
  226.  
  227. C Quel est le type de la premiere donnee (POINT OU LIGNE) ?
  228. C IPOIN1 = premier POINT de la ligne a construire
  229. IPOIN3=0
  230. ITP1=0
  231. ITP2=0
  232. ITTEMP=' '
  233. CALL MESLIR(-168)
  234. CALL QUETYP(ITTEMP,1,IRETOU)
  235. IF (IERR.NE.0) RETURN
  236. CALL MESLIR(-168)
  237. CALL LIROBJ(ITTEMP,IPT1,1,IRETOU)
  238. CALL ACTOBJ(ITTEMP,IPT1,1)
  239. IF (ITTEMP.EQ.'MAILLAGE') THEN
  240. SEGACT IPT1
  241. CALL EXTRPO(IPT1,1,IPOIN1)
  242. IF(IERR .NE. 0)RETURN
  243. CALL EXTRPO(IPT1,2,IPOIN3)
  244. IF(IERR .NE. 0)RETURN
  245. ITP1=1
  246. ELSE
  247. IF (ITTEMP.NE.'POINT ') THEN
  248. MOTERR(1:8)=ITTEMP
  249. CALL ERREUR(39)
  250. RETURN
  251. ENDIF
  252. IPOIN1=IPT1
  253. IPOIN3=IPT1
  254. ENDIF
  255.  
  256. C Lecture des troisieme et quatrieme POINTs si necessaire
  257. IF (INOMB.EQ.2) GOTO 503
  258. CALL MESLIR(-167)
  259. CALL LIROBJ('POINT ',IPCEN,1,IRETOU)
  260. IF (IERR.NE.0) RETURN
  261. IF (INOMB.EQ.3) GOTO 503
  262. CALL LIROBJ('POINT ',IPCEN2,1,IRETOU)
  263. IF (IERR.NE.0) RETURN
  264.  
  265. C Lecture eventuelle du second POINT ou LIGNE
  266. C IPOIN2 = POINT final de la ligne a construire
  267. 503 CALL MESLIR(-166)
  268. CALL LIROBJ('MAILLAGE',IP2,0,IRETOU)
  269. IF (IRETOU.EQ.1) THEN
  270. CALL ACTOBJ('MAILLAGE',IP2,1)
  271. IPT6=IP2
  272. CALL EXTRPO(IP2,2,IPOIN2)
  273. IF(IERR .NE. 0)RETURN
  274. ITP2=1
  275. ELSE
  276. CALL MESLIR(-166)
  277. CALL LIROBJ('POINT ',IPOIN2,0,IRETOU)
  278. IF (IRETOU.EQ.0) THEN
  279. IPOIN2=IPOIN3
  280. IF (IPOIN1.EQ.IPOIN2) CALL ERREUR(20)
  281. ENDIF
  282. IF (IPOIN1.EQ.IPOIN2) CALL ERREUR(432)
  283. ENDIF
  284.  
  285. C Affectation des densites suivant les donnees
  286. IBP1=(IDIM+1)*IPOIN1
  287. IBP2=(IDIM+1)*IPOIN2
  288. IF (IFONC.NE.0) THEN
  289. IF (IMPOI.EQ.0) DEN1=XCOOR(IBP1)
  290. IF (IMPOF.EQ.0) DEN2=XCOOR(IBP2)
  291. ENDIF
  292. IF ((INBR.LE.0).AND.((DEN1*DEN2).LE.XZERO)) CALL ERREUR(17)
  293. IF (IERR.NE.0) RETURN
  294. * A quoi ça sert ?
  295. DENI=DEN1
  296. DECA=DEN2-DEN1
  297.  
  298.  
  299. C-----------------------------------------------------------------------
  300. C BRANCHEMENT SUIVANT LE TYPE DE LIGNE
  301. C-----------------------------------------------------------------------
  302.  
  303. GOTO (11,11,13,13,11,11,13,17,18,18,19,13),ILIG
  304.  
  305.  
  306. C --------------------------
  307. C Segment de DROITE (DROI)
  308. C --------------------------
  309.  
  310. 11 XDIS=XCOOR(IBP2-IDIM)-XCOOR(IBP1-IDIM)
  311. IF (IDIM.GE.2) THEN
  312. YDIS=XCOOR(IBP2-IDIM+1)-XCOOR(IBP1-IDIM+1)
  313. ZDIS=XZERO
  314. IF (IDIM.GE.3) ZDIS=XCOOR(IBP2-IDIM+2)-XCOOR(IBP1-IDIM+2)
  315. DIS=SQRT(MAX(XDIS**2+YDIS**2+ZDIS**2,XZERO))
  316. ELSE
  317. DIS=ABS(XDIS)
  318. ENDIF
  319. C Modif. pour permettre de faire des montagnes
  320. C Si decoupage impose, on accepte points confondus
  321. IF ((INBR.NE.0).AND.(DIS.EQ.XZERO)) DIS=XUN
  322. IF (DIS.EQ.XZERO) CALL ERREUR(21)
  323. IF (IERR.NE.0) RETURN
  324. DLONG=DIS
  325. IF (ABS(DLONG).LE.XZPREC*ABS(DEN1)) THEN
  326. CALL ERREUR(21)
  327. RETURN
  328. ENDIF
  329. DEN1=DEN1/DIS
  330. DEN2=DEN2/DIS
  331. XIN =XCOOR(IBP1-IDIM)
  332. IF (IDIM.EQ.1) THEN
  333. YIN=XZERO
  334. ZIN=XZERO
  335. ELSE IF (IDIM.EQ.2) THEN
  336. YIN=XCOOR(IBP1-IDIM+1)
  337. ZIN=XZERO
  338. ELSE IF (IDIM.EQ.3) THEN
  339. YIN=XCOOR(IBP1-IDIM+1)
  340. ZIN=XCOOR(IBP1-IDIM+2)
  341. ENDIF
  342. GOTO 200
  343.  
  344.  
  345. C --------------------------------------------
  346. C Arc de CERCLE (CERC) ou de PARABOLE (PARA)
  347. C --------------------------------------------
  348.  
  349. 13 CONTINUE
  350. C CERC .. 'PASS' .. (<=> CER3 ..) --> GOTO 19
  351. IF(ICLE2.EQ.2) GOTO 19
  352. C CERC .. 'ROTA' .. --> GOTO 20
  353. IF(ICLE2.EQ.3) GOTO 20
  354.  
  355. IBPC=(IDIM+1)*IPCEN
  356. XCEN=XCOOR(IBPC-IDIM)
  357. YCEN=XCOOR(IBPC-IDIM+1)
  358. ZCEN=XZERO
  359. IF (IDIM.GE.3) ZCEN=XCOOR(IBPC-IDIM+2)
  360. C RECALCUL DU CENTRE DANS LE CAS DE L'ARC DE CERCLE
  361. C On ne recalcule pas le centre des paraboles
  362. IF (ILIG.NE.7) THEN
  363. XVP= XCOOR(IBP2-IDIM) -XCOOR(IBP1-IDIM)
  364. XMP=(XCOOR(IBP2-IDIM) +XCOOR(IBP1-IDIM))*XDEMI
  365. YVP= XCOOR(IBP2-IDIM+1)-XCOOR(IBP1-IDIM+1)
  366. YMP=(XCOOR(IBP2-IDIM+1)+XCOOR(IBP1-IDIM+1))*XDEMI
  367. ZVP=XZERO
  368. ZMP=XZERO
  369. IF (IDIM.GE.3) THEN
  370. ZVP=XCOOR(IBP2-IDIM+2)-XCOOR(IBP1-IDIM+2)
  371. ZMP=(XCOOR(IBP2-IDIM+2)+XCOOR(IBP1-IDIM+2))*XDEMI
  372. ENDIF
  373. SNOR=SQRT(MAX(XVP**2+YVP**2+ZVP**2,XZERO))
  374. IF (SNOR.LT.XPETIT) THEN
  375. CALL ERREUR(21)
  376. RETURN
  377. ENDIF
  378. XVP=XVP/SNOR
  379. YVP=YVP/SNOR
  380. ZVP=ZVP/SNOR
  381. RO=(XMP-XCEN)*XVP+(YMP-YCEN)*YVP+(ZMP-ZCEN)*ZVP
  382. XCEN=XCEN+RO*XVP
  383. YCEN=YCEN+RO*YVP
  384. ZCEN=ZCEN+RO*ZVP
  385. ENDIF
  386. XV1=XCOOR(IBP1-IDIM)-XCEN
  387. XV2=XCOOR(IBP2-IDIM)-XCEN
  388. YV1=XCOOR(IBP1-IDIM+1)-YCEN
  389. YV2=XCOOR(IBP2-IDIM+1)-YCEN
  390. ZV1=XZERO
  391. ZV2=XZERO
  392. IF (IDIM.GE.3) THEN
  393. ZV1=XCOOR(IBP1-IDIM+2)-ZCEN
  394. ZV2=XCOOR(IBP2-IDIM+2)-ZCEN
  395. ENDIF
  396. IF (ILIG.NE.7) THEN
  397. C Arc de CERCLE
  398. c calcul du rayon moyen
  399. RAY1=SQRT(MAX(XV1**2+YV1**2+ZV1**2,XZERO))
  400. RAY2=SQRT(MAX(XV2**2+YV2**2+ZV2**2,XZERO))
  401. IF (RAY1.LT.(0.99D0*RAY2) .OR.
  402. & RAY1.GT.(1.01D0*RAY2)) CALL ERREUR(21)
  403. IF (IERR.NE.0) RETURN
  404. RAY=SQRT(MAX(RAY1*RAY2,XZERO))
  405. c produit scalaire V1*V2 = |V1|.|V2|.cos(ANG)
  406. SCAL=XV1*XV2+YV1*YV2+ZV1*ZV2
  407. c produit vectoriel V1 PVEC V2 = |V1|.|V2|.sin(ANG)
  408. XVECT=YV1*ZV2-YV2*ZV1
  409. YVECT=XV2*ZV1-XV1*ZV2
  410. ZVECT=XV1*YV2-XV2*YV1
  411. VECT=SQRT(MAX(XVECT**2+YVECT**2+ZVECT**2,XZERO))
  412. ANG=ATAN2(VECT,SCAL)
  413. IF (IIMPI.EQ.1) WRITE (IOIMP,9143) ANG
  414. 9143 FORMAT(' ANGLE DE L ARC EN RADIAN ',G12.5)
  415. IF (ABS(ANG).GT.XPI) THEN
  416. CALL ERREUR(313)
  417. RETURN
  418. ENDIF
  419. SINA=SIN(ANG)
  420. c longueur de la ligne
  421. DLONG=ANG*RAY
  422. IF (ABS(DLONG).LE.XZPREC*ABS(DEN1)) THEN
  423. CALL ERREUR(21)
  424. RETURN
  425. ENDIF
  426. DEN1=DEN1/DLONG
  427. DEN2=DEN2/DLONG
  428. c vecteur de base orthogonal a V1 et dans le plan (V1,V2)
  429. XV3=(YVECT*ZV1-YV1*ZVECT)/(RAY*RAY*SINA)
  430. YV3=(XV1*ZVECT-XVECT*ZV1)/(RAY*RAY*SINA)
  431. ZV3=(XVECT*YV1-XV1*YVECT)/(RAY*RAY*SINA)
  432. ELSE
  433. C Arc de PARABOLE
  434. XIN=XCOOR(IBP1-IDIM)
  435. YIN=XCOOR(IBP1-IDIM+1)
  436. XDIS=XDEMI*(XCOOR(IBP2-IDIM)-XIN)
  437. YDIS=XDEMI*(XCOOR(IBP2-IDIM+1)-YIN)
  438. IF (IDIM.GE.3) THEN
  439. ZIN =XCOOR(IBP1-IDIM+2)
  440. ZDIS=XDEMI*(XCOOR(IBP2-IDIM+2)-ZIN)
  441. ELSE
  442. ZIN =XZERO
  443. ZDIS=XZERO
  444. ENDIF
  445. DIS=2*SQRT(XDIS*XDIS+YDIS*YDIS+ZDIS*ZDIS)
  446. XIN=XIN+XDIS
  447. YIN=YIN+YDIS
  448. ZIN=ZIN+ZDIS
  449. DLONG=XDEMI*(DIS+SQRT(XV1*XV1+YV1*YV1+ZV1*ZV1)
  450. . +SQRT(XV2*XV2+YV2*YV2+ZV2*ZV2))
  451. IF (ABS(DLONG).LE.XZPREC*ABS(DEN1)) THEN
  452. CALL ERREUR(21)
  453. RETURN
  454. ENDIF
  455. DEN1=DEN1/DLONG
  456. DEN2=DEN2/DLONG
  457. XV=-XQUAR*(XV1+XV2)
  458. YV=-XQUAR*(YV1+YV2)
  459. ZV=-XQUAR*(ZV1+ZV2)
  460. ENDIF
  461. GOTO 200
  462.  
  463.  
  464. C ----------------
  465. C NON disponible
  466. C ----------------
  467. 17 GOTO 200
  468.  
  469.  
  470. C ---------------------------
  471. C Arc de CUBIQUE (CUBP,CUBT)
  472. C ---------------------------
  473.  
  474. C ::::::::: AJOUT PIFFARD 16/01/86 : POUR CUBIQUES :::::
  475. C ...... RECUPERATION DES 4 PTS DE BASE SI CUBP .....
  476. C ...... DES 2 PTS DE BASE + TGTES SI CUBT .....
  477. 18 XINI=XCOOR(IBP1-IDIM)
  478. YINI=XCOOR(IBP1-IDIM+1)
  479. ZINI=XZERO
  480. IF (IDIM.GE.3) ZINI=XCOOR(IBP1-IDIM+2)
  481. IBP3=(IDIM+1)*IPCEN
  482. XINT1=XCOOR(IBP3-IDIM)
  483. YINT1=XCOOR(IBP3-IDIM+1)
  484. ZINT1=XZERO
  485. IF (IDIM.GE.3) ZINT1=XCOOR(IBP3-IDIM+2)
  486. IBP4=(IDIM+1)*IPCEN2
  487. XINT2=XCOOR(IBP4-IDIM)
  488. YINT2=XCOOR(IBP4-IDIM+1)
  489. ZINT2=XZERO
  490. IF (IDIM.GE.3) ZINT2=XCOOR(IBP4-IDIM+2)
  491. XFIN=XCOOR(IBP2-IDIM)
  492. YFIN=XCOOR(IBP2-IDIM+1)
  493. ZFIN=XZERO
  494. IF (IDIM.GE.3) ZFIN=XCOOR(IBP2-IDIM+2)
  495. C Calcul de la distance entre les 2 ou 4 pts pour determiner
  496. C le nombre de segments a creer si il n'est pas donne
  497. ZDIS=XZERO
  498. IF (ILIG.EQ.9) THEN
  499. XDIS=XFIN-XINT2-XINT1-XINI
  500. YDIS=YFIN-YINT2-YINT1-YINI
  501. IF (IDIM.GE.3) ZDIS=ZFIN-ZINT2-ZINT1-ZINI
  502. ELSE
  503. XDIS=XFIN-XINI
  504. YDIS=YFIN-YINI
  505. IF (IDIM.GE.3) ZDIS=ZFIN-ZINI
  506. ENDIF
  507. DIS=SQRT(XDIS*XDIS+YDIS*YDIS+ZDIS*ZDIS)
  508. IF (DIS.EQ.0) CALL ERREUR(21)
  509. IF (IERR.NE.0) RETURN
  510. DLONG=DIS
  511. DEN1=DEN1/DIS
  512. DEN2=DEN2/DIS
  513. GOTO 200
  514.  
  515.  
  516. C -----------------------
  517. C Arc de CERCLE (CERC3)
  518. C -----------------------
  519.  
  520. C Ajout PIFFARD : JANV 86 POUR CERCLE PASSANT PAR 3 PTS
  521. C Remanie par VERPEAUX on se branche sur le cercle normal
  522. 19 IBPC=(IDIM+1)*IPCEN
  523. XOM1=XCOOR(IBP1-IDIM)
  524. YOM1=XCOOR(IBP1-IDIM+1)
  525. ZOM1=XZERO
  526. IF (IDIM.GE.3) ZOM1=XCOOR(IBP1-IDIM+2)
  527. XOM2=XCOOR(IBPC-IDIM)
  528. YOM2=XCOOR(IBPC-IDIM+1)
  529. ZOM2=XZERO
  530. IF (IDIM.GE.3) ZOM2=XCOOR(IBPC-IDIM+2)
  531. XOM3=XCOOR(IBP2-IDIM)
  532. YOM3=XCOOR(IBP2-IDIM+1)
  533. ZOM3=XZERO
  534. IF (IDIM.GE.3) ZOM3=XCOOR(IBP2-IDIM+2)
  535. XM1=XOM2-XOM1
  536. YM1=YOM2-YOM1
  537. ZM1=ZOM2-ZOM1
  538. SNOR=SQRT(XM1*XM1+YM1*YM1+ZM1*ZM1)
  539. XM1=XM1/SNOR
  540. YM1=YM1/SNOR
  541. ZM1=ZM1/SNOR
  542. XM2=XOM3-XOM2
  543. YM2=YOM3-YOM2
  544. ZM2=ZOM3-ZOM2
  545. SNOR=SQRT(XM2*XM2+YM2*YM2+ZM2*ZM2)
  546. XM2=XM2/SNOR
  547. YM2=YM2/SNOR
  548. ZM2=ZM2/SNOR
  549. XVECT=YM1*ZM2-YM2*ZM1
  550. YVECT=ZM1*XM2-ZM2*XM1
  551. ZVECT=XM1*YM2-XM2*YM1
  552. VECT=SQRT(XVECT*XVECT+YVECT*YVECT+ZVECT*ZVECT)
  553. XVECT=XVECT/VECT
  554. YVECT=YVECT/VECT
  555. ZVECT=ZVECT/VECT
  556. XV1=YM1*ZVECT-YVECT*ZM1
  557. YV1=ZM1*XVECT-ZVECT*XM1
  558. ZV1=XM1*YVECT-XVECT*YM1
  559. XM3=XOM3-XOM1
  560. YM3=YOM3-YOM1
  561. ZM3=ZOM3-ZOM1
  562. RO=(XM2*XM3+YM2*YM3+ZM2*ZM3)/(XM2*XV1+YM2*YV1+ZM2*ZV1)
  563. XCEN=XDEMI*(XOM1+XOM2+RO*XV1)
  564. YCEN=XDEMI*(YOM1+YOM2+RO*YV1)
  565. ZCEN=XDEMI*(ZOM1+ZOM2+RO*ZV1)
  566. IF (IIMPI.EQ.1) WRITE(IOIMP,6032) XCEN,YCEN,ZCEN
  567. 6032 FORMAT(2X,'CENTRE DU CERCLE =',3G13.5)
  568. RAY=SQRT((XOM1-XCEN)**2+(YOM1-YCEN)**2+(ZOM1-ZCEN)**2)
  569. IF (IIMPI.EQ.1) WRITE(IOIMP,6033) RAY
  570. 6033 FORMAT(2X,'RAYON DU CERCLE =',G13.5)
  571. XV1=XCOOR(IBP1-IDIM)-XCEN
  572. XV2=XCOOR(IBP2-IDIM)-XCEN
  573. YV1=XCOOR(IBP1-IDIM+1)-YCEN
  574. YV2=XCOOR(IBP2-IDIM+1)-YCEN
  575. ZV1=XZERO
  576. ZV2=XZERO
  577. IF (IDIM.GE.3) THEN
  578. ZV1=XCOOR(IBP1-IDIM+2)-ZCEN
  579. ZV2=XCOOR(IBP2-IDIM+2)-ZCEN
  580. ENDIF
  581. SCAL=XV1*XV2+YV1*YV2+ZV1*ZV2
  582. XVE=YV1*ZV2-YV2*ZV1
  583. YVE=XV2*ZV1-XV1*ZV2
  584. ZVE=XV1*YV2-XV2*YV1
  585. VE=XVE*XVECT+YVE*YVECT+ZVE*ZVECT
  586. ANG=ATAN2(VE,SCAL)
  587. IF (ANG.LT.XZERO) ANG=ANG+(XDEUX*XPI)
  588. IF (IIMPI.EQ.1) WRITE(IOIMP,9143) ANG
  589. SINA=SIN(ANG)
  590. DLONG=ANG*RAY
  591. IF (ABS(DLONG).LE.XZPREC*ABS(DEN1)) THEN
  592. CALL ERREUR(21)
  593. RETURN
  594. ENDIF
  595. DEN1=DEN1/DLONG
  596. DEN2=DEN2/DLONG
  597. XV3=YVECT*ZV1-YV1*ZVECT
  598. YV3=XV1*ZVECT-XVECT*ZV1
  599. ZV3=XVECT*YV1-XV1*YVECT
  600. ILIG=3
  601. GOTO 200
  602.  
  603.  
  604. C -------------------------------------------------------------
  605. C Arc de CERCLE defini par 1 point + axe + ANGle (CERC 'ROTA')
  606. C -------------------------------------------------------------
  607.  
  608. 20 CONTINUE
  609. c Connus : ANG , points : IPOIN1 (IBP1), IPCEN et IPOIN2(IBP2)
  610. c Reste a calculer : points : vrai CENtre, V1, V3, DLONG
  611. c PC = 1 POINT DE L'AXE
  612. IF(IDIM.LE.2) IPCEN=IPOIN2
  613. IBPC=(IDIM+1)*IPCEN
  614. XC=XCOOR(IBPC-IDIM)
  615. YC=XCOOR(IBPC-IDIM+1)
  616. ZC=XZERO
  617. IF (IDIM.GE.3) ZC=XCOOR(IBPC-IDIM+2)
  618. C CALCUL DE L'AXE DE ROTATION
  619. IF (IDIM.LE.2) THEN
  620. XVECT=XZERO
  621. YVECT=XZERO
  622. ZVECT=XUN
  623. ELSE
  624. XVECT=XCOOR(IBP2-IDIM)-XC
  625. YVECT=XCOOR(IBP2-IDIM+1)-YC
  626. ZVECT=XCOOR(IBP2-IDIM+2)-ZC
  627. NVECT=SQRT(XVECT**2+YVECT**2+ZVECT**2)
  628. IF (NVECT.LE.XPETIT) THEN
  629. CALL ERREUR(277)
  630. RETURN
  631. ENDIF
  632. XVECT=XVECT/NVECT
  633. YVECT=YVECT/NVECT
  634. ZVECT=ZVECT/NVECT
  635. ENDIF
  636. C CALCUL DU VRAI CENTRE = PROJECTION DE P1 SUR L'AXE
  637. XV1=XCOOR(IBP1-IDIM)
  638. YV1=XCOOR(IBP1-IDIM+1)
  639. ZV1=XZERO
  640. IF (IDIM.GE.3) ZV1=XCOOR(IBP1-IDIM+2)
  641. XC1=XV1-XC
  642. YC1=YV1-YC
  643. ZC1=XZERO
  644. IF (IDIM.GE.3) ZC1=ZV1-ZC
  645. SCAL=XC1*XVECT+YC1*YVECT+ZC1*ZVECT
  646. XCEN=XC+SCAL*XVECT
  647. YCEN=YC+SCAL*YVECT
  648. ZCEN=ZC+SCAL*ZVECT
  649. C VECTEUR V1 ET SA NORME (= LE RAYON)
  650. XV1=XV1-XCEN
  651. YV1=YV1-YCEN
  652. ZV1=ZV1-ZCEN
  653. RAY=SQRT(XC1**2+YC1**2+ZC1**2)
  654. c longueur de la ligne
  655. DLONG=ANG*RAY
  656. IF (ABS(DLONG).LE.XZPREC*ABS(DEN1)) THEN
  657. CALL ERREUR(21)
  658. RETURN
  659. ENDIF
  660. DEN1=DEN1/DLONG
  661. DEN2=DEN2/DLONG
  662. c V3 = vecteur de base orthogonal a V1 et dans le plan (V1,VVECT)
  663. XV3=(YVECT*ZV1-YV1*ZVECT)
  664. YV3=(XV1*ZVECT-XVECT*ZV1)
  665. ZV3=(XVECT*YV1-XV1*YVECT)
  666. GOTO 200
  667.  
  668.  
  669. C ----------------------------------------------------------------------
  670. C DECOUPAGE DE LA LIGNE (BASE SUR L'ABSCISSE CURVILIGNE)
  671. C ----------------------------------------------------------------------
  672.  
  673. 200 CONTINUE
  674. CALL DECOUP(INBR,DEN1,DEN2,APROG,NBELEM,DENI,DECA,DLONG)
  675. IF (IERR.NE.0) RETURN
  676.  
  677. IF (IIMPI.EQ.1) WRITE(IOIMP,1000) NBELEM,APROG
  678. 1000 FORMAT (/,' ELEMENT ',I6,' RAISON ',G12.5)
  679. IF ((ILIG.EQ.9).AND.(NBELEM.LT.3)) CALL ERREUR(21)
  680. IF (IERR.NE.0) RETURN
  681.  
  682.  
  683. C ----------------------------------------------------------------------
  684. c CREATION DU MAILLAGE ET SEGMENT DE TRAVAIL
  685. C ----------------------------------------------------------------------
  686.  
  687. NBSOUS=0
  688. NBREF=0
  689. SEGINI MELEME
  690. ITYPEL=ITYPELT
  691. * Le nombre de noeuds à créer
  692. NBNO=NBELEM*(NBNN-1)-1
  693. c cas CERC 'ROTA' : le dernier point n'existe pas, il faut le creer
  694. c sauf si 360 degres + 'ELIM'
  695. IF(ICLE2.EQ.3.AND.IELIM.NE.2) NBNO=NBNO+1
  696. SEGINI TABPAR
  697. C Initialisation segment de travail (arc de CUBIQUE)
  698. IF ((ILIG.EQ.9).OR.(ILIG.EQ.10)) THEN
  699. NSTRAV=NITER+1
  700. SEGINI STRAV
  701. ENDIF
  702. C
  703. C Creation des elements (objet MAILLAGE MELEME)
  704. C Calcul des densites de chacun des points
  705. DIN=DEN1
  706. IPOO=nbpts
  707. DO IELEM=1,NBELEM
  708. * Point initial du segment
  709. IF (IELEM.EQ.1) THEN
  710. NUM(1,IELEM)=IPOIN1
  711. ELSE
  712. NUM(1,IELEM)=IPOO
  713. ENDIF
  714. DIN=DIN*APROG
  715. TABPAR(IELEM)=DIN
  716. IF (NBNN.EQ.3) THEN
  717. * Point milieu du segment (si quadratique)
  718. IPOO=IPOO+1
  719. NUM(2,IELEM)=IPOO
  720. ENDIF
  721. * Point final du segment
  722. IF (IELEM.EQ.NBELEM) THEN
  723. IF (ICLE2.EQ.3.AND.IELIM.NE.2) THEN
  724. IPOO=IPOO+1
  725. NUM(NBNN,IELEM)=IPOO
  726. ELSE
  727. IF (IELIM.EQ.2) THEN
  728. c 360 degres + 'ELIM' -> le dernier point = 1er point
  729. NUM(NBNN,IELEM)=IPOIN1
  730. ELSE
  731. c le dernier point est fourni en entree : on l'utilise
  732. NUM(NBNN,IELEM)=IPOIN2
  733. ENDIF
  734. ENDIF
  735. ELSE
  736. IPOO=IPOO+1
  737. NUM(NBNN,IELEM)=IPOO
  738. ENDIF
  739. ENDDO
  740.  
  741. C ----------------------------------------------------------------------
  742. C CREATION DES POINTS
  743. C ----------------------------------------------------------------------
  744.  
  745. C :::::: AJOUT PIFFARD POUR CUBIQUES ::::::
  746. IF ((ILIG.EQ.9).OR.(ILIG.EQ.10)) THEN
  747. XPOIN(1)=XINI
  748. YPOIN(1)=YINI
  749. ZPOIN(1)=ZINI
  750. DIST(1) =XZERO
  751. ENDIF
  752.  
  753. DPAR=0
  754. IADR=NBPTS
  755. IF (NBNO.GT.0) THEN
  756. NBPTS=IADR+NBNO
  757. SEGADJ MCOORD
  758. ENDIF
  759.  
  760. IF ( (ILIG .eq. 7) .and. (NBELEM .gt. 0) ) THEN
  761. C Cas de la syntaxe GIBIANE "parabole1 = PARA N1 P1 P2 P3"
  762. C avec N1 un entier positif
  763. C avec P1, P2 et P3 des points
  764.  
  765. C on cherche a calculer la longueur curviligne de la parabole 'parabole1' que l'on vas creer :
  766. C x = XIN+PARA*XDIS+(1-PARA*PARA)*XV
  767. C y = YIN+PARA*YDIS+(1-PARA*PARA)*YV
  768. C z = ZIN+PARA*ZDIS+(1-PARA*PARA)*ZV
  769. C dl = sqrt(dx^2 + dy^2 + dz^2)
  770. C => dl = sqrt(para_a*PARA^2 + para_b*PARA + para_c) : simplification
  771. para_a = (4.*XV*XV) + (4.*YV*YV)
  772. para_b = (-4.*XV*XDIS) + (-4.*YV*YDIS)
  773. para_c = (XDIS*XDIS) + (YDIS*YDIS)
  774. IF(IDIM .ge. 3) THEN
  775. para_a = para_a + (4.*ZV*ZV)
  776. para_b = para_b + (-4.*ZV*ZDIS)
  777. para_c = para_c + (ZDIS*ZDIS)
  778. ENDIF
  779. C => dl = sqrt(Beta * (x^2 + 1)) : changement de variable
  780. Beta = -1. * ((para_b**2) - (4.*para_a*para_c)) / (4.*para_a)
  781. Alpha = para_b / (2.*para_a)
  782. C DLONG = integrale de dl entre x0=x(PARA=-1) et x1=x(PARA=1)
  783. P0 = -1.
  784. P1 = 1.
  785. x0 = sqrt(abs(para_a / Beta)) * (P0 + Alpha)
  786. x1 = sqrt(abs(para_a / Beta)) * (P1 + Alpha)
  787.  
  788. dpdx_2 = Beta / sqrt(abs(para_a))
  789.  
  790. dl_x0 = (0.5 * x0 * sqrt((x0*x0) + 1.)) +
  791. & (0.5 * LOG(x0 + sqrt((x0*x0) + 1.)))
  792. dl_x1 = (0.5 * x1 * sqrt((x1*x1) + 1.)) +
  793. & (0.5 * LOG(x1 + sqrt((x1*x1) + 1.)))
  794. DLONG = dpdx_2 * (dl_x1 - dl_x0)
  795. C DENI = longueur curviligne de chaque "morceau" de la parabole
  796. DENI=DLONG/NBELEM
  797. IF(ITYPEL.NE.2) DENI=DENI/2
  798. ER_TOL = DLONG * XZPREC * 1.E4
  799. PARA = P0
  800. ENDIF
  801.  
  802. C---------------------> BOUCLE SUR LES ELEMENTS A CREER
  803. DO 301 IELEM=1,NBELEM
  804. DIN=TABPAR(IELEM)
  805. IF (NBNN.EQ.3) THEN
  806. IMIL=1
  807. NPO=2
  808. ELSE
  809. IMIL=0
  810. NPO=1
  811. ENDIF
  812. IF (IELEM.EQ.NBELEM) THEN
  813. IF (.NOT.(ICLE2.EQ.3.AND.IELIM.NE.2)) NPO=NPO-1
  814. ENDIF
  815. C---------------------> BOUCLE SUR LES NOEUDS A CREER
  816. C Si IMIL=1 point milieu de l'element a creer
  817. C Si IMIL=0 point final de l'element a creer
  818. DO 3010 IPO=1,NPO
  819. IF (NBNN.EQ.3) THEN
  820. DPAR=DPAR+DIN*XDEMI
  821. ELSE
  822. DPAR=DPAR+DIN
  823. ENDIF
  824. GOTO (1001,1001,1003,1003,1001,1001,1007,1008,1009,1009,
  825. . 1010,1003),ILIG
  826. C - Segment de DROITE
  827. 1001 CONTINUE
  828. XCOOR(IADR*(IDIM+1)+1)=XIN+DPAR*XDIS
  829. IF (IDIM.GE.2) THEN
  830. XCOOR(IADR*(IDIM+1)+2)=YIN+DPAR*YDIS
  831. IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZIN+DPAR*ZDIS
  832. ENDIF
  833. GOTO 1020
  834. C - Arc de CERCLE
  835. 1003 CONTINUE
  836. PARA=ANG*DPAR
  837. COSA=COS(PARA)
  838. SINA=SIN(PARA)
  839. IF (ILIG.EQ.12.AND.IMIL.EQ.1) THEN
  840. PAR2=ANG*DIN/XDEUX
  841. XDD=XDEUX-COS(PAR2)-(TAN(PAR2)/RAY/(XDEUX*RAY/SIN(PAR2)
  842. $ -XUN))
  843. ELSE
  844. XDD=XUN
  845. ENDIF
  846. XCOOR(IADR*(IDIM+1)+1)=XCEN+(COSA*XV1+SINA*XV3)*XDD
  847. XCOOR(IADR*(IDIM+1)+2)=YCEN+(COSA*YV1+SINA*YV3)*XDD
  848. IF (IDIM.GE.3) THEN
  849. XCOOR(IADR*(IDIM+1)+3)=ZCEN+(COSA*ZV1+SINA*ZV3)*XDD
  850. ENDIF
  851. GOTO 1020
  852.  
  853. C - Arc de PARABOLE
  854. 1007 continue
  855. IF (INBR .gt. 0) THEN
  856. dl_x1 = ( DENI / dpdx_2 ) + dl_x0
  857. C on cherche x1 tels que "dl_x1 = 0.5 * ( x1*sqrt(x1^2 + 1.) + asinh(x1) )"
  858. C => racine de "x1^4 + x1^2 - asinh(x1)^2 + 4*asinh(x1)*dl_x1 - 4*dl_x1^2"
  859. C => racine de "x1^4 + 4*dl_x1*x1 - 4*dl_x1^2", avec approximation asinh(x1)=x1
  860. NRoot = 0
  861. DO iq=1, 4
  862. Q(iq) = 0
  863. ENDDO
  864. CALL QUARTI(1.d0, 0.d0, 0.d0, 4.d0*dl_x1,-4.d0*dl_x1
  865. $ *dl_x1,Q(1), Q(2), Q(3), Q(4), NRoot)
  866. P0=PARA
  867. do iq=1, 4
  868. x1_ = Q(iq)
  869. P1_ = (x1_ / sqrt(abs(para_a/Beta))) - Alpha
  870. if (iq .eq. 1) then
  871. x1 = x1_
  872. P1 = P1_
  873. endif
  874. if((abs(P1_ - P0).lt.abs(P1 - P0)).and.(P1_.gt.P0))
  875. $ then
  876. x1 = x1_
  877. P1 = P1_
  878. endif
  879. enddo
  880. C boucle de convergence vers la solution "exacte"
  881. i_x1 = 0
  882. 107 continue
  883. C x1' = ((erreur_souhaite - erreur(x1)) / derreur_dx(x1)) + x1
  884. e_x1 = (0.5 * ((x1*sqrt((x1*x1) + 1)) +
  885. & LOG(x1 + sqrt((x1*x1) + 1.)))) - dl_x1
  886. x1 = ((0. - e_x1) / sqrt((x1*x1) + 1)) + x1
  887. i_x1=i_x1+1
  888. if(abs(e_x1) .gt. ER_TOL .and. i_x1 .lt. 10) goto 107
  889. C fin de la boucle de convergence vers la solution "exacte"
  890. dl_x0 = dl_x1
  891. PARA = (x1 / sqrt(abs(para_a/Beta))) - Alpha
  892. ELSE
  893. PARA=2*(DPAR-XDEMI)
  894. ENDIF
  895. XCOOR(IADR*(IDIM+1)+1)=XIN+PARA*XDIS+(1-PARA*PARA)*XV
  896. XCOOR(IADR*(IDIM+1)+2)=YIN+PARA*YDIS+(1-PARA*PARA)*YV
  897. IF(IDIM .ge. 3) THEN
  898. XCOOR(IADR*(IDIM+1)+3)=ZIN+PARA*ZDIS+(1-PARA*PARA)*ZV
  899. ENDIF
  900. GOTO 1020
  901. C - NON disponible
  902. 1008 CONTINUE
  903. GOTO 1020
  904. C - Arc de CUBIQUE (CUBT et CUBP)
  905. 1009 CONTINUE
  906. IF (IMIL.EQ.1) THEN
  907. CALL ERREUR(16)
  908. RETURN
  909. ELSE
  910. UP=IELEM/NBELEM
  911. U(1,i)=UP
  912. F1=2*UP**3-3*UP**2+1
  913. F2=-2*UP**3+3*UP**2
  914. F3=UP**3-2*UP**2+UP
  915. F4=UP**3-UP**2
  916. IF (ILIG.EQ.9) THEN
  917. C ...... CUBIQUES PASSANT PAR 4 POINTS .......
  918. DO J=1,4
  919. TU(J)=+CU(J,1)*XINI+CU(J,2)*XINT1+CU(J,3)*XINT2
  920. $ +CU(J,4)*XFIN
  921. ENDDO
  922. XPOIN(IELEM+1)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4)
  923. DO J=1,4
  924. TU(J)=+CU(J,1)*YINI+CU(J,2)*YINT1+CU(J,3)*YINT2
  925. $ +CU(J,4)*YFIN
  926. ENDDO
  927. YPOIN(IELEM+1)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4)
  928. IF (IDIM.GE.3) THEN
  929. DO J=1,4
  930. TU(J)=+CU(J,1)*ZINI+CU(J,2)*ZINT1+CU(J,3)*ZINT2
  931. $ +CU(J,4)*ZFIN
  932. ENDDO
  933. ZPOIN(IELEM+1)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4)
  934. ENDIF
  935. ELSE
  936. C ...... CUBIQUES PASSANT PAR 2 POINTS + TGTES .......
  937. XPOIN(IELEM+1)=F1*XINI+F2*XFIN+F3*XINT1+F4*XINT2
  938. YPOIN(IELEM+1)=F1*YINI+F2*YFIN+F3*YINT1+F4*YINT2
  939. IF (IDIM.GE.2) ZPOIN(IELEM+1)=F1*ZINI+F2*ZFIN+F3*ZINT1
  940. $ +F4*ZINT2
  941. ENDIF
  942. C ... SI IUNIF = 0 PAS DE CORRECTION SUR LES PARAM DE LA CUBIQUE ....
  943. IF (IUNIF.EQ.1) THEN
  944. DIST(IELEM+1)=DIST(IELEM)+SQRT((XPOIN(IELEM+1)
  945. $ -XPOIN(IELEM))**2+(YPOIN(IELEM+1)-YPOIN(IELEM))
  946. $ **2+(ZPOIN(IELEM+1)-ZPOIN(IELEM))**2)
  947. IF (IIMPI.EQ.1) WRITE(IOIMP,6012) DIST(IELEM+1)
  948. 6012 FORMAT(2X,'DISTANCE CUMULEE ENTRE PTS =',G13.5)
  949. DLONGT(1)=DIST(IELEM+1)
  950. GOTO 301
  951. ELSE
  952. XCOOR(IADR*(IDIM+1)+1)=XPOIN(IELEM+1)
  953. XCOOR(IADR*(IDIM+1)+2)=YPOIN(IELEM+1)
  954. IF(IDIM.GT.2) XCOOR(IADR*(IDIM+1)+3)=ZPOIN(IELEM+1)
  955. ENDIF
  956. GOTO 1020
  957. ENDIF
  958. C - Erreur anormale car ILG=11 --> ILIG=3
  959. 1010 CONTINUE
  960. CALL ERREUR(5)
  961. RETURN
  962. C - Densite du point cree
  963. 1020 CONTINUE
  964. XCOOR((IADR+1)*(IDIM+1))=DENI+DECA*DPAR
  965. IADR=IADR+1
  966. IMIL=1-IMIL
  967. 3010 CONTINUE
  968. 301 CONTINUE
  969. C Cas particulier de l'arc de CUBIQUE (CUBP) avec option UNIF
  970. C :::::::::::::::: AJOUT PIFFARD POUR CUBIQUE PAR 4 POINTS ::::::::::::
  971. C .... SI CUBIQUE ET SI IUNIF = 1 CORRECTION POUR PTS UNIFORMES.
  972. IF ((ILIG.EQ.9).AND.(IUNIF.EQ.1)) THEN
  973. C .... ON RECALCULE LES POINTS AVEC LA CORRECTION SUR LES PARAM U .
  974. C .. ON FAIT 5 ITERATIONS MAXI (ON PEUT SORTIR AVT PAR DMOY)
  975. DLONGT(1)=DLONGT(1)+SQRT((XFIN-XPOIN(NBNO+1))**2
  976. & +(YFIN-YPOIN(NBNO+1))**2
  977. & +(ZFIN-ZPOIN(NBNO+1))**2)
  978. DMOY(1)=DLONGT(1)/NBELEM
  979. IF (IIMPI.EQ.1) WRITE(IOIMP,6003) DMOY(1)
  980. 6003 FORMAT(2X,'DISTANCE MOYENNE =',G13.5)
  981. DO KT=1,NITER
  982. IF (IIMPI.EQ.1) WRITE(IOIMP,6000) KT,NITER
  983. 6000 FORMAT(2X,'ITERATION',I2,' SUR ',I2,' DEMANDEES')
  984. DO I=2,NBNO+1
  985. U(KT+1,I)=U(KT,I)*(I-XUN)/(NBNO+1)*DLONGT(KT)/DIST(I)
  986. UP=U(KT+1,I)
  987. F1=2*UP**3-3*UP**2+1
  988. F2=-2*UP**3+3*UP**2
  989. F3=UP**3-2*UP**2+UP
  990. F4=UP**3-UP**2
  991. DO J=1,4
  992. TU(J)=+CU(J,1)*XINI+CU(J,2)*XINT1
  993. & +CU(J,3)*XINT2+CU(J,4)*XFIN
  994. ENDDO
  995. XPOIN(I)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4)
  996. DO J=1,4
  997. TU(J)=+CU(J,1)*YINI+CU(J,2)*YINT1
  998. & +CU(J,3)*YINT2+CU(J,4)*YFIN
  999. ENDDO
  1000. YPOIN(I)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4)
  1001. ZPOIN(I)=XZERO
  1002. IF (IDIM.GE.3) THEN
  1003. DO J=1,4
  1004. TU(J)=+CU(J,1)*ZINI+CU(J,2)*ZINT1
  1005. & +CU(J,3)*ZINT2+CU(J,4)*ZFIN
  1006. ENDDO
  1007. ZPOIN(I)=F1*TU(1)+F2*TU(2)+F3*TU(3)+F4*TU(4)
  1008. ENDIF
  1009. DIST2(I)=SQRT((XPOIN(I)-XPOIN(I-1))**2
  1010. & +(YPOIN(I)-YPOIN(I-1))**2+(ZPOIN(I)-ZPOIN(I-1))**2)
  1011. IF (IIMPI.EQ.1) WRITE(IOIMP,6001) DIST2(I)
  1012. 6001 FORMAT(2X,'DISTANCE ENTRE 2 PTS =',G13.5)
  1013. DIST(I)=DIST(I-1)+DIST2(I)
  1014. DLONGT(KT+1)=DLONGT(KT+1)+DIST2(I)
  1015. ENDDO
  1016. DLONGT(KT+1)=DLONGT(KT+1)+SQRT((XFIN-XPOIN(NBNO+1))**2
  1017. & +(YFIN-YPOIN(NBNO+1))**2+(ZFIN-ZPOIN(NBNO+1))**2)
  1018. IF (IIMPI.EQ.1) WRITE(IOIMP,6002) DLONGT(KT+1)
  1019. 6002 FORMAT(2X,'DISTANCE TOTALE ENTRE PTS =',G13.5)
  1020. DMOY(KT)=DLONGT(KT+1)/(NBNO+1)
  1021. IF (IIMPI.EQ.1) WRITE(IOIMP,6003) DMOY(KT)
  1022. C6003 FORMAT(2X,'DISTANCE MOYENNE =',G13.5)
  1023. CRIT=0.05D0*DMOY(KT)
  1024. ICHANG=0
  1025. DO IL=2,NBNO+1
  1026. IF (ABS(DIST2(IL)-DMOY(KT)).GT.CRIT) ICHANG=ICHANG+1
  1027. ENDDO
  1028. IF (ICHANG.EQ.0) GOTO 307
  1029. ENDDO
  1030. C .... J'AI FINI LA CORRECTION SI :
  1031. C .. - JE SUIS ARRIVE AU BOUT DES ITERATIONS
  1032. C . - J'AI TROUVE LA BONNE REPARTITION
  1033. C ---- STOCKAGE DES POINTS DS LA BASE GIBI -----
  1034. 307 DPAR=XZERO
  1035. DO i=2,NBNO+1
  1036. DPAR=DPAR+TABPAR(i)
  1037. XCOOR(IADR*(IDIM+1)+1)=XPOIN(i)
  1038. XCOOR(IADR*(IDIM+1)+2)=YPOIN(i)
  1039. IF (IDIM.GE.3) XCOOR(IADR*(IDIM+1)+3)=ZPOIN(i)
  1040. XCOOR((IADR+1)*(IDIM+1))=DENI+DECA+DPAR
  1041. IADR=IADR+1
  1042. ENDDO
  1043. ENDIF
  1044.  
  1045.  
  1046. C ----------------------------------------------------------------------
  1047. C UN PEU DE MENAGE !
  1048. C ----------------------------------------------------------------------
  1049.  
  1050. SEGSUP TABPAR
  1051. IF ((ILIG.EQ.9).OR.(ILIG.EQ.10)) SEGSUP STRAV
  1052.  
  1053. IPT2=MELEME
  1054. DO i=1,NUM(/2)
  1055. ICOLOR(i)=IDCOUL
  1056. ENDDO
  1057.  
  1058. c EVENTUELLE FUSION AVEC LIGNES INITIALES ET FINALES
  1059. IF (ITP1.NE.0) THEN
  1060. ltelq=.false.
  1061. CALL FUSE(IPT1,MELEME,IPT2,ltelq)
  1062. SEGDES IPT1,MELEME
  1063. ENDIF
  1064. MELEME=IPT2
  1065. IF (ITP2.NE.0) THEN
  1066. ltelq=.false.
  1067. CALL FUSE(MELEME,IP2,IPT2,ltelq)
  1068. SEGDES IPT6,MELEME
  1069. ENDIF
  1070.  
  1071. c ECRITURE
  1072. MELEME=IPT2
  1073. SEGDES MELEME
  1074. IF (IERR.NE.0) SEGSUP MELEME
  1075. CALL ACTOBJ('MAILLAGE',MELEME,1)
  1076. CALL ECROBJ('MAILLAGE',MELEME)
  1077.  
  1078. RETURN
  1079. END
  1080.  
  1081.  
  1082.  
  1083.  
  1084.  

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