$$$$ CH_THETX
* CH_THETX  PROCEDUR  JB251061  21/06/14    21:15:01     11034          
* =============================================================================
*PROCEDURE DE CALCUL DU CHAMP THETA/PI EN XFEM
*---------------------------------------------
*
* DESCRIPTION : PROCEDURE ANALOGUE A CH_THETA MAIS POUR LES ELEMENTS XFEM.
*               DETERMINE UN CHPOINT DE NORME CONSTANTE A L'INTERIEUR D'UNE
*               COURONNE ENTOURANT LE FRONT DE FISSURE ET NUL A L'EXTERIEUR DE
*               CETTE COURONNE. LE CHPOINT THETA AINSI DEFINI REPRESENTE LA
*               DIRECTION DE PROPAGATION EVENTUELLE DE LA FISSURE.
* =============================================================================  
DEBP CH_THETX SUPTAB*TABLE OBJUTI*'TABLE' BOOL*'TABLE' ;

* RECUPERATION DE CERTAINES DONNEES
GDIME = VALE 'DIME' ;
GMODE = VALE 'DIME' ;
MOD_MEC = OBJUTI.'MOD_MEC' ;
MUI = EXTR (EXTR MOD_MEC 'DEPL') (LECT 1 PAS 1 GDIME) ;
M_FRONT = OBJUTI.'FRONT' ;


* I - COMPATIBILITE DES DONNEES
* -----------------------------

SI (NON (EXIS SUPTAB 'MAILLAGE')) ;
    MESS 'ERREUR : ON N A PAS TROUVE DANS LA' ;
    MESS '         TABLE L OBJET MAILLAGE' ;
    QUIT CH_THETX ;
SINON ;
    MAILLAGE = SUPTAB.'MAILLAGE' ;            
FINSI ; 
SI (NON (EXIS SUPTAB 'FRONT_FISSURE')) ;
    MESS 'ERREUR : LE FRONT DE LA FISSURE N EST PAS DONNE' ;
    QUIT CH_THETX ;
SINON ;
    PFISS = SUPTAB.'FRONT_FISSURE' ;
    SI (EGA (TYPE PFISS) 'POINT') ;
        PFISS1 = MANU 'POI1' PFISS ;
    SINON ;
        PFISS1 = CHAN 'POI1' PFISS ;
    FINSI ;
FINSI ;               
SI ((NON (EXIS SUPTAB 'PSI')) OU (NON (EXIS SUPTAB 'PHI'))) ;  
    MESS 'ERREUR : CHPOINT PSI et PHI NON FOURNIS' ;
    QUIT CH_THETX ;
SINON ;      
    PSI1 = SUPTAB.'PSI' ;
    PHI1 = SUPTAB.'PHI' ;
FINSI ; 
SI (EXIS SUPTAB 'COUCHE') ;
    COUCHE = SUPTAB.'COUCHE' ;
SINON ;
    SI (EXIS SUPTAB 'CHAMP_THETA') ;
        THETA = SUPTAB.'CHAMP_THETA' ;
    SINON ;
        MESS 'ERREUR : IL FAUT SOIT LE NOMBRE DE COUCHES SOIT UN CHAMP THETA' ;
    FINSI ;
FINSI ;

* ON DETERMINE SI UNE OU DEUX LEVRES SONT MODELISEES
NB_LEVRES = 2 ;
* SI UNE DEMI-EPROUVETTE EST MODELISEE, 
* ON ESPERE QUE PHI(X)>0 OU <0 QQSOIT X
MINY1 = MINI PHI1 ;  MAXY1 = MAXI PHI1 ;
SI ((MAXY1*MINY1) <EG 0.) ;
    MINY1A = MINI MINY1 MAXY1 'ABS' ;
    MAXY1A = MAXI MINY1 MAXY1 'ABS' ;
    SI ((MINI MINY1 MAXY1 'ABS') < ((MAXI MINY1 MAXY1 'ABS')*1.E-3)) ;
        NB_LEVRES = 1 ;
    FINSI ;
