Télécharger la plaquette Cast3M


fr en

Télécharger xpost2.eso

Retour à la liste

Numérotation des lignes : 

xpost2

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

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