Sources Esope | Cast3M

reduri
C REDURI    SOURCE    FANDEUR   22/01/19    21:15:13     11256          
      subroutine reduri(mrigid,moleme,irigre)
      implicit real*8(a-h,o-z)
      implicit integer (i-n)
-INC SMRIGID
-INC SMELEME
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
      segment icpr(nbpts)
      segment inode(ino)
      segment jelnum(imaxel,ino)
      segment izone(imaxel,ino)
      irigre=0
      ipt1=moleme
      segact ipt1
      meleme=ipt1
      nbso1=ipt1.lisous(/1)
      ino=0
 
C   Cas d'un MELEME vide => RIGIDITE vide
      segact mrigid
      IF (nbso1.eq.0) THEN
        IF (IPT1.NUM(/2).EQ.0) THEN
          NRIGEL = 0
          SEGINI,RI1
          RI1.MTYMAT = MTYMAT
          RI1.IFORIG = IFORIG
          irigre = RI1
          segdes mrigid
          RETURN
        ENDIF
      ENDIF
*
*   creation, d'une numerotation locale
*
      segini icpr
      do  1 i =1, max(nbso1,1)
         if(nbso1.ne.0) then
            meleme=ipt1.lisous(i)
            segact meleme
         endif
         do 2 j=1,num(/2)
           do 22 k=1,num(/1)
             ia= num(k,j)
             if(icpr(ia).eq.0) then
               ino=ino+1
               icpr(ia)=ino
             endif
   22      continue
    2   continue
    1 continue
*
* on compte combien d'elements touche un noeud
*
      segini inode
      do 3 i=1,max(nbso1,1)
         if(nbso1.ne.0) then
            meleme=ipt1.lisous(i)
         endif
         do 4 j=1,num(/2)
           do 42 k=1,num(/1)
             ia=num(k,j)
             ib= icpr(ia)
             inode(ib)=inode(ib)+1
   42      continue
    4    continue
    3 continue
      imaxel=0
      do i=1,ino
        imaxel=max(imaxel,inode(i))
        inode(i)=0
      enddo
      segini jelnum,izone
*
* jelnum(i,k) dira le Ieme element qui touche le noeud k
* izone (i,k) dira dans quel lisous le ieme element se trouve
* attention le noeud K est le plus petit de l'élément
      do 5 i=1,max(nbso1,1)
         if(nbso1.ne.0) then
            meleme=ipt1.lisous(i)
         endif
         do 6 j=1,num(/2)
           ipeti=nbpts*(idim+1)
           do  k=1,num(/1)
              ipeti=min(num(k,j),ipeti)
           enddo
             ib= icpr(ipeti)
             inode(ib)=inode(ib)+1
             ic=inode(ib)
             jelnum(ic,ib)=j
             izone(ic,ib)=i
    6    continue
    5 continue
*
*   travail sur les rigidites elementaires
*
      nrigel=coerig(/1)
      nrige=irigel(/1)
      segini ri1
      ri1.mtymat=mtymat
      ri1.iforig=iforig
      irzon=0
      do 10 i=1,irigel(/2)
        ipt3=irigel(1,i)
        segact ipt3
        xmatri=irigel(4,i)
        ielzo=0
        do 11 j=1,ipt3.num(/2)
          ipeti=nbpts*(idim+1)
          do k=1,ipt3.num(/1)
            ipeti=min(ipeti,ipt3.num(k,j))
          enddo
 
* on regarde s'il existe un element de ipt1 ayant ce noeud en plus
* petitre position si non on passe a l'élément suivant
          ib=icpr(ipeti)
          if(ib.eq.0) go to 11
          if(inode(ib).eq.0) go to 11
          do 13 mm=1,inode(ib)
           if(nbso1.ne.0)then
             meleme= ipt1.lisous(izone(mm,ib))
           endif
           if(ipt3.num(/1).ne.num(/1)) go to 13
           iel=jelnum(mm,ib)
           do  in=1,ipt3.num(/1)
             if( ipt3.num(in,j).ne.num(in,iel))go to 13
           enddo
* on a trouver un element a conserver
          if(ielzo.eq.0) then
            segini,ipt4=ipt3
            irzon=irzon+1
            do kk=1,irigel(/1)
              ri1.irigel(kk,irzon)=irigel(kk,i)
            enddo
            ri1.irigel(1,irzon)=ipt4
            ri1.coerig(irzon)=coerig(i)
            xmatri=irigel(4,i)
            segact xmatri
            segini,xmatr1=xmatri
            ri1.irigel(4,irzon)= xmatr1
          endif
          ielzo=ielzo+1
              do io=1,re(/2)
                 do iu=1,re(/1)
                     xmatr1.re(iu,io,ielzo)=re(iu,io,j)
                 enddo
              enddo
*          imatr1.imattt(ielzo)=imattt(j)
          do kk=1,ipt3.num(/1)
            ipt4.num(kk,ielzo)=ipt3.num(kk,j)
          enddo
   13     continue
   11   continue
* on ajuste les longueur si besoin
        if(ielzo.ne.0) then
          if(ielzo.ne.ipt4.num(/2)) then
           nbsous=0
           nbelem= ielzo
           nbnn=ipt4.num(/1)
           nbref=0
           segadj ipt4
           nelrig=ielzo
           nligrp= xmatr1.re(/2)
           nligrd= xmatr1.re(/1)
           segadj xmatr1
           segdes ipt4,xmatr1
          else
           ri1.irigel(4,irzon)=xmatri
           ri1.irigel(1,irzon)=ipt3
           segsup xmatr1,ipt4
          endif
        endif
        segdes xmatri
   10 continue
      if(irzon.eq.0) then
        segdes mrigid
        call erreur(21)
        return
      elseif(irzon.ne.coerig(/1)) then
         nrigel=irzon
         nrige=irigel(/1)
         segadj ri1
         segdes ri1
         irigre=ri1
      else
         do io=1,irzon
          if( ri1.irigel(1,io).ne.irigel(1,io) )go to 20
         enddo
         irigre= mrigid
         go to 21
   20    continue
         irigre=ri1
 
   21    continue
      endif
      do iou=1,irigel(/2)
        ipt3=irigel(1,iou)
        segdes ipt3
      enddo
      segdes mrigid
      if(irigre.eq.mrigid) then
         segsup ri1
      else
       segdes ri1
      endif
      segsup izone,icpr,jelnum,inode
c      return
      end