lammps-sys 0.6.0

/*
//@HEADER
// ************************************************************************
//
//                        Kokkos v. 2.0
//              Copyright (2014) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Christian R. Trott (crtrott@sandia.gov)
//
// ************************************************************************
//@HEADER
*/

#ifndef KOKKOS_HEXEXPLICITFUNCTIONS_HPP
#define KOKKOS_HEXEXPLICITFUNCTIONS_HPP

#include <cmath>

namespace Explicit {

struct Hex8Functions
{
  static const unsigned SpatialDim    = 3 ;
  static const unsigned ElemNodeCount = 8 ;

  // Indices for full 3x3 tensor:

  static const unsigned K_F_XX = 0 ;
  static const unsigned K_F_YY = 1 ;
  static const unsigned K_F_ZZ = 2 ;
  static const unsigned K_F_XY = 3 ;
  static const unsigned K_F_YZ = 4 ;
  static const unsigned K_F_ZX = 5 ;
  static const unsigned K_F_YX = 6 ;
  static const unsigned K_F_ZY = 7 ;
  static const unsigned K_F_XZ = 8 ;
  static const unsigned K_F_SIZE = 9 ;

  //  Indexes into a 3 by 3 symmetric tensor stored as a length 6 vector

  static const unsigned K_S_XX = 0 ;
  static const unsigned K_S_YY = 1 ;
  static const unsigned K_S_ZZ = 2 ;
  static const unsigned K_S_XY = 3 ;
  static const unsigned K_S_YZ = 4 ;
  static const unsigned K_S_ZX = 5 ;
  static const unsigned K_S_YX = 3 ;
  static const unsigned K_S_ZY = 4 ;
  static const unsigned K_S_XZ = 5 ;
  static const unsigned K_S_SIZE = 6 ;

  //  Indexes into a 3 by 3 skew symmetric tensor stored as a length 3 vector

  static const unsigned K_V_XY = 0 ;
  static const unsigned K_V_YZ = 1 ;
  static const unsigned K_V_ZX = 2 ;
  static const unsigned K_V_SIZE = 3 ;

  //--------------------------------------------------------------------------

  template< typename ScalarA , typename ScalarB >
  KOKKOS_INLINE_FUNCTION static
  double dot8( const ScalarA * const a , const ScalarB * const b )
  { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3] +
           a[4] * b[4] + a[5] * b[5] + a[6] * b[6] + a[7] * b[7] ; }

  //--------------------------------------------------------------------------

  template< class ScalarPrecise ,
            class ScalarCompact >
  KOKKOS_INLINE_FUNCTION static
  void grad( const ScalarPrecise x[] ,
             const ScalarPrecise z[] ,
                   ScalarCompact grad_y[] )
  {
    const ScalarCompact R42=(x[3] - x[1]);
    const ScalarCompact R52=(x[4] - x[1]);
    const ScalarCompact R54=(x[4] - x[3]);

    const ScalarCompact R63=(x[5] - x[2]);
    const ScalarCompact R83=(x[7] - x[2]);
    const ScalarCompact R86=(x[7] - x[5]);

    const ScalarCompact R31=(x[2] - x[0]);
    const ScalarCompact R61=(x[5] - x[0]);
    const ScalarCompact R74=(x[6] - x[3]);

    const ScalarCompact R72=(x[6] - x[1]);
    const ScalarCompact R75=(x[6] - x[4]);
    const ScalarCompact R81=(x[7] - x[0]);

    const ScalarCompact t1=(R63 + R54);
    const ScalarCompact t2=(R61 + R74);
    const ScalarCompact t3=(R72 + R81);

    const ScalarCompact t4 =(R86 + R42);
    const ScalarCompact t5 =(R83 + R52);
    const ScalarCompact t6 =(R75 + R31);

    //  Calculate Y gradient from X and Z data

    grad_y[0] = (z[1] *  t1) - (z[2] * R42) - (z[3] *  t5)  + (z[4] *  t4) + (z[5] * R52) - (z[7] * R54);
    grad_y[1] = (z[2] *  t2) + (z[3] * R31) - (z[0] *  t1)  - (z[5] *  t6) + (z[6] * R63) - (z[4] * R61);
    grad_y[2] = (z[3] *  t3) + (z[0] * R42) - (z[1] *  t2)  - (z[6] *  t4) + (z[7] * R74) - (z[5] * R72);
    grad_y[3] = (z[0] *  t5) - (z[1] * R31) - (z[2] *  t3)  + (z[7] *  t6) + (z[4] * R81) - (z[6] * R83);
    grad_y[4] = (z[5] *  t3) + (z[6] * R86) - (z[7] *  t2)  - (z[0] *  t4) - (z[3] * R81) + (z[1] * R61);
    grad_y[5] = (z[6] *  t5) - (z[4] *  t3)  - (z[7] * R75) + (z[1] *  t6) - (z[0] * R52) + (z[2] * R72);
    grad_y[6] = (z[7] *  t1) - (z[5] *  t5)  - (z[4] * R86) + (z[2] *  t4) - (z[1] * R63) + (z[3] * R83);
    grad_y[7] = (z[4] *  t2) - (z[6] *  t1)  + (z[5] * R75) - (z[3] *  t6) - (z[2] * R74) + (z[0] * R54);
  }

