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

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