19 #ifndef GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_IMPLCITSMALLSTRAINQUASISTATIC_IMPL_HPP_
20 #define GEOS_PHYSICSSOLVERS_SOLIDMECHANICS_KERNELS_IMPLCITSMALLSTRAINQUASISTATIC_IMPL_HPP_
27 namespace solidMechanicsLagrangianFEMKernels
31 template<
typename SUBREGION_TYPE,
32 typename CONSTITUTIVE_TYPE,
39 SUBREGION_TYPE
const & elementSubRegion,
40 FE_TYPE
const & finiteElementSpace,
41 CONSTITUTIVE_TYPE & inputConstitutiveType,
47 real64 const (&inputGravityVector)[3] ):
54 inputConstitutiveType,
60 m_X( nodeManager.referencePosition()),
61 m_disp( nodeManager.getField< fields::solidMechanics::totalDisplacement >() ),
62 m_uhat( nodeManager.getField< fields::solidMechanics::incrementalDisplacement >() ),
63 m_gravityVector{ inputGravityVector[0], inputGravityVector[1], inputGravityVector[2] },
64 m_density( inputConstitutiveType.getDensity() )
68 template<
typename SUBREGION_TYPE,
69 typename CONSTITUTIVE_TYPE,
77 m_finiteElementSpace.template setup< FE_TYPE >( k, m_meshData, stack.feStack );
79 localIndex const numSupportPoints = m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
81 stack.numRows = 3 * numSupportPoints;
82 stack.numCols = stack.numRows;
85 for(
localIndex a = 0; a < numSupportPoints; ++a )
87 localIndex const localNodeIndex = m_elemsToNodes( k, a );
90 for(
int i = 0; i < numDofPerTestSupportPoint; ++i )
92 #if defined(CALC_FEM_SHAPE_IN_KERNEL)
93 stack.
xLocal[ a ][ i ] = m_X[ localNodeIndex ][ i ];
95 stack.
u_local[ a ][i] = m_disp[ localNodeIndex ][i];
96 stack.
uhat_local[ a ][i] = m_uhat[ localNodeIndex ][i];
97 stack.localRowDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i;
98 stack.localColDofIndex[a*3+i] = m_dofNumber[localNodeIndex]+i;
103 real64 const stabilizationScaling = computeStabilizationScaling( k );
104 m_finiteElementSpace.template addGradGradStabilizationMatrix
105 < FE_TYPE, numDofPerTrialSupportPoint,
true >( stack.feStack,
107 -stabilizationScaling );
110 template<
typename SUBREGION_TYPE,
111 typename CONSTITUTIVE_TYPE,
113 template<
typename STRESS_MODIFIER >
118 StackVariables & stack,
119 STRESS_MODIFIER && stressModifier )
const
121 real64 dNdX[ numNodesPerElem ][ 3 ];
122 real64 const detJxW = m_finiteElementSpace.template getGradN< FE_TYPE >( k, q, stack.xLocal, stack.feStack, dNdX );
124 real64 strainInc[6] = {0};
127 typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness;
129 FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, strainInc );
131 m_constitutiveUpdate.smallStrainUpdate( k, q, m_dt, strainInc, stress, stiffness );
133 stressModifier( stress );
137 stress[i] *= -detJxW;
140 real64 const gravityForce[3] = { m_gravityVector[0] * m_density( k, q )* detJxW,
141 m_gravityVector[1] * m_density( k, q )* detJxW,
142 m_gravityVector[2] * m_density( k, q )* detJxW };
144 real64 N[numNodesPerElem];
145 FE_TYPE::calcN( q, stack.feStack, N );
146 FE_TYPE::plusGradNajAijPlusNaFi( dNdX,
150 reinterpret_cast< real64 (&)[numNodesPerElem][3]
>(stack.localResidual) );
151 real64 const stabilizationScaling = computeStabilizationScaling( k );
152 m_finiteElementSpace.template addEvaluatedGradGradStabilizationVector< FE_TYPE, numDofPerTrialSupportPoint >( stack.feStack,
154 reinterpret_cast< real64 (&)[numNodesPerElem][3]
>(stack.localResidual),
155 -stabilizationScaling );
156 #if !defined( GEOS_USE_HIP )
157 stiffness.template upperBTDB< numNodesPerElem >( dNdX, -detJxW, stack.localJacobian );
159 stiffness.template BTDB< numNodesPerElem >( dNdX, -detJxW, stack.localJacobian );
164 template<
typename SUBREGION_TYPE,
165 typename CONSTITUTIVE_TYPE,
175 #if !defined( GEOS_USE_HIP )
177 CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps::template fillLowerBTDB< numNodesPerElem >( stack.localJacobian );
179 localIndex const numSupportPoints = m_finiteElementSpace.template numSupportPoints< FE_TYPE >( stack.feStack );
182 for(
int localNode = 0; localNode < numSupportPoints; ++localNode )
185 for(
int dim = 0; dim < numDofPerTestSupportPoint; ++dim )
187 localIndex const dof = LvArray::integerConversion< localIndex >( stack.localRowDofIndex[ numDofPerTestSupportPoint * localNode + dim ] - m_dofRankOffset );
188 if( dof < 0 || dof >= m_matrix.numRows() )
190 m_matrix.template addToRowBinarySearchUnsorted< parallelDeviceAtomic >( dof,
191 stack.localRowDofIndex,
192 stack.localJacobian[ numDofPerTestSupportPoint * localNode + dim ],
195 RAJA::atomicAdd< parallelDeviceAtomic >( &m_rhs[ dof ], stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] );
196 maxForce = fmax( maxForce, fabs( stack.localResidual[ numDofPerTestSupportPoint * localNode + dim ] ) );
204 template<
typename SUBREGION_TYPE,
205 typename CONSTITUTIVE_TYPE,
207 template<
typename POLICY,
208 typename KERNEL_TYPE >
212 KERNEL_TYPE
const & kernelComponent )
214 return Base::template kernelLaunch< POLICY, KERNEL_TYPE >( numElems, kernelComponent );
#define GEOS_HOST_DEVICE
Marks a host-device function.
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
#define GEOS_FORCE_INLINE
Marks a function or lambda for inlining.
This class provides an interface to ObjectManagerBase in order to manage edge data.
The FaceManager class provides an interface to ObjectManagerBase in order to manage face data.
The NodeManager class provides an interface to ObjectManagerBase in order to manage node data.
static real64 kernelLaunch(localIndex const numElems, KERNEL_TYPE const &kernelComponent)
Kernel Launcher.
GEOS_HOST_DEVICE void quadraturePointKernel(localIndex const k, localIndex const q, StackVariables &stack, STRESS_MODIFIER &&stressModifier=NoOpFunctors{}) const
Performs a state update at a quadrature point.
ImplicitSmallStrainQuasiStatic(NodeManager const &nodeManager, EdgeManager const &edgeManager, FaceManager const &faceManager, localIndex const targetRegionIndex, SUBREGION_TYPE const &elementSubRegion, FE_TYPE const &finiteElementSpace, CONSTITUTIVE_TYPE &inputConstitutiveType, arrayView1d< globalIndex const > const inputDofNumber, globalIndex const rankOffset, CRSMatrixView< real64, globalIndex const > const inputMatrix, arrayView1d< real64 > const inputRhs, real64 const inputDt, real64 const (&inputGravityVector)[3])
Constructor.
GEOS_HOST_DEVICE real64 complete(localIndex const k, StackVariables &stack) const
Performs the complete phase for the kernel.
GEOS_HOST_DEVICE void setup(localIndex const k, StackVariables &stack) const
Copy global values from primary field to a local stack array.
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer > CRSMatrixView
Alias for CRS Matrix View.
double real64
64-bit floating point type.
GEOSX_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
GEOSX_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Kernel variables allocated on the stack.
real64 uhat_local[numNodesPerElem][numDofPerTrialSupportPoint]
Stack storage for the element local nodal incremental displacement.
real64 u_local[numNodesPerElem][numDofPerTrialSupportPoint]
Stack storage for the element local nodal displacement.