Sources Esope | Cast3M

dbblx
C DBBLX     SOURCE    PV        21/04/26    21:15:08     10866          
*  dedouble les multiplicateurs de Lagrange
*
      subroutine dbblx(mrigid,lagdua)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMRIGID
-INC SMELEME
*     write(6,*) ' dans dbblx '
      idimp1 = idim + 1
 
      segact mcoord*mod
      nbini = nbpts
 
      nbsous = 0
      nbref  = 0
      nbnn   = 2
      nbelem = 0
      segini,ipt8
      ipt8.itypel = 2
      nbele8 = nbelem
 
      segact mrigid*mod
      nrigel = irigel(/2)
 
C     Boucle d'activation et AJUSTEMENT MCOORD en 1 coup
      do 9 ir = 1, nrigel
        ipt1 = irigel(1,ir)
        segact ipt1
        i_z = ipt1.itypel
        if (i_z .ne. 22) goto 9
        nbpts  = nbpts + ipt1.num(/2)
 9    continue
      segadj,mcoord
 
      nbpts=nbini
      do 10 ir = 1, nrigel
        ipt1 = irigel(1,ir)
        segact ipt1
        i_z = ipt1.itypel
*       write(6,*) ' dbblx itypel ',ipt1.itypel
        if (i_z .ne. 22) goto 10
*       write(6,*) ' dbblx dedoublement ',ir
        segini,meleme=ipt1
*       write(6,*) ' dbblx ipt1 meleme ',ipt1,meleme
        itypel=49
        nbsous = 0
        nbref  = 0
        nbnn   = meleme.num(/1)+1
        nbelem = meleme.num(/2)
        segadj meleme
        nbsup  = nbelem
 
        nbpts0 = nbpts
        nbpts  = nbpts0 + nbsup
 
        do 100 j = 1, nbsup
          do 120 i = nbnn,3,-1
            meleme.num(i,j) = meleme.num(i-1,j)
 120      continue
          meleme.num(2,j) = nbpts0 + j
          ip2 = (meleme.num(2,j)-1) * idimp1
          ip1 = (meleme.num(1,j)-1) * idimp1
          do k = 1, idimp1
            xcoor(ip2+k) = xcoor(ip1+k)
          enddo
 100    continue
        irigel(1,ir) = meleme
**  on garde la liste des noeuds rajoutés
*        nbsous = 0
*        nbref  = 0
        nbnn   = 2
        nbelem = nbele8 + nbsup
        segadj ipt8
        do 130 j = 1, nbsup
          j8 = nbele8 + j
          ipt8.num(1,j8) = meleme.num(1,j)
          ipt8.num(2,j8) = meleme.num(2,j)
 130    continue
        nbele8 = nbelem
*
        des1 = irigel(3,ir)
        segini,descr=des1
        segdes,des1
        nligrp = lisinc(/2)+1
        nligrd = lisdua(/2)+1
        segadj descr
        do 200 i = nligrp, 3, -1
          lisinc(i) = lisinc(i-1)
          noelep(i) = noelep(i-1)+1
 200    continue
        lisinc(2) = 'LX'
        noelep(2) = 2
        do 210 i = nligrd, 3, -1
          lisdua(i) = lisdua(i-1)
          noeled(i) = noeled(i-1)+1
 210    continue
        lisdua(2)='FLX'
        noeled(2)=2
        segdes,descr
        irigel(3,ir) = descr
 
        xmatr1 = irigel(4,ir)
        segini,xmatri=xmatr1
        segdes,xmatr1
        nelrig = re(/3)
        segadj,xmatri
        do 300 im = 1, nelrig
          do     i = nligrp, 1, -1
            do     j = nligrd, 2, -1
*              re(i,j,im) = re(i,j-1,im)
              re(j,i,im) = re(j-1,i,im)
            enddo
          enddo
          do     j = nligrd, 1, -1
            do     i = nligrp, 2, -1
*              re(i,j,im) = re(i-1,j,im)
              re(j,i,im) = re(j,i-1,im)
            enddo
          enddo
* normaliser les nouveaux termes par rapport au max de la relation
          xnorm=0.D0
          do i1=1,re(/1)
          do i2=1,re(/2)
           xnorm=max(abs(re(i1,i2,im)),xnorm)
          enddo
          enddo
          re(1,1,im) = -xnorm+re(1,1,im)
          re(1,2,im) = +xnorm
          re(2,1,im) = +xnorm 
          re(2,2,im) = re(1,1,im)
 300    continue
        segdes,xmatri
        irigel(4,ir) = xmatri
 
  10  continue
*
*
      if (nbele8.eq.0) then
        segsup,ipt8
        ipt8 = 0
      endif
 
      lagdua = ipt8
 
      end