  template< class ScalarPrecise ,
            class ScalarCompact >
  static KOKKOS_INLINE_FUNCTION
  void grad( const ScalarPrecise x[] ,
             const ScalarPrecise y[] ,
             const ScalarPrecise z[] ,
                   ScalarCompact grad_x[] ,
                   ScalarCompact grad_y[] ,
                   ScalarCompact grad_z[] )
  {
    grad( x , z , grad_y );
    grad( z , y , grad_x );
    grad( y , x , grad_z );
  }

  //--------------------------------------------------------------------------

  template< class ScalarPrecise ,
            class ScalarCompact >
  KOKKOS_INLINE_FUNCTION static
  void polar_decomp( const float dt ,
                     const ScalarCompact v_gr[] ,
                           ScalarPrecise stretch[] /* INOUT */ ,
                           ScalarCompact str_ten[] /* OUT */ ,
                           ScalarCompact rot[]     /* OUT */ )
  {
    const float dt_half = 0.5 * dt;

    ScalarCompact vort[ K_V_SIZE ];  // Vorticity

    //  Symmetric part
    str_ten[K_S_XX] = v_gr[K_F_XX];
    str_ten[K_S_YY] = v_gr[K_F_YY];
    str_ten[K_S_ZZ] = v_gr[K_F_ZZ];
    str_ten[K_S_XY] = 0.5 * ( v_gr[K_F_XY] + v_gr[K_F_YX] );
    str_ten[K_S_YZ] = 0.5 * ( v_gr[K_F_YZ] + v_gr[K_F_ZY] );
    str_ten[K_S_ZX] = 0.5 * ( v_gr[K_F_ZX] + v_gr[K_F_XZ] );

    //  Skew Symmetric part
    vort[K_V_XY] = 0.5 * ( v_gr[K_F_XY] - v_gr[K_F_YX] );
    vort[K_V_YZ] = 0.5 * ( v_gr[K_F_YZ] - v_gr[K_F_ZY] );
    vort[K_V_ZX] = 0.5 * ( v_gr[K_F_ZX] - v_gr[K_F_XZ] );

    //   calculate the rates of rotation via gauss elimination.

    ScalarCompact z1 = str_ten[K_S_XY] * stretch[K_S_ZX] -
                       str_ten[K_S_ZX] * stretch[K_S_XY] +
                       str_ten[K_S_YY] * stretch[K_S_YZ] -
                       str_ten[K_S_YZ] * stretch[K_S_YY] +
                       str_ten[K_S_YZ] * stretch[K_S_ZZ] -
                       str_ten[K_S_ZZ] * stretch[K_S_YZ];

    ScalarCompact z2 = str_ten[K_S_ZX] * stretch[K_S_XX] -
                       str_ten[K_S_XX] * stretch[K_S_ZX] +
                       str_ten[K_S_YZ] * stretch[K_S_XY] -
                       str_ten[K_S_XY] * stretch[K_S_YZ] +
                       str_ten[K_S_ZZ] * stretch[K_S_ZX] -
                       str_ten[K_S_ZX] * stretch[K_S_ZZ];

    ScalarCompact z3 = str_ten[K_S_XX] * stretch[K_S_XY] -
                       str_ten[K_S_XY] * stretch[K_S_XX] +
                       str_ten[K_S_XY] * stretch[K_S_YY] -
                       str_ten[K_S_YY] * stretch[K_S_XY] +
                       str_ten[K_S_ZX] * stretch[K_S_YZ] -
                       str_ten[K_S_YZ] * stretch[K_S_ZX];

    {
      //   forward elimination

      const ScalarCompact a1inv  = 1.0 / (stretch[K_S_YY] + stretch[K_S_ZZ]);
      const ScalarCompact a4BYa1 = -1 * stretch[K_S_XY] * a1inv;
      const ScalarCompact a2inv  = 1.0 / (stretch[K_S_ZZ] + stretch[K_S_XX] + stretch[K_S_XY] * a4BYa1);

     const ScalarCompact a5 =  -stretch[K_S_YZ] + stretch[K_S_ZX] * a4BYa1;

      z2 -= z1 * a4BYa1;
      const ScalarCompact a6BYa1 = -1 * stretch[K_S_ZX] * a1inv;
      const ScalarCompact a5BYa2 = a5 * a2inv;
      z3 -= z1 * a6BYa1 - z2 * a5BYa2;

      //   backward substitution -

      z3 /= (stretch[K_S_XX] + stretch[K_S_YY] + stretch[K_S_ZX] * a6BYa1 + a5 * a5BYa2);
      z2 = (z2 - a5 * z3) * a2inv;
      z1 = (z1*a1inv - a6BYa1 * z3 -a4BYa1 * z2);
    }

    //   calculate rotation rates - recall that spin_rate is an asymmetric tensor,
    //   so compute spin rate vector as dual of spin rate tensor,
    //   i.e   w_i = e_ijk * spin_rate_jk

    z1 += vort[K_V_YZ];
    z2 += vort[K_V_ZX];
    z3 += vort[K_V_XY];

    {
      //   update rotation tensor:
      //  1) premultiply old rotation tensor to get right-hand side.

      ScalarCompact r_XX = rot[K_F_XX] + dt_half*( z3 * rot[K_F_YX] - z2 * rot[K_F_ZX] );
      ScalarCompact r_YX = rot[K_F_YX] + dt_half*( z1 * rot[K_F_ZX] - z3 * rot[K_F_XX] );
      ScalarCompact r_ZX = rot[K_F_ZX] + dt_half*( z2 * rot[K_F_XX] - z1 * rot[K_F_YX] );
      ScalarCompact r_XY = rot[K_F_XY] + dt_half*( z3 * rot[K_F_YY] - z2 * rot[K_F_ZY] );
      ScalarCompact r_YY = rot[K_F_YY] + dt_half*( z1 * rot[K_F_ZY] - z3 * rot[K_F_XY] );
      ScalarCompact r_ZY = rot[K_F_ZY] + dt_half*( z2 * rot[K_F_XY] - z1 * rot[K_F_YY] );
      ScalarCompact r_XZ = rot[K_F_XZ] + dt_half*( z3 * rot[K_F_YZ] - z2 * rot[K_F_ZZ] );
      ScalarCompact r_YZ = rot[K_F_YZ] + dt_half*( z1 * rot[K_F_ZZ] - z3 * rot[K_F_XZ] );
      ScalarCompact r_ZZ = rot[K_F_ZZ] + dt_half*( z2 * rot[K_F_XZ] - z1 * rot[K_F_YZ] );


      //  2) solve for new rotation tensor via gauss elimination.
      //   forward elimination -

      const ScalarCompact a12 = - dt_half * z3;
      const ScalarCompact a13 =   dt_half * z2;
            ScalarCompact b32 = - dt_half * z1;
      const ScalarCompact a22inv = 1.0 / (1.0 + a12 * a12);

      const ScalarCompact a13a12 = a13*a12;
      const ScalarCompact a23 = b32 + a13a12;

      r_YX += r_XX * a12;
      r_YY += r_XY * a12;
      r_YZ += r_XZ * a12;

      b32 = (b32 - a13a12) * a22inv;

      r_ZX += r_XX * a13 + r_YX * b32;
      r_ZY += r_XY * a13 + r_YY * b32;
      r_ZZ += r_XZ * a13 + r_YZ * b32;

      //   backward substitution -

      const ScalarCompact a33inv = 1.0 / (1.0 + a13 * a13 + a23 * b32);

      rot[K_F_ZX] = r_ZX * a33inv;
      rot[K_F_ZY] = r_ZY * a33inv;
      rot[K_F_ZZ] = r_ZZ * a33inv;
      rot[K_F_YX] = ( r_YX - rot[K_F_ZX] * a23 ) * a22inv;
      rot[K_F_YY] = ( r_YY - rot[K_F_ZY] * a23 ) * a22inv;
      rot[K_F_YZ] = ( r_YZ - rot[K_F_ZZ] * a23 ) * a22inv;
      rot[K_F_XX] = r_XX - rot[K_F_ZX] * a13 - rot[K_F_YX] * a12;
      rot[K_F_XY] = r_XY - rot[K_F_ZY] * a13 - rot[K_F_YY] * a12;
      rot[K_F_XZ] = r_XZ - rot[K_F_ZZ] * a13 - rot[K_F_YZ] * a12;
    }

    //   update stretch tensor in the new configuration -

    const ScalarCompact a1 = str_ten[K_S_XY] + vort[K_V_XY];
    const ScalarCompact a2 = str_ten[K_S_YZ] + vort[K_V_YZ];
    const ScalarCompact a3 = str_ten[K_S_ZX] + vort[K_V_ZX];
    const ScalarCompact b1 = str_ten[K_S_ZX] - vort[K_V_ZX];
    const ScalarCompact b2 = str_ten[K_S_XY] - vort[K_V_XY];
    const ScalarCompact b3 = str_ten[K_S_YZ] - vort[K_V_YZ];

    const ScalarCompact s_XX = stretch[K_S_XX];
    const ScalarCompact s_YY = stretch[K_S_YY];
    const ScalarCompact s_ZZ = stretch[K_S_ZZ];
    const ScalarCompact s_XY = stretch[K_S_XY];
    const ScalarCompact s_YZ = stretch[K_S_YZ];
    const ScalarCompact s_ZX = stretch[K_S_ZX];

    stretch[K_S_XX] += dt * (str_ten[K_S_XX] * s_XX + ( a1 + z3 ) * s_XY + ( b1 - z2 ) * s_ZX);
    stretch[K_S_YY] += dt * (str_ten[K_S_YY] * s_YY + ( a2 + z1 ) * s_YZ + ( b2 - z3 ) * s_XY);
    stretch[K_S_ZZ] += dt * (str_ten[K_S_ZZ] * s_ZZ + ( a3 + z2 ) * s_ZX + ( b3 - z1 ) * s_YZ);
    stretch[K_S_XY] += dt * (str_ten[K_S_XX] * s_XY + ( a1 )      * s_YY + ( b1      ) * s_YZ - z3 * s_XX + z1 * s_ZX);
    stretch[K_S_YZ] += dt * (str_ten[K_S_YY] * s_YZ + ( a2 )      * s_ZZ + ( b2      ) * s_ZX - z1 * s_YY + z2 * s_XY);
    stretch[K_S_ZX] += dt * (str_ten[K_S_ZZ] * s_ZX + ( a3 )      * s_XX + ( b3      ) * s_XY - z2 * s_ZZ + z3 * s_YZ);
  }