FINSI ;

* ELEMENTS POUR LESQUELS ON VEUT FAIRE LE CALCUL
SI (EGA GDIME 3) ;
    SI (EXIS SUPTAB 'NOEUDS_AVANCES') ;
        NOAV = SUPTAB.'NOEUDS_AVANCES' ;
    SINON ;
        NOAV = CHAN 'POI1' M_FRONT ;
    FINSI ;
    ELAV = VIDE 'MAILLAGE' ;
    REPE INOEU (NBEL NOAV) ;
        POIN1 = (NOAV ELEM &INOEU) POIN 1 ;
        ELAV = ELAV ET (MAILLAGE ELEM 'CONTENANT' POIN1 'TOUS' ) ;
    FIN INOEU ;
FINSI ;
    
    
* II - CREATION DES CHAMPS THETA
* ------------------------------

* II.1 - CAS 2D
* *************

SI (EGA GDIME 2) ;
    
    SI (EXIS SUPTAB 'COUCHE') ;
    
* ELEMENT CONTENANT LE FRONT DE FISSURE
        EFISS0 = MAILLAGE ELEM 'CONTENANT' PFISS ;
      
* TEST POUR VOIR SI ON EST SUR UN NOEUD, DANS UN ELEMENT OU SUR UN BORD
        MOD0 = MODE EFISS0 'MECANIQUE' 'ELASTIQUE' 'XQ4R' ;
        PSI10 = CHAN 'CHAM' MOD0 (REDU PSI1 EFISS0) 'STRESSES' ;
        PHI10 = CHAN 'CHAM' MOD0 (REDU PHI1 EFISS0) 'STRESSES' ;
        LOGPSI10 = ((MAXI PSI10) *  (MINI PSI10)) < 0. ;
        LOGPHI10 = ((MAXI PHI10) *  (MINI PHI10)) < 0. ;
        SI(LOGPSI10 ET LOGPHI10) ;
            MBOUGER = EFISS0 ;
        SINON ;
            MBOUGER = EFISS0 POIN 'PROCH' PFISS ;
            SI((LOGPSI10 OU LOGPHI10) ET (EGA COUCHE 0)) ;
                MESS 'ATTENTION : UTILISATION DE 0 COUCHE + XFEM DELICATE' ;
                MESS '            RESULTAT PROBABLEMENT DE MAUVAISE PRECISION' ;
            FINSI ;
        FINSI ;

* CREATION DU MAILLAGE A BOUGER : MBOUGER
        REPE BCOUCH COUCHE ;
            MBOUGER = MAILLAGE ELEM 'APPUYE' 'LARG' MBOUGER ;
        FIN BCOUCH ;
* CREATION DU MAILLAGE OU THETA N'EST PAS NUL : MAIL
        MAIL = ELEM MAILLAGE 'APPU' 'LARG' MBOUGER ;
* PSI ET PHI DOIVENT ETRE ENTIREMENT DEFINIS SUR MAIL
        MINTER = INTE (CHAN 'POI1' MAIL) (EXTR PSI1 'MAILLAGE') ;
        NBEL1  = NBEL (DIFF (CHAN 'POI1' MAIL) MINTER ) ;
        SI (NBEL1 NEG 0) ;
            MESS NBEL1 'POINTS DES COUCHES OU PSI ET PHI SONT INDEFINIS' ;
            ERRE 21 ;
        FINSI ;
        
* DETERMINATION DU CHAMP THETA :
* ON DEFINIT UN MODELE THERMIQUE POUR CALCULER LES GRADIENTS
        MOD7 = MODE MAIL 'THERMIQUE' ;
        GPSI = GRAD MOD7 (NOMC 'T' PSI1 'NATU' 'DIFFUS') ;
        GPSI = CHAN GPSI 'TYPE' 'SCALAIRE' ;
        GPHI = GRAD MOD7 (NOMC 'T' SUPTAB.'PHI' 'NATU' 'DIFFUS') ;
        GPHI = CHAN GPHI 'TYPE' 'SCALAIRE' ;
