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

Last login: Tue Jun 24 13:42:13 on ttys005

The default interactive shell is now zsh.
To update your account to use zsh, please run `chsh -s /bin/zsh`.
For more details, please visit https://support.apple.com/kb/HT208050.
edr07170:~ sg$ cd FEap
edr07170:FEap sg$ cd ver87
edr07170:ver87 sg$ cd program/
edr07170:program sg$ vi pmacr.f
edr07170:program sg$ cd
edr07170:~ sg$ cd Collaborators/Ali/ElementsandMacros2023-05-08/ElementsandMacrosintel/
edr07170:ElementsandMacrosintel sg$ ls
AAvlce.f		elmt05.o		umacr20.o
AAvlce.o		elmt07.f		umacr49.f
Basics.f90		elmt07.o		umacr49.o
Basics.o		elmt17.f		umacr50.f
CommonTerms.f90		elmt17.o		umacr50.o
CommonTerms.o		feap			umesh9.f
DNRvlce.f		feap.dSYM		umesh9.o
DNRvlce.o		fevlcealt.f		usolve.f
FreeEnergyModels.f90	fevlcealt.o		vlcemod1.f
FreeEnergyModels.o	freeenergymodels.mod	vlcemod1.o
NRres.f			inptvlce.f		vlcemod2.f
NRres.o			inptvlce.o		vlcemod2.o
NRtan.f			kinevlce.f		vlcemod3.f
NRtan.o			kinevlce.o		vlcemod3.o
NRvlce.f		makefile		vlcemod4.f
NRvlce.o		p_paraview.f		vlcemod4.o
PKvlce.f		p_paraview.o		vlcemod5.f
PKvlce.o		pmacr2.f		vlcemod5.o
arclen.f		pmacr2.o		vlcemod6.f90
arclen.o		slcnvlce.f		vlcemod6.o
basics.mod		slcnvlce.o		volvlce.f
commonterms.mod		umacr13.f		volvlce.o
david.zip		umacr13.o		xpardiso.h
elmt05.f		umacr20.f
edr07170:ElementsandMacrosintel sg$ vi elmt
elmt05.f  elmt05.o  elmt07.f  elmt07.o  elmt17.f  elmt17.o  
edr07170:ElementsandMacrosintel sg$ vi elmt17.f
edr07170:ElementsandMacrosintel sg$ vi umesh9.f
edr07170:ElementsandMacrosintel sg$ vi elmt17.f

          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