20 #ifndef GEOSX_PHYSICSSOLVERS_CONTACT_EXPLICITMPM_HPP_
21 #define GEOSX_PHYSICSSOLVERS_CONTACT_EXPLICITMPM_HPP_
23 #include "constitutive/solid/SolidUtilities.hpp"
25 #include "physicsSolvers/solidMechanics/MPMSolverFields.hpp"
30 namespace solidMechanicsMPMKernels
35 void polarDecomposition(
real64 (& R)[3][3],
36 real64 const (&matrix)[3][3] )
39 LvArray::tensorOps::copy< 3, 3 >( R, matrix );
40 real64 RInverse[3][3] = { {0} },
41 RInverseTranspose[3][3] = { {0} },
42 RRTMinusI[3][3] = { {0} };
46 real64 tolerance = 10 * LvArray::NumericLimits< real64 >::epsilon;
48 while( errorSquared > tolerance * tolerance && iter < 100 )
54 LvArray::tensorOps::internal::SquareMatrixOps< 3 >::invert( RInverse, R );
55 LvArray::tensorOps::transpose< 3, 3 >( RInverseTranspose, RInverse );
56 LvArray::tensorOps::add< 3, 3 >( R, RInverseTranspose );
57 LvArray::tensorOps::scale< 3, 3 >( R, 0.5 );
61 LvArray::tensorOps::copy< 3, 3 >( copyR, R );
62 LvArray::tensorOps::Rij_eq_AikBjk< 3, 3, 3 >( RRTMinusI, R, copyR );
63 LvArray::tensorOps::addIdentity< 3 >( RRTMinusI, -1.0 );
64 for( std::ptrdiff_t i = 0; i < 3; i++ )
66 for( std::ptrdiff_t j = 0; j < 3; j++ )
68 errorSquared += RRTMinusI[i][j] * RRTMinusI[i][j];
74 GEOS_LOG_RANK(
"Polar decomposition did not converge in 100 iterations!" );
94 template<
typename POLICY,
typename CONSTITUTIVE_WRAPPER >
96 CONSTITUTIVE_WRAPPER
const & constitutiveWrapper,
114 #if defined(GEOSX_USE_CUDA)
115 LvArray::tensorOps::copy< 6 >( oldStress[p][0], particleStress[p] );
119 real64 strainIncrement[6];
120 strainIncrement[0] = velocityGradient[p][0][0] * dt;
121 strainIncrement[1] = velocityGradient[p][1][1] * dt;
122 strainIncrement[2] = velocityGradient[p][2][2] * dt;
123 strainIncrement[3] = (velocityGradient[p][1][2] + velocityGradient[p][2][1]) * dt;
124 strainIncrement[4] = (velocityGradient[p][0][2] + velocityGradient[p][2][0]) * dt;
125 strainIncrement[5] = (velocityGradient[p][0][1] + velocityGradient[p][1][0]) * dt;
128 real64 fOld[3][3] = { {0} };
129 real64 fNew[3][3] = { {0} };
130 LvArray::tensorOps::copy< 3, 3 >( fNew, deformationGradient[p] );
131 LvArray::tensorOps::copy< 3, 3 >( fOld, deformationGradient[p] );
132 LvArray::tensorOps::scaledAdd< 3, 3 >( fOld, fDot[p], -dt );
135 real64 rotBeginning[3][3] = { {0} };
136 real64 rotEnd[3][3] = { {0} };
137 polarDecomposition( rotBeginning, fOld );
138 polarDecomposition( rotEnd, fNew );
142 constitutive::SolidUtilities::hypoUpdate2_StressOnly( constitutiveWrapper,
152 LvArray::tensorOps::copy< 6 >( particleStress[p], stress );
155 constitutiveWrapper.saveConvergedState( p, 0 );
#define GEOS_HOST_DEVICE
Marks a host-device function.
#define GEOS_LOG_RANK(msg)
Log a message to the rank output stream.
double real64
64-bit floating point type.
LvArray::SortedArrayView< T, localIndex, LvArray::ChaiBuffer > SortedArrayView
A sorted array view of local indices.
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
GEOSX_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
A struct to update particle stresses.
static void launch(SortedArrayView< localIndex const > const indices, CONSTITUTIVE_WRAPPER const &constitutiveWrapper, real64 dt, arrayView3d< real64 const > const deformationGradient, arrayView3d< real64 const > const fDot, arrayView3d< real64 const > const velocityGradient, arrayView2d< real64 > const particleStress)
Launch the kernel function doing constitutive updates.