C CMCT3B    SOURCE    PV090527  26/04/30    21:15:18     12529          
      SUBROUTINE CMCT3B(DESCOB,LSINCO,LTINCO,LSINCD,IRIG2)
*_______________________________________________________________________
c
c      opérateur cmct
c
c entrée
c   ICHP  : champ par point qui stocke la masse inversée M-1
c   IRIGB : rigidité B
c   IRIGC : rigidité C
c
c sortie
c   IRIG2 : rigidité contenant la matrice condensée C M-1 Bt
c
*_______________________________________________________________________

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
*

-INC PPARAM
-INC CCOPTIO
-INC SMRIGID
-INC SMELEME
-INC CCHAMP
-INC SMCOORD
-INC SMCHPOI
*
* Description des objets à traiter
* PNTOB : pointeur de l'objet
* TYPOB : type de l'objet
* PNTCOB : pointeur vers une version locale de l'objet
*          pour la RIGIDITE IMAT
*                                PNTCOB(1,IMAT) pointe vers un CORES1
*                                PNTCOB(2,IMAT) pointe vers un LSINCO
*                                PNTCOB(3,IMAT) pointe vers un MCOEF

      SEGMENT DESCOB
        INTEGER PNTOB(NMAT)
        INTEGER PNTCOB(3,NMAT)
      ENDSEGMENT
*  stockage des noms des composantes duales.
      SEGMENT LSINCD
        CHARACTER*(LOCOMP) LISIND(NLIGD)
      ENDSEGMENT
*
* Noms d'inconnues primales correspondant à LISIND
*
      SEGMENT CORIDP
        CHARACTER*(LOCOMP) LISIDP(NLIGD)
      ENDSEGMENT
*
* tableau pour pointer vers MCOEF à partir du nombre d'inconnues
*
      SEGMENT LSINCO
        INTEGER LESINC(NINC+1,2,NMAT)
      ENDSEGMENT
      SEGMENT LTINCO
        REAL*8  XMAS(NINC)
      ENDSEGMENT
*
*  tableau des coefficient de la matrice C
*  ordonné dans l'ordre des inconnues
      SEGMENT MCOEF
*        numero du noeud support du multiplicateur ligne 1
      INTEGER ICOEF(2,NCOEF)
*         valeur des coefficients
        REAL*8 XCOEF(NCOEF)
      ENDSEGMENT
*
      SEGMENT WORK1
        REAL*8 XDUM(NBNN)
      ENDSEGMENT
*
      segment inddl(ninc)
      segment icat(0)
      segment ivcat(jg)
      segment incat(ncat)
      segment jcat(ncat+1)
      segment kcat(ninc)
      segment khash(ninc)
* dans ce segment, valeurs positives => iinc uniques dans knuniq
*     valeurs negatives => iinc en doublon valeur absolue pointe sur le
*     unique dans knuniq
      segment kuniq(ninc)
* nombre de iinc uniques dans chaque categorie
      segment inuniq(ncat)
* stokage des iinc
      segment knuniq(nuniq)
* nombre de iinc partageant un meme unique
      segment nnuniq(nuniq)
* un segment de correspondance
      segment jnuniq(nuniq)

      LOGICAL LDBG,LFOUND
*
      LDBG=.FALSE.
*_______________________________________________________________________
*
*   il ne reste plus qu' a creer les matrices élémentaires
*
* Construction des noms d'inconnues primales correspondant à LISIND
      NLIGD=LISIND(/2)
      SEGINI CORIDP
      DO ILIGD=1,NLIGD
         CALL PLACE(NOMDU,LNOMDU,idx,LISIND(ILIGD))
         IF (idx.EQ.0) THEN
            LISIDP(ILIGD)=LISIND(ILIGD)
         ELSE
            LISIDP(ILIGD)=NOMDD(idx)
         ENDIF
      ENDDO
*
*     Pour diminuer le nombre de NRIGEL a creer, on va chercher les
*     redondances dans le segment LESINC(I,1,IMAT), d'abord par nombre
*     de ddls (LESINC(I,2,IMAT)) puis par egalite des suites de ddls
*
      NINC=LESINC(/1)-1
      NMAT=LESINC(/3)
      MAXDDL=0
      segini inddl
      do IINC=1,NINC
         do IMAT=NMAT,1,-1
            inddl(iinc)=inddl(iinc)+LESINC(IINC,2,IMAT)
         enddo
         MAXDDL=max(maxddl,inddl(iinc))
      enddo
