/*
//@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_EXPLICITFUNCTORS_HPP
#define KOKKOS_EXPLICITFUNCTORS_HPP

#include <cmath>
#include <Kokkos_Core.hpp>
#include <FEMesh.hpp>

namespace Explicit {

template<typename Scalar , class Device >
struct Fields {

  static const int NumStates     = 2 ;
  static const int SpatialDim    = 3 ;
  static const int ElemNodeCount = 8 ;

  // Indices for full 3x3 tensor:

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

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

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

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

  static const int K_V_XY = 0 ;
  static const int K_V_YZ = 1 ;
  static const int K_V_ZX = 2 ;


  typedef Device                           execution_space ;
  typedef typename execution_space::size_type  size_type ;

  typedef HybridFEM::FEMesh<double,ElemNodeCount,execution_space>  FEMesh ;

  typedef typename FEMesh::node_coords_type    node_coords_type ;
  typedef typename FEMesh::elem_node_ids_type  elem_node_ids_type ;
  typedef typename FEMesh::node_elem_ids_type  node_elem_ids_type ;
  typedef typename Kokkos::ParallelDataMap   parallel_data_map ;

  typedef Kokkos::View< double[][ SpatialDim ][ NumStates ] , execution_space > geom_state_array_type ;
  typedef Kokkos::View< Scalar[][ SpatialDim ] , execution_space > geom_array_type ;
  typedef Kokkos::View< Scalar[] ,               execution_space > array_type ;
  typedef Kokkos::View< Scalar ,                 execution_space >  scalar_type ;

  typedef Kokkos::View< Scalar[][  6 ] ,    execution_space >  elem_sym_tensor_type ;
  typedef Kokkos::View< Scalar[][  9 ] ,    execution_space >  elem_tensor_type ;
  typedef Kokkos::View< Scalar[][  9 ][ NumStates ] , execution_space >  elem_tensor_state_type ;
  typedef Kokkos::View< Scalar[][ SpatialDim ][ ElemNodeCount ] , execution_space > elem_node_geom_type ;

  // Parameters:
  const int num_nodes ;
  const int num_nodes_owned ;
  const int num_elements ;

  const Scalar  lin_bulk_visc;
  const Scalar  quad_bulk_visc;
  const Scalar  two_mu;
  const Scalar  bulk_modulus;
  const Scalar  density;

  // Mesh:
  const elem_node_ids_type  elem_node_connectivity ;
  const node_elem_ids_type  node_elem_connectivity ;
  const node_coords_type    model_coords ;

  // Compute:
  const scalar_type                dt ;
  const scalar_type                prev_dt ;
  const geom_state_array_type      displacement ;
  const geom_state_array_type      velocity ;
  const geom_array_type            acceleration ;
  const geom_array_type            internal_force ;
  const array_type                 nodal_mass ;
  const array_type                 elem_mass ;
  const array_type                 internal_energy ;
  const elem_sym_tensor_type       stress_new ;
  const elem_tensor_state_type     rotation ;
  const elem_node_geom_type        element_force ;
  const elem_tensor_type           vel_grad ;
  const elem_sym_tensor_type       stretch ;
  const elem_sym_tensor_type       rot_stretch ;

  Fields(
      const FEMesh & mesh,
      Scalar arg_lin_bulk_visc,
      Scalar arg_quad_bulk_visc,
      Scalar youngs_modulus,
      Scalar poissons_ratio,
      Scalar arg_density )
    : num_nodes(       mesh.parallel_data_map.count_owned +
                       mesh.parallel_data_map.count_receive )
    , num_nodes_owned( mesh.parallel_data_map.count_owned )
    , num_elements(    mesh.elem_node_ids.dimension_0() )
    , lin_bulk_visc(  arg_lin_bulk_visc )
    , quad_bulk_visc( arg_quad_bulk_visc )
    , two_mu(youngs_modulus/(1.0+poissons_ratio))
    , bulk_modulus(youngs_modulus/(3*(1.0-2.0*poissons_ratio)))
    , density(arg_density)

    // mesh

    , elem_node_connectivity( mesh.elem_node_ids ) // ( num_elements , ElemNodeCount )
    , node_elem_connectivity( mesh.node_elem_ids ) // ( num_nodes , ... )
    , model_coords(  mesh.node_coords )            // ( num_nodes , 3 )

    // compute with input/output

    , dt(              "dt" )
    , prev_dt(         "prev_dt" )
    , displacement(    "displacement" ,   num_nodes )
    , velocity(        "velocity" ,       num_nodes )
    , acceleration(    "acceleration" ,   num_nodes_owned )
    , internal_force(  "internal_force" , num_nodes_owned )
    , nodal_mass(      "nodal_mass" ,     num_nodes_owned )
    , elem_mass(       "elem_mass" ,       num_elements )
    , internal_energy( "internal_energy" , num_elements )
    , stress_new(      "stress_new" ,      num_elements )

    // temporary arrays

    , rotation(      "rotation" ,  num_elements )
    , element_force( "element_force" ,  num_elements )
    , vel_grad(      "vel_grad" , num_elements )
    , stretch(       "stretch" , num_elements )
    , rot_stretch(   "rot_stretch" , num_elements )
  { }
};


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

template< typename Scalar , class DeviceType >
KOKKOS_INLINE_FUNCTION
Scalar dot8( const Scalar * a , const Scalar * 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< typename Scalar , class DeviceType >
KOKKOS_INLINE_FUNCTION
void comp_grad( const Scalar * const x ,
                const Scalar * const y ,
                const Scalar * const z,
                Scalar * const grad_x ,
                Scalar * const grad_y ,
                Scalar * const grad_z )
{
  //  calc X difference vectors

  Scalar R42=(x[3] - x[1]);
  Scalar R52=(x[4] - x[1]);
  Scalar R54=(x[4] - x[3]);

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

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

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

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

  Scalar t4 =(R86 + R42);
  Scalar t5 =(R83 + R52);
  Scalar 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);

  //   calc Z difference vectors

  R42=(z[3] - z[1]);
  R52=(z[4] - z[1]);
  R54=(z[4] - z[3]);

  R63=(z[5] - z[2]);
  R83=(z[7] - z[2]);
  R86=(z[7] - z[5]);

  R31=(z[2] - z[0]);
  R61=(z[5] - z[0]);
  R74=(z[6] - z[3]);

  R72=(z[6] - z[1]);
  R75=(z[6] - z[4]);
  R81=(z[7] - z[0]);

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

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

  //  Calculate X gradient from Y and Z data

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

  //  calc Y difference vectors

  R42=(y[3] - y[1]);
  R52=(y[4] - y[1]);
  R54=(y[4] - y[3]);

  R63=(y[5] - y[2]);
  R83=(y[7] - y[2]);
  R86=(y[7] - y[5]);

  R31=(y[2] - y[0]);
  R61=(y[5] - y[0]);
  R74=(y[6] - y[3]);

  R72=(y[6] - y[1]);
  R75=(y[6] - y[4]);
  R81=(y[7] - y[0]);

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

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

  //  Calculate Z gradient from X and Y data

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

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

template< typename Scalar , class DeviceType >
struct initialize_element
{
  typedef DeviceType     execution_space ;

  typedef Explicit::Fields< Scalar , execution_space > Fields ;

  typename Fields::elem_node_ids_type      elem_node_connectivity ;
  typename Fields::node_coords_type        model_coords ;
  typename Fields::elem_sym_tensor_type    stretch ;
  typename Fields::elem_tensor_state_type  rotation ;
  typename Fields::array_type              elem_mass ;

  const Scalar density ;

  initialize_element( const Fields & mesh_fields )
    : elem_node_connectivity( mesh_fields.elem_node_connectivity )
    , model_coords(           mesh_fields.model_coords )
    , stretch(                mesh_fields.stretch )
    , rotation(               mesh_fields.rotation )
    , elem_mass(              mesh_fields.elem_mass )
    , density(                mesh_fields.density )
    {}

  KOKKOS_INLINE_FUNCTION
  void operator()( int ielem )const
  {
    const int K_XX = 0 ;
    const int K_YY = 1 ;
    const int K_ZZ = 2 ;
    const Scalar ONE12TH = 1.0 / 12.0 ;

    Scalar x[ Fields::ElemNodeCount ];
    Scalar y[ Fields::ElemNodeCount ];
    Scalar z[ Fields::ElemNodeCount ];
    Scalar grad_x[ Fields::ElemNodeCount ];
    Scalar grad_y[ Fields::ElemNodeCount ];
    Scalar grad_z[ Fields::ElemNodeCount ];

    for ( int i = 0 ; i < Fields::ElemNodeCount ; ++i ) {
      const int n = elem_node_connectivity( ielem , i );

      x[i]  = model_coords( n , 0 );
      y[i]  = model_coords( n , 1 );
      z[i]  = model_coords( n , 2 );
    }

    comp_grad<Scalar,execution_space>( x, y, z, grad_x, grad_y, grad_z);

    stretch(ielem,K_XX) = 1 ;
    stretch(ielem,K_YY) = 1 ;
    stretch(ielem,K_ZZ) = 1 ;

    rotation(ielem,K_XX,0) = 1 ;
    rotation(ielem,K_YY,0) = 1 ;
    rotation(ielem,K_ZZ,0) = 1 ;

    rotation(ielem,K_XX,1) = 1 ;
    rotation(ielem,K_YY,1) = 1 ;
    rotation(ielem,K_ZZ,1) = 1 ;

    elem_mass(ielem) = ONE12TH * density *
                                 dot8<Scalar,execution_space>( x , grad_x );
  }

  static void apply( const Fields & mesh_fields )
  {
    initialize_element op( mesh_fields );
    Kokkos::parallel_for( mesh_fields.num_elements , op );
  }
};


template<typename Scalar , class DeviceType >
struct initialize_node
{
  typedef DeviceType     execution_space ;

  typedef Explicit::Fields< Scalar , execution_space > Fields ;

  typename Fields::node_elem_ids_type      node_elem_connectivity ;
  typename Fields::array_type              nodal_mass ;
  typename Fields::array_type              elem_mass ;

  static const int ElemNodeCount = Fields::ElemNodeCount ;

  initialize_node( const Fields & mesh_fields )
    : node_elem_connectivity( mesh_fields.node_elem_connectivity )
    , nodal_mass(             mesh_fields.nodal_mass )
    , elem_mass(              mesh_fields.elem_mass )
    {}


  KOKKOS_INLINE_FUNCTION
  void operator()( int inode )const
  {
    const int begin = node_elem_connectivity.row_map[inode];
    const int end   = node_elem_connectivity.row_map[inode+1];

    Scalar node_mass = 0;

    for(int i = begin; i != end; ++i) {
      const int elem_id = node_elem_connectivity.entries( i , 0 );
      node_mass += elem_mass(elem_id);
    }

    nodal_mass(inode) = node_mass / ElemNodeCount ;
  }

  static void apply( const Fields & mesh_fields )
  {
    initialize_node op( mesh_fields );
    Kokkos::parallel_for( mesh_fields.num_nodes_owned , op );
  }
};

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


template<typename Scalar, class DeviceType >
struct grad
{
  typedef DeviceType execution_space ;

  typedef Explicit::Fields< Scalar , execution_space >  Fields ;

  static const int ElemNodeCount = Fields::ElemNodeCount ;

  static const int K_F_XX = Fields::K_F_XX ;
  static const int K_F_YY = Fields::K_F_YY ;
  static const int K_F_ZZ = Fields::K_F_ZZ ;
  static const int K_F_XY = Fields::K_F_XY ;
  static const int K_F_YZ = Fields::K_F_YZ ;
  static const int K_F_ZX = Fields::K_F_ZX ;
  static const int K_F_YX = Fields::K_F_YX ;
  static const int K_F_ZY = Fields::K_F_ZY ;
  static const int K_F_XZ = Fields::K_F_XZ ;

  // Global arrays used by this functor.

  const typename Fields::elem_node_ids_type     elem_node_connectivity ;
  const typename Fields::node_coords_type       model_coords ;
  const typename Fields::geom_state_array_type  displacement ;
  const typename Fields::geom_state_array_type  velocity ;
  const typename Fields::elem_tensor_type       vel_grad ;
  const typename Fields::scalar_type            dt ;

  const int  current_state;
  const int  previous_state;

  // Constructor on the Host to populate this device functor.
  // All array view copies are shallow.
  grad( const Fields &  fields,
        const int arg_current_state,
        const int arg_previous_state)
    : elem_node_connectivity( fields.elem_node_connectivity)
    , model_coords( fields.model_coords)
    , displacement( fields.displacement)
    , velocity( fields.velocity)
    , vel_grad( fields.vel_grad)
    , dt(  fields.dt)
    , current_state(arg_current_state)
    , previous_state(arg_previous_state)
    { }

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

    //   Calculate Velocity Gradients
    KOKKOS_INLINE_FUNCTION
    void v_grad(  int ielem,
      Scalar * vx,       Scalar * vy,       Scalar * vz,
      Scalar * grad_x,     Scalar * grad_y,     Scalar * grad_z,
      Scalar inv_vol) const
    {
      const int K_F_XX = Fields::K_F_XX ;
      const int K_F_YY = Fields::K_F_YY ;
      const int K_F_ZZ = Fields::K_F_ZZ ;
      const int K_F_XY = Fields::K_F_XY ;
      const int K_F_YZ = Fields::K_F_YZ ;
      const int K_F_ZX = Fields::K_F_ZX ;
      const int K_F_YX = Fields::K_F_YX ;
      const int K_F_ZY = Fields::K_F_ZY ;
      const int K_F_XZ = Fields::K_F_XZ ;

      vel_grad(ielem, K_F_XX) = inv_vol * dot8<Scalar,execution_space>( vx , grad_x );
      vel_grad(ielem, K_F_YX) = inv_vol * dot8<Scalar,execution_space>( vy , grad_x );
      vel_grad(ielem, K_F_ZX) = inv_vol * dot8<Scalar,execution_space>( vz , grad_x );

      vel_grad(ielem, K_F_XY) = inv_vol * dot8<Scalar,execution_space>( vx , grad_y );
      vel_grad(ielem, K_F_YY) = inv_vol * dot8<Scalar,execution_space>( vy , grad_y );
      vel_grad(ielem, K_F_ZY) = inv_vol * dot8<Scalar,execution_space>( vz , grad_y );

      vel_grad(ielem, K_F_XZ) = inv_vol * dot8<Scalar,execution_space>( vx , grad_z );
      vel_grad(ielem, K_F_YZ) = inv_vol * dot8<Scalar,execution_space>( vy , grad_z );
      vel_grad(ielem, K_F_ZZ) = inv_vol * dot8<Scalar,execution_space>( vz , grad_z );
    }

  //--------------------------------------------------------------------------
  // Functor operator() which calls the three member functions.


  KOKKOS_INLINE_FUNCTION
  void operator()( int ielem )const
  {
    const int X = 0 ;
    const int Y = 1 ;
    const int Z = 2 ;
    const Scalar dt_scale = -0.5 * dt();

    //  declare and reuse local data for frequently accessed data to
    //  reduce global memory reads and writes.

    Scalar      x[8],      y[8],      z[8];
    Scalar     vx[8],     vy[8],     vz[8];
    Scalar grad_x[8], grad_y[8], grad_z[8];

    // Read global velocity once and use many times
    // via local registers / L1 cache.
    //  store the velocity information in local memory before using,
    //  so it can be returned for other functions to use

    // Read global coordinates and velocity once and use many times
    // via local registers / L1 cache.
    // load X coordinate information and move by half time step

    for ( int i = 0 ; i < ElemNodeCount ; ++i ) {
      const int n = elem_node_connectivity( ielem , i );

      vx[i] = velocity( n , X , current_state );
      vy[i] = velocity( n , Y , current_state );
      vz[i] = velocity( n , Z , current_state );

      x[i]  = model_coords( n , X ) +
              displacement( n , X , current_state ) +
              dt_scale * vx[i];

      y[i]  = model_coords( n , Y ) +
              displacement( n , Y , current_state ) +
              dt_scale * vy[i];

      z[i]  = model_coords( n , Z ) +
              displacement( n , Z , current_state ) +
              dt_scale * vz[i];
    }

    comp_grad<Scalar,execution_space>( x, y, z, grad_x, grad_y, grad_z);

    //  Calculate hexahedral volume from x model_coords and gradient information

    const Scalar inv_vol = 1.0 / dot8<Scalar,execution_space>( x , grad_x );

    v_grad(ielem, vx, vy, vz, grad_x, grad_y, grad_z, inv_vol);
  }

  static void apply( const Fields & fields ,
                     const int arg_current_state ,
                     const int arg_previous_state )
  {
    grad op( fields, arg_current_state , arg_previous_state );
    Kokkos::parallel_for( fields.num_elements , op );
  }
};

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

template<typename Scalar, class DeviceType >
struct decomp_rotate
{
  typedef DeviceType execution_space ;

  typedef Explicit::Fields< Scalar , execution_space >  Fields ;

  static const int ElemNodeCount = Fields::ElemNodeCount ;

  static const int K_F_XX = Fields::K_F_XX ;
  static const int K_F_YY = Fields::K_F_YY ;
  static const int K_F_ZZ = Fields::K_F_ZZ ;
  static const int K_F_XY = Fields::K_F_XY ;
  static const int K_F_YZ = Fields::K_F_YZ ;
  static const int K_F_ZX = Fields::K_F_ZX ;
  static const int K_F_YX = Fields::K_F_YX ;
  static const int K_F_ZY = Fields::K_F_ZY ;
  static const int K_F_XZ = Fields::K_F_XZ ;

  static const int K_S_XX = Fields::K_S_XX ;
  static const int K_S_YY = Fields::K_S_YY ;
  static const int K_S_ZZ = Fields::K_S_ZZ ;
  static const int K_S_XY = Fields::K_S_XY ;
  static const int K_S_YZ = Fields::K_S_YZ ;
  static const int K_S_ZX = Fields::K_S_ZX ;
  static const int K_S_YX = Fields::K_S_YX ;
  static const int K_S_ZY = Fields::K_S_ZY ;
  static const int K_S_XZ = Fields::K_S_XZ ;

  static const int K_V_XY = Fields::K_V_XY ;
  static const int K_V_YZ = Fields::K_V_YZ ;
  static const int K_V_ZX = Fields::K_V_ZX ;

  // Global arrays used by this functor.

  const typename Fields::elem_tensor_state_type     rotation ;
  const typename Fields::elem_tensor_type           vel_grad ;
  const typename Fields::elem_sym_tensor_type       stretch ;
  const typename Fields::elem_sym_tensor_type       rot_stretch ;
  const typename Fields::scalar_type                dt_value ;

  const int  current_state;
  const int  previous_state;

  decomp_rotate( const Fields & mesh_fields ,
                 const int arg_current_state,
                 const int arg_previous_state)
    : rotation(    mesh_fields.rotation )
    , vel_grad(    mesh_fields.vel_grad )
    , stretch(     mesh_fields.stretch )
    , rot_stretch( mesh_fields.rot_stretch )
    , dt_value(    mesh_fields.dt)
    , current_state( arg_current_state)
    , previous_state(arg_previous_state)
    {}

  static void apply( const Fields & mesh_fields ,
                     const int arg_current_state ,
                     const int arg_previous_state )
  {
    decomp_rotate op( mesh_fields , arg_current_state , arg_previous_state );
    Kokkos::parallel_for( mesh_fields.num_elements , op );
  }


  KOKKOS_INLINE_FUNCTION
  void additive_decomp(int ielem, Scalar * v_gr, Scalar * str_ten) const
  {
    //  In addition to calculating stretching_tensor,
    //  use this as an opportunity to load global
    //  variables into a local space

    for ( int i = 0 ; i < 9 ; ++i ) {
      v_gr[i] = vel_grad( ielem , i );
    }

    //
    //  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]);
  }

  KOKKOS_INLINE_FUNCTION
  void polar_decomp(int ielem, Scalar * v_gr, Scalar * str_ten, Scalar * str, Scalar * vort, Scalar * rot_old, Scalar * rot_new)const
  {
    const Scalar dt = dt_value();
    const Scalar dt_half = 0.5 * dt;

    //  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.
    for ( int i = 0 ; i < 6 ; ++i ) {
      str[i] = stretch(ielem, i);
    }

    Scalar z1 = str_ten[K_S_XY] * str[K_S_ZX] -
                str_ten[K_S_ZX] * str[K_S_XY] +
                str_ten[K_S_YY] * str[K_S_YZ] -
                str_ten[K_S_YZ] * str[K_S_YY] +
                str_ten[K_S_YZ] * str[K_S_ZZ] -
                str_ten[K_S_ZZ] * str[K_S_YZ];

    Scalar z2 = str_ten[K_S_ZX] * str[K_S_XX] -
                str_ten[K_S_XX] * str[K_S_ZX] +
                str_ten[K_S_YZ] * str[K_S_XY] -
                str_ten[K_S_XY] * str[K_S_YZ] +
                str_ten[K_S_ZZ] * str[K_S_ZX] -
                str_ten[K_S_ZX] * str[K_S_ZZ];

    Scalar z3 = str_ten[K_S_XX] * str[K_S_XY] -
                str_ten[K_S_XY] * str[K_S_XX] +
                str_ten[K_S_XY] * str[K_S_YY] -
                str_ten[K_S_YY] * str[K_S_XY] +
                str_ten[K_S_ZX] * str[K_S_YZ] -
                str_ten[K_S_YZ] * str[K_S_ZX];

  //   forward elimination
    const Scalar a1inv = 1.0 / (str[K_S_YY] + str[K_S_ZZ]);

    const Scalar a4BYa1 = -1 * str[K_S_XY] * a1inv;

    const Scalar a2inv = 1.0 / (str[K_S_ZZ] + str[K_S_XX] + str[K_S_XY] * a4BYa1);

    const Scalar a5 =  -str[K_S_YZ] + str[K_S_ZX] * a4BYa1;

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

  //   backward substitution -
    z3 /= (str[K_S_XX] + str[K_S_YY] + str[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.

    for ( int i = 0 ; i < 9 ; ++i ) {
      rot_old[i] = rotation(ielem, i, previous_state);
    }

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


  //  2) solve for new rotation tensor via gauss elimination.
  //   forward elimination -
    Scalar a12 = - dt_half * z3;
    Scalar a13 =   dt_half * z2;
    Scalar b32 = - dt_half * z1;
    Scalar a22inv = 1.0 / (1.0 + a12 * a12);

    Scalar a13a12 = a13*a12;
    Scalar 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 Scalar a33inv = 1.0 / (1.0 + a13 * a13 + a23 * b32);

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

    for ( int i = 0 ; i < 9 ; ++i ) {
      rotation(ielem, i, current_state) = rot_new[i] ;
    }

  //   update stretch tensor in the new configuration -
    const Scalar a1 = str_ten[K_S_XY] + vort[K_V_XY];
    const Scalar a2 = str_ten[K_S_YZ] + vort[K_V_YZ];
    const Scalar a3 = str_ten[K_S_ZX] + vort[K_V_ZX];
    const Scalar b1 = str_ten[K_S_ZX] - vort[K_V_ZX];
    const Scalar b2 = str_ten[K_S_XY] - vort[K_V_XY];
    const Scalar b3 = str_ten[K_S_YZ] - vort[K_V_YZ];

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

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

  }


  KOKKOS_INLINE_FUNCTION
  void rotate_tensor(int ielem, Scalar * str_ten, Scalar * str, Scalar * rot_new)const {

    Scalar t[9];
    Scalar rot_str[6]; // Rotated stretch

    t[0] = str_ten[K_S_XX]*rot_new[K_F_XX] +
           str_ten[K_S_XY]*rot_new[K_F_YX] +
           str_ten[K_S_XZ]*rot_new[K_F_ZX];

    t[1] = str_ten[K_S_YX]*rot_new[K_F_XX] +
           str_ten[K_S_YY]*rot_new[K_F_YX] +
           str_ten[K_S_YZ]*rot_new[K_F_ZX];

    t[2] = str_ten[K_S_ZX]*rot_new[K_F_XX] +
           str_ten[K_S_ZY]*rot_new[K_F_YX] +
           str_ten[K_S_ZZ]*rot_new[K_F_ZX];

    t[3] = str_ten[K_S_XX]*rot_new[K_F_XY] +
           str_ten[K_S_XY]*rot_new[K_F_YY] +
           str_ten[K_S_XZ]*rot_new[K_F_ZY];

    t[4] = str_ten[K_S_YX]*rot_new[K_F_XY] +
           str_ten[K_S_YY]*rot_new[K_F_YY] +
           str_ten[K_S_YZ]*rot_new[K_F_ZY];

    t[5] = str_ten[K_S_ZX]*rot_new[K_F_XY] +
           str_ten[K_S_ZY]*rot_new[K_F_YY] +
           str_ten[K_S_ZZ]*rot_new[K_F_ZY];

    t[6] = str_ten[K_S_XX]*rot_new[K_F_XZ] +
           str_ten[K_S_XY]*rot_new[K_F_YZ] +
           str_ten[K_S_XZ]*rot_new[K_F_ZZ];

    t[7] = str_ten[K_S_YX]*rot_new[K_F_XZ] +
           str_ten[K_S_YY]*rot_new[K_F_YZ] +
           str_ten[K_S_YZ]*rot_new[K_F_ZZ];

    t[8] = str_ten[K_S_ZX]*rot_new[K_F_XZ] +
           str_ten[K_S_ZY]*rot_new[K_F_YZ] +
           str_ten[K_S_ZZ]*rot_new[K_F_ZZ];


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

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

    for ( int i = 0 ; i < 6 ; ++i ) {
      rot_stretch(ielem, i) = rot_str[i] ;
    }

    for ( int i = 0 ; i < 6 ; ++i ) {
      stretch(ielem, i) = str[i] ;
    }
  }

  KOKKOS_INLINE_FUNCTION
  void operator()( int ielem )const {

    //   Local scratch space to avoid multiple
    //   accesses to global memory.
    Scalar str_ten[6]; // Stretching tensor
    Scalar str[6];     // Stretch
    Scalar rot_old[9]; // Rotation old
    Scalar rot_new[9]; // Rotation new
    Scalar vort[3];    // Vorticity
    Scalar v_gr[9];    // Velocity gradient

    additive_decomp(ielem, v_gr, str_ten);

    polar_decomp(ielem, v_gr, str_ten, str, vort, rot_old, rot_new);

    rotate_tensor(ielem, str_ten, str, rot_new);
  }
};

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

template<typename Scalar, class DeviceType >
struct internal_force
{
  typedef DeviceType execution_space ;

  typedef Explicit::Fields< Scalar , execution_space >  Fields ;

  static const int ElemNodeCount = Fields::ElemNodeCount ;

  static const int K_F_XX = Fields::K_F_XX ;
  static const int K_F_YY = Fields::K_F_YY ;
  static const int K_F_ZZ = Fields::K_F_ZZ ;
  static const int K_F_XY = Fields::K_F_XY ;
  static const int K_F_YZ = Fields::K_F_YZ ;
  static const int K_F_ZX = Fields::K_F_ZX ;
  static const int K_F_YX = Fields::K_F_YX ;
  static const int K_F_ZY = Fields::K_F_ZY ;
  static const int K_F_XZ = Fields::K_F_XZ ;

  static const int K_S_XX = Fields::K_S_XX ;
  static const int K_S_YY = Fields::K_S_YY ;
  static const int K_S_ZZ = Fields::K_S_ZZ ;
  static const int K_S_XY = Fields::K_S_XY ;
  static const int K_S_YZ = Fields::K_S_YZ ;
  static const int K_S_ZX = Fields::K_S_ZX ;
  static const int K_S_YX = Fields::K_S_YX ;
  static const int K_S_ZY = Fields::K_S_ZY ;
  static const int K_S_XZ = Fields::K_S_XZ ;

  //--------------------------------------------------------------------------
  // Reduction:

  typedef Scalar value_type;

  KOKKOS_INLINE_FUNCTION
  static void init(value_type &update) {
    update = 1.0e32;
  }

  KOKKOS_INLINE_FUNCTION
  static void join( volatile value_type & update,
                    const volatile value_type & source )
  {
    update = update < source ? update : source;
  }

  // Final serial processing of reduction value:
  KOKKOS_INLINE_FUNCTION
  void final( value_type & result ) const
  {
    prev_dt() = dt() ;
    dt() = result ;
  };

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

  // Global arrays used by this functor.

  const typename Fields::elem_node_ids_type      elem_node_connectivity ;
  const typename Fields::node_coords_type        model_coords ;
  const typename Fields::scalar_type             dt ;
  const typename Fields::scalar_type             prev_dt ;
  const typename Fields::geom_state_array_type   displacement ;
  const typename Fields::geom_state_array_type   velocity ;
  const typename Fields::array_type              elem_mass ;
  const typename Fields::array_type              internal_energy ;
  const typename Fields::elem_sym_tensor_type    stress_new ;
  const typename Fields::elem_node_geom_type     element_force ;
  const typename Fields::elem_tensor_state_type  rotation ;
  const typename Fields::elem_sym_tensor_type    rot_stretch ;

  const Scalar     two_mu;
  const Scalar     bulk_modulus;
  const Scalar     lin_bulk_visc;
  const Scalar     quad_bulk_visc;
  const Scalar     user_dt;
  const int        current_state;

  internal_force( const Fields & mesh_fields,
                  const Scalar arg_user_dt,
                  const int arg_current_state )
    : elem_node_connectivity( mesh_fields.elem_node_connectivity )
    , model_coords(           mesh_fields.model_coords )
    , dt(                     mesh_fields.dt )
    , prev_dt(                mesh_fields.prev_dt )
    , displacement(           mesh_fields.displacement )
    , velocity(               mesh_fields.velocity )
    , elem_mass(              mesh_fields.elem_mass )
    , internal_energy(        mesh_fields.internal_energy )
    , stress_new(             mesh_fields.stress_new )
    , element_force(          mesh_fields.element_force )
    , rotation(               mesh_fields.rotation )
    , rot_stretch(            mesh_fields.rot_stretch )
    , two_mu(                 mesh_fields.two_mu )
    , bulk_modulus(           mesh_fields.bulk_modulus )
    , lin_bulk_visc(          mesh_fields.lin_bulk_visc )
    , quad_bulk_visc(         mesh_fields.quad_bulk_visc )
    , user_dt(       arg_user_dt )
    , current_state( arg_current_state )
  {}

  static void apply( const Fields & mesh_fields ,
                     const Scalar arg_user_dt,
                     const int arg_current_state )
  {
    internal_force  op_force( mesh_fields , arg_user_dt , arg_current_state );

    Kokkos::parallel_reduce( mesh_fields.num_elements, op_force );
  }

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

  KOKKOS_INLINE_FUNCTION
  void rotate_tensor_backward(int ielem ,
    const Scalar * const s_n ,
    Scalar * const rot_stress )const
  {
    const int rot_state = current_state ; // 1 ;

    //   t : temporary variables
    //   s_n : stress_new in local memory space
    //   r_n : rotation_new in local memory space
    Scalar t[9], r_n[9];

    r_n[0] = rotation(ielem, 0, rot_state );
    r_n[1] = rotation(ielem, 1, rot_state );
    r_n[2] = rotation(ielem, 2, rot_state );
    r_n[3] = rotation(ielem, 3, rot_state );
    r_n[4] = rotation(ielem, 4, rot_state );
    r_n[5] = rotation(ielem, 5, rot_state );
    r_n[6] = rotation(ielem, 6, rot_state );
    r_n[7] = rotation(ielem, 7, rot_state );
    r_n[8] = rotation(ielem, 8, rot_state );

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

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

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

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

  KOKKOS_INLINE_FUNCTION
  void comp_force(int ielem,
     const Scalar * const vx ,
     const Scalar * const vy ,
     const Scalar * const vz ,
     const Scalar * const grad_x ,
     const Scalar * const grad_y ,
     const Scalar * const grad_z ,
     Scalar * total_stress12th ) const
  {
    Scalar internal_energy_inc = 0 ;

    for(int inode = 0; inode < 8; ++inode) {

      const Scalar fx =
        total_stress12th[K_S_XX] * grad_x[inode] +
        total_stress12th[K_S_XY] * grad_y[inode] +
        total_stress12th[K_S_XZ] * grad_z[inode] ;

      element_force(ielem, 0, inode) = fx ;

      const Scalar fy =
        total_stress12th[K_S_YX] * grad_x[inode] +
        total_stress12th[K_S_YY] * grad_y[inode] +
        total_stress12th[K_S_YZ] * grad_z[inode] ;

      element_force(ielem, 1, inode) = fy ;

      const Scalar fz =
        total_stress12th[K_S_ZX] * grad_x[inode] +
        total_stress12th[K_S_ZY] * grad_y[inode] +
        total_stress12th[K_S_ZZ] * grad_z[inode] ;

      element_force(ielem, 2, inode) = fz ;

      internal_energy_inc +=
        fx * vx[inode] +
        fy * vy[inode] +
        fz * vz[inode] ;
    }

    internal_energy(ielem) = internal_energy_inc ;
  }

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

  KOKKOS_INLINE_FUNCTION
  void get_stress(int ielem , Scalar * const s_n ) const
    {
      const int kxx = 0;
      const int kyy = 1;
      const int kzz = 2;
      const int kxy = 3;
      const int kyz = 4;
      const int kzx = 5;

      const Scalar e = (rot_stretch(ielem,kxx)+rot_stretch(ielem,kyy)+rot_stretch(ielem,kzz))/3.0;

      s_n[kxx] = stress_new(ielem,kxx) += dt() * (two_mu * (rot_stretch(ielem,kxx)-e)+3*bulk_modulus*e);
      s_n[kyy] = stress_new(ielem,kyy) += dt() * (two_mu * (rot_stretch(ielem,kyy)-e)+3*bulk_modulus*e);
      s_n[kzz] = stress_new(ielem,kzz) += dt() * (two_mu * (rot_stretch(ielem,kzz)-e)+3*bulk_modulus*e);

      s_n[kxy] = stress_new(ielem,kxy) += dt() * two_mu * rot_stretch(ielem,kxy);
      s_n[kyz] = stress_new(ielem,kyz) += dt() * two_mu * rot_stretch(ielem,kyz);
      s_n[kzx] = stress_new(ielem,kzx) += dt() * two_mu * rot_stretch(ielem,kzx);
    }

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


  KOKKOS_INLINE_FUNCTION
  void operator()( int ielem, value_type & update )const
  {
    const Scalar ONE12TH = 1.0 / 12.0 ;

    Scalar x[8], y[8], z[8] ;
    Scalar vx[8], vy[8], vz[8];
    Scalar grad_x[8], grad_y[8], grad_z[8];

    // Position and velocity:

    for ( int i = 0 ; i < ElemNodeCount ; ++i ) {
      const int n = elem_node_connectivity(ielem,i);

      x[i] = model_coords(n, 0) + displacement(n, 0, current_state) ;
      y[i] = model_coords(n, 1) + displacement(n, 1, current_state) ;
      z[i] = model_coords(n, 2) + displacement(n, 2, current_state) ;

      vx[i] = velocity(n, 0, current_state);
      vy[i] = velocity(n, 1, current_state);
      vz[i] = velocity(n, 2, current_state);
    }

    // Gradient:

    comp_grad<Scalar,execution_space>( x , y , z , grad_x , grad_y , grad_z );


    const Scalar mid_vol = dot8<Scalar,execution_space>( x , grad_x );

    const Scalar shr = two_mu ;
    const Scalar dil = bulk_modulus + ((2.0*shr)/3.0);

    const Scalar aspect = 6.0 * mid_vol /
                          ( dot8<Scalar,execution_space>( grad_x , grad_x ) +
                            dot8<Scalar,execution_space>( grad_y , grad_y ) +
                            dot8<Scalar,execution_space>( grad_z , grad_z ) );

    const Scalar dtrial = std::sqrt(elem_mass(ielem) * aspect / dil);
    const Scalar traced = (rot_stretch(ielem, 0) + rot_stretch(ielem, 1) + rot_stretch(ielem, 2));

    const Scalar eps = traced < 0 ? (lin_bulk_visc - quad_bulk_visc * traced * dtrial) : lin_bulk_visc ;

    const Scalar bulkq = eps * dil * dtrial * traced;

    Scalar cur_time_step = dtrial * ( std::sqrt( 1.0 + eps * eps) - eps);

    // force fixed time step if input

    cur_time_step = user_dt > 0 ? user_dt : cur_time_step;

    update = update < cur_time_step ? update : cur_time_step;


    Scalar s_n[ 6 ];

    get_stress( ielem, s_n );

    Scalar total_stress12th[6];

    // Get rotated stress:

    rotate_tensor_backward(ielem, s_n , total_stress12th );

    total_stress12th[0] = ONE12TH*( total_stress12th[ 0 ] + bulkq );
    total_stress12th[1] = ONE12TH*( total_stress12th[ 1 ] + bulkq );
    total_stress12th[2] = ONE12TH*( total_stress12th[ 2 ] + bulkq );
    total_stress12th[3] = ONE12TH*( total_stress12th[ 3 ] );
    total_stress12th[4] = ONE12TH*( total_stress12th[ 4 ] );
    total_stress12th[5] = ONE12TH*( total_stress12th[ 5 ] );

    comp_force(ielem, vx, vy, vz,
                      grad_x, grad_y, grad_z, total_stress12th);
  }
};

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

template<typename Scalar, class DeviceType >
struct nodal_step
{
  typedef DeviceType     execution_space ;
  typedef typename execution_space::size_type  size_type;

  typedef Explicit::Fields< Scalar , execution_space >  Fields ;

  const typename Fields::scalar_type            dt ;
  const typename Fields::scalar_type            prev_dt ;
  const typename Fields::node_elem_ids_type     node_elem_connectivity ;
  const typename Fields::node_coords_type       model_coords ;
  const typename Fields::array_type             nodal_mass ;
  const typename Fields::geom_state_array_type  displacement ;
  const typename Fields::geom_state_array_type  velocity ;
  const typename Fields::geom_array_type        acceleration ;
  const typename Fields::geom_array_type        internal_force ;
  const typename Fields::elem_node_geom_type    element_force ;

  const Scalar   x_bc;
  const int      current_state;
  const int      next_state;


  nodal_step( const Fields  & mesh_fields ,
              const Scalar    arg_x_bc,
              const int       arg_current_state,
              const int       arg_next_state)
   : dt(       mesh_fields.dt )
   , prev_dt(  mesh_fields.prev_dt )
   , node_elem_connectivity( mesh_fields.node_elem_connectivity )
   , model_coords(   mesh_fields.model_coords )
   , nodal_mass(     mesh_fields.nodal_mass )
   , displacement(   mesh_fields.displacement )
   , velocity(       mesh_fields.velocity )
   , acceleration(   mesh_fields.acceleration )
   , internal_force( mesh_fields.internal_force )
   , element_force(  mesh_fields.element_force )
   , x_bc(          arg_x_bc )
   , current_state( arg_current_state )
   , next_state(    arg_next_state )
   {
        //std::cout << "finish_step dt: " << dt << std::endl;
        //std::cout << "finish_step prev_dt: " << prev_dt << std::endl;
   }

  static void apply( const Fields  & mesh_fields ,
                     const Scalar    arg_x_bc ,
                     const int       arg_current_state ,
                     const int       arg_next_state )
  {
    nodal_step op( mesh_fields, arg_x_bc, arg_current_state, arg_next_state );

    // Only update the owned nodes:

    Kokkos::parallel_for( mesh_fields.num_nodes_owned , op );
  }

  KOKKOS_INLINE_FUNCTION
  void operator()(int inode) const
    {
      // Getting count as per 'CSR-like' data structure
      const int begin = node_elem_connectivity.row_map[inode];
      const int end   = node_elem_connectivity.row_map[inode+1];

      double local_force[] = {0.0, 0.0, 0.0};

      // Gather-sum internal force from
      // each element that a node is attached to.

      for ( int i = begin; i < end ; ++i ){

        //  node_elem_offset is a cumulative structure, so
        //  node_elem_offset(inode) should be the index where
        //  a particular row's elem_IDs begin
        const int nelem = node_elem_connectivity.entries( i, 0);

        //  find the row in an element's stiffness matrix
        //  that corresponds to inode
        const int elem_node_index = node_elem_connectivity.entries( i, 1);

        local_force[0] += element_force(nelem, 0, elem_node_index);
        local_force[1] += element_force(nelem, 1, elem_node_index);
        local_force[2] += element_force(nelem, 2, elem_node_index);
      }

      internal_force(inode, 0) = local_force[0];
      internal_force(inode, 1) = local_force[1];
      internal_force(inode, 2) = local_force[2];

      // Acceleration:

      Scalar v_new[3];
      Scalar a_current[3];

      const Scalar tol = 1.0e-7;

      // If not on the boundary then: a = F / m
      if ( tol < fabs(model_coords(inode,0)-x_bc) ) {

        const Scalar m = nodal_mass( inode );

        acceleration(inode,0) = a_current[0] = -local_force[0] / m ;
        acceleration(inode,1) = a_current[1] = -local_force[1] / m ;
        acceleration(inode,2) = a_current[2] = -local_force[2] / m ;
      }
      else { //enforce fixed BC
        acceleration(inode,0) = a_current[0] = 0;
        acceleration(inode,1) = a_current[1] = 0;
        acceleration(inode,2) = a_current[2] = 0;
      }

      // Central difference time integration:

      const Scalar dt_disp = dt() ;
      const Scalar dt_vel = ( dt() + prev_dt() ) / 2.0 ;

      velocity(inode,0,next_state) = v_new[0] =
        velocity(inode,0,current_state) + dt_vel * a_current[0];

      velocity(inode,1,next_state) = v_new[1] =
        velocity(inode,1,current_state) + dt_vel * a_current[1];

      velocity(inode,2,next_state) = v_new[2] =
        velocity(inode,2,current_state) + dt_vel * a_current[2];

      displacement(inode,0,next_state) =
        displacement(inode,0,current_state) + dt_disp * v_new[0];

      displacement(inode,1,next_state) =
        displacement(inode,1,current_state) + dt_disp * v_new[1];

      displacement(inode,2,next_state) =
        displacement(inode,2,current_state) + dt_disp * v_new[2];
    }
};

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

template< typename Scalar , class DeviceType >
struct pack_state
{
  typedef DeviceType     execution_space ;
  typedef typename execution_space::size_type  size_type ;

  typedef Explicit::Fields< Scalar , execution_space >  Fields ;

  typedef typename Fields::geom_state_array_type::value_type  value_type ;
  typedef Kokkos::View< value_type* , execution_space >     buffer_type ;

  static const unsigned value_count = 6 ;

  const typename Fields::geom_state_array_type  displacement ;
  const typename Fields::geom_state_array_type  velocity ;
  const buffer_type  output ;
  const size_type    inode_base ;
  const size_type    state_next ;

  pack_state( const buffer_type & arg_output ,
              const Fields      & mesh_fields ,
              const size_type     arg_begin ,
              const size_type     arg_state )
   : displacement( mesh_fields.displacement )
   , velocity(     mesh_fields.velocity )
   , output(       arg_output )
   , inode_base(   arg_begin )
   , state_next(   arg_state )
   {}

  static void apply( const buffer_type & arg_output ,
                     const size_type     arg_begin ,
                     const size_type     arg_count ,
                     const Fields      & mesh_fields ,
                     const size_type     arg_state )
  {
    pack_state op( arg_output , mesh_fields , arg_begin , arg_state );

    Kokkos::parallel_for( arg_count , op );
  }

  KOKKOS_INLINE_FUNCTION
  void operator()( const size_type i ) const
  {
    const size_type inode = inode_base + i ;

    size_type j = i * value_count ;

    output[j++] = displacement( inode , 0 , state_next );
    output[j++] = displacement( inode , 1 , state_next );
    output[j++] = displacement( inode , 2 , state_next );
    output[j++] = velocity( inode , 0 , state_next );
    output[j++] = velocity( inode , 1 , state_next );
    output[j++] = velocity( inode , 2 , state_next );
  }
};

template< typename Scalar , class DeviceType >
struct unpack_state
{
  typedef DeviceType     execution_space ;
  typedef typename execution_space::size_type  size_type ;

  typedef Explicit::Fields< Scalar , execution_space >  Fields ;

  typedef typename Fields::geom_state_array_type::value_type  value_type ;
  typedef Kokkos::View< value_type* , execution_space >     buffer_type ;

  static const unsigned value_count = 6 ;

  const typename Fields::geom_state_array_type  displacement ;
  const typename Fields::geom_state_array_type  velocity ;
  const buffer_type  input ;
  const size_type    inode_base ;
  const size_type    state_next ;

  unpack_state( const buffer_type & arg_input ,
                const Fields      & mesh_fields ,
                const size_type     arg_begin ,
                const size_type     arg_state )
   : displacement( mesh_fields.displacement )
   , velocity(     mesh_fields.velocity )
   , input(        arg_input )
   , inode_base(   arg_begin )
   , state_next(   arg_state )
   {}

  static void apply( const Fields      & mesh_fields ,
                     const size_type     arg_state ,
                     const buffer_type & arg_input ,
                     const size_type     arg_begin ,
                     const size_type     arg_count )
  {
    unpack_state op( arg_input , mesh_fields , arg_begin , arg_state );

    Kokkos::parallel_for( arg_count , op );
  }

  KOKKOS_INLINE_FUNCTION
  void operator()( const size_type i ) const
  {
    const size_type inode = inode_base + i ;

    size_type j = i * value_count ;

    displacement( inode , 0 , state_next ) = input[j++] ;
    displacement( inode , 1 , state_next ) = input[j++] ;
    displacement( inode , 2 , state_next ) = input[j++] ;
    velocity( inode , 0 , state_next ) = input[j++] ;
    velocity( inode , 1 , state_next ) = input[j++] ;
    velocity( inode , 2 , state_next ) = input[j++] ;
  }
};

} /* namespace Explicit */

#endif /* #ifndef KOKKOS_EXPLICITFUNCTORS_HPP */