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    CB215821  25/04/23    21:15:37     12247          
  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.           segini xmatr1
  100.           nrige1=nrige1+1
  101.           if (nrige1.gt.nrigel) then
  102.            nrigel=nrige1+10000
  103.            segadj ri1
  104.           endif
  105.           ri1.coerig(nrige1)=coerig(ir)
  106.           ri1.irigel(2,nrige1)=irigel(2,ir)
  107.           ri1.irigel(3,nrige1)=des1
  108.           ri1.irigel(4,nrige1)=xmatr1
  109.           do i=5,irigel(/1)
  110.            ri1.irigel(i,nrige1)=irigel(i,ir)
  111.           enddo
  112. * on recopie le lx
  113.           ivc=1
  114.           des1.lisinc(1)=lisinc(1)
  115.           des1.lisdua(1)=lisdua(1)
  116.           des1.noelep(1)=noelep(1)
  117.           des1.noeled(1)=noeled(1)
  118. *  recopie des valeurs en notant ou elles vont
  119.           do iv=2,re(/1)
  120.            if (abs(re(iv,1,iel)).gt.critl) then
  121.             ivc=ivc+1
  122.             garde(iv)=ivc
  123.             xmatr1.re(1,ivc,1)=re(iv,1,iel)
  124.             xmatr1.re(ivc,1,1)=re(iv,1,iel)
  125.             des1.lisinc(ivc)=lisinc(iv)
  126.             des1.lisdua(ivc)=lisdua(iv)
  127.             des1.noelep(ivc)=noelep(iv)
  128.             des1.noeled(ivc)=noeled(iv)
  129.            else
  130.             garde(iv)=0
  131.            endif
  132.           enddo
  133. *  recreer le maillage support
  134.           nbnn=num(/1)
  135.           nbelem=re(/3)
  136.           nbsous=0
  137.           nbref=0
  138.           segini ipt1
  139.           ipt1.itypel=22
  140.           ipt1.num(1,1)=num(1,iel)
  141.           do in=2,nbnn
  142.             ipt1.num(in,1)=-num(in,iel)
  143.           enddo
  144.           do ivc=2,des1.noelep(/1)
  145.             ip=des1.noelep(ivc)
  146.             ipt1.num(ip,1)=abs(ipt1.num(ip,1))
  147.           enddo
  148. *  dans ipt1 en negatif les noeud a supprimer
  149.           inc=0
  150. *  dans eleg la nouvelle position des noeuds
  151.           do in=1,ipt1.num(/1)
  152.             if (ipt1.num(in,1).gt.0) then
  153.                inc  = inc+1
  154.                eleg(in)=inc
  155.                ipt1.num(inc,1)=ipt1.num(in,1)
  156.                do ivc=1,des1.noelep(/1)
  157.                 if (des1.noelep(ivc).eq.in) then
  158.                   des1.noelep(ivc)=inc
  159.                   des1.noeled(ivc)=inc
  160.                 endif
  161.                enddo
  162.             else
  163.                eleg(in)=0
  164.             endif
  165.           enddo
  166.           nbnn=inc
  167.            segadj ipt1
  168. * maillage reconstitue et descr mis a jour
  169.            ri1.irigel(1,nrige1)=ipt1
  170. * comme la relation a ete transferee dans un nouveau paquet on l'annule
  171. * on cree ipt2 et xmatr2 si ce n'a pas deja ete fait
  172.            if (ipt2.eq.meleme) then
  173.             segini,ipt2=meleme
  174.             ri2.irigel(1,ir)=ipt2
  175.            endif
  176.            ipt2.num(1,iel)=0
  177.            travail=.true.
  178. *  maintenant, recherche des relations ayant les memes simplifications
  179. **        write (6,*) 'simplification',iel
  180.           kel=1
  181.           do 110 jel=iel+1,re(/3)
  182.           remax=0.d0
  183. *  au cas ou la relation soit normalisee a autre chose que 1:
  184.           do iv=2,re(/1)
  185.            remax=max(remax,abs(re(iv,1,jel)))
  186.           enddo
  187.           critl = crit*remax
  188. * deja traite?
  189.           if (ipt2.num(1,jel).eq.0) goto 110
  190.           do iv=2,re(/1)
  191. * meme simplification?
  192.            if (abs(re(iv,1,jel)).gt.critl.and.garde(iv).eq.0)goto 110
  193.            if (abs(re(iv,1,jel)).le.critl.and.garde(iv).ne.0)goto 110
  194.           enddo
  195. **        write (6,*) 'meme simplification',jel
  196. * a rajouter dans xmatr1.re
  197.           kel=kel+1
  198.           do iv=2,re(/1)
  199.            if (garde(iv).ne.0) then
  200.             ivc=garde(iv)
  201.             xmatr1.re(1,ivc,kel)=re(iv,1,jel)
  202.             xmatr1.re(ivc,1,kel)=re(iv,1,jel)
  203.            endif
  204.           enddo
  205.           do in=1,num(/1)
  206.             if (eleg(in).ne.0) then
  207.               ipt1.num(eleg(in),kel)=num(in,jel)
  208.             endif
  209.           enddo
  210. * annuler l'element
  211.           ipt2.num(1,jel)=0
  212.  110      continue
  213.           nelrig=kel
  214.           segadj xmatr1
  215.           nbelem=kel
  216.           segadj ipt1
  217.           nbsimp=nbsimp+kel
  218.          endif
  219.  100    continue
  220.          segsup trav
  221. * fin du do iel
  222. *  maintenant, on retasse ri2.irigel(ir)
  223.         if (ipt2.ne.meleme) then
  224.             nligrp=re(/1)
  225.             nligrd=re(/2)
  226.             nelrig=re(/3)-nbsimp
  227.             segini,xmatr2
  228.             ri2.irigel(4,ir)=xmatr2
  229.           nbelem=0
  230.           nbnn=ipt2.num(/1)
  231.           nbsous=0
  232.           nbref=0
  233.           do iel=1,ipt2.num(/2)
  234.            if (ipt2.num(1,iel).ne.0) then
  235.             nbelem=nbelem+1
  236.             do in=1,ipt2.num(/1)
  237.               ipt2.num(in,nbelem)=num(in,iel)
  238.             enddo
  239.             do ip=1,xmatr2.re(/2)
  240.             do id=1,xmatr2.re(/1)
  241.               xmatr2.re(id,ip,nbelem)=re(id,ip,iel)
  242.             enddo
  243.             enddo
  244.            endif
  245.           enddo
  246.           if (nbelem.ne.nelrig) call erreur(5)
  247.           segadj ipt2
  248.         endif
  249.         nbst = nbst + nbsimp
  250. *       write(6,*) ' regroupement de relation ',nbsimp
  251.  200  continue
  252. *  plus qu'a tout assembler dans le resultat.
  253. *  au passage on elimine les matrices vides
  254. *  mais avant, si on n'a rien eu a faire, on sort la matrice initiale qui est toujours dans mrigid
  255.       if (.not.travail) return
  256.  
  257.       nrige2=ri2.irigel(/2)
  258.       if (iimpi.eq.7) write(6,*) ' nombre de relations simplifiees ',
  259.      >      nbst
  260. **    write(6,*) ' simplifications retenues annulees ',nrige1,nbann
  261.       nrigel=nrige1+nrige2
  262.       segadj ri2
  263.       nrigel=0
  264.       do ir=1,nrige2
  265.         meleme=ri2.irigel(1,ir)
  266.         segact meleme
  267.         if (num(/2).ne.0) then
  268.           nrigel=nrigel+1
  269.           do is=1,ri2.irigel(/1)
  270.             ri2.irigel(is,nrigel)=ri2.irigel(is,ir)
  271.           enddo
  272.           ri2.coerig(nrigel)=ri2.coerig(ir)
  273.         endif
  274.       enddo
  275.       do ir=1,nrige1
  276.         meleme=ri1.irigel(1,ir)
  277.         segact meleme
  278.         if (num(/2).ne.0) then
  279.           nrigel=nrigel+1
  280.           do is=1,ri2.irigel(/1)
  281.             ri2.irigel(is,nrigel)=ri1.irigel(is,ir)
  282.           enddo
  283.           ri2.coerig(nrigel)=ri1.coerig(ir)
  284.         endif
  285.       enddo
  286.       segsup ri1
  287.       segadj ri2
  288.  
  289.       mrigid = ri2
  290.  
  291.       return
  292.       end
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  
  302.  
  303.  

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