Télécharger la plaquette Cast3M


fr en

Télécharger ligne.eso

Retour à la liste

Numérotation des lignes : 

ligne

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

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