* THETA EST PERPENDICULAIRE A GRAD(PHI) ET DANS LE SENS DE GRAD(PSI)
        NCOMP = EXTR GPSI 'COMP' ;
        THETA = PVEC GPHI NCOMP MUI ;
* ON UTILISE SIGN POUR S'ASSURER QUE THETA ET GRAD(PSI) SONT DANS LE MEME SENS
        SIGN1 = SIGN (PSCA THETA GPSI MUI NCOMP) ;
        THETA = SIGN1 * THETA ;
* DIRECTION DE PROPAGATION DIR1 SUR LE FRONT AVANT DE CHANGER THETA EN CHPO
        DIR1 = PROI M_FRONT (CHAN 'NOEUD' MOD7 THETA) ;
* ON CHANGE THETA EN CHPO ET ON REDUIT SON SUPPORT A MBOUGER
        THETA = CHAN 'CHPO' MOD7 THETA 'MOYE' ;
        THETA = REDU THETA MBOUGER ;
        
    SINON ;
    
* ON INTERPOLE THETA SUR LE FRONT POUR DIR1
        DIR1 = INT_COMP MAILLAGE THETA PFISS1 ;
    
    FINSI ;
    
* ON NORMALISE DIR1 ET ON CALCULE DIR2
    NDIR1 = (PSCA DIR1 DIR1 MUI MUI)**0.5 ;
    DIR1 = DIR1 / NDIR1 ;
    DIR2 = PVEC DIR1 MUI MUI ;
    
* ON NORMALISE THETA ET ON MULTIPLIE PAR 2. / NB_LEVRES
    NTHETA = (PSCA THETA THETA MUI MUI)**0.5 ;
    THETA = THETA / NTHETA * 2. / NB_LEVRES ;
    
* STOCKAGE
    TTHETA = TABL ;
    FINTER = MANU 'CHPO' PFISS 1 'SCAL' 1. 'NATURE' 'DISCRET' ;
    TTHETA.FINTER = THETA ;
    
FINSI ;


* II.2 - CAS 3D
* *************

SI (EGA GDIME 3) ;
    
* ELEMENTS CONTENANT LE FRONT DE FISSURE
    MAILPSI = EXTR PSI1 'MAILLAGE' ;
    MAILPSI = ELEM MAILLAGE 'APPUYE' 'STRICTEMENT' MAILPSI ;
    PSI1CHAM = CHAN 'CHAM' PSI1 MAILPSI ;
    PHI1CHAM = CHAN 'CHAM' PHI1 MAILPSI ;
    TAB_ELF = (ELEM PSI1CHAM 'EGINFE' 0.) INTE (ELEM PSI1CHAM 'EGSUPE' 0.) ;
    EFISS2 = (ELEM PHI1CHAM 'EGINFE' 0.) INTE (ELEM PHI1CHAM 'EGSUPE' 0.) ;
    EFISS0 = TAB_ELF INTE EFISS2 ;
    DX0 = ((MESU EFISS0 'VOLUME') / (NBEL EFISS0)) ** (1./3.) ;
    EPSP = 1.E-4 * DX0 ;
    EPSN = -1.E-4 * DX0 ;
    TAB_ELF = (ELEM PSI1CHAM 'EGINFE' EPSP) INTE (ELEM PSI1CHAM 'EGSUPE' EPSN) ;
    EFISS2 = (ELEM PHI1CHAM 'EGINFE' EPSP) INTE (ELEM PHI1CHAM 'EGSUPE' EPSN) ;
    EFISS0 = TAB_ELF INTE EFISS2 ;
    
* CREATION DU MAILLAGE A BOUGER : MBOUGER
    MBOUGER = EFISS0 ;
    SI (COUCHE > 1) ;                                                  
        REPE ICOUCH (COUCHE - 1) ;                                            
            MBOUGER = ELEM MAILLAGE 'APPUYE' 'LARG' MBOUGER ;              
        FIN  ICOUCH ;
    FINSI ;
