Element Lagrange Multipliers

Consider wanting to move a node radially in the $x$-$y$ plane by an amount given by a proportional load value. This can be accomplished by enforcing the constraint $(X+u_{x})^{2}+(Y+u_{y})^{2}-\mathrm {prop_{value}} ^{2}=0$. To accomplish this we can create an element with one node and implement the constraint in the element.

Note that when using a mesh with this type of element in it, one needs to let FEAP know in the input file that there will be element Lagrange multipliers in the problem by setting the value of nad on the control record. Since we only have one constraint, there will be only one Lagrange multiplier per point element. Supposing we apply the point constraint element in a 3D mesh with 8-node bricks, then the control record will look like:

feap ** Point constraint example **
0 0 0 3 3 8 1

Or perhaps more clearly as

feap ** Point constraint example **
ndm = 3
ndf = 3
nen = 8
nad = 1

The crucial points in the use of element Lagrange multipliers is to note that what is to be programed as the internal force is $\delta \lambda \,c(u)+\lambda \partial c/\partial u\cdot \delta u$. It is important to note that the location of the equations for the Lagrange multiplier dof are located after the nen x ndf. See the programmer manual section on adding elements on the Manuals page.

Crtical points to note:

• Under isw=1, you must declare that the element contribute nlm Lagrange multiplier equations. The common block hdata.h is needed for this purpose.
• The value of the multiplier will be contained in ule(1,1); if there were two constraints in the element the multiplier for the second constraint would be found in ule(1,2). The leading dimension of ule( , ) takes values 1:3. Slot 2 provides the increment to the multiplier from the prior time step and Slot 3 provides the increment within any current Newton iteration.
• The example shows how to set up element time history plots by including the common block elplot.h.

!$Id:$
      subroutine elmt17(d,ul,xl,ix,tl,s,p,ndf,ndm,nst,isw)

!      * * F E A P * * A Finite Element Analysis Program

!....  Copyright (c) 1984-2023: Regents of the University of California
!                               All rights reserved

!-----[--.----+----.----+----.-----------------------------------------]
!     Modification log                                Date (dd/mm/year)
!       Original version                                    01/11/2006
!    1. Added options to output 'stresses' and tplots       24/05/2023
!-----[--.----+----.----+----.-----------------------------------------]
!     Purpose: Radial point contraint element
!              (X+ux)^2 + (Y+uy)^2 - [prop(k)]^2 == 0

!     Inputs:
!       d(*)          - Material set parameters
!       ul(ndf,nen,*) - Solution parameters
!       xl(ndm,nen)   - Element nodal coordinates
!       ix(nen1)      - Element global node numbers
!       tl(nen)       - Element vector (e.g., nodal temperatures)
!       ndf           - Maximum no dof's/node
!       ndm           - Mesh spatial dimension
!       nst           - Element matrix/vector dimension
!       isw           - Switch parameter (set by input commands)

