diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 6a84eb35..b38e5ece 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -54,7 +54,9 @@ enum InterfaceElementType LINEAR_QUAD, ///! 2D linear quadrilateral LINEAR_TET, ///! 3D linear tetrahedron (volume methods and test mesh class) LINEAR_HEX, ///! 3D linear hexahedron (volume methods and test mesh class) - NUM_CONTACT_ELEMENTS = LINEAR_HEX + PARTICLE_CIRCLE, ///! 2D computational geometry (volume methods and test mesh class) + PARTICLE_SPHERE, ///! 3D computational geometry (volume methods and test mesh class) + NUM_CONTACT_ELEMENTS = PARTICLE_SPHERE }; /*! @@ -86,6 +88,7 @@ enum ContactMode SURFACE_TO_SURFACE_CONFORMING, ///! conforming surface-to-surface interaction SURFACE_TO_VOLUME, ///! surface-to-volume interaction VOLUME_TO_VOLUME, ///! volume-to-volume interaction + PARTICLE_TO_PARTICLE, ///! particle-to-particle interaction NUM_CONTACT_MODES }; @@ -184,6 +187,7 @@ enum KinematicPenaltyCalculation { KINEMATIC_CONSTANT, ///! Constant penalty stiffness applied to all contacting face-pairs KINEMATIC_ELEMENT, ///! Element-wise penalty stiffness calculation + KINEMATIC_PARTICLE, ///! Particle-wise penalty stiffness calculation NUM_KINEMATIC_PENALTY_CALCULATION }; diff --git a/src/tribol/mesh/MeshData.cpp b/src/tribol/mesh/MeshData.cpp index 8d44fd1f..c082dbaf 100644 --- a/src/tribol/mesh/MeshData.cpp +++ b/src/tribol/mesh/MeshData.cpp @@ -219,6 +219,188 @@ bool MeshElemData::isValidRatePenalty( PenaltyEnforcementOptions& pen_options ) } // end MeshElemData::isValidRatePenalty() +/////////////////////////////////////////// +// // +// Routines for the MeshElemData struct // +// // +/////////////////////////////////////////// +MeshParticleData::MeshParticleData() +{ + this->m_mat_mod = nullptr; + this->m_thickness = nullptr; + return; +} + +//------------------------------------------------------------------------------ +MeshParticleData::~MeshParticleData() +{ + // Don't delete any pointers to host code data (e.g. force, velocity, + // penalty data, etc.). Don't nullify pointers as the host code + // could register these pointers at the very beginning of a simulation + // and not every cycle. +} + + +//------------------------------------------------------------------------------ +bool MeshParticleData::isValidKinematicPenalty( PenaltyEnforcementOptions& pen_options ) +{ + // Note, this routine is, and should be called only for non-null meshes + KinematicPenaltyCalculation kin_calc = pen_options.kinematic_calculation; + + // check kinematic penalty calculation data + if ( !in_range( kin_calc, NUM_KINEMATIC_PENALTY_CALCULATION) ) + { + // warning already issued when penalty options were set. If penalty options + // were not set, the following else-if will catch this + return false; + } + // a kinematic penalty should always be set. Right now Tribol does not support + // rate only enforcement + else if ( !pen_options.is_kinematic_calculation_set() ) + { + SLIC_WARNING( "MeshParticleData::isValidKinematic(): kinematic penalty calculation data not set; " << + "call tribol::setPenaltyOptions()." ); + return false; + } + + switch (kin_calc) + { + case KINEMATIC_CONSTANT: + { + if ( !this->m_is_kinematic_constant_penalty_set ) + { + SLIC_WARNING("MeshParticleData::isValidKinematicPenalty(): " << + "single stiffness penalty not set."); + return false; + } + else if ( this->m_penalty_stiffness < pen_options.tiny_penalty ) + { + SLIC_WARNING("MeshParticleData::isValidKinematicPenalty(): " << + "single stiffness penalty less than threshold (" << + pen_options.tiny_penalty << "). Consider increasing " << + "for your problem."); + return false; + } + break; + } // end case KINEMATIC_CONSTANT + case KINEMATIC_PARTICLE: + { + if ( !this->m_is_kinematic_particle_penalty_set ) + { + SLIC_WARNING("MeshParticleData::isValidKinematicPenalty(): " << + "particle-wise penalty data not set."); + return false; + } + + // check for positive material modulus and thickness values + bool isValidMatMod = true; + for (int i=0; im_numParticles; ++i) + { + if (this->m_mat_mod[i] <= 0.) + { + isValidMatMod = false; + } + + SLIC_WARNING_IF(!isValidMatMod, "MeshParticleData::isValidKinematicPenalty(): " << + "invalid nonpositive particle material modulus encountered."); + + if (!isValidMatMod) + { + return false; + } + } // end for loop over particles + break; + } // end case KINEMATIC_PARTICLE + default: + // no-op, quiet compiler + break; + } // end switch over kinematic calculation + + return true; + +} // end MeshParticleData::isValidKinematicPenalty() + +//------------------------------------------------------------------------------ +bool MeshParticleData::isValidRatePenalty( PenaltyEnforcementOptions& pen_options ) +{ + // Note, this method is and should only be called for non-null meshes + RatePenaltyCalculation rate_calc = pen_options.rate_calculation; + + // check rate penalty calculation data + if ( !in_range( rate_calc, NUM_RATE_PENALTY_CALCULATION) ) + { + // warning already issued when penalty options were set. If penalty options + // were not set, the following else-if will catch this + return false; + } + // the rate_calc could be set to NONE and this boolean will be true + else if ( !pen_options.is_rate_calculation_set() ) + { + SLIC_WARNING( "MeshParticleData::isValidRatePenalty(): rate penalty calculation data not set. " << + "call tribol::setPenaltyOptions()." ); + return false; + } + + switch (rate_calc) + { + case NO_RATE_PENALTY: + { + // double check that mesh booleans are consistent with rate penalty + // calculation option + if (this->m_is_rate_constant_penalty_set) + { + this->m_is_rate_constant_penalty_set = false; + } + else if (this->m_is_rate_percent_penalty_set) + { + this->m_is_rate_percent_penalty_set = false; + } + break; + } // end case NONE + case RATE_CONSTANT: + { + if ( !this->m_is_rate_constant_penalty_set ) + { + SLIC_WARNING("MeshParticleData::isValidRatePenalty(): " << + "constant rate penalty data not set."); + return false; + } + else if ( this->m_rate_penalty_stiffness < pen_options.tiny_penalty ) + { + SLIC_WARNING("MeshParticleData::isValidRatePenalty(): " << + "constant rate penalty less than threshold (" << + pen_options.tiny_penalty << "). Consider increasing " << + "for your problem."); + return false; + } + break; + } // end case RATE_CONSTANT + case RATE_PERCENT: + { + if ( !this->m_is_rate_percent_penalty_set ) + { + SLIC_WARNING("MeshParticleData::isValidRatePenalty(): " << + "percent rate penalty data not set."); + return false; + } + else if ( this->m_rate_percent_stiffness < pen_options.tiny_penalty || + this->m_rate_percent_stiffness > (1.-pen_options.tiny_penalty) ) + { + SLIC_WARNING("MeshParticleData::isValidRatePenalty(): " << + "rate percent penalty not in (0,1)."); + return false; + } + break; + } // end case RATE_PERCENT + default: + // no-op, quiet compiler + break; + } // end switch on rate calculation + + return true; + +} // end MeshParticleData::isValidRatePenalty() + //------------------------------------------------------------------------------ ///////////////////////////////// @@ -230,16 +412,19 @@ MeshData::MeshData() : m_elementType( tribol::UNDEFINED_ELEMENT ) , m_lengthNodalData(0) , m_numSurfaceNodes(0) + , m_numParticles(0) , m_numCells(0) , m_numNodesPerCell(0) , m_isValid(true) , m_positionX(nullptr), m_positionY(nullptr), m_positionZ(nullptr) , m_dispX(nullptr), m_dispY(nullptr), m_dispZ(nullptr) , m_forceX(nullptr), m_forceY(nullptr), m_forceZ(nullptr) + , m_momentX(nullptr), m_momentY(nullptr), m_momentZ(nullptr) , m_velX(nullptr), m_velY(nullptr), m_velZ(nullptr) , m_nX(nullptr), m_nY(nullptr), m_nZ(nullptr) , m_cX(nullptr), m_cY(nullptr), m_cZ(nullptr) , m_faceRadius(nullptr) + , m_particleRadius(nullptr) , m_area(nullptr) , m_connectivity(nullptr) , m_sortedSurfaceNodeIds(nullptr) @@ -256,6 +441,7 @@ void MeshData::allocateArrays(int dimension) m_cY = new real[m_numCells]; m_area = new real[m_numCells]; m_faceRadius = new real[m_numCells]; + m_particleRadius = new real[m_numParticles]; m_nZ = nullptr; m_cZ = nullptr; @@ -335,6 +521,12 @@ void MeshData::deallocateArrays() m_faceRadius = nullptr; } + if (m_particleRadius != nullptr) + { + delete[] m_particleRadius; + m_particleRadius = nullptr; + } + if (m_area != nullptr) { delete[] m_area; @@ -902,6 +1094,45 @@ int MeshData::checkPenaltyData( PenaltyEnforcementOptions& p_enfrc_options ) } // end switch over constraint types } // end if-non-null mesh + if (this->m_numParticles>0) + { + PenaltyConstraintType constraint_type = p_enfrc_options.constraint_type; + // switch over penalty enforcement options and check for required data + switch (constraint_type) + { + case KINEMATIC: + { + if (!this->m_particleData.isValidKinematicPenalty( p_enfrc_options )) + { + err = 1; + } + break; + } // end KINEMATIC case + + case KINEMATIC_AND_RATE: + { + if (!this->m_particleData.isValidKinematicPenalty( p_enfrc_options )) + { + err = 1; + } + if (!this->m_particleData.isValidRatePenalty( p_enfrc_options )) + { + err = 1; + } + if (!this->m_nodalFields.m_is_velocity_set) + { + SLIC_WARNING("Nodal velocities not set or null pointers; please set for " << + "use with gap rate penalty enforcement."); + err = 1; + } + break; + } // end case KINEMATIC_AND_RATE + default: + // no-op, quiet compiler + break; + } // end switch over constraint types + } // end if-non-null mesh + return err; } // end MeshData::checkPenaltyData() @@ -916,8 +1147,13 @@ int MeshData::localIdFromConn( const int connectivityId ) //------------------------------------------------------------------------------ void MeshData::print(std::ostream& os) const { - const int num_verts = m_lengthNodalData; + if ( m_lengthNodalData > 0 && m_numParticles > 0 ) + { + SLIC_ERROR("MeshData::print(): Cannot have nodes and particles in same MeshData instance."); + } + const int num_verts = m_lengthNodalData + m_numParticles; const int num_elem = m_numCells; + const int num_particles = m_numParticles; const int dim = m_positionZ == nullptr ? 2 : 3; if(num_verts <= 0) @@ -946,6 +1182,17 @@ void MeshData::print(std::ostream& os) const } } os << "\n }"; + // contact moment + if( m_momentX != nullptr ) + { + os << axom::fmt::format("\n\tfx: {}", axom::fmt::join(m_momentX, m_momentX+num_verts, ", ")); + os << axom::fmt::format("\n\tfy: {}", axom::fmt::join(m_momentY, m_momentY+num_verts, ", ")); + if(dim == 3) + { + os << axom::fmt::format("\n\tfz: {}", axom::fmt::join(m_momentZ, m_momentZ+num_verts, ", ")); + } + } + os << "\n }"; os << axom::fmt::format("\n elems ({}) {{", num_elem); diff --git a/src/tribol/mesh/MeshData.hpp b/src/tribol/mesh/MeshData.hpp index e7940691..a64c501e 100644 --- a/src/tribol/mesh/MeshData.hpp +++ b/src/tribol/mesh/MeshData.hpp @@ -97,6 +97,67 @@ struct MeshElemData }; +struct MeshParticleData +{ + + /// Default constructor + MeshParticleData( ); + + /// Destructor + ~MeshParticleData() { } + + int m_numParticles; + + ///////////////////////// + // PARTICLE FIELDS // + ///////////////////////// + real * m_particle_gap; ///< scalar particle gap + const real * m_particle_pressure; ///< scalar particle pressure + + bool m_isGapComputed {false}; ///< true if the particle gaps have been computed + bool m_is_particle_gap_set {false}; ///< true if particle gap field is set + bool m_is_particle_pressure_set {false}; ///< true if particle pressure field is set + ///////////////////////// + + bool m_is_velocity_set {false}; ///< true if particle velocities have been registered + bool m_is_particle_displacement_set {false}; ///< true if particle displacements have been registered + bool m_is_particle_response_set {false}; ///< true if the particle responses have been registered + + const real * m_mat_mod; ///< Bulk/Young's modulus for contacting particles or particle with mesh + const real * m_thickness; ///< Volume element thickness associated with each contact face ( not used ) + + real m_penalty_stiffness {0.}; ///< single scalar kinematic penalty stiffness for each mesh + real m_penalty_scale {1.}; ///< scale factor applied to kinematic penalty only + real m_rate_penalty_stiffness {0.}; ///< single scalar rate penalty stiffness for each mesh + real m_rate_percent_stiffness {0.}; ///< rate penalty is percentage of gap penalty + + bool m_is_kinematic_constant_penalty_set {false}; ///< True if single kinematic constant penalty is set + bool m_is_kinematic_particle_penalty_set {false}; ///< True if the particle-wise kinematic penalty is set + bool m_is_rate_constant_penalty_set {false}; ///< True if the constant rate penalty is set + bool m_is_rate_percent_penalty_set {false}; ///< True if the rate percent penalty is set + + bool m_is_element_thickness_set {false}; ///< True if element thickness is set + + /*! + * \brief Checks if the kinematic penalty data is valid + * + * \param [in] pen_enfrc penalty enforcement option struct + * + * \Return True if the kinematic penalty option has valid data + */ + bool isValidKinematicPenalty( PenaltyEnforcementOptions& pen_options ); + + /*! + * \brief Checks if the rate penalty data is valid + * + * \param [in] pen_enfrc to penalty enforcement option struct + * + * \Return True if the rate penalty option has valid data + */ + bool isValidRatePenalty( PenaltyEnforcementOptions& pen_options ); ///< True if rate penalty option is valid + +}; // end of struct MeshParticleData + class MeshData { public: @@ -116,8 +177,9 @@ class MeshData InterfaceElementType m_elementType; ///< Type of interface element in mesh int m_lengthNodalData; ///< Total number of elements in data arrays (used to create new arrays an index in using connectivity ids) int m_numSurfaceNodes; ///< Total number of unique node ids in the surface connectivity (computed from MeshData::sortSurfaceNodeIds() ) + int m_numParticles; ///< Total number of unique node ids in the particle mesh int m_numCells; ///< Total number of SURFACE cells in the mesh - int m_numNodesPerCell; ///< Number of nodes per SURFACE cell based on cell type + int m_numNodesPerCell; ///< Number of nodes per SURFACE cell based on cell type int m_dim; ///< Dimension of mesh int m_meshId; ///< Mesh Id associated with this data bool m_isValid; ///< True if the mesh is valid @@ -134,10 +196,18 @@ class MeshData real * m_forceY; ///< Y-component of nodal forces real * m_forceZ; ///< Z-component of nodal forces + real * m_momentX; ///< X-component of nodal moments + real * m_momentY; ///< Y-component of nodal moments + real * m_momentZ; ///< Z-component of nodal moments + const real* m_velX; ///< X-component of nodal velocity const real* m_velY; ///< Y-component of nodal velocity const real* m_velZ; ///< Z-component of nodal velocity + const real* m_angVelX; ///< X-component of particle angular velocity + const real* m_angVelY; ///< Y-component of particle angular velocity + const real* m_angVelZ; ///< Z-component of particle angular velocity + real* m_nX; ///< X-component of outward unit face normals real* m_nY; ///< Y-component of outward unit face normals real* m_nZ; ///< Z-component of outward unit face normals @@ -147,6 +217,7 @@ class MeshData real* m_cZ; ///< Z-component of vertex averaged cell centroids real* m_faceRadius; ///< Face radius used in low level proximity check + real* m_particleRadius; ///< Particle radius used in low level proximity check real* m_area; ///< Cell areas @@ -158,8 +229,9 @@ class MeshData real* m_node_nY; ///< Y-component of outward unit node normals real* m_node_nZ; ///< Z-component of outward unit node normals - MeshNodalData m_nodalFields; ///< method specific nodal fields - MeshElemData m_elemData; ///< method/enforcement specific element data + MeshNodalData m_nodalFields; ///< method specific nodal fields + MeshElemData m_elemData; ///< method/enforcement specific element data + MeshParticleData m_particleData; ///< method/enforcement specific particle data public: