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

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