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    CB215821  25/04/23    21:15:20     12247          
  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 SMCOORD
  48. -INC SMELEME
  49. -INC SMLMOTS
  50. *
  51. * Segment de correspondance : noms d'inconnues
  52. *
  53.       segment trav1
  54.       integer cor1p(nligp2),cor1d(nligd2)
  55.       endsegment
  56.       segment trav2
  57.       integer cor1n(nbnn2),cor2n(nbnn)
  58.       endsegment
  59. *
  60. * Segment extensibles pour ranger les rigidités
  61. * créés
  62. *
  63.       segment imatno(0)
  64.       segment idesno(0)
  65.       segment igeono(0)
  66.       segment coerno(0)
  67.       segment irg8no(0)
  68.       segment irg7no(0)
  69.       segment irg6no(0)
  70.       segment irg5no(0)
  71.       segment irg2no(0)
  72. *
  73.       logical lsym,lrela
  74. *
  75. * Executable statements
  76. *
  77.       SEGACT ri1
  78.       SEGACT MLMOT1,MLMOT2
  79. *
  80.       segini imatno,idesno,igeono,coerno,irg7no,irg5no
  81.       segini irg2no,irg6no,irg8no
  82.       do 1000 irig=1,ri1.irigel(/2)
  83.          des1=ri1.irigel(3,irig)
  84.          segact des1
  85.          nligrp=des1.lisinc(/2)
  86.          nligrd=des1.lisdua(/2)
  87.          nligp2=nligrp
  88.          nligd2=nligrd
  89.          segini trav1
  90. * On regarde le descripteur et les inconnues concernées
  91.          nligp2=0
  92.          do 410 iligrp=1,nligrp
  93.             do 420 incp=1,mlmot1.mots(/2)
  94.                if (mlmot1.mots(incp).eq.des1.lisinc(iligrp)) then
  95.                   nligp2=nligp2+1
  96.                   cor1p(nligp2)=iligrp
  97.                   goto 410
  98.                endif
  99.  420        continue
  100.  410     continue
  101.          if (nligp2.eq.0) goto 999
  102.          nligd2=0
  103.          do 510 iligrd=1,nligrd
  104.             do 520 incd=1,mlmot2.mots(/2)
  105.                if (mlmot2.mots(incd).eq.des1.lisdua(iligrd)) then
  106.                   nligd2=nligd2+1
  107.                   cor1d(nligd2)=iligrd
  108.                   goto 510
  109.                endif
  110.  520        continue
  111.  510     continue
  112. *         WRITE(IOIMP,*) 'irig=',irig
  113. *         WRITE(IOIMP,*) 'cor1p='
  114. *         WRITE (IOIMP,2020) (cor1p(I),I=1,cor1p(/1))
  115. *         WRITE(IOIMP,*) 'cor1d='
  116. *         WRITE (IOIMP,2020) (cor1d(I),I=1,cor1d(/1))
  117. * 2020 FORMAT (20(2X,I4) )
  118.          if (nligd2.eq.0) goto 999
  119.          segadj,trav1
  120. * On fait un cas particulier pour les matrices de relations
  121. * Il faut vérifier qu'on a bien gardé la symétrie.
  122. * Il faut également que la relation ait gardé le multiplicateur
  123. * et porte au moins sur un ddl autre !
  124.          ipt1=ri1.irigel(1,irig)
  125.          segact ipt1
  126.          lrela=.false.
  127.          if (ipt1.itypel.eq.22) then
  128.             if (nligp2.eq.nligd2) then
  129.                do 700 ilig=1,nligp2
  130.                   if (cor1p(ilig).ne.cor1d(ilig)) goto 710
  131.  700           continue
  132.                lrela=.true.
  133.  710           continue
  134.             endif
  135.          endif
  136. * l'ordre de ces deux lignes est important
  137.          if (lrela.and.cor1p(1).ne.1) goto 998
  138.          if (lrela.and.nligp2.eq.1) goto 998
  139.          coerno(**)=ri1.coerig(irig)
  140.          irg2no(**) =ri1.irigel(2,irig)
  141.          irg5no(**) =ri1.irigel(5,irig)
  142.          irg6no(**) =ri1.irigel(6,irig)
  143.          irg8no(**) =ri1.irigel(8,irig)
  144. * Le cas particulier facile, où on garde toute la matrice
  145. * élémentaire : on recopie la matrice originale
  146.          if (nligp2.eq.nligrp.and.nligd2.eq.nligrd) then
  147.             igeono(**)=ri1.irigel(1,irig)
  148.             idesno(**) =des1
  149.             imatno(**) =ri1.irigel(4,irig)
  150.             irg7no(**) =ri1.irigel(7,irig)
  151.             goto 999
  152.          endif
  153. * Sinon, il faut réduire la matrice et son descripteur, voire même
  154. * la géométrie
  155. *    Réduction du descripteur
  156.          nligrp=nligp2
  157.          nligrd=nligd2
  158.          segini des2
  159.          do iligp2=1,nligp2
  160.             iligp1=cor1p(iligp2)
  161.             des2.lisinc(iligp2)=des1.lisinc(iligp1)
  162.             des2.noelep(iligp2)=des1.noelep(iligp1)
  163.          enddo
  164.          do iligd2=1,nligd2
  165.             iligd1=cor1d(iligd2)
  166.             des2.lisdua(iligd2)=des1.lisdua(iligd1)
  167.             des2.noeled(iligd2)=des1.noeled(iligd1)
  168.          enddo
  169.          idesno(**)=des2
  170. *    Réduction de la matrice
  171.          xmatr1=ri1.irigel(4,irig)
  172.          segact xmatr1
  173.          nelrig=xmatr1.re(/3)
  174.          segini xmatr2
  175.          do ielrig=1,nelrig
  176.             do iligp2=1,nligp2
  177.                iligp1=cor1p(iligp2)
  178.                do iligd2=1,nligd2
  179.                   iligd1=cor1d(iligd2)
  180.                   xmatr2.re(iligd2,iligp2,ielrig)=
  181.      &                       xmatr1.re(iligd1,iligp1,ielrig)
  182.                enddo
  183.             enddo
  184.          enddo
  185.          segdes xmatr1
  186.          imatno(**) =xmatr2
  187. *    Réduction éventuelle de la géométrie si des noeuds
  188. *    ne sont pas référencés
  189.          nbnn=ipt1.num(/1)
  190.          nbnn2=ipt1.num(/1)
  191.          segini trav2
  192.          nbnn2=0
  193.          do in=1,des2.noelep(/1)
  194.             ic1=cor2n(des2.noelep(in))
  195.             if (ic1.eq.0) then
  196.                nbnn2=nbnn2+1
  197.                cor1n(nbnn2)=des2.noelep(in)
  198.                cor2n(des2.noelep(in))=nbnn2
  199.             endif
  200.          enddo
  201.          do in=1,des2.noeled(/1)
  202.             ic1=cor2n(des2.noeled(in))
  203.             if (ic1.eq.0) then
  204.                nbnn2=nbnn2+1
  205.                cor1n(nbnn2)=des2.noeled(in)
  206.                cor2n(des2.noeled(in))=nbnn2
  207.             endif
  208.          enddo
  209. *       Si tous les noeuds ont été vus, on peut peut-être
  210. *       garder la géométrie
  211.          if (nbnn2.eq.nbnn) then
  212.             if (ipt1.itypel.eq.22.and.(.not.lrela)) then
  213.                segini,ipt2=ipt1
  214.                ipt2.itypel=28
  215.                segdes ipt2
  216.             else
  217.                ipt2=ipt1
  218.             endif
  219.          else
  220. *       Sinon, on ne garde que les noeuds vus dans le nouveau
  221. *       maillage et on modifie le descripteur
  222.             nbnn=nbnn2
  223.             nbelem=ipt1.num(/2)
  224.             nbsous=0
  225.             nbref=0
  226.             segini ipt2
  227.             ipt2.itypel=28
  228.             if (lrela) ipt2.itypel=22
  229.             do iel=1,nbelem
  230.                do ibnn2=1,nbnn2
  231.                   ibnn=cor1n(ibnn2)
  232.                   ipt2.num(ibnn2,iel)=ipt1.num(ibnn,iel)
  233.                enddo
  234.             enddo
  235.             do iligrp=1,des2.noelep(/1)
  236.                des2.noelep(iligrp)=cor2n(des2.noelep(iligrp))
  237.             enddo
  238.             do iligrd=1,des2.noeled(/1)
  239.                des2.noeled(iligrd)=cor2n(des2.noeled(iligrd))
  240.             enddo
  241.             segdes ipt2
  242.          endif
  243.          segsup trav2
  244.          igeono(**)=ipt2
  245.          segdes des2
  246. *     Garde-t-on la symétrie ou l'antisymétrie ? Par défaut non, mais si
  247. *     les ddl primaux et duaux que l'on garde sont les mêmes, oui !
  248.          lsym=.false.
  249.          if (ri1.irigel(7,irig).ne.2) then
  250.             if (nligp2.eq.nligd2) then
  251.                do 600 ilig=1,nligp2
  252.                   if (cor1p(ilig).ne.cor1d(ilig)) goto 610
  253.  600           continue
  254.                lsym=.true.
  255.  610           continue
  256.             endif
  257.          endif
  258.          if (lsym) then
  259.             irg7no(**)  =ri1.irigel(7,irig)
  260.             xmatr2.symre=ri1.irigel(7,irig)
  261.          else
  262.             irg7no(**)  =2
  263.             xmatr2.symre=2
  264.          endif
  265.          segdes xmatr2
  266.  998     continue
  267.          segdes ipt1
  268.  999     continue
  269.          segdes des1
  270.          segsup trav1
  271.  1000 CONTINUE
  272.       SEGDES MLMOT1,MLMOT2
  273. *
  274. *  il ne reste plus qu'a ranger dans ri2 les nouvelles raideurs engendrées
  275. *
  276.       NRIGEL=imatno(/1)
  277.       SEGINI RI2
  278.       RI2.MTYMAT=RI1.MTYMAT
  279.       RI2.IFORIG=RI1.IFORIG
  280.       SEGDES RI1
  281.       DO irig=1,NRIGEL
  282.          RI2.COERIG(irig)=coerno(irig)
  283.          RI2.IRIGEL(1,irig)=igeono(irig)
  284.          RI2.IRIGEL(2,irig)=irg2no(irig)
  285.          RI2.IRIGEL(3,irig)=idesno(irig)
  286.          RI2.IRIGEL(4,irig)=imatno(irig)
  287.          RI2.IRIGEL(5,irig)=irg5no(irig)
  288.          RI2.IRIGEL(6,irig)=irg6no(irig)
  289.          RI2.IRIGEL(7,irig)=irg7no(irig)
  290.          RI2.IRIGEL(8,irig)=irg8no(irig)
  291.       ENDDO
  292.       SEGDES RI2
  293.       segsup imatno,idesno,igeono,coerno,irg7no,irg5no
  294.       segsup irg2no,irg6no,irg8no
  295. *
  296. * Normal termination
  297. *
  298. *      IRET=0
  299.       RETURN
  300. *
  301. * End of subroutine EXINCR
  302. *
  303.       END
  304.  
  305.  
  306.  
  307.  

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