Télécharger la plaquette Cast3M


fr en

Télécharger xpost1.eso

Retour à la liste

Numérotation des lignes : 

xpost1

  1. C XPOST1    SOURCE    CB215821  24/04/12    21:17:30     11897          
  2. c
  3.       SUBROUTINE XPOST1(IPCHP1,IPMOD1,IPCHP2)
  4. c
  5. c     Construit un Chpoint avec des ddl "physiques" (UX UY) en
  6. c     RECOmbinant les ddl Xfem (UX AX B1X C1X ...)
  7. c     -> Utile pour le tracé de déformées
  8. c
  9. C     SYNTAXE :
  10. c       chpo2 = XFEM 'RECO' chpo1 MODEL_XFEM
  11. C
  12. c
  13.       IMPLICIT INTEGER(I-N)
  14.       IMPLICIT REAL*8 (A-H,O-Z)
  15. C
  16.       PARAMETER   (XTOL1=1.d-7)
  17.       PARAMETER   (NBRMAX=5)
  18.  
  19. C
  20. C  SEGMENTS INCLUDE
  21. -INC CCREEL
  22.  
  23. -INC PPARAM
  24. -INC CCOPTIO
  25. -INC SMCOORD
  26. -INC SMELEME
  27. -INC SMCHPOI
  28. -INC SMCHAML
  29. -INC SMMODEL
  30. -INC SMLREEL
  31. c
  32.       POINTEUR    MCHEX1.MCHELM,MPHI1.MPOVAL,MPSI1.MPOVAL
  33. C
  34.  
  35. c       SEGMENT MRACC
  36. c          INTEGER        TLENR(NBENR1,NBPT)
  37. c       ENDSEGMENT
  38.       SEGMENT MRACC
  39.          INTEGER        TLENR(NBENR1,NBPT)
  40.          INTEGER        IPOCHA(NBPT,NCOMP1)
  41.       ENDSEGMENT
  42. c       SEGMENT MRACC
  43. c          INTEGER        NOD2PT(NBPT)
  44. c          INTEGER        TLENR(NBENR1,NBPT1)
  45. c          INTEGER        IPOCHA(NCOMP1,NBPT1)
  46. c       ENDSEGMENT
  47. C
  48. C
  49. C_____________________________________________________________
  50. C  INITIALISATION DES INCONNUES obligatoires et facultatives
  51.  
  52.       CHARACTER*(LOCOMP) MODDL,MODDL2
  53.  
  54.       PARAMETER (NOBL=4,NFAC=27)
  55.       CHARACTER*(LOCOMP) DDLOBL(NOBL),DDLFAC(NFAC)
  56.       DATA DDLOBL/'UX  ','UY  ','UZ  ','LX  '/
  57.       DATA DDLFAC/'AX  ','AY  ','AZ  ',
  58.      >            'B1X ','B1Y ','B1Z ',
  59.      >            'C1X ','C1Y ','C1Z ',
  60.      >            'D1X ','D1Y ','D1Z ',
  61.      >            'E1X ','E1Y ','E1Z ',
  62.      >            'B2X ','B2Y ','B2Z ',
  63.      >            'C2X ','C2Y ','C2Z ',
  64.      >            'D2X ','D2Y ','D2Z ',
  65.      >            'E2X ','E2Y ','E2Z '/
  66.  
  67.       INTEGER TENR(NFAC),TNI(NFAC),TF2O(NFAC)
  68. c     ddlfac correspond a quel enrichissement?
  69.       DATA TENR/1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3/
  70. c     ddlfac correspond a quel fonction d'enrichissement?
  71.       DATA  TNI/1,1,1,1,1,1,2,2,2,3,3,3,4,4,4,1,1,1,2,2,2,3,3,3,4,4,4/
  72. c     ddlfac correspond a quel ddlobl?
  73.       DATA TF2O/1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3/
  74. c     tables pour mise en concordance des ddl
  75.       INTEGER TOBL(NOBL),TFAC(NOBL,4,NBRMAX)
  76.       INTEGER           TIFAC(NOBL,4,NBRMAX)
  77. c     fonctions de forme
  78.       REAL*8  SHPWRK(4)
  79.  
  80.       nbpt = nbpts
  81.  
  82. C_____________________________________________________________
  83. C  ACTIVATION ...
  84.  
  85. C     ...DU MODELE
  86.       MMODEL=IPMOD1
  87.       NSOUS = KMODEL(/1)
  88.  
  89. C     ...DU CHPOINT avec ddl xfem
  90.       MCHPOI=IPCHP1
  91.       NSOUP0 = IPCHP(/1)
  92.  
  93. C_____________________________________________________________
  94. C  TRAVAIL PRELIMINAIRE SUR L ENRICHISSEMENT'
  95. C  on les transforme en 2 chpoints psi et phi
  96. c  avec le type d enrichissement pour composante
  97.  
  98.       NBENR1= NBRMAX
  99.       NCOMP1=NFAC
  100.       SEGINI,MRACC
  101.       N=NBPT
  102.       NC=NBRMAX
  103.       segini,MPHI1,MPSI1
  104.       NBPT = 0
  105.       NBENR2= 0
  106.  
  107. c ----boucle sur les zones du modele
  108.       DO 0001 ISOUS=1,NSOUS
  109. c         write(6,*) 'zone du modele sISOUS=',ISOUS,'/',NSOUS
  110.          IMODEL = KMODEL(ISOUS)
  111.          NOBMOD=IVAMOD(/1)
  112.          IF(NOBMOD.NE.0) THEN
  113.  
  114. c -------boucle sur les objet ivamod jusqu'a trouver un mchelm d'enrichissement
  115.          DO 0002 iobmo1=1,NOBMOD
  116.            if((TYMODE(iobmo1)).ne.'MCHAML  ')   goto 0002
  117.            MCHEX1 = IVAMOD(iobmo1)
  118.            if((MCHEX1.TITCHE).ne.'ENRICHIS')  goto 0003
  119.  
  120. *          on a trouvé une zone enrichie => on crée une table d'enrichissment nodal
  121.            N1 = MCHEX1.ICHAML(/1)
  122. c ---------boucle sur les maillages elementaire du chaml d'enrichissement
  123.            do 0010 i1=1,N1
  124.               if((MCHEX1.INFCHE(i1,4)).ne.0)  goto 0010
  125.  
  126.               MELEME = MCHEX1.IMACHE(i1)
  127.               MCHAM1 = MCHEX1.ICHAML(i1)
  128.               NBNN = NUM(/1)
  129.               NBELEM = NUM(/2)
  130.               N2=MCHAM1.IELVAL(/1)
  131.               ITYP1 = ITYPEL
  132. c--------     On exclu les sous zone de relation
  133.               if (ITYP1.EQ.48) goto 0010
  134. c ------------boucle sur les enrichissements possibles
  135.               do i2=1,N2
  136.                   MELVAL=MCHAM1.IELVAL(i2)
  137.                   do j=1,NBELEM
  138.                   do i=1,NBNN
  139.                     inoeu = NUM(i,j)
  140.                     NBPT = max(NBPT,inoeu)
  141.                     MLREEL = IELCHE(i,j)
  142. c                  write(6,*) 'i,j,inoeu,mlreel',i,j,inoeu,MLREEL,NBPT
  143.                     if (MLREEL.ne.0) then
  144.                       TLENR(i2,inoeu) = MLREEL
  145.                       NBENR2 = max(NBENR2,i2)
  146.                       if (i2.eq.1)  then
  147.                         MPHI1.VPOCHA(inoeu,i2) = PROG(i)
  148.                       else
  149.                         MPSI1.VPOCHA(inoeu,i2) = PROG(i)
  150.                         MPHI1.VPOCHA(inoeu,i2) = PROG(NBNN+i)
  151.                       endif
  152. c                     avec MPSI1 et MPHI1, on peut fermer MLREEL
  153.                     endif
  154.                   enddo
  155.                   enddo
  156.               enddo
  157. c ------------fin de boucle sur les enrichissements possibles
  158.  0010      continue
  159. c ---------fin de boucle sur les maillages elementaire du chaml d'enrichissement
  160.  
  161.  0003    CONTINUE
  162.  0002    CONTINUE
  163. c -------fin de boucle sur les objet ivamod
  164.  
  165.          ENDIF
  166.  0001 CONTINUE
  167. c ----fin de boucle sur les zones du modele
  168.  
  169.       N=NBPT
  170. *      NC=NBENR1
  171.       NC=NBENR2
  172.       segadj,MPHI1,MPSI1
  173.  
  174. C_____________________________________________________________
  175. C  INITIALISATION
  176. c
  177. c  ...du mchpo2
  178.       NAT=JATTRI(/1)
  179.       NSOUPO=NSOUP0
  180.       SEGINI,MCHPO2=MCHPOI
  181. C      NPOCHA=0
  182.  
  183.  
  184. C_____________________________________________________________
  185. C>>>>>>>>>>>>>>>>>>>>>>>  BOUCLE SUR LES ZONES DU CHPOINT
  186.       DO 1000 ISOUP=1,NSOUP0
  187. c        write(*,*) 'xpost1: ZONE ',ISOUP,' / ',NSOUP0
  188.  
  189. C_____________________________________
  190. C  ACTIVATION DU SOUS CHPOINT
  191.         MSOUPO = IPCHP(ISOUP)
  192.         MELEME = IGEOC
  193.         MPOVAL = IPOVAL
  194.  
  195. c nbre de composante, de points
  196.         NCOMP = NOCOMP(/2)
  197.         NBPT0 = NUM(/2)
  198.         NC    = NCOMP
  199.  
  200. C  INITIALISATION du MSOUP2 (avec NC =celui de départ)
  201.         SEGINI,MSOUP2=MSOUPO
  202.         MCHPO2.IPCHP(ISOUP)=MSOUP2
  203.  
  204. c quelles sont les composantes obligatoires (=physiques) et ou sont elles?
  205. c on en deduit NC
  206.         NC=0
  207.         DO 1001 IOBL=1,NOBL
  208.            MODDL2 = DDLOBL(IOBL)
  209.            DO ICOMP=1,NCOMP
  210.              MODDL = NOCOMP(ICOMP)
  211. *     on a trouvé cette comp obl dans le chpoint d entree
  212.              IF(MODDL.eq.MODDL2) THEN
  213.                NC=NC+1
  214.                MSOUP2.NOCOMP(NC) = MODDL2
  215.                TOBL(NC) = ICOMP
  216.                GOTO 1001
  217.              ENDIF
  218.           ENDDO
  219.  1001   CONTINUE
  220.         IF(NC.lt.NOBL) THEN
  221.         DO IOBL=(NC+1),NOBL
  222.             TOBL(IOBL) = 0
  223.         ENDDO
  224.         ENDIF
  225.  
  226. c ...et les facultatives(enrichissement)?
  227.         do i1=1,NOBL
  228.         do i2=1,4
  229.         do i3=1,NBRMAX
  230.            TFAC(i1,i2,i3) = 0
  231.         enddo
  232.         enddo
  233.         enddo
  234.         NF=0
  235.         IFAC = 0
  236.         DO 1002 IFAC=1,NFAC
  237.            MODDL2 = DDLFAC(IFAC)
  238.            DO ICOMP=1,NCOMP
  239.               MODDL = NOCOMP(ICOMP)
  240. *    on a trouvé une comp fac qui va etre ajouté a la comp obl dans le chpoint de sortie
  241.               IF(MODDL.eq.MODDL2) THEN
  242.                 NF=NF+1
  243.                 IOBL=TF2O(IFAC)
  244.                 INI =TNI (IFAC)
  245.                 IENR=TENR(IFAC)
  246.                 TFAC(IOBL,INI,IENR) = ICOMP
  247.                 TIFAC(IOBL,INI,IENR) = IFAC
  248.                 GOTO 1002
  249.               ENDIF
  250.            ENDDO
  251. *     on n a pas trouvé la composante facultative
  252.                 IOBL=TF2O(IFAC)
  253.                 INI =TNI (IFAC)
  254.                 IENR=TENR(IFAC)
  255.                 TIFAC(IOBL,INI,IENR) = IFAC
  256.  1002   CONTINUE
  257.  
  258.  
  259. C  INITIALISATION du MSOUP2 (avec NC ajusté) et du MPOVA2
  260.         segadj,MSOUP2
  261.         N = VPOCHA(/1)
  262.         SEGINI,MPOVA2
  263.         MSOUP2.IPOVAL = MPOVA2
  264.  
  265. c       actifs a ce stade : MCHPO2,MCHPOI,MELEME,MPOVAL,MPOVA2
  266.  
  267. C_____________________________________
  268. C>>>>>>>>>>>>>  BOUCLE SUR LES POINTS
  269.         DO 2000 J=1,NBPT0
  270. c          write(*,*) 'point ',J,' / ',NBPT0,NC
  271.  
  272. C______________________
  273. C         On commence par recopier les valeurs obligatoires (UX,UY)
  274.           DO IC2=1,NC
  275.               ICOMP = TOBL(IC2)
  276.               MPOVA2.VPOCHA(J,IC2) = VPOCHA(J,ICOMP)
  277.           ENDDO
  278.  
  279. C______________________
  280. C         ce noeud est il IENR2-enrichi?
  281.           inoeu = NUM(1,J)
  282.  
  283.           DO 3000 IENR2=1,NBENR2
  284.  
  285.              MLREEL = TLENR(IENR2,inoeu)
  286.  
  287. C            si ce noeud n est pas enrichi on ne fait rien
  288.              IF(MLREEL.eq.0)  GOTO 3000
  289.  
  290. C            si ce noeud est enrichi,
  291. c            on recupere les infos relatives a l enrichissement
  292.  
  293. c            on calcule les fonctions d enrichissement
  294.  
  295. c------------pour IENR=1, fonction H, ddl AX et AY
  296.              IF (IENR2.eq.1) THEN
  297.               PHIX = MPHI1.VPOCHA(inoeu,IENR2)
  298.               NBNI = 1
  299.               SHPWRK(1) =  SIGN(1.D0,PHIX)
  300.               if (abs(phix).lt.xtol1) then
  301.                SHPWRK(1) =  0.d0
  302.               else
  303.                SHPWRK(1) =  SIGN(1.D0,PHIX)
  304.               endif
  305.              ENDIF
  306. c------------fin du cas IENR=1, fonction H
  307.  
  308. c------------pour IENR>1, 4 fonctions F
  309.              IF (IENR2.ge.2) THEN
  310.               PHIX = MPHI1.VPOCHA(inoeu,IENR2)
  311.               PSIX = MPSI1.VPOCHA(inoeu,IENR2)
  312.               if (ABS(phix).lt.xtol1) then
  313.                HX = 0.d0
  314.               else
  315.                HX = SIGN(1.D0,PHIX)
  316.               endif
  317.               RX   = ( (PSIX**2.) + (PHIX**2.) )**0.5
  318.               IF (RX.LT.XTOL1) THEN
  319.                 THETAX = 0.d0
  320.               ELSE
  321.                 THETAX = HX * ((XPI/2.) - (ATAN2(PSIX,ABS(PHIX))))
  322.               ENDIF
  323.               SIN1T = SIN(THETAX)
  324. C              COS1T = COS(THETAX)
  325.               SIN05T = SIN(THETAX/2.)
  326.               COS05T = COS(THETAX/2.)
  327.               NBNI = 4
  328.               SHPWRK(1) = (RX**0.5) * SIN05T
  329.               SHPWRK(2) = (RX**0.5) * COS05T
  330.               SHPWRK(3) = (RX**0.5) * SIN05T * SIN1T
  331.               SHPWRK(4) = (RX**0.5) * COS05T * SIN1T
  332.              ENDIF
  333. c------------fin du cas IENR>1, fonction F
  334.  
  335. c            on ajoute les fonctions d enrichissement
  336.              DO 3900 IC2=1,NC
  337.              DO 3900 INI=1,NBNI
  338.               ICOMP = TFAC(IC2,INI,IENR2)
  339.               IFAC = TIFAC(IC2,INI,IENR2)
  340. c             1ere fois qu on voit le besoin de (composante + noeud)
  341.               if(IPOCHA(inoeu,IFAC).eq.0) IPOCHA(inoeu,IFAC)=1
  342. c             cas ou on a pas trouvé cette composante dans cette zone du
  343. c             chpoint solution => on saute simplement
  344.               if(ICOMP.eq.0) goto 3900
  345.               MPOVA2.VPOCHA(J,IC2) = MPOVA2.VPOCHA(J,IC2)
  346.      $                + ( SHPWRK(INI) * VPOCHA(J,ICOMP) )
  347.               IPOCHA(inoeu,IFAC) = IPOCHA(inoeu,IFAC)+1
  348.  3900        CONTINUE
  349.  
  350.  3000     CONTINUE
  351. C<<<<<<<<< FIN DE BOUCLE SUR LES enrichissements
  352.  
  353.  
  354.  
  355.  2000   CONTINUE
  356. C<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES POINTS
  357.  
  358.  
  359.  1000 CONTINUE
  360. C<<<<<<<<<<<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES ZONES
  361.  
  362.       segsup,MPHI1,MPSI1
  363.  
  364.  
  365. C_____________________________________________________________
  366. C     PETIT AVERTISSEMENT S'IL MANQUE DES ENRICHISSEMENTS
  367. C
  368. C   --initialisation du nombre d erreur sur le noms de composantes
  369.       NBERR1=0
  370.  
  371. c   --recherche des erreurs
  372.       DO inoeu=1,NBPT
  373.       DO icomp1=1,NCOMP1
  374. c       -cas ou on a pas besoin de cette composante en ce noeud :
  375. c        (IPOCHA(inoeu,icomp1).EQ.0)
  376. c       -Cas ou on a trouvé un noeud enrichi sans la composante associée
  377. c        dans le chpoint => avertissement :
  378.          IF (IPOCHA(inoeu,icomp1).EQ.1) THEN
  379.            NBERR1=NBERR1+1
  380.            if (IIMPI.ge.1) then
  381.              write(IOIMP,991) DDLFAC(icomp1),inoeu
  382.   991        format(2X,'ABSENCE DE LA COMPOSANTE ',A4,' AU NOEUD ',I6,
  383.      $            ' DANS LE CHPOINT FOURNI a XFEM RECO')
  384.            endif
  385.          ENDIF
  386.       ENDDO
  387.       ENDDO
  388. c   --cas ou on a une ou des erreurs--
  389.       IF (NBERR1.gt.0) THEN
  390.         write(IOIMP,*) 'OPERATEUR XFEM RECO : ABSENCE DANS LE CHPOINT ',
  391.      &   'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE'
  392.       ENDIF
  393.  
  394.       segsup,MRACC
  395.  
  396.  
  397. C_____________________________________________________________
  398. C     ON RETOURNE LE CHPOINT
  399.  
  400.       IPCHP2=MCHPO2
  401.  
  402.       END
  403.  
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  
  411.  

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