  //--------------------------------------------------------------------------

  template< typename ScalarCompact >
  static KOKKOS_INLINE_FUNCTION
  void rotate_tensor( const ScalarCompact str_ten[] ,
                      const ScalarCompact rot[] ,
                            ScalarCompact rot_str[] )
  {
    ScalarCompact t[9];

    t[0] = str_ten[K_S_XX]*rot[K_F_XX] + str_ten[K_S_XY]*rot[K_F_YX] + str_ten[K_S_XZ]*rot[K_F_ZX];
    t[1] = str_ten[K_S_YX]*rot[K_F_XX] + str_ten[K_S_YY]*rot[K_F_YX] + str_ten[K_S_YZ]*rot[K_F_ZX];
    t[2] = str_ten[K_S_ZX]*rot[K_F_XX] + str_ten[K_S_ZY]*rot[K_F_YX] + str_ten[K_S_ZZ]*rot[K_F_ZX];

    t[3] = str_ten[K_S_XX]*rot[K_F_XY] + str_ten[K_S_XY]*rot[K_F_YY] + str_ten[K_S_XZ]*rot[K_F_ZY];
    t[4] = str_ten[K_S_YX]*rot[K_F_XY] + str_ten[K_S_YY]*rot[K_F_YY] + str_ten[K_S_YZ]*rot[K_F_ZY];
    t[5] = str_ten[K_S_ZX]*rot[K_F_XY] + str_ten[K_S_ZY]*rot[K_F_YY] + str_ten[K_S_ZZ]*rot[K_F_ZY];

    t[6] = str_ten[K_S_XX]*rot[K_F_XZ] + str_ten[K_S_XY]*rot[K_F_YZ] + str_ten[K_S_XZ]*rot[K_F_ZZ];
    t[7] = str_ten[K_S_YX]*rot[K_F_XZ] + str_ten[K_S_YY]*rot[K_F_YZ] + str_ten[K_S_YZ]*rot[K_F_ZZ];
    t[8] = str_ten[K_S_ZX]*rot[K_F_XZ] + str_ten[K_S_ZY]*rot[K_F_YZ] + str_ten[K_S_ZZ]*rot[K_F_ZZ];


    rot_str[ K_S_XX ] = rot[K_F_XX] * t[0] + rot[K_F_YX] * t[1] + rot[K_F_ZX] * t[2];
    rot_str[ K_S_YY ] = rot[K_F_XY] * t[3] + rot[K_F_YY] * t[4] + rot[K_F_ZY] * t[5];
    rot_str[ K_S_ZZ ] = rot[K_F_XZ] * t[6] + rot[K_F_YZ] * t[7] + rot[K_F_ZZ] * t[8];

    rot_str[ K_S_XY ] = rot[K_F_XX] * t[3] + rot[K_F_YX] * t[4] + rot[K_F_ZX] * t[5];
    rot_str[ K_S_YZ ] = rot[K_F_XY] * t[6] + rot[K_F_YY] * t[7] + rot[K_F_ZY] * t[8];
    rot_str[ K_S_ZX ] = rot[K_F_XZ] * t[0] + rot[K_F_YZ] * t[1] + rot[K_F_ZZ] * t[2];
  }

