Télécharger la plaquette Cast3M


fr en

Télécharger exincr.eso

Retour à la liste

Numérotation des lignes : 

exincr

  1. C EXINCR    SOURCE    FANDEUR   22/03/01    21:15:04     11301          
  2.       SUBROUTINE EXINCR(RI1,MLMOT1,MLMOT2,RI2)
  3.       IMPLICIT REAL*8 (A-H,O-Z)
  4.       IMPLICIT INTEGER (I-N)
  5. C***********************************************************************
  6. C NOM         : EXINCR
  7. C DESCRIPTION : Extrait d'une rigidité la sous-matrice
  8. C               d'inconnues primales et duales celles données
  9. C               en argument CH*4
  10. * Dans le cas d'une relation, il serait prudent
  11. * de créer un nouveau multiplicateur mais alors si
  12. * on recombine les matrices extraites, on ne reobtient
  13. * pas l'originale. D'autre part, si on ne crée pas
  14. * de nouveau multiplicateur reso buggera sur la
  15. * matrice recombinée (présence d'un même multiplicateur
  16. * dans deux matrices différentes) On choisit de ne pas
  17. * créer de multiplicateur ici.
  18. C
  19. C
  20. C LANGAGE     : ESOPE
  21. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  22. C               mél : gounand@semt2.smts.cea.fr
  23. C***********************************************************************
  24. C APPELES          :
  25. C APPELES (E/S)    :
  26. C APPELES (BLAS)   :
  27. C APPELES (CALCUL) :
  28. C APPELE PAR       :
  29. C***********************************************************************
  30. C SYNTAXE GIBIANE    :
  31. C    RIG2  = EXTR  RIG1 LMOT1 LMOT2 ;
  32. C
  33. C ENTREES            :
  34. C ENTREES/SORTIES    :
  35. C SORTIES            :
  36. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  37. C***********************************************************************
  38. C VERSION    : v1, 01/06/2011, version initiale
  39. C HISTORIQUE : v1, 01/06/2011, création
  40. C HISTORIQUE :
  41. C HISTORIQUE :
  42. C***********************************************************************
  43.  
  44. -INC PPARAM
  45. -INC CCOPTIO
  46. -INC SMRIGID
  47. -INC SMELEME
  48. -INC SMLMOTS
  49. *
  50. * Segment de correspondance : noms d'inconnues
  51. *
  52.       segment trav1
  53.       integer cor1p(nligp2),cor1d(nligd2)
  54.       endsegment
  55.       segment trav2
  56.       integer cor1n(nbnn2),cor2n(nbnn)
  57.       endsegment
  58. *
  59. * Segment extensibles pour ranger les rigidités
  60. * créés
  61. *
  62.       segment imatno(0)
  63.       segment idesno(0)
  64.       segment igeono(0)
  65.       segment coerno(0)
  66.       segment irg8no(0)
  67.       segment irg7no(0)
  68.       segment irg6no(0)
  69.       segment irg5no(0)
  70.       segment irg2no(0)
  71. *
  72.       logical lsym,lrela
  73. *
  74. * Executable statements
  75. *
  76.       SEGACT ri1
  77.       SEGACT MLMOT1,MLMOT2
  78. *
  79.       segini imatno,idesno,igeono,coerno,irg7no,irg5no
  80.       segini irg2no,irg6no,irg8no
  81.       do 1000 irig=1,ri1.irigel(/2)
  82.          des1=ri1.irigel(3,irig)
  83.          segact des1
  84.          nligrp=des1.lisinc(/2)
  85.          nligrd=des1.lisdua(/2)
  86.          nligp2=nligrp
  87.          nligd2=nligrd
  88.          segini trav1
  89. * On regarde le descripteur et les inconnues concernées
  90.          nligp2=0
  91.          do 410 iligrp=1,nligrp
  92.             do 420 incp=1,mlmot1.mots(/2)
  93.                if (mlmot1.mots(incp).eq.des1.lisinc(iligrp)) then
  94.                   nligp2=nligp2+1
  95.                   cor1p(nligp2)=iligrp
  96.                   goto 410
  97.                endif
  98.  420        continue
  99.  410     continue
  100.          if (nligp2.eq.0) goto 999
  101.          nligd2=0
  102.          do 510 iligrd=1,nligrd
  103.             do 520 incd=1,mlmot2.mots(/2)
  104.                if (mlmot2.mots(incd).eq.des1.lisdua(iligrd)) then
  105.                   nligd2=nligd2+1
  106.                   cor1d(nligd2)=iligrd
  107.                   goto 510
  108.                endif
  109.  520        continue
  110.  510     continue
  111. *         WRITE(IOIMP,*) 'irig=',irig
  112. *         WRITE(IOIMP,*) 'cor1p='
  113. *         WRITE (IOIMP,2020) (cor1p(I),I=1,cor1p(/1))
  114. *         WRITE(IOIMP,*) 'cor1d='
  115. *         WRITE (IOIMP,2020) (cor1d(I),I=1,cor1d(/1))
  116. * 2020 FORMAT (20(2X,I4) )
  117.          if (nligd2.eq.0) goto 999
  118.          segadj,trav1
  119. * On fait un cas particulier pour les matrices de relations
  120. * Il faut vérifier qu'on a bien gardé la symétrie.
  121. * Il faut également que la relation ait gardé le multiplicateur
  122. * et porte au moins sur un ddl autre !
  123.          ipt1=ri1.irigel(1,irig)
  124.          segact ipt1
  125.          lrela=.false.
  126.          if (ipt1.itypel.eq.22) then
  127.             if (nligp2.eq.nligd2) then
  128.                do 700 ilig=1,nligp2
  129.                   if (cor1p(ilig).ne.cor1d(ilig)) goto 710
  130.  700           continue
  131.                lrela=.true.
  132.  710           continue
  133.             endif
  134.          endif
  135. * l'ordre de ces deux lignes est important
  136.          if (lrela.and.cor1p(1).ne.1) goto 998
  137.          if (lrela.and.nligp2.eq.1) goto 998
  138.          coerno(**)=ri1.coerig(irig)
  139.          irg2no(**) =ri1.irigel(2,irig)
  140.          irg5no(**) =ri1.irigel(5,irig)
  141.          irg6no(**) =ri1.irigel(6,irig)
  142.          irg8no(**) =ri1.irigel(8,irig)
  143. * Le cas particulier facile, où on garde toute la matrice
  144. * élémentaire : on recopie la matrice originale
  145.          if (nligp2.eq.nligrp.and.nligd2.eq.nligrd) then
  146.             igeono(**)=ri1.irigel(1,irig)
  147.             idesno(**) =des1
  148.             imatno(**) =ri1.irigel(4,irig)
  149.             irg7no(**) =ri1.irigel(7,irig)
  150.             goto 999
  151.          endif
  152. * Sinon, il faut réduire la matrice et son descripteur, voire même
  153. * la géométrie
  154. *    Réduction du descripteur
  155.          nligrp=nligp2
  156.          nligrd=nligd2
  157.          segini des2
  158.          do iligp2=1,nligp2
  159.             iligp1=cor1p(iligp2)
  160.             des2.lisinc(iligp2)=des1.lisinc(iligp1)
  161.             des2.noelep(iligp2)=des1.noelep(iligp1)
  162.          enddo
  163.          do iligd2=1,nligd2
  164.             iligd1=cor1d(iligd2)
  165.             des2.lisdua(iligd2)=des1.lisdua(iligd1)
  166.             des2.noeled(iligd2)=des1.noeled(iligd1)
  167.          enddo
  168.          idesno(**)=des2
  169. *    Réduction de la matrice
  170.          xmatr1=ri1.irigel(4,irig)
  171.          segact xmatr1
  172.          nelrig=xmatr1.re(/3)
  173.          segini xmatr2
  174.          do ielrig=1,nelrig
  175.             do iligp2=1,nligp2
  176.                iligp1=cor1p(iligp2)
  177.                do iligd2=1,nligd2
  178.                   iligd1=cor1d(iligd2)
  179.                   xmatr2.re(iligd2,iligp2,ielrig)=
  180.      &                       xmatr1.re(iligd1,iligp1,ielrig)
  181.                enddo
  182.             enddo
  183.          enddo
  184.          segdes xmatr1
  185.          imatno(**) =xmatr2
  186. *    Réduction éventuelle de la géométrie si des noeuds
  187. *    ne sont pas référencés
  188.          nbnn=ipt1.num(/1)
  189.          nbnn2=ipt1.num(/1)
  190.          segini trav2
  191.          nbnn2=0
  192.          do in=1,des2.noelep(/1)
  193.             ic1=cor2n(des2.noelep(in))
  194.             if (ic1.eq.0) then
  195.                nbnn2=nbnn2+1
  196.                cor1n(nbnn2)=des2.noelep(in)
  197.                cor2n(des2.noelep(in))=nbnn2
  198.             endif
  199.          enddo
  200.          do in=1,des2.noeled(/1)
  201.             ic1=cor2n(des2.noeled(in))
  202.             if (ic1.eq.0) then
  203.                nbnn2=nbnn2+1
  204.                cor1n(nbnn2)=des2.noeled(in)
  205.                cor2n(des2.noeled(in))=nbnn2
  206.             endif
  207.          enddo
  208. *       Si tous les noeuds ont été vus, on peut peut-être
  209. *       garder la géométrie
  210.          if (nbnn2.eq.nbnn) then
  211.             if (ipt1.itypel.eq.22.and.(.not.lrela)) then
  212.                segini,ipt2=ipt1
  213.                ipt2.itypel=28
  214.                segdes ipt2
  215.             else
  216.                ipt2=ipt1
  217.             endif
  218.          else
  219. *       Sinon, on ne garde que les noeuds vus dans le nouveau
  220. *       maillage et on modifie le descripteur
  221.             nbnn=nbnn2
  222.             nbelem=ipt1.num(/2)
  223.             nbsous=0
  224.             nbref=0
  225.             segini ipt2
  226.             ipt2.itypel=28
  227.             if (lrela) ipt2.itypel=22
  228.             do iel=1,nbelem
  229.                do ibnn2=1,nbnn2
  230.                   ibnn=cor1n(ibnn2)
  231.                   ipt2.num(ibnn2,iel)=ipt1.num(ibnn,iel)
  232.                enddo
  233.             enddo
  234.             do iligrp=1,des2.noelep(/1)
  235.                des2.noelep(iligrp)=cor2n(des2.noelep(iligrp))
  236.             enddo
  237.             do iligrd=1,des2.noeled(/1)
  238.                des2.noeled(iligrd)=cor2n(des2.noeled(iligrd))
  239.             enddo
  240.             segdes ipt2
  241.          endif
  242.          segsup trav2
  243.          igeono(**)=ipt2
  244.          segdes des2
  245. *     Garde-t-on la symétrie ou l'antisymétrie ? Par défaut non, mais si
  246. *     les ddl primaux et duaux que l'on garde sont les mêmes, oui !
  247.          lsym=.false.
  248.          if (ri1.irigel(7,irig).ne.2) then
  249.             if (nligp2.eq.nligd2) then
  250.                do 600 ilig=1,nligp2
  251.                   if (cor1p(ilig).ne.cor1d(ilig)) goto 610
  252.  600           continue
  253.                lsym=.true.
  254.  610           continue
  255.             endif
  256.          endif
  257.          if (lsym) then
  258.             irg7no(**)  =ri1.irigel(7,irig)
  259.             xmatr2.symre=ri1.irigel(7,irig)
  260.          else
  261.             irg7no(**)  =2
  262.             xmatr2.symre=2
  263.          endif
  264.          segdes xmatr2
  265.  998     continue
  266.          segdes ipt1
  267.  999     continue
  268.          segdes des1
  269.          segsup trav1
  270.  1000 CONTINUE
  271.       SEGDES MLMOT1,MLMOT2
  272. *
  273. *  il ne reste plus qu'a ranger dans ri2 les nouvelles raideurs engendrées
  274. *
  275.       NRIGEL=imatno(/1)
  276.       SEGINI RI2
  277.       RI2.MTYMAT=RI1.MTYMAT
  278.       RI2.IFORIG=RI1.IFORIG
  279.       SEGDES RI1
  280.       DO irig=1,NRIGEL
  281.          RI2.COERIG(irig)=coerno(irig)
  282.          RI2.IRIGEL(1,irig)=igeono(irig)
  283.          RI2.IRIGEL(2,irig)=irg2no(irig)
  284.          RI2.IRIGEL(3,irig)=idesno(irig)
  285.          RI2.IRIGEL(4,irig)=imatno(irig)
  286.          RI2.IRIGEL(5,irig)=irg5no(irig)
  287.          RI2.IRIGEL(6,irig)=irg6no(irig)
  288.          RI2.IRIGEL(7,irig)=irg7no(irig)
  289.          RI2.IRIGEL(8,irig)=irg8no(irig)
  290.       ENDDO
  291.       SEGDES RI2
  292.       segsup imatno,idesno,igeono,coerno,irg7no,irg5no
  293.       segsup irg2no,irg6no,irg8no
  294. *
  295. * Normal termination
  296. *
  297. *      IRET=0
  298.       RETURN
  299. *
  300. * End of subroutine EXINCR
  301. *
  302.       END
  303.  
  304.  
  305.  

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