C ELEMEL    SOURCE    GOUNAND   25/08/04    21:15:06     12340          
      SUBROUTINE ELEMEL(IPT1,IPT2,MELEME,nbltot)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : ELEMEL
C DESCRIPTION : Extrait les elements de IPT1 appuyes sur les
C               elements de IPT2
C               Un element de IPT1 est appuye sur un element de IPT2
C               s'il contient tous les noeuds de ce dernier
C
C
C LANGAGE     : ESOPE
C AUTEUR      : Stephane GOUNAND (CEA/DES/ISAS/DM2S/SEMT/LTA)
C               mel : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          :
C APPELES (E/S)    :
C APPELES (BLAS)   :
C APPELES (CALCUL) :
C APPELE PAR       :
C***********************************************************************
C SYNTAXE GIBIANE    :
C ENTREES            :
C ENTREES/SORTIES    :
C SORTIES            :
C***********************************************************************
C VERSION    : v1, 01/08/2025, version initiale
C HISTORIQUE : v1, 01/08/2025, creation
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC SMELEME
-INC SMCOORD
      segment icpr(nbpts)
      segment ipos(nno2+1)
      segment izone(nlel)
      segment inum(nlel)
      segment linum(imaxel)
      segment lizone(imaxel)
      segment imail(nbsou1)
      logical lfound,ldbg
*
* Executable statements
*
      ldbg=.false.
*
*   Numerotation locale a ipt2
*
      call actobj(ipt2,'MAILLAGE',1)
      segini icpr
      nno2=0
      ipt4=ipt2
      DO IOB=1,MAX(1,IPT2.LISOUS(/1))
         IF (IPT2.LISOUS(/1).NE.0) IPT4=IPT2.LISOUS(IOB)
         if (ldbg) write(ioimp,*) 'IOB=',IOB
         if (ldbg) segprt,IPT4
         DO J=1,IPT4.NUM(/2)
            DO I=1,IPT4.NUM(/1)
               IKI=IPT4.NUM(I,J)
               IF (ICPR(IKI).EQ.0) THEN
                  nno2=nno2+1
                  ICPR(IKI)=nno2
               ENDIF
            ENDDO
         ENDDO
      ENDDO
*     Recopions les sous-maillages de ipt1 dans imail parce qu'on va les
*     tagger, voir les reduire apres
      call actobj(ipt1,'MAILLAGE',1)
      nbsou1=MAX(1,IPT1.LISOUS(/1))
      segini imail
      ipt3=ipt1
      DO isous=1,nbsou1
         IF (IPT1.LISOUS(/1).NE.0) IPT3=IPT1.LISOUS(isous)
         segini,ipt5=ipt3
         imail(isous)=ipt5
      ENDDO
      if (ldbg) write(ioimp,*) 'nno2,nbsou1=',nno2,nbsou1
*
* Correspondance chaque noeud de ipt2, les elements de ipt1 le contenant
*
      segini ipos
      DO isous=1,nbsou1
         ipt3=imail(isous)
         DO J=1,ipt3.NUM(/2)
            DO I=1,ipt3.NUM(/1)
               ikil=ICPR(ipt3.NUM(I,J))
               if (ikil.ne.0) ipos(ikil)=ipos(ikil)+1
            ENDDO
         ENDDO
      ENDDO
      imaxel=ipos(1)
      i_z = ipos(1)
      DO 13 i = 2, nno2
         imaxel=max(imaxel,ipos(i))
         i_z = i_z + ipos(i)
         ipos(i) = i_z
 13   CONTINUE
      nlel = ipos(nno2)
      ipos(nno2+1) = nlel
      segini,izone
      segini,inum
      DO isous=1,nbsou1
         ipt3=imail(isous)
         DO J=1,ipt3.NUM(/2)
            DO I=1,ipt3.NUM(/1)
               ikil=ICPR(ipt3.NUM(I,J))
               if (ikil.ne.0) then
                  ide=ipos(ikil)
                  izone(ide)=isous
                  inum(ide)=j
                  ipos(ikil)=ide-1
               endif
            ENDDO
         ENDDO
      ENDDO
c      segprt,ipos
c      write(ioimp,*) 'imaxel=',imaxel
c      segprt,izone
c      segprt,inum
*
*     Pour chaque element de ipt2, on regarde les elements de ipt1 qui
*     le contiennent. On tagge les elements de IPT1 en changeant le
*     signe du 1er noeud.
*
      segini,lizone
      segini,linum
      call actobj(ipt2,'MAILLAGE',1)
*      segprt,ipt2
      ipt4=ipt2
      DO IOB=1,MAX(1,IPT2.LISOUS(/1))
         IF (IPT2.LISOUS(/1).NE.0) IPT4=IPT2.LISOUS(IOB)
         DO 30 IL2=1,IPT4.NUM(/2)
            if (ldbg) write(ioimp,*) 'ipt2 isou element ',iob,il2