  //--------------------------------------------------------------------------

  template< class ScalarPrecise ,
            class ScalarCompact >
  static KOKKOS_INLINE_FUNCTION
  void rotate_tensor_backward( const ScalarPrecise stress[] ,
                               const ScalarCompact rot[] ,
                                     ScalarCompact rot_stress[] )
  {
    ScalarCompact t[9] ;

    t[0] = stress[K_S_XX]*rot[K_F_XX]+ stress[K_S_XY]*rot[K_F_XY]+ stress[K_S_XZ]*rot[K_F_XZ];
    t[1] = stress[K_S_YX]*rot[K_F_XX]+ stress[K_S_YY]*rot[K_F_XY]+ stress[K_S_YZ]*rot[K_F_XZ];
    t[2] = stress[K_S_ZX]*rot[K_F_XX]+ stress[K_S_ZY]*rot[K_F_XY]+ stress[K_S_ZZ]*rot[K_F_XZ];
    t[3] = stress[K_S_XX]*rot[K_F_YX]+ stress[K_S_XY]*rot[K_F_YY]+ stress[K_S_XZ]*rot[K_F_YZ];
    t[4] = stress[K_S_YX]*rot[K_F_YX]+ stress[K_S_YY]*rot[K_F_YY]+ stress[K_S_YZ]*rot[K_F_YZ];
    t[5] = stress[K_S_ZX]*rot[K_F_YX]+ stress[K_S_ZY]*rot[K_F_YY]+ stress[K_S_ZZ]*rot[K_F_YZ];
    t[6] = stress[K_S_XX]*rot[K_F_ZX]+ stress[K_S_XY]*rot[K_F_ZY]+ stress[K_S_XZ]*rot[K_F_ZZ];
    t[7] = stress[K_S_YX]*rot[K_F_ZX]+ stress[K_S_YY]*rot[K_F_ZY]+ stress[K_S_YZ]*rot[K_F_ZZ];
    t[8] = stress[K_S_ZX]*rot[K_F_ZX]+ stress[K_S_ZY]*rot[K_F_ZY]+ stress[K_S_ZZ]*rot[K_F_ZZ];

    rot_stress[ K_S_XX ] = rot[K_F_XX]*t[0] + rot[K_F_XY]*t[1] + rot[K_F_XZ]*t[2];
    rot_stress[ K_S_YY ] = rot[K_F_YX]*t[3] + rot[K_F_YY]*t[4] + rot[K_F_YZ]*t[5];
    rot_stress[ K_S_ZZ ] = rot[K_F_ZX]*t[6] + rot[K_F_ZY]*t[7] + rot[K_F_ZZ]*t[8];

    rot_stress[ K_S_XY ] = rot[K_F_XX]*t[3] + rot[K_F_XY]*t[4] + rot[K_F_XZ]*t[5];
    rot_stress[ K_S_YZ ] = rot[K_F_YX]*t[6] + rot[K_F_YY]*t[7] + rot[K_F_YZ]*t[8];
    rot_stress[ K_S_ZX ] = rot[K_F_ZX]*t[0] + rot[K_F_ZY]*t[1] + rot[K_F_ZZ]*t[2];
  }

