Sources Esope | Cast3M

courb5
C COURB5    SOURCE    PV        20/03/24    21:16:23     10554          
      SUBROUTINE COURB5 (MCOURB,MABSCI)
************************************************************************
*
*                             C O U R B 5
*                             -----------
*
* FONCTION:
* ---------
*
*     CALCUL DES COORDONNEES DES NOUVEAUX POINTS CREES POUR LA COURBE
*     POLYNOMIALE, COURBE POUR LAQUELLE IL N'AVAIT PAS ETE DEMANDE DE
*     REGULARITE SELON L'ABSCISSE CURVILIGNE.
*
* MODULES UTILISES:
* -----------------
*
      IMPLICIT INTEGER(I-N)
      IMPLICIT real*8 (a-h,o-z)
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SMCOORD
-INC TMCOURB
*
* PARAMETRES:   (E)=ENTREE   (S)=SORTIE   (+ = CONTENU DANS UN COMMUN)
* -----------
*
*     MCOURB  (E)  SEGMENT ACTIF.
*             (S)  ETENDU EN SORTIE (TABLEAU "UCOU" REMPLI).
*     MABSCI  (E)  SEGMENT ACTIF, DEFINI DANS LE S-P "COURB2".
*                  ICI, IL S'AGIT D'ABSCISSES PARAMETRIQUES.
*    +IDIM    (E)  VOIR LE COMMUN "COPTIO".
*    +ILCOUR  (E)  VOIR LE COMMUN "CGEOME".
*    +MCOORD  (E)  SEGMENT ACTIF.
*             (S)  ETENDU EN SORTIE.
*
      SEGMENT,MABSCI
         REAL*8 ABSCIS(NPOIN)
      ENDSEGMENT
*
* VARIABLES:
* ----------
*
*     D1     = DISTANCE DU POINT COURANT AVEC LE POINT PRECEDENT.
*     D3     = DISTANCE DU POINT COURANT AVEC LE POINT SUIVANT.
*     NPOLD  = NOMBRE DE POINTS AVANT AJOUT.
*     NPOIN  = NOMBRE DE POINTS A AJOUTER.
*     U2     = PARAMETRE DU POINT COURANT.
*     UINF, USUP = BORNES ENCADRANT LE PARAMETRE DU POINT-MILIEU
*                  COURANT.
*
*     EN GENERAL, INDICE 1 POUR POINT PRECEDENT, INDICE 2 POUR POINT
*     COURANT ET INDICE 3 POUR POINT SUIVANT.
*
      INTEGER NPOLD,NPOIN
      REAL*8 USUP,UINF,D1,D3
      REAL*8 U2
*
* FONCTIONS:
* ----------
*
      REAL*8 POLYNO
*
* REMARQUES:
* ----------
*
*     CE SOUS-PROGRAMME N'EST PAS PREVU POUR FONCTIONNER AVEC UN NOMBRE
*     DE POINTS NUL.
*
* AUTEUR, DATE DE CREATION:
* -------------------------
*
*     PASCAL MANIGOT     10 SEPTEMBRE 1986
*     P.M.   24/02/87 : REMODELAGE, AVEC UN NOUVEAU SEGMENT "MCOURB".
*
* LANGAGE:
* --------
*
*     ESOPE77 FORTRAN77 + EXTENSION: DECLARATION "REAL*8".
*
************************************************************************
*
      NPOIN = ABSCIS(/1)
      IDIMP1 = IDIM + 1
      segact mcoord*mod
      NPOLD = nbpts 
*
*     DETERMINATION DES COORDONNEES DES NOUVEAUX POINTS.
*
      LONG = NPOIN
      SEGADJ,MCOURB
      NBPTA = NPOLD
      NBPTS=NBPTA+NPOIN
      SEGADJ MCOORD
      MCOFCO = ICOFCO
      SEGACT,MCOFCO
*
      DO 300 IB=1,NPOIN
         U2 = ABSCIS(IB)
         UCOU(IB) = U2
         XCOOR(NBPTA*IDIMP1+1) = POLYNO (COFCOU(1,1),ND1COU,1,U2)
         XCOOR(NBPTA*IDIMP1+2) = POLYNO (COFCOU(1,2),ND1COU,1,U2)
         IF (IDIM .EQ. 3)
     #   XCOOR(NBPTA*IDIMP1+3) = POLYNO (COFCOU(1,3),ND1COU,1,U2)
         NBPTA=NBPTA+1
  300    CONTINUE
*     END DO
*
      IF (ILCOUR .EQ. 3) THEN
*
*        AJUSTEMENT DE LA PLACE DES POINTS-MILIEUX POUR QU'ILS SOIENT
*        BIEN ... AU MILIEU DES ELEMENTS.
*        LES POINTS-MILIEUX SONT AUX PLACES 1, 3, ... , NPOIN  (AU
*        COEFFICIENT MULTIPLICATIF IDIM+1 PRES) DANS LA PARTIE AJOUTEE
*        DE "XCOOR".
*        ON EN PROFITE POUR INSCRIRE LA DENSITE DES POINTS-MILIEUX.
*
         I3 = (PT1COU - 1) * IDIMP1
         X3 = XCOOR(I3+1)
         Y3 = XCOOR(I3+2)
         IF (IDIM .EQ. 3) Z3 = XCOOR(I3+3)
         U3 = U1COU