*
      segini icat
      jg=maxddl
      segini ivcat
      ncat=maxddl
      segini incat
      do IINC=1,NINC
         NDDL=inddl(iinc)
         if (incat(nddl).eq.0) then
            icat(**)=nddl
            ivcat(nddl)=icat(/1)
         endif
         incat(nddl)=incat(nddl)+1
      enddo
      IF (LDBG) THEN
         write(ioimp,*) 'Classement par categorie (nombre de ddls)'
         write(ioimp,*) '  Categories'
         segprt,icat
         segprt,ivcat
         write(ioimp,*) '  Nombre par categorie'
         segprt,incat
      ENDIF
      ncat=icat(/1)
      segini jcat
      jcat(1)=1
      do ic=1,ncat
         jcat(ic+1)=jcat(ic)+incat(icat(ic))
         incat(icat(ic))=0
      enddo
      segsup icat
      segadj incat
* verif
      if (jcat(ncat+1)-1.ne.ninc) then
         write(ioimp,*) 'ncat,ninc=',ncat,ninc
         segprt,jcat
         call erreur(5)
         return
      endif
      segini kcat
      do IINC=1,NINC
         nucat=ivcat(inddl(iinc))
         idx=jcat(nucat)+incat(nucat)
         kcat(idx)=iinc
         incat(nucat)=incat(nucat)+1
      enddo
      segsup incat
      segsup ivcat
      segsup inddl
*
* Classement par categorie (nombre de ddls )
*
      IF (LDBG) THEN
         write(ioimp,*) '  Classement des inconnues duales'
         segprt,jcat
         segprt,kcat
      ENDIF
* Un hash pour trier plus vite les non-identiques
      segini khash
      do IINC=1,NINC
         IHASH=0
         DO IMAT=NMAT,1,-1
            MCOEF=PNTCOB(3,IMAT)
            IDEB=LESINC(IINC,1,IMAT)
            IFIN=IDEB+LESINC(IINC,2,IMAT)-1
            DO IDX=IDEB,IFIN
               IHASH=IHASH*(NLIGD+1)+ICOEF(2,IDX)
            ENDDO
         ENDDO
         khash(IINC)=IHASH
      enddo
      IF (LDBG) THEN
         write(ioimp,*) 'Hash pour les inconnues duales'
         segprt,khash
      ENDIF
      segini kuniq
      segini inuniq
      nuniq=ninc
      segini knuniq
      segini nnuniq
      do ic=1,ncat
         idxdeb=jcat(ic)
         idxfin=jcat(ic+1)-1
         do j=idxdeb,idxfin
* Le premier de la passe est forcement unique
            iinc=kcat(j)
            lfound=.false.
            if (j.ne.idxdeb) then
               ihash=khash(iinc)
               do 32 k=1,inuniq(ic)
                  juniq=idxdeb+k-1
                  iinck=knuniq(juniq)
                  ihashk=khash(iinck)
                  if (ihashk.ne.ihash) goto 32
                  do imat=nmat,1,-1
                     mcoef=pntcob(3,imat)
                     kdeb=lesinc(iinc,1,imat)
                     knb =lesinc(iinc,2,imat)
                     kdebk=lesinc(iinck,1,imat)
                     do 33 inb=1,knb
                        icoefu=icoef(2,kdeb+inb-1)
                        icoefk=icoef(2,kdebk+inb-1)
                        if (icoefu.ne.icoefk) goto 32
 33                  continue
                  enddo
                  lfound=.true.
                  goto 34
 32            continue
 34            continue
            endif
            if (.not.lfound) then
               juniq=idxdeb+inuniq(ic)
               inuniq(ic)=inuniq(ic)+1
               knuniq(juniq)=iinc
               nnuniq(juniq)=1
               kuniq(j)=juniq
*               kuniq(iinc)=juniq
            else
               kuniq(j)=-juniq
*               kuniq(iinc)=-juniq
               nnuniq(juniq)=nnuniq(juniq)+1
            endif
         enddo
      enddo
      segsup khash
      IF (LDBG) THEN
         write(ioimp,*) 'Primales unique par duale'
         write(ioimp,*) '  Nombre de groupes de primales/duale'
         segprt,inuniq
         write(ioimp,*) '  Stockage des duales uniques'
         segprt,knuniq
         write(ioimp,*) '  Combien partagent un meme uniq'
         segprt,nnuniq
         write(ioimp,*) '  Pointeur pour chaque duale vers l''unique'
         segprt,kuniq
      ENDIF
      NRIGEL=0
      do i=1,ncat
         nrigel=nrigel+inuniq(i)
      enddo
      IF (LDBG) THEN
         write(ioimp,2030) 'Nombre de sous-zones =',NRIGEL
      ENDIF