!     Outputs:
!       s(nst,nst,2)  - Element matrix (stiffness, mass, etc.)
!       p(nst,2)      - Element vector (residual, lump mass, etc.)
!-----[--.----+----.----+----.-----------------------------------------]
      implicit none

      include   'bdata.h'               ! head
      include   'cdata.h'               ! nen
      include   'iofile.h'              ! iow
      include   'eldata.h'              ! pstype,mct,n_el,ma
      include   'elplot.h'              ! tt
      include   'hdata.h'               ! nlm
      include   'lmdata.h'              ! ule(*,*)
      include   'prld1.h'               ! prldv(*)
      include   'umac1.h'               ! utx(*)

      integer       ::  ndf             ! Max DOF's/node
      integer       ::  ndm             ! Mesh spatial dimention
      integer       ::  nst             ! Matrix 's' dimension
      integer       ::  isw             ! Switch for solution steps
      integer       ::  ix(*)           ! Element global nodes
      real (kind=8) ::  d(*)            ! Material set parameters
      real (kind=8) ::  ul(ndf,*)       ! Local nodal solutions
      real (kind=8) ::  xl(ndm,*)       ! Local nodal coordinates
      real (kind=8) ::  tl(*)           ! Local nodal load array
      real (kind=8) ::  s(nst,*)        ! Element matrix
      real (kind=8) ::  p(ndf,*)        ! Element vector

      ! Element variables
      logical            :: pcomp, errck, tinput
      character (len=15) :: tdata(1)
      real (kind=8)      :: td(1),ux,uy,la,rr,c
      integer            :: nla

      if(isw.lt.0) then

        utx(1) = 'rpcon'     ! Set element name
                                          
      elseif(isw.eq.0) then    
        if(ior.lt.0) then
          write(*,2000)       
        endif

      elseif(isw.eq.1) then                ! Input material set data

        nlm = 1                            ! Set the additional equation for the LM

        tdata(1) = 'start'
        do while (.not.pcomp(tdata(1),'    ',4))
          errck = tinput(tdata,1,td(1),1)
          if(pcomp(tdata(1),'prop',4)) then
            d(1) = td(1)                   ! Proptable to look in for r value
          endif
        end do
        write(iow,2001) nint(d(1))

      elseif(isw.eq.2) then                ! Check input data
      elseif(isw.eq.3 .or.isw.eq.6) then   ! Compute residual/tangent

        ux = ul(1,1)+xl(1,1)               ! Current position x
        uy = ul(2,1)+xl(2,1)               ! Current position y
        la = ule(1,1)                      ! Lagrange multiplier
        rr = prldv(nint(d(1)))             ! Prop load given radius
        c  = ux*ux+uy*uy-rr*rr             ! constraint

        nla = ndf*nen + 1                  ! row/col for constraint equation

        ! Residual
        p(1,1) = -2.d0*ux*la
        p(2,1) = -2.d0*uy*la
        p(nla,1) = -c

        ! Stiffness
        if(isw.eq.3) then
          s(1,1) = 2.d0*la
          s(1,nla) = 2.d0*ux 
          s(2,2) = 2.d0*la
          s(2,nla) = 2.d0*uy 
          s(nla,1) = 2.d0*ux
          s(nla,2) = 2.d0*uy
        endif

        ! Tplot options
        tt(1) = -2*(xl(1,1)+ul(1,1))*ule(1,1)   ! Fx
        tt(2) = -2*(xl(2,1)+ul(2,1))*ule(1,1)   ! Fy
        tt(3) = prldv(nint(d(1)))               ! radius

      elseif(isw.eq.4 .or.isw.eq.8) then   ! Output/plot element data
        if(isw.eq.4) then

          mct = mct - 2
          if(mct.lt.0) then
            write(iow,3000) head
            if(ior.lt.0) then
              write(iow,3000) head
            endif
            mct = 50
          endif

          write(iow,3010) n_el,ma,xl(1:3,1),ul(1:2,1),ule(1,1),
     &                 -2*(xl(1,1)+ul(1,1))*ule(1,1),
     &                 -2*(xl(2,1)+ul(2,1))*ule(1,1),
     &                 -2*ule(1,1)*prldv(nint(d(1)))
          if(ior.lt.0) then
            write(*,3000) head
            write(*,3010) n_el,ma,xl(1:3,1),ul(1:2,1),ule(1,1),
     &                 -2*(xl(1,1)+ul(1,1))*ule(1,1),
     &                 -2*(xl(2,1)+ul(2,1))*ule(1,1),
     &                 -2*ule(1,1)*prldv(nint(d(1)))
          endif
        endif

      elseif(isw.eq.5) then                ! Compute mass matrix
      elseif(isw.eq.14) then               ! Set non-zero history values
      endif

!     Formats

2000  format('    Elmt 17: Radial motion constraint.')
2001  format(5x,'Radial Constraint:'//
     &      10x,'Proportional Load Number =',1p,1i3/1x)
3000  format(1x,19a4,a3//5x,'Point Constraint'//
     &      10x,'Elmt Mat   1-Coord     2-Coord     3-Coord'/
     &      21x,'1-Displ     2-Displ     LagrangeM'/
     &      21x,'1-Force     2-Force     Mag.Force')
3010  format(i14,i4,1p,3e12.4/(18x,1p,3e12.4))
      end subroutine elmt17