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

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