*
*     Compacter knuniq, nnuniq et mettre à jour kuniq
*     juniq varie au debut entre 1 et ninc puis apres entre 1 et nrigel
*     On utilise jnuniq comme tableau de conversion juniq_avant juniq_apres
*
*      if (.false.) then
      if (nrigel.lt.ninc) then
         segini jnuniq
         juniqp=0
         do juniqv=1,ninc
            if (nnuniq(juniqv).ne.0) then
               juniqp=juniqp+1
               nnuniq(juniqp)=nnuniq(juniqv)
               knuniq(juniqp)=knuniq(juniqv)
               jnuniq(juniqv)=juniqp
            endif
         enddo
         do iinc=1,ninc
            if (kuniq(iinc).ge.0) then
               kuniq(iinc)=jnuniq(kuniq(iinc))
            else
               ixxx=-kuniq(iinc)
               kuniq(iinc)=-jnuniq(ixxx)
            endif
         enddo
         if (juniqp.ne.nrigel) then
            write(ioimp,*) 'juniqp.ne.nrigel',juniqp,nrigel
            call erreur(5)
            return
         endif
         nuniq=nrigel
         segadj,knuniq
         segadj,nnuniq
         segsup jnuniq
         IF (LDBG) THEN
            write(ioimp,*) 'Apres compactage'
            write(ioimp,*) '  Stockage des duales uniques'
            segprt,knuniq
            write(ioimp,*) '  Combien partagent un meme uniq'
            segprt,nnuniq
            write(ioimp,*) '  Pointeur pour chaque duale vers l''unique'
            segprt,kuniq
         ENDIF
      endif
*
      SEGINI MRIGID
      IRIG2 = MRIGID
      MTYMAT = 'CMCT    '
      IFORIG = IFOUR
*     Pour le mode de calcul on devrait s'appuyer sur celui des matrices/champs en entree ?
*      VRAI ! A faire

*
*    boucle sur les categories
*
      NBNN=maxddl
      SEGINI WORK1
*!      segini jnuniq
      iricnt=0
 2030 FORMAT (A,20(2X,I5))

      do ic=1,ncat
*!         idxdeb=jcat(ic)
         nrigun=inuniq(ic)
         do irigun=1,nrigun
            iricnt=iricnt+1
            coerig(iricnt)=1.D0
*!            juniq=idxdeb+irigun-1
            juniq=iricnt
            nbl=nnuniq(juniq)
            iinc=knuniq(juniq)
            if (LDBG) THEN
*!               write(ioimp,2030)
*!     $              'ic,idxdeb,irigun,iricnt,juniq,nbl,iinc=',ic
*!     $              ,idxdeb,irigun,iricnt,juniq,nbl,iinc
               write(ioimp,2030)
     $              'ic,irigun,iricnt,juniq,nbl,iinc=',ic
     $              ,irigun,iricnt,juniq,nbl,iinc
            endif
            NBNN=0
            IDXDUA=0
*
* S'il y a deux rigidités à multiplier, on stocke la primale
* puis la duale dans le MELEME et dans XDUM. A ce moment-là,
* IDXDUA+1 est l'index de départ des noeuds et des valeurs de la duale.
* On parcourt donc DESCOB à l'envers car on avait d'abord mis
* la duale C avant la primale B
*
            DO IMAT=NMAT,1,-1
               NBNN=NBNN+LESINC(iinc,2,IMAT)
               IF (IDXDUA.EQ.0.AND.NMAT.EQ.2) IDXDUA=NBNN
            ENDDO
* creation du maillage et du vecteur des coefficients
            NBELEM = nbl
            NBSOUS = 0
            NBREF = 0
            SEGINI MELEME
            ITYPEL = 32
            IRIGEL(1,iricnt) = MELEME
*
*    segment descripteur DESCR
*
            NLIGRP=NBNN
            IF (NMAT.EQ.1) THEN
               NLIGRP = NBNN
               NLIGRD = NBNN
            ELSE
               NLIGRP = IDXDUA
               NLIGRD = NBNN-IDXDUA
            ENDIF
            SEGINI DESCR
            IRIGEL(3,iricnt) = DESCR
