C ISOVA2    SOURCE    BP208322  11/12/20    21:15:02     7231           
      SUBROUTINE ISOVA2(MELEME,MCHAML,XISO,XTOL,ELEMS)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : ISOVA2
C DESCRIPTION : Construit le maillage d'une isovaleur d'un champ par
C               éléments
C     
C     
C     
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C VERSION    : v1, 17/12/2008, version initiale
C HISTORIQUE : v1, 17/12/2008, création
C HISTORIQUE : 19/12/2011, bp : ajout cub8
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
-INC CCOPTIO
-INC CCGEOME
-INC SMCOORD
-INC SMELEME
-INC SMCHAML
-INC SMLENTI
*     
* Segments ajustables contenant les éléments créés ainsi que leur couleur
* ELEM(1) contient des POI1
* ELEM(2) contient des SEG2
* ELEM(3) contient des TRI3
* ELEM(4) contient des TET4
* ELEM(5) contient des QUA4
*     
      PARAMETER (NTYEL=5)
      SEGMENT ELEMS
      POINTEUR ELEM(NTYEL).MLENTI
      ENDSEGMENT
*     
      PARAMETER (NMAXNO=8)
      REAL*8 VAL(NMAXNO),VALP(NMAXNO)
      LOGICAL LVAL(NMAXNO),LTRI
      INTEGER IPERM(NMAXNO),NUMP(NMAXNO)
      CHARACTER*8 MCHA
      INTEGER IMPR,IRET
      
*     tableaux descriptifs des cas et des aretes 
      PARAMETER (NMAXPT=16,NMAXCA=256,npt2=2,NMAXAR=12)
      INTEGER MATCAS(NMAXPT,NMAXCA),MATARE(npt2,NMAXAR)

*     
* Fonctions locales
*     
      LOGICAL LEGISO
      LEGISO(X)=(ABS(X-XISO).LT.XTOL)
*     
* Executable statements
*     
      IMPR=0
      IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans isova2.eso' 
*     
C***********************************************************************
c     Boucle sur les zones elementaires      
C***********************************************************************
c
      NNO=NUM(/1)
      NEL=NUM(/2)
      N2=IELVAL(/1)
      DO I2=1,N2
      
         MELVAL=IELVAL(I2)
         SEGACT MELVAL
         
* Vérification de la cohérence entre la champ et la géométrie
         N1PTEL=VELCHE(/1)
         N1EL=VELCHE(/2)
         IF ((N1PTEL.NE.1.AND.N1PTEL.NE.NNO).OR.
     $        (N1EL.NE.1.AND.N1EL.NE.NEL)) THEN
* 148 2
*Le champ par élément et l'objet géométrique correspondant n'ont pas le
* 148 2
*même nombre de points par élément : contacter votre support
            CALL ERREUR(148)
            RETURN
         ENDIF
C
c******* Aiguillage selon type d element *******************************
c
         IF(ITYPEL.eq.1)  GOTO 101
         IF(ITYPEL.eq.2)  GOTO 102
         IF(ITYPEL.eq.4)  GOTO 104 
         if(ITYPEL.eq.14)  goto 114
         if(ITYPEL.eq.23)  goto 123
         
******** Cas non traité *********************************************
         MOTERR(1:4)=NOMS(ITYPEL)
*  44 2
*Type d'élément inconnu %m1:4
         CALL ERREUR(44)
         RETURN
         
******** Cas des POI1 ***********************************************
 101     CONTINUE
         DO IEL=1,NEL
               VAL(1)=VELCHE(1,MIN(N1EL,IEL))
               IF (LEGISO(VAL(1))) THEN
                  MLENTI=ELEM(1)
                  LECT(**)=NUM(1,IEL)
                  LECT(**)=ICOLOR(IEL)
               ENDIF
         ENDDO
         goto 900
               
******** Cas des SEG2 ***********************************************
 102     CONTINUE
         DO IEL=1,NEL
               VAL(1)=VELCHE(1,MIN(N1EL,IEL))
               VAL(2)=VELCHE(MIN(N1PTEL,2),MIN(N1EL,IEL))
               VMIN=MIN(VAL(1),VAL(2))
               VMAX=MAX(VAL(1),VAL(2))