*
         DO 310 IB=1,NPOIN,2
*
            U1 = U3
            X1 = X3
            Y1 = Y3
            I2 = (IB + NPOLD - 1) * IDIMP1
            U2 = ABSCIS(IB)
            X2 = XCOOR(I2+1)
            Y2 = XCOOR(I2+2)
            IF (IB .EQ. NPOIN) THEN
               I3 = (PT2COU - 1) * IDIMP1
               U3 = U2COU
            ELSE
               I3 = (IB + NPOLD) * IDIMP1
               U3 = ABSCIS(IB+1)
            END IF
            X3 = XCOOR(I3+1)
            Y3 = XCOOR(I3+2)
            IF (IDIM .EQ. 3) THEN
               Z1 = Z3
               Z2 = XCOOR(I2+3)
               Z3 = XCOOR(I3+3)
               D1 = SQRT( (X2-X1)**2 + (Y2-Y1)**2 + (Z2-Z1)**2 )
               D3 = SQRT( (X2-X3)**2 + (Y2-Y3)**2 + (Z2-Z3)**2 )
            ELSE
               D1 = SQRT( (X2-X1)**2 + (Y2-Y1)**2 )
               D3 = SQRT( (X2-X3)**2 + (Y2-Y3)**2 )
            END IF
            UINF = U1
            USUP = U3
*
            NFOIS = 0
  315       IF (ABS(D3-D1) .GT. (0.1*(D1+D3)) .AND. NFOIS.LT.5) THEN
*
               IF (D1 .GT. D3) THEN
                  USUP = U2
               ELSE
                  UINF = U2
               END IF
               U2 = (UINF + USUP) / 2.
               X2 = POLYNO (COFCOU(1,1),ND1COU,1,U2)
               Y2 = POLYNO (COFCOU(1,2),ND1COU,1,U2)
               IF (IDIM .EQ. 3) THEN
                  Z2 = POLYNO (COFCOU(1,3),ND1COU,1,U2)
                  D1 = SQRT( (X2-X1)**2 + (Y2-Y1)**2 + (Z2-Z1)**2 )
                  D3 = SQRT( (X2-X3)**2 + (Y2-Y3)**2 + (Z2-Z3)**2 )
               ELSE
                  D1 = SQRT( (X2-X1)**2 + (Y2-Y1)**2 )
                  D3 = SQRT( (X2-X3)**2 + (Y2-Y3)**2 )
               END IF
               NFOIS = NFOIS + 1
*
               GOTO 315
            END IF
*
            IF (NFOIS .GT. 0) THEN
*              RECTIFICATION DES COORDONNEES ET DU PARAMETRE:
               XCOOR(I2+1) = X2
               XCOOR(I2+2) = Y2
               IF (IDIM .EQ. 3) XCOOR(I2+3) = Z2
               UCOU(IB) = U2
            END IF
            XCOOR(I2+IDIMP1) = D1 + D3
*           CI-DESSUS, ON NE DIVISE PAS PAR 2 CAR, POUR UN "SEG3", IL NE
*           FAUT PAS CONFONDRE LARGEUR DE MAILLE ET DISTANCE ENTRE 2
*           POINTS.
            IF (IIMPI .EQ. 342) WRITE (IOIMP,*) 'NFOIS = ',NFOIS
*
  310       CONTINUE
*        END DO
*
*        LA BOUCLE 330 VA TRAITER LES POINTS NOUVEAUX N. 2, 4, 6, ...
         IB1 = 2
         IB2 = NPOIN - 1
         IB3 = 2
*
      ELSE
*        ILCOUR = 2
*
*        LA BOUCLE 330 VA TRAITER TOUS LES POINTS NOUVEAUX.
         IB1 = 1
         IB2 = NPOIN
         IB3 = 1
*
      END IF
*
      SEGDES,MCOFCO
*
*     DETERMINATION DES DENSITES DES EXTREMITES DES SEGMENTS:
*
      IF (IB1 .LE. IB2) THEN
         I2 = (PT1COU - 1) * IDIMP1
         I3 = (IB1 + NPOLD - 1) * IDIMP1
         D3 = (XCOOR(I3+1) - XCOOR(I2+1) ) ** 2
     &      + (XCOOR(I3+2) - XCOOR(I2+2) ) ** 2
         IF (IDIM .EQ. 3) D3 = D3 + (XCOOR(I3+3) - XCOOR(I2+3) )**2
         D3 = SQRT(D3)
         DO 330 IB=IB1,IB2,IB3
            I2 = I3
            IF (IB .EQ. IB2) THEN
               I3 = (PT2COU - 1) * IDIMP1
            ELSE
               I3 = (IB + NPOLD - 1 + IB3) * IDIMP1
            END IF
            D1 = D3
            D3 = (XCOOR(I3+1) - XCOOR(I2+1) ) ** 2
     &         + (XCOOR(I3+2) - XCOOR(I2+2) ) ** 2
            IF (IDIM .EQ. 3) D3 = D3 + (XCOOR(I3+3) - XCOOR(I2+3) )**2
            D3 = SQRT(D3)
            XCOOR(I2+IDIMP1) = (D1 + D3) / 2.
  330       CONTINUE
*        END DO
      END IF
*
      END