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

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