diff --git a/src/BODY/pair_body_rounded_polygon.cpp b/src/BODY/pair_body_rounded_polygon.cpp
index 9f9cdea31b1abd7e5c7b21e826aa64a3d6db6355..62a6d77a01b94d57cc34fc06e6827df3080fecaf 100644
--- a/src/BODY/pair_body_rounded_polygon.cpp
+++ b/src/BODY/pair_body_rounded_polygon.cpp
@@ -848,6 +848,7 @@ int PairBodyRoundedPolygon::vertex_against_edge(int i, int j,
// done with the edges from body j,
// given that vertex ni interacts with only one vertex from one edge of body j
+ // comment out this break to take into account concave shapes
// break;
@@ -955,6 +956,7 @@ int PairBodyRoundedPolygon::vertex_against_edge(int i, int j,
// done with the edges from body j,
// given that vertex ni interacts with only one edge from body j
+ // comment out this break to take into account concave shapes
// break;
@@ -1080,12 +1082,12 @@ int PairBodyRoundedPolygon::compute_distance_to_vertex(int ibody,
// check if x0 (the queried vertex) and xmi (the body's center of mass)
// are on the different sides of the edge
- int m = 1;//opposite_sides(xi1, xi2, x0, xmi);
+ int m = opposite_sides(xi1, xi2, x0, xmi);
if (m == 0) {
// x0 and xmi are on not the opposite sides of the edge
- // leave xpi for another edge to detect -- for convex shapes only
+ // leave xpi for another edge to detect
mode = NONE;
diff --git a/src/BODY/pair_body_rounded_polyhedron.cpp b/src/BODY/pair_body_rounded_polyhedron.cpp
index bc94b9cc805cc70cbb2f2f87a2ee015fcb9d10c2..02c2abe1e337edf619f310a583b8faa48d4a3f1e 100644
--- a/src/BODY/pair_body_rounded_polyhedron.cpp
+++ b/src/BODY/pair_body_rounded_polyhedron.cpp
@@ -47,7 +47,7 @@ using namespace MathConst;
#define DELTA 10000
#define EPSILON 1e-3
#define MAX_FACE_SIZE 4 // maximum number of vertices per face (same as BodyRoundedPolyhedron)
-#define MAX_CONTACTS 16 // for 3D models
+#define MAX_CONTACTS 32 // for 3D models
//#define _POLYHEDRON_DEBUG
@@ -308,36 +308,8 @@ void PairBodyRoundedPolyhedron::compute(int eflag, int vflag)
// also consider point contacts and line contacts
if (num_contacts > 0) {
- double contact_area;
- if (num_contacts == 1) {
- contact_area = 0;
- } else if (num_contacts == 2) {
- contact_area = num_contacts * A_ua;
- } else {
- int m;
- double xc[3],dx,dy,dz;
- xc[0] = xc[1] = xc[2] = 0;
- for (m = 0; m < num_contacts; m++) {
- xc[0] += contact_list[m].xi[0];
- xc[1] += contact_list[m].xi[1];
- xc[2] += contact_list[m].xi[2];
- }
-
- xc[0] /= (double)num_contacts;
- xc[1] /= (double)num_contacts;
- xc[2] /= (double)num_contacts;
-
- contact_area = 0.0;
- for (m = 0; m < num_contacts; m++) {
- dx = contact_list[m].xi[0] - xc[0];
- dy = contact_list[m].xi[1] - xc[1];
- dz = contact_list[m].xi[2] - xc[2];
- contact_area += (dx*dx + dy*dy + dz*dz);
- }
- contact_area *= (MY_PI/(double)num_contacts);
- }
rescale_cohesive_forces(x, f, torque, contact_list, num_contacts,
- contact_area, k_nij, k_naij, facc);
+ k_nij, k_naij, facc);
}
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,
@@ -383,6 +355,8 @@ void PairBodyRoundedPolyhedron::settings(int narg, char **arg)
mu = force->numeric(FLERR,arg[2]);
A_ua = force->numeric(FLERR,arg[3]);
cut_inner = force->numeric(FLERR,arg[4]);
+
+ if (A_ua < 0) A_ua = 1;
}
/* ----------------------------------------------------------------------
@@ -1381,7 +1355,17 @@ int PairBodyRoundedPolyhedron::interaction_face_to_edge(int ibody,
#endif
} else if (interact == EF_INTERSECT_INSIDE) {
-
+ // need to do something here to resolve overlap!!
+ // p is the intersection between the edge and the face
+ int jflag = 1;
+ if (d1 < d2)
+ pair_force_and_torque(jbody, ibody, xpj1, hi1, d1, contact_dist,
+ k_n, k_na, shift, x, v, f, torque, angmom,
+ jflag, energy, facc);
+ else
+ pair_force_and_torque(jbody, ibody, xpj2, hi2, d2, contact_dist,
+ k_n, k_na, shift, x, v, f, torque, angmom,
+ jflag, energy, facc);
}
return interact;
@@ -1661,6 +1645,7 @@ void PairBodyRoundedPolyhedron::contact_forces(int ibody, int jbody,
ft[1] = -c_t * vt2;
ft[2] = -c_t * vt3;
+ // these are contact forces (F_n, F_t and F_ne) only
// cohesive forces will be scaled by j_a after contact area is computed
// mu * fne = tangential friction deformation during gross sliding
// see Eq. 4, Fraige et al.
@@ -1697,13 +1682,49 @@ void PairBodyRoundedPolyhedron::contact_forces(int ibody, int jbody,
void PairBodyRoundedPolyhedron::rescale_cohesive_forces(double** x,
double** f, double** torque, Contact* contact_list, int &num_contacts,
- double contact_area, double k_n, double k_na, double* facc)
+ double k_n, double k_na, double* facc)
{
int m,ibody,jbody;
- double delx,dely,delz,fx,fy,fz,R,fpair,r,shift;
+ double delx,dely,delz,fx,fy,fz,R,fpair,r,shift,contact_area;
+
+ int num_unique_contacts = 0;
+ if (num_contacts == 1) {
+ num_unique_contacts = 1;
+ contact_area = 0;
+ } else if (num_contacts == 2) {
+ num_unique_contacts = 2;
+ contact_area = num_contacts * A_ua;
+ } else {
+ find_unique_contacts(contact_list, num_contacts);
+
+ double xc[3],dx,dy,dz;
+ xc[0] = xc[1] = xc[2] = 0;
+ num_unique_contacts = 0;
+ for (int m = 0; m < num_contacts; m++) {
+ if (contact_list[m].unique == 0) continue;
+ xc[0] += contact_list[m].xi[0];
+ xc[1] += contact_list[m].xi[1];
+ xc[2] += contact_list[m].xi[2];
+ num_unique_contacts++;
+ }
+
+ xc[0] /= (double)num_unique_contacts;
+ xc[1] /= (double)num_unique_contacts;
+ xc[2] /= (double)num_unique_contacts;
+
+ contact_area = 0.0;
+ for (int m = 0; m < num_contacts; m++) {
+ if (contact_list[m].unique == 0) continue;
+ dx = contact_list[m].xi[0] - xc[0];
+ dy = contact_list[m].xi[1] - xc[1];
+ dz = contact_list[m].xi[2] - xc[2];
+ contact_area += (dx*dx + dy*dy + dz*dz);
+ }
+ contact_area *= (MY_PI/(double)num_unique_contacts);
+ }
+
double j_a = contact_area / (num_contacts * A_ua);
if (j_a < 1.0) j_a = 1.0;
-
shift = k_na * cut_inner;
for (m = 0; m < num_contacts; m++) {
@@ -2302,6 +2323,41 @@ void PairBodyRoundedPolyhedron::total_velocity(double* p, double *xcm,
vi[2] = omega[0]*r[1] - omega[1]*r[0] + vcm[2];
}
+/* ----------------------------------------------------------------------
+ Determine the length of the contact segment, i.e. the separation between
+ 2 contacts, should be extended for 3D models.
+------------------------------------------------------------------------- */
+
+double PairBodyRoundedPolyhedron::contact_separation(const Contact& c1,
+ const Contact& c2)
+{
+ double x1 = 0.5*(c1.xi[0] + c1.xj[0]);
+ double y1 = 0.5*(c1.xi[1] + c1.xj[1]);
+ double z1 = 0.5*(c1.xi[2] + c1.xj[2]);
+ double x2 = 0.5*(c2.xi[0] + c2.xj[0]);
+ double y2 = 0.5*(c2.xi[1] + c2.xj[1]);
+ double z2 = 0.5*(c2.xi[2] + c2.xj[2]);
+ double rsq = (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1) + (z2 - z1)*(z2 - z1);
+ return rsq;
+}
+
+/* ----------------------------------------------------------------------
+ find the number of unique contacts
+------------------------------------------------------------------------- */
+
+void PairBodyRoundedPolyhedron::find_unique_contacts(Contact* contact_list, int& num_contacts)
+{
+ int n = num_contacts;
+ for (int i = 0; i < n - 1; i++) {
+
+ for (int j = i + 1; j < n; j++) {
+ if (contact_list[i].unique == 0) continue;
+ double d = contact_separation(contact_list[i], contact_list[j]);
+ if (d < EPSILON) contact_list[j].unique = 0;
+ }
+ }
+}
+
/* ---------------------------------------------------------------------- */
void PairBodyRoundedPolyhedron::sanity_check()
diff --git a/src/BODY/pair_body_rounded_polyhedron.h b/src/BODY/pair_body_rounded_polyhedron.h
index d7ee1f013e66a837948bdd6bddfb5464630e50fe..05985ef813a82790aebd56740525c102f2749130 100644
--- a/src/BODY/pair_body_rounded_polyhedron.h
+++ b/src/BODY/pair_body_rounded_polyhedron.h
@@ -35,12 +35,13 @@ class PairBodyRoundedPolyhedron : public Pair {
double init_one(int, int);
struct Contact {
- int ibody, jbody; // body (i.e. atom) indices (not tags)
- int type; // 0 = VERTEX-FACE; 1 = EDGE-EDGE
- double fx,fy,fz; // unscaled cohesive forces at contact
- double xi[3]; // coordinates of the contact point on ibody
- double xj[3]; // coordinates of the contact point on jbody
+ int ibody, jbody; // body (i.e. atom) indices (not tags)
+ int type; // 0 = VERTEX-FACE; 1 = EDGE-EDGE
+ double fx,fy,fz; // unscaled cohesive forces at contact
+ double xi[3]; // coordinates of the contact point on ibody
+ double xj[3]; // coordinates of the contact point on jbody
double separation; // contact surface separation
+ int unique;
};
protected:
@@ -124,7 +125,11 @@ class PairBodyRoundedPolyhedron : public Pair {
int jflag, double& energy, double* facc);
void rescale_cohesive_forces(double** x, double** f, double** torque,
Contact* contact_list, int &num_contacts,
- double contact_area, double k_n, double k_na, double* facc);
+ double k_n, double k_na, double* facc);
+
+ double contact_separation(const Contact& c1, const Contact& c2);
+
+ void find_unique_contacts(Contact* contact_list, int& num_contacts);
void sum_torque(double* xm, double *x, double fx, double fy, double fz, double* torque);
int opposite_sides(double* n, double* x0, double* a, double* b);