UMESH1.F - - Example

umesh[0-9].f files in the user directory allow one to create user mesh commands. To set up a user mesh command one needs to edit the file umesh1.f (for example) to include the code one wishes to have executed when the user mesh command is called from inside FEAP. The interface can be programed in fortran or it can be programed in C. In the example that follows a fortran interface is utilzed.

The main steps in creating a user mesh command whether in C or Fortran are

Set the name of the mesh command to whatever 4-letter name you wish it to have when FEAP calls it.
Set up the source code you wish to have executed when the command is called from FEAP

The mesh command name is set in the common variable uct in the common block umac1. When uct is equal to 'mes1' (for umesh1.f or .c) it should be set to whatever you wish for the name of your mesh command.

A example umesh1.f is given below. The name of the user macro is set to tors. The mesh command itself is used for applying a rotation to a set of boundary displacements to the nodes on a given x, y, or z plane in accordance with a rotation about either the x, y, or z axis. The input is done in a polling fashion and looks for the key words height <x,y,z> value to define the x, y, or z coordinate value of the plane of interest. It then reads a twist rate value for applying the rotation.

Source file:

c$Id: umesh1.f,v 1.1 1998/08/28 18:19:29 rlt Exp $
subroutine umesh1(prt)

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

c-----[--.----+----.----+----.-----------------------------------------]
c Purpose: Dummy user input routine

c Inputs:
c prt - Flag, output results if true

c      Outputs:
c         none   - Users are responsible for generating outputs
c                  through common blocks, etc. See programmer
c                  manual for example.
c-----[--.----+----.----+----.-----------------------------------------]

implicit none

      real*8       td,cvalue,rate,tol
      real*8       x(3),c,s,y(3)
      integer      coord,i,j,k
      character    text(2)*15
      logical      pcomp,prt,tinput,errck

      include 'cdata.h'
      include 'iofile.h'
      include 'pointer.h'
      include 'sdata.h'
      include 'umac1.h'

include 'comblk.h'

data tol /1.d-5/

save

c Set command word

      if(pcomp(uct,'mes1',4)) then
        uct = 'tors' ! Set the name you wish for the mesh command and return
        return
      endif

errck = tinput(text,2,td,1)

do while (.not.pcomp(' ',text(1),4)) ! Polling input loop for plane height and twist rate

         if( pcomp(text(1),'heig',4)) then
           if(     pcomp(text(2),'x',1)) then
               coord = 1
           elseif( pcomp(text(2),'y',1)) then
               coord = 2
           elseif( pcomp(text(2),'z',1)) then
               coord = 3
           endif
           cvalue = td
         elseif( pcomp(text(1),'twis',4)) then
           rate = td
         endif

errck = tinput(text,2,td,1)
end do

      if(ior.lt.0) then
          write(*,1001)
      endif
      write(iow,1001)

      if(ior.lt.0) then
          write(*,1002) coord,cvalue,rate
      endif
      write(iow,1002) coord,cvalue,rate

c = cos(cvalue*rate)
s = sin(cvalue*rate)

do i = 1, numnp ! loop over nodes to find those on the plane of interest

         j = np(43)+(i-1)*3 ! np(43) is the blank common pointer to the coordinates
        x(1) = hr(j)
        x(2) = hr(j+1)
        x(3) = hr(j+2)

if(abs(x(coord)-cvalue).lt.tol*cvalue) then ! only find those nodes near the plane of interest

          if(coord.eq.1) then   ! Determine the nodal displacements desired
            y(1) = 0.d0
            y(2) = c*x(2) - s*x(3) - x(2)
            y(3) = s*x(2) + c*x(3) - x(3)
          elseif(coord.eq.2) then
            y(2) = 0.d0
            y(3) = c*x(3) - s*x(1) - x(3)
            y(1) = s*x(3) + c*x(1) - x(1)
          elseif(coord.eq.3) then
            y(3) = 0.d0
            y(1) = c*x(1) - s*x(2) - x(1)
            y(2) = s*x(1) + c*x(2) - x(2)
          endif

          k = np(27)+ndf*numnp+(i-1)*3 ! Set the displacement values in np(27)+ndf*numnp
          hr(k) = y(1)
          hr(k+1) = y(2)
          hr(k+2) = y(3)

          if(ior.lt.0) then
              write(*,1003) i,y
          endif
          write(iow,1003) i,y

endif
end do

1001 format(/5x,'Torsional Boundary Condition'/)
1002 format(5x,'Longitudinal axis = ',i5/
& 5x,'Application height = ',f8.5/
& 5x,'Twist Rate = ',f8.5//)
1003 format(5x,'Node ',i5,'; Disp ',3f8.5)

end