*     Y a-t-il une isovaleur ?
               IF (((VMIN-XISO).LE.XTOL).AND.((VMAX-XISO).GE.-XTOL))
     $              THEN
                  NPEQ1=0
                  DO IVAL=1,2
                     IF (LEGISO(VAL(IVAL))) THEN
                        NPEQ1=NPEQ1+1
                     ENDIF
                  ENDDO
                  IF (NPEQ1.EQ.2) THEN
* Le segment entier
                     MLENTI=ELEM(2)
                     LECT(**)=NUM(1,IEL)
                     LECT(**)=NUM(2,IEL)
                     LECT(**)=ICOLOR(IEL)
                  ELSEIF (NPEQ1.EQ.1) THEN
* Une des extrémités
                     MLENTI=ELEM(1)
                     IF (LEGISO(VAL(1))) THEN
* Le premier point
                        LECT(**)=NUM(1,IEL)
                     ELSE
* Le deuxième point
                        LECT(**)=NUM(2,IEL)
                     ENDIF
                     LECT(**)=ICOLOR(IEL)
                  ELSE
* Un point milieu
                     MLENTI=ELEM(1)
                     VAL1=VAL(1)
                     VAL2=VAL(2)
                     NUM1=NUM(1,IEL)
                     NUM2=NUM(2,IEL)
                     CALL ISOVA3(XISO,VAL1,VAL2,NUM1,NUM2,MLENTI)
                     LECT(**)=ICOLOR(IEL)
                  ENDIF
               ENDIF
         ENDDO
         goto 900
               
******** Cas des TRI3 ***********************************************
 104     CONTINUE
         DO IEL=1,NEL
               VAL(1)=VELCHE(1,MIN(N1EL,IEL))
               VAL(2)=VELCHE(MIN(N1PTEL,2),MIN(N1EL,IEL))
               VAL(3)=VELCHE(MIN(N1PTEL,3),MIN(N1EL,IEL))
               VMIN=MIN(VAL(1),VAL(2),VAL(3))
               VMAX=MAX(VAL(1),VAL(2),VAL(3))
*     Y a-t-il une isovaleur ?
               IF (((VMIN-XISO).LE.XTOL).AND.((VMAX-XISO).GE.-XTOL))
     $              THEN
* Permutation permettant d'avoir VAL en ordre croissant
* VAL(IPERM(1)).LE.VAL(IPERM(2)).LE.VAL(IPERM(3))
                  DO IVAL=1,3
                     IPERM(IVAL)=IVAL
                  ENDDO
                  CALL ISHELR(3,IPERM,3,VAL,IMPR,IRET)
                  NPEQ1=0
                  DO IVAL=1,3
                     VALP(IVAL)=VAL(IPERM(IVAL))
                     NUMP(IVAL)=NUM(IPERM(IVAL),IEL)
                     LVAL(IVAL)=(LEGISO(VALP(IVAL)))
                     IF (LVAL(IVAL)) THEN
                        NPEQ1=NPEQ1+1
                     ENDIF
                  ENDDO
C                        WRITE(IOIMP,*) 'IEL=',IEL
C                        DO I=1,3
C                          WRITE(IOIMP,*) I,' ',NUM(I,IEL),' ',VAL(I),
C     $                          ' ',IPERM(I),' ',LVAL(I)
C                        ENDDO
                  IF (NPEQ1.EQ.3) THEN
* Le triangle entier
                     MLENTI=ELEM(3)
                     LECT(**)=NUMP(1)
                     LECT(**)=NUMP(2)
                     LECT(**)=NUMP(3)
                     LECT(**)=ICOLOR(IEL)
                  ELSEIF (NPEQ1.EQ.2) THEN
* Un des côtés
                     MLENTI=ELEM(2)
                     IF (LVAL(1)) THEN
                        LECT(**)=NUMP(1)
                     ELSE
                        LECT(**)=NUMP(3)
                     ENDIF
                     LECT(**)=NUMP(2)
                     LECT(**)=ICOLOR(IEL)
                  ELSEIF (NPEQ1.EQ.1) THEN
* Un des points ou un segment
                     IF (LVAL(1)) THEN
