diff --git a/modules/core/MarmotFiniteElementCore/include/Marmot/MarmotElement.h b/modules/core/MarmotFiniteElementCore/include/Marmot/MarmotElement.h index 14fd4817..86a89e3f 100644 --- a/modules/core/MarmotFiniteElementCore/include/Marmot/MarmotElement.h +++ b/modules/core/MarmotFiniteElementCore/include/Marmot/MarmotElement.h @@ -148,6 +148,26 @@ class MarmotElement { double dT, double& pNewdT ) = 0; + /** + * @brief Perform element computations for explicit time integration. + * @param[in] QTotal Total dof vector. + * @param[in] dQ Incremental dof vector. + * @param[out] Pint Internal force vector. + * @param[in] time Current time. + * @param[in] dT Time step size. + * @param[out] pNewdT Suggested new time step size. + * + * @note Default implementation throws an exception. + */ + virtual void computeYourselfExplicit( const double* QTotal, + const double* dQ, + double* Pint, + const double* time, + double dT, + double& pNewdT ) + { + throw std::invalid_argument( MakeString() << __PRETTY_FUNCTION__ << " not yet implemented" ); + }; /** * @brief Compute contribution from distributed surface loads. * @param[in] loadType Type of load. @@ -194,6 +214,16 @@ class MarmotElement { throw std::invalid_argument( MakeString() << __PRETTY_FUNCTION__ << " not yet implemented" ); }; + /** + * @brief Compute lumped damping matrix. + * @param[out] C damping as diagonal matrix. + * @note Default implementation throws an exception. + */ + virtual void computeLumpedDamping( double* C ) + { + throw std::invalid_argument( MakeString() << __PRETTY_FUNCTION__ << " not yet implemented" ); + }; + /** * @brief Compute consistent inertia matrix. * @param[out] I Inertia matrix. diff --git a/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h b/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h index c20c2ea9..dedd701c 100644 --- a/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h +++ b/modules/core/MarmotFiniteStrainMechanicsCore/include/Marmot/MarmotMaterialFiniteStrain.h @@ -123,8 +123,26 @@ class MarmotMaterialFiniteStrain { * */ virtual void computeStress( ConstitutiveResponse< 3 >& response, AlgorithmicModuli< 3 >& tangents, - const Deformation< 3 >&, - const TimeIncrement& ) const = 0; + const Deformation< 3 >& deformation, + const TimeIncrement& timeIncrement ) const = 0; + + /** + * @brief Explicit version of computeStress for use in explicit time integration schemes. + * @param[inout] response ConstitutiveResponse instance + * @param[in] deformation Deformation instance + * @param[in] timeIncrement TimeIncrement instance + * + * @note The default implementation calls computeStress and ignores the algorithmic tangent. + * @note Derived classes may override this method for efficiency reasons. + */ + virtual void computeStressExplicit( ConstitutiveResponse< 3 >& response, + const Deformation< 3 >& deformation, + const TimeIncrement& timeIncrement ) const + { + AlgorithmicModuli< 3 > tangents; + computeStress( response, tangents, deformation, timeIncrement ); + } + /** * @brief Computes the Kirchhoff stress given the deformation, time increment, and eigen deformation. * @param[inout] response ConstitutiveResponse instance @@ -154,6 +172,19 @@ class MarmotMaterialFiniteStrain { AlgorithmicModuli< 3 >& algorithmicModuli, const Deformation< 3 >& deformation, const TimeIncrement& timeIncrement ) const; + + /** + * @brief Explicit version of computePlaneStrain for use in explicit time integration schemes. + * @note The default implementation calls computePlaneStrain and ignores the algorithmic tangent. + */ + virtual void computePlaneStrainExplicit( ConstitutiveResponse< 3 >& response, + const Deformation< 3 >& deformation, + const TimeIncrement& timeIncrement ) const + { + AlgorithmicModuli< 3 > algorithmicModuli; + computePlaneStrain( response, algorithmicModuli, deformation, timeIncrement ); + } + /** * @brief Compute stress under plane strain conditions with eigen deformation. * @param[inout] response ConstitutiveResponse instance @@ -170,6 +201,20 @@ class MarmotMaterialFiniteStrain { const Deformation< 3 >& deformation, const TimeIncrement& timeIncrement, const std::tuple< double, double, double >& eigenDeformation ) const; + + /** + * @brief Explicit version of computePlaneStrain with eigen deformation for use in explicit time integration schemes. + * @note The default implementation calls computePlaneStrain and ignores the algorithmic tangent. + */ + virtual void computePlaneStrainExplicit( ConstitutiveResponse< 3 >& response, + const Deformation< 3 >& deformation, + const TimeIncrement& timeIncrement, + const std::tuple< double, double, double >& eigenDeformation ) const + { + AlgorithmicModuli< 3 > algorithmicModuli; + computePlaneStrain( response, algorithmicModuli, deformation, timeIncrement, eigenDeformation ); + } + /** * @brief Compute stress under plane stress conditions. * @param[inout] response ConstitutiveResponse instance @@ -185,6 +230,18 @@ class MarmotMaterialFiniteStrain { const Deformation< 2 >& deformation, const TimeIncrement& timeIncrement ) const; + /** + * @brief Explicit version of computePlaneStress for use in explicit time integration schemes. + * @note The default implementation calls computePlaneStress and ignores the algorithmic tangent. + */ + virtual void computePlaneStressExplicit( ConstitutiveResponse< 2 >& response, + const Deformation< 2 >& deformation, + const TimeIncrement& timeIncrement ) const + { + AlgorithmicModuli< 2 > algorithmicModuli; + computePlaneStress( response, algorithmicModuli, deformation, timeIncrement ); + } + /** * @brief Find the eigen deformation that corresponds to a given eigen stress. * @param initialGuess Initial guess for the eigen deformation. @@ -236,4 +293,23 @@ class MarmotMaterialFiniteStrain { stateVars[i] = 0.0; } } + + /** + * @brief Get the mass density of the material. + * @return Mass density + * + * @note The default implementation throws a std::runtime_error. + */ + virtual double getDensity() const { throw std::runtime_error( "getDensity() not implemented for this material" ); } + + /** + * @brief Get the damping coefficient of the material. + * @return Damping coefficient + * + * @note The default implementation throws a std::runtime_error. + */ + virtual double getDampingCoefficient() const + { + throw std::runtime_error( "getDampingCoefficient() not implemented for this material" ); + } }; diff --git a/modules/core/MarmotMechanicsCore/include/Marmot/MarmotMaterialHypoElastic.h b/modules/core/MarmotMechanicsCore/include/Marmot/MarmotMaterialHypoElastic.h index 0e9d517a..152f783f 100644 --- a/modules/core/MarmotMechanicsCore/include/Marmot/MarmotMaterialHypoElastic.h +++ b/modules/core/MarmotMechanicsCore/include/Marmot/MarmotMaterialHypoElastic.h @@ -128,6 +128,22 @@ class MarmotMaterialHypoElastic { const double* dStrain, const timeInfo& timeInfo ) const = 0; + /** + * Explicit version of @ref computeStress for use in explicit time integration schemes. + * The algorithmic tangent is not needed in explicit schemes and will therefore not be computed. + * @param[in,out] state A state3D instance carrying stress, strain energy, and state variables + * @param[in] dStrain linearized strain increment + * @param[in] timeInfo Structure carrying time information + * + * @note The default implementation calls @ref computeStress and ignores the algorithmic tangent. + * @note Derived classes may override this method for efficiency reasons. + */ + virtual void computeStressExplicit( state3D& state, const double* dStrain, const timeInfo& timeInfo ) const + { + double dStress_dStrain[36]; + computeStress( state, dStress_dStrain, dStrain, timeInfo ); + } + /** * Plane stress implementation of @ref computeStress. */ @@ -182,5 +198,10 @@ class MarmotMaterialHypoElastic { } } - virtual double getDensity() { return -1; } + virtual double getDensity() { throw std::runtime_error( "getDensity() not implemented for this material." ); } + + virtual double getDampingCoefficient() + { + throw std::runtime_error( "getDampingCoefficient() not implemented for this material." ); + } }; diff --git a/modules/elements/DisplacementFiniteElement/include/Marmot/DisplacementFiniteElement.h b/modules/elements/DisplacementFiniteElement/include/Marmot/DisplacementFiniteElement.h index 46456be1..6374a339 100644 --- a/modules/elements/DisplacementFiniteElement/include/Marmot/DisplacementFiniteElement.h +++ b/modules/elements/DisplacementFiniteElement/include/Marmot/DisplacementFiniteElement.h @@ -282,6 +282,27 @@ namespace Marmot::Elements { double dT, double& pNewdT ); + /** # + * @brief Compute internal force only (no tangent stiffness). + * @details Uses the small-strain relation \f$\Delta\boldsymbol{\varepsilon}=\mathbf{B}\,\Delta\mathbf{u}\f$ and + * integrates + * \f[ + * \mathbf{P}_e = \sum_{qp} \mathbf{B}^\mathsf{T} \boldsymbol{\sigma}\, J_0 w. + * \f] + * If pNewdT<1, the routine returns early to signal time step reduction. + * @param QTotal Total displacement vector. + * @param dQ Incremental displacement. + * @param Pe Internal force vector (accumulated). + * @param time Time data forwarded to materials. + * @param dT Time increment. + * @param pNewdT Suggested scaling of dT by the material; if reduced (<1), the routine returns early. + */ + void computeYourselfExplicit( const double* QTotal, + const double* dQ, + double* Pe, + const double* time, + double dT, + double& pNewdT ); /** * @brief Compute consistent mass matrix using material density. * @details \f$\mathbf{M}_e = \sum_{qp} \rho\, \mathbf{N}^\mathsf{T}\mathbf{N}\, J_0 w\f$. @@ -289,11 +310,19 @@ namespace Marmot::Elements { void computeConsistentInertia( double* M ); /** - * @brief Compute lumped mass vector via row-sum of consistent mass. - * @details \f$\mathbf{m}_e = \mathrm{rowsum}(\mathbf{M}_e)\f$. + * @brief Compute the lumped (diagonal) mass matrix. + * @details Uses the manifold-based approach according to + * Yang et al. (2017) "A rigorous and unified mass lumping scheme for higher-order elements", CMAME */ void computeLumpedInertia( double* M ); + /** + * @brief Compute lumped damping matrix using a damping coefficient. + * @details Uses the manifold-based approach according to + * Yang et al. (2017) "A rigorous and unified mass lumping scheme for higher-order elements", CMAME + */ + void computeLumpedDamping( double* C ); + /** * @brief Access a named state view at a quadrature point. * @note Using "sdv" returns the raw material state vector and is deprecated. @@ -589,6 +618,144 @@ namespace Marmot::Elements { } } + template < int nDim, int nNodes > + void DisplacementFiniteElement< nDim, nNodes >::computeYourselfExplicit( const double* QTotal_, + const double* dQ_, + double* Pe_, + const double* time, + double dT, + double& pNewDT ) + { + using namespace Marmot; + using namespace ContinuumMechanics::VoigtNotation; + + Map< const RhsSized > QTotal( QTotal_ ); + Map< const RhsSized > dQ( dQ_ ); + Map< RhsSized > Pe( Pe_ ); + + Voigt S, dE; + + for ( QuadraturePoint& qp : qps ) { + + const BSized& B = qp.B; + dE = B * dQ; + + if constexpr ( nDim == 1 ) { + +/* MarmotMaterialHypoElastic::state1D state; */ +/* MarmotMaterialHypoElastic::timeInfo timeInfo; */ + +/* // set state info */ +/* state.stress = reduce3DVoigt< ParentGeometryElement::voigtSize >( qp.managedStateVars->stress )( 0 ); */ +/* state.strainEnergyDensity = 0.0; */ +/* state.stateVars = qp.managedStateVars->materialStateVars.data(); */ + +/* // set time info */ +/* timeInfo.time = time[1]; */ +/* timeInfo.dT = dT; */ +/* try { */ +/* qp.material->computeUniaxialStress( state, C.data(), dE.data(), timeInfo ); */ +/* } */ +/* catch ( const std::runtime_error& e ) { */ +/* pNewDT = 0.5; */ +/* return; */ +/* } */ +/* Eigen::VectorXd stress1D( 1 ); */ +/* stress1D( 0 ) = state.stress; */ +/* qp.managedStateVars->stress = make3DVoigt< ParentGeometryElement::voigtSize >( stress1D ); */ # + + throw std::runtime_error( "Explicit uniaxial stress not implemented yet" ); + } + + else if constexpr ( nDim == 2 ) { + + if ( sectionType == SectionType::PlaneStress ) { + + /* MarmotMaterialHypoElastic::state2D state; */ + /* MarmotMaterialHypoElastic::timeInfo timeInfo; */ + + /* // set state info */ + /* state.stress = reduce3DVoigt< ParentGeometryElement::voigtSize >( qp.managedStateVars->stress + * ); */ + /* state.strainEnergyDensity = 0.0; */ + /* state.stateVars = qp.managedStateVars->materialStateVars.data(); */ + + /* // set time info */ + /* timeInfo.time = time[1]; */ + /* timeInfo.dT = dT; */ + /* try { */ + /* qp.material->computePlaneStress( state, C.data(), dE.data(), timeInfo ); */ + /* } */ + /* catch ( const std::runtime_error& e ) { */ + /* pNewDT = 0.5; */ + /* return; */ + /* } */ + /* qp.managedStateVars->stress = make3DVoigt< ParentGeometryElement::voigtSize >( state.stress ); */ + /* S = state.stress; */ + throw std::runtime_error( "Explicit plane stress not implemented yet" ); + } + + else if ( sectionType == SectionType::PlaneStrain ) { + + Vector6d dE6 = planeVoigtToVoigt( dE ); + Matrix6d C66; + + // Vector6d S6 = qp.managedStateVars->stress; + MarmotMaterialHypoElastic::state3D state; + MarmotMaterialHypoElastic::timeInfo timeInfo; + + // set state info + state.stress = qp.managedStateVars->stress; + state.strainEnergyDensity = 0.0; + state.stateVars = qp.managedStateVars->materialStateVars.data(); + + // set time info + timeInfo.time = time[1]; + timeInfo.dT = dT; + try { + qp.material->computeStressExplicit( state, dE6.data(), timeInfo ); + } + catch ( const std::runtime_error& e ) { + pNewDT = 0.5; + return; + } + qp.managedStateVars->stress = state.stress; + + S = reduce3DVoigt< ParentGeometryElement::voigtSize >( state.stress ); + } + } + + else if constexpr ( nDim == 3 ) { + if ( sectionType == SectionType::Solid ) { + + MarmotMaterialHypoElastic::state3D state; + MarmotMaterialHypoElastic::timeInfo timeInfo; + + // set state info + state.stress = qp.managedStateVars->stress; + state.strainEnergyDensity = 0.0; + state.stateVars = qp.managedStateVars->materialStateVars.data(); + + // set time info + timeInfo.time = time[1]; + timeInfo.dT = dT; + try { + qp.material->computeStressExplicit( state, dE.data(), timeInfo ); + } + catch ( const std::runtime_error& e ) { + pNewDT = 0.5; + return; + } + qp.managedStateVars->stress = state.stress; + S = state.stress; + } + } + + qp.managedStateVars->strain += make3DVoigt< ParentGeometryElement::voigtSize >( dE ); + + Pe -= B.transpose() * S * qp.J0xW; + } + } template < int nDim, int nNodes > void DisplacementFiniteElement< nDim, nNodes >::setInitialConditions( StateTypes state, const double* values ) { @@ -701,14 +868,49 @@ namespace Marmot::Elements { template < int nDim, int nNodes > void DisplacementFiniteElement< nDim, nNodes >::computeLumpedInertia( double* M ) { - Map< RhsSized > Me( M ); - Me.setZero(); + Map< RhsSized > LMM( M ); + LMM.setZero(); - KeSizedMatrix CMM; - CMM.setZero(); - computeConsistentInertia( CMM.data() ); + constexpr int nNodesLinear = pow( 2, nDim ); + auto linGeometryEl = MarmotGeometryElement< nDim, nNodesLinear >(); + for ( const auto& qp : qps ) { + const auto N_ = this->N( qp.xi ); + const auto N_lin = linGeometryEl.N( qp.xi ); + + VectorXd N_weighted = 0.5 * ( N_ ); + N_weighted.head( nNodesLinear ) += 0.5 * N_lin; - Me = CMM.rowwise().sum(); + const double rho = qp.material->getDensity(); + VectorXd m_ = N_weighted * qp.detJ * qp.weight * rho; + for ( int i = 0; i < nNodes; i++ ) { + for ( int d = 0; d < nDim; d++ ) + LMM( i * nDim + d ) += m_( i ); + } + } + } + + template < int nDim, int nNodes > + void DisplacementFiniteElement< nDim, nNodes >::computeLumpedDamping( double* C ) + { + Map< RhsSized > LCM( C ); + LCM.setZero(); + + constexpr int nNodesLinear = pow( 2, nDim ); + auto linGeometryEl = MarmotGeometryElement< nDim, nNodesLinear >(); + for ( const auto& qp : qps ) { + const auto N_ = this->N( qp.xi ); + const auto N_lin = linGeometryEl.N( qp.xi ); + + VectorXd N_weighted = 0.5 * ( N_ ); + N_weighted.head( nNodesLinear ) += 0.5 * N_lin; + + const double eta = qp.material->getDampingCoefficient(); + VectorXd m_ = N_weighted * qp.detJ * qp.weight * eta; + for ( int i = 0; i < nNodes; i++ ) { + for ( int d = 0; d < nDim; d++ ) + LCM( i * nDim + d ) += m_( i ); + } + } } template < int nDim, int nNodes > diff --git a/modules/elements/DisplacementFiniteStrainULElement/include/Marmot/DisplacementFiniteStrainULElement.h b/modules/elements/DisplacementFiniteStrainULElement/include/Marmot/DisplacementFiniteStrainULElement.h index 40be8574..29338f05 100644 --- a/modules/elements/DisplacementFiniteStrainULElement/include/Marmot/DisplacementFiniteStrainULElement.h +++ b/modules/elements/DisplacementFiniteStrainULElement/include/Marmot/DisplacementFiniteStrainULElement.h @@ -345,6 +345,46 @@ namespace Marmot::Elements { double dT, double& pNewdT ); + /** @brief Compute the negative element residual vector (internal force only, no tangent stiffness) + * + * For a given displacement \f$\mathbf{q}^{(n+1)}\f$ at the current time step \f$t^{(n+1)} = t^{(n)} + \Delta\,t\f$, + * compute the internal work contribution for the negative element residual vector (right hand side of global + * newton) \f$-\int_{V_0}\,\mathbf{N}_{A,i}\,\tau_{ij}\,dV_0\f$. + * + * @param QTotal[in] Pointer to the total element displacement vector at the current time step + * @param dQ[in] Pointer to the increment of the element displacement vector at the current time step + * @param Pe[in,out] Pointer to the negative element residual vector (right hand side of global newton) + * @param time[in] Pointer to the time at the beginning of the current time step + * @param dT[in] Length of the current time step + * @param pNewdT[in,out] Suggested length of the next time step + */ + void computeYourselfExplicit( const double* QTotal, + const double* dQ, + double* Pe, + const double* time, + double dT, + double& pNewdT ); + + /** + * @brief Compute consistent mass matrix using material density. + * @details \f$\mathbf{M}_e = \sum_{qp} \rho\, \mathbf{N}^\mathsf{T}\mathbf{N}\, J_0 w\f$. + */ + void computeConsistentInertia( double* M ); + + /** + * @brief Compute lumped (diagonal) mass matrix using material density. + * @details Using the manifold based approach according to + * Yang et al. (2017) "A rigorous and unified mass lumping scheme for higher-order elements", CMAME + */ + void computeLumpedInertia( double* M ); + + /** + * @brief Compute lumped damping matrix using a damping coefficient. + * @details Using the manifold based approach according to + * Yang et al. (2017) "A rigorous and unified mass lumping scheme for higher-order elements", CMAME + */ + void computeLumpedDamping( double* C ); + /** @brief Get a view to a state variable at a specific quadrature point of the element * @param stateName[in] Name of the state variable * @param qpNumber[in] Number of the quadrature point where the state variable is stored @@ -631,6 +671,107 @@ namespace Marmot::Elements { K.template block< bsU, bsU >( idxU, idxU ) += Map< Matrix< double, bsU, bsU > >( torowmajor( k_UU ).data() ); } + template < int nDim, int nNodes > + void DisplacementFiniteStrainULElement< nDim, nNodes >::computeYourselfExplicit( const double* qTotal, + const double* dQ, + double* rightHandSide, + const double* time, + double dT, + double& pNewDT ) + { + using namespace Fastor; + + const static Tensor< double, nDim, nDim > I( + ( Eigen::Matrix< double, nDim, nDim >() << Eigen::Matrix< double, nDim, nDim >::Identity() ).finished().data() ); + + // in ... + const auto qU_np = TensorMap< const double, nNodes, nDim >( qTotal ); + + // ... and out: residuals + TensorMap< double, nNodes, nDim > r_U( rightHandSide ); + + Eigen::Map< Eigen::VectorXd > rhs( rightHandSide, sizeLoadVector ); + + for ( auto& qp : qps ) { + + using namespace Marmot::FastorIndices; + + const auto& dNdX_ = qp.dNdX; + + const auto dNdX = Tensor< double, nDim, nNodes >( dNdX_.data(), ColumnMajor ); + + const auto F_np = evaluate( einsum< Ai, jA >( qU_np, dNdX ) + I ); + + const Material::Deformation< nDim > deformation = { F_np }; + + const Material::TimeIncrement timeIncrement{ time[1], dT }; + + Material::ConstitutiveResponse< nDim > response = { 0, 0, 0, nullptr }; + try { + if constexpr ( nDim == 2 ) { + + if ( sectionType == SectionType::PlaneStrain ) { + + using namespace Marmot; + + Material::ConstitutiveResponse< 3 > + response3D{ FastorStandardTensors::Tensor33d( qp.managedStateVars->stress.data(), Fastor::ColumnMajor ), + -1.0, + -1.0, + qp.managedStateVars->materialStateVars.data() }; + + Material::Deformation< 3 > deformation3D{ + expandTo3D( deformation.F ), + }; + + deformation3D.F( 2, 2 ) = 1.0; + + if ( hasEigenDeformation ) + qp.material->computePlaneStrainExplicit( response3D, + deformation3D, + timeIncrement, + { qp.managedStateVars->F0_XX, + qp.managedStateVars->F0_YY, + qp.managedStateVars->F0_ZZ } ); + else + qp.material->computePlaneStrainExplicit( response3D, deformation3D, timeIncrement ); + + response = { reduceTo2D< U, U >( response3D.tau ), + response3D.rho, + response3D.elasticEnergyDensity, + qp.managedStateVars->materialStateVars.data() }; + + qp.managedStateVars->stress = Marmot::mapEigenToFastor( response3D.tau ).reshaped(); + } + } + else { + response = { Marmot::FastorStandardTensors::Tensor33d( qp.managedStateVars->stress.data(), ColumnMajor ), + -1.0, + -1.0, + qp.managedStateVars->materialStateVars.data() }; + + qp.material->computeStressExplicit( response, deformation, timeIncrement ); + + // implicit conversion to col major + qp.managedStateVars->stress = Marmot::mapEigenToFastor( response.tau ).reshaped(); + } + } + catch ( const std::runtime_error& ) { + pNewDT = 0.25; + return; + } + const auto dNdx = evaluate( einsum< ji, jA >( inv( F_np ), dNdX ) ); + + const double& J0xW = qp.J0xW; + + const auto& tau = response.tau; + + // r[ node, dim ] (swap to abuse directly colmajor layout) + // directly operate via TensorMap + r_U -= ( +einsum< iA, ij >( dNdx, tau ) ) * J0xW; + } + } + template < int nDim, int nNodes > void DisplacementFiniteStrainULElement< nDim, nNodes >::computeDistributedLoad( MarmotElement::DistributedLoadTypes loadType, @@ -758,6 +899,67 @@ namespace Marmot::Elements { r_U += this->NB( this->N( qp.xi ) ).transpose() * f * qp.J0xW; } + template < int nDim, int nNodes > + void DisplacementFiniteStrainULElement< nDim, nNodes >::computeConsistentInertia( double* M ) + { + Eigen::Map< KSizedMatrix > Me( M ); + Me.setZero(); + + for ( const auto& qp : qps ) { + const auto N_ = this->NB( this->N( qp.xi ) ); + const double rho = qp.material->getDensity(); + Me += N_.transpose() * N_ * qp.J0xW * rho; + } + } + + template < int nDim, int nNodes > + void DisplacementFiniteStrainULElement< nDim, nNodes >::computeLumpedInertia( double* M ) + { + Eigen::Map< RhsSized > LMM( M ); + LMM.setZero(); + + constexpr int nNodesLinear = std::pow( 2, nDim ); + auto linGeometryEl = MarmotGeometryElement< nDim, nNodesLinear >(); + for ( const auto& qp : qps ) { + const auto N_ = this->N( qp.xi ); + const auto N_lin = linGeometryEl.N( qp.xi ); + + Eigen::VectorXd N_weighted = 0.5 * ( N_ ); + N_weighted.head( nNodesLinear ) += 0.5 * N_lin; + + const double rho = qp.material->getDensity(); + Eigen::VectorXd m_ = N_weighted * qp.J0xW * rho; + for ( int i = 0; i < nNodes; i++ ) { + for ( int d = 0; d < nDim; d++ ) + LMM( i * nDim + d ) += m_( i ); + } + } + } + + template < int nDim, int nNodes > + void DisplacementFiniteStrainULElement< nDim, nNodes >::computeLumpedDamping( double* C ) + { + Eigen::Map< RhsSized > LCM( C ); + LCM.setZero(); + + constexpr int nNodesLinear = std::pow( 2, nDim ); + auto linGeometryEl = MarmotGeometryElement< nDim, nNodesLinear >(); + for ( const auto& qp : qps ) { + const auto N_ = this->N( qp.xi ); + const auto N_lin = linGeometryEl.N( qp.xi ); + + Eigen::VectorXd N_weighted = 0.5 * ( N_ ); + N_weighted.head( nNodesLinear ) += 0.5 * N_lin; + + const double eta = qp.material->getDampingCoefficient(); + Eigen::VectorXd m_ = N_weighted * qp.J0xW * eta; + for ( int i = 0; i < nNodes; i++ ) { + for ( int d = 0; d < nDim; d++ ) + LCM( i * nDim + d ) += m_( i ); + } + } + } + template < int nDim, int nNodes > std::vector< double > DisplacementFiniteStrainULElement< nDim, nNodes >::getCoordinatesAtCenter() { @@ -804,6 +1006,13 @@ namespace Marmot::Elements { const double* time, double dT, double& pNewdT ); + + void computeYourselfExplicit( const double* QTotal, + const double* dQ, + double* Pe, + const double* time, + double dT, + double& pNewdT ); }; template < int nNodes > @@ -958,4 +1167,100 @@ namespace Marmot::Elements { torowmajor( k_UU ).data() ); } + template < int nNodes > + void AxiSymmetricDisplacementFiniteStrainULElement< nNodes >::computeYourselfExplicit( const double* qTotal, + const double* dQ, + double* rightHandSide, + const double* time, + double dT, + double& pNewDT ) + { + constexpr int nDim = 2; + + using Parent = DisplacementFiniteStrainULElement< 2, nNodes >; + const auto sizeLoadVector = DisplacementFiniteStrainULElement< 2, nNodes >::sizeLoadVector; + using Material = MarmotMaterialFiniteStrain; + + using namespace Fastor; + + const static Tensor< double, nDim, nDim > I( + ( Eigen::Matrix< double, nDim, nDim >() << Eigen::Matrix< double, nDim, nDim >::Identity() ).finished().data() ); + + // in ... + const auto qU_np = TensorMap< const double, nNodes, nDim >( qTotal ); + + // ... and out: residuals + TensorMap< double, nNodes, nDim > r_U( rightHandSide ); + + Eigen::Map< Eigen::VectorXd > rhs( rightHandSide, sizeLoadVector ); + + for ( auto& qp : Parent::qps ) { + + using namespace Marmot::FastorIndices; + + auto N_ = this->N( qp.xi ); + const auto& dNdX_ = qp.dNdX; + + Eigen::Vector2d coords = this->NB( N_ ) * this->coordinates; + const double r = coords[0]; + + const auto N = Tensor< double, nNodes >( N_.data() ); + const auto dNdX = Tensor< double, nDim, nNodes >( dNdX_.data(), ColumnMajor ); + + const auto u_np = evaluate( einsum< A, Ai >( N, qU_np ) ); + + const auto F_np = evaluate( einsum< Ai, jA >( qU_np, dNdX ) + I ); + + const Material::Deformation< nDim > deformation = { + F_np, + }; + + const Material ::TimeIncrement timeIncrement{ time[0], dT }; + + Material::ConstitutiveResponse< nDim > response; + + using namespace Marmot; + Material::ConstitutiveResponse< 3 > + response3D{ FastorStandardTensors::Tensor33d( qp.managedStateVars->stress.data(), Fastor::ColumnMajor ), + -1.0, + -1.0, + qp.managedStateVars->materialStateVars.data() }; + + Material::Deformation< 3 > deformation3D{ expandTo3D( deformation.F ) }; + + deformation3D.F( 2, 2 ) = 1 + u_np[0] / r; + + try { + qp.material->computePlaneStrainExplicit( response3D, deformation3D, timeIncrement ); + } + catch ( const std::runtime_error& ) { + pNewDT = 0.25; + return; + } + response = { reduceTo2D< U, U >( response3D.tau ), + response3D.rho, + response3D.elasticEnergyDensity, + qp.managedStateVars->materialStateVars.data() }; + + qp.managedStateVars->stress = Marmot::mapEigenToFastor( response3D.tau ).reshaped(); + + const auto dNdx = evaluate( einsum< ji, jA >( inv( F_np ), dNdX ) ); + + const double J0xWxRx2Pi = qp.J0xW * 2 * Constants::Pi * r; + + const auto& tau = response.tau; + + // r[ node, dim ] (swap to abuse directly colmajor layout) + // directly operate via TensorMap + r_U -= ( +einsum< iA, ij >( dNdx, tau ) ) * J0xWxRx2Pi; + + const double F33 = 1 + u_np[0] / r; + const double invF33 = 1. / F33; + + for ( int A = 0; A < nNodes; A++ ) { + r_U( A, 0 ) -= ( +N( A ) * invF33 * response3D.tau( 2, 2 ) / r ) * J0xWxRx2Pi; + } + } + } + } // namespace Marmot::Elements diff --git a/modules/materials/ADLinearElastic/include/Marmot/ADLinearElastic.h b/modules/materials/ADLinearElastic/include/Marmot/ADLinearElastic.h index 13be6b3b..b2c79273 100644 --- a/modules/materials/ADLinearElastic/include/Marmot/ADLinearElastic.h +++ b/modules/materials/ADLinearElastic/include/Marmot/ADLinearElastic.h @@ -50,6 +50,9 @@ namespace Marmot::Materials { ADLinearElastic( const double* materialProperties, int nMaterialProperties, int materialNumber ); + double getDensity() override; + double getDampingCoefficient() override; + protected: void computeStressAD( state3DAD& state, const autodiff::dual* dStrain, const timeInfo& timeInfo ) const; diff --git a/modules/materials/ADLinearElastic/src/ADLinearElastic.cpp b/modules/materials/ADLinearElastic/src/ADLinearElastic.cpp index cdc70338..17ce75c1 100644 --- a/modules/materials/ADLinearElastic/src/ADLinearElastic.cpp +++ b/modules/materials/ADLinearElastic/src/ADLinearElastic.cpp @@ -17,7 +17,23 @@ namespace Marmot::Materials { E( materialProperties[0] ), nu( materialProperties[1] ) { - assert( nMaterialProperties == 2 ); + assert( nMaterialProperties >= 2 ); + } + + double ADLinearElastic::getDensity() + { + if ( nMaterialProperties >= 3 ) + return materialProperties[2]; + throw std::runtime_error( + std::string( MakeString() << __PRETTY_FUNCTION__ << ": Density not specified for ADLinearElastic." ) ); + } + + double ADLinearElastic::getDampingCoefficient() + { + if ( nMaterialProperties >= 4 ) + return materialProperties[3]; + throw std::runtime_error( std::string( + MakeString() << __PRETTY_FUNCTION__ << ": Damping coefficient not specified for ADLinearElastic." ) ); } void ADLinearElastic::computeStressAD( state3DAD& state, const autodiff::dual* dStrain, diff --git a/modules/materials/B4/include/Marmot/B4.h b/modules/materials/B4/include/Marmot/B4.h index ac7f035e..d2564126 100644 --- a/modules/materials/B4/include/Marmot/B4.h +++ b/modules/materials/B4/include/Marmot/B4.h @@ -159,6 +159,9 @@ namespace Marmot::Materials { B4( const double* materialProperties, int nMaterialProperties, int materialLabel ); + double getDensity() override; + double getDampingCoefficient() override; + void computeStress( state3D& state, double* dStressDDStrain, const double* dStrain, diff --git a/modules/materials/B4/src/B4.cpp b/modules/materials/B4/src/B4.cpp index 31190b71..325fbf76 100644 --- a/modules/materials/B4/src/B4.cpp +++ b/modules/materials/B4/src/B4.cpp @@ -1,6 +1,7 @@ #include "Marmot/B4.h" #include "Marmot/B4Shrinkage.h" #include "Marmot/MarmotElasticity.h" +#include "Marmot/MarmotJournal.h" #include "Marmot/MarmotMaterialHypoElastic.h" #include "Marmot/MarmotTypedefs.h" #include @@ -169,4 +170,19 @@ namespace Marmot::Materials { return; } + double B4::getDensity() + { + if ( nMaterialProperties >= 23 ) + return materialProperties[22]; + throw std::runtime_error( std::string( MakeString() << __PRETTY_FUNCTION__ << ": Density not specified for B4." ) ); + } + + double B4::getDampingCoefficient() + { + if ( nMaterialProperties >= 24 ) + return materialProperties[23]; + throw std::runtime_error( + std::string( MakeString() << __PRETTY_FUNCTION__ << ": Damping coefficient not specified for B4." ) ); + } + } // namespace Marmot::Materials diff --git a/modules/materials/CompressibleNeoHooke/include/Marmot/CompressibleNeoHooke.h b/modules/materials/CompressibleNeoHooke/include/Marmot/CompressibleNeoHooke.h index f3082d7b..45f0aff1 100644 --- a/modules/materials/CompressibleNeoHooke/include/Marmot/CompressibleNeoHooke.h +++ b/modules/materials/CompressibleNeoHooke/include/Marmot/CompressibleNeoHooke.h @@ -28,7 +28,6 @@ #pragma once #include "Marmot/MarmotJournal.h" #include "Marmot/MarmotMaterialFiniteStrain.h" -#include namespace Marmot::Materials { @@ -80,6 +79,32 @@ namespace Marmot::Materials { const Deformation< 3 >&, const TimeIncrement& ) const override; + /** + * @brief Get material density. + * @return Density value. + */ + double getDensity() const override + { + if ( this->nMaterialProperties < 3 ) { + throw std::runtime_error( + std::string( MakeString() << __PRETTY_FUNCTION__ << ": No density given! nMaterialProperties < 3." ) ); + } + return this->materialProperties[2]; + } + + /** + * @brief Get damping coefficient. + * @return Damping coefficient. + */ + double getDampingCoefficient() const override + { + if ( this->nMaterialProperties < 4 ) { + throw std::runtime_error( std::string( + MakeString() << __PRETTY_FUNCTION__ << ": No damping coefficient given! nMaterialProperties < 4." ) ); + } + return this->materialProperties[3]; + } + /** * @brief Initialize the state layout (no state variables here). */ diff --git a/modules/materials/FiniteStrainJ2Plasticity/include/Marmot/FiniteStrainJ2Plasticity.h b/modules/materials/FiniteStrainJ2Plasticity/include/Marmot/FiniteStrainJ2Plasticity.h index 5c1f0779..aae1d3ef 100644 --- a/modules/materials/FiniteStrainJ2Plasticity/include/Marmot/FiniteStrainJ2Plasticity.h +++ b/modules/materials/FiniteStrainJ2Plasticity/include/Marmot/FiniteStrainJ2Plasticity.h @@ -30,10 +30,12 @@ #include "Marmot/MarmotEnergyDensityFunctions.h" #include "Marmot/MarmotFastorTensorBasics.h" #include "Marmot/MarmotFiniteStrainPlasticity.h" +#include "Marmot/MarmotJournal.h" #include "Marmot/MarmotMaterialFiniteStrain.h" #include "Marmot/MarmotMath.h" #include "Marmot/MarmotTypedefs.h" #include +#include #include #include @@ -178,6 +180,32 @@ namespace Marmot::Materials { const Deformation< 3 >& deformation, const TimeIncrement& timeIncrement ) const; + /** + * @brief Get material density. + * @return Density value. + */ + double getDensity() const override + { + if ( this->nMaterialProperties < 8 ) { + throw std::runtime_error( + std::string( MakeString() << __PRETTY_FUNCTION__ << ": No density given! nMaterialProperties < 8." ) ); + } + return this->density; + } + + /** + * @brief Get damping coefficient. + * @return Damping coefficient. + */ + double getDampingCoefficient() const override + { + if ( this->nMaterialProperties < 9 ) { + throw std::runtime_error( std::string( + MakeString() << __PRETTY_FUNCTION__ << ": No damping coefficient given! nMaterialProperties < 9." ) ); + } + return this->materialProperties[8]; + } + void initializeStateLayout() override { stateLayout.add( "Fp", 9 ); // plastic deformation gradient diff --git a/modules/materials/LinearElastic/include/Marmot/LinearElastic.h b/modules/materials/LinearElastic/include/Marmot/LinearElastic.h index 02da6f52..7d4afb9f 100644 --- a/modules/materials/LinearElastic/include/Marmot/LinearElastic.h +++ b/modules/materials/LinearElastic/include/Marmot/LinearElastic.h @@ -44,6 +44,7 @@ namespace Marmot::Materials { LinearElastic( const double* materialProperties, int nMaterialProperties, int materialNumber ); double getDensity() override; + double getDampingCoefficient() override; protected: /// @brief Type of isotropic and anisotropic behavior. diff --git a/modules/materials/LinearElastic/src/LinearElastic.cpp b/modules/materials/LinearElastic/src/LinearElastic.cpp index f531e6b0..a473d66d 100644 --- a/modules/materials/LinearElastic/src/LinearElastic.cpp +++ b/modules/materials/LinearElastic/src/LinearElastic.cpp @@ -1,6 +1,5 @@ #include "Marmot/LinearElastic.h" #include "Marmot/MarmotElasticity.h" -#include "Marmot/MarmotJournal.h" #include "Marmot/MarmotMath.h" #include "Marmot/MarmotTypedefs.h" #include "Marmot/MarmotUtility.h" @@ -15,7 +14,7 @@ namespace Marmot::Materials { LinearElastic::LinearElastic( const double* materialProperties, int nMaterialProperties, int materialNumber ) : MarmotMaterialHypoElastic::MarmotMaterialHypoElastic( materialProperties, nMaterialProperties, materialNumber ), // clang-format off - anisotropicType( nMaterialProperties == 2 || nMaterialProperties == 11 || nMaterialProperties == 15 ? static_cast< Type >( nMaterialProperties ) : static_cast< Type >( nMaterialProperties - 1 ) ), + anisotropicType( nMaterialProperties >= 15 ? Type::Orthotropic : (nMaterialProperties >= 11 ? Type::TransverseIsotropic : Type::Isotropic) ), E1( materialProperties[0] ), E2( anisotropicType == Type::Isotropic ? E1 : materialProperties[1] ), E3( anisotropicType == Type::Orthotropic ? materialProperties[2] : E2 ), @@ -90,9 +89,23 @@ namespace Marmot::Materials { double LinearElastic::getDensity() { - if ( nMaterialProperties == 3 || nMaterialProperties == 12 || nMaterialProperties == 16 ) - return materialProperties[nMaterialProperties - 1]; - else - throw std::runtime_error( "Density not specified for this material" ); + int base_props = static_cast< int >( anisotropicType ); + if ( nMaterialProperties >= base_props + 1 ) + return materialProperties[base_props]; + else { + throw std::runtime_error( + std::string( MakeString() << __PRETTY_FUNCTION__ << ": Density not specified for this material." ) ); + } + } + + double LinearElastic::getDampingCoefficient() + { + int base_props = static_cast< int >( anisotropicType ); + if ( nMaterialProperties >= base_props + 2 ) + return materialProperties[base_props + 1]; + else { + throw std::runtime_error( std::string( + MakeString() << __PRETTY_FUNCTION__ << ": Damping coefficient not specified for this material." ) ); + } } } // namespace Marmot::Materials diff --git a/modules/materials/VonMises/include/Marmot/VonMises.h b/modules/materials/VonMises/include/Marmot/VonMises.h index 5ae26e33..4b70fe0f 100644 --- a/modules/materials/VonMises/include/Marmot/VonMises.h +++ b/modules/materials/VonMises/include/Marmot/VonMises.h @@ -47,6 +47,7 @@ namespace Marmot::Materials { const double* dStrain, const timeInfo& timeInfo ) const override; + void computeStressExplicit( state3D& state, const double* dStrain, const timeInfo& timeInfo ) const override; /** * @brief Get material density. * @return Density value. @@ -54,6 +55,8 @@ namespace Marmot::Materials { */ double getDensity() override; + double getDampingCoefficient() override; + void initializeStateLayout() override { stateLayout.add( "kappa", 1 ); diff --git a/modules/materials/VonMises/src/VonMises.cpp b/modules/materials/VonMises/src/VonMises.cpp index b6d4256b..2908eff7 100644 --- a/modules/materials/VonMises/src/VonMises.cpp +++ b/modules/materials/VonMises/src/VonMises.cpp @@ -2,7 +2,6 @@ #include "Marmot/MarmotConstants.h" #include "Marmot/MarmotElasticity.h" #include "Marmot/MarmotTypedefs.h" -#include "Marmot/MarmotVoigt.h" #include "Marmot/VonMisesConstants.h" namespace Marmot::Materials { @@ -12,11 +11,22 @@ namespace Marmot::Materials { double VonMisesModel::getDensity() { - if ( this->nMaterialProperties < 7 ) - throw std::runtime_error( MakeString() << __PRETTY_FUNCTION__ << ": No density given! nMaterialProperties < 7" ); + if ( this->nMaterialProperties < 7 ) { + throw std::runtime_error( + std::string( MakeString() << __PRETTY_FUNCTION__ << ": No density given! nMaterialProperties < 7." ) ); + } return this->materialProperties[6]; } + double VonMisesModel::getDampingCoefficient() + { + if ( this->nMaterialProperties < 8 ) { + throw std::runtime_error( std::string( + MakeString() << __PRETTY_FUNCTION__ << ": No damping coefficient given! nMaterialProperties < 8." ) ); + } + return this->materialProperties[7]; + } + void VonMisesModel::computeStress( state3D& state, double* dStress_dStrain, const double* dStrain, @@ -118,4 +128,97 @@ namespace Marmot::Materials { } } + void VonMisesModel::computeStressExplicit( state3D& state, const double* dStrain, const timeInfo& timeInfo ) const + + { + // elasticity parameters + const double& E = this->materialProperties[0]; + const double& nu = this->materialProperties[1]; + // plasticity parameters + const double& yieldStress = this->materialProperties[2]; + const double& HLin = this->materialProperties[3]; + const double& deltaYieldStress = this->materialProperties[4]; + const double& delta = this->materialProperties[5]; + + // map to stress, strain and tangent + mVector6d S( state.stress.data() ); + const auto dE = Map< const Vector6d >( dStrain ); + + // compute elastic stiffness + const auto Cel = ContinuumMechanics::Elasticity::Isotropic::stiffnessTensor( E, nu ); + + // handle zero strain increment + if ( dE.isZero( 1e-14 ) ) { + return; + } + + // get current hardening variable + double& kappa = stateLayout.getAs< double& >( state.stateVars, "kappa" ); + + // isotropic hardening law + auto fy = [&]( double kappa_ ) { + return yieldStress + HLin * kappa_ + deltaYieldStress * ( 1. - std::exp( -delta * kappa_ ) ); + }; + + // derivative of fy wrt dKappa + auto dfy_ddKappa = [&]( double kappa_ ) { return HLin + deltaYieldStress * delta * std::exp( -delta * kappa_ ); }; + + // yield function + auto f = [&]( double rho_, double kappa_ ) { return rho_ - Constants::sqrt2_3 * fy( kappa_ ); }; + + // compute elastic predictor + const Vector6d trialStress = S + Cel * dE; + + using namespace ContinuumMechanics::VoigtNotation; + const double rhoTrial = std::sqrt( 2. * Invariants::J2( trialStress ) ); + + if ( f( rhoTrial, kappa ) >= 0.0 ) { + // plastic step + const double G = E / ( 2. * ( 1. + nu ) ); + + auto g = [&]( double deltaKappa ) { + return rhoTrial - Constants::sqrt6 * G * deltaKappa - Constants::sqrt2_3 * fy( kappa + deltaKappa ); + }; + + // variables for return mapping + int counter = 0; + double dKappa = 0; + double dLambda = 0; + double dg_ddKappa = 0; + + // compute return mapping direction + Vector6d n = ContinuumMechanics::VoigtNotation::IDev * trialStress / rhoTrial; + + double g_val = g( dKappa ); + + // always perform 5 iterations for explicit scheme + // this allows to avoid if-conditions inside the time stepping loop + // and improves performance + while ( counter < 5 ) { + + // compute derivative of g wrt kappa + dg_ddKappa = -Constants::sqrt6 * G - Constants::sqrt2_3 * dfy_ddKappa( kappa + dKappa ); + + // update dKappa and iteration counter + dKappa -= g_val / dg_ddKappa; + g_val = g( dKappa ); + counter += 1; + } + + if ( std::abs( g_val ) > VonMisesConstants::innerNewtonTol ) { + throw std::runtime_error( "return mapping failed to converge in VonMisesModel::computeStressExplicit" ); + } + + dLambda = Constants::sqrt3_2 * dKappa; + + // update material state + S = trialStress - 2. * G * dLambda * n; + kappa = kappa + dKappa; + } + else { + // elastic step + S = trialStress; + } + } + } // namespace Marmot::Materials