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  26/04/30    21:15:41     12529          
  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.        RIGREL=0
  210.        segini xmatri
  211. *
  212.        segini descr
  213.         lisinc(1)=des1.lisinc(1)
  214.         lisdua(1)=des1.lisdua(1)
  215.         noelep(1)=1
  216.         noeled(1)=1
  217.         do inc=2,nligrd
  218.          lisinc(inc)=des1.lisinc(mod(inc-2,idim)+2)
  219.          lisdua(inc)=des1.lisdua(mod(inc-2,idim)+2)
  220.          noelep(inc)=(inc-2)/idim+2
  221.          noeled(inc)=(inc-2)/idim+2
  222.          enddo
  223.  
  224.        nrigel=nrigel+1
  225.        segadj mrigid
  226.        irigel(1,nrigel)=meleme
  227.        irigel(3,nrigel)=descr
  228.        irigel(4,nrigel)=xmatri
  229.        irigel(6,nrigel)=1
  230.        coerig(nrigel)=1.d0
  231.        kel=0
  232.        do 260 jel=iel,ielem(/2)
  233.         if (nbnl(jel).ne.nbnn) goto 260
  234.         kel=kel+1
  235.         do 265 in=1,nbnn
  236.          num(in,kel)=ielem(in,jel)
  237.  265    continue
  238.          if (icpr(num(1,kel)).lt.0) icolor(kel)=1
  239. ***     if (idim.eq.2.or.nbnn.eq.2) then
  240. ***      xnorm2=(abs(val(1,jel))+abs(val(3,jel)))**2+
  241. ***  >         (abs(val(2,jel))+abs(val(4,jel)))**2
  242. ***     else
  243. ***      xnorm2=(abs(val(1,jel))+abs(val(4,jel))+abs(val(7,jel)))**2+
  244. ***  >          (abs(val(2,jel))+abs(val(5,jel))+abs(val(8,jel)))**2+
  245. ***  >          (abs(val(3,jel))+abs(val(6,jel))+abs(val(9,jel)))**2
  246. ***     endif
  247. ***     xnorm=sqrt(xnorm2)
  248. ***     if (mchpoi.eq.0) xnorm=1.d0
  249.         xnorm=1.d0
  250.         if (mchpoi.ne.0) xnorm=dnorm(abs(icpr(num(1,kel))))
  251. C
  252.         do 270 in=1,idim*(nbnn-1)
  253.          re(in+1,1,kel)=val(in,jel)/xnorm
  254.          re(1,in+1,kel)=val(in,jel)/xnorm
  255.  270    continue
  256.         nbnl(jel)=0
  257.         if (kel.ge.nblim) goto 261
  258.  260   continue
  259.  261   continue
  260.  250   continue
  261.        segsup des1
  262. **    call prrigi(mrigid,0)
  263.  
  264. *  et maintenant le second membre
  265.        if (mchpoi.ne.0) then
  266.          nc=MSOUP1.NOCOMP(/2)
  267.          n=nblag
  268.          segini,msoupo,mpoval
  269.          nocomp(1)='FLX '
  270.          nocomp(2)='TAIL '
  271.          if (nc.eq.3) then
  272.            nocomp(3)='FADH'
  273.          endif
  274.          nbnn=1
  275.          nbelem=nblag
  276.          segini meleme
  277.          itypel=1
  278. CCCC     mpova1=msoup1.ipoval
  279. CCCC     segact mpova1
  280.          ipt1=msoup1.igeoc
  281.          segact ipt1
  282.          do i=1,ipt1.num(/2)
  283.            lag=ipt1.num(1,i)
  284.            kel=abs(icpr(lag))
  285. C
  286.            num(1,kel)=lag
  287.            vpocha(kel,1)=mpova1.vpocha(i,1)/dnorm(kel)+vpocha(kel,1)
  288. C
  289.            if (idim.ne.3) then
  290.              vpocha(kel,2)=mpova1.vpocha(i,2)+vpocha(kel,2)
  291.            else
  292.              vpocha(kel,2)=sqrt(2.D0*mpova1.vpocha(i,2))+vpocha(kel,2)
  293.            endif
  294. C
  295.            if (nc.eq.3) then
  296.              vpocha(kel,3)=mpova1.vpocha(i,3)/dnorm(kel)+vpocha(kel,3)
  297.            endif 
  298.          enddo
  299.          segsup,msoup1,mpova1,ipt1
  300.          igeoc=meleme
  301.          ipoval=mpoval
  302.          ipchp(1)=msoupo
  303.        endif
  304. ***
  305. **    call ecchpo(mchpoi,0)
  306.        segsup val,ielem,nbnl,dnorm
  307.        return
  308.        end
  309.  
  310.  
  311.  
  312.  

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