$$$$ VARIHC
* VARIHC    PROCEDUR  MB234859  26/08/27    21:15:09     12630          
****************************************************
* Procedure de tracking par resolution d'un probleme
* de conduction thermique stationnaire
* Application au modele E-FEM
* Version  : 
*     Variable online global tracking 
* Creation : 
*     Francesco RICCARDI 22/08/2016
*     Benjamin RICHARD
* Contact :  
*     Francesco.Riccardi[at]cea.fr
*     Benjamin.Richard[at]cea.fr 
* Institution : 
*     CEA\DEN\DANS\DM2S\SEMT\EMSI
****************************************************
* Voir notice TRACKING pour plus d'informations
****************************************************
*
*******----DEBUT Procedure VARIHC----*******
*
DEBP VARIHC TAB2*'TABLE';


* Indice d'evolution
NPAS = dime tab2. deplacements;
NPAS = NPAS - 1;
*LIST NPAS;

* Modeles
MODEL1 = TAB2. 'MODELE';
MODELS = EXTR MODEL1 'ZONE';

* Modele et materiau EFEM
MODEL11 = MODELS. 1;
SURFACE1 = MODELS. 2;
MATr1 = TAB2. 'WTABLE'. 'CARACTERISTIQUES';
MATr11 = REDU MATr1 SURFACE1;
E11 = EXCO MATr11 YOUN;
NU11 = EXCO MATr11 NU;
P11 = EXCO MATr11 FT;
P12 = EXCO MATR11 XNX;
P13 = EXCO MATR11 XNY;
P14 = EXCO MATr11 IND1;
P15 = EXCO MATR11 XIE;
P16 = EXCO MATR11 YIE;

* Coordonnees de la geometrie
X0 = COOR 1 SURFACE1;
Y0 = COOR 2 SURFACE1;
X0S = CHAN 'CHAM' X0 MODEL11 'STRESSES';
X0S = CHAN 'TYPE' X0S 'CONTRAINTES PRINCIPALES';
X0S = CHAN 'COMP' 'SI11' X0S;
Y0S = CHAN 'CHAM' Y0 MODEL11 'STRESSES';
Y0S = CHAN 'TYPE' Y0S 'CONTRAINTES PRINCIPALES';
Y0S = CHAN 'COMP' 'SI11' Y0S;

* Nombre de root elements
NROOTS = TAB2. 'NROOTS'; 
*LIST NROOTS;

* Nombre maximum de fissures
MAXR = TAB2. 'NRMAX'; 

* Variables internes : NFLA
FLAG = EXCO (TAB2. 'VARIABLES_INTERNES'. NPAS) 'NFLA';
*TRAC FLAG MODEL11;
VFLAG = MAXI FLAG;
ELACT1 = FLAG ELEM 'EGALE' 1.;

* Maximum des contraintes principales
SIG1 = TAB2. 'CONTINUATION'. 'CONTRAINTES';
PRIN1 = PRIN SIG1 MODEL11;
SIGMA11 = EXCO 'SI11' PRIN1;
COS1X = NOMC 'SI11' (EXCO PRIN1 'COX1');
COS1Y = NOMC 'SI11' (EXCO PRIN1 'COY1');
SIGMAX = MAXI SIGMA11;

* 
*----- Debut algorithme de tracking --------------------------- 
* 
SI ((SIGMAX > 0.) ET (NPAS < (TAB2. 'NPAS_TRACKING')));


* Elements qui violent le critere de resistance
***dans la phase elastique lineaire VFLAG=0***
SI (VFLAG EGA 0.);
    ELSP11 = SIGMA11 ELEM 'MAXI';
SINON;
    MESH0 = TAB2. 'MESH'. (NPAS - 1);
    SURFRED = SURFACE1 DIFF MESH0;
    SIGRED = REDU SIGMA11 SURFRED;
