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

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