Sources Esope | Cast3M

psiphi
C PSIPHI    SOURCE    CB215821  20/11/25    13:37:58     10792          
      SUBROUTINE PSIPHI
***********************************************************************
c OPERATEUR : PSIP
c
c APPEL :
c CHP1 (CHP2) (CHP3) = PSIP MAIL1 MAIL2 (CRIT1) (|'DEUX' | | P1 (P2) | );
c                                                |'TROI' | | MAIL3   |
c
c FONCTION : calcule les fonctions distances signées (level set)
c           relatives aux maillages MAIL2 (surface de fissure)
c           et MAIL3 en 3D (front de fissure)
c           ou P1 en 2D (pointe de fissure)
c           aux noeuds de MAIL1 (pas trop éloignés).
c
c CREATION : chat (16/07/2007)
c MODIFS   : bp (2009 -> 2012) : diverses corrections/ameliorations
c               repertoriees par les fiches d'anomalies/developpement
c            bp 14/03/2012 : on split psiphi en psiphi psip2d psip3d et
c                            zonag2 (evolué seulement le 18/12/2013)
c
c TO DO : cas 3d avec plusieurs fronts (ou 1 front discontinu)
c
***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMELEME
-INC SMCHPOI
-INC CCREEL
      PARAMETER(NMOCLE=2)
      CHARACTER*4 MOCLE(NMOCLE)
      DATA MOCLE/'DEUX','TROI'/
 
      if(iimpi.ge.1) write(ioimp,*)'==== ENTREE DANS PSIPHI ====='
      segact mcoord
 
***********************************************************************
*     INITIALISATIONs et LECTUREs
***********************************************************************
*
* calcule t'on phi? phi ET psi? phi, psi et tau?
      ideux=0
      call lirmot (mocle,nmocle,ideux,0)
      if (ideux.eq.2.and.idim.lt.3) then
        write(*,*) ' ON NE PEUT CALCULER TROIS LEVEL SET QU EN 3D !'
        write(*,*) ' ON CONTINUE AVEC LE CALCUL DE 2 LEVEL SET ...'
        ideux=1
      endif
 
* lecture des points pour lesquels on veut calculer phi et psi
      call LIROBJ('MAILLAGE',ipt1,1,iretou)
      melmai=ipt1
      if(ierr.ne.0) return
      call change(ipt1,1)
c       segact ipt1
c rem: inutile car change laisse ipt1 actif
 
* lecture du maillage de la fissure
      call LIROBJ('MAILLAGE',meleme,1,iretou)
      melfis=meleme
      if(ierr.ne.0) return
 
* lecture de la pointe de fissure (objet de type point) ou du front
      ip1=0
      ip2=0
      melfro=0
    1 continue
      if (ideux.ge.1) then
      if (idim.eq.2) then
        call LIROBJ('POINT',ipt,0,iretou)
        if (iretou.ne.0) then
          if (ip1.eq.0) then
             ip1 = ipt
             go to 1
          else
             ip2 = ipt
          endif
        endif
      else
        call LIROBJ('MAILLAGE', melfro ,1,iretou)
        if(iretou.eq.0) write(ioimp,*) 'Il manque le maillage du front'
        if(IERR.NE.0) return
      endif
      endif
 
* lecture facultative d'une longueur max (=plus grande taille des elements
* concernés par la fissure)
      xcrit=0.D0
      call lirree(xcrit,0,ircrit)
 
 
 
***********************************************************************
* CRÉATION DES MCHPOI DE SORTIE (1 POUR PSI , 2 POUR PHI, 3 POUR TAU)
***********************************************************************
      mpova1=0
      mpova2=0
      mpova3=0
 
*-----PHI-----
c       if(idebug.eq.1) write(6,*) '----Creation de PHI----'
      nat=1
      nsoupo=1
      nc=1
      n=ipt1.num(/2)
      segini,mchpo2
      segini,msoup2
      mchpo2.jattri(1)=1
      mchpo2.ipchp(1)=msoup2
      mchpo2.ifopoi=ifour
      segdes mchpo2
      msoup2.igeoc=ipt1
      msoup2.nocomp(1)='PHI'
      msoup2.noharm(1)=nifour
      segini,mpova2
      msoup2.ipoval=mpova2
c       segdes,msoup2
 
*-----PSI-----
      if (ideux.ge.1) then
c       if(idebug.eq.1) write(6,*) '----Creation de PSI----'
        segini,mchpo1
        segini,msoup1
        mchpo1.jattri(1)=1
        mchpo1.ipchp(1)=msoup1
        mchpo1.ifopoi=ifour
        segdes,mchpo1
        msoup1.igeoc=ipt1
        msoup1.nocomp(1)='PSI'
        msoup1.noharm(1)=nifour
        segini,mpova1
        msoup1.ipoval=mpova1
c         segdes,msoup1
      else
        msoup1 = 0
      endif
 
*-----TAU-----
      if (ideux.ge.2) then
c       if(idebug.eq.1) write(6,*) '----Creation de TAU----'
        segini,mchpo3
        segini,msoup3
        mchpo3.jattri(1)=1
        mchpo3.ipchp(1)=msoup3
        mchpo3.ifopoi=ifour
        segdes,mchpo3
        msoup3.igeoc=ipt1
        msoup3.nocomp(1)='TAU'
        msoup3.noharm(1)=nifour
        segini,mpova3
        msoup3.ipoval=mpova3
c         segdes,msoup3
      else
        msoup3 = 0
      endif
 
 
***********************************************************************
*     PSIPHI 3D/2D
***********************************************************************
      if (idim.eq.3) then
c       write(ioimp,*)'appel PSIP3D (',ideux,ipt1,melfis,melfro,xcrit
        call PSIP3D(ideux,ipt1,melfis,melfro,xcrit,
     &              msoup1,msoup2,msoup3)
      else
c       write(ioimp,*)'appel PSIP3D (',ideux,ipt1,melfis,ip1,ip2,xcrit
        call PSIP2D(ideux,ipt1,melfis,ip1,ip2,xcrit,
     &              msoup1,msoup2)
      endif
 
 
***********************************************************************
*     ECRITURE ET FIN DU PROGRAMME
***********************************************************************
 
c-----ecriture de (PSI) PHI ((TAU)) --------------
      if(ideux.ge.2) then
        call ACTOBJ('CHPOINT ',mchpo3,1)
        call ECROBJ('CHPOINT ',mchpo3)
      endif
      call ACTOBJ('CHPOINT ',mchpo2,1)
      call ECROBJ('CHPOINT ',mchpo2)
      if(ideux.ge.1) then
        call ACTOBJ('CHPOINT',mchpo1,1)
        call ECROBJ('CHPOINT',mchpo1)
      endif
 
      end