Télécharger la plaquette Cast3M


fr en

Télécharger relasi.eso

Retour à la liste

Numérotation des lignes : 

relasi

  1. C RELASI    SOURCE    PV        19/02/27    21:15:06     10128          
  2. *  simplification des relations
  3. *  on elimine dans les relations les noeuds dont toutes les contributions sont < crit
  4. *
  5.       subroutine relasi(mrigid)
  6. *
  7.       implicit real*8 (a-h,o-z)
  8.       logical travail
  9.  
  10. -INC PPARAM
  11. -INC CCOPTIO
  12. -INC SMRIGID
  13. -INC SMELEME
  14.  
  15. *  segment de travail pour le regroupement des relations
  16. *  garde indique en quel position du nouveau re on va
  17. *  eleg indique la nouvelle position dans l'element d'un noeud garde
  18.       segment trav
  19.        integer garde(nligrd)
  20.        integer eleg(nbnn)
  21.       endsegment
  22.       character*4 inco1,inco2
  23.  
  24.       if (iimpi.eq.7) write(6,*) 'entree dans relasi'
  25.       crit = 1d-10
  26.       travail=.false.
  27.       segact mrigid
  28.       nrigel=0
  29.       nrige1=0
  30. *  ri1 va contenir les relations  simplifiees
  31. *  ri2 va contenir celles que l'on garde
  32.       segini ri1
  33.       segini,ri2=mrigid
  34.       nbst=0
  35.       do 200 ir=1,irigel(/2)
  36. **     write (6,*) ' boucle 200 ',ir
  37.         meleme=irigel(1,ir)
  38.         descr=irigel(3,ir)
  39.         xmatri=irigel(4,ir)
  40.         segact meleme,descr
  41.         segact xmatri*mod
  42.         if (itypel.ne.22) then
  43.           goto 200
  44.         endif
  45. **   fait au vol si besoin
  46.         ipt2=meleme
  47. **      segini,ipt2=meleme
  48.           nligrp=re(/1)
  49.           nligrd=re(/2)
  50.           nelrig=re(/3)
  51.         xmatr2=xmatri
  52. **      segini,xmatr2
  53.         ri2.irigel(1,ir)=ipt2
  54.         ri2.irigel(4,ir)=xmatr2
  55.         nbnn=num(/1)
  56.         nbref=0
  57.         nbsous=0
  58.         nbsimp=0
  59.         nbann=0
  60.         nbnn=num(/1)
  61.         segini trav
  62.         do 100 iel=1,re(/3)
  63. *  test deja traite
  64.          if (ipt2.num(1,iel).eq.0) goto 100
  65.          iok=1
  66.          remax=0.d0
  67. *  au cas ou la relation soit normalisee a autre chose que 1:
  68.          do iv=2,re(/1)
  69.           remax=max(remax,abs(re(iv,1,iel)))
  70.          enddo
  71.          critl = crit*remax
  72.  
  73. *  test (et correction) qu'une inconnue apparait plusieurs fois dans la relation
  74.          do iv1=1,nligrd
  75.            ip1=num(noeled(iv1),iel)
  76.            inco1=lisdua(iv1)
  77.          do iv2=iv1+1,nligrd
  78.            ip2=num(noeled(iv2),iel)
  79.            inco2=lisdua(iv2)
  80.            if (ip1.eq.ip2.and.inco1.eq.inco2) then 
  81.             re(1,iv1,iel)=re(1,iv1,iel)+re(1,iv2,iel)
  82.             re(1,iv2,iel)=0.d0
  83.             re(iv1,1,iel)=re(iv1,1,iel)+re(iv2,1,iel)
  84.             re(iv2,1,iel)=0.d0
  85.            endif
  86.          enddo
  87.          enddo
  88. *  y a t il des simplifications dans cette relation?
  89.          do iv=2,re(/1)
  90.           if (abs(re(iv,1,iel)).gt.critl) iok=iok+1
  91.          enddo
  92.          if (iok.ne.re(/2)) then
  93. *  il y a des termes a eliminer. on cree donc un nouveau paquet
  94.           nligrp=iok
  95.           nligrd=iok
  96.           nelrig=re(/3)
  97.           segini des1
  98.           segini xmatr1
  99.           nrige1=nrige1+1
  100.           if (nrige1.gt.nrigel) then
  101.            nrigel=nrige1+10000
  102.            segadj ri1
  103.           endif
  104.           ri1.coerig(nrige1)=coerig(ir)
  105.           ri1.irigel(2,nrige1)=irigel(2,ir)
  106.           ri1.irigel(3,nrige1)=des1
  107.           ri1.irigel(4,nrige1)=xmatr1
  108.           do i=5,irigel(/1)
  109.            ri1.irigel(i,nrige1)=irigel(i,ir)
  110.           enddo
  111. * on recopie le lx
  112.           ivc=1
  113.           des1.lisinc(1)=lisinc(1)
  114.           des1.lisdua(1)=lisdua(1)
  115.           des1.noelep(1)=noelep(1)
  116.           des1.noeled(1)=noeled(1)
  117. *  recopie des valeurs en notant ou elles vont
  118.           do iv=2,re(/1)
  119.            if (abs(re(iv,1,iel)).gt.critl) then
  120.             ivc=ivc+1
  121.             garde(iv)=ivc
  122.             xmatr1.re(1,ivc,1)=re(iv,1,iel)
  123.             xmatr1.re(ivc,1,1)=re(iv,1,iel)
  124.             des1.lisinc(ivc)=lisinc(iv)
  125.             des1.lisdua(ivc)=lisdua(iv)
  126.             des1.noelep(ivc)=noelep(iv)
  127.             des1.noeled(ivc)=noeled(iv)
  128.            else
  129.             garde(iv)=0
  130.            endif
  131.           enddo
  132. *  recreer le maillage support
  133.           nbnn=num(/1)
  134.           nbelem=re(/3)
  135.           nbsous=0
  136.           nbref=0
  137.           segini ipt1
  138.           ipt1.itypel=22
  139.           ipt1.num(1,1)=num(1,iel)
  140.           do in=2,nbnn
  141.             ipt1.num(in,1)=-num(in,iel)
  142.           enddo
  143.           do ivc=2,des1.noelep(/1)
  144.             ip=des1.noelep(ivc)
  145.             ipt1.num(ip,1)=abs(ipt1.num(ip,1))
  146.           enddo
  147. *  dans ipt1 en negatif les noeud a supprimer
  148.           inc=0
  149. *  dans eleg la nouvelle position des noeuds
  150.           do in=1,ipt1.num(/1)
  151.             if (ipt1.num(in,1).gt.0) then
  152.                inc  = inc+1
  153.                eleg(in)=inc
  154.                ipt1.num(inc,1)=ipt1.num(in,1)
  155.                do ivc=1,des1.noelep(/1)
  156.                 if (des1.noelep(ivc).eq.in) then
  157.                   des1.noelep(ivc)=inc
  158.                   des1.noeled(ivc)=inc
  159.                 endif
  160.                enddo
  161.             else
  162.                eleg(in)=0
  163.             endif
  164.           enddo
  165.           nbnn=inc
  166.            segadj ipt1
  167. * maillage reconstitue et descr mis a jour
  168.            ri1.irigel(1,nrige1)=ipt1
  169. * comme la relation a ete transferee dans un nouveau paquet on l'annule
  170. * on cree ipt2 et xmatr2 si ce n'a pas deja ete fait
  171.            if (ipt2.eq.meleme) then
  172.             segini,ipt2=meleme
  173.             ri2.irigel(1,ir)=ipt2
  174.            endif
  175.            ipt2.num(1,iel)=0
  176.            travail=.true.
  177. *  maintenant, recherche des relations ayant les memes simplifications
  178. **        write (6,*) 'simplification',iel
  179.           kel=1
  180.           do 110 jel=iel+1,re(/3)
  181.           remax=0.d0
  182. *  au cas ou la relation soit normalisee a autre chose que 1:
  183.           do iv=2,re(/1)
  184.            remax=max(remax,abs(re(iv,1,jel)))
  185.           enddo
  186.           critl = crit*remax
  187. * deja traite?
  188.           if (ipt2.num(1,jel).eq.0) goto 110
  189.           do iv=2,re(/1)
  190. * meme simplification?
  191.            if (abs(re(iv,1,jel)).gt.critl.and.garde(iv).eq.0)goto 110
  192.            if (abs(re(iv,1,jel)).le.critl.and.garde(iv).ne.0)goto 110
  193.           enddo
  194. **        write (6,*) 'meme simplification',jel
  195. * a rajouter dans xmatr1.re
  196.           kel=kel+1
  197.           do iv=2,re(/1)
  198.            if (garde(iv).ne.0) then
  199.             ivc=garde(iv)
  200.             xmatr1.re(1,ivc,kel)=re(iv,1,jel)
  201.             xmatr1.re(ivc,1,kel)=re(iv,1,jel)
  202.            endif
  203.           enddo
  204.           do in=1,num(/1)
  205.             if (eleg(in).ne.0) then
  206.               ipt1.num(eleg(in),kel)=num(in,jel)
  207.             endif
  208.           enddo
  209. * annuler l'element
  210.           ipt2.num(1,jel)=0
  211.  110      continue
  212.           nelrig=kel
  213.           segadj xmatr1
  214.           nbelem=kel
  215.           segadj ipt1
  216.           nbsimp=nbsimp+kel
  217.          endif
  218.  100    continue
  219.          segsup trav
  220. * fin du do iel
  221. *  maintenant, on retasse ri2.irigel(ir)
  222.         if (ipt2.ne.meleme) then
  223.             nligrp=re(/1)
  224.             nligrd=re(/2)
  225.             nelrig=re(/3)-nbsimp
  226.             segini,xmatr2
  227.             ri2.irigel(4,ir)=xmatr2
  228.           nbelem=0
  229.           nbnn=ipt2.num(/1)
  230.           nbsous=0
  231.           nbref=0
  232.           do iel=1,ipt2.num(/2)
  233.            if (ipt2.num(1,iel).ne.0) then
  234.             nbelem=nbelem+1
  235.             do in=1,ipt2.num(/1)
  236.               ipt2.num(in,nbelem)=num(in,iel)
  237.             enddo
  238.             do ip=1,xmatr2.re(/2)
  239.             do id=1,xmatr2.re(/1)
  240.               xmatr2.re(id,ip,nbelem)=re(id,ip,iel)
  241.             enddo
  242.             enddo
  243.            endif
  244.           enddo
  245.           if (nbelem.ne.nelrig) call erreur(5)
  246.           segadj ipt2
  247.         endif
  248.         nbst = nbst + nbsimp
  249. *       write(6,*) ' regroupement de relation ',nbsimp
  250.  200  continue
  251. *  plus qu'a tout assembler dans le resultat.
  252. *  au passage on elimine les matrices vides
  253. *  mais avant, si on n'a rien eu a faire, on sort la matrice initiale qui est toujours dans mrigid
  254.       if (.not.travail) return
  255.  
  256.       nrige2=ri2.irigel(/2)
  257.       if (iimpi.eq.7) write(6,*) ' nombre de relations simplifiees ',
  258.      >      nbst
  259. **    write(6,*) ' simplifications retenues annulees ',nrige1,nbann
  260.       nrigel=nrige1+nrige2
  261.       segadj ri2
  262.       nrigel=0
  263.       do ir=1,nrige2
  264.         meleme=ri2.irigel(1,ir)
  265.         segact meleme
  266.         if (num(/2).ne.0) then
  267.           nrigel=nrigel+1
  268.           do is=1,ri2.irigel(/1)
  269.             ri2.irigel(is,nrigel)=ri2.irigel(is,ir)
  270.           enddo
  271.           ri2.coerig(nrigel)=ri2.coerig(ir)
  272.         endif
  273.       enddo
  274.       do ir=1,nrige1
  275.         meleme=ri1.irigel(1,ir)
  276.         segact meleme
  277.         if (num(/2).ne.0) then
  278.           nrigel=nrigel+1
  279.           do is=1,ri2.irigel(/1)
  280.             ri2.irigel(is,nrigel)=ri1.irigel(is,ir)
  281.           enddo
  282.           ri2.coerig(nrigel)=ri1.coerig(ir)
  283.         endif
  284.       enddo
  285.       segsup ri1
  286.       segadj ri2
  287.  
  288.       mrigid = ri2
  289.  
  290.       return
  291.       end
  292.  
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  

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