*  Le point P1
                        MLENTI=ELEM(1)
                        LECT(**)=NUMP(1)
                        LECT(**)=ICOLOR(IEL)
                     ELSEIF (LVAL(2)) THEN
* Un segement reliant P2 à l'intersection sur (P1, P3)
                        MLENTI=ELEM(2)
                        LECT(**)=NUMP(2)
                        CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1),NUMP(3)
     $                       ,MLENTI)
                        LECT(**)=ICOLOR(IEL)
                     ELSE
*  Le point P3
                        MLENTI=ELEM(1)
                        LECT(**)=NUMP(3)
                        LECT(**)=ICOLOR(IEL)
                     ENDIF
                  ELSE
* Un segment reliant (P1,P2) à (P1,P3)
*                 ou (P2,P3) à (P1,P3)
                     MLENTI=ELEM(2)
                     CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1),NUMP(3)
     $                    ,MLENTI)
                     IF (VALP(2).GT.XISO) THEN
* Un segment reliant (P1,P2) à (P1,P3)
                        CALL ISOVA3(XISO,VALP(1),VALP(2),NUMP(1),NUMP(2)
     $                       ,MLENTI)
                        
                     ELSE
* Un segment reliant (P2,P3) à (P1,P3)
                        CALL ISOVA3(XISO,VALP(2),VALP(3),NUMP(2),NUMP(3)
     $                       ,MLENTI)
                     ENDIF
                     LECT(**)=ICOLOR(IEL)
                  ENDIF
               ENDIF
         ENDDO
         goto 900
               
******** Cas des CUB8 ***********************************************
 114     CONTINUE
*   On remplit les tableaux MATCAS et MATARE decrivant les cas et aretes 
c          write(6,*) 'ISOVA2: cub8: appel a ISOVA4'
         CALL ISOVA4(ITYPEL,MATCAS,NMAXCA,NMAXPT,MATARE,2,NMAXAR)
c          write(6,*) 'ISOVA2: MATARE=',(MATARE(iou,1),iou=1,2)
c          write(6,*) '               ',(MATARE(iou,2),iou=1,2)
c          write(6,*) '                ...'
c          write(6,*) 'ISOVA2: MATCAS=',(MATCAS(iou,1),iou=1,16)
c          write(6,*) '               ',(MATCAS(iou,2),iou=1,16)
c          write(6,*) '                ...'
c          write(6,*) '               ',(MATCAS(iou,255),iou=1,16)
c          write(6,*) '               ',(MATCAS(iou,256),iou=1,16)
         
* REM (BP) : 
* - on ne fait pas de test avec la tolerance (cf LEGISO autres itypel) et
*   on espere que VAL(IVAL).le.XISO au lieu de VAL(IVAL).lt.XISO suffira...
* - on ne construit que des triangles TRI3 en se basant sur l algo :
*   "the marching cubes algorithm" (domaine public)  
         DO 214 IEL=1,NEL
            ICAS = 1
            DO IVAL=1,8
               VAL(IVAL)=VELCHE(MIN(N1PTEL,IVAL),MIN(N1EL,IEL))
c                if(VAL(IVAL).le.XISO) ICAS=ICAS+(2**(IVAL-1))
            ENDDO
c           on teste les noeuds (attention: numerotation cub8 cast3m 
c           differente de celle du marching cubes algorithm)            
            if(VAL(1).le.XISO) ICAS=ICAS+1
            if(VAL(4).le.XISO) ICAS=ICAS+2
            if(VAL(3).le.XISO) ICAS=ICAS+4
            if(VAL(2).le.XISO) ICAS=ICAS+8
            if(VAL(5).le.XISO) ICAS=ICAS+16
            if(VAL(8).le.XISO) ICAS=ICAS+32
            if(VAL(7).le.XISO) ICAS=ICAS+64
            if(VAL(6).le.XISO) ICAS=ICAS+128
c             write(6,*) 'ICAS=',ICAS
c             write(6,*) 'MATCAS=',(MATCAS(iou,ICAS),iou=1,16)
*           Y a-t-il une isovaleur ?
            IF (ICAS.gt.1.or.ICAS.lt.256) THEN
                MLENTI=ELEM(3)
                do iar1=1,15
                  NUMAR1=MATCAS(iar1,ICAS)
