Télécharger la plaquette Cast3M


fr en

Télécharger xpost3.eso

Retour à la liste

Numérotation des lignes : 

xpost3

  1. C XPOST3    SOURCE    CB215821  24/04/12    21:17:31     11897          
  2. c
  3. C XPOST2    SOURCE    CB215821  19/07/30    21:18:43     10273          
  4.       SUBROUTINE XPOST3(IPCHP1,IPMOD1,IPT2in,IPCHP2)
  5. c
  6. c     Construit 2 Chpoints avec des ddl "physiques" (UX UY)
  7. c     sur un maillage quelconque d'interpolation
  8. c     en RECOmbinant les ddl Xfem (UX AX B1X C1X ...)
  9. c     -> Utile pour projeter
  10. c
  11. C     SYNTAXE :
  12. c       chpo2 chpo3 = XFEM 'RECO' chpo1 MODEL_XFEM geo2
  13. C
  14. c
  15.       IMPLICIT INTEGER(I-N)
  16.       IMPLICIT REAL*8 (A-H,O-Z)
  17. C
  18.       PARAMETER   (NBRMAX=5)
  19.       PARAMETER   (XTOL=1.D-5)
  20. C
  21. C  SEGMENTS INCLUDE
  22.  
  23. -INC PPARAM
  24. -INC CCOPTIO
  25. -INC CCREEL
  26. -INC SMCOORD
  27. -INC SMELEME
  28. -INC SMCHPOI
  29. -INC SMCHAML
  30. -INC SMMODEL
  31. -INC SMLREEL
  32. c
  33.       POINTEUR    MCHEX1.MCHELM
  34. C
  35.       SEGMENT IDEJVU(NBPT2)
  36. c
  37.       segment wrk4
  38.         real*8 xel(3,nbnn),shpp(6,nbnn),qsi(3),XPT2(3)
  39.       endsegment
  40. C
  41. C_____________________________________________________________
  42. C  INITIALISATION DES INCONNUES obligatoires et facultatives
  43.       PARAMETER (NOBL=3,NFAC=27)
  44.       CHARACTER*(LOCOMP) DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2
  45.       DATA DDLOBL/'UX  ','UY  ','UZ  '/
  46.       DATA DDLFAC/'AX  ','AY  ','AZ  ',
  47.      >            'B1X ','B1Y ','B1Z ',
  48.      >            'C1X ','C1Y ','C1Z ',
  49.      >            'D1X ','D1Y ','D1Z ',
  50.      >            'E1X ','E1Y ','E1Z ',
  51.      >            'B2X ','B2Y ','B2Z ',
  52.      >            'C2X ','C2Y ','C2Z ',
  53.      >            'D2X ','D2Y ','D2Z ',
  54.      >            'E2X ','E2Y ','E2Z '/
  55.       INTEGER TENR(NFAC),TNI(NFAC),TF2O(NFAC)
  56. c     ddlfac correspond a quel enrichissement?
  57.       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/
  58. c     ddlfac correspond a quel fonction d'enrichissement?
  59.       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/
  60. c     ddlfac correspond a quel ddlobl?
  61.       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/
  62. c     tables pour mise en concordance des ddl
  63.       INTEGER TOBL(NOBL),TFAC(NOBL,4,NBRMAX)
  64.       INTEGER           TIFAC(NOBL,4,NBRMAX)
  65. c     fonctions de forme
  66.       REAL*8  SHPWRK(4),SHPWR2(4)
  67.  
  68. C_____________________________________________________________
  69. C
  70.       NBERR1=0
  71.  
  72. C  LECTURE 
  73. c       CALL LIROBJ('MAILLAGE',IPT2in,1,IRETOU)
  74. c       CALL LIROBJ('CHPOINT ',MCHPOI,1,IRETOU)
  75. c       CALL LIROBJ('MMODEL  ',MMODEL,1,IRETOU)
  76. c       CALL ACTOBJ('MAILLAGE',IPT2in,1)
  77. c       CALL ACTOBJ('CHPOINT ',MCHPOI,1)
  78. c       CALL ACTOBJ('MMODEL  ',MMODEL,1)
  79. c       IF(IERR.NE.0) RETURN
  80. C  LECTURE faite par xpost
  81.       MCHPOI=IPCHP1
  82.       MMODEL=IPMOD1
  83.       IPT2  = IPT2in      
  84.  
  85.       SEGACT,MCOORD
  86.  
  87. C_____________________________________________________________
  88. C     TRANSFO EN MCHELM DU CHPOINT avec ddl xfem
  89.       call CHAME1(0,MMODEL,MCHPOI,' ',IPCHEL,1)
  90.  
  91. C_____________________________________________________________
  92. C  CHANGEMENT DE SUPPORT POUR LE CHPOINT DE SORTIE
  93.       call change(ipt2,1)
  94. c       write(6,*) 'xpost2: retour de change (ipt2 doit etre actif)',ipt2
  95.  
  96. C_____________________________________________________________
  97. C  ACTIVATION DU MODELE
  98.       NSOUS = KMODEL(/1)
  99.  
  100. C_____________________________________________________________
  101. c       write(6,*) 'C  INITIALISATION du mchpo2'
  102. C
  103.       NAT=JATTRI(/1)
  104.       NSOUPO=1
  105.       SEGINI,MCHPO2=MCHPOI
  106.       segadj,MCHPO2
  107.       SEGINI,MCHPO3=MCHPO2
  108. c     on sait que l'on a seulement besoin de UX et UY en 2D
  109. c     (+léger que dans xpost1.eso)
  110.       NC=IDIM
  111.       segini,MSOUP2,MSOUP3
  112.       MCHPO2.IPCHP(1) = MSOUP2
  113.       MCHPO3.IPCHP(1) = MSOUP3
  114.       do ic2=1,IDIM
  115.          MSOUP2.NOCOMP(ic2) = DDLOBL(ic2)
  116.          MSOUP3.NOCOMP(ic2) = DDLOBL(ic2)
  117.          MSOUP2.NOHARM(ic2) = NIFOUR
  118.          MSOUP3.NOHARM(ic2) = NIFOUR
  119.       enddo
  120.       MSOUP2.IGEOC = IPT2
  121.       MSOUP3.IGEOC = IPT2
  122.       NBNN2 = IPT2.NUM(/1)
  123.       NBELE2 = IPT2.NUM(/2)
  124.       NBPT2 = NBELE2
  125.       N = NBPT2
  126. C       segini,MPOVA2,MPOVA3
  127. C     bp: on ajoute IDEJVU pour ne caluler le champ de depl qu'une seule fois si le point arive sur le bord d un element
  128.       segini,MPOVA2,MPOVA3,IDEJVU
  129.       MSOUP2.IPOVAL = MPOVA2
  130.       MSOUP3.IPOVAL = MPOVA3
  131.  
  132. C     CHAMP DE DEPLACEMENT X-FEM
  133.       MCHELM = IPCHEL
  134.  
  135. C_____________________________________________________________
  136. C>>>>>>>>>>>>>>>>>>>>>>>  BOUCLE SUR LES ZONES DU Modele
  137.       DO 1000 ISOUS=1,NSOUS
  138.         IMODEL = KMODEL(ISOUS)
  139. c         write(6,*) '____________________________________'
  140. c         write(6,*) 'xpost2: ZONE',ISOUS,' / ',NSOUS,'  IMODEL=',IMODEL
  141.  
  142. C_____________________________________
  143. C  RECHERCHE DU MCHAM1 issu du MCHEX1 D ENRICHISSEMENT
  144.         MCHAM1=0
  145.         NBENR2=0
  146.         isouX=0
  147.         NOBMOD=IVAMOD(/1)
  148. c         write(6,*) 'IMAMOD,NEFMOD,NOBMOD',IMAMOD,NEFMOD,NOBMOD
  149. C Si sous modele de sure on passe au sous modele suivant
  150.         if (NEFMOD.eq.259) goto 1000
  151.         IF (NOBMOD.NE.0) THEN
  152. c        boucle sur les objets contenus dans ivamod
  153.          DO 1002 iobmo1=1,NOBMOD
  154. c            write(6,*) '-->1002: objet=',iobmo1,'/',NOBMOD
  155. c            write(*,*) '   TYMODE(iobmo1)=',TYMODE(iobmo1)
  156.            if((TYMODE(iobmo1)).ne.'MCHAML')   goto 1002
  157.            MCHEX1 = IVAMOD(iobmo1)
  158. c            write(*,*) '   MCHEX1,MCHEX1.TITCHE=',MCHEX1,MCHEX1.TITCHE
  159.            if((MCHEX1.TITCHE).ne.'ENRICHIS')  goto 1003
  160. c REM. IMPORTANTE: a ce jour 1 seule zone XFEM, donc on prend ICHAML(isouX=1)
  161. c mais on ecrit ci dessous un prototype de la recherche de la zone du MCHEX1
  162. c associée a la zone du sous modele en cours via le pointeur du maillage
  163.            nsouX = MCHEX1.ICHAML(/1)
  164.            do isouX=1,nsouX
  165. c               write(*,*) 'isouX=',isouX,'/',nsouX,' : ',IMACHE(isouX)
  166.               if(IMAMOD.eq.MCHEX1.IMACHE(isouX)) goto 1004
  167.            enddo
  168.            isouX=0
  169.            goto 1003
  170.  1004      continue
  171.            MCHAM1 = MCHEX1.ICHAML(isouX)
  172.            NBENR2 = MCHAM1.IELVAL(/1)
  173.  1003      continue
  174.  1002    CONTINUE
  175.         ENDIF
  176. c         if(isouX.ne.0) then
  177. c           write(6,*) 'isouX/n',isouX,nsouX
  178. c           write(6,*) 'IMAMOD,IMACHE=',IMAMOD,IMACHE(isouX)
  179. c         endif
  180. c         write(6,*) 'MCHAM1,NBENR2=',MCHAM1,NBENR2
  181.  
  182. C_____________________________________
  183. C  ACTIVATION DU SOUS CHPOINT devenu MCHELM  ET DES COMPOSANTES MELVAL
  184.         MCHAML = ICHAML(ISOUS)
  185.         N2= IELVAL(/1)
  186.  
  187. c nbre de composante, de points
  188.         NCOMP = N2
  189.  
  190.  
  191. C_____________________________________
  192. c quelles sont les composantes obligatoires (=physiques) et ou sont elles?
  193. c on en deduit NC
  194.         NC=0
  195.         DO 1011 IOBL=1,NOBL
  196.            MODDL2 = DDLOBL(IOBL)
  197.            DO ICOMP=1,NCOMP
  198.              MODDL = NOMCHE(ICOMP)
  199. *     on a trouvé cette comp obl dans le chpoint d entree
  200.              IF (MODDL.eq.MODDL2) THEN
  201.                NC=NC+1
  202.                TOBL(NC) = ICOMP
  203.                GOTO 1011
  204.              ENDIF
  205.           ENDDO
  206.  1011   CONTINUE
  207.         IF (NC.lt.NOBL) THEN
  208.         DO IOBL=(NC+1),NOBL
  209.             TOBL(IOBL) = 0
  210.         ENDDO
  211.         ENDIF
  212.  
  213. c ...et les facultatives(enrichissement)?
  214.         do i1=1,NOBL
  215.         do i2=1,4
  216.         do i3=1,NBRMAX
  217.            TFAC(i1,i2,i3) = 0
  218.         enddo
  219.         enddo
  220.         enddo
  221.         NF=0
  222.         IFAC = 0
  223.         DO 1012 IFAC=1,NFAC
  224.            MODDL2 = DDLFAC(IFAC)
  225.            DO ICOMP=1,NCOMP
  226.               MODDL = NOMCHE(ICOMP)
  227. *    on a trouvé une comp fac qui va etre ajouté a la comp obl dans le chpoint de sortie
  228.               IF (MODDL.eq.MODDL2) THEN
  229.                 NF=NF+1
  230.                 IOBL=TF2O(IFAC)
  231.                 INI =TNI (IFAC)
  232.                 IENR=TENR(IFAC)
  233.                 TFAC(IOBL,INI,IENR) = ICOMP
  234.                 TIFAC(IOBL,INI,IENR) = IFAC
  235.                 GOTO 1012
  236.               ENDIF
  237.            ENDDO
  238. *     on n a pas trouvé la composante facultative
  239.                 IOBL=TF2O(IFAC)
  240.                 INI =TNI (IFAC)
  241.                 IENR=TENR(IFAC)
  242.                 TIFAC(IOBL,INI,IENR) = IFAC
  243.  1012   CONTINUE
  244. c         write(*,*) 'TOBL et TFAC remplis'
  245.  
  246. c       recup du maillage
  247.         MELEME = IMAMOD
  248.         NBNN0=NUM(/1)
  249.         NBELE0=NUM(/2)
  250.         MELGEO = ITYPEL
  251.  
  252. c       initialisation du segment de travail
  253.         NBNN = NBNN0
  254.         segini,wrk4
  255.  
  256.  
  257. C_____________________________________
  258. C>>>>>>>>>>>>>  BOUCLE SUR LES ELEMENTS DU MODELE
  259.         DO 2000 J0=1,NBELE0
  260.  
  261. c         XEL = coor des noeuds de l'ef J0 de MELEME
  262.           call DOXE(XCOOR,idim,NBNN0,NUM,J0,XEL)
  263.  
  264.  
  265. C_____________________________________
  266. C>>>>>>>>>>>>>  BOUCLE SUR LES NOEUD DU CHP2
  267.           DO 2002 J2=1,NBPT2
  268.  
  269. *           point J2 IDEJaVU dans un autre element => on saute
  270.             if(IDEJVU(J2).eq.1) goto 2002
  271.  
  272. C            inoeu2 = IPT2.NUM(1,J2)
  273. c           XPT2 = coor du pt J2 de IPT2
  274.             call DOXE(XCOOR,idim,1,IPT2.NUM,J2,XPT2)
  275.  
  276. c calcul iteratif des coord QSI du pt XPT2 dans le repere local de l'ef defini par XEL
  277.             call qsijs(XEL,MELGEO,NBNN,idim,XPT2,SHPP,QSI,iret)
  278.             if(iret.ne.0) goto 2002
  279.             do i0=1,NBNN0
  280. *              if (SHPP(1,i0).lt.0.) goto 2002
  281. *               a revoir car pb lorsqu'on tombe "pile" sur un bord
  282.               if(SHPP(1,i0).lt.(-1.E-7)) goto 2002
  283.             enddo
  284.  
  285. *           On a trouvé que le point J2 appartenait a l'EF J0
  286.             IDEJVU(J2) = 1
  287. c             write(6,*) '_______________________________'
  288. c             write(*,*) 'ELEM ',J0,'/',NBELE0,'  point ',J2,'/',NBPT2
  289. c             write(6,*) 'QSI= ',(QSI(iou),iou=1,3)
  290. c             write(6,*) 'SHPP=',(SHPP(1,iou),iou=1,4)
  291.  
  292. C______________________
  293. C         On commence par interpoler la partie relative aux ddls (UX,UY)
  294.           DO IC2=1,NC
  295.              MPOVA2.VPOCHA(J2,IC2) = 0.d0
  296.              ICOMP = TOBL(IC2)
  297.              MELVAL = IELVAL(ICOMP)
  298.              DO I0=1,NBNN0
  299.                MPOVA2.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2)
  300.      $                      + ( SHPP(1,I0) * VELCHE(I0,J0) )
  301.              ENDDO
  302.              MPOVA3.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2)
  303.           ENDDO
  304.  
  305. C______________________
  306. C         le noeud inoeu (I0^ieme noeud de l ef J0) est il IENR2-enrichi?
  307.           if(NBENR2.eq.0) goto 2002
  308.           DO 3000 I0=1,NBNN0
  309.           inoeu = NUM(I0,J0)
  310.  
  311.           DO 3001 IENR2=1,NBENR2
  312.  
  313.              MELVA1=MCHAM1.IELVAL(IENR2)
  314.              MLREEL=MELVA1.IELCHE(I0,J0)
  315.  
  316. C            si ce noeud n est pas IENR2-enrichi on ne fait rien
  317.              IF(MLREEL.eq.0)  GOTO 3001
  318.  
  319. C            si ce noeud est enrichi,...
  320. c                write(6,*) 'I0,inoeu,IENR2,MLREEL=',
  321. c      $                     I0,inoeu,IENR2,MLREEL
  322.  
  323. c            ...on calcule les fonctions d enrichissement :
  324.  
  325. c------------pour IENR=1, fonction H, ddl AX et AY
  326.              IF (IENR2.eq.1) THEN
  327.                NBNI = 1              
  328. c               SHPWRK(1) =  1.d0
  329. c               SHPWR2(1) =  -1.d0 
  330.                 PHIX = 0.d0
  331.                 PHIMAX = 1.D-16
  332.                 do iii0 = 1,NBNN0
  333.                   XSHPP = SHPP(1,iii0)
  334.                   XPHI1 = PROG(iii0)
  335.                   PHIX  = PHIX + (XSHPP * XPHI1)
  336.                   PHIMAX=MAX(PHIMAX,ABS(XPHI1))
  337.                 enddo
  338.                 IF(ABS(PHIX).LE.XTOL*PHIMAX) THEN
  339.                   SHPWRK(1) =  1.d0
  340.                   SHPWR2(1) =  -1.d0 
  341.                 ELSE
  342.                   HX = SIGN(1.D0,PHIX)
  343.                   SHPWRK(1) =  HX
  344.                   SHPWR2(1) =  HX               
  345.                 ENDIF
  346.              ENDIF
  347. c------------fin du cas IENR=1, fonction H
  348.  
  349. c------------pour IENR>1, 4 fonctions F
  350.              IF (IENR2.ge.2) THEN
  351.               NBNI = 4
  352. c             psi = \sum Ni(x) \psi_i
  353. c             phi = \sum Ni(x) \phi_i
  354.               PSIX = 0.d0
  355.               PHIX = 0.d0
  356.               PHIMAX = 1.D-16
  357. c               WRITE(6,*) 'XPSI1=',(PROG(IOU),IOU=1,4)
  358.               do iii0 = 1,NBNN0
  359.                XSHPP = SHPP(1,iii0)
  360.                XPSI1 = PROG(iii0)
  361.                XPHI1 = PROG(NBNN0+iii0)
  362.                PSIX  = PSIX + (XSHPP * XPSI1)
  363.                PHIX  = PHIX + (XSHPP * XPHI1)
  364.                PHIMAX=MAX(PHIMAX,ABS(XPHI1))
  365.               enddo
  366. *              write(*,*) 'PHIX,PSIX=',PHIX,PSIX
  367.               RX   = ( (PSIX**2) + (PHIX**2) )**0.5
  368.               IF(RX.LT.XTOL*PHIMAX) THEN
  369.                 THETAX = 0.
  370.                 THETA2 = 0.
  371.               ELSEIF(ABS(PHIX).LE.XTOL*PHIMAX .AND. PSIX.LE.0.) THEN
  372.                 THETAX  = XPI
  373.                 THETA2 = -1.*XPI
  374.               ELSE
  375.                 HX = SIGN(1.D0,PHIX)
  376.                 THETAX = HX * ((XPI/2.) - (ATAN(PSIX/ABS(PHIX))))
  377.                 THETA2 = THETAX
  378.               ENDIF
  379.               RX05= (RX**0.5)
  380. c               write(6,*) 'PSIX,PHIX,RX05=',PSIX,PHIX,RX05
  381.               SIN1T = SIN(THETAX)
  382. C              COS1T = COS(THETAX)
  383.               SIN05T = SIN(THETAX/2.)
  384.               COS05T = COS(THETAX/2.)
  385. c               if (IENR2.lt.2) then
  386. c                SHPWRK(1) = 0.d0
  387. c               else
  388.               SHPWRK(1) = (RX05) * SIN05T
  389. c               endif
  390.               SHPWRK(2) = (RX05) * COS05T
  391.               SHPWRK(3) = (RX05) * SIN05T * SIN1T
  392.               SHPWRK(4) = (RX05) * COS05T * SIN1T
  393.               SIN12 = SIN(THETA2)
  394. C              COS12 = COS(THETA2)
  395.               SIN052 = SIN(THETA2/2.)
  396.               COS052 = COS(THETA2/2.)
  397. c             if (IENR2.lt.2) then
  398. c               SHPWR2(1) = 0.d0
  399. c             else
  400.               SHPWR2(1) = (RX05) * SIN052
  401. c             endif
  402.               SHPWR2(2) = (RX05) * COS052
  403.               SHPWR2(3) = (RX05) * SIN052 * SIN12
  404.               SHPWR2(4) = (RX05) * COS052 * SIN12
  405.              ENDIF
  406. c------------fin du cas IENR>1, fonction F
  407.  
  408. c            on ajoute les fonctions d enrichissement
  409.              DO 3900 IC2=1,NC
  410.              DO 3901 INI=1,NBNI
  411.               ICOMP = TFAC(IC2,INI,IENR2)
  412. c             cas ou on a pas trouvé cette composante dans cette zone du
  413. c             chpoint solution => AVERTISSEMENT puis on saute
  414.               if (ICOMP.eq.0) then
  415.                 NBERR1=NBERR1+1
  416.                 if (IIMPI.ge.1) then
  417.                  write(IOIMP,991) DDLFAC(icomp),inoeu
  418.   991            format(2X,'ABSENCE DE LA COMPOSANTE ',A4,' AU NOEUD ',
  419.      $                  I6,' DANS LE CHPOINT FOURNI a XFEM FISS')
  420.                 endif
  421.                 goto 3900
  422.               endif
  423.               MELVAL = IELVAL(ICOMP)
  424.               MPOVA2.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2)
  425.      $           + ( (SHPP(1,I0) * SHPWRK(INI)) * VELCHE(I0,J0) )
  426.               MPOVA3.VPOCHA(J2,IC2) = MPOVA3.VPOCHA(J2,IC2)
  427.      $            + ( (SHPP(1,I0) * SHPWR2(INI)) * VELCHE(I0,J0) )
  428.  3901        CONTINUE
  429.  3900        CONTINUE
  430.  
  431.  
  432.  3001     CONTINUE
  433. C<<<<<<<<< FIN DE BOUCLE SUR LES enrichissements
  434.  3000     CONTINUE
  435. C<<<<<<<<< FIN DE BOUCLE SUR LES noeuds
  436.  
  437.  
  438.  
  439.  2002   CONTINUE
  440. C<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES noeud du chp2
  441.  2000   CONTINUE
  442. C<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES elements
  443.  
  444.       segsup,wrk4
  445.  
  446.  1000 CONTINUE
  447. C<<<<<<<<<<<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES ZONES
  448.  
  449.       SEGSUP,IDEJVU
  450.  
  451. c   --cas ou on a une ou des erreurs--
  452.       IF (NBERR1.gt.0) THEN
  453.         write(IOIMP,*) 'OPERATEUR XFEM FISS : ABSENCE DANS LE CHPOINT ',
  454.      &   'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE'
  455.       ENDIF
  456.  
  457. C_____________________________________________________________
  458. C
  459. C  ECRITURE
  460.       CALL ACTOBJ('CHPOINT ',MCHPO2,1)
  461.       CALL ACTOBJ('CHPOINT ',MCHPO3,1)
  462.       CALL ECROBJ('CHPOINT ',MCHPO2)
  463.       CALL ECROBJ('CHPOINT ',MCHPO3)
  464.  
  465.       END
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  

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