***pour assurer que les points test n'appartiennent pas a MESH0***  
***cela peut se passer si on a des bloccages***   
    ELSP11 = SIGRED ELEM 'EGSUPE' (MAXI P11);
FINSI;
NTEST = NBEL ELSP11;
SI (NTEST > 0);
    SP11 = REDU SIGMA11 ELSP11;
FINSI;

TABISO0 = TABLE;
TABROOT = TABLE;
TABEROOT = TABLE;
SI (VFLAG EGA 0);
    TABEROOT. 1 = SIGMA11 ELEM 'MAXI';    
    TABROOT. 1 = BARY (TABEROOT. 1);
    SI (NROOTS EGA 1);
        TABROOT. 1 = TAB2. 'PROOTS'. 1;
        TABEROOT. 1 = TAB2. 'EROOTS'. 1;
    FINSI;
SINON;
    i = 0;
***NROOTS est surement >0***
    REPETER BOUCLE1 NROOTS;
        indice = i * 4;
        LXIN0 = TAB2. 'TABISO'. (indice + 1);
        LYIN0 = TAB2. 'TABISO'. (indice + 2);
        LXOUT0 = TAB2. 'TABISO'. (indice + 3);
        LYOUT0 = TAB2. 'TABISO'. (indice + 4);
        MESH0 = TAB2. 'MESH_ISO'. (i + 1);
* Nombre de points qui composent la fissure  
        DIMM0 = DIME LXIN0;     
***si tous les elements sont actives on touche rien*** 
        FLAGISO = REDU FLAG MESH0;
        MESHACT = FLAGISO ELEM 'SUPERIEUR' 0.;
        NACT = NBEL MESHACT;   
        SI ((NACT > 0) ET (NACT < DIMM0));
            SI ((NACT + 1) EGA DIMM0);
                LISTEN0 = DIMM0;
            SINON;
                LISTEN0 = ENTI (PROG (NACT + 1) PAS 1 DIMM0);
            FINSI;
* Partie activee de la fissure            
            LXIN0 = ENLE LXIN0 LISTEN0;
            LYIN0 = ENLE LYIN0 LISTEN0;
            LXOUT0 = ENLE LXOUT0 LISTEN0;
            LYOUT0 = ENLE LYOUT0 LISTEN0;
        FINSI;
        LXG0 = (LXIN0 + LXOUT0) / 2.; 
        LYG0 = (LYIN0 + LYOUT0) / 2.;  
        TABROOT. (i + 1) = TAB2. 'PROOTS'. (i + 1);
        TABEROOT. (i + 1) = TAB2. 'EROOTS'. (i + 1);
***le point root est le point output du dernier element active***
***si NACT=0 il coincide avec celui du pas precedent***           
        SI (NACT > 0);
            TABROOT. (i + 1) = (EXTR LXOUT0 NACT) (EXTR LYOUT0 NACT);
            TABEROOT. (i + 1) = (MESH0 ELEM 'CONTENANT'
            ((EXTR LXG0 NACT) (EXTR LYG0 NACT))) ELEM 1;
        FINSI;
        TABISO0. (indice + 1) = LXIN0;
        TABISO0. (indice + 2) = LYIN0;
        TABISO0. (indice + 3) = LXOUT0;
        TABISO0. (indice + 4) = LYOUT0;
        i = i + 1;
    FIN BOUCLE1;
FINSI;    

SI (NTEST > 0);
* Points root potentiels ordonnes par contrainte decroissante
    TABTEST = TABLE;
    i = 0;
    LSP11  = 'PROG' ;
    LPTEST = 'PROG' ;
    REPETER BOUCLE2 NTEST;
        i = i + 1;
        STEST = EXTR SP11 'SI11' 1 i 1;
        LSP11 = LSP11 ET STEST ;
        LPTEST = LPTEST ET i ;
    FIN BOUCLE2;
    TABTEST. 1 = LPTEST;
    TABTEST. 2 = LSP11;
    TABTEST = ORDO TABTEST 2 'DECR';
FINSI;

