subroutine forcem(press,xp1,xp2,mnn,mn)
c
c CONTACT: ap@marc.de
c      
      implicit real*8 (a-h,o-z)                                               dp
      dimension mn(7),xp1(3),xp2(3)
c
c PURPOSE
c  Apply an angular accelleration to an element.
c
c USAGE
c  -define rotation axis with ROTATION A model definition option
c  -specify the DIST LOADS type for magnitude and arbitrary direction
c  -specify the desired angular accelleration as the DIST LOADS 
c   pressure value
c
c  To define the model with mentat one can proceed as follows:
c  -Define a centrifugal load as a boundary condition.
c   Here one also defines the rotation axis.
c  -Write the MARC input file when all other definitions are done.
c  -Edit the MARC input file and modify the DIST LOADS block.
c   Change the distribited load type 100 to the appropriate
c   type which flags a nonuniform load in arbitrary direction.
c   Change the value of the distributed load to the desired
c   angular accelleration.
c
c PROCEDURE
c  The subroutine will compute the accelleration components at
c  each integration point which result from the specified angular
c  accelleration about the specified rotation axis. The magnitude
c  of this accelleration vector is simply angular accelleration
c  times distance to the rotation axis.
c  This acceleration vector will be scaled with the mass density
c  at the integration point.
c  For shell and beam elements the vector will also be scaled with
c  shell thickness or beam cross-sectional area.
c  The subroutine returns the magnitude of this scaled accelleration
c  vector and its normalized components.
c
c* * * * * *
c
c INPUT
c  press    distributed load magnitude given in input
c  xp1      integration point coordinates
c  mnn      integration point number
c  mn       element number,etc
c
c OUTPUT
c  press    magnitude of accelleration vector
c  xp2      normalized accelleration vector components
c
c* * * * * *
c
c     x1     point on rotation axis
c     xl     direction cosines of rotation axis
c     xp     direction cosines of vector from x1 to xp1
c
c* * * * * *
c
      dimension x1(3),xl(3),xp(3)
c
      include '../common/matdat'
      include '../common/elmcom'
      include '../common/heat'
      include '../common/space'
      include '../common/arrays'
      include '../common/lass'
c
      logical debug
      data debug/.true./
c
c... data base offset (here use internal element number)
      lofr = (n-1)*nelstr
c
c check for different element classes
      if (lclass.eq.8) then
c volume elements
       area = 1d0
      else if (lclass.eq.2) then
c shell elements (get thickness)
       la1  = ithick + lofr + nn - 1
       area = vars(la1)
      else if (lclass.eq.1) then
c truss & elastic beam elements (get cross sectional area)
       la1  = igeomr + lofr
       area = vars(la1)
      else
c remaining classes are not accounted for
       area = 1d0
       write(6,1000) mn(1),lclass
      end if
      if (debug) write(0,'(a20,2i14)') 'ELEMENT: ',mn(1),mnn
      if (debug) write(0,'(a20,i14,2e14.5)') 'lclass: ',lclass,area,rho
c
 1000 format('Element',i7,': this element class is not accounted for
     * (class',i3,')')
c
c define angular accelleration (value from input)
      ang_acc = press
c
c get point on the rotation axis
      x1(1) = vars(icentr)
      x1(2) = vars(icentr+1)
      x1(3) = vars(icentr+2)
      if (debug) write(0,'(a20,3e14.5)') 'centre: ',(x1(ijk),ijk=1,3)
c
c get direction cosines of rotation axis
      xl(1) = vars(icenta)
      xl(2) = vars(icenta+1)
      xl(3) = vars(icenta+2)
      if (debug) write(0,'(a20,3e14.5)') 'direction cos: ',
     *           (xl(ijk),ijk=1,3)
c
c compute direction vector from x1 to integration point xp1
      xp(1) = xp1(1) - x1(1)
      xp(2) = xp1(2) - x1(2)
      xp(3) = xp1(3) - x1(3)
c
c accelleration direction is vector product xl X xp
      call vecprd(xl,xp,xp2)
      call x2norm(xp2,3,rp)
      do i1=1,3
       xp2(i1) = xp2(i1)/rp
      enddo
      if (debug) write(0,'(a20,3e14.5)') 'acc dir cos:',
     *           (xp2(ijk),ijk=1,3)
c
c magnitude is product of angular accelleration and radius
      press = ang_acc * rp * rho * area
      if (debug) write(0,'(a20,2e14.5)') 'ang_acc, radius: ',ang_acc,rp
c
      return
      end
      subroutine vecprd(x,y,z)
c
c* * * * *
c
c PURPOSE:
c  compute vector out product z = x X y
c
c INPUT:
c  x,y      input vectors
c OUTPUT
c  z        output vector
c
c* * * * *
c
      implicit double precision (a-h,o-z)
c
      dimension x(3),y(3),z(3)
c
      z(1) = x(2)*y(3) - x(3)*y(2)
      z(2) = x(3)*y(1) - x(1)*y(3)
      z(3) = x(1)*y(2) - x(2)*y(1)
c
      return
      end
      subroutine x2norm(x,n,a)
c
c* * * * *
c
c PURPOSE
c  compute vector 2-norm
c
c INPUT:
c  x       input vector
c  n       vector dimension
c
c OUTPUT
c  a       2-norm
c
c* * * * *
c
      implicit double precision (a-h,o-z)
c
      dimension x(1)
c
      a = 0d0
c
      do i1=1,n
       a = a + x(i1)*x(i1)
      enddo
c
      a = dsqrt(a)
c
      return
      end