*
            NELRIG =nbl
            rigrel=0
            SEGINI XMATRI
            IRIGEL(4,iricnt)=XMATRI
*     S'il y a deux matrices en entrée, le résultat
*     n'est sans doute pas symétrique
            IF (NMAT.EQ.2) then
               IRIGEL(7,iricnt)=2
               xmatri.symre=2
            endif
*! On recycle nnuniq pour stocker iricnt
*!            nnuniq(juniq)=iricnt
* On recycle nnuniq comme compteur
            nnuniq(juniq)=0
         enddo
*
*     Remplissage des objets
*
         idxdeb=jcat(ic)
         idxfin=jcat(ic+1)-1
         do jc=idxdeb,idxfin
            iinc=kcat(jc)
*     lfirst=kuniq(iinc).gt.0
            juniq=abs(kuniq(jc))
            iricou=juniq
*!            iricou=nnuniq(juniq)
            meleme=irigel(1,iricou)
            descr=irigel(3,iricou)
            xmatri=irigel(4,iricou)
            iincu=knuniq(juniq)
*!            jnuniq(juniq)=jnuniq(juniq)+1
*!            ibl=jnuniq(juniq)
            nnuniq(juniq)=nnuniq(juniq)+1
            ibl=nnuniq(juniq)
            if (LDBG) THEN
               write(ioimp,2030)
     $              'ic,jc,iinc,juniq,iincu,ibl,iricou='
     $              ,ic,jc,iinc,juniq,iincu,ibl,iricou
            endif
*     remplissage MELEME
            INOEU = 0
            DO IMAT=NMAT,1,-1
               MCOEF=PNTCOB(3,IMAT)
               DO 1200 J=0,LESINC(iinc,2,IMAT)-1
                  INOEU = INOEU + 1
                  NUM(INOEU,ibl) = ICOEF(1,J+LESINC(iinc,1,IMAT))
                  XDUM(INOEU) = XCOEF(J+LESINC(iinc,1,IMAT))
 1200          CONTINUE
            ENDDO

*     remplissage DESCR (une seule fois, on a tout fait pour)
            if (ibl.eq.1) then
               IPASS=0
               DO IMAT=NMAT,1,-1
                  MCOEF=PNTCOB(3,IMAT)
                  IPASS=IPASS+1
                  INOEU=0
                  DO 1300 J=0,LESINC(iincu,2,IMAT)-1
                     INOEU = INOEU + 1
                     IF (IPASS.EQ.1) THEN
                        LISINC(INOEU) = LISIDP(ICOEF(2,J+LESINC(iincu,1
     $                       ,IMAT)))
                        NOELEP(INOEU)=INOEU
                        IF (NMAT.EQ.1) THEN
                           LISDUA(INOEU) = LISIND(ICOEF(2,J+LESINC(iincu
     $                          ,1,IMAT)))
                           NOELED(INOEU)=INOEU
                        ENDIF
                     ENDIF
                     IF (IPASS.EQ.2) THEN
                        LISDUA(INOEU)
     $                       = LISIND(ICOEF(2,J+LESINC(iincu,1,IMAT)))
                        NOELED(INOEU)=IDXDUA+INOEU
                     ENDIF
 1300             CONTINUE
               ENDDO
            endif
*     remplissage XMATRI
            IF (LTINCO.NE.0) THEN
               XCOEF=XMAS(iinc)
            ELSE
               XCOEF=1.D0
            ENDIF
            nligrd=re(/1)
            nligrp=re(/2)
            if (nmat.eq.2) then
               idxdua=nligrp
            else
               idxdua=0
            endif
            DO 600 J=1,NLIGRP
               DO 500 K=1,NLIGRD
                  RE(K,J,ibl)=XDUM(IDXDUA+K)*XDUM(J)*XCOEF
 500           CONTINUE
 600        CONTINUE
* recyclage XDUM
            INOEU = 0
            DO IMAT=NMAT,1,-1
               MCOEF=PNTCOB(3,IMAT)
               DO 2200 J=0,LESINC(iinc,2,IMAT)-1
                  INOEU = INOEU + 1
                  XDUM(INOEU) = 0.D0
 2200          CONTINUE
            ENDDO
         enddo
      enddo
*!      segsup jnuniq
      segsup nnuniq
      segsup knuniq
      segsup inuniq
      segsup kuniq
      segsup jcat
      segsup kcat
      SEGSUP WORK1
      SEGSUP CORIDP
*
*_______________________________________________________________________
      RETURN
      END