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

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