Télécharger la plaquette Cast3M


fr en

Télécharger baryc5.eso

Retour à la liste

Numérotation des lignes : 

baryc5

  1. C BARYC5    SOURCE    GOUNAND   25/11/24    21:15:02     12406          
  2.       SUBROUTINE BARYC5(MELEMX,NKPVIR,TRAVK,KGRAV)
  3.       IMPLICIT REAL*8(A-H,O-Z)
  4.       IMPLICIT INTEGER(I-N)
  5. C***********************************************************************
  6. C NOM         : BARYC5
  7. C DESCRIPTION :
  8. C     CALCUL LE BARYCENTRE D'UN MAILLAGE ET LA METRIQUE MOYENNE ASSOCIEE
  9. C     Repris de baryce.eso.
  10. C     Par rapport à baryc2, ignore le noeud virtuel KPVIRT
  11. C     Gère le noeud virtuel KPVIRT
  12. C
  13. C     Par rapport à baryc3, on gère le MCOORD différemment, en vue de
  14. C     faire moins de SEGADJ à l'aide du segment TRAVK
  15. C
  16. C     Comme baryc4 mais avec un MELEMX
  17. C
  18. C     KGRAV est le numéro du nouveau noeud créé
  19. C
  20. C
  21. C LANGAGE     : ESOPE
  22. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SEMT/LTA)
  23. C               mél : gounand@semt2.smts.cea.fr
  24. C***********************************************************************
  25. C APPELES          :
  26. C APPELES (E/S)    :
  27. C APPELES (BLAS)   :
  28. C APPELES (CALCUL) :
  29. C APPELE PAR       :
  30. C***********************************************************************
  31. C SYNTAXE GIBIANE    :
  32. C ENTREES            :
  33. C ENTREES/SORTIES    :
  34. C SORTIES            :
  35. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  36. C***********************************************************************
  37. C VERSION    : v1, 26/09/2017, version initiale
  38. C HISTORIQUE : v1, 26/09/2017, création
  39. C HISTORIQUE :
  40. C HISTORIQUE :
  41. C***********************************************************************
  42. -INC PPARAM
  43. -INC CCOPTIO
  44. -INC CCREEL
  45. -INC TMATOP2
  46. -INC TMATOP1
  47. *-INC SMELEMX
  48. -INC SMCOORD
  49.       POINTEUR KCOORD.MCOORD
  50. *-INC SMETRIQ
  51.       POINTEUR KCMETR.METRIQ
  52. *-INC STRAVJ
  53.       POINTEUR TRAVK.TRAVJ
  54.       LOGICAL lchang
  55.       PARAMETER(NGRAV=3)
  56.       DIMENSION XGRAV(NGRAV)
  57.       PARAMETER(NMET=6)
  58.       DIMENSION XMET(NMET)
  59.       DIMENSION A(3,3)
  60.       DIMENSION IDXSYM(3,3,3)
  61.       DIMENSION INVSYM(2,6,3)
  62.       DIMENSION XSUM(3),TLNSUM(3,3)
  63.       DIMENSION DXGRAV(3)
  64.       DIMENSION DLMGRA(3,3)
  65.       LOGICAL LDBG
  66. *
  67.       DATA ((A(I,J),I=1,3),J=1,3) /9*0.D0/
  68.       DATA (((IDXSYM(I,J,K),I=1,1),J=1,1),K=1,1) /1/
  69.       DATA (((IDXSYM(I,J,K),I=1,2),J=1,2),K=2,2) /1,2,2,3/
  70.       DATA (((IDXSYM(I,J,K),I=1,3),J=1,3),K=3,3) /1,2,4,2,3,5,4,5,6/
  71. *
  72.       DATA (((INVSYM(I,J,K),I=1,2),J=1,1),K=1,1) /1,1/
  73.       DATA (((INVSYM(I,J,K),I=1,2),J=1,3),K=2,2) /1,1,2,1,2,2/
  74.       DATA (((INVSYM(I,J,K),I=1,2),J=1,6),K=3,3)
  75.      $     /1,1,2,1,2,2,3,1,3,2,3,3/
  76. *
  77. *
  78. * Executable statements
  79. *
  80. *     SEGACT,MCOORD
  81. *      LDBG=.TRUE.
  82.       LDBG=impr.ge.3
  83. *      LDBG=.FALSE.
  84.       IDIMP1=IDIM+1
  85.       KCOORD=TRAVK.COORD
  86.       KCMETR=TRAVK.CMETR
  87. *     mis dans topv2      SEGACT KCOORD*MOD
  88.       DO I=1,NGRAV
  89.          XGRAV(I)=0.D0
  90.       ENDDO
  91.       DO I=1,NMET
  92.          XMET(I)=0.D0
  93.       ENDDO
  94.       NPOIN=0
  95. *      SEGACT,MELEMX
  96. *      DO 3 J=1,IPT1.NUM(/2)
  97. *     DO 31 I=1,IPT1.NUM(/1)
  98.       IF (ldbg) write(ioimp,*) 'BARYC5: IMET,NLCOU,NNCOU,NKPVIR='
  99.      $     ,IMET,NLCOU,NNCOU,NKPVIR
  100.       DO 3 ILCOU=1,NLCOU
  101.          DO 31 INCOU=1,NNCOU
  102.             INO=MELEMX.NUMX(INCOU,ILCOU)
  103.             IF (NKPVIR.GT.0) THEN
  104. *               ldbg=.true.
  105.                IF (INO.LE.NKPVIR) GOTO 31
  106.             ENDIF
  107.             NPOIN=NPOIN+1
  108.             IREF=IDIMP1*(INO-1)
  109.             DO 5 L=1,IDIM
  110.                XGRAV(L)=XGRAV(L)+KCOORD.XCOOR(IREF+L)
  111.  5          CONTINUE
  112.             IF (KCMETR.NE.0) THEN
  113.                if (ldbg) then
  114.                   write(ioimp,*) 'ilcou,incou,ino=',ilcou,incou,ino
  115.                   write(ioimp,*) '  coo=',(KCOORD.XCOOR(IREF+L),L=1
  116.      $                 ,IDIM)
  117.                   write(ioimp,*) '  met=',(KCMETR.XIN(L,INO),L=1
  118.      $                 ,KCMETR.XIN(/1))
  119.                endif
  120.                DO 6 ININ=1,KCMETR.XIN(/1)
  121.                   XMET(ININ)=XMET(ININ)+KCMETR.XIN(ININ,INO)
  122.  6             CONTINUE
  123.             ENDIF
  124.  31      CONTINUE
  125.  3    CONTINUE
  126.       if (ldbg) write(ioimp,*) 'NPOIN=',NPOIN
  127. *      SEGDES,MELEMX
  128. C     ON MET LE CENTRE DE GRAVITE DANS LA TABLE DES POINTS
  129. C     ET LA METRIQUE ASSOCIEE LE CAS ECHANT
  130.       NPCOUN=TRAVK.NPCOU+1
  131.       CALL TOPADP(TRAVK,NPCOUN,1,lchang,'baryc5 : TRAVK')
  132.       if (ierr.ne.0) return
  133.       if (iveri.ge.2) then
  134.          call vetopi(travk,'baryc5 : Apres extension travk 1')
  135.          if (ierr.ne.0) return
  136.       endif
  137.       do i=1,idim
  138.          xgrav(i)=xgrav(i)/NPOIN
  139.       enddo
  140.       IF (KCMETR.NE.0) THEN
  141.          DO 12 ININ=1,KCMETR.XIN(/1)
  142.             KCMETR.XIN(ININ,NPCOUN)=XMET(ININ)/NPOIN
  143.  12      CONTINUE
  144.       ENDIF
  145.       if (ldbg) then
  146.          write(ioimp,*) 'Centre de gravite'
  147.          write(ioimp,*) '  coo=',(xgrav(L),L=1,IDIM)
  148.          write(ioimp,*) '  met=',(KCMETR.XIN(L,NPCOUN),L=1,KCMETR.XIN(
  149.      $        /1))
  150.       endif
  151. *
  152. * On essaie de mieux placer le point milieu
  153. *
  154.       if (kcmetr.ne.0.and.imomet.eq.1.and.imobar.eq.1) then
  155. *
  156.          tsum=0.d0
  157.          do j=1,idim
  158.             xsum(j)=0.d0
  159.             dxgrav(j)=0.d0
  160.             do i=1,idim
  161.                dlmgra(i,j)=0.d0
  162.                tlnsum(i,j)=0.d0
  163.             enddo
  164.          enddo
  165.          if (imet.eq.4) then
  166.             imx=idim
  167.          else
  168.             imx=1
  169.          endif
  170.  
  171.          DO 4 INL=1,NLCOU
  172.             DO 41 INN=1,NNCOU
  173.                INO=MELEMX.NUMX(INN,INL)
  174.                IF (NKPVIR.GT.0) THEN
  175.                   IF (INO.LE.NKPVIR) GOTO 41
  176.                ENDIF
  177.                IREF=IDIMP1*(INO-1)
  178. * +1 car c'est le log de la metrique inverse
  179.                DO J=1,imx
  180.                   DO I=1,imx
  181. *                     A(I,J)=KCMETR.XIN(IDXSYM(I,J,imx),INO)
  182.                      A(I,J)=(KCMETR.XIN(IDXSYM(I,J,imx),INO)
  183.      $                    +KCMETR.XIN(IDXSYM(I,J,imx),NPCOUN))/2
  184. *                     A(I,J)=KCMETR.XIN(IDXSYM(I,J,imx),NPCOUN)
  185.                   enddo
  186.                ENDDO
  187.                if (ldbg) then
  188.                   write(ioimp,*) 'Log Tenseur metrique inverse voulu'
  189.                   write(ioimp,*) '  a=',((a(L,K),L=1,imx),K=1,imx)
  190.                endif
  191. * Projection suivant la direction centre de gravite-noeud puis exponentielle
  192.                if (imet.eq.4) then
  193.                   aproj=0.d0
  194.                   do j=1,idim
  195.                      do i=1,idim
  196.                         aproj=aproj+(KCOORD.XCOOR(IREF+i)-XGRAV(i))*a(i
  197.      $                       ,j)*(KCOORD.XCOOR(IREF+j)-XGRAV(j))
  198.                      enddo
  199.                   enddo
  200.                   dx2=0.d0
  201.                   do i=1,idim
  202.                      dx2=dx2+(KCOORD.XCOOR(IREF+i)-XGRAV(i))**2
  203.                   enddo
  204.                   aproj=aproj/dx2
  205.                   xmii=exp(aproj)
  206. *    Imet = 3
  207.                ELSE
  208.                   xmii=exp(a(1,1))
  209.                ENDIF
  210.                if (ldbg) then
  211.                   write(ioimp,*) 'Tenseur metrique inverse voulu'
  212.                   write(ioimp,*) '  xmii=',xmii
  213.                endif
  214.                tsum=tsum+xmii
  215.                DO J=1,imx
  216.                   DO I=1,imx
  217.                      tlnsum(i,j)=tlnsum(i,j)+xmii*(KCMETR.XIN(IDXSYM(I,J
  218.      $                    ,imx),INO)-KCMETR.XIN(IDXSYM(I,J,imx),NPCOUN))
  219.                   ENDDO
  220.                ENDDO
  221.                do i=1,idim
  222.                   xsum(i)=xsum(i)+xmii*(KCOORD.XCOOR(IREF+i)
  223.      $                 -XGRAV(i))
  224.                enddo
  225.  41         CONTINUE
  226.  4       CONTINUE
  227.          DO J=1,imx
  228.             DO I=1,imx
  229.                tlnsum(i,j)=tlnsum(i,j)/NPOIN
  230.             ENDDO
  231.          ENDDO
  232.          do i=1,idim
  233.             xsum(i)=xsum(i)/NPOIN
  234.          enddo
  235.          tsum=tsum/NPOIN
  236.          if (ldbg) then
  237.             write(ioimp,*) 'Moyennes'
  238.             write(ioimp,*) '  xsum=',(xsum(L),L=1,IDIM)
  239.             write(ioimp,*) '  tsum=',tsum
  240.             write(ioimp,*) '  tlnsum=',((tlnsum(L,K),L=1,imx),K=1,imx)
  241.          endif
  242.          xrii=1.d0/tsum
  243.          DO J=1,imx
  244.             DO I=1,imx
  245.                dlmgra(i,j)=dlmgra(i,j)+xrii*tlnsum(i,j)
  246.             ENDDO
  247.          ENDDO
  248.          do i=1,idim
  249.             dxgrav(i)=dxgrav(i)+xrii*xsum(i)
  250.          enddo
  251.          if (ldbg) then
  252.             write(ioimp,*) 'Depl gravite ameliore ?'
  253.             write(ioimp,*) ' dxgrav=',(dxgrav(L),L=1,IDIM)
  254.             write(ioimp,*) 'DLog metrique gravite ameliore ?'
  255.             write(ioimp,*) '  dlmgra=',((dlmgra(L,K),L=1,imx),K=1,imx)
  256.          endif
  257.          do l=1,idim
  258.             xgrav(l)=xgrav(l)+dxgrav(l)
  259.          enddo
  260.          DO ININ=1,KCMETR.XIN(/1)
  261.             I=INVSYM(1,ININ,IDIM)
  262.             J=INVSYM(2,ININ,IDIM)
  263.             KCMETR.XIN(ININ,NPCOUN)=KCMETR.XIN(ININ,NPCOUN)+DLMGRA(I,J)
  264.          ENDDO
  265. *
  266.          if (ldbg) then
  267.             write(ioimp,*) 'Centre de gravite ameliore ?'
  268.             write(ioimp,*) '  coo=',(xgrav(L),L=1,IDIM)
  269.             write(ioimp,*) '  met=',(KCMETR.XIN(L,NPCOUN),L=1
  270.      $           ,KCMETR.XIN(/1))
  271.          endif
  272.       endif
  273. *      write(ioimp,*) 'npcoun,npcou,npmax',npcoun,travk.npcou,travk.npmax
  274. *      NBPTS=XCOOR(/1)/IDIMP1+1
  275. *      SEGADJ,MCOORD
  276. *      IREF=(NBPTS-1)*IDIMP1
  277.       IREF=(NPCOUN-1)*IDIMP1
  278.       DO 11 I=1,IDIMP1
  279.          KCOORD.XCOOR(IREF+I)=XGRAV(I)
  280.  11   CONTINUE
  281. *     KGRAV=XCOOR(/1)/IDIMP1
  282. *      if (ldbg) then
  283. *         write(ioimp,*) 'Stop !'
  284. *         ierr=1
  285. *         return
  286. *      endif
  287.       KGRAV=NPCOUN
  288.       if (iveri.ge.2.and.lchang) then
  289.          call vetopi(travk,'baryc5 : Apres extension travk 2')
  290.          if (ierr.ne.0) return
  291.       endif
  292.  
  293.       RETURN
  294.       END
  295.  
  296.  

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