*
*\\\\ Resolution du probleme de conduction ////*
*
* Cosinus des tangentes
COS2X = EXCO PRIN1 'COX2';
COS2Y = EXCO PRIN1 'COY2';
COS2Z = EXCO PRIN1 'COZ2';

* Champ vectoriel des tangentes
*SIGMA22 = MANU 'CHML' MODEL11 'SI22' 1. 'TYPE' 
*'CONTRAINTES PRINCIPALES' 'STRESSES';
*PRIN12 = (SIGMA22 ET COS2X ET COS2Y ET COS2Z);
*VECTAN = VECT PRIN12 MODEL11 0.001 'SI22';    

* Modele thermique
MOTH = MODE SURFACE1 THERMIQUE ANISOTROPE;
* Perturbation conductivite numerique
TTPETIT = 1.e-6;
* Tenseur conductivite anisotrope
TTX = NOMC 'UX' COS2X;
TTY = NOMC 'UX' COS2Y;
KXX = (TTX ** 2) + ((TTX ** 0) * TTPETIT);
KYY = (TTY ** 2) + ((TTY ** 0) * TTPETIT);
KXY = TTX * TTY;
PX = 1.0 0.;
MATH = MATE MOTH 'DIRECTION' PX K11 KXX K21 KXY K22 KYY;
* Matrice de rigidite
KTH = COND MOTH MATH;
* Conditions aux limites
PCL1 = CHAN 'POI1' (TAB2. 'LBC'); 
CLT1 = BLOQ PCL1 'T';
CLTH = CLT1;
* Temperature imposee
VBC = TAB2. 'BCSTH';
DE1 = DEPI CLT1 VBC;
DETH = DE1;
* Assemblage
KTHTOT = KTH ET CLTH;
* Resolution
RESTH = RESO KTHTOT DETH;
T0 = CHAN 'CHAM' RESTH MOTH 'NOEUD';
T0 = EXCO T0 'T';
TAB2. 'RESTH'. NPAS = RESTH;   

*    
*\\\\ Extraction des isovaleurs ////*   
*
SI (NROOTS EGA MAXR);
    NTEST = 0;
FINSI;
TISOPAS = TABLE;
TABMESH1 = TABLE;
MESH1 = VIDE 'MAILLAGE';
***on choisit un point du domaine en faisant l'hypothese***
***que aucune fissure passe par ceci***
P0 = MANU 'POI1' ((X0 * Y0) POIN 'MINI');
i = 0;
REPETER BOUCLE3 (NROOTS + NTEST);          
    i = i + 1;
* Fissure activees aux pas precedents   
    SI (i <EG NROOTS);    
        PPA = TABROOT. i;
        MPPA = MANU 'POI1' PPA;
        T0P = PROI T0 MPPA; 
        TG = EXTR T0P 'T' PPA;
        ISOVG = ISOV T0 TG;
*        TRAC (ISOVG ET (CONT SURFACE1));
***cette phase est critique***
***si on perd l'unicite de la solution ORDO ne marche pas***
***arreter la procedure de tracking***
        ISOVG = ORDO ISOVG;         
                T0G = MANU 'CHML' MOTH 'T' TG 'NOEUD';
                T01 = T0 - T0G;
                T02 = T01 MASQ 'SUPERIEUR' 0.;
                MAILA = T02 ELEM 'EGALE' 1.;
                MAILB = T02 ELEM 'EGALE' 0.; 
                MESHISO = SURFACE1 DIFF (MAILA ET MAILB);
* Extraction du maillage  
                MM1 = (CHAN 'POI1' MESHISO) ET P0;
                MMT = (CHAN 'POI1' MESH1) ET P0;        
                CHECK = MMT INTE MM1;
                SI ((NBEL CHECK) > 1);
            ISOVG = TAB2. 'ISOTEST'. i;
                FINSI;
        MESH1 = MESH1 ET MESHISO;
        TISOPAS. i = ISOVG;
        TABMESH1. i = MESHISO;
        TAB2. 'ISOTEST'. i = ISOVG;       
