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

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