Hi,
I am wondering if there is a way to extract each individual contact force acting on a sphere clump (a multisphere)? I know there is a method called pair/gran/local that is able to get per-sphere contact forces. However this command does not offer an option indicating the id_multisphere that each specific sphere belongs to and computing stress through contact force and branch vector can be very difficult since we need identify the center of mass coordinates of two multisphere bodies that are in contact. Any suggestion is greatly appreciated!
Best,
Liuchi
ckloss | Mon, 01/02/2017 - 13:45
Hi Liuchi,
Hi Liuchi,
unfortunately such an option is not yet available. However, adding the multisphere ID to compute pair/gran/local would be trivial. I'll add that to the "nice-to-have list" for one of the next releases
Christoph
richti83 | Wed, 01/04/2017 - 13:48
I've allready made a little
I've allready made a little modification to compute_pair_gran local to sum all pair-forces acting on one multisphere:
//line 62
#include "fix_multisphere.h"
//line 126
fix_ms_ = static_cast(modify->find_fix_style_strict("multisphere", 0));
if (fix_ms_)
{
ms_ = &fix_ms_->data();
printf("Detected multisphere\n");
}
else error->warning(FLERR, "Did not find fix multisphere");
//line 477
if(fflag)
{
array[ipair][n++] = fx;
array[ipair][n++] = fy;
array[ipair][n++] = fz;
double frc[3] = { -1 * fx,-1 * fy,-1 * fz };
const int ibodyi = fix_ms_ ? ((fix_ms_->belongs_to(i) > -1) ? (ms_->map(fix_ms_->belongs_to(i))) : -1) : -1;
if (ibodyi >= 0)
ms_->add_pf(frc, ibodyi);
const int ibodyj = fix_ms_ ? ((fix_ms_->belongs_to(j) > -1) ? (ms_->map(fix_ms_->belongs_to(j))) : -1) : -1;
if (ibodyj >= 0)
ms_->add_pf(frc, ibodyj);
}
comute_pair_gran_local.h:
class FixMultisphere* fix_ms_;
class Multisphere *ms_;
multisphere.cpp
//line 101
temp_(*customValues_.addElementProperty< ScalarContainer >("temp","comm_exchange_borders","frame_invariant","restart_yes")),
temp_old_(*customValues_.addElementProperty< ScalarContainer >("temp_old","comm_exchange_borders","frame_invariant","restart_yes")),
pf_(*customValues_.addElementProperty< VectorContainer >("pf", "comm_exchange_borders", "frame_invariant", "restart_yes")) //new
{
and in multisphere.h:
inline void add_pf(double *frc, int ibody_local) {
vectorAdd3D(pf_(ibody_local), frc, pf_(ibody_local));
//printf("add_pf %d (%f %f %f) at step %ld\n", ibody_local, frc[0], frc[1], frc[2],update->ntimestep);
}
the value can be accesed with compute rigid out of the box
compute fc all pair/gran/local force
...
compute bid all rigid property id_multisphere # ID
compute btype all rigid property clumptype # type
compute xcm all rigid property xcm # center of mass ("COM")
compute fcm all rigid property fcm # force at COM
compute quat all rigid property quat # rotation arround COM
compute vel all rigid property vcm # velocity of COM
compute pf all rigid property pf # pair forces
..
dump rigids all local 1000 post/rigids*.dump c_xcm[1] c_xcm[2] c_xcm[3] c_quat[1] c_quat[2] c_quat[3] c_quat[4] c_vel[1] c_vel[2] c_vel[3] c_pf[1] c_pf[2] c_pf[3] c_fcm[1] c_fcm[2] c_fcm[3]
dump dmp all custom 1000 post/dump*.liggghts id type mol x y z vx vy vz fx fy fz omegax omegay omegaz radius