Télécharger la plaquette Cast3M


fr en

Télécharger mucpr1.eso

Retour à la liste

Numérotation des lignes : 

mucpr1

  1. C MUCPR1    SOURCE    GOUNAND   25/05/05    21:15:07     12259          
  2.       subroutine mucpr1(mchpo1,mrigid,mchpoi)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8(A-H,O-Z)
  5. *  multiplication rapide d'une matrice par un champoint dans le cas
  6. *  ou on connait deja la structure du chanpoint resultant, qui a par
  7. *  exemple ete obtenu par un appel precedent a mucpri
  8. *
  9. *  entree : mrigid, mchpo1
  10. *  sortie : mchpoi
  11. *
  12. -INC CCREEL
  13.  
  14. -INC PPARAM
  15. -INC CCOPTIO
  16. -INC SMCHPOI
  17. -INC SMELEME
  18. -INC SMCOORD
  19. -INC SMRIGID
  20.       segment itrav1
  21. *  liste des inconnues primales dans le champ etendu
  22.         character*(LOCOMP) nocop(nctop)
  23.         integer ifop(nctop)
  24.       endsegment
  25.       segment itrav2
  26. *  liste des inconnues duales dans le champ etendu
  27.         character*(LOCOMP) nocod(nctod)
  28.         integer ifod(nctod),ifodh(nctod),ifodp(nctod)
  29.       endsegment
  30.       segment itrav3
  31. *  champ primal  etendu
  32. **      real*8 xvalp(nctop,ncpr )
  33.         real*8 xvalp(ncpr,nctop )
  34.       endsegment
  35.       segment itrav4
  36. *  champ dual  etendu
  37. **      real*8 xvald(nctod,ncpr )
  38.         real*8 xvald(ncpr,nctod)
  39.       endsegment
  40.       segment itrav5
  41. *  position inconnues primales et duales de la raideur dans les champs etendus
  42.         integer icposp(nbpma),icposd(nbdma)
  43. *  buffers primal et dual
  44.         real*8 xbuffp(nbpma),xbuffd(nbdma)
  45.         character*(LOCOMP) charm
  46.       endsegment
  47.       segment icpr(nbpts)
  48.  
  49.       character*4 cnoha
  50.       integer*4 inoha
  51.       equivalence(cnoha,inoha)
  52.       data cnoha/'NOHA'/
  53. *
  54. *  creation du champ resultat
  55. *
  56. *  traitement ctrlC
  57.       if(ierr.ne.0) return
  58.       segini icpr
  59.       ncpr=0
  60.       segact mrigid
  61.       nrigel=irigel(/2)
  62.       mchpo2=imgeo2
  63.       segact mchpo2
  64.       if (mchpo2.ifopoi.ne.iforig) then
  65.         moterr(1:8)='CHPO-RIG'
  66.         interr(1)=mchpo2.ifopoi
  67.         interr(2)=iforig
  68.         interr(3)=IFOUR
  69. c-dbg      write(ioimp,*) '1132 mucpr1',mchpo2,mrigid
  70.         call erreur(1132)
  71.       endif
  72.       segini,mchpoi=mchpo2
  73.       nbz=ipchp(/1)
  74.       do 1  isous=1,nbz
  75.        msoup2=ipchp(isous)
  76.        segini,msoupo=msoup2
  77.        ipchp(isous)=msoupo
  78.        meleme=igeoc
  79.        segact meleme
  80.        nelt=num(/2)
  81. *       nnoe=num(/1)
  82.        do im=1,nelt
  83. *       do in=1,nnoe
  84. *        ip=num(in,im)
  85.         ip=num(1,im)
  86.         if (icpr(ip).eq.0) then
  87.          ncpr=ncpr+1
  88.          icpr(ip)=ncpr
  89.         endif
  90. *       enddo
  91.        enddo
  92.        n=nelt
  93.        nc=noharm(/1)
  94.        segini mpoval
  95.        ipoval=mpoval
  96.    1  continue
  97. *
  98. *  constitution de la liste des composantes du chant primal
  99. *
  100.       nctop=0
  101.       segact mchpo1
  102.       if (mchpo1.ifopoi.ne.iforig) then
  103.         moterr(1:8)='CHPO-RIG'
  104.         interr(1)=mchpo1.ifopoi
  105.         interr(2)=iforig
  106.         interr(3)=ifour
  107. c-dbg      write(ioimp,*) '1132 mucpr1 (2)',mchpo1,mrigid
  108.         call erreur(1132)
  109.       endif
  110.       do 10 isoupo=1,mchpo1.ipchp(/1)
  111.        msoup1=mchpo1.ipchp(isoupo)
  112.        segact msoup1
  113.        mpova1=msoup1.ipoval
  114.        segact mpova1
  115.        ipt1=msoup1.igeoc
  116.        segact ipt1
  117.        do im=1,ipt1.num(/2)
  118. *       do in=1,ipt1.num(/1)
  119. *        ip=ipt1.num(in,im)
  120.         ip=ipt1.num(1,im)
  121.         if (icpr(ip).eq.0) then
  122.          ncpr=ncpr+1
  123.          icpr(ip)=ncpr
  124.         endif
  125. *       enddo
  126.        enddo
  127.        nctop=nctop+msoup1.nocomp(/2)
  128.   10  continue
  129.       segini itrav1
  130.       nctop=0
  131.       do 11 isous=1,mchpo1.ipchp(/1)
  132.        msoup1=mchpo1.ipchp(isous)
  133.        do 12 i=1,msoup1.nocomp(/2)
  134.          do 13 j=1,nctop
  135.           if (ifop(j).ne.msoup1.noharm(i)) goto 13
  136.           if (nocop(j).eq.msoup1.nocomp(i)) goto 12
  137.   13     continue
  138.          nctop=nctop+1
  139.          nocop(nctop)=msoup1.nocomp(i)
  140.          ifop(nctop)=msoup1.noharm(i)
  141. c       write(6,*) 'la',nctop,nocop(nctop),ifop(nctop),i,isous
  142.   12   continue
  143.   11   continue
  144. *     write (6,*) ' mucpr1 nbpts,ncpr,nctop ',
  145. *    >  nbpts,ncpr,nctop
  146. *
  147. *  constitution de la liste des composantes du champ dual
  148. *
  149.       nctod=0
  150.       do 15 isoupo=1,ipchp(/1)
  151.        msoupo=ipchp(isoupo)
  152.        mpoval=ipoval
  153.        segact mpoval*mod
  154.        nctod=nctod+nocomp(/2)
  155.   15  continue
  156. *     write (6,*) ' dans mucpr1 nctod-1 ',nctod,ipchp(/1)
  157.       segini itrav2
  158.       nctod=0
  159.       do 16 isous=1,ipchp(/1)
  160.        msoupo=ipchp(isous)
  161.        do 17 i=1,nocomp(/2)
  162.          do 18 j=1,nctod
  163.           if (ifod(j).ne.noharm(i)) goto 18
  164.           if (nocod(j).eq.nocomp(i)) goto 17
  165.   18     continue
  166.          nctod=nctod+1
  167.          nocod(nctod)=nocomp(i)
  168.          ifod(nctod)=noharm(i)
  169.   17   continue
  170.   16   continue
  171. *     write (6,*) ' dans mucpr1 nctod-2 ',nctod
  172. *
  173. *  expansion du champ primal
  174. *
  175.       segini itrav3,itrav4
  176.       do 20 isous=1,mchpo1.ipchp(/1)
  177.        msoup1=mchpo1.ipchp(isous)
  178.        mpova1=msoup1.ipoval
  179.        ipt1=msoup1.igeoc
  180.        do 22 i=1,msoup1.nocomp(/2)
  181.         do 23 j=1,nctop
  182.           if (ifop(j).ne.msoup1.noharm(i)) goto 23
  183.           if (nocop(j).eq.msoup1.nocomp(i)) goto 24
  184.   23    continue
  185. c        write(6,*) 'mucpr1 - 1'
  186.         call erreur(5)
  187.   24    continue
  188.          do 25 iel=1,mpova1.vpocha(/1)
  189. **         xvalp(j,icpr(ipt1.num(1,iel)))=mpova1.vpocha(iel,i)
  190.            xvalp(icpr(ipt1.num(1,iel)),j)=mpova1.vpocha(iel,i)
  191.    25    continue
  192.   22   continue
  193.   20  continue
  194. *
  195. *  creation de itrav5
  196. *
  197.       nbpma=0
  198.       nbdma=0
  199.       do 27 irige=1,nrigel
  200.        descr=irigel(3,irige)
  201.        segact descr
  202.        nbpma=max(nbpma,lisinc(/2))
  203.        nbdma=max(nbdma,lisdua(/2))
  204.   27  continue
  205.       segini itrav5
  206.        if (nbdma.eq.0) goto 51
  207. *
  208. *  travail effectif : boucle sur la matrice.
  209. *
  210.       do 50 irige=1,nrigel
  211.         write (charm,fmt='(a4)') irigel(5,irige)
  212. cjk        write (charm,fmt='(a4)') 'NOHA'
  213.        descr=irigel(3,irige)
  214.        segact descr
  215. *  remplissage de icposp et icposd
  216.        jsauv=0
  217.        do 30 i=1,lisinc(/2)
  218.          do 35 j=1,nctop
  219. cjk           if (ifop(j).ne.irigel(5,irige).and.charm.ne.'NOHA') goto 35
  220.            if (ifop(j).ne.irigel(5,irige)) goto 35
  221.            if (lisinc(i).eq.nocop(j)) goto 37
  222.   35     continue
  223. *  inconnu du chp primal pas dans la raideur
  224. *          write(6,*) 'inconnu du chp primal pas dans la raideur'
  225.          j=0
  226.   37     continue
  227.          icposp(i)=j
  228.          if (jsauv*j.ne.0) then
  229.           if (ifop(jsauv).ne.ifop(j)) then
  230.             call erreur(5)
  231.           endif
  232.          endif
  233.          if (j.ne.0) jsauv=j
  234.   30   continue
  235.        do 40 i=1,lisdua(/2)
  236.          do 45 j=1,nctod
  237.            if (ifod(j).ne.irigel(5,irige).and.charm.ne.'NOHA') goto 45
  238.            if (lisdua(i).eq.nocod(j)) goto 47
  239.   45     continue
  240.          icposd(i)=0
  241.          goto 40
  242. *          write (6,*) ' dans mucpr1 ',lisdua(/2),nctod
  243. *          write (6,*)  irigel(5,irige),i,(lisdua(j),j=1,lisdua(/2))
  244. *          write (6,*)  (nocod(j),ifod(j),j=1,nctod)
  245. *         write(6,*) 'mucpr1 - 3'
  246.          call erreur(5)
  247.   47     continue
  248.            if (charm.eq.'NOHA'.and.jsauv.ne.0) then
  249.              ifodh(j)=ifop(jsauv)+2**17
  250.            endif
  251.          icposd(i)=j
  252.   40   continue
  253. *  transfert des valeurs dans le buffer primal
  254.        ipt2=irigel(1,irige)
  255.        segact ipt2
  256. *      write (6,*) ' ipt2 ',ipt2.num(/1),ipt2.num(/2)
  257.        xmatri=irigel(4,irige)
  258.        isyme=irigel(7,irige)
  259.        if (ipt2.itypel.eq.22) then
  260.          segact xmatri*mod
  261.          irela=1
  262.        else
  263.          segact xmatri
  264.          irela=0
  265.        endif
  266.       nod=noeled(/1)
  267.       nop=noelep(/1)
  268. *     if (nod.ne.nop) isyme=2
  269.        do 55 iel=1,ipt2.num(/2)
  270.         ibfp=0
  271.         do 60 ip=1,nop
  272.          if (icposp(ip).eq.0.or.icpr(ipt2.num(noelep(ip),iel)).eq.0)
  273.      >   then
  274.            xbuffp(ip)=0.d0
  275.          else
  276. *          write (6,*) ' icposp icposd ',icposp(ip),icposd(ip)
  277. *          write (6,*) ' ip xbuffp ',ip,xbuffp(ip)
  278. *          write (6,*) ' icposp noelep ',icposp(ip),noelep(ip)
  279. *          write (6,*) ' irige coerig ',irige,coerig(irige)
  280. *          write (6,*) ' ipt2 ',ipt2.num(noelep(ip),iel)
  281. *          write (6,*) ' icpr ',icpr(ipt2.num(noelep(ip),iel))
  282. *          write (6,*) ' xvalp ',xvalp(icposp(ip),
  283. *    >                  icpr(ipt2.num(noelep(ip),iel)))
  284. **           xbuffp(ip)=xvalp(icposp(ip),icpr(ipt2.num(noelep(ip),iel)))*
  285.            xbuffp(ip)=xvalp(icpr(ipt2.num(noelep(ip),iel)),icposp(ip))*
  286.      >                coerig(irige)
  287.            if (abs(xbuffp(ip)).gt.xpetit) ibfp=ibfp+1
  288.          endif
  289.   60    continue
  290.  
  291. *  operation
  292. *  on impose l'egalite des multiplicateurs de Lagrange dans mucpr2
  293.         if (nop*nod.ne.0.and.ibfp.ne.0) then
  294.           call mucpr2(nop,nod,re(1,1,iel),xbuffp,
  295.      $    xbuffd,isyme,irela)
  296. *  transfert du buffer dual dans le vecteur dual etendu
  297.          do 65 id=1,nod
  298.           i1d=icposd(id)
  299.           if (i1d.eq.0) then
  300. *           write (6,*) ' i1d nul '
  301.            goto 65
  302.           endif
  303.           i2d=icpr(ipt2.num(noeled(id),iel))
  304.           if (i2d.eq.0) goto 65
  305. **        xvald(i1d,i2d)=xvald(i1d,i2d)+xbuffd(id)
  306.           xvald(i2d,i1d)=xvald(i2d,i1d)+xbuffd(id)
  307.    65    continue
  308. **      else
  309. **        do id=1,nod
  310. **         xbuffd(id)=0.D0
  311. **        enddo
  312.         endif
  313.   55   continue
  314. *       if (xmatri.ne.0) segdes xmatri
  315.        segdes descr,xmatri
  316.   50   continue
  317.  
  318.   51   continue
  319. *
  320. *  transfert final dans le vecteur dual
  321. *
  322.       isous=0
  323.       do 80 isouso=1,ipchp(/1)
  324.        isous=isous+1
  325.        msoupo=ipchp(isouso)
  326.         mpoval=ipoval
  327.         meleme=igeoc
  328.         segact meleme
  329.        do 82 i=1,nocomp(/2)
  330.         do 83 j=1,nctod
  331.          if (ifod(j).ne.noharm(i).and.ifodh(j).ne.noharm(i)) goto 83
  332.          if (nocod(j).eq.nocomp(i)) goto 84
  333.   83    continue
  334.         ifodp(i)=0
  335.         goto 82
  336. *        write(6,*) 'mucpr1 - 4'
  337.         call erreur(5)
  338.   84    continue
  339.         if (ifodh(j).ne.0) noharm(i)=ifodh(j)-2**17
  340. *        write (charm,fmt='(a4)') noharm(i)
  341.          if (charm.eq.'NOHA') noharm(i)=nifour
  342.         ifodp(i)=j
  343.   82   continue
  344.        ieln=0
  345.        ipt1=0
  346.        do 85 iel=1,vpocha(/1)
  347.        ieln=ieln+1
  348.        if (ieln.ne.iel) num(1,ieln)=num(1,iel)
  349.        ig=0
  350.        do 86 i=1,vpocha(/2)
  351.          if (ifodp(i).eq.0) then
  352. *          write (6,*) ' ifodp nul '
  353.           goto 86
  354.          endif
  355. **       xv=xvald(ifodp(i),icpr(num(1,iel)))
  356.          xv=xvald(icpr(num(1,iel)),ifodp(i))
  357.          vpocha(ieln,i)=xv
  358. *         if (abs(xv).gt.1e-35) ig=1
  359.          if (abs(xv).gt.xpetit) ig=1
  360.   86   continue
  361.        if (ig.eq.0) then
  362. *  il y a lieu de simplifier
  363.         if (ipt1.eq.0) then
  364.          segini,ipt1=meleme
  365.          meleme=ipt1
  366.         endif
  367.         ieln=ieln-1
  368.        endif
  369.   85   continue
  370.       if (ieln.ne.vpocha(/1)) then
  371.        nbnn  =1
  372.        nbelem=ieln
  373.        nbref =0
  374.        nbsous=0
  375.        segadj,meleme
  376.        n=ieln
  377.        nc=vpocha(/2)
  378. *      write (6,*) ' chpoin simplifie ',vpocha(/1),ieln
  379.        segadj mpoval
  380.       endif
  381.       if (meleme.ne.igeoc) call crech1(meleme,1)
  382.       igeoc=meleme
  383.       segact,mpoval*nomod,meleme*nomod,msoupo*nomod
  384.       if (ieln.eq.0) then
  385.        isous=isous-1
  386.        segsup mpoval,msoupo
  387.       else
  388.        ipchp(isous)=ipchp(isouso)
  389. *      write(6,*) 'isous ipchp ',isous,ipchp(isous),ipchp(/1)
  390.       endif
  391.   80  continue
  392. *OLD      mchpoi.ifopoi = mchpo1.ifopoi
  393.       if (isous.ne.ipchp(/1)) then
  394.        nat=jattri(/1)
  395.        nsoupo=isous
  396.        segadj mchpoi
  397.       endif
  398.       mchpoi.ifopoi = mrigid.iforig
  399.       segact,mchpoi*nomod
  400.       segdes,mrigid
  401.       segsup itrav1,itrav2,itrav3,itrav4,itrav5,icpr
  402.       end
  403.  
  404.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales