* fichier : isosurf.dgibi
************************************************************************
* Section : Mathematiques Fonctions
************************************************************************
************************ CAS TEST : isosurf.dgibi ******************
GRAPH = FAUX  ;
'SAUT' 'PAGE' ;
*-------------------------------------------------------------------
* TEST @ISOSURF
*                ISOSURFACES POUR MAILLAGE DE TETRAHEDRES
********************************************************************
*  Test de la procedure qui extrait les isosurfaces dont les valeurs 
*  sont listées dans une liste (LIS1) d'un champoint (HANA1) appuyé 
*  sur un maillage (MASSIF0) exclusivement constitué de tetrahèdres 
*  (TET4).
*
*  Le résultat final est constitué du maillage surfacique (TRI3)
*  regroupant l'ensemble des isosurface MAIF1 et du champoint CHPF1 
*  des isovaleurs LIS1 appuyées sur MAIF1.
*
*  Si GRAPH est faux, seule l'isosurface de valeur 1500 est extraite
*  et un test de non regression est realise sur la base de 
*  l'apartenance d'un point de coordonnees analytiques a la surface 
*  extraite.
*
*  Si GRAPH est vrai, le meme test de non regression est realise,
*  mais 6 isosurfaces dont les valeurs sont listees dans LIS1 sont
*  tracees. 
********************************************************************
*
'OPTION' 'ECHO' 1 ;
'TITRE' 'Isosurface Tetra'   ;
OPTI DIME 3 ELEM TET4   ;
OPTI ISOV SURFACE ;
*
*-------------------------------------------------------------------
*--------------------- Création du maillage 3D ---------------------
*
EPSI1 = 0.000001 ;
*----------------------------------------- 
*          Dimensions de base
*-----------------------------------------
LX1 = 10.D0 ;
LY1 = 10.D0 ;
LZ1 = 10.D0 ;
*-----------------------------------------
*               Points
*-----------------------------------------
* Points tetra
P1 = 0.D0 0.D0 0.D0 ;
P2 = LX1 0.D0 0.D0 ;
P3 = LX1 LY1 0.D0 ;
P4 = (LX1 / 2.D0) (LY1 / 4.D0) LZ1 ;
*-----------------------------------------
*             Lignes base
*-----------------------------------------
ND1 = 13 ;
LXY1 = DROIT ND1 P1 P2 ;
LXY2 = DROIT ND1 P2 P3 ;
LXY3 = DROIT ND1 P3 P1 ;
COT1 = LXY1 ET LXY2 ET LXY3 ;
*
LZZ1 = DROIT ND1 P1 P4 ;
LZZ2 = DROIT ND1 P2 P4 ;
LZZ3 = DROIT ND1 P3 P4 ;
COT2 = LXY1 ET LZZ2 ET (INVE LZZ1) ;
COT3 = LXY2 ET LZZ3 ET (INVE LZZ2) ;
COT4 = LXY3 ET LZZ1 ET (INVE LZZ3) ;
*-----------------------------------------------------------
*                         Surfaces
*-----------------------------------------------------------
SUB1 = SURF COT1 'PLANE' ; 
SUZ1 = SURF COT2 'PLANE' ;
SUZ2 = SURF COT3 'PLANE' ;
SUZ3 = SURF COT4 'PLANE' ;
SUT1 = SUB1 ET SUZ1 ET SUZ2 ET SUZ3 ;
ELIM EPSI1 SUT1 ;
*-------------------------------------------------------------------
*                           -- VOLUMES --
*-------------------------------------------------------------------
MASSIF0 = COUL ROUG (VOLU SUT1) ;
SI GRAPH ;
   TRAC CACH MASSIF0 ;
FINSI ;
*-------------------------------------------------------------------
*=====================================================================
*                             -- MAILLAGES QUAF --
*=====================================================================
QFTOT = CHANGE MASSIF0 QUAF ;
*
*=====================================================================
*                                 -- MODELE--
*=====================================================================
MODHYB = 'MODELE' QFTOT 'DARCY' 'ANISOTROPE'     ;
HYSOM = 'DOMA' MODHYB 'SOMMET'       ;
XXS YYS ZZS = 'COOR' HYSOM ;
*
*=====================================================================
* Valeurs charges imposees
GRA1  = 2000.D0  ;
*
* Champoint
PI2 = 2.0D0 * 3.1414D0 ;
HANA1 = (COS (PI2 * (XXS - (LX1 / 2.D0))))
      * (COS (PI2 * (XXS - (LX1 / 2.D0))))
      * (COS (PI2 * (YYS - (LY1 / 2.D0))))
      * (COS (PI2 * (YYS - (LY1 / 2.D0))))
      * (COS (PI2 * (ZZS - (LZ1 / 2.D0)))) 
      * (COS (PI2 * (ZZS - (LZ1 / 2.D0))))
      * GRA1  ;
*
*------------------Tracer champoint-----------------------------
SI GRAPH ;
   TRAC CACH HANA1 MASSIF0 ;
FINSI ;
* ===========================================================
* Liste des isovaleurs a extraire
LIS1 = PROG 1500.0 ;
SI GRAPH ;
   LIS1 = PROG 1900.0 1700.0 1500.0 800.0 1040.0 1200.0 ;
FINSI ;
*
* Procedure @isosurf
MAIF1 CHPF1 = @ISOSURF MASSIF0 LIS1 HANA1 ;
*
SI GRAPH ;
   TRAC CHPF1 MAIF1 ;
FINSI ;
*
CONT1 = COT1 ET COT2 ET COT3 ET COT4 ;
VAL1 = MINI HANA1 ;
CHP2 = MANU 'CHPO' CONT1 1 'SCAL' VAL1 'NATURE' 'DISCRET' ;
MAT1 = MAIF1 ET CONT1 ;
CHT1 = CHPF1 ET CHP2 ;
SI GRAPH ;
   TRAC CACH FACE CHT1 MAT1 ;
FINSI ;
*
* Verification de l'algoritme
* Un point de coordonnees XT1, YT1, ZT1 identifie comme
* appartenant a la surface d'isovaleur 1500.0D0 definie 
* dans LIS1 est il reelement contenu dans l'isosurface extraite ?
VAL0 = (COS (PI2 * ((LY1 / 4.D0) - (LY1 / 2.D0))))
      * (COS (PI2 * ((LY1 / 4.D0) - (LY1 / 2.D0)))) ;
VAL1 = (1500.0D0 / (GRA1 * VAL0)) ** 0.5D0 ;
VAL2 = ((1.D0 - (VAL1 * VAL1)) ** 0.5D0) / VAL1 ;
VAL3 = ATG (VAL2) ;
ZT1 = (VAL3 / PI2 ) + (LZ1 / 2.D0) ;
XT1 = LX1 / 2.D0 ;
YT1 = LY1 / 4.D0 ;
PV1 = XT1 YT1 ZT1 ;
*
OPTI ERREUR IGNORE ;
GEO1 = (MAIF1 ELEM 'CONTENANT' PV1) ELEM 1  ;
OPTI ERREUR NORMAL ;
MOT1 = TYPE GEO1 ;
SI (EGA MOT1 'ANNULE') ;
   MESS 'Point de reference hors isosurface' ;
   'ERREUR' 5 ;
FINSI ;
*
FIN ;
************************************************************************