* CREATION DU MAILLAGE OU THETA N'EST PAS NUL : MAIL
    MAIL = ELEM MAILLAGE 'APPU' 'LARG' MBOUGER ;
* ON VERIFIE QUE PSI ET PHI SONT ENTIEREMENT DEFINIS SUR MAIL
    MINTER = (CHAN 'POI1' MAIL) INTE (EXTR PSI1 'MAILLAGE') ;
    NBEL1 = NBEL (DIFF (CHAN 'POI1' MAIL) MINTER ) ;
    SI (NBEL1 NEG 0) ; 
        MESS NBEL1 'POINTS DES COUCHES OU PSI ET PHI SONT INDEFINIS' ;
        ERRE 21 ;
    FINSI ;
* ON SE RESTREINT AU DOMAINE DE DEFINITION DE PSI ET PHI
    MAIL = MAIL INTE MAILPSI ;


* ON COMMENCE PAR CREER LE CHAMP THETA GLOBAL :
* *********************************************

* ON DEFINIT UN MODELE THERMIQUE POUR CALCULER LES GRADIENTS
    MOD7 = MODE MAIL 'THERMIQUE' ;
    GPSI = GRAD MOD7 (NOMC 'T' PSI1 'NATU' 'DIFFUS') ;
    GPSI = CHAN GPSI 'TYPE' 'SCALAIRE' ;
    GPHI = GRAD MOD7 (NOMC 'T' SUPTAB.'PHI' 'NATU' 'DIFFUS') ;
    GPHI = CHAN GPHI 'TYPE' 'SCALAIRE' ;
    NCOMP = EXTR GPSI 'COMP' ;
* TAW = GRAD(PSI)^GRAD(PHI)
    TAW1 = PVEC GPSI GPHI NCOMP NCOMP NCOMP ;
* THETA = GRAD(PHI)^TAW
    THETA = PVEC GPHI TAW1 NCOMP NCOMP MUI ;
* ON CHANGE THETA EN CHPO
    THETA =  CHAN 'CHPO' MOD7 THETA 'MOYE' ;
* ON NORME THETA, EN FAISANT ATTENTION QUE LE CHAMP PEUT ETRE QUASI-NUL
* PAR ENDROITS => ON UTILISE MASQ
    NTHETA = (PSCA THETA THETA MUI MUI)**0.5 ;
    MTHETA = MASQ NTHETA 'SUPERIEUR' 1.E-2 ;
    BTHETA = BORN NTHETA 'SCAL' 'MINIMUM' 1.E-2 ;
    THETA = CHAN (THETA / BTHETA * MTHETA) 'ATTRIBUT' 'NATURE' 'DIFFUS' ;
    OBJUTI.'V1' = THETA ;
* ON REDUIT THETA A MBOUGER
    THETA = REDU THETA MBOUGER ;
* ON FAIT DE MEME AVEC GRAD(PHI), MAIS IL N'EST A PRIORI JAMAIS NUL DONC PAS
* BESOIN DE MASQ
    GPHI = CHAN 'CHPO' MOD7 GPHI 'MOYE' ;
    NGPHI = (PSCA GPHI GPHI NCOMP NCOMP)**0.5 ;
    GPHI = CHAN (GPHI / NGPHI) 'ATTRIBUT' 'NATURE' 'DIFFUS' ;
    OBJUTI.'V2' = NOMC NCOMP MUI GPHI ;
    GPHI = REDU GPHI MBOUGER ;
* ON FAIT DE MEME AVEC TAW, LUI AUSSI PEUT ETRE QUASI-NUL
    TAW1 = CHAN 'CHPO' MOD7 TAW1 'MOYE' ;
    NTAW1 = (PSCA TAW1 TAW1 NCOMP NCOMP)**0.5 ;
    MTAW1 = MASQ NTAW1 'SUPERIEUR' 1.E-2 ;
    BTAW1 = BORN NTAW1 'SCAL' 'MINIMUM' 1.E-2 ;
    TAW1 = CHAN (TAW1 / BTAW1 * MTAW1) 'ATTRIBUT' 'NATURE' 'DIFFUS' ;
    OBJUTI.'V3' = NOMC NCOMP MUI TAW1 ;
    TAW1 = REDU TAW1 MBOUGER ;
