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

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