Question about contact model API with custom property value

Raheem Sterling_22160

Hi

I have introduced a new contact model, the force was calculated as style F_n = Kn*(rij_thisstep-rij_0), rij_thisstep(double) means the element1.position-element2.position at the current time, and rij_0 means the value in the initial state which stored like double rij_0 = m_customvalue[0]. When I generated a packing without any external force(like gravity, particle body force) and internal force, just this new contact model. But with the passage of time, rij_thisstep-rij_0 accumulated from 1e-9 to 0.1. I am not sure if this caused by the claculation in C++. Are there any method to prevent this?

the code is as below, and the information in the debugging process, rij_0 and rij had a little difference as figure below.

if(timeStepData.time >= m_requestedBondTime  && m_BondStatus[0] == 0.0 && nBondSeperation < 0.0){             m_BondStatusDelta[0] += 1.0;             m_InitialVectorIDelta[0] += (element2.position - element1.position).getX(); //ni1_0_x             m_InitialVectorIDelta[1] += (element2.position - element1.position).getY(); //ni1_0_y             m_InitialVectorIDelta[2] += (element2.position - element1.position).getZ(); //ni1_0_z  } if(m_BondStatus[0] == 1.0){             rij_0 = CSimple3DVector(m_InitialVectorI[0], m_InitialVectorI[1], m_InitialVectorI[2]);             rij_thisstep = element2.position - element1.position; }

Thanks

Raheem Sterling_22160

Hi Stephen,

I have edited the question and sorry about your answer is missing, when I debugged the code, rij_0 was a constant vector where the first "if" wouldn't be accessed anymore when the bond was created because bonstatus changes from 0 to 1, and rij_thisstep wasn't a custom value.

Regards!

Stefan Pantaleev_21979

Raheem Sterling_22160 said:

Hi Stephen,

I have edited the question and sorry about your answer is missing, when I debugged the code, rij_0 was a constant vector where the first "if" wouldn't be accessed anymore when the bond was created because bonstatus changes from 0 to 1, and rij_thisstep wasn't a custom value.

Regards!

Hi Raheem,

I deleted my original reply as I noticed you were setting rij_0 only once using Bond status.

The accumulation of rij_thisstep-rij_0 might be due to multi-threading - have you tried debugging on a single CPU thread?

Are we certain that rij_thisstep is not being accessed anywhere outside of this code snippet?

Best regards,

Stefan

Raheem Sterling_22160

Stefan Pantaleev_21979 said:

Hi Raheem,

I deleted my original reply as I noticed you were setting rij_0 only once using Bond status.

The accumulation of rij_thisstep-rij_0 might be due to multi-threading - have you tried debugging on a single CPU thread?

Are we certain that rij_thisstep is not being accessed anywhere outside of this code snippet?

Best regards,

Stefan

Hi Stephen,

I debugged the code with only one processor to avoid the problem you mentioned, because my code works well in the particles line on global x, y, z axis.

And rij_thisstep was only changed here. I have test like predefine double ndistance  = element1.position-element2.position which minus the value first, then assign the ndistance to custom property value and rij_thisstep, but it was still not work, it seems something happeend when the value passing from m_InitialVectorIDelta[0] to m_InitialVector[0].

Then I want to define the rij_thisstep through custom property also, as below, but the return value of rij_thisstep was zero all the time.

if(m_BondStatus[0] == 1.0)     {          const double* m_CurrentRij = contactCustomProperties->getValue(iCURRENT_RIJ);         double* m_CurrentRijDelta = contactCustomProperties->getDelta(iCURRENT_RIJ);             rij_thisstep = CSimple3DVector(m_CurrentRij[0], m_CurrentRij[1], m_CurrentRij[2]);          m_CurrentRijDelta[0] = RIJ.getX();         m_CurrentRijDelta[1] = RIJ.getY();         m_CurrentRijDelta[2] = RIJ.getZ();                  double problemValue = rij_thisstep.length()-rij_0.length();         //Rest custom property value         m_CurrentRijDelta[0] = -m_CurrentRij[0];         m_CurrentRijDelta[1] = -m_CurrentRij[1];         m_CurrentRijDelta[2] = -m_CurrentRij[2]; } 

Best Regards!

Raheem

Stefan Pantaleev_21979

Raheem Sterling_22160 said:

Hi Stephen,

I debugged the code with only one processor to avoid the problem you mentioned, because my code works well in the particles line on global x, y, z axis.

And rij_thisstep was only changed here. I have test like predefine double ndistance  = element1.position-element2.position which minus the value first, then assign the ndistance to custom property value and rij_thisstep, but it was still not work, it seems something happeend when the value passing from m_InitialVectorIDelta[0] to m_InitialVector[0].

Then I want to define the rij_thisstep through custom property also, as below, but the return value of rij_thisstep was zero all the time.

if(m_BondStatus[0] == 1.0)     {          const double* m_CurrentRij = contactCustomProperties->getValue(iCURRENT_RIJ);         double* m_CurrentRijDelta = contactCustomProperties->getDelta(iCURRENT_RIJ);             rij_thisstep = CSimple3DVector(m_CurrentRij[0], m_CurrentRij[1], m_CurrentRij[2]);          m_CurrentRijDelta[0] = RIJ.getX();         m_CurrentRijDelta[1] = RIJ.getY();         m_CurrentRijDelta[2] = RIJ.getZ();                  double problemValue = rij_thisstep.length()-rij_0.length();         //Rest custom property value         m_CurrentRijDelta[0] = -m_CurrentRij[0];         m_CurrentRijDelta[1] = -m_CurrentRij[1];         m_CurrentRijDelta[2] = -m_CurrentRij[2]; } 

Best Regards!

Raheem

Hi Raheem,

The deltas defined in the current timestep are applied in the next timestep.

In the above code it looks like you are setting the delta twice and the final value is the second value -m_CurrentRij[n] so you are effectively zeroing m_CurrentRij every timestep. This is why you are getting the zeros.

Why rij_thisstep is accumulating with the original code is still not clear to me however. I can't see an obvious reason from the code snipped apart from the particles actually separating by that cumulative distance.

Best regards,

Stefan

Raheem Sterling_22160

Stefan Pantaleev_21979 said:

Hi Raheem,

The deltas defined in the current timestep are applied in the next timestep.

In the above code it looks like you are setting the delta twice and the final value is the second value -m_CurrentRij[n] so you are effectively zeroing m_CurrentRij every timestep. This is why you are getting the zeros.

Why rij_thisstep is accumulating with the original code is still not clear to me however. I can't see an obvious reason from the code snipped apart from the particles actually separating by that cumulative distance.

Best regards,

Stefan

Hi Stefan,

I have used the custom property value for both rij_0 and rij_thisstep, and they are exactly the same, I guess it because the calculation in C++ for float, because other values in my code also have the same question. For example:

ni1 = ni1_0.operator*(Orienation_i0_inv).operator*(Orienation_i);

where the Orienation_i0_inv is the matrix at initial state, and orienation_i is the matrix at the current time, even I capped the velocity and angvel of particles, ni1 and ni1_0 was different. 

It has a significant effect on my simulation if I want to reach a stress free state.

Best Regards！

Raheem