Télécharger la plaquette Cast3M


fr en

Télécharger pro2.eso

Retour à la liste

Numérotation des lignes : 

pro2

  1. C PRO2      SOURCE    MB234859  25/09/08    21:15:59     12358          
  2.       SUBROUTINE PRO2(IPGEA,IPCHEL,isort,IPOUT,ilphmo)
  3. C=====================================================================
  4. C
  5. C SOUS-PROGRAMME DE L'OPERATEUR PROIET
  6. C
  7. C ENTREES :
  8. C
  9. C   IPGEA  = POINTEUR DE UN OBJET MAILLAGE.
  10. C   IPCHEL = POINTEUR DE UN OBJET CHAMELEM
  11. C  isort  0  on veut un chpo brutal a 1 sous zone non repartitionné
  12. C             en sortie pointuer su un segment ipca
  13. C         1  on veut un champ  conforme  partitionne
  14. c ilphmo  pointeur sur un listmot phases du modele cible
  15. C SORTIES  :
  16. C
  17. C   IPOUT = POINTEUR SUR UN OBJET DE TYPE CHPOINT
  18. C
  19. C=====================================================================
  20. C
  21. C   TRES INSPIRE DE LA FONCTION ACCRO DEVELOPPEE PAR JMB
  22. C
  23. C=====================================================================
  24.  
  25.       IMPLICIT INTEGER(I-N)
  26.       IMPLICIT REAL*8 (A-H,O-Z)
  27. *-
  28. -INC CCREEL
  29.  
  30. -INC PPARAM
  31. -INC CCOPTIO
  32. -INC SMCOORD
  33. -INC SMELEME
  34. -INC SMCHPOI
  35. -INC SMLCHPO
  36. -INC SMLMOTS
  37. -INC SMCHAML
  38. -INC TMTRAV
  39. -INC CCGEOME
  40. -INC SMMODEL
  41. -INC SMLENTI
  42.       SEGMENT nomiii
  43.            character*(LOCOMP)  NOMINC(0)
  44.       endsegment
  45.       SEGMENT MDATA
  46.          INTEGER IPCA(NBSZ)
  47.       ENDSEGMENT
  48.       SEGMENT IPEX(NODES)
  49.       SEGMENT ISOM(NBSOM(IPT3.ITYPEL))
  50.       SEGMENT MTR1
  51.        CHARACTER*4 IPCOM(0)
  52.       ENDSEGMENT
  53.       SEGMENT SEPHAS
  54.        CHARACTER*8 LIPHAS(NPHAS)
  55.       ENDSEGMENT
  56.       SEGMENT LLIZA(NLIZA6)
  57.       SEGMENT/MTR4/(IPHAR(0))
  58.       segment inomil(0)
  59.       SEGMENT/MTRA/(NOPOIN(nbpts))
  60. *
  61.       LOGICAL drphas
  62. C
  63.       SEGMENT ICOUNT
  64.       INTEGER IEINT(2,NODES),IDEJVU(NODES),IACCRO
  65. * ient (1,i)= j, ient(2,i)=k  veut dire que le ieme point de ipt2 est dans
  66. * le k ieme element de la jeme sous zone.
  67. * idejvu(i)=1 veut dire que le ieme noeud de ipt2 est deja traité
  68.       REAL*8 QQQ(3,NODES),CKRIT(NODES)
  69. * CKRIT est le meilleur des criteres déja rencontré
  70.       ENDSEGMENT
  71.       SEGMENT ICOSOM
  72.         INTEGER IPTT(NBSZ),ICENV(NBSZ),IVOL(NBSZ)
  73.       ENDSEGMENT
  74.       SEGMENT WORK1
  75.         REAL*8 CORDEL(3,NBNN),QSI(3),XPO(3),SHPP(6,NBNN),V1(3),V2(3)
  76.         REAL*8 XEL(3,NBNN),BPSS(3,3),XPT(3),cored(3,NBNN),V3(3)
  77.         REAL*8 R(3),U(3),W(3,3)
  78.         INTEGER II(3)
  79.       ENDSEGMENT
  80.       segment siele
  81.         integer iiele(2,nsous)
  82.       endsegment
  83.       external shape
  84.  
  85. C    lecture optionnelle du modele associe au MCHAML origine
  86.       IPMOD=0
  87.       SIELE= 0
  88.       call LIROBJ('MMODEL  ',IPMOD,0,iret)
  89. *     write(6,*) ' dans pro2 lecture de ipmod ',ipmod
  90.       if(iret.eq.1) then
  91.          call ACTOBJ('MMODEL  ',IPMOD,1)
  92.          mmodel = ipmod
  93.          NSOUS=KMODEL(/1)
  94.          segini siele
  95.          DO 30 ISOUS = 1, NSOUS
  96.             IMODEL = KMODEL(ISOUS)
  97.             MELEME = IMAMOD
  98.             MELE = NEFMOD
  99. C
  100.             iiele(1,isous)=IMAMOD
  101.             iiele(2,isous)= INFELE(14)
  102.  30      CONTINUE
  103.       endif
  104. C   lecture optionnelle du critere de rattrapage de points exterieurs
  105.       ccrit=1.D-5
  106.       call lirree(ccrit,0,iretou)
  107.       ccrit = -ABS(ccrit)
  108. *       write(6,*) ' dans pro2 lecture de ccrit', ccrit
  109. C
  110. C     ON VA CHANGER LE MAILLAGE RECEPTEUR EN ELEMENTS POI1
  111. C         SI IL NE L EST PAS DEJA
  112.  
  113.       IPT2 = IPGEA
  114.       SEGACT IPT2
  115.       IF(IPT2.ITYPEL.NE.1) CALL CHANGE(IPGEA,1)
  116.       IF(IERR.NE.0) RETURN
  117. C
  118.       IPT2 = IPGEA
  119.       SEGACT IPT2
  120.       NODES = IPT2.NUM(/2)
  121.       SEGINI ICOUNT
  122.       IACCRO  = 0
  123.       do 3259 iu=1,nodes
  124.          ckrit(iu)=-1.d10
  125.  3259 continue
  126. C
  127. C   LE MAILLAGE SUPPORT DU CHAMELEM DOIT ETRE MASSIF
  128. C
  129.       MCHEL1 = IPCHEL
  130.       SEGACT MCHEL1
  131.       NBSZ = MCHEL1.IMACHE(/1)
  132.  
  133. *     Si MCHAML vide : erreur 
  134.       if (NBSZ.eq.0) then
  135.          moterr(1:8) = 'MCHAML  '
  136.          call erreur(1027)
  137.          return
  138.       endif
  139.  
  140. *     write(6,*) ' nbsz',nbsz
  141. C
  142.       nbnnma=0
  143.       DO 10 IOB=1, NBSZ
  144. C
  145.          IPT1=MCHEL1.IMACHE(IOB)
  146.          SEGACT IPT1
  147. C
  148.          ITY = IPT1.ITYPEL
  149. C   En dimension 1 : EF massifs s'appuient sur SEG2 ou SEG3
  150.          IF (IDIM.EQ.1) THEN
  151.             IF (ITY.NE.2.AND.ITY.NE.3) THEN
  152.                CALL ERREUR(16)
  153.                RETURN
  154.             ENDIF
  155.          ELSE
  156. C   En dimension 2,3 :
  157. cbp           IF (ITY.LT.4) THEN
  158. C            IF (ITY.LT.4.OR.(IDIM.eq.3.and.ITY.LT.14)) THEN
  159. C               write(ioimp,*) 'IDIM,ITYPEL=',IDIM,ITY
  160. C               CALL ERREUR(16)
  161. C               RETURN
  162. C            ENDIF
  163.          ENDIF
  164.          nbnnma=max(nbnnma,ipt1.num(/1))
  165. C
  166.  10   CONTINUE
  167. C      SEGACT MCOORD
  168.  
  169.       SEGINI MDATA
  170. Ckich on regarde les phases du MCHAML a projeter
  171.       nphas = nbsz
  172.       segini sephas
  173.       liphas(1) = mchel1.conche(1)(17:24)
  174.       nphas = 1
  175.       if ( mchel1.conche(/2).gt.1) then
  176.          do 11 jpha =2, mchel1.conche(/2)
  177.             do k = 1,nphas
  178.                if (mchel1.conche(jpha)(17:24).eq.liphas(k)) goto 11
  179.             enddo
  180.             nphas = nphas + 1
  181.             liphas(nphas) = mchel1.conche(jpha)(17:24)
  182.  11      continue
  183.          segadj sephas
  184.       endif
  185. Ckich on fabrique un maillage contenant tous les maillages du champ
  186. c de pointeurs distincts
  187.       nbsous=nbsz
  188.       nbelem=0
  189.       nbnn=0
  190.       nbref=0
  191.       segini ipt6
  192.  
  193.       nbsous = 1
  194.       ipt6.lisous(1) = mchel1.imache(1)
  195.       DO 20 IOB=2,NBSZ
  196.  
  197.          ipt5 = mchel1.imache(iob)
  198.          segact ipt5
  199.          if   (ipt5.itypel.eq.48) goto 20
  200.          do is= 1, nbsous
  201.             if (ipt6.lisous(is).eq.mchel1.imache(iob)) goto 20
  202.          enddo
  203. c verifie que les maillages sont disjoints : avertit
  204.          do 22 is= 1, nbsous
  205.             ip1 = ipt6.lisous(is)
  206.             ip2 = mchel1.imache(iob)
  207.             call interb(ip1,ip2,iret,iob1)
  208.             if (iret.eq.1) goto 22
  209.             write(ioimp,*) 'maillages non disjoints dans le MCHAML'
  210. * kich : on peut faire mieux, pardon.
  211.             goto 23
  212.  22      continue
  213.  
  214.  23      nbsous = nbsous + 1
  215.          IPT6.lisous(nbsous)=MCHEL1.IMACHE(IOB)
  216.  20   continue
  217.       segadj ipt6
  218. C
  219. C       ICI LE  ZONAGE
  220. C============================================
  221.       CALL  ZONEG2(IPT6,IPT2,ICount,1)
  222. c       write(6,*) ' dans pro2 sortie de zoneg2'
  223. C============================================
  224. C
  225. C-------------------------
  226. C            TABLEAU  des POINT PAS  VUS
  227. C     en previson de rectifications ulterieures
  228.       NBOUB  = NODES - IACCRO
  229. *       write(6,*) ' nboubl ' , nboub
  230.       SEGACT IPT2
  231.  
  232. *      SEGINI IOUBLI
  233. *     SEGACT IEXCLU
  234.       IUOO = 0
  235.       DO 21 I=1,NODES
  236. *         IP = IPT2.NUM(1,I)
  237.          IF(IDEJVU(I).EQ.1) THEN
  238.             go to 210
  239.          ELSE
  240.             IF(CKRIT(I).GT.CCRIT) THEN
  241.                IDEJVU(I)=1
  242.                IUOO= IUOO + 1
  243. *      IPTOU(IUO,1)=IP
  244.          ELSE
  245. *            IP = IPT2.NUM(1,I)
  246. *            write(6,*) ' ip ckrit (i)'
  247. *            write(6,*) ip, ckrit(i), qqq(1,i),qqq(2,i),qqq(3,i)
  248.                GO TO 21
  249.             ENDIF
  250.          endif
  251.  210     IPCA(IEINT(1,I))=1
  252.  21   CONTINUE
  253.       iaccro=iaccro+iuoo
  254. *       write(6,*) ' nodes iaccro',nodes,iaccro
  255. CC   on donne le nombre de points non projetés
  256.  
  257.       if (NBOUB.ne.iuoO) then
  258.          inonac=NBOUB-IUOO
  259.          INTERR(1)=inonac
  260. *       write(6,*) ' nombre de points non accrochés ',interr(1)
  261.          CALL ERREUR(-322)
  262.          CALL SOUCIS(inonac)
  263.       endif
  264.  
  265.       segact ipt6
  266. Ckich : on fabrique 1 chpoint par phase presente dans le MMODEL
  267. c stocke dans un listchpo, dans l ordre de MOTS (posibilite terme nul)
  268.       if (ilphmo.lt.0) then
  269.          jgn = 8
  270.          jgm = 1
  271.          segini mlmots
  272.          mots(1) = liphas(1)
  273.       else
  274.          mlmots = ilphmo
  275. c     n1 = nphas
  276.          n1 = mots(/2)
  277.          segini mlchpo
  278.       endif
  279.  
  280.       ipout = 0
  281.       DO 700 IMP = 1,MOTS(/2)
  282.          drphas = .false.
  283.          do iphas = 1,NPHAS
  284.            if(mots(imp).eq.liphas(iphas)) then
  285.              drphas = .true.
  286.              goto 703
  287.            endif
  288.          enddo
  289.  703     if(.not.drphas) then
  290.             write(6,*) ' pas de projection phase ', mots(imp)
  291.             goto 700
  292.          endif
  293.  
  294. C** recherche du noms de composantes existantes dans le chpoint
  295.          segini nomiii
  296.          NA=0
  297.          DO 220 i=1,nbsz
  298. ckich maintenant maillage disjoints        if(ipca(i).eq.0) go to 220
  299.             if (mchel1.conche(i)(17:24).ne.liphas(iphas)) goto 220
  300.             mchaml = mchel1.ichaml(i)
  301.             segact mchaml
  302.             do 221 j=1,nomche(/2)
  303.                IF( TYPCHE(J).NE.'REAL*8') GO TO 221
  304.                IF(NA.NE.0) THEN
  305.                   do 222 k=1,nominc(/2)
  306.                      if( nominc(k).eq.nomche(j)) go to 221
  307.  222              continue
  308.                   nominc(**)=nomche(j)
  309.                ELSE
  310.                   nominc(**)=nomche(j)
  311.                   na=1
  312.                ENDIF
  313.  221        continue
  314.  220     continue
  315.          nnin=nominc(/2)
  316.          nnnoe=iaccro
  317. *       write(6,*) ' nnin nnoe ',nnin,nnnoe,iphas
  318.          segini mtrav
  319. C
  320. C****  remplissage de mtrav en calculant les valeurs
  321. C
  322.          do 311 iu = 1 , nominc(/2)
  323.             inco(iu)=nominc(iu)
  324.  311     continue
  325.          segsup nomiii
  326.  
  327. c recense les zones du mchel1 compatibles avec le maillage de ipt6 et la phase
  328.          nliza6 = ipt6.lisous(/1)
  329.          segini lliza
  330.          do iza6 = 1,nliza6
  331.             do ips = 1,nbsz
  332.                IF (ipt6.lisous(iza6).eq.mchel1.imache(ips).AND.
  333.      &              liphas(iphas).eq.mchel1.conche(ips)(17:24)) THEN
  334.                   if (lliza(iza6).eq.0) then
  335.                      jg = 1
  336.                      segini mlenti
  337.                      LLIZA(IZA6) = mlenti
  338.                      lect(1) = ips
  339.                   else
  340.                      mlenti = lliza(iza6)
  341.                      jg = lect(/1)+1
  342.                      segadj mlenti
  343.                      lect(jg) = ips
  344.                   endif
  345.                ENDIF
  346.             enddo
  347.             mlenti = lliza(iza6)
  348.          enddo
  349.  
  350.          segact ipt2
  351.          ipla=0
  352.          nbnn = nbnnma
  353.          SEGINI ,WORK1
  354.          do 224 I=1,nodes
  355.             if(idejvu(i).eq.1) then
  356.                ipla=ipla+1
  357.                igeo(ipla)=ipt2.num(1,i)
  358.                ip=ipt2.num(1,i)
  359.                igeo(ipla)=ipt2.num(1,i)
  360.                iza6=ieint(1,i)
  361.                iel=ieint(2,i)
  362. *
  363.                IF (LLIZA(iza6).EQ.0)  goto 302
  364.                mlenti = lliza(iza6)
  365.                nza = lect(/1)
  366.  
  367.                do 300 iz = 1,nza
  368.                   IZA=lect(iz)
  369.  
  370.                   mchaml=mchel1.ichaml(iza)
  371.                   ipt1=mchel1.imache(iza)
  372.                   segact ipt1, mchaml
  373.  
  374.  
  375.                   if(ipmod.ne.0) then
  376. C    recherche  element fini  pour fonction de forme dans qsijs
  377.                      do lk=1,nsous
  378.                         if(ipt1.eq.iiele(1,lk)) then
  379.                            kel=iiele(2,lk)
  380.                            goto  31
  381.                         endif
  382.                      enddo
  383.  31                  continue
  384.                   else
  385. C  par defaut on utilise les fonctions de transformation geometriques
  386.                      kel =IPT1.ITYPEL
  387.                   endif
  388. C
  389. C         Coordonnées des sommets de l'élements IEL
  390. C
  391.                   nbnn=ipt1.num(/1)
  392.                   CALL DOXE(XCOOR,IDIM,NBNN,IPT1.NUM,IEL,CORDEL)
  393. C
  394.                   IDI2 = IDIM
  395.                   do ih=1,idim
  396.                      qsi(ih)=qqq(ih,i)
  397.                   enddo
  398. C
  399. C         Coordonnées  du point dans l'element de référence
  400. C
  401.                   CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,IRET)
  402. *         if ( qsi(1).lt.-1.00001.or.qsi(1).gt.1.00001) then
  403. *         if ( qsi(2).lt.-1.00001.or.qsi(2).gt.1.00001) then
  404. *         if ( qsi(3).lt.-1.00001.or.qsi(3).gt.1.00001) then
  405. *         write(6,*) ip,iza,iel,kel,(qsi(ih),ih=1,idim)
  406. *         write(6,*) ( shpp(1,ioup),ioup=1,nbnn)
  407. *         endif
  408. *        endif
  409. *        endif
  410. C
  411. C         Calcul des valeurs pour chaque type de composantes
  412. C
  413.                   do 200 il=1,ielval(/1)
  414.                      IF (TYPCHE(il).EQ.'REAL*8') THEN
  415. C
  416.                         do 225 kcol=1,nnin
  417. *        write(6,2323)il,nomche(il)(1:4),inco(1)
  418. * 2323   format('IL ', I2, 'nomche', A4, ' inco' , A4)
  419.                            if(nomche(il).eq.inco(kcol)) go to 226
  420.  225                    continue
  421.  226                    continue
  422. *             if(kcol.ne.1) then
  423. *                      write(6,*) 'kcol nnin',kcol,nnin
  424. *                      write(6,*)' iza il ipla',iza,il,ipla
  425. *          return
  426. *             endif
  427.                         MELVAL = IELVAL(il)
  428.                         SEGACT MELVAL
  429.                         IEMN = MIN(IEL,VELCHE(/2))
  430.                         NHAR(kcol) = MCHEL1.INFCHE(iza,3)
  431.                         DO 150, NP = 1, NBNN
  432.                            IBMN = MIN(NP,velche(/1))
  433.                            bb(kcol,ipla)=SHPP(1,NP)*VELCHE(IBMN,IEMN)
  434.      $                          +bb(kcol,ipla)
  435. C
  436.  150                    CONTINUE
  437.                         ibin(kcol,ipla)=ibin(kcol,ipla)+1
  438. C
  439.                      ENDIF
  440.  200              CONTINUE
  441.  300           continue
  442.  302           continue
  443.             endif
  444.  224     CONTINUE
  445.  
  446.          do iza6 = 1,nliza6
  447.             mlenti = lliza(iza6)
  448.             segsup mlenti
  449.          enddo
  450.          segsup lliza
  451.  
  452.          ITRAV=MTRAV
  453.          DO 4370 IB=1,NNNOE
  454.             DO 43701 IA=1,NNIN
  455. *      write(6,fmt='('' ia ib bb ibin '',2i4,e12.5,i5)')ia,ib,bb(ia,ib),
  456. *     * ibin(ia,ib)
  457.                IF (IBIN(IA,IB).GT.1)  THEN
  458. *        write(6,*) ' ibin gt 1  ib ,ia' , ib ,ia, ibin(ia,ib)
  459.                   BB(IA,IB)=BB(IA,IB) / IBIN(IA,IB)
  460.                   IBIN(IA,IB)=1
  461.                ENDIF
  462. 43701       CONTINUE
  463.  4370    CONTINUE
  464. C
  465. C     RECONSTUCTION DE LA PARTITION
  466. C
  467.          CALL CRECHP(ITRAV,ICHPOI1)
  468. C
  469.          IPOUT=ICHPOI1
  470.          MCHPOI=ICHPOI1
  471. *   le champ est forcement diffus
  472.          JATTRI(1)=1
  473.          MTYPOI='PROJECTI'
  474.          IFOPOI=IFOUR
  475.          SEGSUP MTRAV
  476. ckich dans le listchpoint
  477. *      write(6,*) 'pro2',ilphmo,mchpoi,ipout
  478.          if (ilphmo.gt.0) then
  479. *      write(6,*) 'pro2',mchpoi,ipout,'iphas',liphas(iphas),iphas,imp
  480.             ichpoi(imp) = MCHPOI
  481.          endif
  482.  700  CONTINUE
  483.       if (ipout.le.0) then
  484.          call erreur(-79)
  485.          return
  486.       endif
  487.       if (ilphmo.gt.0) ipout = mlchpo
  488.       SEGSUP ICOUNT,MDATA
  489.  
  490.       END
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  

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