* On cherche le noeud de IPT2 avec le plus petit nombre de voisins
* dans IPT1, on saute l'element si un des noeuds n'est pas dans IPT1
            NVOIS=IGRAND
            IPMIN=0
            DO IP2=1,IPT4.NUM(/1)
               ik=ICPR(IPT4.NUM(IP2,IL2))
               nik=ipos(ik+1)-ipos(ik)
               if (ldbg) write(ioimp,*) '  ipt2 noeud ',ip2,IPT4.NUM(IP2
     $              ,IL2),nik
               if (nik.lt.nvois) then
                  nvois=nik
                  ipmin=ip2
               endif
            ENDDO
            if (ipmin.eq.0) goto 30
            ik=ICPR(IPT4.NUM(IPMIN,IL2))
            ideb=ipos(ik)+1
            ifin=ipos(ik+1)
            nik=ifin-ideb+1
            do iik=1,nik
               lizone(iik)=izone(ideb+iik-1)
               linum(iik) =inum(ideb+iik-1)
            enddo
            if (ldbg) write(ioimp,*) '   ipt2 ipmin,nik=',ipmin,nik
            do 300 iik=1,nik
               iz=lizone(iik)
               in=linum(iik)
               if (ldbg) write(ioimp,*) '   iik,iz,in=',iik,iz,in
               DO IP2=1,IPT4.NUM(/1)
                  if (ldbg) then
                     write(ioimp,*) '   ip2,ipmin=',ip2,ipmin
                     write(ioimp,*) '     lizone :',(lizone(ibb),ibb=1
     $                    ,nik)
                     write(ioimp,*) '     linum  :',(linum(ibb),ibb=1
     $                    ,nik)
                  endif
                  if (ip2.ne.ipmin) then
                     lfound=.false.
                     ik=ICPR(IPT4.NUM(IP2,IL2))
                     ideb=ipos(ik)+1
                     ifin=ipos(ik+1)
                     do ie=ideb,ifin
                        iz2=izone(ie)
                        in2=inum(ie)
                        if (iz2.eq.iz.and.in2.eq.in) then
                           lfound=.true.
                           goto 3000
                        endif
                     enddo
 3000                continue
                     if (ldbg) write(ioimp,*) '   iik,iz,in,lfound=',iik
     $                    ,iz,in,lfound
                     if (.not.lfound) then
                        lizone(iik)=0
                        linum(iik)=0
                        goto 300
                     endif
                  endif
               ENDDO
 300        continue
            do iik=1,nik
               iz=lizone(iik)
               in=linum(iik)
               if (iz.ne.0.and.in.ne.0) then
                  ipt3=imail(iz)
                  ipt3.num(1,in)=-abs(ipt3.num(1,in))
               endif
            enddo
 30      CONTINUE
      ENDDO
      segsup lizone
      segsup linum
      segsup izone
      segsup inum
      segsup ipos
      segsup icpr
*
* On reduit de maniere adequate les maillages de imail
*
      nsou2=nbsou1
      nbltot=0
c      segprt,imail
      do isous=1,nbsou1
         ipt3=imail(isous)
         nel=ipt3.num(/2)
         nno=ipt3.num(/1)
         nbref=ipt3.lisref(/1)
         iel2=0
         do iel=1,nel
            if (ipt3.num(1,iel).lt.0) then
               iel2=iel2+1
               ipt3.num(1,iel2)=-ipt3.num(1,iel)
               do ino=2,nno
                  ipt3.num(ino,iel2)=ipt3.num(ino,iel)
               enddo
            endif
         enddo
         if (iel2.eq.0) then
            segsup ipt3
            imail(isous)=0
            nsou2=nsou2-1
         elseif (iel2.le.nel) then
            nbnn=nno
            nbelem=iel2
            nbsous=0
            segadj ipt3
         endif
         nbltot=nbltot+iel2
      enddo
c      write(ioimp,*) 'nsou2=',nsou2
      if (nsou2.eq.0) then
         call melvid(ilcour, meleme)
      elseif (nsou2.eq.1) then
         do isous=1,nbsou1
            if (imail(isous).ne.0) then
               meleme=imail(isous)
               goto 45
            endif
         enddo
 45      continue
      else
         nbnn=0
         nbelem=0
         nbsous=nsou2
         nbref=0
         segini meleme
         isou2=0
         do isous=1,nbsou1
            if (imail(isous).ne.0) then
               isou2=isou2+1
               lisous(isou2)=imail(isous)
               if (isou2.eq.nsou2) goto 55
            endif
         enddo
 55      continue
      endif
      segsup imail
*
* Normal termination
*
      RETURN
*
* End of subroutine ELEMEL
*
      END
 
