Télécharger la plaquette Cast3M


fr en

Télécharger pron.eso

Retour à la liste

Numérotation des lignes : 

pron

  1. C PRON      SOURCE    OF166741  24/10/07    21:15:42     12016          
  2.  
  3.       SUBROUTINE PRON(IPMODE,IPCHE,IPCHT,ISUP,IPOUT)
  4.  
  5.       IMPLICIT INTEGER(I-N)
  6.       IMPLICIT REAL*8 (A-H,O-Z)
  7.  
  8. *-------------------------------------------------------------------- *
  9. *                                                                     *
  10. *   Sous-programme appelé par  l'opérateur   PROI                     *
  11. *                                                                     *
  12. *   Projection d'un chamelem aux noeuds ou aux points d'integration   *
  13. *   des elements support du modele                                    *
  14. *   Evaluation  faite au point courant par utilisation des fonctions  *
  15. *   de formes du modele associe au champ d origine (lu en option      *
  16. *   dans PRO2) ou par defaut a partir des fonctions geometriques      *
  17. *   associées à l'élément d origine                                   *
  18. *                                                                     *
  19. *   __________________________________________________________________*
  20. *                                                                     *
  21. *   ENTREE :                                                          *
  22. *          IPMODE  modele sur lequel on va projeter                   *
  23. *          IPCHE   MCHAML de caracteristiques (si le modele comprend  *
  24. *                  des coques)                                        *
  25. *          IPCHT  MCHAML a projeter (obligatoirement defini aux noeuds)
  26. *          ISUP   support du resultat  1 noeuds                       *
  27. *                                      2 pts d integration contraintes*
  28. *   SORTIE :                                                          *
  29. *          IPOUT   MCHAML resultat                                    *
  30. *                                                                     *
  31. *   Projection d'un chamelem sur le support d un modele               *
  32. *                                                                     *
  33. *---------------------------------------------------------------------*
  34. *
  35. -INC PPARAM
  36. -INC CCOPTIO
  37. -INC CCREEL
  38. *
  39. -INC SMCOORD
  40. -INC SMMODEL
  41. -INC SMCHAML
  42. -INC SMELEME
  43. -INC SMCHPOI
  44. -INC SMLCHPO
  45. -INC SMLMOTS
  46. -INC TMTRAV
  47. -INC SMINTE
  48.         character*(LOCOMP) nomm
  49.         segment snomip
  50.           character*(LOCOMP) nomip(0)
  51.         endsegment
  52.         segment ipcb(nsschp)
  53.         segment  icpr(2,nbpts)
  54. C
  55.       segment info
  56.         integer infell(JG)
  57.       endsegment
  58. C
  59.         segment sicoq
  60.         integer icocoq(nbsous),imfr(nbsous),iiele(nbsous)
  61.         integer iminte(nbsous),imint1(nbsous),ibnap(nbsous)
  62.         integer inbg(nbsous),ibbas(nbsous)
  63.         endsegment
  64. C
  65.       SEGMENT WRK4
  66.        REAL*8 BPSS(3,3),XEL(3,NBN1) ,XU(3),SHP(6,NBN1),S(6)
  67.        REAL*8 TXR(3,3,NBN1),XE(3,NBG),XEG(3,NBG3)
  68.       ENDSEGMENT
  69.  
  70.       SEGMENT VECT
  71.         REAL*8  VEC1(IDIM)
  72.         REAL*8  VEC2(IDIM)
  73.         REAL*8  VECN(IDIM,NBN1*2)
  74.       ENDSEGMENT
  75.  
  76.       segment phachm
  77.         integer iphach(nsouz)
  78.         integer mphach(nsouz)
  79.       endsegment
  80.       LOGICAL logach
  81.  
  82. c     Compteur du nombre de sous modeles de sure
  83.       ISURE = 0
  84.  
  85.       IPOUT = 0
  86.  
  87.       MMODEL=IPMODE
  88.       NSOUS=KMODEL(/1)
  89. *
  90.       if(IPCHE.NE.0) THEN
  91.         MCHELM=IPCHE
  92.         NCHAM=ICHAML(/1)
  93.       endif
  94. *
  95.       segact mcoord*mod
  96.       NBNOE = nbpts       
  97. *
  98. *     Création du maillage provisoire support pour interpolation
  99. *
  100. C ipt2  contiendra les super couche au point de gauss
  101.       NBSOUS = NSOUS
  102.       NBREF = 0
  103.       NBELEM = 0
  104.       NBNN =  0
  105.       SEGINI IPT2
  106.       SEGINI SICOQ
  107. c listmots des phases
  108.       ilphmo = -1
  109.       jgn = 8
  110.       jgm = nsous
  111.       segini mlmots
  112.       ilphmo = mlmots
  113.       jgm = 1
  114. *
  115. * Boucle sur l'ensemble des sous zones du modeles sur lequel on va
  116. *   projetter   pour   initialisations
  117.       nponou=0
  118.       npcoq =0
  119. *
  120.       DO 30 ISOUS = 1, NSOUS
  121. *
  122.         IMODEL = KMODEL(ISOUS)
  123. *
  124.         MELE = NEFMOD
  125.         MELEME = IMAMOD
  126.         NBELE1= NUM(/2)
  127.         NBN1 = NUM(/1)
  128.       IF(INFMOD(/1).NE.0)THEN
  129.       IF(INFMOD(1).NE.0)THEN
  130.           NPINT=INFMOD(1)
  131.           IBNAP(isous) =NPINT
  132.       ELSE
  133.           NPINT=0
  134.           ibnap(isous)=1
  135.       ENDIF
  136.       ELSE
  137.           NPINT=0
  138.           ibnap(isous)=1
  139.       ENDIF
  140. C
  141.       if(INFMOD(/1).lt.2+isup) then
  142.       call elquoi(mele,npint,isup,ipinf,imodel)
  143.       if (ierr.ne.0) return
  144.         info=ipinf
  145.         mfr  =infell(13)
  146.         iele =infell(14)
  147.         minte=infell(11)
  148.         minte1=infELL(12)
  149.         nbg=nbn1
  150.         if(isup.eq.2) then
  151.         nbg=shptot(/3)
  152.         elseif(isup.eq.3) then
  153.         nbg=infell(6)
  154.         elseif(isup.eq.4) then
  155.         nbg=infell(3)
  156.         elseif(isup.eq.5) then
  157.         nbg=infell(4)
  158.          endif
  159.         segsup info
  160.         else
  161.         mfr  =infele(13)
  162.         iele =infele(14)
  163.         minte=infmod(2+isup)
  164.         minte1=infmod(8)
  165.         nbg=nbn1
  166.         if(isup.eq.2) then
  167.         nbg=shptot(/3)
  168.         elseif(isup.eq.3) then
  169.         nbg=infele(6)
  170.         elseif(isup.eq.4) then
  171.         nbg=infele(3)
  172.         elseif(isup.eq.5) then
  173.         nbg=infele(4)
  174.         endif
  175.         endif
  176. C
  177.         imfr(isous)=mfr
  178.         iiele(isous)=iele
  179.         iminte(isous)=minte
  180.         imint1(isous)=minte1
  181.         inbg(isous)=nbg
  182.         NBSOUS = 0
  183.         NBREF = 0
  184.         NBELEM = NBELE1
  185.       IF (mfr.EQ.3.OR.mfr.EQ.5.OR.mfr.EQ.9) THEN
  186. C
  187.       if(IPCHE.eq.0) then
  188.         call erreur(404)
  189.         return
  190.       endif
  191. C
  192. C----------------------------coque ---------------------
  193.         icocoq(isous) = 1
  194.         npcoq = npcoq + 1
  195. C  nombre de points par nappes
  196.         if(isup.GT.2) then
  197.                 if(mfr.eq.5) then
  198. C -----------------      cas des  coques epaisses
  199.                         if(nbg.eq.8.or.nbg.eq.6) then
  200.                         IBNAP(isous)=2
  201.                         ibbas(isous)=nbg/2
  202.                         else
  203.                         IBNAP(isous)=nbg/nbn1*2
  204.                         endif
  205.                 else
  206.                 ibbas(isous) = nbg/ibnap(isous)
  207.                 endif
  208.         elseif(isup.eq.1) then
  209.         ibnap(isous)=1
  210.         ibbas(isous) = nbn1
  211.         elseif(isup.eq.2) then
  212.         ibnap(isous)= 1
  213.         ibbas(isous) = 1
  214.         endif
  215. *
  216. *      write(6,*) 'nombre de noeuds par nappe ',
  217. *     &       ibbas(isous),' nappes ',ibnap(isous)
  218. *
  219. * création du maillage  des points supports pour cette sz
  220. *
  221.         NBNN = ibbas(isous)*ibnap(isous)
  222.         nponou = nponou + nbnn*nbele1
  223.         SEGINI IPT1
  224.         IPT1.ITYPEL = 28
  225.         ipt2.lisous(isous) = ipt1
  226. *write(6,*) 'sz ',isous ,'nbnn nbelem  ipt1' ,nbnn,nbelem,ipt1
  227.       ELSE
  228. C-----------------------------massif ------------------------
  229.         if(isup.GT.2) then
  230.         NBNN = NBG
  231.         elseif(isup.eq.2) then
  232.         NBNN=1
  233.         else
  234.         NBNN= nbn1
  235.         endif
  236. C
  237.         SEGINI IPT1
  238.         if(isup.eq.2) then
  239.         else
  240.         IPT1.ITYPEL = 28
  241.         endif
  242.         ipt2.lisous(isous) = ipt1
  243.         nponou=nponou+nbnn*nbele1
  244. C
  245.       ENDIF
  246.  
  247.         if (isous.eq.1) then
  248.          mots(1) = conmod(17:24)
  249.         else
  250.          do ipl = 1,jgm
  251.           if (mots(ipl).eq.conmod(17:24)) goto 27
  252.          enddo
  253.          jgm = jgm + 1
  254.          mots(jgm) = conmod(17:24)
  255.  27      continue
  256.         endif
  257. C
  258. *        segsup info
  259. 30    CONTINUE
  260.  
  261.       segadj mlmots
  262. C-----------------------------------------------------------
  263. C
  264. C     si il y a des coques ou si l'on veut le champ aux points
  265. C  d integration, on va fabriquer des points support provisoires
  266. C
  267. C----------------------------------------------------------
  268. C         write(6,*) ' nombre de points crees ' ,nponou
  269. C (fdp) ajustement de MCOORD pour les noeuds provisoires
  270. C       il sera remis a sa taille initiale a la fin du programme
  271.          NBPTS = NBNOE + nponou
  272.          SEGADJ MCOORD
  273. *-----------------------------------------------------------
  274. *     Boucle sur l'ensemble des sous zones du modeles
  275. *------------------------------------------------------------
  276.         inupo =NBNOE
  277.       DO 100 ISOUS = 1,NSOUS
  278.         IMODEL = KMODEL(ISOUS)
  279.         MELEME = IMAMOD
  280.         IF (itypel.eq.48) then
  281.             ISURE = ISURE+1
  282.             goto 100
  283.         ENDIF
  284.         nbele1= num(/2)
  285.         nbn1 = num(/1)
  286.         ipt1=ipt2.lisous(isous)
  287.         segact ipt1*mod
  288.         NNB = ibbas(isous)
  289.         NBG =inbg(isous)
  290.         NBG3=NBG*3
  291.         SEGINI WRK4
  292.         SEGINI VECT
  293.       IF(ICOCOQ(ISOUS).NE.0) THEN
  294. C-----------------------------------------------------------
  295. C  -------------------------------LA SZ EST UNE COQUE ------
  296. C-----------------------------------------------------------
  297.            NUCHA = 0
  298.        DO 15, NUCH = 1, NCHAM
  299. *
  300.             IF ( CONCHE(NUCH).EQ.CONMOD.AND.
  301.      &           IMACHE(NUCH).EQ.IMAMOD) NUCHA = NUCH
  302. *
  303.   15      CONTINUE
  304. *
  305.           IF (NUCHA.NE.0) THEN
  306.             MCHAML=ICHAML(NUCHA)
  307. *
  308.             MELVA1 = 0
  309.             MELVA2 = 0
  310.             NCOMP = IELVAL(/1)
  311.             DO 20, I = 1, NCOMP
  312.               IF (NOMCHE(I).EQ.'EPAI') THEN
  313.                 MELVA1 = IELVAL(I)
  314.               ELSEIF (NOMCHE(I).EQ.'EXCE') THEN
  315.                 MELVA2 = IELVAL(I)
  316.               ENDIF
  317.   20        CONTINUE
  318.           ENDIF
  319. C
  320.       MELE = NEFMOD
  321. c   coque integree ou pas ?
  322.         npint=infmod(1)
  323.         nbnap=ibnap(isous)
  324.         mfr  =imfr(isous)
  325.         iele =iiele(isous)
  326.         minte=iminte(isous)
  327.         minte1=imint1(isous)
  328.         nbg = inbg(isous)
  329.         nnb = ibbas(isous)
  330. *-------------------------------------------------------------------
  331. *  Boucle sur les elements de la sous zone du modele  (c est une coque)
  332. *-------------------------------------------------------------------
  333.         DO 95 iel=1, NBELE1
  334.         ipo=0
  335.         call doxe(xcoor,idim,nbn1,num,iel,xel)
  336. *    write(6,*) ' coordonnees des noeuds '
  337. *    write(6,2000) ((xel(ik,jk),ik=1,idim),jk=1,nbn1)
  338. *    write(6,*) '-------------------- '
  339. 2000     format(3E12.5)
  340. * on veut les coordonnees aux point d integration et normales
  341. *
  342.       IF(ISUP.GT.1) THEN
  343. C  --- projection des points d'integration sur la surface moyenne--
  344.       DO 601 ip=1,nnb
  345.       do 600 ik=1,idim
  346.       xe(ik,ip)=0.D0
  347.       DO 602 jk = 1,NBN1
  348.       XE(Ik,ip) = XE(Ik,ip) + SHPTOT(1,JK,ip)*(XEL(IK,jk))
  349. 602   CONTINUE
  350. 600   CONTINUE
  351. 601   CONTINUE
  352.       do i=1,nnb
  353.       do j=1,idim
  354.       xeg(j,i)=xe(j,i)
  355.       enddo
  356.       enddo
  357. C      ELSEIF(isup.eq.2) THEN
  358. C----------------------------          on travaille aux  CDG
  359. C      do  ik=1,idim
  360. C      xeg(ik,1)=0.D0
  361. C      do  ip=1,nbg
  362. C      DO  jk = 1,NBN1
  363. C      XEG(Ik,1) = XEG(Ik,jk) + SHPTOT(1,JK,ip)*(XEL(IK,jk))
  364. C      enddo
  365. C      enddo
  366. C      enddo
  367.       ELSE
  368. C  -------   on travaille aux noeuds
  369.       do i=1,nbn1
  370.       do j=1,idim
  371.       xeg(j,i)=xel(j,i)
  372.       enddo
  373.       enddo
  374.       ENDIF
  375. C      write(6,*) '  coordonnes des points de  base '
  376. C     write(6,2576) (ip,(xeg(il,ip),il=1,3),ip=1,nnb)
  377. C      write(6,*) ' ----------------- '
  378. 2576  format(I4,2X,3E12.5)
  379.  
  380. C----------- les normales aux points ad hoc ------------------
  381. C
  382. *---------------    coque minces   3 noeuds ---------------
  383.        if(iele.eq.4)  then
  384. C    --------coq3   dkt  dst   ----------------
  385.       do 98 i=1,idim
  386.       VEC1(i)=xel(i,2)-xel(i,1)
  387. 98    VEC2(i)=xel(i,3)-xel(i,1)
  388.       vecn(1,1) = vec1(2)*vec2(3)-vec1(3)*vec2(2)
  389.       vecn(2,1) = vec1(3)*vec2(1)-vec1(1)*vec2(3)
  390.       vecn(3,1) = vec1(1)*vec2(2)-vec1(2)*vec2(1)
  391.       vnor =sqrt(vecn(1,1)*vecn(1,1)+vecn(2,1)*vecn(2,1)
  392.      & +vecn(3,1)*vecn(3,1))
  393.       do 91 i=1,idim
  394. 91    vecn(i,1)=vecn(i,1)/vnor
  395.       do 99 j=2,nbn1
  396.       do 99 i=1,idim
  397. 99    vecn(i,j)=vecn(i,1)
  398. C  petite rectif necessaire pour le tri3
  399.       if(isup.eq.2) then
  400.        do j=1,idim
  401.        xeg(j,1) = XZERO
  402.                do ip=1,3
  403.                xeg(j,1) = xeg(j,1) + xel(j,ip)
  404.                enddo
  405.        xeg(j,1)=xeg(j,1)/3.D0
  406.        enddo
  407.        endif
  408.        elseif (iele.eq.8) then
  409. C    --------------- coq4 ----------------
  410.       call cq4loc(xel,xe,BPSS,irrt,1)
  411. C  ici c est une estimation du plan moyen
  412. *     write(6,fmt='(3E12.5)') ((bpss(i,j),i=1,3),j=1,3)
  413.       do ip=1,ibbas(isous)
  414.       do i=1,3
  415.       vecn(i,ip)=bpss(3,i)
  416.       enddo
  417.       enddo
  418.         elseif (iele.eq.2) then
  419. C   ---------------- coq2 --------------
  420.         vnor=0.D0
  421.         do 92 i=1,idim
  422.         VEC1(i)=xel(i,2)-xel(i,1)
  423. 92      vnor=vnor + vec1(i)*vec1(i)
  424.         vnor = sqrt(vnor)
  425.         vecn(1,1)=-vec1(2)/vnor
  426.         vecn(2,1)=vec1(1)/vnor
  427.         vecn(1,2)=vecn(1,1)
  428.         vecn(2,2)=vecn(2,1)
  429. C
  430.         elseif (iele.eq.5.or.iele.eq.10) then
  431. *-----------    coques epaisses  coq8 coq6
  432.        if(isup.eq.1) then
  433.       call cq8loc(xel,nbn1,MINTE1.SHPTOT,TXR,irr)
  434.        else
  435.       call cq8loc(xel,nbn1,SHPTOT,TXR,irr)
  436.        endif
  437.  
  438. C      write(6,*) ' normales aux points ad hoc  coques epaisses'
  439.         do i=1,nbn1
  440.         do k=1,3
  441.         vecn(k,i)=txr(k,3,i)
  442.         enddo
  443.         enddo
  444. 2222  format(I4,3E12.5)
  445. 2001  format(3E10.3)
  446.         endif
  447. C--------------------------------------------------
  448. *     write(6,*)  ' normales   a l element  '
  449. *     write(6,fmt='(3E12.5)') ((vecn(l,ip),l=1,3),ip=1,nbn1)
  450. C--------------------------------------------------
  451. C
  452.         do 603 ipoi=1,nnb
  453.         if(melva1.ne.0) then
  454.         ibmn =MIN(iel,melva1.velche(/2))
  455.         igmn =MIN(ipoi,melva1.velche(/1))
  456.         epai = melva1.velche(igmn,ibmn)
  457.         else
  458.         epai=0.D0
  459.         endif
  460. C
  461.         if(melva2.ne.0) then
  462.         ibmn =MIN(iel,melva2.velche(/2))
  463.         igmn =MIN(ipoi,melva2.velche(/1))
  464.         exce = melva2.velche(igmn,ibmn)
  465.         else
  466.         exce =0.D0
  467.         endif
  468. C  si on projette aux points d integration  calcul des coordonnes
  469. C                sur l element reel
  470. *  on range les points sur les normales et non   pas par nappes
  471. *  pour ne pas chercher n fois exce et epai attention au MCHAML
  472. *
  473.         do 603 inap =1,ibnap(isous)
  474.         ipo=ipo+1
  475.         inupo = inupo+1
  476.                 if(ibnap(isous).gt.1)  then
  477. C                if(npint.ne.0.or.mfr.eq.5) then
  478.                 igau=(inap-1)*inbg(isous)/ibnap(isous)+1
  479.                 zzz = dzegau(igau)
  480.                 else
  481. C               zzz = (1.D0*inap-2.D0)
  482.                 zzz = XZERO
  483.                 endif
  484.         dnor=exce+epai*zzz/2.D0
  485.         ipt1.num(ipo,iel) =inupo
  486.         do 97 i=1,idim
  487. 97      XCOOR((inupo-1)*(IDIM+1)+i) = xeg(i,ipoi)+vecn(i,ipoi)*dnor
  488. C       write(6,2003) iel,ipoi,inupo-nbnoe,
  489. C     & ( XCOOR((inupo-1)*(IDIM+1)+i),i=1,3),dnor,igau
  490. 603     continue
  491. C        write(6,*) 'ipt1 iel ipo nbp crees ' ,
  492. C     & ipt1,iel,inap,(inupo-nbnoe)
  493.  
  494. 95     CONTINUE
  495. 2003   format(3I4,2X,4e12.5,I4)
  496. C
  497. C
  498.       ELSEIF(ICOCOQ(ISOUS).EQ.0) THEN
  499. C---------------------------------------------------------
  500. C ---------------------------- LA SZ  EST UN MASSIF ------
  501. C---------------------------------------------------------
  502.         IF(ISUP.GT.1) THEN
  503. C ----------------     on travaille aux points d integration
  504.         minte=iminte(isous)
  505. 2234    format(6e12.5)
  506. C
  507.       DO 195 iel=1,nbele1
  508.        call doxe(xcoor,idim,nbn1,num,iel,xel)
  509. C  ------ PROJECTION
  510.       DO 700 ip=1,shptot(/3)
  511.                 DO ik=1,idim
  512.                 XU(ik)=0.D0
  513.                 DO  jk = 1,NBN1
  514.                 XU(Ik) = XU(Ik) + SHPTOT(1,JK,ip)*(XEL(IK,jk))
  515.                 ENDDO
  516.                 ENDDO
  517.        inupo=inupo+1
  518.         ipt1.num(ip,iel) = inupo
  519.                 do  i=1,idim
  520.                 XCOOR((inupo-1)*(IDIM+1)+i) = xu(i)
  521.                 enddo
  522. C       write(6,2223) inupo,(xu(i),i=1,idim)
  523. 700   CONTINUE
  524. 195   CONTINUE
  525. 2223  format(I6,2X,3E12.5)
  526.        ELSE
  527. C     on travaille aux noeuds
  528.        do 197 iel=1,nbele1
  529.        call doxe(xcoor,idim,nbn1,num,iel,xel)
  530.                 do j=1,nbn1
  531.                 inupo=inupo+1
  532.                 ipt1.num(j,iel) = inupo
  533.                 do i=1,idim
  534.                 XCOOR((inupo-1)*(IDIM+1)+i) = xel(i,j)
  535.                 enddo
  536.                 enddo
  537. 197    continue
  538.        ENDIF
  539. C    fin de  l'alternative   coques ou massifs
  540.       ENDIF
  541. C--------------------------------
  542.       SEGSUP WRK4,VECT
  543. *      segsup info
  544. *
  545.  100  CONTINUE
  546. C   fin de creation des noeuds supports provisoires
  547. C        ENDIF
  548. C--------------//////////////////////////////////------------------
  549. C  maintenant il faut un meleme poi1 de tous les points sur lesquels
  550. C          on doit interpoler une valeur
  551.          segini,ipt4=ipt2
  552.         ipgeom=ipt4
  553.         if( isup.ne.2) then
  554.         call change(ipgeom,1)
  555.         segsup ipt4
  556.         endif
  557. C       call ecrobj('MAILLAGE',ipgeom)
  558. C----------------------------------------------------------------
  559. C----------------------------------------------------------------
  560. *-------------on est pret a faire l interpolation sur ipgeom ----
  561. C
  562. *write(6,*) ' maillage apres changer poi1 ' ,ipgeom
  563. C       call ecmail(ipgeom)
  564.         isort=0
  565.  
  566. *        write(6,*) ' dans pron appel a pro2'
  567.           CALL PRO2(IPGEOM,IPCHT,isort,IPOUT,ilphmo)
  568. *      write(6,*) ' dans pron sortie de pro2',nsous,ilphmo
  569.          if (ierr.ne.0) return
  570. C----------------------------------------------------------------
  571. C----------------------------------------------------------------
  572. C----------------------------------------------------------------
  573. C  ---      il faut maintenant reconstituer les MCHAML
  574. C  ---       a partir du chpo  construit  sur ipgeom ----------------
  575. C------------ INITIALISATION   du MCHAML RESULTAT ------------------
  576.        mlchpo = ipout
  577.        nsouz  = nsous
  578.        segini phachm
  579.         N1=NSOUS - ISURE
  580.         L1=12
  581.         N3=6
  582.         SEGINI MCHEL1
  583.         MCHEL1.TITCHE='SCALAIRE'
  584.         MCHEL1.IFOCHE=IFOUR
  585.  
  586. * boucle phases modele
  587.        DO 7000 IPHAS = 1,JGM
  588. *         write(6,*) 'iphas ',jgm,iphas, ' mots ',mots(iphas),ichpoi(/1)
  589.         segini icpr,snomip
  590. *        mdata=isort
  591.         mchpoi=ichpoi(iphas)
  592.          if (mchpoi.le.0) goto 7000
  593.         nsschp=ipchp(/1)
  594.         segini ipcb
  595. C
  596.         do i=1,nsschp
  597. *        mchpoi=ipca(i)
  598. C       call ecchpo( mchpoi)
  599.         msoupo=ipchp(i)
  600.         ipcb(i)=msoupo
  601. C       write(6,*) i,ipca(i),msoupo
  602.         mpoval=ipoval
  603.         nc = nocomp(/2)
  604. C
  605. C tableau des general des composantes
  606.         do ic = 1,nc
  607.         nomm= nocomp(ic)
  608.         if(nomip(/2).eq.0) then
  609.           nomip(**)=nomm
  610.         else
  611.           do  k=1,nomip(/2)
  612.           if(nomm.eq.nomip(k)) goto 4322
  613.         enddo
  614.         nomip(**)=nomm
  615. 4322    continue
  616.         endif
  617.         enddo
  618. *       write(6,*) (nomip(l),l=1,nomip(/2))
  619. *    reperage de la position  ses points dans le chpo
  620.         ipt5=igeoc
  621.         do j=1,ipt5.num(/2)
  622.         icpr(1,ipt5.num(1,j))=j
  623.         icpr(2,ipt5.num(1,j))= i
  624.         enddo
  625.         enddo
  626. C  i=icpr(2,k)   point k venant du msoupo du msoupo i  dans ipcb
  627.  
  628.         kphach = 0
  629.         inupo=0
  630.         DO 2100 isous = 1,nsous
  631.         imodel= KMODEL(isous)
  632.          if (conmod(17:24).ne.mots(iphas)) goto 2100
  633.         meleme = imamod
  634.          if (itypel.eq.48) goto 2100; 
  635. C
  636.         ipt5=ipt2.lisous(isous)
  637.         N2 = nomip(/2)
  638.         mfr=imfr(isous)
  639. * pour une phase donnée ne cree le mchaml qu une fois par maillage
  640. * distinct
  641.  
  642.         if (kphach.gt.0) then
  643.           do kpha = 1,kphach
  644.             if (imamod.eq.iphach(kpha)) then
  645.               mchaml = mphach(kpha)
  646.               logach = .true.
  647.               goto 2105
  648.             endif
  649.           enddo
  650. *             création du nouveau segment MCHAML
  651.               SEGINI MCHAML
  652.               kphach = kphach + 1
  653.               iphach(kphach) = imamod
  654.               mphach(kphach) = mchaml
  655.               logach = .false.
  656.         else
  657. *             création du nouveau segment MCHAML
  658.               SEGINI MCHAML
  659.               kphach = kphach + 1
  660.               iphach(kphach) = imamod
  661.               mphach(kphach) = mchaml
  662.               logach = .false.
  663.         endif
  664.  
  665.  2105         MCHEL1.IMACHE(isous)=MELEME
  666.               MCHEL1.ICHAML(isous)=MCHAML
  667.               MCHEL1.CONCHE(isous)=CONMOD
  668.               MCHEL1.INFCHE(isous,1)=0
  669.               MCHEL1.INFCHE(isous,2)=0
  670.               MCHEL1.INFCHE(isous,3)=NIFOUR
  671.           IF(isup.eq.1) then
  672.               MCHEL1.INFCHE(isous,4)=0
  673.               MCHEL1.INFCHE(isous,5)=0
  674.               MCHEL1.INFCHE(isous,6)=1
  675.           ELSE
  676.               MCHEL1.INFCHE(isous,4)=iminte(isous)
  677.               MCHEL1.INFCHE(isous,5)=0
  678.               MCHEL1.INFCHE(isous,6)=isup
  679.           ENDIF
  680. *
  681.         if (logach) goto 2100
  682.  
  683.         N1EL=NUM(/2)
  684.         N2PTEL=0
  685.         N2EL=0
  686. C
  687.         IF(ICOCOQ(ISOUS).EQ.1) THEN
  688. C--------------------------------------------------------
  689. C       boucle sur les composantes   et c est une coque
  690. C--------------------------------------------------------
  691.         inocom=0
  692.         DO  2200 icomp=1,nomip(/2)
  693.         inocom=inocom+1
  694. *write(6,*) 'composante ', nomip(icomp)
  695.         idebco = inupo
  696.  
  697. C         N1PTEL= NUM(/1)*nbnap
  698.           N1PTEL= nnb*nbnap
  699.           SEGINI MELVAL
  700.           IELVAL(inocom)= MELVAL
  701.           NOMCHE(inocom)= NOCOMP(icomp)
  702.           TYPCHE(inocom)= 'REAL*8'
  703. C
  704.           DO 162 NUEL=1,num(/2)
  705.           DO 162 NUPT=1,nnb
  706.              DO 172 IPOS = 1,nbnap
  707.              ipop = (ipos-1)*nnb+nupt
  708.              inupo = ipt5.num(ipop,nuel)
  709.        jh=icpr (1,inupo)
  710.        ipa=icpr(2,inupo)
  711.        if(jh.eq.0) go to 172
  712. C  il faut verifier si vpocha existe pour ce point
  713.        msoupo = ipcb(ipa)
  714.        mpoval=ipoval
  715.        do l=1,nocomp(/2)
  716.        if(nocomp(l).eq.nomip(icomp)) then
  717.        vvv=vpocha(jh,icomp)
  718.        goto  4557
  719.        endif
  720.        enddo
  721.        vvv=XZERO
  722. 4557   continue
  723. C       write(6,2004) NUEL,inupo,jh,ipop,
  724. C     & ( XCOOR((inupo+nbnoe-1)*(IDIM+1)+i),i=1,3),vvv
  725.               VELCHE(ipop,NUEL) = vvv
  726.  172      CONTINUE
  727.  162      CONTINUE
  728.  
  729.         if(icomp.lt.nocomp(/2)) inupo =idebco
  730. 2200    continue
  731.  
  732.         ELSE
  733. C--------------------------------------------------------
  734. C       boucle sur les composantes   et c est un massif
  735. C--------------------------------------------------------
  736.  
  737.         if(isup.eq.1) then
  738.         N1PTEL=NUM(/1)
  739.         elseif(isup.eq.2) then
  740.         N1PTEL=1
  741.         else
  742.         N1PTEL= inbg(isous)
  743.         endif
  744. C
  745.         DO  2220 icomp=1,NOMIP(/2)
  746. *        write(6,*) 'composante ', nomip(icomp)
  747.         idebco = inupo
  748. C
  749.  
  750.         SEGINI MELVAL
  751.         NOMCHE(icomp)=nomip(icomp)
  752.         TYPCHE(icomp)='REAL*8'
  753.         IELVAL(icomp)=MELVAL
  754. C
  755.           DO 163 NUEL=1,N1EL
  756.           DO 163 NUPT=1,N1PTEL
  757.               inupo=ipt5.num(nupt,nuel)
  758.        jh=icpr(1,inupo)
  759.        ipa=icpr(2,inupo)
  760.        if(jh.eq.0) go to 163
  761. C  il faut verifier si vpocha existe pour ce point
  762.        msoupo = ipcb(ipa)
  763.        mpoval=ipoval
  764.        do l=1,nocomp(/2)
  765.        if(nocomp(l).eq.nomip(icomp)) then
  766.        vvv=vpocha(jh,icomp)
  767.        goto  4558
  768.        endif
  769.        enddo
  770.        vvv=XZERO
  771. 4558   continue
  772.  
  773. C        write(6,2003) NUEL,inupo,jh,
  774. C     & ( XCOOR((inupo+nbnoe-1)*(IDIM+1)+i),i=1,3),vvv
  775.               VELCHE(NUPT,NUEL) = vvv
  776.  163      CONTINUE
  777.         ipcham= mchaml
  778.         ipc1 = melval
  779.         call comred(ipc1)
  780.         IELVAL(icomp)=ipc1
  781.         melval = ipc1
  782.         if(icomp.lt.nocomp(/2)) inupo =idebco
  783. C  ----  fin de la boucle sur les composantes
  784. 2220    continue
  785.  
  786. C----------------  fin du traitement des massifs
  787.         ENDIF
  788.         segsup ipt5
  789. C    fin de la boucle sur les sous zones du modele
  790. 2100    continue
  791. C   destruction des chpo intermediaires
  792.         do i=1,ipcb(/1)
  793.           msoupo=ipcb(i)
  794. cgf          ipt5= igeoc (correction 7284)
  795.           mpoval= ipoval
  796. *          mchpoi=ipca(i)
  797.           segsup mpoval,msoupo
  798. cgf          segsup ipt5 (correction 7284)
  799.         enddo
  800.         segsup mchpoi
  801.         segsup ipcb,snomip
  802.         segsup icpr
  803.  7000   CONTINUE
  804. C (fdp) re-ajustement de MCOORD a sa taille initiale
  805.         NBPTS = NBNOE
  806.         SEGADJ MCOORD
  807. C       retrait des maillages temporaires du pre-conditionnement
  808. c       (leurs numero de noeuds depasse la taille de MCOOR)
  809.         call crech1b
  810. C (fdp) suppression du maillage temporaire IPT1
  811.         ipt1=ipgeom
  812.         segsup ipt1
  813.         segsup ipt2
  814.         segsup phachm
  815.         IPOUT=MCHEL1
  816.         segsup sicoq
  817.  
  818.         return
  819.       END
  820.  
  821.  
  822.  
  823.  
  824.  
  825.  

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