  //--------------------------------------------------------------------------

  template< class ScalarPrecise ,
            class ScalarCompact >
  KOKKOS_INLINE_FUNCTION static
  void update_stress( const float dt ,
                      const float two_mu ,
                      const float bulk_modulus ,
                      const ScalarCompact rot_str[] ,
                            ScalarPrecise stress[] )
  {
    const ScalarCompact e = rot_str[ K_S_XX ] + rot_str[ K_S_YY ] + rot_str[ K_S_ZZ ] ;
    const ScalarCompact eb = e * bulk_modulus ;
    const ScalarCompact e3 = e / 3.0 ;

    stress[K_S_XX] += dt * ( two_mu * ( rot_str[K_S_XX] - e3 ) + eb );
    stress[K_S_YY] += dt * ( two_mu * ( rot_str[K_S_YY] - e3 ) + eb );
    stress[K_S_ZZ] += dt * ( two_mu * ( rot_str[K_S_ZZ] - e3 ) + eb );

    stress[K_S_XY] += dt * two_mu * rot_str[K_S_XY];
    stress[K_S_YZ] += dt * two_mu * rot_str[K_S_YZ];
    stress[K_S_ZX] += dt * two_mu * rot_str[K_S_ZX];
  }

  //--------------------------------------------------------------------------

