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

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