HexExplicitFunctions.hpp
0541d52c17a1lammm-devel

HexExplicitFunctions.hpp
View Options

			/*
			//@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 = (z1a1inv - 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 */

lammps/lib/kokkos/example/multi_fem/HexExplicitFunctions.hpp0541d52c17a1lammm-devel

HexExplicitFunctions.hppView Options

lammps/lib/kokkos/example/multi_fem/HexExplicitFunctions.hpp
0541d52c17a1lammm-devel

HexExplicitFunctions.hpp
View Options