  template< class ScalarPrecise ,
            class ScalarCompact >
  static KOKKOS_INLINE_FUNCTION
  void comp_force( const ScalarPrecise vx[] ,
                   const ScalarPrecise vy[] ,
                   const ScalarPrecise vz[] ,
                   const ScalarCompact grad_x[] ,
                   const ScalarCompact grad_y[] ,
                   const ScalarCompact grad_z[] ,
                   const ScalarCompact total_stress12th[] ,
                         ScalarCompact force[][ SpatialDim ] ,
                         ScalarCompact & energy )
  {
    ScalarPrecise internal_energy = 0 ;

    for ( unsigned inode = 0; inode < ElemNodeCount ; ++inode ) {

      force[inode][0] = total_stress12th[K_S_XX] * grad_x[inode] +
                        total_stress12th[K_S_XY] * grad_y[inode] +
                        total_stress12th[K_S_XZ] * grad_z[inode] ;

      force[inode][1] = total_stress12th[K_S_YX] * grad_x[inode] +
                        total_stress12th[K_S_YY] * grad_y[inode] +
                        total_stress12th[K_S_YZ] * grad_z[inode] ;

      force[inode][2] = total_stress12th[K_S_ZX] * grad_x[inode] +
                        total_stress12th[K_S_ZY] * grad_y[inode] +
                        total_stress12th[K_S_ZZ] * grad_z[inode] ;

      internal_energy += force[inode][0] * vx[inode] +
                         force[inode][1] * vy[inode] +
                         force[inode][2] * vz[inode] ;
    }

    energy = internal_energy ;
  }

  //--------------------------------------------------------------------------
};

} // namespace Explicit

#endif /* #ifndef KOKKOS_HEXEXPLICITFUNCTIONS_HPP */