Skip to content

Commit 7aecb9b

Browse files
committed
refactor: ModelSBKelvin
1 parent afd5782 commit 7aecb9b

File tree

3 files changed

+82
-85
lines changed

3 files changed

+82
-85
lines changed

Models/SolidBonds/BondModelKelvin/ModelSBKelvin.cpp

Lines changed: 42 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -14,77 +14,79 @@ CModelSBKelvin::CModelSBKelvin()
1414
m_hasGPUSupport = true;
1515
}
1616

17-
void CModelSBKelvin::CalculateSB(double _time, double _timeStep, size_t _iLeft, size_t _iRight, size_t _iBond, SSolidBondStruct& _bonds, unsigned* _pBrokenBondsNum) const
17+
void CModelSBKelvin::CalculateSB(double _time, double _timeStep, size_t _iLeft, size_t _iRight, size_t _iBond, SSolidBondStruct& _bonds, unsigned* _brokenBondsNum) const
1818
{
1919
// relative angle velocity of contact partners
20-
CVector3 relAngleVel = Particles().AnglVel(_iLeft) - Particles().AnglVel(_iRight);
20+
const CVector3 relAngleVel = Particles().AnglVel(_iLeft) - Particles().AnglVel(_iRight);
2121

22-
double dMu = Bonds().Viscosity(_iBond);
22+
const double mu = Bonds().Viscosity(_iBond);
2323
// the bond in the global coordinate system
24-
CVector3 currentBond = GetSolidBond(Particles().Coord(_iRight), Particles().Coord(_iLeft), m_PBC);
25-
double dDistanceBetweenCenters = currentBond.Length();
26-
CVector3 rAC = currentBond*0.5;
24+
const CVector3 currentBond = GetSolidBond(Particles().Coord(_iRight), Particles().Coord(_iLeft), m_PBC);
25+
const double distanceBetweenCenters = currentBond.Length();
26+
const CVector3 rAC = currentBond * 0.5;
2727

2828
// !! do not delete source equation
2929
//relativeVelocity = (m_ParticleVelocities[ nLeftParticleID ]-m_ParticleAnglVel[ nLeftParticleID ]*rAC)
3030
// - ( m_ParticleVelocities[ nRightParticleID ] + m_ParticleAnglVel[ nRightParticleID ]*rAC);
3131

3232
// optimized
33-
CVector3 sumAngleVelocity = Particles().AnglVel(_iLeft) + Particles().AnglVel(_iRight);
34-
CVector3 relativeVelocity = Particles().Vel(_iLeft) - Particles().Vel(_iRight) - sumAngleVelocity*rAC;
33+
const CVector3 sumAngleVelocity = Particles().AnglVel(_iLeft) + Particles().AnglVel(_iRight);
34+
const CVector3 relativeVelocity = Particles().Vel(_iLeft) - Particles().Vel(_iRight) - sumAngleVelocity * rAC;
3535

36-
CVector3 currentContact = currentBond/dDistanceBetweenCenters;
37-
CVector3 tempVector = Bonds().PrevBond(_iBond)*currentBond;
36+
const CVector3 currentContact = currentBond / distanceBetweenCenters;
37+
const CVector3 tempVector = Bonds().PrevBond(_iBond) * currentBond;
3838

39-
CVector3 Phi = currentContact*(DotProduct(sumAngleVelocity, currentContact)*_timeStep*0.5);
39+
const CVector3 phi = currentContact * (DotProduct(sumAngleVelocity, currentContact) * _timeStep * 0.5);
4040

41-
CMatrix3 M( 1 + tempVector.z*Phi.z + tempVector.y*Phi.y, Phi.z - tempVector.z - tempVector.y*Phi.x, -Phi.y - tempVector.z*Phi.x + tempVector.y,
42-
tempVector.z - Phi.z - tempVector.x*Phi.y, tempVector.z*Phi.z + 1 + tempVector.x*Phi.x, -tempVector.z*Phi.y + Phi.x - tempVector.x,
43-
-tempVector.y - tempVector.x*Phi.z + Phi.y, -tempVector.y*Phi.z + tempVector.x - Phi.x, tempVector.y*Phi.y + tempVector.x*Phi.x + 1);
41+
const CMatrix3 M(1 + tempVector.z * phi.z + tempVector.y * phi.y, phi.z - tempVector.z - tempVector.y * phi.x, -phi.y - tempVector.z * phi.x + tempVector.y,
42+
tempVector.z - phi.z - tempVector.x * phi.y, tempVector.z * phi.z + 1 + tempVector.x * phi.x, -tempVector.z * phi.y + phi.x - tempVector.x,
43+
-tempVector.y - tempVector.x * phi.z + phi.y, -tempVector.y * phi.z + tempVector.x - phi.x, tempVector.y * phi.y + tempVector.x * phi.x + 1);
4444

45-
CVector3 normalVelocity = currentContact * DotProduct(currentContact, relativeVelocity);
46-
CVector3 tangentialVelocity = relativeVelocity - normalVelocity;
45+
const CVector3 normalVelocity = currentContact * DotProduct(currentContact, relativeVelocity);
46+
const CVector3 tangentialVelocity = relativeVelocity - normalVelocity;
4747

4848
// normal angle velocity
49-
CVector3 normalAngleVel = currentContact*DotProduct(currentContact, relAngleVel);
50-
CVector3 tangAngleVel = relAngleVel - normalAngleVel;
49+
const CVector3 normalAngleVel = currentContact * DotProduct(currentContact, relAngleVel);
50+
const CVector3 tangAngleVel = relAngleVel - normalAngleVel;
5151

5252
// calculate the force
53-
double dStrainTotal = (dDistanceBetweenCenters-Bonds().InitialLength(_iBond)) / Bonds().InitialLength(_iBond);
53+
const double strainTotal = (distanceBetweenCenters - Bonds().InitialLength(_iBond)) / Bonds().InitialLength(_iBond);
5454

55-
CVector3 vNormalForce = currentContact*(-1*Bonds().CrossCut(_iBond)*Bonds().NormalStiffness(_iBond)*dStrainTotal);
55+
const CVector3 normalForce = currentContact * (-1 * Bonds().CrossCut(_iBond) * Bonds().NormalStiffness(_iBond) * strainTotal);
5656

57-
CVector3 vDampingForce = -dMu *normalVelocity*Bonds().CrossCut(_iBond)*Bonds().NormalStiffness(_iBond)*fabs(dStrainTotal);
58-
if (vDampingForce.Length() > 0.5* vNormalForce.Length())
59-
vDampingForce *= 0.5*vNormalForce.Length() / vDampingForce.Length();
57+
CVector3 dampingForceNorm = -mu * normalVelocity * Bonds().CrossCut(_iBond) * Bonds().NormalStiffness(_iBond) * fabs(strainTotal);
58+
if (dampingForceNorm.Length() > 0.5 * normalForce.Length())
59+
dampingForceNorm *= 0.5 * normalForce.Length() / dampingForceNorm.Length();
6060

61-
_bonds.TangentialOverlap(_iBond) = M* Bonds().TangentialOverlap(_iBond) - tangentialVelocity*_timeStep;
62-
_bonds.TangentialForce(_iBond) = Bonds().TangentialOverlap(_iBond)*(Bonds().TangentialStiffness(_iBond)*Bonds().CrossCut(_iBond) / Bonds().InitialLength(_iBond));
61+
_bonds.TangentialOverlap(_iBond) = M * Bonds().TangentialOverlap(_iBond) - tangentialVelocity * _timeStep;
62+
_bonds.TangentialForce(_iBond) = Bonds().TangentialOverlap(_iBond) * (Bonds().TangentialStiffness(_iBond) * Bonds().CrossCut(_iBond) / Bonds().InitialLength(_iBond));
6363

64-
CVector3 vDampingTangForce = -dMu * tangentialVelocity*Bonds().TangentialOverlap(_iBond).Length()*(Bonds().TangentialStiffness(_iBond)*Bonds().CrossCut(_iBond) / Bonds().InitialLength(_iBond));
65-
if (vDampingTangForce.Length() > 0.5*Bonds().TangentialForce(_iBond).Length())
66-
vDampingTangForce *= 0.5*Bonds().TangentialForce(_iBond).Length() / vDampingTangForce.Length();
64+
CVector3 dampingForceTang = -mu * tangentialVelocity * Bonds().TangentialOverlap(_iBond).Length() * (Bonds().TangentialStiffness(_iBond) * Bonds().CrossCut(_iBond) / Bonds().InitialLength(_iBond));
65+
if (dampingForceTang.Length() > 0.5 * Bonds().TangentialForce(_iBond).Length())
66+
dampingForceTang *= 0.5 * Bonds().TangentialForce(_iBond).Length() / dampingForceTang.Length();
6767

68-
_bonds.NormalMoment(_iBond) = M* Bonds().NormalMoment(_iBond)- normalAngleVel*(_timeStep* 2 * Bonds().AxialMoment(_iBond)*Bonds().TangentialStiffness(_iBond)/ Bonds().InitialLength(_iBond));
69-
_bonds.TangentialMoment(_iBond) = M* Bonds().TangentialMoment(_iBond) - tangAngleVel*(_timeStep*Bonds().NormalStiffness(_iBond)*Bonds().AxialMoment(_iBond)/ Bonds().InitialLength(_iBond));
70-
_bonds.TotalForce(_iBond) = vNormalForce + Bonds().TangentialForce(_iBond) + vDampingForce + vDampingTangForce;
68+
_bonds.NormalMoment(_iBond) = M * Bonds().NormalMoment(_iBond) - normalAngleVel * (_timeStep * 2 * Bonds().AxialMoment(_iBond) * Bonds().TangentialStiffness(_iBond) / Bonds().InitialLength(_iBond));
69+
_bonds.TangentialMoment(_iBond) = M * Bonds().TangentialMoment(_iBond) - tangAngleVel * (_timeStep * Bonds().NormalStiffness(_iBond) * Bonds().AxialMoment(_iBond) / Bonds().InitialLength(_iBond));
70+
_bonds.TotalForce(_iBond) = normalForce + Bonds().TangentialForce(_iBond) + dampingForceNorm + dampingForceTang;
7171

72-
_bonds.UnsymMoment(_iBond) = currentBond * 0.5*Bonds().TangentialForce(_iBond);
72+
_bonds.UnsymMoment(_iBond) = currentBond * 0.5 * Bonds().TangentialForce(_iBond);
7373
_bonds.PrevBond(_iBond) = currentBond;
7474

75-
if (m_parameters[0].value == 0 ) return; // consider breakage
75+
if (m_parameters[0].value == 0.0) return; // consider breakage
76+
7677
// check the bond destruction
77-
double dForceLength = vNormalForce.Length();
78-
if (dStrainTotal <= 0) // compression
79-
dForceLength *= -1;
80-
double dMaxStress = dForceLength / Bonds().CrossCut(_iBond) + Bonds().TangentialMoment(_iBond).Length()*Bonds().Diameter(_iBond) / (2 * Bonds().AxialMoment(_iBond));
81-
double dMaxTorque = Bonds().TangentialForce(_iBond).Length() / Bonds().CrossCut(_iBond) + Bonds().NormalMoment(_iBond).Length()*Bonds().Diameter(_iBond) / (2 * 2 * Bonds().AxialMoment(_iBond));
78+
double forceLength = normalForce.Length();
79+
if (strainTotal <= 0) // compression
80+
forceLength *= -1;
81+
82+
const double maxStress = forceLength / Bonds().CrossCut(_iBond) + Bonds().TangentialMoment(_iBond).Length() * Bonds().Diameter(_iBond) / (2 * Bonds().AxialMoment(_iBond));
83+
const double maxTorque = Bonds().TangentialForce(_iBond).Length() / Bonds().CrossCut(_iBond) + Bonds().NormalMoment(_iBond).Length() * Bonds().Diameter(_iBond) / (2 * 2 * Bonds().AxialMoment(_iBond));
8284

83-
if ( ( dMaxStress >= Bonds().NormalStrength(_iBond) ) || ( dMaxTorque >= Bonds().TangentialStrength(_iBond)) )
85+
if (maxStress >= Bonds().NormalStrength(_iBond) || maxTorque >= Bonds().TangentialStrength(_iBond))
8486
{
8587
_bonds.Active(_iBond) = false;
8688
_bonds.EndActivity(_iBond) = _time;
87-
*_pBrokenBondsNum += 1;
89+
*_brokenBondsNum += 1;
8890
}
8991
}
9092

