Procédures | Cast3M

* @ISOSURF  PROCEDUR  GOUNAND   15/09/16    21:15:00     8625           
************************************************************************
*    Procedure qui extrait les isosurfaces dont les valeurs sont
*    listées dans une liste de réels (LIS1) d'un champoint (HANA1)
*    appuyé sur un maillage (MASSIF0).
*
*    Le résultat final est constitué du maillage surfacique
*    regroupant l'ensemble des isosurface MAIF1 et du champoint
*    CHPF1 des isovaleurs LIS1 appuyées sur MAIF1.
*
*    Postraitement TRAC CACH MAIF1 CHPF1 ;
*
*  Syntaxe :
*  ---------
*            MAIL1 CHPF1 = @ISOSURF MASSIF0 LIS1 HANA1 ;
*
*    Entrée  :
*    ---------
*    MASSIF0  : Maillage support du champoint
*
*    LIS1     : Liste (LISTREEL) d'isovaleurs à rechercher
*
*    HANA1    : Champoint appuye sur MASSIF0
*
*    Sortie  :
*    ---------
*    MAIF1    : Maillage de l'ensemble des isosurfaces
*
*    CHPF1    : Champoint des isovaleurs LIS1 appuyées sur MAIF1
*
************************************************************************
* Remarque 1 : Attention la procédure utilise une élimination des points
*              doubles des isosurfaces extraites avec un epsilon
*              calculé aussi intelligemment que possible.
*              (Gounand: l'algo en Gibiane utilise ELIM mais l'algo en
*               Esope n'en a plus besoin, et fait l'élimination en
*               interne dans l'opérateur ISOV depuis la fiche 8152).
*
* Remarque 2 : pour l'algo Gibiane :
*              6 isovaleurs sur  1016 tetra : cpu =    3.7 s
*              6 isovaleurs sur  6950 tetra : cpu =   39 s
*              6 isovaleurs sur 52000 tetra : cpu = 1315 s
*   ---> A reprogrammer en fortran (Gounand : c'est fait !)
*************************************************************************
'DEBPROC' @ISOSURF MASSIF0*'MAILLAGE' LIS1*'LISTREEL' HANA1*'CHPOINT' ;
*
* Paramètres
* lgibi = FAUX : on utilise l'opérateur ISOV pour construire
*                les isosurfaces
* lgibi = VRAI : ancien algorithme en Gibiane pour mémoire
lgibi = FAUX ;
* Garde-fou liste vide
NIS1 = DIME LIS1 ;
SI ( NIS1 '<' 1 )          ;
     'ERREUR' 'Liste vide' ;
     'QUITTER' @ISOSURF    ;
FINSI ;
*
* Garde-fou isovaleur non comprise dans le champoint
MAXC1 = 'MAXIMUM' HANA1 ;
MINC1 = 'MINIMUM' HANA1 ;
MAXL1 = 'MAXIMUM' LIS1 ;
MINL1 = 'MINIMUM' LIS1 ;
vmax = 'MAXIMUM' ('PROG' maxc1 maxl1 minc1 minl1) 'ABS' ;
vref = '+' ('*' vmax 1.D-14) 1.D-100 ;
SI ( (MAXL1 '>' ('+' MAXC1 vref)) 'OU' (MINL1 '<' ('-' MINC1 vref)) )  ;
     'ERREUR' 'Isovaleur de la liste hors champs' ;
     'QUITTER' @ISOSURF    ;
FINSI ;
*
* Changement du maillage en tétraèdres
*
chp = hana1 ; mail = massif0 ;
mailt = 'CHANGER' mail 'TET4' ;
*
* S'il y a une différence de nombres de noeuds entre mail et mailt
* on va faire une projection
*
'SI' ('>' ('NBNO' mailt) ('NBNO' mail)) ;
   mailp = 'CHANGER' 'POI1' mail ;
   mailtp = 'CHANGER' 'POI1' mailt ;
   dmp = 'DIFF' mailp mailtp ;
   cha = 'CHANGER' 'CHAM' chp mail ;
   chp2 = 'PROI' dmp cha ;
   chp = chp '+' chp2 ;
'FINSI' ;
massif0 = mailt ;
hana1   = chp ;
*
* Garde-fou maillage non uniquement constitué de TET4
*gounand Inutile normalement vu les lignes précédentes
*LMOT1 = MASSIF0 ELEM 'TYPE' ;
*NMT1 = DIME LMOT1 ;
*SI (NMT1 > 1) ;
*   'ERREUR' 'Maillage comprenant plusieurs types' ;
*   'QUITTER' @ISOSURF    ;
*SINON ;
*   MOT1 = EXTR LMOT1 1 ;
*   SI ( NEG MOT1 'TET4' ) ;
*     'ERREUR' 'Maillage non constitué de TET4' ;
*     'QUITTER' @ISOSURF    ;
*   FINSI ;
*FINSI ;
*
LIS1 = ORDO LIS1 ;
*
* Calcul du paramètre d'élimination
*
lmm = 'PROG' 1.D-30 ;
'REPETER' idim ('VALEUR' 'DIME') ;
   cc = 'COORDONNEE' &idim MASSIF0 ;
   dmm = '-' ('MAXIMUM' cc) ('MINIMUM' cc) ;
   lmm = 'ET' lmm ('PROG' dmm) ;
'FIN' idim ;
EPSI1 = 0.000001 '*' ('MAXIMUM' lmm) ;
*
'SI' lgibi ;
*
* Début de l'algorithme codé en gibiane
*
* nombre elements du maillage
NB1 = NBEL MASSIF0 ;
* Boucle sur les isovaleurs
REPETER BOUS1 NIS1 ;
ISV1 = EXTR LIS1 &BOUS1 ;
*
* Nombre elements avec triangle isovaleur
NBIV1 = 0 ;
* Boucle sur les elements constitutifs du maillage
REPETER BOUC1 NB1 ;
* On extrait le nieme element
 TET1 = MASSIF0 ELEM &BOUC1 ;
 TET2 = CHANGER 'POI1' TET1 ;
* On extrait les valeurs du champoint et les points
 P1 = TET2 POIN 1 ;
 P2 = TET2 POIN 2 ;
 P3 = TET2 POIN 3 ;
 P4 = TET2 POIN 4 ;
 VAL1 = EXTR HANA1 SCAL P1 ;
 VAL2 = EXTR HANA1 SCAL P2 ;
 VAL3 = EXTR HANA1 SCAL P3 ;
 VAL4 = EXTR HANA1 SCAL P4 ;
 LIT1 = PROG VAL1 VAL2 VAL3 VAL4 ;
 MAX1 = MAXI LIT1 ;
 MIN1 = MINI LIT1 ;
* Rentrer dans element si ISV1 (isovaleur) comprise
SI ((ISV1 > MAX1) OU (ISV1 < MIN1)) ;
* Pas d'isovaleur dans cet element
 STRI1 = FAUX ;
SINON ;
* Il y a une isovaleur dans l'element
 EXP0 = FAUX ;
 STRI1 = FAUX ;
* On ordonne les valeurs aux points (et les points)
 REPETER BOUC2 ;
   SI (VAL1 &lt;EG VAL2) ;
      SI (VAL2 &lt;EG VAL3) ;
         SI (VAL3 &lt;EG VAL4) ;
            QUITTER BOUC2 ;
         SINON ;
            VAX1 = VAL4 ;
            PAX1 = P4   ;
            VAL4 = VAL3 ;
            P4   = P3 ;
            VAL3 = VAX1 ;
            P3   = PAX1 ;
         FINSI ;
      SINON ;
         VAX1 = VAL3 ;
         PAX1 = P3   ;
         VAL3 = VAL2 ;
         P3   = P2 ;
         VAL2 = VAX1 ;
         P2   = PAX1 ;
      FINSI ;
   SINON ;
      VAX1 = VAL2 ;
      PAX1 = P2   ;
      VAL2 = VAL1 ;
      P2   = P1 ;
      VAL1 = VAX1 ;
      P1   = PAX1 ;
   FINSI ;
 FIN BOUC2 ;
* On a fini d'ordonner les valeurs
* On teste si l'isovaleur correspond à un
* noeud du tetrahèdre
 NPEQ1 = 0 ;
 SI (ISV1 EGA VAL1) ;
  NPEQ1 = NPEQ1 + 1 ;
 FINSI ;
 SI (ISV1 EGA VAL2) ;
  NPEQ1 = NPEQ1 + 1 ;
 FINSI ;
 SI (ISV1 EGA VAL3) ;
  NPEQ1 = NPEQ1 + 1 ;
 FINSI ;
*gounand SI (ISV1 EGA VAL3) ;
 SI (ISV1 EGA VAL4) ;
  NPEQ1 = NPEQ1 + 1 ;
 FINSI ;
* NPEQ1 indique le nombre de noeuds du tetrahèdre
* qui correspondent à l'isovaleur
 SI (NPEQ1 EGA 4) ;
*  Isosurface = 4 faces du tetrahedre
   NBIV1 = NBIV1 + 1 ;
   STRI1 = VRAI ;
   MAIL1 = CHANGER 'FACE' TET1 ;
 SINON ;
   SI (NPEQ1 EGA 3) ;
*  Isosurface = 1 face du tetrahedre
     STRI1 = VRAI ;
     NBIV1 = NBIV1 + 1 ;
     SI (ISV1 EGA VAL1) ;
*  Point P1
       PIN1 = P1 ;
     SINON ;
*  Point P4 ;
       PIN1 = P4 ;
     FINSI ;
*  Toujours P3 et P2 inclus
     PIN2 = P2 ;
     PIN3 = P3 ;
   SINON ;
     SI (NPEQ1 EGA 2) ;
*  Isosurface = 1 segment ou 1 triangle
*   appuye sur 2 points du tetrahedre
       SI ((ISV1 EGA VAL1) OU (ISV1 EGA VAL4)) ;
*  Isosurface = 1 segment (on ne fait rien)
       SINON ;
*  Isosurface = 1 triangle appuye sur 2 points tetra
          NBIV1 = NBIV1 + 1 ;
          STRI1 = VRAI ;
*  point intersection (P1,P4)
          PIN1 = P2 ;
          PIN2 = P3 ;
          VEC1 = P4 MOINS P1 ;
          DVE1 = (ISV1 - VAL1) / (VAL1 - VAL4) ;
          VEC2 = DVE1 * VEC1 ;
          PIN3 = P1 MOINS VEC2 ;
       FINSI ;
     SINON ;
       SI (NPEQ1 EGA 1) ;
*  Isosurface = 1 point ou 1 triangle
*   appuye sur 1 point du tetrahedre
         SI ((ISV1 EGA VAL1) OU (ISV1 EGA VAL4)) ;
*  Isosurface = 1 point (on ne fait rien)
         SINON ;
*  Isosurface = 1 triangle appuye sur 1 point
*   du tetrahedre
           NBIV1 = NBIV1 + 1 ;
           STRI1 = VRAI ;
           SI (ISV1 EGA VAL2) ;
*  Le point appuye est P2
*  les autres points intersection (P1,P3)
             PIN1 = P2 ;
             VEC1 = P3 MOINS P1 ;
             DVE1 = (ISV1 - VAL1) / (VAL1 - VAL3) ;
             VEC2 = DVE1 * VEC1 ;
             PIN2 = P1 MOINS VEC2 ;
           SINON ;
*  Le point appuye est P3
*  les autres points intersection (P2,P4)
             PIN1 = P3 ;
             VEC1 = P4 MOINS P2 ;
             DVE1 = (ISV1 - VAL2) / (VAL2 - VAL4) ;
             VEC2 = DVE1 * VEC1 ;
             PIN2 = P2 MOINS VEC2 ;
           FINSI ;
*  Il y a toujours un point intersection (P1,P4)
           VEC1 = P4 MOINS P1 ;
           DVE1 = (ISV1 - VAL1) / (VAL1 - VAL4) ;
           VEC2 = DVE1 * VEC1 ;
           PIN3 = P1 MOINS VEC2 ;
         FINSI ;
       SINON ;
*  Isosurface = 1 triangle quelconque section
*   du tetrahedre ou un quadrangle (= 2 triangles)
         NBIV1 = NBIV1 + 1 ;
         STRI1 = VRAI ;
         SI (ISV1 < VAL2) ;
*  Isosurface entre V1 et V2
*  Points intersection (P1,P2) (P1,P3)
*  Un seul triangle
           VEC1 = P2 MOINS P1 ;
           DVE1 = (ISV1 - VAL1) / (VAL1 - VAL2) ;
           VEC2 = DVE1 * VEC1 ;
           PIN1 = P1 MOINS VEC2 ;
           VEC1 = P3 MOINS P1 ;
           DVE1 = (ISV1 - VAL1) / (VAL1 - VAL3) ;
           VEC2 = DVE1 * VEC1 ;
           PIN2 = P1 MOINS VEC2 ;
         SINON ;
           SI (ISV1 > VAL3) ;
*  Isosurface entre V3 et V4
*  Points intersection (P3,P4) (P2,P4)
*  Un seul triangle
             VEC1 = P4 MOINS P3 ;
             DVE1 = (ISV1 - VAL3) / (VAL3 - VAL4) ;
             VEC2 = DVE1 * VEC1 ;
             PIN1 = P3 MOINS VEC2 ;
             VEC1 = P4 MOINS P2 ;
             DVE1 = (ISV1 - VAL2) / (VAL2 - VAL4) ;
             VEC2 = DVE1 * VEC1 ;
             PIN2 = P2 MOINS VEC2 ;
           SINON ;
*  Isosurface entre V2 et V3 ;
*  Points intersection (P1,P3) (P2,P3) (P2,P4)
*  Un quadrangle = 2 triangles
             EXP0 = VRAI ;
             VEC1 = P3 MOINS P1 ;
             DVE1 = (ISV1 - VAL1) / (VAL1 - VAL3) ;
             VEC2 = DVE1 * VEC1 ;
             PIN0 = P1 MOINS VEC2 ;
             VEC1 = P3 MOINS P2 ;
             DVE1 = (ISV1 - VAL2) / (VAL2 - VAL3) ;
             VEC2 = DVE1 * VEC1 ;
             PIN1 = P2 MOINS VEC2 ;
             VEC1 = P4 MOINS P2 ;
             DVE1 = (ISV1 - VAL2) / (VAL2 - VAL4) ;
             VEC2 = DVE1 * VEC1 ;
             PIN2 = P2 MOINS VEC2 ;
           FINSI ;
         FINSI ;
*  Il y toujours un point intersection
*  entre P1 et P4
         VEC1 = P4 MOINS P1 ;
         DVE1 = (ISV1 - VAL1) / (VAL1 - VAL4) ;
         VEC2 = DVE1 * VEC1 ;
         PIN3 = P1 MOINS VEC2 ;
       FINSI ;
     FINSI ;
   FINSI ;
 FINSI ;
FINSI ;
*
* Fabrication de la surface si elle existe
SI STRI1 ;
* Premier triangle
  LS1 = DROIT 1 PIN1 PIN2 ;
  LS2 = DROIT 1 PIN2 PIN3 ;
  LS3 = DROIT 1 PIN3 PIN1 ;
  CONT1 = LS1 ET LS2 ET LS3 ;
  MAIL1 = COUL VERT (SURF CONT1 'PLANE') ;
  SI EXP0 ;
* Second triangle
    LS1 = DROIT 1 PIN3 PIN0 ;
    LS2 = DROIT 1 PIN0 PIN1 ;
    LS3 = DROIT 1 PIN1 PIN3 ;
    CONT1 = LS1 ET LS2 ET LS3 ;
    MAIL2 = COUL VERT (SURF CONT1 'PLANE') ;
    MAIL1 = MAIL1 ET MAIL2 ;
  FINSI ;
  SI (NBIV1 EGA 1) ;
    MAIT1 = MAIL1 ;
  SINON ;
    MAIT1 = MAIT1 ET MAIL1 ;
  FINSI ;
FINSI ;
FIN BOUC1 ;
* Elimination points multiples du maillage de surface
* obtenu par somme de triangles independants
ELIM EPSI1 MAIT1 ;
CHP1 = MANU 'CHPO' MAIT1 1 'SCAL' ISV1 'NATURE' 'DISCRET' ;
SI (NIS1 > 1) ;
   SI (&BOUS1 EGA 1) ;
     MAIF1 = MAIT1 ;
     CHPF1 = CHP1 ;
   SINON ;
     MAIF1 = MAIF1 ET MAIT1 ;
     CHPF1 = CHPF1 ET CHP1 ;
   FINSI ;
SINON ;
  MAIF1 = MAIT1 ;
  CHPF1 = CHP1 ;
FINSI ;
FIN BOUS1 ;
*
* Fin de l'algorithme en gibiane
* et début de la version fortran
'SINON' ;
   hana2 = 'CHANGER' 'CHAM' hana1 massif0 ;
   lmail = vrai ;
   'REPETER' BOUS1 NIS1 ;
      ISV1 = EXTR LIS1 &BOUS1 ;
*dbg      'MESSAGE' ('CHAINE' 'isv1=' isv1) ;
      MAIT1 = 'ISOV' hana2 ISV1 ;
*dbg      ymin = 'MINIMUM' ('COORDONNEE' 2 MAIT1) ;
*dbg      'MESSAGE' ('CHAINE' 'ymin=' ymin) ;
      'SI' ('>' ('NBEL' mait1) 0) ;
*sg: normalement inutile depuis fiche 8152   'ELIM' EPSI1 MAIT1 ;
         CHP1 = 'MANU' 'CHPO' MAIT1 1 'SCAL' ISV1 'NATURE' 'DISCRET' ;
         'SI' lmail ;
            MAIF1 = MAIT1 ;
            CHPF1 = CHP1 ;
            lmail = FAUX ;
         'SINON' ;
            MAIF1 = MAIF1 ET MAIT1 ;
            CHPF1 = CHPF1 ET CHP1 ;
         'FINSI' ;
      'FINSI' ;
   'FIN' BOUS1 ;
'FINSI' ;
*
'FINPROC' MAIF1 CHPF1 ;
*
********************************************************************