Sources Esope | Cast3M

mucpri
C MUCPRI    SOURCE    PV090527  23/03/21    21:15:09     11610          
      subroutine mucpri(mchpo1,mrigid,mchpoi)
C
C  ****  multiplication d'une matrice(mrigid) par un champoin (mchpo1)
C  ****  le resultat est un champoin (mchpoi).
C  ****  iret=ire2*ire1
C  ****  le champpoint resultat a la meme dimension que la matrice.
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC SMCHPOI
-INC SMRIGID
-INC SMCOORD
      segment itrav1
*  liste des inconnues duales dans le champ etendu
        character*(LOCOMP) nocod(nctod)
        integer ifod(nctod)
      endsegment
      segment itrav2
*  position des comps du chpo cree dans le nocod, nb pts par paquet du chpo
        integer liscomp(nctod,nsoupo),nbp(nsoupo)
      endsegment
      segment itrav4
*  champ dual  etendu
        real*8 ivald(nctod,nbpts)
      endsegment
      character*4 cnoha
      integer*4   inoha
      equivalence(cnoha,inoha)
      data cnoha/'NOHA'/
*
*  constitution de la liste des composantes duales
*
      segact mrigid
      if (imgeo2.ne.0) goto 1000
*  creation d'un champoin prototype
*  attention les mpoval ne sont pas crees
      segact mrigid*mod
      nrigel=irigel(/2)
      nctod=0
      do 27 irige=1,nrigel
       descr=irigel(3,irige)
       segact descr
       nctod=nctod+lisdua(/2)
  27  continue
      segini itrav1
      nctod=0
      do 30 irige=1,nrigel
       descr=irigel(3,irige)
       do 32 i=1,lisdua(/2)
         do 33 j=1,nctod
          if (ifod(j).ne.irigel(5,irige)) goto 33
          if (nocod(j).eq.lisdua(i)) goto 31
  33     continue
         nctod=nctod+1
         nocod(nctod)=lisdua(i)
         ifod(nctod)=irigel(5,irige)
  31   continue
  32   continue
  30  continue
*
*  expansion du chpo dual de la rigidite
*
      segini itrav4
      do 40 irige=1,nrigel
        descr=irigel(3,irige)
        ipt2=irigel(1,irige)
        segact ipt2
        do 45 ic=1,lisdua(/2)
           do 47 ir=1,nctod
             if (ifod(ir).ne.irigel(5,irige)) goto 47
             if (nocod(ir).eq.lisdua(ic)) goto 50
  47       continue
           goto 45
  50       continue
           if (irigel(5,irige).eq.inoha) then
              ivava=2**17+2**16
           else
              ivava=2**17+ifod(ir)
           endif
           do 55 iel=1,ipt2.num(/2)
              ivald(ir,ipt2.num(noeled(ic),iel))=ivava
  55       continue
  45    continue
  40  continue
*
*
*  creation du champoin resultant (sans les valeurs)
*
      nat=1
      nsoupo=10
      nsoupr=0
      segini mchpoi,itrav2
      segact mchpo1
*     jattri(1) = (mchpo1.jattri(1) * 4 )/ 2
*  par defaut les chpts issus d'un produit seront discrets !  PV
      jattri(1) = 2
      mochde='créé par mucpri'
      mtypoi=' '
*  Le(s) chpoint(s) resultat(s) ont IFOPOI = IOFRIG
*OLD      ifopoi = mchpo1.ifopoi
      ifopoi = iforig
      if (iforig .ne. mchpo1.ifopoi) then
        interr(1)=mchpo1.ifopoi
        interr(2)=iforig
        interr(3)=ifour
c-dbg      write(ioimp,*) '1132 mucpri',mchpo1,mrigid
        call erreur(1132)
        ifopoi = ifour
      end if
*  on sauve de suite dans la raideur le prototype
      imgeo2=mchpoi
*
      do 100 i=1,nbpts
        do 110 ir=1,nctod
         if (ivald(ir,i).ne.0) goto 120
 110     continue
        goto 100
 120  continue
*  point a garder test liste des composantes
      do 130 isoupo=nsoupr,1,-1
        msoupo=ipchp(isoupo)
*  regarder si meme composantes et harmonique que isoupo
        do 140 ic=1,nocomp(/2)
         if (ivald(liscomp(ic,isoupo),i).ne.2**17+noharm(ic)) goto 130
 140    continue
*  verifier pas d'autres composantes
        indtot=0
        do 145 ir=1,nctod
          if (ivald(ir,i).ne.0) indtot=indtot+1
 145    continue
        if (indtot.ne.nocomp(/2)) goto 130
*  ok on ajoute le pt dans sa liste
        nbp(isoupo)=nbp(isoupo)+1
        meleme=igeoc
        num(1,nbp(isoupo))=i
        goto 100
 130  continue
*  creation d'un nouveau paquet
 150  continue
      nsoupr=nsoupr+1
      if (nsoupr.gt.nsoupo) then
        nsoupo=nsoupr+10
        segadj mchpoi,itrav2
      endif
      nc=0
      do 160 ir=1,nctod
       if (ivald(ir,i).eq.0) goto 160
       nc=nc+1
       liscomp(nc,nsoupr)=ir
 160  continue
      segini msoupo
      ipchp(nsoupr)=msoupo
      do 165 ic=1,nc
       nocomp(ic)=nocod(liscomp(ic,nsoupr))
       noharm(ic)=ifod(liscomp(ic,nsoupr))
 165  continue
      nbelem=nbpts
      nbnn=1
      nbsous=0
      nbref=0
      segini meleme
      itypel=1
      igeoc=meleme
      num(1,1)=i
      nbp(nsoupr)=1
      goto 100
 100  continue
      do 170 isoup=1,nsoupr
        msoupo=ipchp(isoup)
        meleme=igeoc
        nbelem=nbp(isoup)
        if(meleme.ne.num(/2)) then
        nbnn=1
        nbsous=0
        nbref=0
        segini ipt1
        ipt1.itypel=itypel
        do i=1,nbelem
         ipt1.num(1,i)=num(1,i)
        enddo
        segsup meleme
        meleme=ipt1
        endif
* si possible ne pas creer de nouveau meleme
        call crech1(meleme,1)
        igeoc=meleme
 170  continue
      segsup itrav1,itrav2,itrav4
      nsoupo=nsoupr
      segadj mchpoi  
 1000 continue
      call mucpr1(mchpo1,mrigid,mchpoi)
      return
      end