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( ). Also note that we only consider active nodes in the computation. One can determine if a node is active by looking at its NDTYP which is stored in the integer array located at np(190), thus ndtyp(1:numnp) => mr(np(190):np(190)+numnp-1). Active nodes have an ndtyp greater equal to zero, inactive nodes have negative values.

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
      
          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