| std::enable_if< I< sizeof...(BF), void >::type recursiveBuildStress(const double &fij, const Vector3d &ra, const Vector3d &rA, const Vector3d &rb, const Vector3d &rB, const Vector3d &rab, const Vector3d &rAB, const int &i_gridPoint, const int &i_stress, std::tuple< Stress< BF, stressType > &... > t){ if(I==i_stress &&stressType==Cauchy) { assert(rab.squaredNorm()>epsilon);std::get< I >(t).field[i_gridPoint]=std::get< I >(t).field[i_gridPoint]+std::get< I >(t).method.bondFunction(ra, rb) *fij *rab.transpose() *rab/rab.norm();} else if(I==i_stress &&stressType==Piola) { assert(rab.squaredNorm()>epsilon);std::get< I >(t).field[i_gridPoint]=std::get< I >(t).field[i_gridPoint]+std::get< I >(t).method.bondFunction(rA, rB) *fij *rAB.transpose() *rab/rab.norm();} else recursiveBuildStress< I+1 >(fij, ra, rA, rb, rB, rab, rAB, i_gridPoint, i_stress, t);}template< std::size_t I=0, typename ...TStress >inline typename std::enable_if< I==sizeof...(TStress), void >::typerecursiveBuildContinuumFields(const double &mass, const Vector3d &velocity, const Vector3d &position, const int &i_gridPoint, const int &i_stress, std::tuple< TStress &... > t){ if(sizeof...(TStress)!=0) assert(0);}template< std::size_t I=0, typename ...BF >inline typename std::enable_if< I< sizeof...(BF), void >::typerecursiveBuildContinuumFields(const double &mass, const Vector3d &velocity, const Vector3d &position, const int &i_gridPoint, const int &i_stress, std::tuple< Stress< BF, Cauchy > &... > t){ if(I==i_stress) { auto &stress=std::get< I >(t);if(position.norm() >=stress.method.getAveragingDomainSize()) return;double weight=stress.method(position);stress.massDensityField[i_gridPoint]+=mass *weight;stress.momentumDensityField[i_gridPoint]+=mass *weight *velocity;} else recursiveBuildContinuumFields< I+1 >(mass, velocity, position, i_gridPoint, i_stress, t);}template< std::size_t I=0, typename ...TStress >inline typename std::enable_if< I==sizeof...(TStress), void >::typerecursiveFinalizeContinuumVelocity(const int &i_gridPoint, const int &i_stress, std::tuple< TStress &... > t){ if(sizeof...(TStress)!=0) assert(0);}template< std::size_t I=0, typename ...BF >inline typename std::enable_if< I< sizeof...(BF), void >::typerecursiveFinalizeContinuumVelocity(const int &i_gridPoint, const int &i_stress, std::tuple< Stress< BF, Cauchy > &... > t){ if(I==i_stress) { auto &stress=std::get< I >(t);if(stress.massDensityField[i_gridPoint] > epsilon) stress.velocityField[i_gridPoint]=stress.momentumDensityField[i_gridPoint]/stress.massDensityField[i_gridPoint];else stress.velocityField[i_gridPoint]=Vector3d::Zero();} else recursiveFinalizeContinuumVelocity< I+1 >(i_gridPoint, i_stress, t);}template< std::size_t I=0, typename ...TStress >inline typename std::enable_if< I==sizeof...(TStress), Vector3d >::typerecursiveGetContinuumVelocity(const int &i_gridPoint, const int &i_stress, std::tuple< TStress &... > t){ if(sizeof...(TStress)!=0) assert(0);return Vector3d::Zero();}template< std::size_t I=0, typename ...BF >inline typename std::enable_if< I< sizeof...(BF), Vector3d >::typerecursiveGetContinuumVelocity(const int &i_gridPoint, const int &i_stress, std::tuple< Stress< BF, Cauchy > &... > t){ if(I==i_stress) return std::get< I >(t).velocityField[i_gridPoint];return recursiveGetContinuumVelocity< I+1 >(i_gridPoint, i_stress, t);}template< std::size_t I=0, typename ...TStress >inline typename std::enable_if< I==sizeof...(TStress), void >::typerecursiveBuildKineticStress(const double &mass, const Vector3d &velocity, const Vector3d &position, const int &i_gridPoint, const int &i_stress, std::tuple< TStress &... > t){ if(sizeof...(TStress)!=0) assert(0);}template< std::size_t I=0, typename ...BF >inline typename std::enable_if< I< sizeof...(BF), void >::typerecursiveBuildKineticStress(const double &mass, const Vector3d &velocity, const Vector3d &position, const int &i_gridPoint, const int &i_stress, std::tuple< Stress< BF, Cauchy > &... > t){ if(I==i_stress) { auto &stress=std::get< I >(t);if(position.norm() >=stress.method.getAveragingDomainSize()) return;stress.field[i_gridPoint]=stress.field[i_gridPoint] - amuAngstromSquaredPerPicosecondSquaredToEv *stress.method(position) *mass *velocity.transpose() *velocity;} else recursiveBuildKineticStress< I+1 >(mass, velocity, position, i_gridPoint, i_stress, t);}template< std::size_t I=0, typename ...TStress >inline typename std::enable_if< I==sizeof...(TStress), void >::typerecursiveNullifyStress(std::tuple< TStress &... > t){}template< std::size_t I=0, StressType stressType, typename ...BF >inline typename std::enable_if< I< sizeof...(BF), void >::typerecursiveNullifyStress(std::tuple< Stress< BF, stressType > &... > t){ auto &stress=std::get< I >(t);std::fill(stress.field.begin(), stress.field.end(), Matrix3d::Zero());if constexpr(stressType==Cauchy) { std::fill(stress.momentumDensityField.begin(), stress.momentumDensityField.end(), Vector3d::Zero());std::fill(stress.massDensityField.begin(), stress.massDensityField.end(), 0.0);std::fill(stress.velocityField.begin(), stress.velocityField.end(), Vector3d::Zero());} recursiveNullifyStress< I+1 >(t);}double averagingDomainSize_max(const std::tuple<> t){ return 0;}template< std::size_t I=0, typename ...TStress >inline typename std::enable_if< I==sizeof...(TStress) -1, double >::type averagingDomainSize_max(const std::tuple< TStress &... > t){ return std::get< I >(t).method.getAveragingDomainSize();}template< std::size_t I=0, typename ...TStress >inline typename std::enable_if< I< sizeof...(TStress) -1, double >::type averagingDomainSize_max(const std::tuple< TStress &... > t){ return std::max(std::get< I >(t).method.getAveragingDomainSize(), averagingDomainSize_max< I+1 >(t));}std::map< Grid< Reference > *, double > recursiveGridMaxAveragingDomainSizeMap(const std::tuple<> &t){std::map< Grid< Reference > *, double > map;return map;}template< std::size_t I=0, StressType stressType, typename TGrid=typename std::conditional< stressType==Piola, Grid< Reference >, Grid< Current > >::type, typename ...BF >inline typename std::enable_if< I==sizeof...(BF), std::map< TGrid *, double > >::type recursiveGridMaxAveragingDomainSizeMap(const std::tuple< Stress< BF, stressType > &... > t){ std::map< TGrid *, double > map;return map;}template< std::size_t I=0, StressType stressType, typename TGrid=typename std::conditional< stressType==Piola, Grid< Reference >, Grid< Current > >::type, typename ...BF >inline typename std::enable_if< I< sizeof...(BF), std::map< TGrid *, double > >::type recursiveGridMaxAveragingDomainSizeMap(const std::tuple< Stress< BF, stressType > &... > t){ std::map< TGrid *, double > map;const std::map< TGrid *, double > &rmap=recursiveGridMaxAveragingDomainSizeMap< I+1 >(t);map[std::get< I >(t).pgrid]=std::get< I >(t).method.getAveragingDomainSize();for(const auto &pair :rmap) if(std::get< I >(t).pgrid==pair.first) { map[pair.first]=std::max(map.at(std::get< I >(t).pgrid), pair.second);} else { auto result=map.insert(pair);assert(result.second==true);} return map;}inline std::vector< std::pair< Grid< Current > *, double > > getTGridDomainSizePairs(const std::tuple<> &emptyTuple){ std::vector< std::pair< Grid< Current > *, double > > emptyVectorPair;return emptyVectorPair;}inline std::vector< std::pair< Grid< Reference > *, double > > getTGridDomainSizePairs(std::tuple<> &&emptyTuple){ std::vector< std::pair< Grid< Reference > *, double > > emptyVectorPair;return emptyVectorPair;}template< std::size_t I=0, StressType stressType, typename TGrid=typename std::conditional< stressType==Piola, Grid< Reference >, Grid< Current > >::type, typename ...BF >inline typename std::enable_if< I==sizeof...(BF) -1, std::vector< std::pair< TGrid *, double > > >::type getTGridDomainSizePairs(const std::tuple< Stress< BF, stressType > &... > t){ std::vector< std::pair< TGrid *, double > > vectorPair;vectorPair.push_back({std::get< I >(t).pgrid, std::get< I >(t).method.getAveragingDomainSize()});return vectorPair;}template< std::size_t I=0, StressType stressType, typename TGrid=typename std::conditional< stressType==Piola, Grid< Reference >, Grid< Current > >::type, typename ...BF >inline typename std::enable_if< I< sizeof...(BF) -1, std::vector< std::pair< TGrid *, double > > >::type getTGridDomainSizePairs(const std::tuple< Stress< BF, stressType > &... > t){ std::vector< std::pair< TGrid *, double > > vectorPair;vectorPair.push_back({std::get< I >(t).pgrid, std::get< I >(t).method.getAveragingDomainSize()});std::vector< std::pair< TGrid *, double > > next=getTGridDomainSizePairs< I+1 >(t);vectorPair.insert(vectorPair.end(), next.begin(), next.end());return vectorPair;}std::vector< GridBase * > getBaseGridList(const std::tuple<> t){ std::vector< GridBase * > pgridBaseVector;return pgridBaseVector;}template< std::size_t I=0, typename ...TStress >inline typename std::enable_if< I==sizeof...(TStress) -1, std::vector< GridBase * > >::type getBaseGridList(const std::tuple< TStress &... > t){ std::vector< GridBase * > pgridBaseVector;pgridBaseVector.push_back(std::get< I >(t).pgrid);return pgridBaseVector;}template< std::size_t I=0, typename ...TStress > std::enable_if< I< sizeof...(TStress) -1, std::vector< GridBase * > >::type getBaseGridList(const std::tuple< TStress &... > t){ std::vector< GridBase * > pgridBaseVector;pgridBaseVector.push_back(std::get< I >(t).pgrid);std::vector< GridBase * > | next = getBaseGridList<I+1>(t) |
