Télécharger la plaquette Cast3M


fr en

Télécharger prom.eso

Retour à la liste

Numérotation des lignes : 

prom

  1. C PROM      SOURCE    CB215821  24/04/12    21:16:57     11897          
  2.  
  3.       SUBROUTINE PROM(IPMODE,IPCHE,IPCHT,ISUP,IPOUT)
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6.       EXTERNAL SHAPE
  7. *-------------------------------------------------------------------- *
  8. *                                                                     *
  9. *   Sous-programme appelé par l'opérateur   PROI                      *
  10. *   Projection d'un chamelem sur le support d un modele  par          *
  11. *   minimisation de l integrale   (Ue-Us)**2 ds  sur chaque element   *
  12. *    Ue evalue a partir des fonctions de formes du modele associe     *
  13. *        au champ d origine si il est  fourni ( sinon a partir des    *
  14. *        fonctions de transformation geometriques de l element)       *
  15. *                                                                     *
  16. *    L integrale est exprime à partir des fonctions de forme associees*
  17. *    au modèle support recepteur                                      *
  18. *---------------------------------------------------------------------*
  19. *                                                                     *
  20. *   ENTREE :                                                          *
  21. *          IPMODE  modele sur lequel on va projeter                   *
  22. *          IPCHE   MCHAML de caracteristique (si le maodele comprends *
  23. *                  des coques)                                        *
  24. *         IPCHT   MCHAML a projeter (obligatoirement defini aux noeuds)
  25. *         ISUP    remis  a 1  pour le moment (on travaille aux noeuds)*
  26. *                                                                     *
  27. *   SORTIE :                                                          *
  28. *          IPOUT   MCHAML resultat                                    *
  29. *                                                                     *
  30. *---------------------------------------------------------------------*
  31. *
  32.  
  33. -INC PPARAM
  34. -INC CCOPTIO
  35. -INC SMCOORD
  36. -INC SMMODEL
  37. -INC SMCHAML
  38. -INC SMELEME
  39. -INC SMCHPOI
  40. -INC SMLCHPO
  41. -INC SMLMOTS
  42. -INC TMTRAV
  43. -INC SMINTE
  44. -INC CCREEL
  45. *
  46.         PARAMETER ( O1=1.D0,O2=2.D0)
  47.         character*(LOCOMP) nomm
  48.         segment snomip
  49.           character*(LOCOMP) nomip(0)
  50.         endsegment
  51. *        SEGMENT MDATA
  52. *         INTEGER IPCA(*),IPEXI(*)
  53. *        ENDSEGMENT
  54.         segment ipcb(nsschp)
  55.         segment icpr(2,nbpts)
  56.         segment icpr1(nbpts)
  57. C
  58.       segment info
  59.       integer infell(JG)
  60.       endsegment
  61. C
  62.         segment sicoq
  63.         integer icocoq(nbsous),imfr(nbsous),iiele(nbsous)
  64.         integer iminte(nbsous),imint1(nbsous),ibnap(nbsous)
  65.         integer inbg(nbsous),ibbas(nbsous)
  66.         endsegment
  67. C
  68.        segment WRK5
  69.         REAL*8 W(npt),A(NBN1,NBN1),B(NBN1,NBN1),XPU(3),RHS(NBN1)
  70.         REAL*8 SHPXXX(6,NBNN,npt),XX(3,npt),VAL(NBN1),POIDS(npt)
  71.         INTEGER NDD
  72.         endsegment
  73. C
  74.       SEGMENT WRK4
  75.        REAL*8 BPSS(3,3),XEL(3,NBN1),XU(3),SHP(6,NBN1),S(6)
  76.        REAL*8 TXR(3,3,NBN1),XE(3,NBG)
  77.       ENDSEGMENT
  78.       SEGMENT VECT
  79.         REAL*8  VEC1(IDIM)
  80.         REAL*8  VEC2(IDIM)
  81.         REAL*8  VECN(IDIM,ibbas(isous))
  82.       ENDSEGMENT
  83.  
  84. C= Logique a TRUE si mode de calcul est axisymetrique 1D/2D ou Fourier
  85.       LOGICAL LOGAXI
  86. C= Quelques constantes (2.Pi et 4.Pi)
  87.       PARAMETER (X2Pi =6.283185307179586476925286766559D0)
  88.       PARAMETER (X4Pi=12.566370614359172953850573533118D0)
  89.  
  90. C   on   remet  isup a 1  pour le moment
  91.       ISUP=1
  92.       MMODEL=IPMODE
  93.       SEGACT,MMODEL
  94.       NSOUS=KMODEL(/1)
  95. C
  96.       ndec=0
  97.       call lirent(ndec,0,iret)
  98.       if(ndec.eq.0) ndec = 3
  99.       if(ndec.lt.0.or.ndec.gt.10) then
  100.         moterr='integration'
  101.         call erreur(645)
  102.         return
  103.       endif
  104. *
  105.       segact mcoord*mod
  106.       NBNOE = nbpts 
  107.  
  108.       if(ipche.NE.0) then
  109.       MCHELM=IPCHE
  110.       segact MCHELM
  111.       NCHAM=ICHAML(/1)
  112.       endif
  113. *
  114. *     Création du maillage provisoire support pour interpolation
  115. *
  116.       NBSOUS = NSOUS
  117.       NBREF = 0
  118.       NBELEM = 0
  119.       NBNN =  0
  120. C ipt2  contiendra les super couche ou points de gauss
  121.       SEGINI IPT2
  122.       SEGINI SICOQ
  123. c listmots des phases
  124.       ilphmo = -1
  125.       jgn = 8
  126.       jgm = nsous
  127.       segini mlmots
  128.       ilphmo = mlmots
  129.       jgm = 1
  130. *
  131. * Boucle sur l'ensemble des sous zones du modeles sur lequel on va
  132. *   projetter   pour   initialisations
  133.         nponou=0
  134.         npcoq =0
  135.       DO 30 ISOUS = 1, NSOUS
  136. *
  137.         IMODEL = KMODEL(ISOUS)
  138.         SEGACT,IMODEL
  139.         MELEME = IMAMOD
  140.         SEGACT MELEME
  141.         MELE = NEFMOD
  142.         NBELE1= NUM(/2)
  143.         NBN1 = NUM(/1)
  144.         ibnap(isous)=1
  145. C
  146. C   COQUE INTEGREE OU PAS ?
  147. C
  148.       IF(INFMOD(/1).NE.0)THEN
  149.           NPINT=INFMOD(1)
  150.           if( npint.ne.0)IBNAP(isous) =NPINT
  151.       ELSE
  152.           NPINT=0
  153. C          ibnap(isous)=1
  154.       ENDIF
  155.  
  156.       if( infmod(/1).lt.2+isup) then
  157.       call elquoi(mele,npint,isup,ipinf,imodel)
  158.       if (ierr.ne.0) then
  159.         segsup mchelm
  160.         iret=0
  161.         return
  162.       endif
  163.         info=ipinf
  164.         mfr  =infell(13)
  165.         iele =infell(14)
  166.         minte=infell(11)
  167.         minte1=infell(12)
  168.         nbg=infell(6)
  169.         segsup info
  170.          else
  171.         minte=infmod(isup+2)
  172.         minte1=infMOD(12)
  173.         nbg=infele(6)
  174.         mfr  =infele(13)
  175.         iele =infele(14)
  176.         endif
  177.         imfr(isous)=mfr
  178.         iiele(isous)=iele
  179.         iminte(isous)=minte
  180.         imint1(isous)=minte1
  181.         inbg(isous)=nbg
  182. *       write(6,*) 'mfr iele mele minte minte1 nbg',mfr,iele,
  183. *     & mele,minte,minte1,nbg
  184.         NBSOUS = 0
  185.         NBREF = 0
  186.         NBELEM = NBELE1
  187.  
  188.       IF(mfr.eq.3.or.mfr.eq.5.or.mfr.eq.9) THEN
  189. C------------------coques
  190.       IF(IPCHE.eq.0) THEN
  191. C  pour le coques il faut les caracteristiques
  192.       call erreur(404)
  193. *      segsup info
  194.       return
  195.       ENDIF
  196. C
  197.       icocoq(isous)=1
  198.       nbnn = ndec*ndec
  199.       if(iele.lt.4) nbnn=ndec
  200.       ibbas(isous)= nbnn
  201.        if(idim.eq.3.or.(iele.lt.4.and.idim.eq.2))  then
  202.        ibnap(isous)=ndec
  203.        nbnn= nbnn*ndec
  204.        inbg(isous)=nbnn
  205.        endif
  206.  
  207.       ELSE
  208. C--------------- massifs
  209.  
  210.         IF (mfr.EQ.1) THEN
  211.           IF (IDIM.EQ.3) THEN
  212.             NBNN=ndec*ndec*ndec
  213.           ELSE IF (IDIM.EQ.2) THEN
  214.             NBNN=ndec*ndec
  215. C*        ELSE IF (IDIM.EQ.1) THEN
  216.           ELSE
  217.             NBNN=ndec
  218.           ENDIF
  219.         ELSE
  220.           NBNN=ndec*ndec
  221.         ENDIF
  222.  
  223. C    petite rectif utile pour les elements speciaux
  224.         if(iele.eq.32) then
  225. C   les    polygones
  226.         nbnn =nbnn*nbn1
  227.         elseif(iele.eq.25.or.iele.eq.26) then
  228. C  les pyramides
  229.         nbnn = nbnn*4
  230.         endif
  231. C
  232.         inbg(isous)= NBNN
  233.       ENDIF
  234. *
  235.         SEGINI IPT1
  236.         IPT1.ITYPEL = 28
  237.         ipt2.lisous(isous) = ipt1
  238.         nponou=nponou+nbnn*nbelem
  239.  
  240.         if (isous.eq.1) then
  241.          mots(1) = conmod(17:24)
  242.         else
  243.          do ipl = 1,jgm
  244.           if (mots(ipl).eq.conmod(17:24)) goto 27
  245.          enddo
  246.          jgm = jgm + 1
  247.          mots(jgm) = conmod(17:24)
  248.  27      continue
  249.         endif
  250.  
  251. *        segsup info
  252. 30    CONTINUE
  253.       segadj mlmots
  254. C-----------------------------------------------------------
  255. C
  256. C
  257. C   on va fabriquer des points support provisoires a l interieur
  258. C        de chaque element des maillages du modele  recepteur
  259. C
  260. C
  261. C----------------------------------------------------------
  262. C        write(6,*) ' nombre de points crees ' ,nponou
  263. C (DIP) ajustement de MCOORD pour les noeuds provisoires
  264. C       il sera remis a sa taille initiale a la fin du programme
  265.          NBPTS = NBNOE + nponou
  266.          SEGADJ MCOORD
  267. *-----------------------------------------------------------
  268. *     Boucle sur l'ensemble des sous zones du modeles
  269. *------------------------------------------------------------
  270.         inupo =NBNOE
  271.       DO 100 ISOUS = 1,NSOUS
  272.         IMODEL = KMODEL(ISOUS)
  273.         segact imodel
  274.         MELEME = IMAMOD
  275.         segact meleme
  276.         nbele1= num(/2)
  277.         nbn1 = num(/1)
  278.         ipt1=ipt2.lisous(isous)
  279.         segact ipt1*mod
  280.         NNB = ibbas(isous)
  281.         NBG = inbg(isous)
  282.         SEGINI WRK4
  283.             MELVA1 = 0
  284.             MELVA2 = 0
  285.       IF(icocoq(isous).ne.0) THEN
  286. C  -------------------------------la sz est une coque ------
  287.  
  288.            NUCHA = 0
  289.        DO 15, NUCH = 1, NCHAM
  290. *
  291.             IF ( CONCHE(NUCH).EQ.CONMOD.AND.
  292.      &           IMACHE(NUCH).EQ.IMAMOD) NUCHA = NUCH
  293. *
  294.   15      CONTINUE
  295. *
  296.           IF (NUCHA.NE.0) THEN
  297.             MCHAML=ICHAML(NUCHA)
  298.             SEGACT,MCHAML
  299. *
  300.             MELVA1 = 0
  301.             MELVA2 = 0
  302.             NCOMP = IELVAL(/1)
  303.             DO 20, I = 1, NCOMP
  304.               IF (NOMCHE(I).EQ.'EPAI') THEN
  305.                 MELVA1 = IELVAL(I)
  306.                 SEGACT, MELVA1
  307.               ELSEIF (NOMCHE(I).EQ.'EXCE') THEN
  308.                 MELVA2 = IELVAL(I)
  309.                 SEGACT, MELVA2
  310.               ENDIF
  311.   20        CONTINUE
  312.           ENDIF
  313.       ENDIF
  314. C
  315.       MELE = NEFMOD
  316. c   coque integree ou pas ?
  317.         if(infmod(/1).ne.0)then
  318.         npint=infmod(1)
  319.         else
  320.         npint=0
  321.         endif
  322. C
  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. C
  330.         segact minte
  331.         if(mfr.eq.5) segact minte1
  332.  
  333.         nnb = ibbas(isous)
  334. C
  335.        ipo=0
  336. *-------------------------------------------------------------------
  337. *   Pavage de l element de reference
  338. *-------------------------------------------------------------------
  339.        npt=1000
  340.        NBNN = inbg(isous)
  341.        segini WRK5
  342.        iwkr= wrk5
  343.        iele=iiele(isous)
  344.        CALL PAVPOI(WRK5,iele,ndec,nbn1)
  345.        segini vect
  346. *-------------------------------------------------------------------
  347. *  Boucle sur les elements de la sous zone du modele
  348. *-------------------------------------------------------------------
  349.  
  350.       DO IEL=1,NBELE1
  351.       CALL DOXE(XCOOR,IDIM,NBN1,NUM,IEL,XEL)
  352.  
  353. C   si c est une coque pouvant etre excentree
  354. C     recherche des normales  aux noeuds
  355.       if(icocoq(isous).ne.0)  then
  356. C
  357.        IF(IELE.EQ.4) THEN
  358. C    --------coq3   dkt  dst   ----------------
  359.       do  i=1,idim
  360.       VEC1(i)=xel(i,2)-xel(i,1)
  361.       VEC2(i)=xel(i,3)-xel(i,1)
  362.       enddo
  363.       vecn(1,1) = vec1(2)*vec2(3)-vec1(3)*vec2(2)
  364.       vecn(2,1) = vec1(3)*vec2(1)-vec1(1)*vec2(3)
  365.       vecn(3,1) = vec1(1)*vec2(2)-vec1(2)*vec2(1)
  366.       vnor =sqrt(vecn(1,1)*vecn(1,1)+vecn(2,1)*vecn(2,1)
  367.      & +vecn(3,1)*vecn(3,1))
  368.       do    i=1,idim
  369.       vecn(i,1)=vecn(i,1)/vnor
  370.       enddo
  371.       do    j=2,ibbas(isous)
  372.       do    i=1,idim
  373.       vecn(i,j)=vecn(i,1)
  374.       enddo
  375.       enddo
  376.                 ELSEIF (IELE.EQ.8) THEN
  377. C    --------------- coq4 ----------------
  378.       call cq4loc(xel,xe,BPSS,irrt,1)
  379. C  ici c est une estimation du plan moyen
  380. *     write(6,fmt='(3E12.5)') ((bpss(i,j),i=1,3),j=1,3)
  381.       do i=1,3
  382.       do ip=1,ibbas(isous)
  383.       vecn(i,ip)=bpss(3,i)
  384.       enddo
  385.       enddo
  386.         ELSEIF(IELE.eq.2) THEN
  387. C   ---------------- coq2 --------------
  388.         vnor=0.D0
  389.         do    i=1,idim
  390.         VEC1(i)=xel(i,2)-xel(i,1)
  391.         vnor=vnor + vec1(i)*vec1(i)
  392.         enddo
  393.         vnor = sqrt(vnor)
  394.         vecn(1,1)=-vec1(2)/vnor
  395.         vecn(2,1)=vec1(1)/vnor
  396.         vecn(1,2)=vecn(1,1)
  397.         vecn(2,2)=vecn(2,1)
  398.         ELSEIF (iele.eq.5.or.iele.eq.10) THEN
  399. C  preparation pour coques integrees
  400.          SURFX=XZERO
  401.          SURFY=XZERO
  402.          SURFZ=XZERO
  403.                  ENDIF
  404.       endif
  405. C -------------------  boucle sur les points d integration
  406.       inuo=inupo
  407. C      write(6,*) ' ndd et inupo ' ,ndd,inupo
  408.       do igau=1,ndd
  409. C
  410.       do j=1,idim
  411.       xpu(j)=xx(j,igau)
  412.       enddo
  413. C      write(6,*)(xpu(J),j=1,idim)
  414. C
  415.       if(mele.GE.111.and.mele.le.122) then
  416. C  les plolygones
  417.       call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  418.       else
  419.       CALL SHAPE(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  420.       endif
  421.  
  422. C       PROJECTION  DANS L ELEMENT REEL
  423.       DO     IK = 1,IDIM
  424.       XU(IK)= XZERO
  425.       DO     JK = 1,nbn1
  426.       XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK))
  427.       enddo
  428.       enddo
  429. C
  430.       if ( mfr.eq.5) then
  431. C             coques integrees ds l epaisseur
  432.             DO I=1,6
  433.                S(I)=XZERO
  434.             CONTINUE
  435.             ENDDO
  436.             DO INOE=1,NBN1
  437.                S(1)=S(1)+SHP(2,INOE)*XEL(2,INOE)
  438.                S(2)=S(2)+SHP(3,INOE)*XEL(3,INOE)
  439.                S(3)=S(3)+SHP(3,INOE)*XEL(2,INOE)
  440.                S(4)=S(4)+SHP(2,INOE)*XEL(3,INOE)
  441.                S(5)=S(5)+SHP(3,INOE)*XEL(1,INOE)
  442.                S(6)=S(6)+SHP(2,INOE)*XEL(1,INOE)
  443.             ENDDO
  444.             SURFX=S(1)*S(2)-S(3)*S(4)
  445.             SURFY=S(4)*S(5)-S(2)*S(6)
  446.             SURFZ=S(6)*S(3)-S(5)*S(1)
  447.             SURF=SQRT(SURFX**2+SURFY**2+SURFZ**2)
  448. C      write(6,*) ' normales au pt igau du plan de base '
  449.         vecn(1,igau)=surfx/surf
  450.         vecn(2,igau)=surfy/surf
  451.         vecn(3,igau)=surfz/surf
  452.       endif
  453. C
  454.       IF(ICOCOQ(ISOUS).EQ.0) THEN
  455. C   elements  standards
  456.       inupo=inupo+1
  457.       ipt1.num(igau,iel)=inupo
  458.       do i=1,idim
  459.       xcoor((inupo-1)*(idim+1)+i)=XU(i)
  460.       enddo
  461. C      write(6,3335) (xu(i),i=1,idim)
  462. C3335  format(3e12.5)
  463.       ELSE
  464. C   traitement special des coques
  465.         if(melva1.ne.0) then
  466.         ibmn =MIN(iel,melva1.velche(/2))
  467.         igmn =MIN(igau,melva1.velche(/1))
  468.         epai = melva1.velche(igmn,ibmn)
  469.         else
  470.         epai =xzero
  471.         endif
  472. C
  473.         if(melva2.ne.0) then
  474.         ibmn =MIN(iel,melva2.velche(/2))
  475.         igmn =MIN(igau,melva2.velche(/1))
  476.         exce = melva2.velche(igmn,ibmn)
  477.         else
  478.         exce =xzero
  479.         endif
  480.  
  481.       dp = epai/ibnap(isous)
  482.       dnor = exce -(epai + dp)/2.D0
  483. C
  484.       do inap=1,ibnap(isous)
  485.       dnor=dnor+dp
  486.       ipos = (inap-1)*ibbas(isous)+igau
  487.       ipt1.num(ipos,iel)=inuo+ipos
  488.       do i=1,idim
  489.       xcoor((inuo+ipos-1)*(idim+1)+i)= XU(i)+dnor*vecn(i,igau)
  490.       enddo
  491.       enddo
  492.       inupo=inupo+ibnap(isous)
  493.       ENDIF
  494.  
  495. C ---------------fin de la boucle sur les points de d integration
  496.       enddo
  497.  
  498. C  fin de la boucle sur les elements
  499.       enddo
  500.        segsup wrk5,vect
  501.       SEGSUP WRK4
  502. *      segsup info
  503. *
  504.  100  CONTINUE
  505. C   fin de creation des noeuds supports provisoires
  506. 2003   format(3I4,2X,4e12.5,I4)
  507. C--------------//////////////////////////////////------------------
  508. C  maintenant in faut un meleme poi1  de tous les points sur lesquels
  509. C          on doit interpoler une valeur
  510.         segini icpr1
  511.         ipt5=ipt2
  512.          nbelem=0
  513.         do io=1,max(1,ipt2.lisous(/1))
  514.          if (ipt2.lisous(/1).ne.0) ipt5=ipt2.lisous(io)
  515.          segact ipt5
  516.          do ip=1,ipt5.num(/1)
  517.          do il=1,ipt5.num(/2)
  518.           ipt=ipt5.num(ip,il)
  519.           if (icpr1(ipt).eq.0) then
  520.            nbelem=nbelem+1
  521.            icpr1(ipt)=nbelem
  522.           endif
  523.          enddo
  524.          enddo
  525.         enddo
  526.         nbnn=1
  527.         nbsous=0
  528.         nbref=0
  529. **      write (6,*) ' nbelem ',nbelem
  530.         segini ipt4
  531.         ipgeom=ipt4
  532.         nbelem=0
  533.         do i=1,icpr1(/1)
  534.          if (icpr1(i).ne.0) then
  535.           nbelem=nbelem+1
  536. **        write (6,*) ' i nbelem icpr ',i,nbelem,icpr1(i)
  537.           ipt4.num(1,nbelem)=i
  538.          endif
  539.         enddo
  540.         segsup icpr1
  541. C----------------------------------------------------------------
  542. C----------------------------------------------------------------
  543. *-------------on est pret a faire l interpolation sur ipgeom ----
  544. C
  545. *write(6,*) ' maillage apres changer poi1 ' ,ipgeom
  546. C          call ecmail(ipgeom)
  547.          isort=0
  548.           CALL PRO2(IPGEOM,IPCHT,isort,IPOUT,ilphmo)
  549.          if (ierr.ne.0) return
  550. C----------------------------------------------------------------
  551. C----------------------------------------------------------------
  552. C
  553. C Definition de "constantes" selon le mode de calcul
  554. C Modes de calcul 1D/2D axisymetrique
  555.       IF (IFOMOD.EQ.0.OR.IFOMOD.EQ.4) THEN
  556.         LOGAXI=.TRUE.
  557.         CteAxi=X2Pi
  558. C Mode de calcul 2D Fourier
  559.       ELSE IF (IFOMOD.EQ.1) THEN
  560.         LOGAXI=.TRUE.
  561.         IF (NIFOUR.EQ.0) THEN
  562.           CteAxi=X2Pi
  563.         ELSE
  564.           CteAxi=XPi
  565.         ENDIF
  566.       ELSE
  567.         LOGAXI=.FALSE.
  568.       ENDIF
  569. C----------------------------------------------------------------
  570. C  ---      il faut maintenant reconstituer les MCHAML
  571. C  ---       a partir du chpo  construit  sur ipgeom ----------------
  572.        mlchpo = ipout
  573.        segact mlchpo
  574.         mchel5= ipcht
  575.         segact mchel5
  576.         N1=NSOUS
  577.         L1=12
  578.         N3=6
  579.         SEGINI MCHEL1
  580. *        MCHEL1.TITCHE='SCALAIRE'
  581.         MCHEL1.TITCHE = mchel5.titche
  582.         MCHEL1.IFOCHE=IFOUR
  583. *//////////////////////////////
  584.  
  585. * boucle phases modele
  586.        DO 7000 IPHAS = 1,JGM
  587. *         write(6,*) 'iphas ',iphas
  588.         segini icpr,snomip
  589. *        mdata=isort
  590. *        segact mdata
  591.         mchpoi=ichpoi(iphas)
  592.         segact mchpoi
  593.         nsschp=ipchp(/1)
  594.         segini ipcb
  595.  
  596. C        ip=0
  597.         do i=1,nsschp
  598. *        mchpoi=ipca(i)
  599. C        call ecchpo(mchpoi)
  600. *        segact mchpoi
  601.         msoupo=ipchp(i)
  602.         ipcb(i)=msoupo
  603.         segact msoupo
  604.         mpoval=ipoval
  605.         segact mpoval
  606.         nc = nocomp(/2)
  607. C
  608. C tableau des general des composantes
  609.         do ic = 1,nc
  610.         nomm= nocomp(ic)
  611.         if(nomip(/2).eq.0) then
  612.         nomip(**)=nomm
  613.         else
  614.         do  k=1,nomip(/2)
  615.         if(nomm.eq.nomip(k)) goto 4322
  616.         enddo
  617.         nomip(**)=nomm
  618. 4322    continue
  619.         endif
  620.         enddo
  621. *       write(6,*) (nomip(l),l=1,nomip(/2))
  622. *    reperage de la position  ses points dans le chpo
  623.         ipt5=igeoc
  624.         segact ipt5
  625. *        ip=0
  626.         do j=1,ipt5.num(/2)
  627. *        ip=ip+1
  628.         icpr(1,ipt5.num(1,j))=j
  629.         icpr(2,ipt5.num(1,j))= i
  630.         enddo
  631.         enddo
  632. C  i=icpr(2,k)   point k venant du msoupo du msoupo i  dans ipcb
  633.  
  634.         inupo=0
  635.         segact ipt2
  636.         DO 2100 isous = 1,nsous
  637.         imodel= KMODEL(isous)
  638.         segact imodel
  639.          if (conmod(17:24).ne.mots(iphas)) goto 2100
  640.         meleme = imamod
  641.         segact meleme
  642.         idi2=idim
  643.         ipt5=ipt2.lisous(isous)
  644.         segact ipt5
  645.         if(icocoq(isous).ne.0) idi2 = idim-1
  646.         mfr=imfr(isous)
  647. *             création du nouveau segment MCHAML
  648.           N2=nomip(/2)
  649.               SEGINI MCHAML
  650.               MCHEL1.IMACHE(isous)=MELEME
  651.               MCHEL1.ICHAML(isous)=MCHAML
  652.               MCHEL1.CONCHE(isous)=CONMOD
  653.               MCHEL1.INFCHE(isous,1)=0
  654.               MCHEL1.INFCHE(isous,2)=0
  655.               MCHEL1.INFCHE(isous,3)=NIFOUR
  656.               MCHEL1.INFCHE(isous,4)=0
  657.               MCHEL1.INFCHE(isous,5)=0
  658.               MCHEL1.INFCHE(isous,6)=1
  659. *
  660.         N1EL=NUM(/2)
  661.         N2PTEL=0
  662.         N2EL=0
  663. C--------------------------------------------------------
  664. C       boucle sur les composantes
  665. C--------------------------------------------------------
  666.         NBNN=num(/1)
  667.         N1PTEL= NBNN
  668.         NBN1=NBNN
  669.  
  670. C
  671.         SEGINI WRK5,wrk4
  672.       iele = iiele(isous)
  673.       iwrk= wrk5
  674.       CALL PAVPOI(WRK5,iele,ndec,nbnn)
  675. C  -----  remplissage de shpxxx  ( shptot)
  676.       DO IGAU=1,NDD
  677.         do j=1,idim
  678.         xpu(j)=xx(j,igau)
  679.         enddo
  680.         if(mele.GE.111.and.mele.le.122) then
  681.         call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  682.         else
  683.         call shape(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  684.         endif
  685.       do  id=1,6
  686.       do  no=1,nbnn
  687.       shpxxx(id,no,igau)=shp(id,no)
  688.       enddo
  689.       enddo
  690.       ENDDO
  691.       N1PTEL=NBNN
  692. C ----------- creation des melvals recepteurs -----
  693.       DO   icomp=1,NOMIP(/2)
  694.         SEGINI MELVAL
  695.         NOMCHE(icomp)=nomip(icomp)
  696.         TYPCHE(icomp)='REAL*8'
  697.         IELVAL(icomp)=MELVAL
  698.       ENDDO
  699. C--------------------------------------------------
  700. C--------  boucle sur les elements       ----------
  701.         AXIS=O1
  702.         DO 2210 iel=1,num(/2)
  703.         CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XEL)
  704.         CALL ZERO(A,NBNN,NBNN)
  705. C       SSOM1=XZERO
  706. C
  707. C
  708.         DO inap=1,ibnap(isous)
  709.         DO 40 IGAU=1,ndd
  710.          do i=1,idim
  711.          xpu(i)=xx(i,igau)
  712.          enddo
  713.          if(mele.GE.111.and.mele.le.122) then
  714. C  les plolygones
  715.          call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret)
  716.          else
  717.          call shape(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  718.          endif
  719.          CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC)
  720. C Modes de calcul 1D/2D axisymetrique et 2D Fourier
  721.          IF (LOGAXI) THEN
  722.            CALL DISTRR(XEL,SHP,NBNN,RR)
  723.            AXIS=CteAxi*RR
  724. C Mode de calcul 1D spherique
  725.          ELSE IF (IFOMOD.EQ.5) THEN
  726.            CALL DISTRR(XEL,SHP,NBNN,RR)
  727.            AXIS=X4Pi*RR*RR
  728.          ENDIF
  729. C
  730. C   on veut calculer la fonction au point d integration  projete
  731. C           PROJECTION
  732.       DO IK = 1,IDIM
  733.       XU(IK)= XZERO
  734.       DO  JK = 1,NBNN
  735.       XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK))
  736.       ENDDO
  737.       ENDDO
  738. C
  739. **          precis=0.d0
  740.             CPZA=poids(igau)*DJAC*AXIS
  741.             DO 50 INO1=1,NBNN
  742.             DO 60 INO2=1,NBNN
  743.             Q = SHPXXX(1,INO1,IGAU)*SHPXXX(1,INO2,IGAU)*CPZA
  744. CCC         Q = SHPXXX(1,INO1,IGAU)*SHPXXX(1,INO2,IGAU)
  745. CCC  1              *poids(igau)*DJAC*AXIS
  746.             A(INO1,INO2)=A(INO1,INO2)+Q
  747. **          precis=max(precis,abs(q))
  748.  60         CONTINUE
  749.  50         CONTINUE
  750.  40         CONTINUE
  751.          enddo
  752. *  la precision de invalm est maintenant relative
  753.          precis = 1d-10
  754. C
  755. C  ----  fin preparation matrice 'masse ' sur element courant
  756. C         do i=1,NBNN
  757. C        do j=1,NBNN     garder eventuellement pour calculer
  758. C        B(i,j)=A(i,j)   la difference sur les integrales
  759. C        enddo
  760. C        enddo
  761. C
  762.          call invalm(A,NBNN,NBNN,KERRE,precis)
  763.          if (kerre.ne.0) then
  764.           moterr(1:8)='transfert'
  765.           call erreur(700)
  766.          endif
  767.  
  768.  
  769.  
  770. C  ------------------  boucle  sur les composantes
  771.         DO 2220 icomp=1,nomip(/2)
  772.         idebco = inupo
  773.         do i=1,nbnn
  774.         RHS(i) =XZERO
  775.         enddo
  776. C  ---------  boucle  sur les points d integration
  777. C
  778.         AXIS = O1
  779.         DO 42 inap=1,ibnap(isous)
  780.         DO 41 IGAU=1,ndd
  781.          do i=1,idim
  782.          xpu(i)=xx(i,igau)
  783.          enddo
  784.             CALL SHAPE(xpu(1),xpu(2),xpu(3),iele,shp,iret)
  785.          CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC)
  786. C Modes de calcul 1D/2D axisymetrique et 2D Fourier
  787.          IF (LOGAXI) THEN
  788.            CALL DISTRR(XEL,SHP,NBNN,RR)
  789.            AXIS=CteAxi*RR
  790. C Mode de calcul 1D spherique
  791.          ELSE IF (IFOMOD.EQ.5) THEN
  792.            CALL DISTRR(XEL,SHP,NBNN,RR)
  793.            AXIS=X4Pi*RR*RR
  794.          ENDIF
  795.  
  796.        ipop=(inap-1)*ndd+igau
  797.        inupo=ipt5.num(ipop,iel)
  798.  
  799.        jh=icpr (1,inupo)
  800.        ipa=icpr(2,inupo)
  801. C  il faut verifier si vpocha existe pour ce point
  802.        msoupo = ipcb(ipa)
  803.        mpoval=ipoval
  804.                 do l=1,nocomp(/2)
  805.        if(nocomp(l).eq.nomip(icomp)) then
  806.        vvv=vpocha(jh,icomp)
  807.        goto  4555
  808.        endif
  809.                 enddo
  810.        vvv=XZERO
  811. 4555   continue
  812. *       write(6,7777) inupo,jh,poids(igau),vvv
  813. C
  814.         melval=ielval(icomp)
  815.         CPZA=vvv*poids(igau)*DJAC*AXIS
  816.          DO INO1 = 1,NBNN
  817.            RHS(ino1)=RHS(ino1)+CPZA*SHPXXX(1,ino1,igau)
  818. CCC      RHS(ino1)=RHS(ino1)+vvv*poids(igau)*
  819. CCC  1          DJAC*SHPXXX(1,INO1,IGAU)*AXIS
  820.          ENDDO
  821. 41       continue
  822. 42       continue
  823.  
  824. 7777  format(2I4,2X,2E12.5)
  825. C           SSOM1 =SSOM1+SURF*poids(igau)*vvv
  826. C
  827.          do i=1,NBNN
  828.          val(i)=XZERO
  829.          do j=1,NBNN
  830.          val(i)=val(i)+A(i,j)*RHS(J)
  831. *        write(6,*) ' A RHS ',a(i,j),rhs(j)
  832.          enddo
  833.          enddo
  834. C
  835.          do I=1,NBNN
  836.          velche(i,iel)=val(i)
  837.          enddo
  838. C
  839. *         if(icomp.lt.nomip(/2)) inupo =idebco
  840. C  ----  fin de la boucle sur les composantes
  841. 2220    continue
  842. C    fin de la boucle sur les elements
  843. 2210     continue
  844. C
  845.         do i=1,ielval(/1)
  846.           ipc1 = ielval(i)
  847.           call comred(ipc1)
  848.           ielval(i)=ipc1
  849.           melval = ipc1
  850.         enddo
  851.  
  852. C----------------  fin du traitement des massifs
  853.         segsup wrk5,wrk4,ipt5
  854.  
  855. C       segsup ipt3
  856. C    fin de la boucle sur les sous zones du modele
  857. 2100    continue
  858. C  destruction des chpo intermediaires
  859.         do i=1,ipcb(/1)
  860.           msoupo=ipcb(i)
  861. cgf          ipt5= igeoc (correction 7284)
  862.           mpoval= ipoval
  863. *          mchpoi=ipca(i)
  864.           segsup mpoval,msoupo
  865. cgf          segsup ipt5 (correction 7284)
  866.         enddo
  867.         segsup mchpoi
  868.         segsup ipcb,snomip
  869.         segsup icpr
  870.  7000   CONTINUE
  871. C (fdp) re-ajustement de MCOORD a sa taille initiale
  872.         NBPTS = NBNOE
  873.         SEGADJ MCOORD
  874. C       retrait des maillages temporaires du pre-conditionnement
  875. c       (leurs numero de noeuds depasse la taille de MCOOR)
  876.         call crech1b
  877. C (fdp) suppression du maillage temporaire IPT1
  878.         ipt1=ipgeom
  879.         segsup ipt1
  880.         segact ipt2
  881.         do io=1,ipt2.lisous(/1)
  882.           ipt5=ipt2.lisous(io)
  883. *         segsup ipt5
  884.         enddo
  885.         segsup ipt2
  886.         segsup ipt4
  887.         if(ipche.ne.0) then
  888.           mchelm=ipche
  889.         endif
  890.         IPOUT=MCHEL1
  891.         segsup sicoq
  892.  
  893.         return
  894.       END
  895.  
  896.  
  897.  
  898.  
  899.  
  900.  
  901.  
  902.  

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