Models/SolidBonds/BondModelKelvin/ModelSBKelvin.cu

Lines changed: 39 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -87,72 +87,67 @@ void __global__ CUDA_CalcSBForce_Kelvin_kernel(
8787
if (!_bondActivities[i]) continue;
8888

8989
// relative angle velocity of contact partners
90-
CVector3 relAngleVel = _partAnglVels[_bondLeftIDs[i]] - _partAnglVels[_bondRightIDs[i]];
90+
const CVector3 relAngleVel = _partAnglVels[_bondLeftIDs[i]] - _partAnglVels[_bondRightIDs[i]];
9191

9292
// the bond in the global coordinate system
93-
CVector3 currentBond = GetSolidBond(_partCoords[_bondRightIDs[i]], _partCoords[_bondLeftIDs[i]], PBC);
94-
double dDistanceBetweenCenters = currentBond.Length();
95-
96-
double dBondInitLength = _bondInitialLengths[i];
93+
const CVector3 currentBond = GetSolidBond(_partCoords[_bondRightIDs[i]], _partCoords[_bondLeftIDs[i]], PBC);
94+
double distanceBetweenCenters = currentBond.Length();
9795

9896
// optimized
99-
CVector3 sumAngleVelocity = _partAnglVels[_bondLeftIDs[i]] + _partAnglVels[_bondRightIDs[i]];
100-
CVector3 relativeVelocity = _partVels[_bondLeftIDs[i]] - _partVels[_bondRightIDs[i]] - sumAngleVelocity*currentBond*0.5;
97+
const CVector3 sumAngleVelocity = _partAnglVels[_bondLeftIDs[i]] + _partAnglVels[_bondRightIDs[i]];
98+
const CVector3 relativeVelocity = _partVels[_bondLeftIDs[i]] - _partVels[_bondRightIDs[i]] - sumAngleVelocity * currentBond * 0.5;
10199

102-
CVector3 currentContact = currentBond / dDistanceBetweenCenters;
103-
CVector3 tempVector = _bondPrevBonds[i] * currentBond;
100+
const CVector3 currentContact = currentBond / distanceBetweenCenters;
101+
const CVector3 tempVector = _bondPrevBonds[i] * currentBond;
104102

105-
CVector3 Phi = currentContact*(DotProduct(sumAngleVelocity, currentContact)*_timeStep*0.5);
103+
const CVector3 phi = currentContact * (DotProduct(sumAngleVelocity, currentContact) * _timeStep * 0.5);
106104

107-
CMatrix3 M( 1 + tempVector.z*Phi.z + tempVector.y*Phi.y, Phi.z - tempVector.z - tempVector.y*Phi.x, -Phi.y - tempVector.z*Phi.x + tempVector.y,
108-
tempVector.z - Phi.z - tempVector.x*Phi.y, tempVector.z*Phi.z + 1 + tempVector.x*Phi.x, -tempVector.z*Phi.y + Phi.x - tempVector.x,
109-
-tempVector.y - tempVector.x*Phi.z + Phi.y, -tempVector.y*Phi.z + tempVector.x - Phi.x, tempVector.y*Phi.y + tempVector.x*Phi.x + 1);
105+
const CMatrix3 M(1 + tempVector.z * phi.z + tempVector.y * phi.y, phi.z - tempVector.z - tempVector.y * phi.x, -phi.y - tempVector.z * phi.x + tempVector.y,
106+
tempVector.z - phi.z - tempVector.x * phi.y, tempVector.z * phi.z + 1 + tempVector.x * phi.x, -tempVector.z * phi.y + phi.x - tempVector.x,
107+
-tempVector.y - tempVector.x * phi.z + phi.y, -tempVector.y * phi.z + tempVector.x - phi.x, tempVector.y * phi.y + tempVector.x * phi.x + 1);
110108

111-
CVector3 normalVelocity = currentContact * DotProduct(currentContact, relativeVelocity);
112-
CVector3 tangentialVelocity = relativeVelocity - normalVelocity;
109+
const CVector3 normalVelocity = currentContact * DotProduct(currentContact, relativeVelocity);
110+
const CVector3 tangentialVelocity = relativeVelocity - normalVelocity;
113111

114112
// normal angle velocity
115-
CVector3 normalAngleVel = currentContact*DotProduct(currentContact, relAngleVel);
116-
CVector3 tangAngleVel = relAngleVel - normalAngleVel;
113+
const CVector3 normalAngleVel = currentContact * DotProduct(currentContact, relAngleVel);
114+
const CVector3 tangAngleVel = relAngleVel - normalAngleVel;
117115

118116
// calculate the force
119-
double dStrainTotal = (dDistanceBetweenCenters-dBondInitLength) / dBondInitLength;
120-
CVector3 vNormalForce = currentContact * (-1*_bondCrossCuts[i] * _bondNormalStiffnesses[i] * dStrainTotal);
121-
122-
double dMu = _bondViscosities[i];
117+
const double strainTotal = (distanceBetweenCenters - _bondInitialLengths[i]) / _bondInitialLengths[i];
118+
const CVector3 normalForce = currentContact * (-1 * _bondCrossCuts[i] * _bondNormalStiffnesses[i] * strainTotal);
123119

124-
CVector3 vDampingForce = -dMu * normalVelocity*_bondCrossCuts[i] * _bondNormalStiffnesses[i] *fabs(dStrainTotal);
125-
if (vDampingForce.Length() > 0.5*vNormalForce.Length())
126-
vDampingForce *= 0.5*vNormalForce.Length() / vDampingForce.Length();
120+
CVector3 dampingForceNorm = -_bondViscosities[i] * normalVelocity * _bondCrossCuts[i] * _bondNormalStiffnesses[i] * fabs(strainTotal);
121+
if (dampingForceNorm.Length() > 0.5 * normalForce.Length())
122+
dampingForceNorm *= 0.5 * normalForce.Length() / dampingForceNorm.Length();
127123

128124
_bondTangentialOverlaps[i] = M * _bondTangentialOverlaps[i] - tangentialVelocity * _timeStep;
129-
const CVector3 vTangentialForce = _bondTangentialOverlaps[i] * (_bondTangentialStiffnesses[i] * _bondCrossCuts[i] / dBondInitLength);
130-
131-
CVector3 vDampingTangForce = -dMu * tangentialVelocity*_bondTangentialOverlaps[i].Length() *(_bondTangentialStiffnesses[i] * _bondCrossCuts[i] / dBondInitLength);
132-
if (vDampingTangForce.Length() > 0.5*vTangentialForce.Length())
133-
vDampingTangForce *= 0.5*vTangentialForce.Length() / vDampingTangForce.Length();
125+
const CVector3 tangentialForce = _bondTangentialOverlaps[i] * (_bondTangentialStiffnesses[i] * _bondCrossCuts[i] / _bondInitialLengths[i]);
134126

135-
const CVector3 vBondNormalMoment = M * _bondNormalMoments[i] - normalAngleVel * (_timeStep * 2 * _bondAxialMoments[i] * _bondTangentialStiffnesses[i] / dBondInitLength);
136-
const CVector3 vBondTangentialMoment = M * _bondTangentialMoments[i] - tangAngleVel * (_timeStep * _bondNormalStiffnesses[i] * _bondAxialMoments[i] / dBondInitLength);
127+
CVector3 dampingForceTang = -_bondViscosities[i] * tangentialVelocity * _bondTangentialOverlaps[i].Length() * (_bondTangentialStiffnesses[i] * _bondCrossCuts[i] / _bondInitialLengths[i]);
128+
if (dampingForceTang.Length() > 0.5 * tangentialForce.Length())
129+
dampingForceTang *= 0.5 * tangentialForce.Length() / dampingForceTang.Length();
137130

138-
_bondNormalMoments[i] = vBondNormalMoment;
139-
_bondTangentialMoments[i] = vBondTangentialMoment;
131+
const CVector3 bondNormalMoment = M * _bondNormalMoments[i] - normalAngleVel * (_timeStep * 2 * _bondAxialMoments[i] * _bondTangentialStiffnesses[i] / _bondInitialLengths[i]);
132+
const CVector3 bondTangentialMoment = M * _bondTangentialMoments[i] - tangAngleVel * (_timeStep * _bondNormalStiffnesses[i] * _bondAxialMoments[i] / _bondInitialLengths[i]);
133+
const CVector3 bondUnsymMoment = currentBond * 0.5 * tangentialForce;
140134

141-
const CVector3 vUnsymMoment = currentBond*0.5 * vTangentialForce;
135+
_bondNormalMoments[i] = bondNormalMoment;
136+
_bondTangentialMoments[i] = bondTangentialMoment;
142137
_bondPrevBonds[i] = currentBond;
143-
_bondTotalForces[i] = vNormalForce + vTangentialForce + vDampingForce + vDampingTangForce;
138+
_bondTotalForces[i] = normalForce + tangentialForce + dampingForceNorm + dampingForceTang;
144139

145140
if (m_vConstantModelParameters[0])
146141
{
147-
double dForceLength = vNormalForce.Length();
148-
if (dStrainTotal <= 0) // compression
149-
dForceLength *= -1;
142+
double forceLength = normalForce.Length();
143+
if (strainTotal <= 0) // compression
144+
forceLength *= -1;
150145

151146
// check the bond destruction
152-
double dMaxStress = dForceLength / _bondCrossCuts[i] + vBondTangentialMoment.Length() * _bondDiameters[i] / (2 * _bondAxialMoments[i]);
153-
double dMaxTorque = vTangentialForce.Length() / _bondCrossCuts[i] + vBondNormalMoment.Length() * _bondDiameters[i] / (4 * _bondAxialMoments[i]);
147+
const double maxStress = forceLength / _bondCrossCuts[i] + bondTangentialMoment.Length() * _bondDiameters[i] / (2 * _bondAxialMoments[i]);
148+
const double maxTorque = tangentialForce.Length() / _bondCrossCuts[i] + bondNormalMoment.Length() * _bondDiameters[i] / (4 * _bondAxialMoments[i]);
154149

155-
if (dMaxStress >= _bondNormalStrengths[i] || dMaxTorque >= _bondTangentialStrengths[i])
150+
if (maxStress >= _bondNormalStrengths[i] || maxTorque >= _bondTangentialStrengths[i])
156151
{
157152
_bondActivities[i] = false;
158153
_bondEndActivities[i] = _time;
@@ -161,9 +156,9 @@ void __global__ CUDA_CalcSBForce_Kelvin_kernel(
161156
}
162157

163158
// apply forces and moments directly to particles, only if bond is not broken
164-
const CVector3 partForce = vNormalForce + vTangentialForce + vDampingForce + vDampingTangForce;
165-
const CVector3 partMoment1 = vBondNormalMoment + vBondTangentialMoment - vUnsymMoment;
166-
const CVector3 partMoment2 = vBondNormalMoment + vBondTangentialMoment + vUnsymMoment;
159+
const CVector3 partForce = normalForce + tangentialForce + dampingForceNorm + dampingForceTang;
160+
const CVector3 partMoment1 = bondNormalMoment + bondTangentialMoment - bondUnsymMoment;
161+
const CVector3 partMoment2 = bondNormalMoment + bondTangentialMoment + bondUnsymMoment;
167162
CUDA_VECTOR3_ATOMIC_ADD(_partForces[_bondLeftIDs[i]], partForce);
168163
CUDA_VECTOR3_ATOMIC_ADD(_partMoments[_bondLeftIDs[i]], partMoment1);
169164
CUDA_VECTOR3_ATOMIC_SUB(_partForces[_bondRightIDs[i]], partForce);

Models/SolidBonds/BondModelKelvin/ModelSBKelvin.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@ class CModelSBKelvin : public CSolidBondModel
1111
public:
1212
CModelSBKelvin();
1313

14-
void CalculateSB(double _time, double _timeStep, size_t _iLeft, size_t _iRight, size_t _iBond, SSolidBondStruct& _bonds, unsigned* _pBrokenBondsNum) const override;
14+
void CalculateSB(double _time, double _timeStep, size_t _iLeft, size_t _iRight, size_t _iBond, SSolidBondStruct& _bonds, unsigned* _brokenBondsNum) const override;
1515
void ConsolidatePart(double _time, double _timeStep, size_t _iBond, size_t _iPart, SParticleStruct& _particles) const override;
1616

1717
void SetParametersGPU(const std::vector<double>& _parameters, const SPBC& _pbc) override;

0 commit comments

Comments
 (0)