Element Lagrange Multipliers
Jump to navigation
Jump to search
Consider wanting to move a node radially in the - plane by an amount given by a proportional load value. This can be accomplished by enforcing the constraint . 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 .
!$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