User macro to flip the boundary codes
Consider the problem of mimicking the behavior of a standard load/unload experiment that is displacement controlled up to the point where the load cell reads zero force and the machine then switches to load control. In a FEAP computation the boundary codes on the driven nodes can be dynamically changed by changing the values in the ID(ndf,numnp,1:2) array. In particular if for a given dof and node number the value of ID( , , 2) is set to zero then FEAP will treat the dof as active and look for forced values as set by the FORCe command in the input file. If these have been omitted, then the force will be taken as zero.
Here is an example of a MACRo command that will achieve this purpose. It allows one to identify the nodes to be switched using a simple search on a given plane. Then it allows one to switch the boundary codes once the reactions on the selected nodes drop below a specified tolerance for times greater than some given value. Note that after the boundary codes have been changed, it necessary to recompute the equation profile as there now additional active equations in the problem; this can be done with a call to pnewpr( ).
Here is an input file that will use the MACRo called ULNM
feap ** UNLM user macro test, viscoelastic bar unloading ** ndf = 2 ndm = 2 nen = 4 mate 1 ! Viscoelastic material solid elastic isotropic 10 0.3 viscoelastic deviatoric 0.5 1.0 block ! Strip geometry cart 5 10 1 0.0 0.0 2 5.0 0.0 3 5.0 10.0 4 0.0 10.0 eboun ! Clamped top and bottom 2 0 1 1 2 10 1 1 edisp ! Driven displacement at the top, No defined forces, they will default to zero 2 10 0 1.0 end batch ! prop curve for the displacements prop,,1 end 2 1 0.0 0.0 1.0 1.0 2.0 0.0 batch ! Get time history of the diplacement at the top-left node, and top surface reactions tplot end disp,,2,0.0,10.0 rsum,2,61,66 show batch plot,defo,1 dt,,.01 loop,,100 ! Extend to displacement of 1.0 time loop,,10 tang,,1 next plot,wipe plot,mesh plot,boun next UNLM,edir,2,10.0 ! Set-up nodes to be released loop,,500 time loop,,10 tang,,1 next plot,wipe plot,mesh plot,boun reac ! compute reactions which are needed for the UNLM macro UNLM,forc,0.01,1.5 ! Check if cross-head reaction force is below 0.01 and that time is above 1.5, then only flip next end stop
After the computation one can plot the time histories that have been saved to see that when the reaction force dropped below 0.01, the boundary conditions where flipped.
The source code for the MARCo UNLM is given as follows. Note that the logic here is quite simple and one could create a much better control algorithm to find more carefully on the zero crossing, but this example is really meant to illustrate the use of FEAP data structures and provide a non-trivial programming example. As a side note, ID( , , 1) contains the equations numbers (necessarily greater than zero for active equations).
!$Id:$
subroutine umacr49(lct,ctl)
! * * F E A P * * A Finite Element Analysis Program
!.... Copyright (c) 1984-2022: Regents of the University of California
! All rights reserved
!-----[--.----+----.----+----.-----------------------------------------]
! Modification log Date (dd/mm/year)
! Original version 01/11/2006
! 1. Remove 'prt' from argument list 09/07/2009
! 2. Add 'uhelpfl' comment option 19/10/2017
! 3. Add 'help' comment option - replaces 'uhelpfl' 07/01/2019
!-----[--.----+----.----+----.-----------------------------------------]
! Purpose: Unload to set level and flip bcs to zero load
! 1. UNLM EDIR VAL TOL
! 2. In time loop: REAC; UNLM FVAL TVAL (command applies
! once REAC < FVAL and T > TVAL
!
! Inputs:
! lct - Command character parameters
! ctl(3) - Command numerical parameters
! Outputs:
! N.B. Users are responsible for command actions. See
! programmers manual for example.
!-----[--.----+----.----+----.-----------------------------------------]
implicit none
include 'cdata.h'
include 'comfil.h'
include 'endata.h'
include 'iofile.h'
include 'pdata0.h'
include 'pointer.h'
include 'sdata.h'
include 'tdata.h'
include 'umac1.h'
include 'comblk.h'
character (len=15) :: lct
logical :: pcomp
real (kind=8) :: ctl(3)
! umacro variables
logical :: flipflg,pass1
integer :: edir, cnt, i, n
integer, allocatable :: bn(:)
integer, pointer :: id(:,:,:)
real (kind=8) :: xval, fval, tval, xtol, cforce
real (kind=8), pointer :: x(:,:), reac(:,:)
save
! Set command word
if(pcomp(uct,'ma49',4)) then ! Usual form
uct = 'unlm' ! Specify 'name'
flipflg = .false. ! BCs have not yet been flipped
pass1 = .true. ! First pass through the macro
elseif(pcomp(lct,'help',4)) then ! Write help information
write(*,*) 'COMMAND: UNLM unload macro'
elseif(urest.eq.1) then ! Read restart data
elseif(urest.eq.2) then ! Write restart data
else ! Perform user operation
! Pointers into global arrays
x(1:ndm,1:numnp) => hr(np(43):np(43)+ndm*numnp-1)
reac(1:ndf,1:numnp) => hr(np(26):np(26)+nneq-1)
id(1:ndf,1:numnp,1:2) => mr(np(31):np(31)+2*nneq-1)
if(pcomp(lct,'edir',4)) then
! Set up edge to look at
pass1 = .false.
edir = nint(ctl(1))
xval = ctl(2)
xtol = ctl(3)
! Set tolerance if not set
if(xtol.eq.0.0d0) then
xtol = 1.0d-6*dxmsh(edir)
endif
write(iow,1000) edir,xval,xtol
! Count the number of nodes on edge
cnt = 0
do n = 1,numnp
if(mr(np(190)-1+n).ge.0 .and. abs(x(edir,n)-xval).le.xtol)
& then
cnt = cnt + 1
write(iow,*) n
endif
end do ! n
! Store the node numbers
allocate(bn(cnt))
cnt = 0
do n = 1,numnp
if(mr(np(190)-1+n).ge.0 .and. abs(x(edir,n)-xval).le.xtol)
& then
cnt = cnt + 1
bn(cnt) = n
endif
end do ! n
write(iow,1003)
write(iow,1004) (bn(i),i=1,cnt)
elseif(pcomp(lct,'forc',4)) then
if(pass1) write(*,*) 'WARNING: Call UNLM with EDIR first'
! Force and Time values for the switch
fval = ctl(1) ! force value
tval = ctl(2) ! time value
if(.not.rfl) write(*,*) 'Call REAC before UNLM'
! Get edge reaction summed
cforce = 0.d0
do n = 1,cnt
cforce = cforce - reac(edir,bn(n))
end do ! n
rfl = .false.
if(ior.lt.0) write(*,*) ttim,cforce,flipflg
write(iow,*) ttim,cforce,flipflg
if(.not.flipflg) then ! Only check if BCs not yet flipped
if(ttim.ge.tval .and. cforce.lt.fval) then
flipflg = .true.
! Flip Boundary Codes
do n = 1,cnt
id(edir,bn(n),2) = 0
end do
write(iow,*) 'No. nodes released ',cnt
if(ior.lt.0) write(*,*) 'No. nodes released ',cnt
call pnewpr()
endif
endif
endif
endif
1000 format(/'Unloading command data'/,' Direction ',i4/,
& ' Coor. val ',e12.4/,' Tol (if def) ',e12.4)
1003 format(/'Identified nodes by UNLM'/)
1004 format(i4)
end subroutine umacr49