C IMPOFU    SOURCE    PV090527  26/04/30    21:15:41     12529          
*  reunion des relations portant sur le meme multiplicateur de lagrange
*  si mchpoi est nul, on se contente de nettoyer les termes petits
*
      subroutine impofu(mrigid,mchpoi)
      implicit real*8 (a-h,o-z)
-INC SMRIGID
-INC SMCHPOI
-INC SMELEME

-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC CCGEOME
-INC CCREEL
*  nombre max d'elements par paquet
      parameter (nblim=16)
*
      segment icpr(nbpts)
      segment val((nbpo-1)*idim,nblag)
      segment ielem(nbpo,nblag)
      segment nbnl(nblag)
      segment dnorm(nblag)
      segini icpr
**    call prrigi(mrigid,0)
**    call ecchpo(mchpoi,0)
      segact mrigid
      nblag=0
*  regroupement des raideurs dans un grand tableau dimensionne au max
* on commence par reperer les multiplicateurs de lagrange et leur destination
      do 10 i=1,irigel(/2)
       meleme=irigel(1,i)
       segact meleme
       do 20 iel=1,num(/2)
        lag=num(1,iel)
        if (icpr(lag).eq.0) then
         nblag=nblag+1
         icpr(lag)=nblag
*  noter si formulation faible
         if (icolor(iel).eq.1) icpr(lag)=-nblag
        endif
  20   continue
  10  continue
* sauvegarde d'un descr pour avoir les noms d'inconnues
      des1=irigel(3,1)
      segact des1
* creation melme et xmatri a la taille max
      nbpo=7
      segini val,ielem,nbnl,dnorm
* remplissage, d'abort le mult
      do i=1,icpr(/1)
       if (icpr(i).ne.0) then
        ielem(1,abs(icpr(i)))=i
       endif
      enddo
C
      if (mchpoi.ne.0) then
        segact mchpoi*mod
        msoup1=ipchp(1)
        segact msoup1
        mpova1=msoup1.ipoval
        segact mpova1
      endif
C
      do 100 i=1,irigel(/2)
       meleme=irigel(1,i)
       xmatri=irigel(4,i)
       segact xmatri,meleme
       imod=0
       do 110 iel=1,num(/2)
        lag=num(1,iel)
        kel=abs(icpr(lag))
        remax=-1.
        do iv=2,re(/1)
          remax=max(remax,abs(re(iv,1,iel)))
        enddo
* recherche du point significatif le plus haut
        isig=0
        do 115 in=2,num(/1)
         ip=num(in,iel)
*  la contribution n'est pas significative on saute le noeud
         do iv=1,idim
           if (abs(re(1+(in-2)*idim+iv,1,iel)).gt.1d-10*remax) goto 116
         enddo
*  pour assurer que la somme des termes est nulle, on corrige le point isig
         if (imod.eq.0) segact xmatri*mod
         imod=1
         isign=isig
         if (isig.eq.0) isign=in+1
         if (isign.gt.num(/1)) call erreur(5)
         do iv=1,idim
           re(1+(isign-2)*idim+iv,1,iel)=
     >        re(1+(isign-2)*idim+iv,1,iel)+
     >        re(1+(in-2)*idim+iv,1,iel)
              re(1+(in-2)*idim+iv,1,iel)=0.d0
         enddo
         if (iimpi.ne.0) write (6,*) ' noeud elimine dans relation',
     >         in,ip,re(1+(in-2)*idim+1 ,1,iel),remax
         goto 115
 116     continue
         if (isig.eq.0) isig=in
 115    continue
        do 120 in=2,num(/1)
         ip=num(in,iel)
*  la contribution n'est pas significative on saute le noeud
         do iv=1,idim
           if (abs(re(1+(in-2)*idim+iv,1,iel)).gt.1d-10*remax) goto 160
         enddo
         goto 120
 160     continue
         do 130 ir=2,nbpo
          if (ielem(ir,kel).eq.0) goto 150
          if (ielem(ir,kel).ne.ip) goto 130
*  le noeud est deja dans l'element, on ajoute les valeurs
          do iv=1,idim
           val((ir-2)*idim+iv,kel)=val((ir-2)*idim+iv,
     >        kel)+re(1+(in-2)*idim+iv,1,iel)
          enddo
          goto 120
 130     continue
*  le noeud n'est pas dans l'element, on ajoute le noeud et les valeurs
 150    continue
        if (ir.gt.nbpo) then
         nbpo=ir
         segadj ielem,val
        endif
        ielem(ir,kel)=ip
        do iv=1,idim
         val((ir-2)*idim+iv,kel)=re(1+(in-2)*idim+iv,1,iel)
        enddo
 120    continue
C
        if (mchpoi.ne.0) then
          dnorm(kel) = dnorm(kel) + mpova1.vpocha(iel,2)
        endif
C
 110   continue
      segsup meleme,xmatri
      descr=irigel(3,i)
      if (i.ne.1) segsup descr
 100  continue
      segsup mrigid
*
*  eclatement en paquets de meme nb de noeuds
*  et renormalisation
*

*  nb noeuds par element
      do iel=1,ielem(/2)
        do in=1,ielem(/1)
         if (ielem(in,iel).eq.0) goto 200
        enddo
 200    continue
        nbnl(iel)=in-1
      enddo