* ON INTERPOLE LES 3 DIRECTIONS SUR LE FRONT
    DIR1 = INT_COMP MBOUGER THETA PFISS1 ;
    DIR2 = NOMC NCOMP MUI (INT_COMP MBOUGER GPHI PFISS1) ;
    DIR3 = NOMC NCOMP MUI (INT_COMP MBOUGER TAW1 PFISS1) ;
* ON LES NORME ET ON LES CHANGE EN DIFFUS
    NDIR1 = (PSCA DIR1 DIR1 MUI MUI)**0.5 ;
    NDIR2 = (PSCA DIR2 DIR2 MUI MUI)**0.5 ;
    NDIR3 = (PSCA DIR3 DIR3 MUI MUI)**0.5 ;
    DIR1 = CHAN (DIR1 / NDIR1) 'ATTRIBUT' 'NATURE' 'DIFFUS' ;
    DIR2 = CHAN (DIR2 / NDIR2) 'ATTRIBUT' 'NATURE' 'DIFFUS' ;
    DIR3 = CHAN (DIR3 / NDIR3) 'ATTRIBUT' 'NATURE' 'DIFFUS' ;
    
    
* ENSUITE ON IDENTIFIE CHAQUE ELEMENT AUTOUR DU FRONT A BOUGER :
* **************************************************************

    NPFISS1 = NBEL PFISS1 ; 
    TAB_ELF = TABL ;
    REPE IPFISS1 NPFISS1 ;
        PFI = PFISS1 POIN &IPFISS1 ;
* EXTRACTION DES ELEMENTS CONTENANT LE POINT COURANT
        ELI = MBOUGER ELEM 'CONTENANT' PFI 'TOUS' ;
        SI (EGA &IPFISS1 1) ;
            TAB_ELF.(1) = ELI ;
        SINON ;
            SI ((DIME TAB_ELF) > 1) ;
                INTE1 = ELI INTE TAB_ELF.((DIME TAB_ELF) - 1) 'NOVERIF' ;
                SI (NEG (NBEL INTE1) 0) ;
* S'IL Y A DES ELEMENTS EN COMMUN AVEC L'AVANT DERNIER INDICE DE LA TABLE
* ALORS ON LES ENLEVE
                    ELI = DIFF ELI INTE1 ;
                    SI (EGA (NBEL ELI) 0) ;
                        ITER IPFISS1 ;
                    FINSI ;
                FINSI ;
            FINSI ;
        FINSI ;
* ON REGARDE ENSUITE L'INTERSECTION AVEC LE DERNIER INDICE DE LA TABLE
        INTE2 = ELI INTE TAB_ELF.(DIME TAB_ELF) 'NOVERIF' ;
        SI (EGA (NBEL INTE2) 0) ;
* S'IL N'Y A PAS D'ELEMENT EN COMMUN ON A TERMINE
            TAB_ELF.((DIME TAB_ELF) + 1) = ELI ;
        SINON ;
* SINON IL FAUT SUPPRIMER/AJOUTER DES ELEMENTS
            ETEMP = DIFF TAB_ELF.(DIME TAB_ELF) INTE2 ;
            SI (NEG (NBEL ETEMP) 0) ;
* DANS CE CAS IL Y AVAIT DES ELEMENTS EN TROP DANS LE DERNIER INDICE
* => ON LES ENLEVE ET ON MET CE QUI RESTE DANS UN NOUVEL INDICE
                TAB_ELF.(DIME TAB_ELF) = ETEMP ;
                INTE2 = ELI INTE ETEMP 'NOVERIF' ;
                ELI = DIFF ELI INTE2 ;
                SI (NEG (NBEL ELI) 0) ;
                    TAB_ELF.((DIME TAB_ELF) + 1) = ELI ;
                FINSI ;
            SINON ;