* Nouvelles fissures
***on fait l'hypothese d'avoir une seule nouvelle fissure***
***par pas de temps***        
    SINON;
        iTEST = ENTI (EXTR (TABTEST. 1) (i - NROOTS));
        ETEST = ELSP11 ELEM iTEST;
        PPP = BARY ETEST;
        MPPP = MANU 'POI1' PPP;
        T0P = PROI T0 MPPP; 
        TG = EXTR T0P 'T' PPP;
        ISOVG = ISOV T0 TG;        
* Extraction du maillage 
                T0G = MANU 'CHML' MOTH 'T' TG 'NOEUD';
                T01 = T0 - T0G;
                T02 = T01 MASQ 'SUPERIEUR' 0.;
                MAILA = T02 ELEM 'EGALE' 1.;
                MAILB = T02 ELEM 'EGALE' 0.;
                MESHISO = SURFACE1 DIFF (MAILA ET MAILB);        
                MM1 = (CHAN 'POI1' MESHISO) ET P0;
                MMT = (CHAN 'POI1' MESH1) ET P0;        
                CHECK = MMT INTE MM1;        
* On verifie s'il y a des intersection entre maillages
***dans ce cas la on ecarte l'isovaleur***        
        SI ((NBEL CHECK) EGA 1);
            ISOVG = ORDO ISOVG;
*            TRAC (ISOVG ET (CONT SURFACE1));
            T0RED = REDU T0 ETEST;
            ISOROOT = ISOV T0RED TG;
            PROOT = ISOROOT POIN 1;
            CHECK1 = SURFACE1 ELEM 'CONTENANT' PROOT ;
***on fait l'hypothese que le root element appartient au contour***            
            SI ((NBEL CHECK1) > 1);
                PROOT = ISOROOT POIN 2;
            FINSI;
            TABROOT. (NROOTS + 1) = PROOT;
                        TABEROOT. (NROOTS + 1) = ETEST; 
            TISOPAS. (NROOTS + 1) = ISOVG;
            TAB2. 'ISOTEST'. (NROOTS + 1) = ISOVG;
            TABMESH1. (NROOTS + 1) = MESHISO;
            TAB2. 'NROOTS' = (TAB2. 'NROOTS') + 1;
            QUIT BOUCLE3;
                FINSI;
    FINSI;    
FIN BOUCLE3;    
TAB2. 'PROOTS' = TABROOT;
TAB2. 'EROOTS' = TABEROOT;

***cette partie sert a creer un systeme d'abscisses curvilignes***
***avec sense positif selon celui de propagation des fissures***
***l'origine est en pointe de fissure***
TABISO1 = TABLE;
NN = DIME TABROOT;
*LIST NN;
i = 0;
REPETER BOUCLE4 NN;
    i = i + 1;
    indice = (i - 1) * 4;
    PROOT = TABROOT. i;    
    EROOT = TABEROOT. i;
    ISOVG = TISOPAS. i;
    MESHISO = TABMESH1. i;
    NUMEL = (NBNO ISOVG) - 1;
    NACT = 0;
    SI ((VFLAG > 0.) ET (i <EG NROOTS));
        FLAGISO = REDU FLAG (TAB2. 'MESH_ISO'. i);
        MESHACT = FLAGISO ELEM 'SUPERIEUR' 0.;
        NACT = NBEL MESHACT;
    FINSI;
    XISOVG = COOR 1 ISOVG;
    YISOVG = COOR 2 ISOVG;
    EVOLX = EVOL 'CHPO' XISOVG ISOVG;
    EVOLY = EVOL 'CHPO' YISOVG ISOVG;
    LISX = EXTR EVOLX 'SCAL';
    LISY = EXTR EVOLY 'SCAL';
    LISABSC = EXTR EVOLX 'ABSC';
    NUMP = NBNO ISOVG;
    LISP = PROG 1 PAS 1 NUMP;    
    j = 0;
    REPETER B_INT NUMP;
        j = j + 1;
        Pj = ISOVG POIN j;
        SI ('<EG' ('DIST' Pj PROOT) 1.E-10);
            jroot = j;