*  elimination elem en double si point ligne
      if (idim.eq.3) then
      do 300 iel=1,ielem(/2)
       if (nbnl(iel).ne.4) goto 300
        do 301 jel=iel+1,ielem(/2)
        if (nbnl(jel).ne.4) goto 301
        if (ielem(4,iel).ne.ielem(4,jel)) goto 301
         if (ielem(2,iel)*ielem(3,iel).ne.ielem(2,jel)*ielem(3,jel)) 
     >         goto 301
         if (ielem(2,iel)+ielem(3,iel).ne.ielem(2,jel)+ielem(3,jel)) 
     >         goto 301
      if (iimpi.ne.0) write (6,*) 'elimination elem seg en double'
        nbnl(jel)=0
 301    continue
 300  continue
      endif
*  elimination elem en double si point point
      do 310 iel=1,ielem(/2)
       if (nbnl(iel).ne.3) goto 310
        do 311 jel=iel+1,ielem(/2)
        if (nbnl(jel).ne.3) goto 311
        if (ielem(3,iel).ne.ielem(3,jel)) goto 311
        if (ielem(2,iel).ne.ielem(2,jel)) goto 311
        if (iimpi.ne.0) write(6,*) 'elimination elem poin en double'
        nbnl(jel)=0
 311    continue
 310  continue

      nrigel=0
      segini mrigid
      ichole = 0
      imgeo1 = 0
      imgeo2 = 0
      isupeq = 0
      iforig = ifour
      mtymat='RIGIDITE'
*
      nbsous=0
      nbref=0
      do 250 iel=1,ielem(/2)
       if (nbnl(iel).eq.0) goto 250
       nbnn=nbnl(iel)
       nbelem=0
       do 255 jel=iel,ielem(/2)
        if (nbnl(jel).eq.nbnn) nbelem=nbelem+1
        if (nbelem.ge.nblim) goto 256
 255   continue
 256   continue
       segini meleme
       itypel=22
       nligrd=(nbnn-1)*idim+1
       nligrp=nligrd
       nelrig=nbelem
       RIGREL=0
       segini xmatri
*
       segini descr
        lisinc(1)=des1.lisinc(1)
        lisdua(1)=des1.lisdua(1)
        noelep(1)=1
        noeled(1)=1
        do inc=2,nligrd
         lisinc(inc)=des1.lisinc(mod(inc-2,idim)+2)
         lisdua(inc)=des1.lisdua(mod(inc-2,idim)+2)
         noelep(inc)=(inc-2)/idim+2
         noeled(inc)=(inc-2)/idim+2
         enddo

       nrigel=nrigel+1
       segadj mrigid
       irigel(1,nrigel)=meleme
       irigel(3,nrigel)=descr
       irigel(4,nrigel)=xmatri
       irigel(6,nrigel)=1
       coerig(nrigel)=1.d0
       kel=0
       do 260 jel=iel,ielem(/2)
        if (nbnl(jel).ne.nbnn) goto 260
        kel=kel+1
        do 265 in=1,nbnn
         num(in,kel)=ielem(in,jel)
 265    continue
         if (icpr(num(1,kel)).lt.0) icolor(kel)=1
***     if (idim.eq.2.or.nbnn.eq.2) then
***      xnorm2=(abs(val(1,jel))+abs(val(3,jel)))**2+
***  >         (abs(val(2,jel))+abs(val(4,jel)))**2
***     else
***      xnorm2=(abs(val(1,jel))+abs(val(4,jel))+abs(val(7,jel)))**2+
***  >          (abs(val(2,jel))+abs(val(5,jel))+abs(val(8,jel)))**2+
***  >          (abs(val(3,jel))+abs(val(6,jel))+abs(val(9,jel)))**2
***     endif
***     xnorm=sqrt(xnorm2)
***     if (mchpoi.eq.0) xnorm=1.d0
        xnorm=1.d0
        if (mchpoi.ne.0) xnorm=dnorm(abs(icpr(num(1,kel))))
C
        do 270 in=1,idim*(nbnn-1)
         re(in+1,1,kel)=val(in,jel)/xnorm
         re(1,in+1,kel)=val(in,jel)/xnorm
 270    continue
        nbnl(jel)=0
        if (kel.ge.nblim) goto 261
 260   continue
 261   continue
 250   continue
       segsup des1
**    call prrigi(mrigid,0)

*  et maintenant le second membre
       if (mchpoi.ne.0) then
         nc=MSOUP1.NOCOMP(/2)
         n=nblag
         segini,msoupo,mpoval
         nocomp(1)='FLX '
         nocomp(2)='TAIL '
         if (nc.eq.3) then
           nocomp(3)='FADH'
         endif
         nbnn=1
         nbelem=nblag
         segini meleme
         itypel=1
CCCC     mpova1=msoup1.ipoval
CCCC     segact mpova1
         ipt1=msoup1.igeoc
         segact ipt1
         do i=1,ipt1.num(/2)
           lag=ipt1.num(1,i)
           kel=abs(icpr(lag))
C
           num(1,kel)=lag
           vpocha(kel,1)=mpova1.vpocha(i,1)/dnorm(kel)+vpocha(kel,1)
C
           if (idim.ne.3) then
             vpocha(kel,2)=mpova1.vpocha(i,2)+vpocha(kel,2)
           else
             vpocha(kel,2)=sqrt(2.D0*mpova1.vpocha(i,2))+vpocha(kel,2)
           endif
C
           if (nc.eq.3) then
             vpocha(kel,3)=mpova1.vpocha(i,3)/dnorm(kel)+vpocha(kel,3)
           endif 
         enddo
         segsup,msoup1,mpova1,ipt1
         igeoc=meleme
         ipoval=mpoval
         ipchp(1)=msoupo
       endif
***
**    call ecchpo(mchpoi,0)
       segsup val,ielem,nbnl,dnorm
       return
       end