* DANS CE CAS IL Y A DES ELEMENTS EN TROP DANS ELI
* => ON LES ENLEVE ET ON STOCKE CE QUI RESTE DANS UN NOUVEL INDICE
                ETEMP = DIFF ELI INTE2 ;
                SI (NEG (NBEL ETEMP) 0) ;
                    TAB_ELF.((DIME TAB_ELF) + 1) = ETEMP ;
                FINSI ;
            FINSI ;
        FINSI ;
    FIN IPFISS1 ;
    
    
* CREATION ET REMPLISSAGE DE LA TABLE TTHETA
* ******************************************

* FONCTIONS D'INTERPOLATION
    TINTER = TABL ;
* ON BOUCLE SUR LES ELEMENTS ENTOURANT LE FRONT
    REPE IELF (DIME TAB_ELF) ;
        MBOUGI = TAB_ELF.&IELF ;
        MINTER = MBOUGI INTE ELAV ;
        SI (EGA (NBEL MINTER) 0) ;
            ITER IELF ;
        FINSI ;
        FINTER = MANU 'CHPO' MBOUGI 1 'SCAL' 1. ;
* ON BOUCLE SUR LE NOMBRE DE COUCHES
        REPE ICOUCH (COUCHE - 1) ;
            MBOUGI = ELEM MAILLAGE 'APPUYE' 'LARG' MBOUGI ;
            FINTER = FINTER + (MANU 'CHPO' MBOUGI 1 'SCAL' 1.) ;
        FIN  ICOUCH ;
        TINTER.(TAB_ELF.&IELF) = FINTER ;
    FIN IELF ;
* TABLE TTHETA
    TTHETA = TABL ;
    TIND = INDE TINTER ;
    REPE IIND (DIME TIND) ;
        FINTER = TINTER.(TIND.&IIND) ;
        THETI = FINTER * THETA ;
        TTHETA.(TIND.&IIND) = THETI ;
    FIN IIND ;
    TTHETA.'GLOBAL' = THETA ;
    
    
* NORMALISATION
* *************

* PETIT FLOTTANT POUR QUE LE MOUVEMENT SOIT INFINITESIMAL
    EPS1 = 1.E-4 ;
    TIND = INDE TTHETA ;
    REPE ICHAM (DIME TIND) ;
        THETI = TTHETA.(TIND.&ICHAM) ;
* ON CREE UN CHAMP THETI DONT LA NORME EST PROCHE DE LA DIMENSION DES ELEMENTS
        NTHETI = (PSCA THETI THETI MUI MUI)**0.5 ;
        THETI = THETI * DX0 / (MAXI NTHETI) ;
        THETA2 = INT_COMP EFISS0 THETI PFISS ;
        PFISS2 = PFISS PLUS (EPS1*THETA2) ;
        AFISS2 = PFISS2 REGL PFISS 1 ;
        DELTA_A = NB_LEVRES * (MESU AFISS2 'SURF') / EPS1 ;
        THETI = THETI * 2. / DELTA_A ;
        TTHETA.(TIND.&ICHAM) = THETI ;
    FIN ICHAM ;

* ON REMPLACE THETA PAR LE CHAMP THETA GLOBAL
    THETA = TTHETA.'GLOBAL' ;

FINSI ;


* II.3 - STOCKAGE DU CHAMP THETA ET DES DIRECTIONS
* ************************************************

SI (EXIS SUPTAB 'COUCHE') ;
    SUPTAB.'CHAMP_THETA' = THETA ;
FINSI ;
SUPTAB.'TAB_THETA' = TTHETA ;

OBJUTI.'DIRECTION1' = DIR1 ;
OBJUTI.'DIRECTION2' = DIR2 ;
SI (EGA GDIME 3) ;
    OBJUTI.'DIRECTION3' = DIR3 ;
FINSI ;

FINP ;