Télécharger ligne.eso

Retour à la liste

Numérotation des lignes :

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

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