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 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....  Copyright (c) 1984-1998: Robert L. Taylor & Juan C. Simo

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.

      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/


c     Set command word

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

      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
           cvalue = td
         elseif( pcomp(text(1),'twis',4)) then
           rate = td

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

      if( then

      if( then
          write(*,1002) coord,cvalue,rate
      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)

          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( then
              write(*,1003) i,y
          write(iow,1003) i,y

      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)