Télécharger la plaquette Cast3M


fr en

Télécharger impofu.eso

Retour à la liste

Numérotation des lignes : 

impofu

  1. C IMPOFU    SOURCE    PV090527  23/02/20    21:15:03     11602          
  2. *  reunion des relations portant sur le meme multiplicateur de lagrange
  3. *  si mchpoi est nul, on se contente de nettoyer les termes petits
  4. *
  5.       subroutine impofu(mrigid,mchpoi)
  6.       implicit real*8 (a-h,o-z)
  7. -INC SMRIGID
  8. -INC SMCHPOI
  9. -INC SMELEME
  10.  
  11. -INC PPARAM
  12. -INC CCOPTIO
  13. -INC SMCOORD
  14. -INC CCGEOME
  15. -INC CCREEL
  16. *  nombre max d'elements par paquet
  17.       parameter (nblim=16)
  18. *
  19.       segment icpr(nbpts)
  20.       segment val((nbpo-1)*idim,nblag)
  21.       segment ielem(nbpo,nblag)
  22.       segment nbnl(nblag)
  23.       segment dnorm(nblag)
  24.       segini icpr
  25. **    call prrigi(mrigid,0)
  26. **    call ecchpo(mchpoi,0)
  27.       segact mrigid
  28.       nblag=0
  29. *  regroupement des raideurs dans un grand tableau dimensionne au max
  30. * on commence par reperer les multiplicateurs de lagrange et leur destination
  31.       do 10 i=1,irigel(/2)
  32.        meleme=irigel(1,i)
  33.        segact meleme
  34.        do 20 iel=1,num(/2)
  35.         lag=num(1,iel)
  36.         if (icpr(lag).eq.0) then
  37.          nblag=nblag+1
  38.          icpr(lag)=nblag
  39. *  noter si formulation faible
  40.          if (icolor(iel).eq.1) icpr(lag)=-nblag
  41.         endif
  42.   20   continue
  43.   10  continue
  44. * sauvegarde d'un descr pour avoir les noms d'inconnues
  45.       des1=irigel(3,1)
  46.       segact des1
  47. * creation melme et xmatri a la taille max
  48.       nbpo=7
  49.       segini val,ielem,nbnl,dnorm
  50. * remplissage, d'abort le mult
  51.       do i=1,icpr(/1)
  52.        if (icpr(i).ne.0) then
  53.         ielem(1,abs(icpr(i)))=i
  54.        endif
  55.       enddo
  56. C
  57.       if (mchpoi.ne.0) then
  58.         segact mchpoi*mod
  59.         msoup1=ipchp(1)
  60.         segact msoup1
  61.         mpova1=msoup1.ipoval
  62.         segact mpova1
  63.       endif
  64. C
  65.       do 100 i=1,irigel(/2)
  66.        meleme=irigel(1,i)
  67.        xmatri=irigel(4,i)
  68.        segact xmatri,meleme
  69.        imod=0
  70.        do 110 iel=1,num(/2)
  71.         lag=num(1,iel)
  72.         kel=abs(icpr(lag))
  73.         remax=-1.
  74.         do iv=2,re(/1)
  75.           remax=max(remax,abs(re(iv,1,iel)))
  76.         enddo
  77. * recherche du point significatif le plus haut
  78.         isig=0
  79.         do 115 in=2,num(/1)
  80.          ip=num(in,iel)
  81. *  la contribution n'est pas significative on saute le noeud
  82.          do iv=1,idim
  83.            if (abs(re(1+(in-2)*idim+iv,1,iel)).gt.1d-10*remax) goto 116
  84.          enddo
  85. *  pour assurer que la somme des termes est nulle, on corrige le point isig
  86.          if (imod.eq.0) segact xmatri*mod
  87.          imod=1
  88.          isign=isig
  89.          if (isig.eq.0) isign=in+1
  90.          if (isign.gt.num(/1)) call erreur(5)
  91.          do iv=1,idim
  92.            re(1+(isign-2)*idim+iv,1,iel)=
  93.      >        re(1+(isign-2)*idim+iv,1,iel)+
  94.      >        re(1+(in-2)*idim+iv,1,iel)
  95.               re(1+(in-2)*idim+iv,1,iel)=0.d0
  96.          enddo
  97.          if (iimpi.ne.0) write (6,*) ' noeud elimine dans relation',
  98.      >         in,ip,re(1+(in-2)*idim+1 ,1,iel),remax
  99.          goto 115
  100.  116     continue
  101.          if (isig.eq.0) isig=in
  102.  115    continue
  103.         do 120 in=2,num(/1)
  104.          ip=num(in,iel)
  105. *  la contribution n'est pas significative on saute le noeud
  106.          do iv=1,idim
  107.            if (abs(re(1+(in-2)*idim+iv,1,iel)).gt.1d-10*remax) goto 160
  108.          enddo
  109.          goto 120
  110.  160     continue
  111.          do 130 ir=2,nbpo
  112.           if (ielem(ir,kel).eq.0) goto 150
  113.           if (ielem(ir,kel).ne.ip) goto 130
  114. *  le noeud est deja dans l'element, on ajoute les valeurs
  115.           do iv=1,idim
  116.            val((ir-2)*idim+iv,kel)=val((ir-2)*idim+iv,
  117.      >        kel)+re(1+(in-2)*idim+iv,1,iel)
  118.           enddo
  119.           goto 120
  120.  130     continue
  121. *  le noeud n'est pas dans l'element, on ajoute le noeud et les valeurs
  122.  150    continue
  123.         if (ir.gt.nbpo) then
  124.          nbpo=ir
  125.          segadj ielem,val
  126.         endif
  127.         ielem(ir,kel)=ip
  128.         do iv=1,idim
  129.          val((ir-2)*idim+iv,kel)=re(1+(in-2)*idim+iv,1,iel)
  130.         enddo
  131.  120    continue
  132. C
  133.         if (mchpoi.ne.0) then
  134.           dnorm(kel) = dnorm(kel) + mpova1.vpocha(iel,2)
  135.         endif
  136. C
  137.  110   continue
  138.       segsup meleme,xmatri
  139.       descr=irigel(3,i)
  140.       if (i.ne.1) segsup descr
  141.  100  continue
  142.       segsup mrigid
  143. *
  144. *  eclatement en paquets de meme nb de noeuds
  145. *  et renormalisation
  146. *
  147.  
  148. *  nb noeuds par element
  149.       do iel=1,ielem(/2)
  150.         do in=1,ielem(/1)
  151.          if (ielem(in,iel).eq.0) goto 200
  152.         enddo
  153.  200    continue
  154.         nbnl(iel)=in-1
  155.       enddo
  156. *  elimination elem en double si point ligne
  157.       if (idim.eq.3) then
  158.       do 300 iel=1,ielem(/2)
  159.        if (nbnl(iel).ne.4) goto 300
  160.         do 301 jel=iel+1,ielem(/2)
  161.         if (nbnl(jel).ne.4) goto 301
  162.         if (ielem(4,iel).ne.ielem(4,jel)) goto 301
  163.          if (ielem(2,iel)*ielem(3,iel).ne.ielem(2,jel)*ielem(3,jel)) 
  164.      >         goto 301
  165.          if (ielem(2,iel)+ielem(3,iel).ne.ielem(2,jel)+ielem(3,jel)) 
  166.      >         goto 301
  167.       if (iimpi.ne.0) write (6,*) 'elimination elem seg en double'
  168.         nbnl(jel)=0
  169.  301    continue
  170.  300  continue
  171.       endif
  172. *  elimination elem en double si point point
  173.       do 310 iel=1,ielem(/2)
  174.        if (nbnl(iel).ne.3) goto 310
  175.         do 311 jel=iel+1,ielem(/2)
  176.         if (nbnl(jel).ne.3) goto 311
  177.         if (ielem(3,iel).ne.ielem(3,jel)) goto 311
  178.         if (ielem(2,iel).ne.ielem(2,jel)) goto 311
  179.         if (iimpi.ne.0) write(6,*) 'elimination elem poin en double'
  180.         nbnl(jel)=0
  181.  311    continue
  182.  310  continue
  183.  
  184.       nrigel=0
  185.       segini mrigid
  186.       ichole = 0
  187.       imgeo1 = 0
  188.       imgeo2 = 0
  189.       isupeq = 0
  190.       iforig = ifour
  191.       mtymat='RIGIDITE'
  192. *
  193.       nbsous=0
  194.       nbref=0
  195.       do 250 iel=1,ielem(/2)
  196.        if (nbnl(iel).eq.0) goto 250
  197.        nbnn=nbnl(iel)
  198.        nbelem=0
  199.        do 255 jel=iel,ielem(/2)
  200.         if (nbnl(jel).eq.nbnn) nbelem=nbelem+1
  201.         if (nbelem.ge.nblim) goto 256
  202.  255   continue
  203.  256   continue
  204.        segini meleme
  205.        itypel=22
  206.        nligrd=(nbnn-1)*idim+1
  207.        nligrp=nligrd
  208.        nelrig=nbelem
  209.        segini xmatri
  210. *
  211.        segini descr
  212.         lisinc(1)=des1.lisinc(1)
  213.         lisdua(1)=des1.lisdua(1)
  214.         noelep(1)=1
  215.         noeled(1)=1
  216.         do inc=2,nligrd
  217.          lisinc(inc)=des1.lisinc(mod(inc-2,idim)+2)
  218.          lisdua(inc)=des1.lisdua(mod(inc-2,idim)+2)
  219.          noelep(inc)=(inc-2)/idim+2
  220.          noeled(inc)=(inc-2)/idim+2
  221.          enddo
  222.  
  223.        nrigel=nrigel+1
  224.        segadj mrigid
  225.        irigel(1,nrigel)=meleme
  226.        irigel(3,nrigel)=descr
  227.        irigel(4,nrigel)=xmatri
  228.        irigel(6,nrigel)=1
  229.        coerig(nrigel)=1.d0
  230.        kel=0
  231.        do 260 jel=iel,ielem(/2)
  232.         if (nbnl(jel).ne.nbnn) goto 260
  233.         kel=kel+1
  234.         do 265 in=1,nbnn
  235.          num(in,kel)=ielem(in,jel)
  236.  265    continue
  237.          if (icpr(num(1,kel)).lt.0) icolor(kel)=1
  238. ***     if (idim.eq.2.or.nbnn.eq.2) then
  239. ***      xnorm2=(abs(val(1,jel))+abs(val(3,jel)))**2+
  240. ***  >         (abs(val(2,jel))+abs(val(4,jel)))**2
  241. ***     else
  242. ***      xnorm2=(abs(val(1,jel))+abs(val(4,jel))+abs(val(7,jel)))**2+
  243. ***  >          (abs(val(2,jel))+abs(val(5,jel))+abs(val(8,jel)))**2+
  244. ***  >          (abs(val(3,jel))+abs(val(6,jel))+abs(val(9,jel)))**2
  245. ***     endif
  246. ***     xnorm=sqrt(xnorm2)
  247. ***     if (mchpoi.eq.0) xnorm=1.d0
  248.         xnorm=1.d0
  249.         if (mchpoi.ne.0) xnorm=dnorm(abs(icpr(num(1,kel))))
  250. C
  251.         do 270 in=1,idim*(nbnn-1)
  252.          re(in+1,1,kel)=val(in,jel)/xnorm
  253.          re(1,in+1,kel)=val(in,jel)/xnorm
  254.  270    continue
  255.         nbnl(jel)=0
  256.         if (kel.ge.nblim) goto 261
  257.  260   continue
  258.  261   continue
  259.  250   continue
  260.        segsup des1
  261. **    call prrigi(mrigid,0)
  262.  
  263. *  et maintenant le second membre
  264.        if (mchpoi.ne.0) then
  265. CCCC     segact mchpoi
  266. CCCC     msoup1=ipchp(1)
  267. CCCC     segact msoup1*mod
  268. CCC      nc=MSOUP1.NOCOMP(/2) - 1
  269.          nc=MSOUP1.NOCOMP(/2)
  270.          n=nblag
  271.          segini,msoupo,mpoval
  272.          nocomp(1)='FLX '
  273.          nocomp(2)='TAIL '
  274.          if (nc.eq.3) then
  275.            nocomp(3)='FADH'
  276.          endif
  277.          nbnn=1
  278.          nbelem=nblag
  279.          segini meleme
  280.          itypel=1
  281. CCCC     mpova1=msoup1.ipoval
  282. CCCC     segact mpova1
  283.          ipt1=msoup1.igeoc
  284.          segact ipt1
  285.          do i=1,ipt1.num(/2)
  286.            lag=ipt1.num(1,i)
  287.            num(1,abs(icpr(lag)))=lag
  288.            vpocha(abs(icpr(lag)),1)=mpova1.vpocha(i,1)/
  289.      >       dnorm(abs(icpr(lag)))+vpocha(abs(icpr(lag)),1)
  290. C
  291.            IF (IDIM.NE.3) THEN
  292.              vpocha(abs(icpr(lag)),2)=mpova1.vpocha(i,2) +
  293.      >        vpocha(abs(icpr(lag)),2)
  294.            ELSE
  295.              vpocha(abs(icpr(lag)),2)=sqrt(2.D0*mpova1.vpocha(i,2)) +
  296.      >        vpocha(abs(icpr(lag)),2)
  297.            ENDIF
  298. C
  299.            IF (nc.eq.3) THEN
  300.              vpocha(abs(icpr(lag)),2)=mpova1.vpocha(i,3)/
  301.      >         dnorm(abs(icpr(lag)))+vpocha(abs(icpr(lag)),2)
  302.            ENDIF 
  303.          enddo
  304.          segsup,msoup1,mpova1,ipt1
  305.          igeoc=meleme
  306.          ipoval=mpoval
  307.          ipchp(1)=msoupo
  308.        endif
  309. ***
  310. **    call ecchpo(mchpoi,0)
  311.        segsup val,ielem,nbnl,dnorm
  312.        return
  313.        end
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  

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