Télécharger la plaquette Cast3M


fr en

Télécharger testma.eso

Retour à la liste

Numérotation des lignes : 

testma

  1. C TESTMA    SOURCE    PV        20/03/30    21:25:08     10567          
  2.       SUBROUTINE TESTMA(IPCHE,IMEL,COND,CONST,IRET,IMODI)
  3. *______________________________________________________________________
  4. *
  5. *             test sur les maillages
  6. *
  7. *    entrees :
  8. *    ---------
  9. *       ipche   chamelem dont il faut vérifier le maillage (type mchaml)
  10. *       imel    maillage du modèle: sert de référence (type meleme)
  11. *               Maillage actif en Entree (?) et en Sortie (oui)
  12. *       cond    vrai si il faut tenir compte de const
  13. *       const   nom de constituant
  14. *
  15. *
  16. *    sortie :
  17. *    --------
  18. *       iret    chamelem réordonné si nécessaire
  19. *               = 0 si pb
  20. *       imodi   =1 il a fallu decouper le maillage imel
  21. *
  22. *______________________________________________________________________
  23. *
  24. *                        declarations
  25. *
  26.       IMPLICIT REAL*8(A-H,O-Z)
  27.       IMPLICIT INTEGER(I-N)
  28. *
  29.  
  30. -INC PPARAM
  31. -INC CCOPTIO
  32. *
  33. -INC SMCHAML
  34. -INC SMELEME
  35. -INC SMCOORD
  36. -INC SMINTE
  37. *
  38.       SEGMENT ICPR(NPT)
  39.       SEGMENT ICPEL(ICPELD)
  40.       SEGMENT ITRAV(NNNT)
  41.       SEGMENT NIMEL(NNNT)
  42.       SEGMENT IPERM(NNNT)
  43.       SEGMENT NZONE(NSOUS)
  44.       SEGMENT ICOM(NBEL)
  45. *
  46.       SEGMENT  ITRA
  47.          INTEGER ICC (IA),IRE(IA,IMA)
  48.       ENDSEGMENT
  49. *
  50.       SEGMENT WTRAV
  51.         REAL*8 QSIREF(NBPGAU),QSICHAM(NBPGAU),ETAREF(NBPGAU)
  52.         REAL*8 ETACHAM(NBPGAU),DZEREF(NBPGAU),DZECHAM(NBPGAU)
  53.         REAL*8 XECHAM(3,NP1),XEREF(3,NP1)
  54.         INTEGER TABOK(NBPGAU),TAB(NBPGAU)
  55.       ENDSEGMENT
  56.       LOGICAL INIT,COND
  57.       CHARACTER*(*) CONST
  58. *
  59. *                        executable
  60. *
  61.       NPT = nbpts
  62.       ICPR  = 0
  63.       ICPEL = 0
  64. *
  65.       IRET  = 0
  66.       IMODI = 0
  67. *
  68.       IPT1 = IMEL
  69.       SEGACT IPT1
  70.       NBSOU1=IPT1.LISOUS(/1)
  71. *  La routine ne travaille qu'avec un maillage elementaire
  72.       IF ( NBSOU1 .GT. 1) THEN
  73. C        SEGDES,IPT1
  74.         call erreur(25)
  75.         RETURN
  76.       ENDIF
  77.       MELEME = IPT1
  78. *
  79. *     champ servant à stocker provisoirement le resultat
  80. *
  81.       MCHELM = IPCHE
  82.       SEGACT MCHELM
  83.       SEGINI,MCHEL1=MCHELM
  84. *
  85. *     segment servant a stocker l'ordre de stockage dans mchel1
  86.       NS01=0
  87.       NSOUS=ICHAML(/1)
  88.       SEGINI NZONE
  89. *
  90. *     creation de wtrav : on l'ajustera ensuite si necessaire
  91. *     mais on sort la creation destruction de la boucle
  92. *     la definition doit etre coherente avec celle dans tabgau
  93.       nbpgau=0
  94.       np1=0
  95.       wtrav=0
  96.       segini wtrav
  97. *
  98. *     boucle conditionnelle sur les sous zones du maillage imel
  99. *
  100.       IBOU1 = 0
  101. *
  102.       SEGINI ICPR
  103.   1   CONTINUE
  104.          IBOU1 = IBOU1 + 1
  105. *
  106. * Initialisations de ICPR, ICPEL la premiere fois puis
  107. * remises a zero pour les fois suivantes faite apres 2
  108. *
  109.          IF (NBSOU1.NE.0) THEN
  110.             MELEME = IPT1.LISOUS(IBOU1)
  111.          ENDIF
  112.          SEGACT MELEME
  113. *
  114. *        nnnt :nb d'elements de imel pour la sous zone consideree
  115. *
  116.          NNNT=NUM(/2)
  117. *
  118.          IA=0
  119.          IMA=0
  120. *
  121. *  fabrication de icpel qui dira combien de fois un point apparait
  122. *              de icpr qui donne une numerotation locale
  123. *
  124.          DO J=1,NNNT
  125.             DO 3 K=1,NUM(/1)
  126.                ID=NUM(K,J)
  127.                IF(ICPR(ID).NE.0) GO TO 3
  128. *              numerotation locale
  129.                IA=IA+1
  130.                ICPR(ID)=IA
  131.    3     CONTINUE
  132.          enddo
  133.          ICPELD=IA
  134.          IF(ICPEL.EQ.0) THEN 
  135.            SEGINI ICPEL
  136.          ELSE
  137.            IF(ICPELD.GT.ICPEL(/1)) SEGADJ ICPEL
  138.            call ooozmr(icpel(1),icpeld)
  139.          ENDIF
  140.  
  141.          DO  J=1,NNNT
  142.             DO  K=1,NUM(/1)
  143.                ID=NUM(K,J)
  144.                ICPEL(ICPR(ID))=ICPEL(ICPR(ID))+1
  145. *              on retient le max des occurences des points
  146.                IMA=MAX(ICPEL(ICPR(ID)),IMA)
  147.          enddo
  148.          enddo
  149. *
  150. *   fabrication de itra :
  151. *       icc    donne le nombre d'elements  a partir du numero local
  152. *       ire    donne les numeros d'elements de imel
  153. *              a partir du numero local et de l'occurence
  154. *
  155.          SEGINI ITRA
  156.          DO J=1,NNNT
  157.             DO K=1,NUM(/1)
  158.                ID=ICPR(NUM(K,J))
  159.                ICC(ID)=ICC(ID)+1
  160.                IA=ICC(ID)
  161.                IRE(ID,IA)=J
  162.          enddo
  163.          enddo
  164. *
  165. *        il faut maintenant regarder si dans un imache il existe
  166. *        les elements de la sous zone du meleme imel
  167. *        boucle sur les sous zone du chamelem
  168. *
  169.          IFLAG=0
  170. *  initialisons ipt2
  171.          IPT2=-1
  172.          DO 10 I=1,ICHAML(/1)
  173.             IF (COND) THEN
  174.                  IF (CONCHE(I) .NE. CONST) GOTO 10
  175.             ENDIF
  176.             MCHAML=ICHAML(I)
  177.             IPT2  =IMACHE(I)
  178.             IF (IPT2 .EQ. MELEME) THEN
  179.                   IFLAG=1
  180.                   NS01=NS01+1
  181.                   IF(NS01 .GT. ICHAML(/1))THEN
  182.                     N1=NS01
  183.                     N3=INFCHE(/2)
  184.                     L1=TITCHE(/1)
  185.                     SEGADJ,MCHEL1
  186.                     DO N33=1,N3
  187.                        MCHEL1.INFCHE(NS01,N33)=MCHEL1.INFCHE(I,N33)
  188.                     ENDDO
  189.                   ENDIF
  190.                   NSOUS=NSOUS+1
  191.                   SEGADJ,NZONE
  192.                   NZONE(NS01)=I
  193.                   MCHEL1.ICHAML(NS01)=MCHEL1.ICHAML(I)
  194.                   MCHEL1.IMACHE(NS01)=MELEME
  195.                   MCHEL1.CONCHE(NS01)=MCHEL1.CONCHE(I)
  196.                   GOTO 10
  197.             ENDIF
  198.             SEGACT,MELEME,IPT2
  199.             IF(ITYPEL.NE.IPT2.ITYPEL) GO TO 11
  200. *
  201.             NBEL = NUM(/2)
  202.             SEGINI ICOM
  203. *
  204.             INIT=.FALSE.
  205. *
  206.             ICO=0
  207.             NP1=IPT2.NUM(/1)
  208.             DO 12 K=1,IPT2.NUM(/2)
  209.                DO 13 L=1,NP1
  210.                   ID=IPT2.NUM(L,K)
  211.                   IDD = ICPR(ID)
  212.                   IF(IDD.EQ.0) GO TO 12
  213.    13          CONTINUE
  214. *
  215. *              ok l'element k possede ses noeuds dans imel
  216. *
  217.                IF(ITYPEL.EQ.1) THEN
  218.                   ID=IPT2.NUM(1,K)
  219.                   IDD=ICPR(ID)
  220.                   IRE1=IRE(IDD,1)
  221.                   GOTO 20
  222.                ENDIF
  223. *
  224. *              ces noeuds correspondent-t-ils a un meme element ire1
  225. *
  226.                ID = IPT2.NUM(1,K)
  227.                IDD=ICPR(ID)
  228.                IDE=ICC(IDD)
  229. c          boucle sur les occurences de idd dans les elements
  230.                DO 14 L=1,IDE
  231.                   IRE1=IRE(IDD,L)
  232.                   DO 15 M=2,NP1
  233.                      IDD2=ICPR(IPT2.NUM(M,K))
  234.                      IF(IDD2.EQ.0) GO TO 12
  235.                      IDE2=ICC(IDD2)
  236.                      DO 16 N=1,IDE2
  237.                         IF(IRE(IDD2,N).EQ.IRE1) GO TO 15
  238.    16                CONTINUE
  239.                      GO TO 14
  240.    15             CONTINUE
  241.                   GOTO 20
  242.    14          CONTINUE
  243. *
  244. *        l'element n'appartient pas a meleme
  245.                GOTO 12
  246. *
  247. *              itrav donne le numero de l'element dans le chamelem
  248. *              nimel donne le numero de l'element dans la sous zone
  249. *                    de imel
  250. cbp            IPERM = 0 si pas de permutation des noeuds
  251. *
  252.    20          CONTINUE
  253. c
  254.                IF(.NOT.INIT) THEN
  255.                   INIT=.TRUE.
  256.                   SEGINI ITRAV
  257.                   SEGINI NIMEL
  258.                   SEGINI IPERM
  259.                ENDIF
  260.                ICO=ICO+1
  261.                IF (ICO .GT. NNNT) THEN
  262.                  CALL ERREUR(472)
  263.                  RETURN
  264.                ENDIF
  265.                ITRAV(ICO)=K
  266.                NIMEL(ICO)=IRE1
  267.                ICOM(IRE1)=1
  268.                IPERM(ICO)=0
  269.                do iii=1,NP1
  270.                  if(IPT2.NUM(iii,K).ne.NUM(iii,IRE1)) then
  271.                    IPERM(ICO)=iii
  272.                    goto 12
  273.                  endif
  274.                enddo
  275.    12       CONTINUE
  276. *
  277.             IF(INIT) THEN
  278. *
  279. *              le nb d'elements est-il egal pour la sous zone imel et
  280. *              la sous zone mchaml
  281. *
  282.                IF (ICO.EQ.NNNT) THEN
  283.                   IFLAG=1
  284.                   NS01=NS01+1
  285.                   IF(NS01 .GT. ICHAML(/1))THEN
  286.                     N1=NS01
  287.                     N3=INFCHE(/2)
  288.                     L1=TITCHE(/1)
  289.                     SEGADJ,MCHEL1
  290.                     DO N33=1,N3
  291.                        MCHEL1.INFCHE(NS01,N33)=MCHEL1.INFCHE(I,N33)
  292.                     ENDDO
  293.                   ENDIF
  294.                   NSOUS=NSOUS+1
  295.                   SEGADJ,NZONE
  296.                   NZONE(NS01)=I
  297.                   MINTE=INFCHE(I,4)
  298.                   IPMINT=MINTE
  299.                   IF (IPMINT.NE.0) SEGACT,MINTE
  300.                   MCHAML=ICHAML(I)
  301.                   SEGINI,MCHAM1=MCHAML
  302.                   MCHEL1.ICHAML(NS01)=MCHAM1
  303.                   MCHEL1.IMACHE(NS01)=MELEME
  304.                   MCHEL1.CONCHE(NS01)=MCHEL1.CONCHE(I)
  305. *
  306.                   DO 24 ICOMP=1,MCHAM1.IELVAL(/1)
  307.                      MELVAL=MCHAM1.IELVAL(ICOMP)
  308.                      SEGACT MELVAL
  309.                      IMEM=0
  310.                      IF (MCHAM1.TYPCHE(ICOMP).EQ.'REAL*8') THEN
  311.                         N1PTEL=VELCHE(/1)
  312.                         N1EL=ICO
  313.                         IF (VELCHE(/2).EQ.1) THEN
  314.                            N1EL=1
  315.                            IMEM=1
  316.                         ENDIF
  317.                         NEL=N1EL
  318.                         N2PTEL=0
  319.                         N2EL  =0
  320.                      ELSE
  321.                         N2PTEL=IELCHE(/1)
  322.                         N2EL=ICO
  323.                         IF (IELCHE(/2).EQ.1) THEN
  324.                            N2EL=1
  325.                            IMEM=1
  326.                         ENDIF
  327.                         NEL=N2EL
  328.                         N1PTEL=0
  329.                         N1EL  =0
  330.                      ENDIF
  331. *
  332.                      SEGINI MELVA1
  333. *
  334.                      DO 21 J=1,NEL
  335.                         K   =ITRAV(J)
  336.                         IRE1=NIMEL(J)
  337.                         KTEST=K
  338.                         ITEST=IRE1
  339.                         IF (NEL.EQ.1) THEN
  340.                            IRE1=1
  341.                            IF (IMEM.EQ.1) THEN
  342.                               K=1
  343.                            ENDIF
  344.                         ENDIF
  345. *
  346. *    dans le cas des champs qui ne sont pas définis aux noeuds
  347. *    vérification de la bonne orientation des maillages
  348. *    sinon création d'un tableau de correspondance
  349. *
  350.                         NBPGAU=MAX(N1PTEL,N2PTEL)
  351. cbp   -cas ou l'on doit identifier la position des pts d'integration
  352. cbp    (ne marche pas avec des coques epaisses MFR=5,
  353. cbp     ni avec DKT dont nombre de points >1 dans l'épaisseur)
  354.                         IF((IPMINT.NE.0).and.(IPERM(J).ne.0)) THEN
  355.  
  356.                           if (qsiref(/1).ne.nbpgau .or.
  357.      &                        xeref(/2).ne.np1) then
  358.                             segadj,wtrav
  359.                           endif
  360.                           CALL TABGAU(IPMINT,ITEST,KTEST,NP1,MELEME,
  361.      +                           IPT2,NBPGAU,IRET1,wtrav)
  362.                           IF (IRET1.EQ.0) THEN
  363.                             CALL ERREUR(757)
  364.                             SEGSUP MCHEL1,ITRA
  365. C                            if (imel.ne.ipt2) SEGDES IPT2
  366. C                            if (imel.ne.meleme) segdes meleme
  367.                             GOTO 900
  368.                           ENDIF
  369.                           IF (N1PTEL.EQ.0) THEN
  370.                            DO 22 IGAU=1,N2PTEL
  371.                               IGAU1=TAB(IGAU)
  372.                               MELVA1.IELCHE(IGAU,IRE1)=IELCHE(IGAU1,K)
  373.    22                      CONTINUE
  374.                           ELSE
  375.                            DO 23 IGAU=1,N1PTEL
  376.                               IGAU1=TAB(IGAU)
  377.                               MELVA1.VELCHE(IGAU,IRE1)=VELCHE(IGAU1,K)
  378.    23                      CONTINUE
  379.                           ENDIF
  380.  
  381. cbp   -cas ou on recopie directement sans chercher la permutation
  382.                         ELSE
  383.                          IF (N1PTEL.EQ.0) THEN
  384.                            DO 25 IGAU=1,N2PTEL
  385.                               MELVA1.IELCHE(IGAU,IRE1)=IELCHE(IGAU,K)
  386.    25                      CONTINUE
  387.                          ELSE
  388.                            DO 26 IGAU=1,N1PTEL
  389.                               MELVA1.VELCHE(IGAU,IRE1)=VELCHE(IGAU,K)
  390.    26                      CONTINUE
  391.                          ENDIF
  392.  
  393.                         ENDIF
  394. *
  395.    21                CONTINUE
  396.                      MCHAM1.IELVAL(ICOMP)=MELVA1
  397. C                     SEGDES MELVA1,MELVAL
  398.    24             CONTINUE
  399. C                  SEGDES MCHAM1,MCHAML
  400.                   SEGSUP ITRAV,NIMEL,IPERM
  401.                ELSE
  402. *
  403. *                 le maillage est a cheval sur plusieurs sous zones
  404. *                 du mchaml
  405. *                 il va falloir scinder le maillage en deux
  406. *                 et recommencer le processus
  407. *
  408.                   IMODI = 1
  409. *       creation des deux maillages dont l'union est le maillage meleme
  410.                   NBNN = NUM(/1)
  411.                   NBELEM = ICO
  412.                   NBREF = 0
  413.                   NBSOUS = 0
  414.                   SEGINI IPT3
  415.                   NBELEM = NUM(/2) - ICO
  416.                   IF (NBELEM.LE.0) THEN
  417.                     CALL ERREUR(26)
  418.                     GOTO 900
  419.                   ENDIF
  420.                   SEGINI IPT4
  421.                   IPT3.ITYPEL = ITYPEL
  422.                   IPT4.ITYPEL = ITYPEL
  423.                   I3 = 0
  424.                   I4 = 0
  425.                   DO 45 II=1,NUM(/2)
  426.                     IF (ICOM(II) .EQ. 1) THEN
  427. *                    l'element est commun aux deux maillages
  428.                      I3=I3+1
  429.                      IPT3.ICOLOR(I3)=ICOLOR(II)
  430.                      DO 43 JJ=1,NUM(/1)
  431.                         IPT3.NUM(JJ,I3)=NUM(JJ,II)
  432.    43                CONTINUE
  433.                     ELSE
  434. *                    l'element appartient seulement a meleme
  435.                      I4=I4+1
  436.                      IPT4.ICOLOR(I4)=ICOLOR(II)
  437.                      DO 44 JJ=1,NUM(/1)
  438.                         IPT4.NUM(JJ,I4)=NUM(JJ,II)
  439.    44                CONTINUE
  440.                     ENDIF
  441.    45            CONTINUE
  442. *
  443. *        modification de ipt1
  444. *
  445.                  IF (IPT1.LISOUS(/1) .LE. 1) THEN
  446.                      NBSOUS=2
  447.                      NBREF =0
  448.                      NBELEM=0
  449.                      NBNN  =0
  450.                      SEGINI IPT5
  451.                      IPT5.LISOUS(1)=IPT3
  452.                      IPT5.LISOUS(2)=IPT4
  453. C                     if (imel.ne.ipt1) SEGDES IPT1
  454.                      IPT1=IPT5
  455.                  ELSE
  456.                      NBSOUS=IPT1.LISOUS(/1)+1
  457.                      NBREF =0
  458.                      NBELEM=0
  459.                      NBNN  =0
  460.                      SEGADJ,IPT1
  461.                      IPT1.LISOUS(IBOU1) =IPT3
  462.                      IPT1.LISOUS(NBSOUS)=IPT4
  463. C                     if (imel.ne.meleme) SEGDES MELEME
  464.                  ENDIF
  465.                  IBOU1 =IBOU1-1
  466.                  NBSOU1=NBSOUS
  467.                  SEGSUP ITRAV,NIMEL,IPERM
  468.                  GOTO 2
  469.                ENDIF
  470.             ENDIF
  471.             SEGSUP ICOM
  472.    11       CONTINUE
  473. C            if (imel.ne.ipt2) SEGDES IPT2
  474.    10    CONTINUE
  475. *
  476. *        on n'a pas trouve de sous zone dans le mchaml qui
  477. *        correspondent a une sous zone de imel
  478. *
  479.          IF (IFLAG.EQ.0) THEN
  480. C            if (imel.ne.meleme) segdes meleme
  481.             SEGSUP MCHEL1
  482.             IF(IPT2.LE.0) THEN
  483.                SEGSUP,ITRA
  484.             ELSE
  485. C              if (imel.ne.ipt2) SEGDES IPT2
  486.               MOTERR(1:8)='MAILLAGE'
  487.               MOTERR(9:16)='MCHAML'
  488.               CALL ERREUR(135)
  489.             ENDIF
  490.             GOTO 900
  491.          ELSE
  492.             SEGSUP,ITRA
  493. C            if (imel.ne.meleme) segdes meleme
  494. C            if (imel.ne.ipt2) SEGDES IPT2
  495.          ENDIF
  496.  
  497.    2     continue
  498.          if (ibou1.lt.nbsou1) then
  499. * remise a zero de icpr
  500.           DO J=1,NNNT
  501.              DO  K=1,NUM(/1)
  502.                 ID=NUM(K,J)
  503.                 ICPR(ID)=0
  504.              enddo
  505.           enddo
  506.           goto 1
  507.          endif
  508. *
  509. *     fin de la boucle sur les sous zones du maillage
  510. *
  511.       IF(NS01.NE.NSOUS) THEN
  512.          N1=NS01
  513.          N3=MCHEL1.INFCHE(/2)
  514.          L1=MCHEL1.TITCHE(/1)
  515.          SEGINI MCHEL2
  516.          MCHEL2.TITCHE=TITCHE
  517.          MCHEL2.IFOCHE=IFOCHE
  518.          DO 31 ISOUS=1,NS01
  519.             MCHEL2.ICHAML(ISOUS) = MCHEL1.ICHAML(ISOUS)
  520.             MCHEL2.IMACHE(ISOUS) = MCHEL1.IMACHE(ISOUS)
  521.             MCHEL2.CONCHE(ISOUS) = MCHEL1.CONCHE(ISOUS)
  522.             NSS=NZONE(ISOUS)
  523.             DO 33 N33=1,N3
  524.                MCHEL2.INFCHE(ISOUS,N33) = MCHEL1.INFCHE(NSS,N33)
  525.   33        CONTINUE
  526.   31     CONTINUE
  527.          SEGSUP MCHEL1
  528.          MCHEL1=MCHEL2
  529.       ENDIF
  530. C     SEGDES,MCHEL1
  531.       IRET=MCHEL1
  532.  
  533.  900  CONTINUE
  534. C      SEGDES MCHELM
  535.       SEGSUP,ICPR,ICPEL,NZONE,WTRAV
  536.  
  537.       END
  538.  
  539.  
  540.  
  541.  
  542.  
  543.  

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