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

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