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    JK148537  24/11/05    21:15:08     12067          
  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.        else
  260.          segact xmatri
  261.        endif
  262.       nod=noeled(/1)
  263.       nop=noelep(/1)
  264. *     if (nod.ne.nop) isyme=2
  265.        do 55 iel=1,ipt2.num(/2)
  266.         ibfp=0
  267.         do 60 ip=1,nop
  268.          if (icposp(ip).eq.0.or.icpr(ipt2.num(noelep(ip),iel)).eq.0)
  269.      >   then
  270.            xbuffp(ip)=0.d0
  271.          else
  272. *          write (6,*) ' icposp icposd ',icposp(ip),icposd(ip)
  273. *          write (6,*) ' ip xbuffp ',ip,xbuffp(ip)
  274. *          write (6,*) ' icposp noelep ',icposp(ip),noelep(ip)
  275. *          write (6,*) ' irige coerig ',irige,coerig(irige)
  276. *          write (6,*) ' ipt2 ',ipt2.num(noelep(ip),iel)
  277. *          write (6,*) ' icpr ',icpr(ipt2.num(noelep(ip),iel))
  278. *          write (6,*) ' xvalp ',xvalp(icposp(ip),
  279. *    >                  icpr(ipt2.num(noelep(ip),iel)))
  280. **           xbuffp(ip)=xvalp(icposp(ip),icpr(ipt2.num(noelep(ip),iel)))*
  281.            xbuffp(ip)=xvalp(icpr(ipt2.num(noelep(ip),iel)),icposp(ip))*
  282.      >                coerig(irige)
  283.            if (abs(xbuffp(ip)).gt.xpetit) ibfp=ibfp+1
  284.          endif
  285.   60    continue
  286.  
  287. *  operation
  288. *  on impose l'egalite des multiplicateurs de Lagrange dans mucpr2
  289.         if (nop*nod.ne.0.and.ibfp.ne.0) then
  290.           call mucpr2(nop,nod,re(1,1,iel),xbuffp,
  291.      $    xbuffd,isyme)
  292. *  transfert du buffer dual dans le vecteur dual etendu
  293.          do 65 id=1,nod
  294.           i1d=icposd(id)
  295.           if (i1d.eq.0) then
  296. *           write (6,*) ' i1d nul '
  297.            goto 65
  298.           endif
  299.           i2d=icpr(ipt2.num(noeled(id),iel))
  300.           if (i2d.eq.0) goto 65
  301. **        xvald(i1d,i2d)=xvald(i1d,i2d)+xbuffd(id)
  302.           xvald(i2d,i1d)=xvald(i2d,i1d)+xbuffd(id)
  303.    65    continue
  304. **      else
  305. **        do id=1,nod
  306. **         xbuffd(id)=0.D0
  307. **        enddo
  308.         endif
  309.   55   continue
  310. *       if (xmatri.ne.0) segdes xmatri
  311.        segdes descr,xmatri
  312.   50   continue
  313.  
  314.   51   continue
  315. *
  316. *  transfert final dans le vecteur dual
  317. *
  318.       isous=0
  319.       do 80 isouso=1,ipchp(/1)
  320.        isous=isous+1
  321.        msoupo=ipchp(isouso)
  322.         mpoval=ipoval
  323.         meleme=igeoc
  324.         segact meleme
  325.        do 82 i=1,nocomp(/2)
  326.         do 83 j=1,nctod
  327.          if (ifod(j).ne.noharm(i).and.ifodh(j).ne.noharm(i)) goto 83
  328.          if (nocod(j).eq.nocomp(i)) goto 84
  329.   83    continue
  330.         ifodp(i)=0
  331.         goto 82
  332. *        write(6,*) 'mucpr1 - 4'
  333.         call erreur(5)
  334.   84    continue
  335.         if (ifodh(j).ne.0) noharm(i)=ifodh(j)-2**17
  336. *        write (charm,fmt='(a4)') noharm(i)
  337.          if (charm.eq.'NOHA') noharm(i)=nifour
  338.         ifodp(i)=j
  339.   82   continue
  340.        ieln=0
  341.        ipt1=0
  342.        do 85 iel=1,vpocha(/1)
  343.        ieln=ieln+1
  344.        if (ieln.ne.iel) num(1,ieln)=num(1,iel)
  345.        ig=0
  346.        do 86 i=1,vpocha(/2)
  347.          if (ifodp(i).eq.0) then
  348. *          write (6,*) ' ifodp nul '
  349.           goto 86
  350.          endif
  351. **       xv=xvald(ifodp(i),icpr(num(1,iel)))
  352.          xv=xvald(icpr(num(1,iel)),ifodp(i))
  353.          vpocha(ieln,i)=xv
  354. *         if (abs(xv).gt.1e-35) ig=1
  355.          if (abs(xv).gt.xpetit) ig=1
  356.   86   continue
  357.        if (ig.eq.0) then
  358. *  il y a lieu de simplifier
  359.         if (ipt1.eq.0) then
  360.          segini,ipt1=meleme
  361.          meleme=ipt1
  362.         endif
  363.         ieln=ieln-1
  364.        endif
  365.   85   continue
  366.       if (ieln.ne.vpocha(/1)) then
  367.        nbnn  =1
  368.        nbelem=ieln
  369.        nbref =0
  370.        nbsous=0
  371.        segadj,meleme
  372.        n=ieln
  373.        nc=vpocha(/2)
  374. *      write (6,*) ' chpoin simplifie ',vpocha(/1),ieln
  375.        segadj mpoval
  376.       endif
  377.       if (meleme.ne.igeoc) call crech1(meleme,1)
  378.       igeoc=meleme
  379.       segact,mpoval*nomod,meleme*nomod,msoupo*nomod
  380.       if (ieln.eq.0) then
  381.        isous=isous-1
  382.        segsup mpoval,msoupo
  383.       else   
  384.        ipchp(isous)=ipchp(isouso)
  385. *      write(6,*) 'isous ipchp ',isous,ipchp(isous),ipchp(/1)
  386.       endif
  387.   80  continue
  388. *OLD      mchpoi.ifopoi = mchpo1.ifopoi
  389.       if (isous.ne.ipchp(/1)) then
  390.        nat=jattri(/1)
  391.        nsoupo=isous 
  392.        segadj mchpoi
  393.       endif
  394.       mchpoi.ifopoi = mrigid.iforig
  395.       segact,mchpoi*nomod
  396.       segdes,mrigid
  397.       segsup itrav1,itrav2,itrav3,itrav4,itrav5,icpr
  398.       end
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  

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