*            LIST jroot;
*            LIST NUMP;
            QUIT B_INT;
        FINSI;
    FIN B_INT;  
    SROOT = EXTR LISABSC jroot;
    LISABSC1 = LISABSC - ((LISP / LISP) * SROOT);
    MM = MAXI LISABSC1;
    TABABS = TABLE;
    TABABS. 1 = LISABSC1;
    TABABS. 2 = LISP;
    TABABS. 3 = LISX;
    TABABS. 4 = LISY;
    SI (NACT EGA 0);
***si MM=0 il faut changer le signe***    
        SI ((MM EGA 0.));
            TABABS. 1 = -1. * LISABSC1;
            TABABS = ORDO TABABS 1 'CROI';
        FINSI;
***le dernier point est un point de sortie***        
        LXIN = ENLE (TABABS. 3) (DIME LISX);
        LYIN = ENLE (TABABS. 4) (DIME LISY);
***le premier point est un point d'entree***        
        LXOUT = ENLE (TABABS. 3) 1;
        LYOUT = ENLE (TABABS. 4) 1;
        TABISO1. (indice + 1) = LXIN;
        TABISO1. (indice + 2) = LYIN;
        TABISO1. (indice + 3) = LXOUT;
        TABISO1. (indice + 4) = LYOUT;        
    SINON;
        TABISO1. (indice + 1) = TABISO0. (indice + 1);
        TABISO1. (indice + 2) = TABISO0. (indice + 2); 
        TABISO1. (indice + 3) = TABISO0. (indice + 3);
        TABISO1. (indice + 4) = TABISO0. (indice + 4);   
        SI (NACT < NUMEL);
***on cherche le point correspondant a la pointe de fissure***        
            PTEST = ISOVG POIN (jroot + 1);   
            TEST = NBEL (EROOT ELEM 'CONTENANT' PTEST 'NOVERIF');
            SI (TEST > 0);   
                TABABS. 1 = -1. * LISABSC1;
                jroot = (NUMP - jroot) + 1; 
            FINSI;  
            TABABS = ORDO TABABS 1 'CROI';
            SI (jroot EGA 2);
                LISTEN = 1;
            SINON;
                LISTEN = ENTI (PROG 1 PAS 1 (jroot - 1));
            FINSI;
            LISX1 = ENLE (TABABS. 3) LISTEN;
            LISY1 = ENLE (TABABS. 4) LISTEN;
* Partie potentielle de la fissure       
            LXIN1 = ENLE LISX1 (DIME LISX1);
            LYIN1 = ENLE LISY1 (DIME LISY1);  
            LXOUT1 = ENLE LISX1 1;
            LYOUT1 = ENLE LISY1 1;
* On integre avec la partie activee         
            TABISO1. (indice + 1) = 
            TABISO0. (indice + 1) ET LXIN1; 
            TABISO1. (indice + 2) = 
            TABISO0. (indice + 2) ET LYIN1;
            TABISO1. (indice + 3) = 
            TABISO0. (indice + 3) ET LXOUT1;
            TABISO1. (indice + 4) = 
            TABISO0. (indice + 4) ET LYOUT1;
        FINSI;
    FINSI;        
FIN BOUCLE4;
TAB2. 'TABISO' = TABISO1;

