Télécharger la plaquette Cast3M


fr en

Télécharger mesu.eso

Retour à la liste

Numérotation des lignes : 

mesu

  1. C MESU      SOURCE    PASCAL    21/03/19    21:15:10     10924          
  2.  
  3. C   fournit la mesure d'un maillage
  4.  
  5.       SUBROUTINE mesu
  6.  
  7.       IMPLICIT INTEGER(I-N)
  8.       IMPLICIT REAL*8 (A-H,O-Z)
  9.  
  10.  
  11. -INC PPARAM
  12. -INC CCOPTIO
  13. -INC CCREEL
  14. -INC CCGEOME
  15. -INC SMCOORD
  16. -INC SMELEME
  17.  
  18.       PARAMETER (NCLE=4)
  19.       SEGMENT /FER/ (NFI(ITT),MAI(IPP),ITOUR)
  20.       CHARACTER*4 mcle(NCLE)
  21.       DATA mcle / 'LONG','SURF','VOLU','DENS' /
  22.  
  23. *dbg      write(ioimp,*) 'coucou mesu'
  24.       segact mcoord
  25.       icle=0
  26.       CALL LIRMOT(mcle,ncle,icle,0)
  27.       IF (IERR.NE.0) RETURN
  28.       CALL LIROBJ('MAILLAGE',meleme,1,iOK)
  29.       CALL ACTOBJ('MAILLAGE',meleme,1)
  30.       IF (IERR.NE.0) RETURN
  31.  
  32.       IF (ICLE.EQ.4) THEN
  33.         CALL MESUDE(meleme)
  34.         IF (IERR.NE.0) RETURN
  35.         RETURN
  36.       ENDIF
  37.  
  38.       idimp1=IDIM+1
  39.       SEGACT,meleme
  40.       ipt1=meleme
  41.       itp=0
  42.       xmesu=xzero
  43.       nsous=lisous(/1)
  44.       DO isou=1,MAX(nsous,1)
  45.          IF (nsous.NE.0) THEN
  46.             ipt1=lisous(isou)
  47.             SEGACT,ipt1
  48.          ENDIF
  49.          jtyp=3
  50.          IF (ksurf(ipt1.itypel).EQ.ipt1.itypel) jtyp=2
  51.          IF (kdegre(ipt1.itypel).EQ.ipt1.itypel) jtyp=1
  52.          IF (jtyp.EQ.1.AND.nsous.GT.1) CALL ERREUR(25)
  53.          IF (itp.EQ.0) itp=jtyp
  54.          IF (itp.NE.jtyp) CALL ERREUR(25)
  55.          IF (icle.EQ.0) icle=itp
  56.          IF (icle.EQ.1.AND.itp.NE.1) CALL ERREUR(21)
  57.          IF (icle.EQ.2.AND.itp.EQ.3) CALL ERREUR(21)
  58.          IF (icle.EQ.3.AND.itp.EQ.1) CALL ERREUR(21)
  59.          IF (IERR.NE.0) RETURN
  60.          GOTO (100,200,300),icle
  61. C  calcul de la longueur d'une ligne
  62.  100     IF (IIMPI.NE.0) WRITE(6,*) ' Mesure d''une ligne'
  63.          nbnn=ipt1.num(/1)
  64.          IF (IDIM.EQ.1) THEN
  65.             DO iel=1,ipt1.num(/2)
  66.                ip1=idimp1*(ipt1.num(1,iel)-1)
  67.                ip2=idimp1*(ipt1.num(nbnn,iel)-1)
  68.                xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
  69.                xmesu=xmesu+ABS(xp12)
  70.             ENDDO
  71.          ELSE IF (IDIM.EQ.2) THEN
  72.             DO iel=1,ipt1.num(/2)
  73.                ip1=idimp1*(ipt1.num(1,iel)-1)
  74.                ip2=idimp1*(ipt1.num(nbnn,iel)-1)
  75.                xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
  76.                yp12=XCOOR(ip2+2)-XCOOR(ip1+2)
  77.                xmesu=xmesu+SQRT(xp12*xp12+yp12*yp12)
  78.             ENDDO
  79.          ELSE IF (IDIM.EQ.3) THEN
  80.             DO iel=1,ipt1.num(/2)
  81.                ip1=idimp1*(ipt1.num(1,iel)-1)
  82.                ip2=idimp1*(ipt1.num(nbnn,iel)-1)
  83.                xp12=XCOOR(ip2+1)-XCOOR(ip1+1)
  84.                yp12=XCOOR(ip2+2)-XCOOR(ip1+2)
  85.                zp12=XCOOR(ip2+3)-XCOOR(ip1+3)
  86.                xmesu=xmesu+SQRT(xp12*xp12+yp12*yp12+zp12*zp12)
  87.             ENDDO
  88.          ENDIF
  89.          GOTO 500
  90. C  calcul de la surface d'une aire
  91.  200     GOTO (220,240),jtyp
  92. C     cas d'une ligne
  93.  220     IF (IIMPI.NE.0) WRITE(6,*) ' Mesure de surface dans un contour'
  94.          CALL AVTRSF(ipt1,fer,ipt5)
  95.          SEGSUP,ipt5
  96.          ib=0
  97.          xb=xzero
  98.          yb=xzero
  99.          zb=xzero
  100.          DO it=1,itour
  101.             imdeb=mai(it)+1
  102.             imfin=mai(it+1)
  103.             DO im=imdeb,imfin
  104.                ib=ib+1
  105.                ip=idimp1*(nfi(im)-1)
  106.                xb=xb+XCOOR(ip+1)
  107.                yb=yb+XCOOR(ip+2)
  108.                IF (IDIM.EQ.3) zb=zb+XCOOR(ip+3)
  109.             ENDDO
  110.          ENDDO
  111.          xb=xb/ib
  112.          yb=yb/ib
  113.          zb=zb/ib
  114.          IF (IIMPI.NE.0) WRITE(6,*) ' Barycentre ',xb,yb,zb
  115.          xs=xzero
  116.          ys=xzero
  117.          zs=xzero
  118.          zv1=xzero
  119.          zv2=xzero
  120.          DO it=1,itour
  121.             imdeb=mai(it)+1
  122.             imfin=mai(it+1)
  123.             DO im=imdeb,imfin
  124.                im1=im+1
  125.                IF (im1.GT.imfin) im1=imdeb
  126.                ip1=idimp1*(nfi(im)-1)
  127.                ip2=idimp1*(nfi(im1)-1)
  128.                xv1=XCOOR(ip1+1)-xb
  129.                yv1=XCOOR(ip1+2)-yb
  130.                xv2=XCOOR(ip2+1)-xb
  131.                yv2=XCOOR(ip2+2)-yb
  132.                IF (IDIM.GE.3) THEN
  133.                   zv1=XCOOR(ip1+3)-zb
  134.                   zv2=XCOOR(ip2+3)-zb
  135.                ENDIF
  136.                xs=xs+(yv1*zv2-zv1*yv2)
  137.                ys=ys+(zv1*xv2-xv1*zv2)
  138.                zs=zs+(xv1*yv2-yv1*xv2)
  139.             ENDDO
  140.          ENDDO
  141.          xmesu=xmesu+0.5D0*SQRT(xs*xs+ys*ys+zs*zs)
  142.          SEGSUP,fer
  143.          GOTO 500
  144. C     cas d'une surface
  145.  240     IF (IIMPI.NE.0) WRITE(6,*) ' Mesure d''une surface'
  146.          iad=LTEL(2,ipt1.itypel)-1
  147.          jad=LDEL(2,iad+1)-1
  148.          ityp=LDEL(1,iad+1)
  149.          npfac=KDFAC(1,ityp)
  150.          idep=KDFAC(2,ityp)
  151.          ifep=idep+3*(KDFAC(3,ityp)-1)
  152.          xs=xzero
  153.          ys=xzero
  154.          zs=xzero
  155.          Cc=xzero
  156.          DO iel=1,ipt1.num(/2)
  157.             DO i=idep,ifep,3
  158.                iafa=LFAC(jad+KFAC(i))
  159.                ibfa=LFAC(jad+KFAC(i+1))
  160.                icfa=LFAC(jad+KFAC(i+2))
  161.                ia=ipt1.num(iafa,iel)
  162.                ib=ipt1.num(ibfa,iel)
  163.                ic=ipt1.num(icfa,iel)
  164.                IF (IIMPI.EQ.5) WRITE(6,*) ' Triangle ',ia,ib,ic
  165.                ia=idimp1*(ia-1)
  166.                ib=idimp1*(ib-1)
  167.                ic=idimp1*(ic-1)
  168.                xb=XCOOR(ib+1)
  169.                yb=XCOOR(ib+2)
  170.                xv1=XCOOR(ia+1)-xb
  171.                yv1=XCOOR(ia+2)-yb
  172.                xv2=XCOOR(ic+1)-xb
  173.                yv2=XCOOR(ic+2)-yb
  174.                IF (IDIM.GE.3) THEN
  175.                   zb=XCOOR(ib+3)
  176.                   zv1=XCOOR(ia+3)-zb
  177.                   zv2=XCOOR(ic+3)-zb
  178.                   xs=yv1*zv2-zv1*yv2
  179.                   ys=zv1*xv2-xv1*zv2
  180.                ENDIF
  181.                zs=xv1*yv2-yv1*xv2
  182.                Cc=Cc+SQRT(xs*xs+ys*ys+zs*zs)
  183.             ENDDO
  184.          ENDDO
  185.          xmesu=xmesu+0.5D0*Cc
  186.          GOTO 500
  187. C  calcul d'un volume
  188.  300     GOTO (320,320,340),jtyp
  189. C  cas du volume interieur a une surface
  190. C  du fait du probleme d'orientation des elements on ne peut pas
  191. C  traiter le probleme sous-maillage par sous maillage
  192.  320     CALL MESUVO(meleme,xmesu)
  193.          GOTO 11
  194. C  cas d'un volume. On va utiliser le decoupage en triangle de la
  195. C  surface des elements mais en le dedoublant de facon a etre sur
  196. C  de correspondre d'un element au voisin
  197. C     SG 20160712 : pour les faces TRI7, QUA9 (ityp=7,8), il n'y a pas
  198. C     besoin de dedoubler car le decoupage de la face est symetrique
  199. C     (barycentrique) 
  200.  340     nbfac=LTEL(1,ipt1.itypel)
  201.          IF (nbfac.EQ.0) GOTO 500
  202.          iad=LTEL(2,ipt1.itypel)-1
  203.          kdeg=KDEGRE(ipt1.itypel)-1
  204.          Cc=xzero
  205.          DO iel=1,ipt1.num(/2)
  206. C     centre de gravite du volume
  207.             xb=xzero
  208.             yb=xzero
  209.             zb=xzero
  210.             nbnn=ipt1.num(/1)
  211.             DO i=1,nbnn
  212.                ipt=idimp1*(ipt1.num(i,iel)-1)
  213.                xb=xb+XCOOR(ipt+1)
  214.                yb=yb+XCOOR(ipt+2)
  215.                zb=zb+XCOOR(ipt+3)
  216.             ENDDO
  217.             xb=xb/nbnn
  218.             yb=yb/nbnn
  219.             zb=zb/nbnn
  220.             DO ifac=1,nbfac
  221.                v=xzero
  222.                ityp=LDEL(1,iad+ifac)
  223.                npfac=KDFAC(1,ityp)
  224.                jad=LDEL(2,iad+ifac)-1
  225.                idep=KDFAC(2,ityp)
  226.                ifep=idep+3*(KDFAC(3,ityp)-1)
  227.                DO i=idep,ifep,3
  228.                   iafa=LFAC(jad+KFAC(i))
  229.                   ibfa=LFAC(jad+KFAC(i+1))
  230.                   icfa=LFAC(jad+KFAC(i+2))
  231.                   ia=ipt1.num(iafa,iel)
  232.                   ib=ipt1.num(ibfa,iel)
  233.                   ic=ipt1.num(icfa,iel)
  234.                   v=v+VOLTRI(xb,yb,zb,ia,ib,ic)
  235.                ENDDO
  236. * SG
  237.                IF (ITYP.NE.7.AND.ITYP.NE.8) THEN
  238.                   DO i=idep,ifep,3
  239.                      ifa=KFAC(i)+kdeg
  240.                      IF (ifa.GT.npfac) ifa=ifa-npfac
  241.                      ifb=KFAC(i+1)+kdeg
  242.                      IF (ifb.GT.npfac) ifb=ifb-npfac
  243.                      ifc=KFAC(i+2)+kdeg
  244.                      IF (ifc.GT.npfac) ifc=ifc-npfac
  245.                      iafa=LFAC(jad+ifa)
  246.                      ibfa=LFAC(jad+ifb)
  247.                      icfa=LFAC(jad+ifc)
  248.                      ia=ipt1.num(iafa,iel)
  249.                      ib=ipt1.num(ibfa,iel)
  250.                      ic=ipt1.num(icfa,iel)
  251.                      v=v+VOLTRI(xb,yb,zb,ia,ib,ic)
  252.                   ENDDO
  253.                ELSE
  254.                   v=v*2
  255.                ENDIF
  256.                IF (IIMPI.EQ.1) WRITE(6,*) 'v-2 ',v
  257.                Cc=Cc+ABS(v)
  258.             ENDDO
  259.          ENDDO
  260.          xmesu=xmesu+Cc/12.D0
  261.  500     CONTINUE
  262.       ENDDO
  263.  
  264.  11   CONTINUE
  265.       CALL ECRREE(xmesu)
  266.  
  267.       END
  268.  
  269.  
  270.  
  271.  
  272.  

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