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  24/08/02    21:15:04     11974          
  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.             call qsijs(XEL,MELGEO,NBNN,idim,XPT2,SHPP,QSI,iret)
  271.             if(iret.ne.0) goto 2002
  272.             do i0=1,NBNN0
  273. *              if (SHPP(1,i0).lt.0.) goto 2002
  274. *               a revoir car pb lorsqu'on tombe "pile" sur un bord
  275.               if(SHPP(1,i0).lt.(-1.E-7)) goto 2002
  276.             enddo
  277.  
  278. *           On a trouvé que le point J2 appartenait a l'EF J0
  279.             IDEJVU(J2) = 1
  280. c             write(6,*) '_______________________________'
  281. c             write(*,*) 'ELEM ',J0,'/',NBELE0,'  point ',J2,'/',NBPT2
  282. c             write(6,*) 'QSI= ',(QSI(iou),iou=1,3)
  283. c             write(6,*) 'SHPP=',(SHPP(1,iou),iou=1,4)
  284.  
  285. C______________________
  286. C         On commence par interpoler la partie relative aux ddls (UX,UY)
  287.           DO IC2=1,NC
  288.              MPOVA2.VPOCHA(J2,IC2) = 0.d0
  289.              ICOMP = TOBL(IC2)
  290.              MELVAL = IELVAL(ICOMP)
  291.              DO I0=1,NBNN0
  292.                MPOVA2.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2)
  293.      $                      + ( SHPP(1,I0) * VELCHE(I0,J0) )
  294.              ENDDO
  295.              MPOVA3.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2)
  296. c              write(*,*) 'U_',IC2,'=',MPOVA2.VPOCHA(J2,IC2)
  297.           ENDDO
  298.  
  299. C______________________
  300. C         le noeud inoeu (I0^ieme noeud de l ef J0) est il IENR2-enrichi?
  301.           if(NBENR2.eq.0) goto 2002
  302.           DO 3000 I0=1,NBNN0
  303.           inoeu = NUM(I0,J0)
  304.  
  305.           DO 3001 IENR2=1,NBENR2
  306.  
  307.              MELVA1=MCHAM1.IELVAL(IENR2)
  308.              MLREEL=MELVA1.IELCHE(I0,J0)
  309.  
  310. C            si ce noeud n est pas IENR2-enrichi on ne fait rien
  311.              IF(MLREEL.eq.0)  GOTO 3001
  312.  
  313. C            si ce noeud est enrichi,...
  314. c                write(6,*) 'I0,inoeu,IENR2,MLREEL=',
  315. c      $                     I0,inoeu,IENR2,MLREEL
  316.  
  317. c            ...on calcule les fonctions d enrichissement :
  318.  
  319. c------------pour IENR=1, fonction H, ddl AX et AY
  320.              IF (IENR2.eq.1) THEN
  321.                NBNI = 1
  322.               SHPWRK(1) =  1.d0
  323.               SHPWR2(1) =  -1.d0
  324.              ENDIF
  325. c------------fin du cas IENR=1, fonction H
  326.  
  327. c------------pour IENR>1, 4 fonctions F
  328.              IF (IENR2.ge.2) THEN
  329.               NBNI = 4
  330.               PHIX = 0.d0
  331.               PSIX = 0.d0
  332. c               WRITE(6,*) 'XPSI1=',(PROG(IOU),IOU=1,4)
  333.               do iii0 = 1,NBNN0
  334.                XSHPP = SHPP(1,iii0)
  335.                XPSI1 = PROG(iii0)
  336.                PSIX  = PSIX + (XSHPP * XPSI1)
  337.               enddo
  338. *              write(6,*) 'PHIX,PSIX=',PHIX,PSIX
  339.               RX   = ( (PSIX**2) + (PHIX**2) )**0.5
  340. *              IF(RX.LT.XTOL) THEN
  341.               IF (PSIX.ge.0.) THEN
  342.                 THETAX = 0.
  343.                 THETA2 = 0.
  344.               ELSE
  345. *                THETAX = HX * ((XPI/2.) - (ATAN(PSIX/ABS(PHIX))))
  346.                 THETAX  = XPI
  347.                 THETA2 = -1.*XPI
  348.               ENDIF
  349.               RX05= (RX**0.5)
  350. c               write(6,*) 'PSIX,PHIX,RX05=',PSIX,PHIX,RX05
  351.               SIN1T = SIN(THETAX)
  352. C              COS1T = COS(THETAX)
  353.               SIN05T = SIN(THETAX/2.)
  354.               COS05T = COS(THETAX/2.)
  355. c               if (IENR2.lt.2) then
  356. c                SHPWRK(1) = 0.d0
  357. c               else
  358.               SHPWRK(1) = (RX05) * SIN05T
  359. c               endif
  360.               SHPWRK(2) = (RX05) * COS05T
  361.               SHPWRK(3) = (RX05) * SIN05T * SIN1T
  362.               SHPWRK(4) = (RX05) * COS05T * SIN1T
  363.               SIN12 = SIN(THETA2)
  364. C              COS12 = COS(THETA2)
  365.               SIN052 = SIN(THETA2/2.)
  366.               COS052 = COS(THETA2/2.)
  367. c             if (IENR2.lt.2) then
  368. c               SHPWR2(1) = 0.d0
  369. c             else
  370.               SHPWR2(1) = (RX05) * SIN052
  371. c             endif
  372.               SHPWR2(2) = (RX05) * COS052
  373.               SHPWR2(3) = (RX05) * SIN052 * SIN12
  374.               SHPWR2(4) = (RX05) * COS052 * SIN12
  375.              ENDIF
  376. c------------fin du cas IENR>1, fonction F
  377.  
  378. c            on ajoute les fonctions d enrichissement
  379.              DO 3900 IC2=1,NC
  380.              DO 3901 INI=1,NBNI
  381.               ICOMP = TFAC(IC2,INI,IENR2)
  382. c             cas ou on a pas trouvé cette composante dans cette zone du
  383. c             chpoint solution => AVERTISSEMENT puis on saute
  384.               if (ICOMP.eq.0) then
  385.                 NBERR1=NBERR1+1
  386.                 if (IIMPI.ge.1) then
  387.                  write(IOIMP,991) DDLFAC(icomp),inoeu
  388.   991            format(2X,'ABSENCE DE LA COMPOSANTE ',A4,' AU NOEUD ',
  389.      $                  I6,' DANS LE CHPOINT FOURNI a XFEM FISS')
  390.                 endif
  391.                 goto 3900
  392.               endif
  393.               MELVAL = IELVAL(ICOMP)
  394.               MPOVA2.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2)
  395.      $           + ( (SHPP(1,I0) * SHPWRK(INI)) * VELCHE(I0,J0) )
  396.               MPOVA3.VPOCHA(J2,IC2) = MPOVA3.VPOCHA(J2,IC2)
  397.      $            + ( (SHPP(1,I0) * SHPWR2(INI)) * VELCHE(I0,J0) )
  398.  3901        CONTINUE
  399.  3900        CONTINUE
  400.  
  401.  
  402.  3001     CONTINUE
  403. C<<<<<<<<< FIN DE BOUCLE SUR LES enrichissements
  404.  3000     CONTINUE
  405. C<<<<<<<<< FIN DE BOUCLE SUR LES noeuds
  406.  
  407. c           do ic2=1,NC
  408. c           write(*,*) 'U_tot',IC2,'=',MPOVA2.VPOCHA(J2,IC2)
  409. c           enddo
  410.  
  411.  
  412.  2002   CONTINUE
  413. C<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES noeud du chp2
  414.  2000   CONTINUE
  415. C<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES elements
  416.  
  417.       segsup,wrk4
  418.  
  419.  1000 CONTINUE
  420. C<<<<<<<<<<<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES ZONES
  421.  
  422.       SEGSUP,IDEJVU
  423.  
  424. c   --cas ou on a une ou des erreurs--
  425.       IF (NBERR1.gt.0) THEN
  426.         write(IOIMP,*) 'OPERATEUR XFEM FISS : ABSENCE DANS LE CHPOINT ',
  427.      &   'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE'
  428.       ENDIF
  429.  
  430. C_____________________________________________________________
  431. C
  432. C  ECRITURE
  433.       CALL ACTOBJ('CHPOINT ',MCHPO2,1)
  434.       CALL ACTOBJ('CHPOINT ',MCHPO3,1)
  435.       CALL ECROBJ('CHPOINT ',MCHPO2)
  436.       CALL ECROBJ('CHPOINT ',MCHPO3)
  437.  
  438.       END
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  
  445.  
  446.  

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