*
*\\\\ Creation des chamelem ////*       
*
***CHM vaut sur les isovaleurs et 0 ailleurs***
CHM = SIGMA11 * 0.;
* Abscisses des points d'entree des fissures
CHIX = CHM;
* Ordonnees des points d'entree des fissures
CHIY = CHM;
* Composante x des normales
CHNX = CHM;
* Composante y des normales
CHNY = CHM;
TAB2. 'MESH'. NPAS = VIDE 'MAILLAGE';
i = 0;
REPETER BOUCLE5 NN;
    i = i + 1;
    indice = (i - 1) * 4;
    MESHISO = TABMESH1. i;
    NACT = 0;
    SI ((VFLAG > 0.) ET (i <EG NROOTS));
        FLAGISO = REDU FLAG (TAB2. 'MESH_ISO'. i);
        MESHACT = FLAGISO ELEM 'SUPERIEUR' 0.;
        NACT = NBEL MESHACT;
    FINSI;
    LXIN = TABISO1. (indice + 1);
    LYIN = TABISO1. (indice + 2);
    LXOUT = TABISO1. (indice + 3);
    LYOUT = TABISO1. (indice + 4);
    LXG = (LXIN + LXOUT) / 2.; 
    LYG = (LYIN + LYOUT) / 2.;
    SI (NACT EGA 0);
        MTEST = MESHISO;
    SINON;
        MDIFF = MESHISO DIFF MESHACT;
        MTEST = MDIFF ET MESHACT; 
    FINSI;
    j = 0;
    REPETER B_INT (DIME LXIN);
        j = j + 1;
* Coordonnees des points entree et sortie
        XIj = EXTR LXIN j;
        YIj = EXTR LYIN j;
        XOj = EXTR LXOUT j;
        YOj = EXTR LYOUT j;
* Coordonnees des points moyens   
        XGj = EXTR LXG j;
        YGj = EXTR LYG j;
        PGj = XGj YGj;    
* Calcul de la tangente et de la normale a la fissure   
        LL = (((XIj - XOj) ** 2) + ((YIj - YOj) ** 2)) ** 0.5;
        SI (YIj > YOj);
            TY = (YIj - YOj) / LL;
            TX = (XIj - XOj) / LL;
        SINON;
            TY = (YOj - YIj) / LL;
            TX = (XOj - XIj) / LL;
        FINSI;
        SI (YIj EGA YOj);
            TY = 0.;
            TX = -1.;
        FINSI;
        NX = TY;
        NY = -1. * TX;    
* Element contenant le point moyen
        ELj = (MTEST ELEM 'CONTENANT' PGj) ELEM 1;
        VALXj = MAXI (REDU X0S ELj);
        VALYj = MAXI (REDU Y0S ELj);
* Champ nul sauf dans ELi 
        CHPX = X0S MASQ 'EGALE' VALXj;
        CHPY = Y0S MASQ 'EGALE' VALYj;
        CHP = CHPX * CHPY;

* Champs par elements pour le model EFEM
        CHM = CHM + CHP;    
        CHIX = CHIX + (CHP * XIj);
        CHIY = CHIY + (CHP * YIj);
        CHNX = CHNX + (CHP * NX);
        CHNY = CHNY + (CHP * NY);

* Extraction de la mesh de la fissure potentielle
        SI(j EGA 1);
            LISO = DROI 1 (XIj YIj) (XOj YOj); 
            MESHi = ELj;
        SINON;
            LISO = LISO ET (DROI 1 (XIj YIj) (XOj YOj));
            MESHi = MESHi ET ELj;
        FINSI;        
    FIN B_INT; 
    TAB2. 'ISOTOT'. i = LISO;
    TAB2. 'MESH_ISO'. i = MESHi; 
    TAB2. 'MESH'. NPAS = (TAB2. 'MESH'. NPAS) ET MESHi;     
    
FIN BOUCLE5;
*TRAC CHM MODEL11;

*
*\\\\ Mise a jour des parametres materiaux ////*        
*
MATr2 = MATE MODEL11 YOUN E11 NU NU11 FT P11 XNX CHNX
XNY CHNY IND1 CHM XIE CHIX YIE CHIY;
MATr1 = MATr2;
TAB2. 'WTABLE'. 'CARACTERISTIQUES' = MATR1;


FINSI;


FINP MATR1 ;