c                   write(6,*) 'iar1,NUMAR1=',iar1,NUMAR1
                  if(NUMAR1.eq.0) goto 214
                  ipoin1 = MATARE(1,NUMAR1)
                  ipoin2 = MATARE(2,NUMAR1)
                  VAL1=VAL(ipoin1)
                  VAL2=VAL(ipoin2)
                  NUM1=NUM(ipoin1,IEL)
                  NUM2=NUM(ipoin2,IEL)
c                   write(6,*) '  ipoin1 et 2 =',ipoin1,ipoin2,NUM1,NUM2
                  CALL ISOVA3(XISO,VAL1,VAL2,NUM1,NUM2,MLENTI)
c                   write(6,*) '  VAL1 et 2=',VAL1,VAL2
                  if(mod(iar1,3).eq.0) LECT(**)=ICOLOR(IEL)
                enddo
            ENDIF
 214     CONTINUE
c          write(6,*) '=> LECT=',(LECT(iou),iou=1,LECT(/1))
         goto 900

******** Cas des TET4 ***********************************************
 123     CONTINUE
         DO IEL=1,NEL
               DO IVAL=1,4
                  VAL(IVAL)=VELCHE(MIN(N1PTEL,IVAL),MIN(N1EL,IEL))
               ENDDO
               VMIN=MIN(VAL(1),VAL(2),VAL(3),VAL(4))
               VMAX=MAX(VAL(1),VAL(2),VAL(3),VAL(4))
*     Y a-t-il une isovaleur ?
               IF (((VMIN-XISO).LE.XTOL).AND.((VMAX-XISO).GE.-XTOL))
     $              THEN
* Permutation permettant d'avoir VAL en ordre croissant
* VAL(IPERM(1)).LE.VAL(IPERM(2)).LE.VAL(IPERM(3)).LE.VAL(IPERM(4))
                  DO IVAL=1,4
                     IPERM(IVAL)=IVAL
                  ENDDO
                  CALL ISHELR(4,IPERM,4,VAL,IMPR,IRET)
                  NPEQ1=0
                  DO IVAL=1,4
                     VALP(IVAL)=VAL(IPERM(IVAL))
                     NUMP(IVAL)=NUM(IPERM(IVAL),IEL)
                     LVAL(IVAL)=(LEGISO(VALP(IVAL)))
                     IF (LVAL(IVAL)) THEN
                        NPEQ1=NPEQ1+1
                     ENDIF
                  ENDDO
C                        WRITE(IOIMP,*) 'IEL=',IEL
C                        DO I=1,3
C                          WRITE(IOIMP,*) I,' ',NUM(I,IEL),' ',VAL(I),
C     $                          ' ',IPERM(I),' ',LVAL(I)
C                        ENDDO
                  IF (NPEQ1.EQ.4) THEN
*  Isovaleur = tetrahedre
                     MLENTI=ELEM(4)
                     DO IVAL=1,4
                        LECT(**)=NUMP(IVAL)
                     ENDDO
                     LECT(**)=ICOLOR(IEL)
                  ELSEIF (NPEQ1.EQ.3) THEN
*  Isovaleur = 1 face du tetrahedre
                     MLENTI=ELEM(3)
                     IF (LVAL(1)) THEN
                        LECT(**)=NUMP(1)
                     ELSE
                        LECT(**)=NUMP(4)
                     ENDIF
                     LECT(**)=NUMP(2)
                     LECT(**)=NUMP(3)
                     LECT(**)=ICOLOR(IEL)
                  ELSEIF (NPEQ1.EQ.2) THEN
*  Isovaleur = 1 segment ou 1 triangle
*   appuye sur 2 points du tetrahedre
                     IF (LVAL(1).OR.LVAL(4)) THEN
*  Isovaleur = 1 segment
                        MLENTI=ELEM(2)
                        DO IVAL=1,4
                           IF (LVAL(IVAL)) THEN
                              LECT(**)=NUMP(IVAL)
                           ENDIF
                        ENDDO
                        LECT(**)=ICOLOR(IEL)
                     ELSE
