444 lines
19 KiB
C++
444 lines
19 KiB
C++
/*
|
|
//@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 H. Carter Edwards (hcedwar@sandia.gov)
|
|
//
|
|
// ************************************************************************
|
|
//@HEADER
|
|
*/
|
|
|
|
#ifndef KOKKOS_HEXEXPLICITFUNCTIONS_HPP
|
|
#define KOKKOS_HEXEXPLICITFUNCTIONS_HPP
|
|
|
|
#include <math.h>
|
|
|
|
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 */
|
|
|