Updated pair body rounded/polygon and rounded/polyhedron (4bd4b2a1) · Commits · 郑智淋 / lammps

src/BODY/pair_body_rounded_polygon.cpp

+4 −2

Original line number	Diff line number	Diff line
		@@ -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;

src/BODY/pair_body_rounded_polyhedron.cpp

+90 −34

Original line number	Diff line number	Diff line
		@@ -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 += (dxdx + dydy + 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 += (dxdx + dydy + 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()

src/BODY/pair_body_rounded_polyhedron.h

+11 −6

Original line number	Diff line number	Diff line
		@@ -41,6 +41,7 @@ class PairBodyRoundedPolyhedron : public Pair {
		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);

Admin message