*  Isovaleur = 1 triangle appuye sur 2 points tetra
                        MLENTI=ELEM(3)
                        LECT(**)=NUMP(2)
                        LECT(**)=NUMP(3)
                        CALL ISOVA3(XISO,VALP(1),VALP(4),NUMP(1),NUMP(4)
     $                       ,MLENTI)
                        LECT(**)=ICOLOR(IEL)
                     ENDIF
                  ELSEIF (NPEQ1.EQ.1) THEN
*  Isovaleur = 1 point ou 1 triangle
*   appuye sur 1 point du tetrahedre
                     IF (LVAL(1).OR.LVAL(4)) THEN
*  Isovaleur = 1 point
                        MLENTI=ELEM(1)
                        DO IVAL=1,4
                           IF (LVAL(IVAL)) THEN
                              LECT(**)=NUMP(IVAL)
                           ENDIF
                        ENDDO
                        LECT(**)=ICOLOR(IEL)
                     ELSE
*  Isosurface = 1 triangle appuye sur 1 point
*   du tetrahedre
                        MLENTI=ELEM(3)
                        IF (LVAL(2)) THEN
*  Le point appuye est P2
*  les autres points intersection (P1,P3)
                           LECT(**)=NUMP(2)
                           CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1)
     $                          ,NUMP(3),MLENTI)
                        ELSE
*  Le point appuye est P3
*  les autres points intersection (P2,P4)
                           LECT(**)=NUMP(3)
                           CALL ISOVA3(XISO,VALP(2),VALP(4),NUMP(2)
     $                          ,NUMP(4),MLENTI)
                        ENDIF
*  Il y a toujours un point intersection (P1,P4)
                        CALL ISOVA3(XISO,VALP(1),VALP(4),NUMP(1)
     $                       ,NUMP(4),MLENTI)
                        LECT(**)=ICOLOR(IEL)
                     ENDIF
                  ELSE
*  Isosurface = 1 triangle quelconque section
*   du tetrahedre ou un quadrangle (= 2 triangles)
                     IF (VALP(2).GT.XISO) THEN
*  Isosurface entre V1 et V2
*  Points intersection (P1,P2) (P1,P3)
*  Un seul triangle
                        MLENTI=ELEM(3)
                        CALL ISOVA3(XISO,VALP(1),VALP(2),NUMP(1)
     $                       ,NUMP(2),MLENTI)
                        CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1)
     $                       ,NUMP(3),MLENTI)
                     ELSEIF (VALP(3).LT.XISO) THEN
*  Isosurface entre V3 et V4
*  Points intersection (P3,P4) (P2,P4)
*  Un seul triangle
                        MLENTI=ELEM(3)
                        CALL ISOVA3(XISO,VALP(3),VALP(4),NUMP(3)
     $                       ,NUMP(4),MLENTI)
                        CALL ISOVA3(XISO,VALP(2),VALP(4),NUMP(2)
     $                       ,NUMP(4),MLENTI)
                     ELSE
*  Isosurface entre V2 et V3 ;
*  Points intersection (P1,P3) (P2,P3) (P2,P4)
*  Un quadrangle = 2 triangles
                        MLENTI=ELEM(5)
                        CALL ISOVA3(XISO,VALP(1),VALP(3),NUMP(1)
     $                       ,NUMP(3),MLENTI)
                        CALL ISOVA3(XISO,VALP(2),VALP(3),NUMP(2)
     $                       ,NUMP(3),MLENTI)
                        CALL ISOVA3(XISO,VALP(2),VALP(4),NUMP(2)
     $                       ,NUMP(4),MLENTI)
                     ENDIF
*  Il y toujours un point intersection
*  entre P1 et P4
                     CALL ISOVA3(XISO,VALP(1),VALP(4),NUMP(1)
     $                    ,NUMP(4),MLENTI)
                     LECT(**)=ICOLOR(IEL)
                  ENDIF
               ENDIF
         ENDDO
         goto 900
         
               
******** on termine proprement **************************************
 900     CONTINUE
         SEGDES MELVAL
         
      ENDDO
c     fin de Boucle sur les zones elementaires      
C***********************************************************************
      
      RETURN
*     
* End of subroutine ISOVA2
*     
      END