subroutine elevar(nm,nnm,layer,gstran,gstres,stress,pstran,
     1cstran,vstran,cauchy,eplas,equivc,swell,krtyp,prang,dt,
     2gsv,ngens1,ngen1,nstats,nstass,thmstr)
c
c CONTACT: ap@marc.de
c      
c* * * * *
c
c User routine to compute model inertia quantities.
c It will print:
c   1) Total mass of model.
c   2) The coordinates of the center of mass.
c   3) The inertia tensor about axis through the center of mass
c      which are parallel to the global coordinate axes.
c   4) The principal moments of inertia and its direction cosines.
c   5) The moment of inertia about an axis specified in the
c      "ROTATION A" model definition option.
c
c NOTE: The routine is restricted to volume elements (like 7, 21),
c       shell elements (like 75, 22), membrane elements (like 18),
c       elastic beam elements (like 52, 98) and truss elements 
c       (like 9, 64).
c       This subroutine will only furnish useful results if the
c       model is a real 3D model (i.e. #coordinates=3).
c       This subroutine must be used in conjunction with subroutine
c       'ubginc' which is attached below.
c       The routine is compliant with marc versions k62, k71 and k72.
c
c TEST STATUS: No tests were done for higher order elements.
c
c LAST MODIFICATION DATE: May 28, 1998
c
c* * * * *
      implicit real*8 (a-h,o-z)                                               dp
      dimension gstran(ngens),gstres(ngens),
     1stress(ngen1),pstran(ngen1),cstran(ngen1),vstran(ngen1),
     2cauchy(ngen1),dt(nstats),gsv(1)
     3,thmstr(1),crang(3,3),krtyp(4)
c
      logical debug
      data debug/.false./
c
      dimension xip(3),xcm(3),xv(6),xp(3),rp(3,3)
c* * * * *
c xip = integration point coordinates
c xcm = coordinates of center of mass
c xv  = workspace vector
c xp  = principal inertia values (xp(1) > xp(2) > xp(3))
c rp  = direction cosines for principal inertia values
c       rp(ij,i1) direction cosine vector of principal value i1
c* * * * *
c
c The common "inertia_com" saves data over calls of elevar.
      common/inertia_com/sm(3),xa(3,3),xl1(3),xl2(3),xl,xm,ncrc,nmel
c* * * * *
c sm   = static moment vector
c xa   = inertia tensor about X,Y,Z-axis through center of mass
c xl1  = 1st point of desired inertia axis
c xl2  = 2nd point of desired inertia axis
c xl   = moment of inertia about desired axis
c xm   = mass of structure
c ncrc = number of space coordinates
c nmel = element counter
c* * * * *
c
      include '../common/strnen'
      include '../common/arrays'
      include '../common/array4'
      include '../common/heat'
      include '../common/dimen'
      include '../common/elmcom'
      include '../common/matdat'
      include '../common/space'
      include '../common/lass'
c
c... return if layer number not equal to one
      if (layer.gt.1) goto 99999
c
c... update the element counter
      if (nnm.eq.1) nmel = nmel + 1
c
c... initialize the inertia value
      if (nmel.eq.1.and.nnm.eq.1.and.layer.eq.1) then
       xm = 0.0d0
       xl = 0.0d0
       do i1=1,3
        sm(i1) = 0.0d0
        do i2=1,3
         xa(i1,i2) = 0.0d0
        enddo
       enddo
c
c... set the space dimension
      ncrc=ncrd
      if (ncrc.gt.3) ncrc=3
      if (ncrc.lt.3) then
       write(6,'(/,15x,a,/)') 'WARNING: cannot handle 2D; skipped'
       goto 99999
      endif
c
c... xl1, xl2 define 1st and 2nd point of desired inertia axis
c
c... get point on the rotation axis as the 1st point
       xl1(1) = vars(icentr)
       xl1(2) = vars(icentr+1)
       xl1(3) = vars(icentr+2)
       if (debug)
     *   write(0,'(a20,3e14.5)') '1st point: ',(xl1(ijk),ijk=1,3)
c
c... get direction cosines of rotation axis and define the 2nd point
       xl2(1) = xl1(1) + vars(icenta)
       xl2(2) = xl1(2) + vars(icenta+1)
       xl2(3) = xl1(3) + vars(icenta+2)
       if (debug)
     *   write(0,'(a20,3e14.5)') '2nd point: ',(xl2(ijk),ijk=1,3)
c
c... check if 1st and 2nd point coincide (if they do then define the
c... 2nd point in global z-direction above the 1st point)
       dd = (xl1(1)-xl2(1))*(xl1(1)-xl2(1)) +
     *      (xl1(2)-xl2(2))*(xl1(2)-xl2(2)) +
     *      (xl1(3)-xl2(3))*(xl1(3)-xl2(3)) 
       if (dd.lt.1d-10) then
        xl2(1) = xl1(1)
        xl2(2) = xl1(2)
        xl2(3) = xl1(3) + 1d0
       endif
      endif
c
c... data base offset (here use internal element number)
      lofr = (n-1)*nelstr
c
c... get the integration point coordinates
      la1  = icrxpt + lofr + (nn-1)*ncrc
      do 100,i1=1,ncrc
       xip(i1) = vars(la1+i1-1)
 100  continue
c
c... d2 = distance squared to desired inertia axis
      p1 = 0.0d0
      p2 = 0.0d0
      p3 = 0.0d0
      do 200, i1=1,ncrc
       p1 = p1 + (xip(i1)-xl1(i1))**2
       p2 = p2 + (xip(i1)-xl1(i1))*(xl2(i1)-xl1(i1))
       p3 = p3 + (xl2(i1)-xl1(i1))**2
 200  continue
      d2 = p1 - p2*p2/p3
c
c... dv = the jacobian (the integration point volume factor)
      la1   = ijacoc  + lofr + nn - 1
      dv    = vars(la1)
c
c... check for different element classes to compute area or thickness
      if (lclass.eq.8) then
c... volume elements
       thick = 1d0
      else if (lclass.eq.2) then
c... shell elements (get thickness)
       la1  = ithick + lofr + nn - 1
       thick = vars(la1)
      else if (lclass.eq.1) then
c... truss & elastic beam elements (get cross sectional area)
       la1  = igeomr + lofr
       thick = vars(la1)
      else
c... remaining classes are not accounted for
       thick = 0d0
       write(6,1111) nm,lclass
      end if
      if (debug) then
       write(0,'(a20,2i14)') 'ELEMENT, IP: ',nm,nnm
       write(0,'(a20,i14,2e14.5)') 'lclass: ',lclass,thick,rho
       write(0,'(a20,3e14.5)') 'IP coordinates',(xip(ijk),ijk=1,3)
      endif
c
 1111 format('Element',i7,': this element class is not accounted for
     * (class',i3,')')
c
c... scale jacobian with shell thickness or beam area
      dv = dv*thick
c
c... sum the integration point contributions
      xl      = xl + rho*d2*dv
      xm      = xm + rho*dv
      sm(1)   = sm(1) + rho*xip(1)*dv
      sm(2)   = sm(2) + rho*xip(2)*dv
      sm(3)   = sm(3) + rho*xip(3)*dv
      xa(1,1) = xa(1,1) + rho*(xip(2)**2+xip(3)**2)*dv
      xa(2,2) = xa(2,2) + rho*(xip(1)**2+xip(3)**2)*dv
      xa(3,3) = xa(3,3) + rho*(xip(1)**2+xip(2)**2)*dv
      xa(1,2) = xa(1,2) - rho*xip(1)*xip(2)*dv
      xa(2,3) = xa(2,3) - rho*xip(2)*xip(3)*dv
      xa(1,3) = xa(1,3) - rho*xip(1)*xip(3)*dv
c
c... print the sum when processing last element and integration point
      if (nmel.eq.numel.and.nnm.eq.intel) then
       xcm(1) = sm(1)/xm
       xcm(2) = sm(2)/xm
       xcm(3) = sm(3)/xm
       xa(1,1) = xa(1,1) - xm*(xcm(2)**2+xcm(3)**2)
       xa(2,2) = xa(2,2) - xm*(xcm(1)**2+xcm(3)**2)
       xa(3,3) = xa(3,3) - xm*(xcm(1)**2+xcm(2)**2)
       xa(1,2) = xa(1,2) + xm*xcm(1)*xcm(2) 
       xa(2,3) = xa(2,3) + xm*xcm(2)*xcm(3) 
       xa(1,3) = xa(1,3) + xm*xcm(1)*xcm(3) 
       xa(2,1) = xa(1,2)
       xa(3,2) = xa(2,3)
       xa(3,1) = xa(1,3)
c... compute principal values
       xv(1) = xa(1,1)
       xv(2) = xa(2,2)
       xv(3) = xa(3,3)
       xv(4) = xa(1,2)
       xv(5) = xa(2,3)
       xv(6) = xa(1,3)
       call princv(xp,rp,xv,3,3,0,0,0,0) 
c... sort principal values according to size
       i1 = 1
       i2 = 2
       i3 = 3
       if (xp(i2).gt.xp(i1)) then
        ii = i1
        i1 = i2
        i2 = ii
       end if
       if (xp(i3).gt.xp(i1)) then
        ii = i1
        i1 = i3
        i3 = ii
       end if
       if (xp(i3).gt.xp(i2)) then
        ii = i2
        i2 = i3
        i3 = ii
       end if
       write(6,5000)
c... print the total mass
       write(6,9000) xm
c... print the center of mass
       write(6,6000) (xcm(j1),j1=1,3)
c... print the inertia tensor w.r.t. center of mass
       write(6,8000) ((xa(j1,j2),j1=1,3),j2=1,3)
c... print the principal values
       write(6,1000) xp(i1)
       write(6,1100) (rp(j1,i1),j1=1,3)
       write(6,2000) xp(i2)
       write(6,2100) (rp(j1,i2),j1=1,3)
       write(6,3000) xp(i3)
       write(6,3100) (rp(j1,i3),j1=1,3)
c... print moment of inertia about desired axis
       write(6,7000) (xl1(j1),j1=1,3),(xl2(j1),j1=1,3),xl
       write(6,5000)
      endif
c
 9000 format(/,15x,'Total mass of the model: ',e14.6)
 6000 format(/,15x,'Center of mass X:',e14.6,' Y:',e14.6,' Z:',e14.6)
 8000 format(/,15x,'Inertia tensor w.r.t. center of mass:',/
     *17x,'Ixx:',e14.6,' Ixy:',e14.6,' Ixz:',e14.6,/,
     *17x,'Iyx:',e14.6,' Iyy:',e14.6,' Iyz:',e14.6,/,
     *17x,'Izx:',e14.6,' Izy:',e14.6,' Izz:',e14.6)
 1000 format(/,15x,'First  principal moment of inertia = ',e14.6)
 2000 format(/,15x,'Second principal moment of inertia = ',e14.6)
 3000 format(/,15x,'Third  principal moment of inertia = ',e14.6)
 1100 format(15x,'First  principal direction cosine vector: ', 3e14.6)
 2100 format(15x,'Second principal direction cosine vector: ', 3e14.6)
 3100 format(15x,'Third  principal direction cosine vector: ', 3e14.6)
 5000 format(/)
 7000 format(/,15x,'Moment of inertia about axis through points',/,
     *17x,'X1:',e14.6,' Y1:',e14.6,' Z1:',e14.6,' and',/,
     *17x,'X2:',e14.6,' Y2:',e14.6,' Z2:',e14.6,/,
     *15x,'is: ',e14.6)
c
99999 continue
      return
c
      end
      subroutine ubginc(inc,incsub)
      implicit real*8 (a-h,o-z)
c
c     subroutine that gets called at the beginning of each increment.
c
      common/inertia_com/sm(3),xa(3,3),xl1(3),xl2(3),xl,xm,ncrc,nmel
c
c... initialize the element counter at the start of the increment.
      nmel = 0
c
      return
      end