Télécharger ligne.eso

Retour à la liste

Numérotation des lignes :

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

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