Télécharger la plaquette Cast3M


fr en

Télécharger scnd3.eso

Retour à la liste

Numérotation des lignes : 

scnd3

  1. C SCND3     SOURCE    MB234859  25/07/21    21:15:02     12330          
  2. C
  3.       subroutine scnd3(mrigid,ri1,ri4)
  4. C----------------------------------------------------------------------
  5. C Eliminer les inconnues dependantes des matrices elementaires
  6. C
  7. C Entrees :
  8. C ---------
  9. C   MRIGID : pointeur MRIGID sur les rigidites elementaires a conserver.
  10. C   ri1    : pointeur MRIGID sur les rigidites elementaires reecrites
  11. C            sous forme de dependance : indique les noeuds a eliminer.
  12. C
  13. C Sortie  :
  14. C ---------
  15. C   ri4    : pointeur MRIGID sur les rigidites elementaires a conserver
  16. C            ou ont ete eliminees les inconnues dependantes.
  17. C
  18. C Remarque : Le(s) ddl(s) dont dependantes les ddls a eliminer ne sont
  19. C            pas eliminables dans la meme passe.
  20. C----------------------------------------------------------------------
  21.       IMPLICIT INTEGER(I-N)
  22.       IMPLICIT REAL*8 (A-H,O-Z)
  23. C
  24. -INC PPARAM
  25. -INC CCOPTIO
  26. -INC CCHAMP
  27. -INC CCREEL
  28. -INC SMCOORD
  29. -INC SMRIGID
  30. -INC SMELEME
  31. C
  32.       segment icpr(nbno)
  33.       segment trav1
  34.         character*4 compd(nbincl),compdd(nbincl)
  35.       endsegment
  36.       segment trav2
  37.         integer irigd(nbnoc,nbincd),nelrigd(nbnoc,nbincd)
  38.       endsegment
  39.       segment trav3
  40.         integer cor1p(nligrp5),cor1d(nligrd5)
  41.       endsegment
  42.       pointeur trav3b.trav3
  43. C
  44.       segment lisde(2,nbd)
  45.       segment ldesc(6,ndes)
  46.       segment lrigi(nbrig)
  47.       segment lelem(3,nelem)
  48.       segment ltrav
  49.         integer iporig(ncomp),ipnouv(ncomp),iligne(nlr,ncomp)
  50.         integer ipinco(ncomp),ipnode(ncomp),innode(ncomp),nbrel(ncomp)
  51.         real*8  xcoeff(nlr,ncomp)
  52.       endsegment
  53.       pointeur ltravp.ltrav,ltravd.ltrav
  54.       segment linco(nnoeud)
  55.       segment xrela(nligrm5)
  56.       segment lincp(nligrp5)
  57.       segment lincd(nligrd5)
  58.       character*(LOCHPO) coelim,coelid
  59.       logical bvu,bnew,bnsym,bok,ldiag
  60.       integer nnhash(4)
  61.       data nnhash/4*0/
  62. C     =================================================================
  63. C     Segments trav1 (noms de composantes) et icpr (numero de noeuds)
  64. C     =================================================================
  65.       ldiag=.false.
  66.       xzref=5.d-12
  67.       nde0=8
  68.       nbdeg=0
  69.       nbincd=0
  70.       nbincl=8
  71.       segini trav1
  72.       nbno=nbpts
  73.       segini icpr
  74.       nbnoc=0
  75.       nrigto=0
  76.       nligp1=0
  77. C
  78. C     Rigidite contenant les inconnues a eliminer
  79.       segact ri1
  80.       nri1=ri1.irigel(/2)
  81.       do 10 irig=1,nri1
  82.         descr=ri1.irigel(3,irig)
  83.         segact descr
  84. C
  85. C       Nombre max d'inco en relation avec le ddl elimine
  86.         lpri=lisinc(/2)
  87.         nligp1=max(nligp1,lpri)
  88. C
  89. C       Composante a eliminer (lisdua(/2)=1)
  90.         coelim=lisdua(1)
  91.         do 50 inc=1,nbincd
  92.           if (compd(inc).eq.coelim) goto 40
  93.  50     continue
  94. C
  95. C       Ajout de la composante
  96.         nbincd=nbincd+1
  97.         if(nbincd.gt.nbincl) then
  98.           nbincl=nbincd+8
  99.           segadj trav1
  100.         endif
  101.         compd(nbincd)=coelim
  102. C
  103.  40     continue
  104. C
  105. C       Identifier les noeuds presents
  106.         meleme=ri1.irigel(1,irig)
  107.         segact meleme
  108.         do j=1,num(/2)
  109.           do k=1,num(/1)
  110.             ip=num(k,j)
  111.             if(icpr(ip).eq.0) then
  112.               nbnoc=nbnoc+1
  113.               icpr(ip)=nbnoc
  114.             endif
  115.           enddo
  116.         enddo
  117. C
  118. C       Boucle sur les composantes dont peut dependre l'inc eliminee
  119.         do 45 iligrd=1,lisinc(/2)
  120.           do 55 inc=1,nbincd
  121.             if (compd(inc).eq.lisinc(iligrd)) goto 45
  122.  55       continue
  123. C
  124. C         Ajout de la composante
  125.           nbincd=nbincd+1
  126.           inc=nbincd
  127.           if(nbincd.gt.nbincl) then
  128.             nbincl=nbincd+8
  129.             segadj trav1
  130.           endif
  131.           compd(nbincd)=lisinc(iligrd)
  132. C
  133.  45     continue
  134.  10   continue
  135. *     write (6,*) ' scnd3 - 1 ',(compd(i),i=1,nbincd)
  136. C
  137. C     Rigidite contenant les matrices elementaires a conserver
  138.       segact mrigid
  139.       nbrig=mrigid.irigel(/2)
  140.       nnoeud=0
  141.       nligrp5=0
  142.       nligrd5=0
  143.       do 110 irig=1,nbrig
  144.         descr=irigel(3,irig)
  145.         segact descr
  146.         lpri=lisinc(/2)
  147.         ldua=lisdua(/2)
  148.         nligrp5=max(nligrp5,lpri)
  149.         nligrd5=max(nligrd5,ldua)
  150.         do 140 iprim=1,lpri
  151.           do 150 ico=1,nbincd
  152.             if (compd(ico).eq.lisinc(iprim)) goto 140
  153.  150      continue
  154. C
  155. C         Ajout de la composante
  156.           nbincd=nbincd+1
  157.           ico=nbincd
  158.           if(nbincd.gt.nbincl) then
  159.             nbincl=nbincd+8
  160.             segadj trav1
  161.           endif
  162.           compd(nbincd)=lisinc(iprim)
  163. C
  164.  140    continue
  165. C
  166. C       Identifier les noeuds presents
  167.         meleme=irigel(1,irig)
  168.         segact meleme
  169.         nnoeud=max(nnoeud,num(/1))
  170.         do j=1,num(/2)
  171.           do i=1,num(/1)
  172.             ip=num(i,j)
  173.             if(icpr(ip).eq.0) then
  174.               nbnoc=nbnoc+1
  175.               icpr(ip)=nbnoc
  176.             endif
  177.           enddo
  178.         enddo
  179.  110  continue
  180.       nligrm5=max(nligrp5,nligrd5)
  181.       segini,linco,xrela,lincp,lincd,trav3
  182. *
  183. *     La matrice ou l'on a elimine les inconnues peut etre plus grosse
  184.       nligrp5=(nligrm5+nligp1)*2
  185.       nligrd5=(nligrm5+nligp1)*2
  186.       segini,trav3b
  187. *     write (6,*) ' scnd3 - 2 ',(compd(i),i=1,nbincd)
  188. C
  189. C     Noms des composantes duales : recuperer dans nomdu
  190. C     eventuellement surcharger/completer par l'utilisateur (opti inco)
  191.       do 300 i=1,nbincd
  192.         do 310 j=1,lnomdd
  193.           if (compd(i).eq.nomdd(j)) goto 320
  194.  310    continue
  195. C
  196. C       write (6,*) 'Inco. non definie ',compd(i),' ',compdd(i)
  197.         moterr(1:LOCHPO)=compd(i)
  198.         call erreur(1150)
  199.         return
  200. C
  201.  320    continue
  202.         compdd(i)=nomdu(j)
  203. **      write (6,*) ' correspondance ', compd(i), compdd(i)
  204.  300  continue
  205. *     write (6,*) ' scnd3 - 3 ',compd(/2),(compd(i),i=1,compd(/2))
  206. *     write (6,*) ' scnd3 - 3 ',compdd(/2),(compdd(i),i=1,compdd(/2))
  207. C     =================================================================
  208. C     Segment trav2 : identifier les numeros des noeuds a eliminer
  209. C     =================================================================
  210.       segini trav2
  211.       do 200 irig=1,nri1
  212.         meleme=ri1.irigel(1,irig)
  213.         descr =ri1.irigel(3,irig)
  214.         do 250 ico=1,nbincd
  215.           if (compd(ico).eq.lisdua(1)) goto 260
  216.  250    continue
  217.         call erreur(5)
  218.  260    continue
  219.         do 240 i=1,num(/2)
  220.           ip=num(noeled(1),i)
  221.           irigd(icpr(ip),ico)=irig
  222.           nelrigd(icpr(ip),ico)=i
  223.  240    continue
  224.  200  continue
  225. C
  226. C     Nombre max de relations dans lesquelles intervient un ddl donne
  227.       nlr=0
  228.       do 201 irig=1,nri1
  229.         meleme=ri1.irigel(1,irig)
  230.         descr =ri1.irigel(3,irig)
  231.         lpri=noelep(/1)
  232.         do 241 i=1,num(/2)
  233.           do 242 j=1,lpri
  234.             ip=num(noelep(j),i)
  235.             do 252 ico=1,nbincd
  236.               if (compd(ico).eq.lisinc(j)) goto 262
  237.  252        continue
  238.             call erreur(5)
  239.  262        continue
  240.             if (nelrigd(icpr(ip),ico).gt.0) goto 242
  241.             nelrigd(icpr(ip),ico)=nelrigd(icpr(ip),ico)-1
  242.             nlr=max(nlr,abs(nelrigd(icpr(ip),ico)))
  243.  242      continue
  244.  241    continue
  245.  201  continue
  246. C     =================================================================
  247. C     Eliminer des matrices elementaires les inc qui doivent l'etre
  248. C     =================================================================
  249.       ndes=nde0
  250.       nbd=nde0
  251.       segini,lrigi,ldesc,lisde
  252.       krigel=0
  253.       kdescr=0
  254. C     ----------------------------
  255. C     Boucle sur les sous-matrices
  256.       do 1000 irig=1,nbrig
  257.         ipt4  =mrigid.irigel(1,irig)
  258.         des4  =mrigid.irigel(3,irig)
  259.         xmatr4=mrigid.irigel(4,irig)
  260.         segact,ipt4,des4,xmatr4
  261. C
  262. C       La matrice est-elle symetrique
  263.         bnsym=(xmatr4.symre.ne.0)
  264. C
  265.         nelrig=xmatr4.re(/3)
  266.         nligrp4=des4.lisinc(/2)
  267.         nligrd4=des4.lisdua(/2)
  268. C
  269. C       Correspondance lisinc, lisdua, numero ico dans trav1
  270.         do 410 iprim=1,nligrp4
  271.           cor1p(iprim)=0
  272.           do 420 ico=1,nbincd
  273.             if (compd(ico).eq.des4.lisinc(iprim)) then
  274.               cor1p(iprim)=ico
  275.               goto 410
  276.             endif
  277.  420      continue
  278. **        write(6,*) 'cor1p pas de correspondance'
  279. **        call erreur(5)
  280.  410    continue
  281. C
  282.         do 510 idua=1,nligrd4
  283.           cor1d(idua)=0
  284.           do 520 ico=1,nbincd
  285.             if (compdd(ico).eq.des4.lisdua(idua)) then
  286.               cor1d(idua)=ico
  287.               goto 510
  288.             endif
  289.  520      continue
  290. **        call erreur(5)
  291.  510    continue
  292. C
  293.         nnoeu=ipt4.num(/1)
  294.         nelem=ipt4.num(/2)
  295.         segini,lelem
  296.         lrigi(irig)=lelem
  297.         bvu=.false.
  298. C
  299. C       ===============================================================
  300. C       Boucle sur les elements
  301. C       ===============================================================
  302.         do 1040 iel=1,nelem
  303. C
  304. C         Reinitialiser linco, lincp et lincd
  305.           do inn=1,nnoeu
  306.             linco(inn)=0
  307.           enddo
  308.           do inn=1,nligrp4
  309.             lincp(inn)=0
  310.           enddo
  311.           if (bnsym) then
  312.             do inn=1,nligrd4
  313.               lincd(inn)=0
  314.             enddo
  315.           endif
  316.           ipass=1
  317. C
  318.           icpt=0
  319.           kelimp=0
  320.           ncor=0
  321.           do 1009 iligri=1,nligrp4
  322.             incp=cor1p(iligri)
  323.             ip=ipt4.num(des4.noelep(iligri),iel)
  324.             inn=des4.noelep(iligri)
  325.             linco(inn)=linco(inn)+1
  326. C
  327. C           Ce ddl est-il a eliminer
  328.             if (nelrigd(icpr(ip),incp).le.0) then
  329.               icpt=icpt+1
  330.               lincp(iligri)=icpt
  331.               goto 1009
  332.             endif
  333. C
  334.             linco(inn)=linco(inn)-1
  335.             irig1=irigd(icpr(ip),incp)
  336.             i1   =nelrigd(icpr(ip),incp)
  337.             des1 =ri1.irigel(3,irig1)
  338.             ncor=ncor+des1.lisinc(/2)
  339.             kelimp=kelimp+1
  340.  1009     continue
  341. C
  342. C         Matrice non symetrique : verifier si la duale de l'inconnue
  343. C         eliminee n'est pas presente
  344.           kelimd=kelimp
  345.           if (bnsym) then
  346.             icpt=0
  347.             kelimd=0
  348.             do 1010 iligrd=1,nligrd4
  349.               incd=cor1d(iligrd)
  350.               ip=ipt4.num(des4.noeled(iligrd),iel)
  351.               inn=des4.noeled(iligrd)
  352.               linco(inn)=linco(inn)+1
  353. C
  354. C             L'inconnue primale associee est-il eliminee
  355.               if (nelrigd(icpr(ip),incd).le.0) then
  356.                 icpt=icpt+1
  357.                 lincd(iligrd)=icpt
  358.                 goto 1010
  359.               endif
  360. C
  361.               linco(inn)=linco(inn)-1
  362.               irig1=irigd(icpr(ip),incd)
  363.               i1   =nelrigd(icpr(ip),incd)
  364.               des1 =ri1.irigel(3,irig1)
  365.               ncor=ncor+des1.lisinc(/2)
  366.               kelimd=kelimd+1
  367.  1010       continue
  368.           endif
  369. C
  370. C         -------------------------------------------------------------
  371. C         Pas de modification a apporter, passer a l'element suivant
  372.           if ((kelimp+kelimd).eq.0) then
  373.             if (.not.bvu) then
  374. C
  375.               nnhash(1)=nnhash(1)+1
  376.               call deshas(nbincd,trav3,nnoeu,des4,idhash)
  377. C
  378. C             Tentative de reutilisation d'un descripteur
  379.               mdes=0
  380.               do 2030 ides=kdescr,1,-1
  381.                 idhtmp=lisde(2,ides)
  382.                 if (idhash.ne.idhtmp) then
  383.                   nnhash(2)=nnhash(2)+1
  384.                   goto 2030
  385.                 endif
  386. C
  387.                 des3=lisde(1,ides)
  388.                 if (nligrd4.ne.des3.noeled(/1)) goto 2029
  389.                 if (nligrp4.ne.des3.noelep(/1)) goto 2029
  390.                 do id1=1,nligrd4
  391.                   if (des4.noeled(id1).ne.des3.noeled(id1)) goto 2029
  392.                   if (des4.lisdua(id1).ne.des3.lisdua(id1)) goto 2029
  393.                 enddo
  394.                 do id1=1,nligrp4
  395.                   if (des4.noelep(id1).ne.des3.noelep(id1)) goto 2029
  396.                   if (des4.lisinc(id1).ne.des3.lisinc(id1)) goto 2029
  397.                 enddo
  398. C
  399. C               Reutilisation d'un descripteur existant
  400.                 nnhash(4)=nnhash(4)+1
  401.                 mdes=ides
  402.                 goto 2032
  403. C
  404.  2029           continue
  405.                 nnhash(3)=nnhash(3)+1
  406.  2030         continue
  407. C
  408. C             Nouveau descripteur
  409.               kdescr=kdescr+1
  410.               mdes=kdescr
  411.               if (lisde(/2).lt.kdescr) then
  412.                 nbd=kdescr+nde0
  413.                 segadj,lisde
  414.               endif
  415.               lisde(1,mdes)=des4
  416.               lisde(2,mdes)=idhash
  417.               goto 2033
  418. C
  419.  2032         continue
  420. C
  421. C             Rigidite avec memes descripteur et caracteristiques?
  422.               iptd=lisde(1,mdes)
  423.               do 2031 irie=krigel,1,-1
  424. C
  425.                 if (ldesc(1,irie).ne.iptd)                  goto 2031
  426.                 if (ldesc(3,irie).ne.mrigid.irigel(2,irig)) goto 2031
  427.                 if (ldesc(4,irie).ne.mrigid.irigel(5,irig)) goto 2031
  428.                 if (ldesc(5,irie).ne.mrigid.irigel(6,irig)) goto 2031
  429.                 if (ldesc(6,irie).ne.mrigid.irigel(7,irig)) goto 2031
  430. C
  431. C               Nouvel element pour la rigidite
  432.                 mrig=irie
  433.                 goto 2034
  434. C
  435.  2031         continue
  436. C
  437.  2033         continue
  438. C
  439. C             Nouvelle rigidite
  440.               krigel=krigel+1
  441.               mrig=krigel
  442.               if (ldesc(/2).lt.mrig) then
  443.                 ndes=mrig+nde0
  444.                 segadj,ldesc
  445.               endif
  446.               ldesc(1,mrig)=lisde(1,mdes)
  447.               ldesc(2,mrig)=0
  448.               ldesc(3,mrig)=mrigid.irigel(2,irig)
  449.               ldesc(4,mrig)=mrigid.irigel(5,irig)
  450.               ldesc(5,mrig)=mrigid.irigel(6,irig)
  451.               ldesc(6,mrig)=mrigid.irigel(7,irig)
  452. C
  453.  2034         continue
  454.               bvu=.true.
  455.             endif
  456. C
  457.             lelem(1,iel)=ipt4
  458.             lelem(2,iel)=mrig
  459.             lelem(3,iel)=-1*xmatr4
  460.             ldesc(2,mrig)=ldesc(2,mrig)+1
  461. C
  462.             goto 1040
  463.           endif
  464. C
  465. C         -------------------------------------------------------------
  466. C         Modifications de la matrice elementraire
  467. C
  468. C         Les noeuds sont-ils toujours presents
  469.           nbn1=0
  470.           do ind=1,nnoeu
  471.             if (linco(ind).gt.0) then
  472.               nbn1=nbn1+1
  473.               linco(ind)=nbn1
  474.             endif
  475.           enddo
  476. C
  477.           ncomp=nligrp4+ncor
  478.           segini,ltravp
  479.           ltrav=ltravp
  480.           nbcomp=nligrp4
  481.           ncompp=ncomp
  482.           kelim=kelimp
  483.           kno=0
  484. C
  485.  1111     continue
  486. C
  487.           if (ipass.eq.2) then
  488.             ncomp=nligrd4+ncor
  489.             segini,ltravd
  490.             ltrav=ltravd
  491.             nbcomp=nligrd4
  492.             ncompd=ncomp
  493.             kelim=kelimd
  494.           endif
  495. C
  496.           do icc=1,nbcomp
  497.             ltrav.iporig(icc)=icc
  498.             if (ipass.eq.1) then
  499.               ipn=des4.noelep(icc)
  500.               incp=cor1p(icc)
  501.             else
  502.               ipn=des4.noeled(icc)
  503.               incp=cor1d(icc)
  504.             endif
  505.             if (linco(ipn).ne.0) then
  506.               ip=ipt4.num(ipn,iel)
  507.               ltrav.innode(icc)=ip
  508.               ltrav.ipnode(icc)=linco(ipn)
  509.               ltrav.ipinco(icc)=incp
  510. C             Ce ddl est-il a eliminer
  511.               if (nelrigd(icpr(ip),incp).le.0) then
  512.                 if (ipass.eq.1) ltrav.ipnouv(icc)=lincp(icc)
  513.                 if (ipass.eq.2) ltrav.ipnouv(icc)=lincd(icc)
  514.               endif
  515.             endif
  516.           enddo
  517. C
  518.           kligrp4=0
  519.           kcomp=nbcomp
  520.           do 1011 iligri=1,nbcomp
  521. C
  522.             if (ipass.eq.1) then
  523.               ipn=des4.noelep(iligri)
  524.               incp=cor1p(iligri)
  525.             else
  526.               ipn=des4.noeled(iligri)
  527.               incp=cor1d(iligri)
  528.             endif
  529.             ip=ipt4.num(ipn,iel)
  530. C
  531. C           Ce ddl est-il a eliminer
  532.             if (nelrigd(icpr(ip),incp).le.0) goto 1011
  533. C
  534. C           Le(s) ddl(s) lie(s) au ddl elimine exite(nt)-il(s)
  535.             irig1=irigd(icpr(ip),incp)
  536.             i1    =nelrigd(icpr(ip),incp)
  537.             ipt1  =ri1.irigel(1,irig1)
  538.             des1  =ri1.irigel(3,irig1)
  539.             xmatr1=ri1.irigel(4,irig1)
  540.             segact,xmatr1
  541. C
  542.             xvalm=0.d0
  543.             do icoel=1,des1.lisinc(/2)
  544.               xvalm=max(xvalm,abs(xmatr1.re(1,icoel,i1)))
  545.             enddo
  546.             xvalm=xvalm*xzref
  547. C
  548. C           -----------------------------------------------------------
  549.             do 1012 icoel=1,des1.lisinc(/2)
  550. C
  551. C             Coefficient de la relation non significatif : iterer
  552.               if (abs(xmatr1.re(1,icoel,i1)).le.xvalm) goto 1012
  553. C
  554.               ipl=ipt1.num(des1.noelep(icoel),i1)
  555. C
  556.               do 1420 ico=1,nbincd
  557.                 if (compd(ico).eq.des1.lisinc(icoel)) then
  558.                   iprim=ico
  559.                   goto 1410
  560.                 endif
  561.  1420         continue
  562.  1410         continue
  563. C
  564. C             Noeuds deja presents dans l'element
  565.               ino=0
  566.               ipo=0
  567.               do 1013 idq=1,ncomp
  568.                 if (ltrav.innode(idq).ne.ipl) goto 1013
  569.                 ino=ltrav.ipnode(idq)
  570.                 if (ltrav.ipinco(idq).ne.iprim) goto 1013
  571.                 ipo=idq
  572.                 goto 1014
  573.  1013         continue
  574. C
  575. C             Il faut ajouter un noeud au segment descripteur
  576.               if (ino.eq.0) then
  577.                 kno=kno+1
  578.                 ino=nbn1+kno
  579.               endif
  580. C
  581. C             Il faut completer le segment descripteur
  582.               kligrp4=kligrp4+1
  583. C
  584.  1014         continue
  585. C
  586.               if (ipo.eq.0) then
  587. C               Conserver les informations du nouveau ddl
  588.                 kcomp=kcomp+1
  589.                 ltrav.innode(kcomp)=ipl
  590.                 ltrav.ipnode(kcomp)=ino
  591.                 ltrav.ipinco(kcomp)=iprim
  592.                 ltrav.ipnouv(kcomp)=nbcomp-kelim+kligrp4
  593.                 ipo=kcomp
  594.               endif
  595. C
  596.               ltrav.nbrel(ipo)=ltrav.nbrel(ipo)+1
  597.               ltrav.xcoeff(ltrav.nbrel(ipo),ipo)=xmatr1.re(1,icoel,i1)
  598.               ltrav.iligne(ltrav.nbrel(ipo),ipo)=ltrav.iporig(iligri)
  599. C
  600.  1012       continue
  601. C           -----------------------------------------------------------
  602. C
  603.  1011     continue
  604. C
  605. C         Cas d'une matrice vide
  606.           if (ipass.eq.1) then
  607.             nligrp=nligrp4-kelimp+kligrp4
  608.             if (nligrp.eq.0) goto 1040
  609. C
  610. C           Matrice non-symetrique : refaire le travail sur les duales
  611.             if (bnsym) then
  612.               ipass=2
  613.               goto 1111
  614.             endif
  615. C
  616. C           Matrice symetrique : infos duals identiques
  617.             ltravd=ltravp
  618.             ncompd=ncompp
  619.             nligrd=nligrp
  620.           else
  621.             nligrd=nligrd4-kelimd+kligrp4
  622.             if (nligrd.eq.0) goto 1040
  623.           endif
  624. C
  625. C         Matrice de relation vide car tout a ete elimine
  626.           ityp4=28
  627.           if ((ipt4.itypel.eq.49).or.(ipt4.itypel.eq.22)) then
  628.             iid=1
  629.             if (ipt4.itypel.eq.49) iid=2
  630.             if (nligrp.le.iid) then
  631.               nbdeg=nbdeg+1
  632.               goto 1040
  633.             endif
  634.             ityp4=ipt4.itypel
  635.           endif
  636. C
  637. C         Creation du nouveau maillage, descripteur, matrice
  638. C
  639.           nbnn=nbn1+kno
  640.           nbelem=1
  641.           nbref=0
  642.           nbsous=0
  643.           segini,ipt2
  644.           ipt2.itypel=ityp4
  645. C
  646.           segini,des2
  647. C
  648.           nelrig=1
  649.           segini,xmatr2
  650.           xmatr2.symre=xmatr4.symre
  651. C
  652. C         Renseigner le maillage et le descripteur
  653.           do 1021 ic=1,ncompp
  654.             inew=ltravp.ipnouv(ic)
  655.             if (inew.ne.0) then
  656.               des2.noelep(inew)=ltravp.ipnode(ic)
  657.               des2.lisinc(inew)=compd(ltravp.ipinco(ic))
  658.               trav3b.cor1p(inew)=ltravp.ipinco(ic)
  659.               ipt2.num(ltravp.ipnode(ic),1)=ltravp.innode(ic)
  660.             endif
  661.  1021     continue
  662.           do 1022 ic=1,ncompd
  663.             inew=ltravd.ipnouv(ic)
  664.             if (inew.ne.0) then
  665.               des2.noeled(inew)=ltravd.ipnode(ic)
  666.               des2.lisdua(inew)=compdd(ltravd.ipinco(ic))
  667.               trav3b.cor1d(inew)=ltravd.ipinco(ic)
  668.               ipt2.num(ltravd.ipnode(ic),1)=ltravd.innode(ic)
  669.             endif
  670.  1022     continue
  671. C
  672. C         Renseigner la matrice
  673.           do 1024 ic=1,ncompp
  674.             iori=ltravp.iporig(ic)
  675.             inew=ltravp.ipnouv(ic)
  676.             if (inew.ne.0) then
  677. C
  678. C             Reinitialiser xrela
  679.               if (iori.ne.0) then
  680.                 do il=1,nligrd4
  681.                   xrela(il)=xmatr4.re(il,iori,iel)
  682.                 enddo
  683.               else
  684.                 do il=1,nligrd4
  685.                   xrela(il)=0.D0
  686.                 enddo
  687.               endif
  688. C
  689.               do inr=1,ltravp.nbrel(ic)
  690.                 xtmp=ltravp.xcoeff(inr,ic)
  691.                 ilig=ltravp.iligne(inr,ic)
  692.                 do il=1,nligrd4
  693.                   xval=xmatr4.re(il,ilig,iel)*xtmp
  694.                   xrela(il)=xrela(il)+xval
  695.                 enddo
  696.               enddo
  697. C
  698.               do ie=1,ncompd
  699.                 knew=ltravd.ipnouv(ie)
  700.                 if (knew.ne.0) then
  701.                   xval=0.D0
  702.                   if (ie.le.nligrd4) xval=xrela(ie)
  703.                   do inr=1,ltravd.nbrel(ie)
  704.                     xtmp=ltravd.xcoeff(inr,ie)
  705.                     ilig=ltravd.iligne(inr,ie)
  706.                     xval=xval+xtmp*xrela(ilig)
  707.                   enddo
  708.                   xmatr2.re(knew,inew,1)=xmatr2.re(knew,inew,1)+xval
  709.                 endif
  710.               enddo
  711. C
  712.             endif
  713.  1024     continue
  714. C
  715. C         Suppression des termes negligeables de la relation
  716.           if (ipt2.itypel.eq.49.or.ipt2.itypel.eq.22) then
  717. C
  718.             renorm=0.d0
  719.             do ir=1,nligrp4
  720.               do jr=1,nligrd4
  721.                 re1=xmatr4.re(ir,jr,1)
  722.                 renorm=renorm+abs(re1)
  723.               enddo
  724.             enddo
  725.             renorm=max(xpetit,renorm*xzref)
  726. C
  727.             iid=2
  728.             if (ipt2.itypel.eq.49) iid=3
  729. C
  730.             idec=0
  731.             do ii=iid,nligrd
  732.               if (abs(xmatr2.re(ii,1,1)).lt.renorm.and.
  733.      >            abs(xmatr2.re(1,ii,1)).lt.renorm) then
  734.                 idec=idec+1
  735.               elseif(idec.ne.0) then
  736.                 xmatr2.re(ii-idec,1,1)=xmatr2.re(ii,1,1)
  737.                 if(iid.eq.3) xmatr2.re(ii-idec,2,1)=xmatr2.re(ii,2,1)
  738.                 xmatr2.re(1,ii-idec,1)=xmatr2.re(1,ii,1)
  739.                 if(iid.eq.3) xmatr2.re(2,ii-idec,1)=xmatr2.re(2,ii,1)
  740.                 des2.lisinc(ii-idec)=des2.lisinc(ii)
  741.                 des2.lisdua(ii-idec)=des2.lisdua(ii)
  742.                 des2.noelep(ii-idec)=des2.noelep(ii)
  743.                 des2.noeled(ii-idec)=des2.noeled(ii)
  744.                 trav3b.cor1p(ii-idec)=trav3b.cor1p(ii)
  745.                 trav3b.cor1d(ii-idec)=trav3b.cor1d(ii)
  746.               endif
  747.             enddo
  748. C
  749.             nligrd=nligrd-idec
  750.             if(nligrd.lt.iid) then
  751. **        write (6,*) 'scnd3 matrice supprimee LX',nligrd,itypel,meleme
  752.               nbdeg=nbdeg+1
  753.               segsup,ipt2,des2,xmatr2
  754.               goto 1040
  755.             endif
  756. C
  757. C           Ajuster les segments et retirer les noeuds non references
  758.             if (idec.ne.0) then
  759. C
  760.               nligrp=nligrp-idec
  761.               nelrig=1
  762.               segadj,des2,xmatr2
  763. C
  764.               do ii=1,nligrp
  765.                 inp=des2.noelep(ii)
  766.                 ipt2.num(inp,1)=-1*abs(ipt2.num(inp,1))
  767.               enddo
  768.               do ii=1,nligrd
  769.                 ind=des2.noeled(ii)
  770.                 ipt2.num(ind,1)=-1*abs(ipt2.num(ind,1))
  771.               enddo
  772. C
  773.               iden=0
  774.               do ii=1,nbnn
  775.                 if (ipt2.num(ii,1).gt.0) then
  776.                   iden=iden+1
  777.                 else
  778.                   if (iden.ne.0) then
  779.                     do no=1,nligrp
  780.                       if (des2.noelep(no).eq.ii) then
  781.                         des2.noelep(no)=des2.noelep(no)-iden
  782.                       endif
  783.                     enddo
  784.                     do no=1,nligrd
  785.                       if (des2.noeled(no).eq.ii) then
  786.                         des2.noeled(no)=des2.noeled(no)-iden
  787.                       endif
  788.                     enddo
  789.                   endif
  790.                   ipt2.num(ii-iden,1)=abs(ipt2.num(ii,1))
  791.                 endif
  792.               enddo
  793.               nbnn=nbnn-iden
  794.               nbelem=1
  795.               nbsous=0
  796.               nbref=0
  797.               segadj,ipt2
  798.             endif
  799. C
  800.           endif
  801. C
  802. C         Matrice symetrique : forcer la symetrie de xmatri
  803.           if(xmatr2.symre.eq.0) then
  804.             ninc=32
  805.             do ib=1,nligrd,ninc
  806.               do jr=1,nligrd
  807.                 do ir=ib,min(jr-1,ib+ninc-1)
  808.                   re1=xmatr2.re(ir,jr,1)
  809.                   re2=xmatr2.re(jr,ir,1)
  810.                   rem=(re1+re2)/2.d0
  811.                   xmatr2.re(ir,jr,1)=rem
  812.                   xmatr2.re(jr,ir,1)=rem
  813.                 enddo
  814.               enddo
  815.             enddo
  816.           endif
  817. C
  818. C         Supprimer les segments de travail
  819.           segsup,ltravp
  820.           if (bnsym) segsup,ltravd
  821. C
  822. C         Tentative de reutilisation d'un descripteur
  823.           nnhash(1)=nnhash(1)+1
  824.           call deshas(nbincd,trav3b,nbnn,des2,idhash)
  825.           jdes=0
  826.           do 1030 ides=kdescr,1,-1
  827.             idhtmp=lisde(2,ides)
  828.             if (idhash.ne.idhtmp) then
  829.               nnhash(2)=nnhash(2)+1
  830.               goto 1030
  831.             endif
  832. C
  833.             des3=lisde(1,ides)
  834.             if (nligrd.ne.des3.noeled(/1)) goto 1029
  835.             if (nligrp.ne.des3.noelep(/1)) goto 1029
  836.             do id1=1,nligrd
  837.               if (des2.noeled(id1).ne.des3.noeled(id1)) goto 1029
  838.               if (des2.lisdua(id1).ne.des3.lisdua(id1)) goto 1029
  839.             enddo
  840.             do id1=1,nligrp
  841.               if (des2.noelep(id1).ne.des3.noelep(id1)) goto 1029
  842.               if (des2.lisinc(id1).ne.des3.lisinc(id1)) goto 1029
  843.             enddo
  844. C
  845. C           Descripteur existant
  846.             nnhash(4)=nnhash(4)+1
  847.             jdes=ides
  848.             segsup,des2
  849.             goto 1032
  850. C
  851.  1029       continue
  852.             nnhash(3)=nnhash(3)+1
  853.  1030     continue
  854. C
  855. C         Nouveau descripteur
  856.           nnhash(3)=nnhash(3)+1
  857.           kdescr=kdescr+1
  858.           jdes=kdescr
  859.           if (lisde(/2).lt.kdescr) then
  860.             nbd=kdescr+nde0
  861.             segadj,lisde
  862.           endif
  863.           lisde(1,jdes)=des2
  864.           lisde(2,jdes)=idhash
  865.           goto 1033
  866. C
  867.  1032     continue
  868. C
  869. C         Rigidite avec memes descripteur et caracteristiques?
  870.           iptd=lisde(1,jdes)
  871.           do 1031 irie=krigel,1,-1
  872. C
  873.             if (ldesc(1,irie).ne.iptd)                  goto 1031
  874.             if (ldesc(3,irie).ne.mrigid.irigel(2,irig)) goto 1031
  875.             if (ldesc(4,irie).ne.mrigid.irigel(5,irig)) goto 1031
  876.             if (ldesc(5,irie).ne.mrigid.irigel(6,irig)) goto 1031
  877.             if (ldesc(6,irie).ne.mrigid.irigel(7,irig)) goto 1031
  878. C
  879.             jrig=irie
  880.             goto 1034
  881. C
  882.  1031     continue
  883. C
  884.  1033     continue
  885. C
  886. C         Nouvelle rigidite
  887.           krigel=krigel+1
  888.           jrig=krigel
  889.           if (ldesc(/2).lt.jrig) then
  890.             ndes=jrig+nde0
  891.             segadj,ldesc
  892.           endif
  893.           ldesc(1,jrig)=lisde(1,jdes)
  894.           ldesc(2,jrig)=0
  895.           ldesc(3,jrig)=mrigid.irigel(2,irig)
  896.           ldesc(4,jrig)=mrigid.irigel(5,irig)
  897.           ldesc(5,jrig)=mrigid.irigel(6,irig)
  898.           ldesc(6,jrig)=mrigid.irigel(7,irig)
  899. C
  900.  1034     continue
  901. C
  902.           lelem(1,iel)=ipt2
  903.           lelem(2,iel)=jrig
  904.           lelem(3,iel)=xmatr2
  905.           ldesc(2,jrig)=ldesc(2,jrig)+1
  906. C
  907.  1040   continue
  908. C       ===============================================================
  909. C
  910.  1000 continue
  911. C     Fin de boucle sur les sous-matrices
  912. C     -----------------------------------
  913.       segsup,linco,xrela,lincp,lincd,trav3,trav3b
  914.       if (ldiag) then
  915.         write(ioimp,*)'Nb calculs hash          =',nnhash(1)
  916.         write(ioimp,*)'Nb comparaisons raccourcies hash=',nnhash(2)
  917.         write(ioimp,*)'Nb comparaisons raccourcies nohash=',nnhash(3)
  918.         write(ioimp,*)'Nb comparaisons completes=',nnhash(4)
  919.         write(ioimp,*)'Nb sous rigidite initiale=',mrigid.irigel(/2)
  920.         write(ioimp,*)'Nb sous rigidite finale  =',krigel
  921.         write(ioimp,*)'Nb de descripteurs       =',kdescr
  922.       endif
  923. C     =================================================================
  924. C     Construire la rigidite ou les inconnues ont ete eliminees
  925. C     =================================================================
  926.       nrigel=krigel
  927.       segini,ri4
  928.       ri4.mtymat='TEMPORAI'
  929.       ri4.iforig=mrigid.iforig
  930.       do ibr=1,nbrig
  931.         lelem=lrigi(ibr)
  932.         ipt4 =mrigid.irigel(1,ibr)
  933. C
  934. C       Renseigner les matrices elementaires
  935.         do 1042 iel=1,ipt4.num(/2)
  936.           ipt6=lelem(1,iel)
  937.           if (ipt6.eq.0) goto 1042
  938.           kdes=lelem(2,iel)
  939.           xmatr6=lelem(3,iel)
  940. C
  941.           itmp=1
  942.           bnew=.true.
  943.           if (xmatr6.lt.0) then
  944.             bnew=.false.
  945.             xmatr6=abs(xmatr6)
  946.             itmp=iel
  947.           endif
  948. C
  949.           des3=ldesc(1,kdes)
  950.           nelt=ldesc(2,kdes)
  951.           ipt5=ri4.irigel(1,kdes)
  952. C
  953. C         Creation de la sous-matrice
  954.           if (ipt5.eq.0) then
  955.             ldesc(2,kdes)=1
  956. C
  957.             nligrp=des3.noelep(/1)
  958.             nligrd=des3.noeled(/1)
  959.             nelrig=nelt
  960. C
  961.             bok=(nelrig.eq.xmatr6.re(/3).and.bnew)
  962.             if (bok) then
  963.               xmatr5=xmatr6
  964.               ipt5=ipt6
  965.               goto 1041
  966.             endif
  967. C
  968.             segini,xmatr5
  969.             xmatr5.symre=xmatr6.symre
  970. C
  971.             nbnn=ipt6.num(/1)
  972.             nbelem=nelt
  973.             nbref=0
  974.             nbsous=0
  975.             segini,ipt5
  976.             if (bnew) then
  977.               ipt5.itypel=28
  978.               if ((ipt6.itypel.eq.49).or.(ipt6.itypel.eq.22)) then
  979.                 ipt5.itypel=ipt6.itypel
  980.               endif
  981.             else
  982.               ipt5.itypel=ipt6.itypel
  983.             endif
  984. C
  985.  1041       continue
  986. C
  987.             ri4.coerig(kdes)=mrigid.coerig(ibr)
  988.             ri4.irigel(1,kdes)=ipt5
  989.             ri4.irigel(2,kdes)=mrigid.irigel(2,ibr)
  990.             ri4.irigel(3,kdes)=des3
  991.             ri4.irigel(4,kdes)=xmatr5
  992.             ri4.irigel(5,kdes)=mrigid.irigel(5,ibr)
  993.             ri4.irigel(6,kdes)=mrigid.irigel(6,ibr)
  994.             ri4.irigel(7,kdes)=mrigid.irigel(7,ibr)
  995.             ri4.irigel(8,kdes)=mrigid.irigel(8,ibr)
  996. C
  997.             if (bok) goto 1042
  998.           endif
  999. C
  1000.           ielmt=ldesc(2,kdes)
  1001.           do inn=1,ipt6.num(/1)
  1002.             ipt5.num(inn,ielmt)=ipt6.num(inn,itmp)
  1003.           enddo
  1004. C
  1005.           xcoer=mrigid.coerig(ibr)/ri4.coerig(kdes)
  1006.           xmatr5=ri4.irigel(4,kdes)
  1007.           do ilc=1,xmatr5.re(/2)
  1008.             do ili=1,xmatr5.re(/1)
  1009.               xmatr5.re(ili,ilc,ielmt)=xcoer*xmatr6.re(ili,ilc,itmp)
  1010.             enddo
  1011.           enddo
  1012. C
  1013.           ldesc(2,kdes)=ldesc(2,kdes)+1
  1014.           if (bnew) segsup,ipt6,xmatr6
  1015.  1042   continue
  1016.         segsup,lelem
  1017.       enddo
  1018.       segsup,icpr,trav1,trav2,lrigi,ldesc
  1019. C
  1020. CCC   WRITE(*,*) 'MRIGID'
  1021. CCC   call prrigi(MRIGID,0)
  1022. CCC   WRITE(*,*) 'RI1'
  1023. CCC   call prrigi(ri1,0)
  1024. CCC   WRITE(*,*) 'RI4'
  1025. CCC   call prrigi(ri4,0)
  1026. C     =================================================================
  1027. C     Desactivation generale
  1028. C     =================================================================
  1029.       segact mrigid
  1030.       do irig=1,nbrig
  1031.         xmatri=irigel(4,irig)
  1032.         meleme=irigel(1,irig)
  1033.         descr=irigel(3,irig)
  1034.         segdes xmatri,descr
  1035.       enddo
  1036.       segact ri4
  1037.       do irig=1,ri4.irigel(/2)
  1038.         xmatri=ri4.irigel(4,irig)
  1039.         meleme=ri4.irigel(1,irig)
  1040.         descr=ri4.irigel(3,irig)
  1041.         segdes xmatri,descr
  1042.       enddo
  1043.       segact ri1
  1044.       do irig=1,nri1
  1045.         xmatri=ri1.irigel(4,irig)
  1046.         meleme=ri1.irigel(1,irig)
  1047.         descr=ri1.irigel(3,irig)
  1048.         segdes xmatri,descr
  1049.       enddo
  1050.       if(iimpi.ne.0.and.nbdeg.ne.0)
  1051.      >  write (6,*) ' nombre de relations degenerees'//
  1052.      >  ' apres elimination ',nbdeg
  1053. *     write(6,*) 'nbincl nligrp5 nnoeud ndes nlr',
  1054. *    >            nbincl,nligrp5,nnoeud,ndes,nlr
  1055.       end
  1056.  
  1057.  
  1058.  

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