Télécharger ligne.eso

Retour à la liste

Numérotation des lignes :

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

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